{'fraction_non_alphanumeric': 0.028123783573374854, 'fraction_numerical': 0.03474114441416894, 'mean_word_length': 4.878861556064074, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling novelty detection in the thalamocortical loop

In complex natural environments, sensory systems are constantly exposed to a large stream of inputs. Novel or rare stimuli, which are often associated with behaviorally important events, are typically processed differently than the steady sensory background, which has less relevance. Neural signatures of such differential processing, commonly referred to as novelty detection, have been identified on the level of EEG recordings as mismatch negativity and the level of single neurons as stimulus-specific adaptation. Here, we propose a multi-scale recurrent network with synaptic depression to explain how novelty detection can arise in the whisker-related part of the somatosensory thalamocortical loop. The architecture and dynamics of the model presume that neurons in cortical layer 6 adapt, via synaptic depression, specifically to a frequently presented stimulus, resulting in reduced population activity in the corresponding cortical column when compared with the population activity evoked by a rare stimulus. This difference in population activity is then projected from the cortex to the thalamus and amplified through the interaction between neurons of the primary and reticular nuclei of the thalamus, resulting in spindle-like, rhythmic oscillations. These differentially activated thalamic oscillations are forwarded to cortical layer 4 as a late secondary response that is specific to rare stimuli that violate a particular stimulus pattern. Model results show a strong analogy between this late single neuron activity and EEG-based mismatch negativity in terms of their common sensitivity to presentation context and timescales of response latency, as observed experimentally. Our results indicate that adaptation in L6 can establish the thalamocortical dynamics that produce signatures of SSA and MMN and suggest a mechanistic model of novelty detection that could generalize to other sensory modalities. Author summary Cortical sensory neurons have been shown to be capable of novelty detection, that is they respond more vigorously when a novel, unexpected stimulus is presented, and less so when the stimulus is part of a predictable sequence. However, the neural mechanism underlying this capability is not yet fully understood. Here, we developed a thalamocortical network model that accounts for novelty detection and reproduces physiologically observed neural response patterns in the somatosensory cortex. Specifically, our results demonstrate that the novelty signal arises from the complex recurrent interplay between thalamic neurons and cortical neurons in layers 4 and 6. This work therefore provides a concrete mechanism that can serve as a starting point for further investigating the neural circuit mechanisms underlying novelty detection.

Sensory cortices are adept at identifying regularities and patterns obscured in the 2 ever-changing input stream to dynamically generate predictions for the forthcoming 3 stimulus. Their responses are then updated according to differences between the 4 anticipated sensory inputs and the actual ones [1]. Mismatch negativity (MMN), 5 observed in human electroencephalograph (EEG) recordings [2], and stimulus-specific 6 adaptation (SSA), evident in single-cell recordings in animals [3], are two prominent 7 instances of this type of short-term plasticity acting over timescales of seconds to 8 minutes. For example, when consecutively presented with an identical stimulus over a 9 short period, the cortical sensory neuron will decline/adapt its response and 10 predict/expect another same one for the next stimulus. If the next stimulus violates the 11 expectation, a larger/unadapted response will be evoked than by the expected stimulus. 12 Therefore, both SSA and MMN conceptually represent a form of novelty detection: 13 unexpected or infrequent stimuli are made more salient by eliciting stronger responses 14 than expected or frequent ones. 15 The phenomenon of attenuated neuronal responses to common (standard) stimuli 16 without generalizing to other, rare (deviant) stimuli has been a long-standing topic of 17 interest, especially in the auditory system, which is referred to as SSA [4][5][6][7][8]. It is 18 typically tested using the oddball paradigm that contains a repetitive standard sequence 19 occasionally interrupted by a deviant [9]. SSA is a robust and widespread phenomenon 20 found in multiple sensory modalities. Ulanovsky et al [4] first observed this phenomenon 21 in the auditory cortex, and following this pioneering work, SSA has also been identified 22 at subcortical stages of the auditory pathway such as the thalamus [10][11][12] and the 23 inferior colliculus [13,14], as well as in other sensory systems such as the somatosensory 24 cortex [5,8] and the visual system [15,16]. Furthermore, SSA was demonstrated in the 25 auditory cortex of anaesthetized [7,17], awake [18] and freely moving rats [19], 26 indicating that this phenomenon is not significantly affected by brain state and likely a 27 hardwired component of sensory processing. 28 MMN, like SSA, is usually studied by applying the oddball paradigm and manifested 29 in human EEG recordings as an additional negative deviation elicited by a deviant 30 stimulus that breaks the regularity established by repeated presentations of the 31 standard [9]. Both cortical SSA and MMN demonstrate a "true" deviance-detecting 32 property, where the deviant response is dependent on its presentation context and 33 beyond mere sensory input depression by applying several control paradigms [6,8,20]. 34 In addition to SSA, animal models also show MMN-like responses that are homologous 35 with human MMN in terms of their electrophysiological, pharmacological and functional 36 properties [21,22]. Despite significant similarities, SSA is unlikely to be the direct cortex (S1) occurring hundreds of milliseconds post-stimulus, similar to MMN latency. 49 The late sensory response has also been tested to be context-specific. These evidence 50 strongly implied that the late response possibly lies in the physiological origin of MMN. 51 Over the past decade, a large number of computational models have been proposed 52 to account for cortical SSA, with a focus on the auditory pathway. One simple type of 53 SSA model is based on a feedforward network with synaptic depression [17,[28][29][30][31] 54 which, unlike general neuronal fatigue, is stimulus-specific. Prominent synaptic 55 depression is found in the thalamocortical pathway, and thus cortical SSA emerges due 56 to the differentially adapted thalamocortical inputs [17,28]. In this class of models, the 57 standard's regularity cannot affect the deviant response because the responses to each 58 stimulus feature solely depend on the amount of adaptation load accumulated in its 59 input channel. Hence the model fails to exhibit true deviance detection in which the 60 presentation context formed by the standards also contributes to the deviant response. 61 Nevertheless, this flaw of the model can be addressed through some modifications, such 62 as employing the model in a multi-layer configuration [30] or assuming the width of the 63 tuning curve of the input channel is history-dependent [3]. On top of that, recurrent 64 networks have been considered an alternative way to generate cortical SSA [32][33][34][35][36]. 65 Predictive coding based neuronal network of the auditory cortex was studied by 66 Wacongne et al. to account for deviance detection in the context of MMN [32]. Yarden 67 et al. proposed a generic SSA mechanism mediated via the propagation of synchronous 68 population activity across a local neural circuit [35]. Recently a laminar network of 69 morphologically plausible, multi-compartmental neuron units was designed to capture 70 SSA in the form of local field potentials [36]. 71 However, there is still no computational model accounting for the biphasic sensory 72 responses in S1 that can be either sensitive to the rarity and presentation context of 73 stimuli or not, depending on the response latency [8]. In particular, the latency of the 74 secondary cortical response cannot be accounted for by the dynamics of a recurrent 75 network involving short-term synaptic depression that is commonly regarded as a 76 indispensable constituent of SSA models. Thus the latency is possibly produced by 77 some subcortical source(s) other than the intra-cortical dynamics. Inspired by the idea 78 that SSA arises from the interaction between differentially adapted populations of 79 neurons tuned for specific feature stimuli [35], here we developed a hierarchal 80 thalamocortical model for the whisker-related region of S1 (also known as the barrel  , which is somatotopically arranged into discrete clusters called "barreloids", relays peripheral sensory information to the L4 barrels (purple), with the relative strength of connections represented by the size of the arrows. Each L4 cluster projects its output vertically to the corresponding L6 infrabarrel (green), which in turn provides feedback excitation to its somatotopic barreloid and thalamic reticular nucleus (TRN, in blue shade), with higher cortical drive on TRN than VPM (specified by the size of the diamonds). Each L4 barrel and L6 infrabarrel also connects to its neighbouring barrels. (B) Population spike (PS) can be generated in the cortical cluster, represented by the transient increase of the cluster's mean firing rate (left). Mean synaptic resources of the neuronal group are depleted by the PS and later recover gradually (right). Higher levels of initial resources evoke a more substantial PS. The time 0 represents the onset of a stimulus whose duration is marked in grey shading. (C) Two classes of cells with different synaptic receptors are used in the thalamic circuit: thalamocortical (TC) relay cells with excitatory AMPA-mediated synapses in the VPM and reticular (RE) cells with inhibitory GABA A -mediated synapses in the TRN. TC and RE neuronal populations are mutually coupled, and the RE population is also recurrently connected. The thalamus can be activated by bottom-up sensory input (grey) and top-down cortical drive (green). (D) A minimal thalamic network of 2 TC and 2 RE cells driven by the cortex is capable of inducing spindle oscillations. (E) Voltage traces of two types of thalamic cells simulated by the izhikevich model in response to a step in direct current (bottom of each panel): initial burst tonic firing in TC neuron (-70 mV at rest, hyperpolarizing pulse of 10 pA) and rebound burst in RE neuron (-62.5mV at rest, depolarizing pulse of -10 pA lasting 120 ms).
(F) The time course of membrane potentials for the 4 thalamic cells. Cortical stimulation is marked by red bars under the traces (0.08 pA for TC 1 and 20 pA for RE 1 ).

98
The architecture of the somatosensory thalamocortical network is organized in a loop as 99 shown in Fig 1A. Based on the stereotypic somatotopic map for representing rodent 100 whiskers, the barrel cortex is modeled as a grid of interconnected vertical columns, each 101 of which primarily mediates stimuli from its principal whisker. Also, we only focus on 102 neuronal clusters within layers 4 (barrels) and 6 (infrabarrels) of each cortical column, 103 which respectively act as the afferent and efferent layers of the cortex from and to the 104 thalamus [37][38][39]. The barreloids in the ventral posterior medial (VPM) nucleus of the 105 thalamus show an identical layout to that of the barrels. Each barreloid relays the 106 information from its principal whisker through an ascending sensory pathway to L4 of 107 its principal as well as surrounding barrels, with stronger coupling to the principal than 108 the surrounding ones. We hypothesize that L2/3 and L5 contribute little to the cortical 109 SSA but merely transmit columnar information from L4 to L6. Therefore, a direct 110 shortcut connection from L4 to L6 within every cortical column is added, which is 111 sufficient in the simplified model proposed here to demonstrate SSA in L6, although in 112 reality the flow of excitation from L4 will go through L3, L2 and L5 in sequence before 113 arriving at L6 [40][41][42]. To close the thalamocortical loop, the outputs of layer 6 are 114 projected back to their somatotopically aligned barreloids and thalamic reticular 115 population. Here we assume that neurons in the thalamic reticular nucleus (TRN) are 116 also whisker-specific, but they are not necessarily arranged in spatially discrete clusters 117 as barreloids. 118 We first focus only on the cortical network that is capable of generating SSA and 119 true deviance detection by modulation of cross-barrel synchronous activities. All 120 cortical intra-and inter-column synapses exhibit short-term depression. Only excitatory 121 populations are considered in the cortical circuitry to simplify the mean-field analysis, 122 and the level of excitability can be regulated by the threshold of neuronal gain function. 123 Every layer-specific population in columns is modelled by a mean-field recurrent 124 network coupled with depressing synapses, which can generate the so-called population 125 spike (PS) in response to external stimulation. The PS is characterized by nearly 126 synchronous firing of a group of neurons within a short time window [43,44]. In the 127 cortical mean-field rate populations we studied here, a population spike is described as a 128 sharp increase in population activity, leading to fast depletion of mean synaptic resource 129 that will recover gradually after the removal of external input. Higher initial synaptic 130 resources cause more substantial PSs to be triggered in response to the same stimulus 131 (Fig 1B).

132
The thalamic circuitry is set up to induce spindle oscillations characterized by 133 November 6, 2021 5/29 intermittent bursting activity at 7-14 Hz, which is postulated as the origin of the late 134 rhythmic responses in L4 of S1 [8], since the frequency of cortical oscillation (∼10Hz) 135 falls into the spindle range. Early studies have shown that network and intrinsic 136 mechanisms act in combination to generate spindle oscillations [45]. Two interacting 137 populations of spiking thalamic neurons are considered to model spindle oscillation: 138 excitatory thalamocortical (TC) neurons and inhibitory thalamic reticular (RE) 139 neurons. As illustrated in Fig 1C, the TC population is excited by whisker stimulation 140 and then provides AMPA-mediated excitatory synaptic input to the RE population and 141 neurons in cortical L4. In turn, the RE neurons suppress the activity of the TC neurons 142 and themselves via GABA A -mediated inhibition. The top-down modulation from L6 143 exerts more effective excitation on RE than TC population due to more substantial 144 corticothalamic conductance [38,45,46]. Unlike the mean-field dynamics of the cortex 145 represented by the rate network, the intrinsic dynamics of individual spiking thalamic spike-bursts upon release from hyperpolarization [45,47] (Fig 1E). 150 We elucidate the mechanism of spindle oscillation triggered by the corticothalamic 151 feedback excitation in a minimal thalamic circuit consisting of 2 RE and 2 TC cells [46] 152 ( Fig 1D). When the cortical drive is exerted on both TC 1 and RE 1 neurons, only the 153 RE 1 neuron is able to elicit spikes due to the considerably higher cortical conductance 154 on RE than TC neurons. The bursting activity of RE 1 cell provide inhibitory 155 postsynaptic potentials (IPSPs) to both TC cells and then induce rebound bursting 156 spikes in TC 2 cells when its potential is released from hyperpolarization and almost 157 reaches its resting state. The firing of TC 2 cell reciprocally excite RE cells, leading to 158 the next cycle of spindle oscillations ( Fig 1F). Minute bump in the potential of TC 1 cell 159 evoked by initial corticothalamic stimulation prevents the cell from generating rebound 160 bursts in the first cycle, which results in the TC cells fire every two cycles while RE 161 cells fire every cycle in the minimal circuit. With the thalamocortical pathway, spindle 162 oscillation is synchronized over cortical areas as hypothesized to be the late oscillatory 163 rhythms found in L4 of S1 [8].

164
Propagation of population activities in the thalamocortical loop 165 The PS evoked in the principal L4 barrel will propagate to the neighbouring L4 barrels 166 via inter-column coupling, as well as to the column-aligned L6 infrabarrel via the 167 L4-to-L6 projection [40]. The L4 responses are largely localized around the principal 168 barrel and substantially attenuated in the surrounding barrels [48,49] ( after stimulus onset, in accordance with experimental findings [50] (Fig 2B, top). The 171 topographical spread of PS in L6 is similar to that in L4 but with broader temporal 172 profiles of the PSs and slightly longer response latencies [50] (Fig 2A and B, bottom individual whisker from the deflection probability on SSA [8].

201
In accordance with experimental findings [8], the L4 network exhibits prominent accordance with the sporadic standard-triggered oscillation recorded in L4 of the 226 cortex [8]. Here these late rhythmic responses occur largely because the deviant did not 227 deplete as many resources in L6 of D2 as another standard would have done, which 228 occasionally allows strong enough cortical drive on the thalamus to generate a spindle 229 oscillation during the presentation of the next standard stimulation.

230
Context-dependent deviance detection 231 It has been suggested that the larger responses to the stimuli features, when deviant 232 than when standard, not only reflect the rarity of the deviant causing less 233 use-dependent adaptation [3], but also the sensitivity to the violation of the expectation 234 established by the repeated standard, indicating a dynamic prediction mechanism for 235 the next stimulus based on the short-term memory of the previous statistical sensory 236 stream [52]. Such sensitivity to contextual information is termed true deviance 237 detection/sensitivity. A common test for true deviance detection uses many-standards 238 control protocol, where besides the standard-deviant stimuli pair in the oddball 239 sequence, many other different stimuli are presented as the standards [20]. In this 240 control, each of these distinct stimuli occurs with the same probability as the deviant 241 stimulus, which consequently eliminates any potential expectation for the next stimulus 242 but preserves the rarity of the deviant. If the response to the deviant stimulus in the 243 oddball sequence exceeds the response to the same deviant one embedded in the control 244 sequence, then the existence of true deviance detection is confirmed [20]. Here in the 245 many-standards condition for the somatosensory domain, the standard stimuli are 246 equiprobably distributed over three whiskers in a row that are adjacent to the deviant 247 whisker, therefore each whisker is deflected with 25% probability [8] (Fig 4A).  activity is elicited in the many-standards than the oddball condition owing to the lower 287 level of L6 resources reduced by the last D3 stimulation in the many-standards condition 288 (Fig 5C). In consequence, only the corticothalamic feedback in the oddball protocol is 289 large enough to trigger spindle oscillations in TC cells (Fig 5E), which in turn cause 290 synchronized activity in the cortical L4 barrel (Fig 5B, top). In a similar vein, the   317 Furthermore, our model predicts the dependence of the novelty-predicting effect on 318 the inter-stimulus interval (ISI) across different cortical layers and thalamus (Fig 7).

319
Due to the feature of biphasic dynamics in our network, CSI is separately computed for 320 the early and late responses. For the early response, the novelty-predicting effect is 321 demonstrated in the cortex over a wide range of ISI and particularly prominent in L6 322 that as we expected to be the novelty detector. Early novelty-predicting responses 323 originated in L6 are subsequently revealed in the thalamic oscillation and secondary L4 324 responses via the L6-to-thalamus-to-L4 pathway.

325
To eliminate the potential effect of sensory input uncertainties on CSI, we removed 326 the stimulus-by-stimulus variation of the thalamocortical drive that was previously 327 presented and reran both protocols on the model. In this case, each external stimulus 328 evokes the same amplitude of thalamic responses. One can clearly see from Fig 8A that 329 the novelty-predicting effect is greatly attenuated near both ends of the range allowing 330 for it because either of the cases will lead to the steady-state of unrecovered/recovered 331 synaptic resources in L6, which are independent of the presentation context. We also 332 investigate the individual role of three types of deviant in the sequenced paradigm in 333 detecting novelty across different ISIs, where the C1 deviant preceded by the B1, C2 334 and D1 standard is defined as type 1, 2 and 3, respectively (Fig 8B). The 335 novelty-predicting effect is the most pronounced for the type 2 deviant and is not 336 influenced much by the increase in ISI after the first several values, which indicates the 337 capacity of novelty prediction in the cortex depends on the specific sequential context in 338 the recent past in randomized paradigms, i.e., the precise information about the order of 339 standard presentations over hundreds of milliseconds before the arrival of deviant.

341
We presented a cortico-thalamo-cortical circuit of somatosensory processing with 342 recurrent dynamics and depressing synapses that is capable of producing both SSA and 343 the contextual sensitivity over different latencies across cortical layers and thalamus, in 344 line with previous experimental observations [8]. The computational results suggest that 345 the early context-specific deviant responses in S1 L6 could be a strong candidate for the 346 electrophysiological correlates of the mid-latency response observed in EEG 347 recordings [22,26], because of their similar response latencies and shared sensitivity to 348 presentation context. Besides, we show that the late context-dependent responses in S1 349 L4, which temporally and functionally resemble MMN, can be produced simply by the 350 interaction between different cortical layers and the thalamus, without involving 351 top-down modulation by higher brain areas.  [54,55], where the mean firing rates show a temporary peak in the principal 360 column upon the onset of stimulation, followed by weaker activity with short latency in 361 the surrounding columns [55]. Previous modelling work has demonstrated that this 362 population burst could be a consequence of a recurrent network with synaptic 363 depression [43,44,56,57]. Here, such a network is the essential ingredient of our 364 hierarchical model to produce SSA.  work has demonstrated that such a mean-field representation is equivalent to a 374 homogeneous pool of leaky integrate-and-fire neurons in terms of its population 375 activities [58]. Second, this specification allows us to focus on the collective population 376 behaviour within the grid network, which is sufficient to generate cortical SSA and true 377 deviance detection. SSA and true deviance sensitivity cannot be accounted for solely by 378 the intrinsic dynamics of individual neurons, but relies to a large extent on network 379 dynamics. It should be noted that although this mean-field network of the barrel cortex 380 ignores some details, such as the potential contribution of additional layers, distinct and their interplay to elicit synchronous oscillations. Specifically, the spiking dynamics 392 of each thalamic neuron is simulated by an Izhikevich model that reproduces similar 393 firing patterns to the conductance-based model used in [46], but is computationally as 394 efficient as the leaky integrate-and-fire model. Our modelling results implicate the slow 395 recovery variable in these thalamocortical neurons, which may describe the inactivation 396 of low-threshold Ca 2+ T-currents [59], as the main factor driving the latency of the 397 secondary response observed in S1 L4 and therefore potentially MMN.

398
Finally, it is worth mentioning that the current implementation of our model still 399 has several shortcomings. A major limitation is the sensitivity of SSA to cortical 400 network parameters, which stems from the different layer-specific response patterns in 401 the feedforward L4-to-L6 network. In other words, the parameter regime allowing for 402 SSA and the true deviance sensitivity in initial L6 responses is constrained by not only 403 the spatiotemporal pattern of whisker-evoked transient activity in the L6 network but 404 also responses in L4 that directly delivers its output to L6. Primarily, in the L6 circuit, 405 the PS evoked by the oddball standard should reduce with repeated stimuli while not 406 generalize to the oddball deviant. Besides, the propagation of activity across the L6 407 network is supposed to be robust enough to occasionally pass the PSs triggered by 408 standards in the many-standards condition into the deviant column, yet too weak for 409 PSs evoked by the oddball standard to do so. Another possible shortcoming of the 410 model is that it omitted cortical layer 2/3 and layer 5 that could modulate the laminar 411 pathway from L4 to L6. Future work incorporating these missing layers could generate 412 richer dynamics that expand the narrow region of SSA existence in the parameter space, 413 leading to a robust SSA phenomenon. Lastly, it would be instructive to test the NMDA 414 dependency of the late response in a spiking implementation of the cortical dynamics 415 explored here to further investigate SSA's connection with MMN.

417
Derivation of spiking to rate network 418 We derive the rate equations for a network of interacting populations mediated by 419 depressing synapses assuming that each population is large and homogeneous, using 420 standard techniques [58,[60][61][62]. By homogenous, we mean that (i) all neurons have the 421 same properties; (ii) each neuron receives identical external input; and (iii) the coupling 422 strength between any pair of neurons within the population is statistically uniform. The 423 general conclusions developed here will be applied in the next section to formalize the 424 mean-field cortical network composed of a grid arrangement of coupled populations. 425 We define the mean firing rate as the population activity A(t) averaged over a pool 426 where N is the number of neurons within the population, n act (t; t + ∆t) denotes the 428 number of spikes of all neurons in the population, which appear betweem time t and 429 t + ∆t. Since ∆t is a very short time interval, each neuron fires at most one spike during 430 the interval and correspondingly n act (t; t + ∆t) also represents the number of neurons 431 that fire between time t and t + ∆t. δ(t) is the Dirac delta function and t f j denotes the 432 firing time of the j-th neuron in the population. This notion of the rate as population 433 activity averaged over many neurons will be widely used in the following derivation.

434
The dynamics of short-term depression for the connection between pre-and 435 post-synaptic neuron pair j, i is given by the phenomenological model [56,63]: where each pre-synaptic spike emitted at t f j utilizes a certain rate U of the available 437 synaptic resource x ij that is a fraction of full resources. The depleted available resource 438 then returns to its baseline value with a time constant τ rec . For the current-based type 439 of post-synaptic input, which is independent of neuronal membrane potential, the 440 increment in the amplitude of post-synaptic current I i (t) triggered by a pre-synaptic 441 spike arriving at t f j is given by where J is the absolute synaptic efficacy and x − ij (t f j ) denotes the value of x ij 443 immediately before the arrival of the spike. 444 We define x i (t) ≡ 1/N N j=1 x ij as the mean available resource in a bunch of 445 synapses from a population of pre-synaptic neurons j to a post-synaptic neuron i. By 446 taking the average over the pool of pre-synaptic neurons on both sides of Eq 2, we get 447 The second term on the right-hand side of Eq 4 can be further transformed, assuming 448 uncorrelated Poisson spiking from a large population of pre-synaptic neurons, to where N is the number of neurons in the pre-synaptic population. Ω denotes the 450 subpopulation of pre-synaptic neurons that fire between time t and t + ∆t, and 451 n act (t; t + ∆t) is the size of the subpopulation. The approximation becomes an equality 452 in the limit of N → ∞. A(t) represents the population activity of pre-synaptic neurons 453 as defined in Eq 1.

454
Using Eq 5, Eq 4 can be reformulated as November 6, 2021 13/29 We find that the dynamics of the mean synaptic resource x i (t) is independent of the 456 neuronal identity and therefore is the same for all neurons in the population.

457
In addition, we assume that each pre-synaptic spike emitted at t f j triggers an 458 instantaneous post-synaptic current with the temporal profile α(t − t f j ) and the efficacy 459 for each synaptic connection in an all-to-all coupled population of N neurons is scaled 460 as J 0 /N . The input current to a neuron i in the population is generated by all spikes of 461 all neurons within the same population, Making use of Eq 5 with t replaced by t − s to substitute the average quantity on 463 the right-hand side of Eq 7, we obtain 464 We dropped the neuronal index i since the input current at time t depends on past 465 mean synaptic resource and population activity and is identical for all neurons.
where N k is the size of population k and Γ k represents the set of neurons that belongs 470 to population k. 471 We assume that each neuron i in population k receives input from all neurons j in Correspondingly the dynamics of mean synaptic resource x i (t) is generalized as By assuming that U kl and τ rec,kl are homogenous across different pools, we drop the 480 subscript kl and reformulate the above equation as Here we change the index of mean synaptic resource from post-synaptic neuron i in 482 population k to pre-synaptic population l since the evolution of mean synaptic resource 483 is governed by the activity of pre-synaptic population and is the same for all neurons 484 within post-synaptic populations.

485
In this case, the input current, which is independent of the neuronal index i but the post-synaptic population index k. 487 Finally, we characterize the dynamics of the input current h k (t) of population k with 488 the differential equation of passive membrane and employ for each population the rate 489 model A l (t) = F l (h l (t)) in which F l is the stationary gain function of neurons in 490 population l. The input current h k (t) takes into account both synaptic coupling I k (t) 491 and external drive I ext where all neurons in the network have the same membrane time constant τ m .

493
Particularly, in the case of instantaneous post-synaptic current pulse α kn (s) = δ(s), 494 Eq 14 can be reduced to a first-order differential equation Modelling the thalamocortical network of the barrel cortex 496 It has been observed that a cortico-thalamo-cortical loop structure exists in the rodent

504
To mimic the somatotopic organization of the barrel cortex, the circuit of each layer 505 is composed of an array of cortical columns (barrels). Each column is modelled by a 506 fully connected excitatory population mediated by synaptic depression, with its 507 mean-field dynamics as analyzed in the last subsection. Every column is also connected 508 to its nearest vertical and horizontal neighbours by inter-column depressing synapses.

509
Furthermore, the L4 columns receive thalamic inputs with cross-whisker tuning and 510 send outputs to their aligned L6 columns.

511
Based on Eqs 15 and 12, the mean-field dynamics of cortical layer 4 network is described by the following equations (see Fig 9 for a schematic representation of the network architecture): Similarly, the mean-field cortical layer 6 network is defined as: The dynamics of each barrel is described by its population activity A m,n and mean 513 fraction of synaptic resource available x m,n , where m and n denote the row and arc where λ is the radius of tuning curve and d is the floored Euclidian distance between 537 barrel m, n and whisker channel p, q on the evenly-spaced gird where the separation 538 distance is 1 where · denotes the floor function that rounds its argument to the nearest integer less 540 than or equal to that argument.

541
A p,q TC describes the thalamocortical input delivered through the whisker channel p, q. 542 It is defined as the population activity of TC cells organized in barreloid m, n, which in 543 practice is evaluated by Eq 1 with a time window ∆t = 2 ms. It is worth noting that 544 the time bin of 2 ms is only used to compute the population activity of TC cells from 545 their spikes, and the cortical population activity is obtained by numerically solving 546 those differential equations of the cortical dynamics with an integration time step of 0.1 547 ms. Specifically, the A p,q TC is updated every time window (2 ms) and held constant for 20 548 integration time steps (each 0.1 ms). 549 We assume that the late cortical rhythmic activity (roughly 10 Hz) that lies in the 550 frequency range of the spindle oscillation (7-14 Hz) [45] emerges from the 551 thalamo-cortical system. Earlier studies have shown that spindle oscillation is initiated 552 in the thalamus and synchronized over the cortex and its minimal neural substrate lies 553 on the interaction between thalamic RE cells and TC relay cells [45,46].

554
Motivated by the anatomy of the VPM and reticular nucleus of the thalamus, we 555 modelled each barreloid in the VPM nucleus by a cluster of 100 TC cells mediated by a 556 pool of 100 RE cells, whose connectivity is illustrated in Fig 9B. The spiking dynamics 557 of each thalamic neuron simulated by the Izhikevich model, which uses the quadratic 558 integrated-and-fire model for the membrane potential equation, is expressed as with after-spike resetting where the time t is in milliseconds; v and u respectively denote the membrane potential 561 (in mV) and recovery variable that represents the difference of all influx and efflux of The thalamic background noise I noise (t) is subject to uniform distribution bound 574 between -0.5 and 0. The conductance-based synaptic current is modelled as where E syn is the synaptic reversal potential. Synaptic conductivity g syn (t) has a time 576 course of exponential decay: where Θ(t) is the Heaviside step function and t f represent the arrival time of pre-578 synaptic spikes. We simulated AMPA and GABA A receptors with the decay time 579 constant τ = 5 and 6 ms; E syn = 0 and -75 mV, respectively. The maximal conductance 580 g syn for each type of connection is scaled asḡ syn =Ḡ syn /C, where C = 100 × 0.6 = 60 is 581 the number of randomly chosen presynaptic partners for each neuron andḠ syn = 2, 0.01 582 and 0.2 nS for AMPA-mediated TC-to-RE, GABA A -mediated RE-to-TC and 583 GABA A -mediated RE-to-RE synapses, respectively.

584
In addition, half of TC cells in each barreloid and their coupled RE cells were 585 randomly picked receiving corticothalamic feedback current, which is expressed as where A L6 (t) denotes the population activity of L6 neurons in the somatotopic 587 infrabarrel. We respectively assigned the corticothalamic coupling strength 588 w CTh = 0.001 and 0.4 for TC and RE neurons.

589
Finally, we respectively stimulated randomly 10% and 50% of TC cells in the two 590 complementary subgroups of each barreloid with sensory input that is described as where B and ξ p,q (t) respectively represent the maximum amplitude and temporal context-dependent deviance detection. In the whisker oddball condition, deviant stimuli are applied to the C2 whisker with 25% probability of occurrence, while standard stimuli are applied to the D2 whisker with 75% probability. In the many-standards condition, the deviant stimulus is presented with the same probability of appearance as in the oddball condition, but standard stimuli are equiprobably distributed over three whiskers D1-D3, each of which is stimulated 25% of the time. The bar height represents the probability of a specific whisker being stimulated. (B) A subset of multi-unit (MU) recordings in the granular layer of rat S1 shows late deviant-specific responses in the whisker oddball and many-standards protocols. Average early spike (top row, left) and late deviant-specific oscillatory responses (top row, right) are shown for standard (blue curves, activity recorded in D2) and deviant (red, activity recorded in C2) stimuli in oddball condition. The late response (bottom row) also demonstrates a distinct difference between two deviant types (both responses recorded in C2) in oddball (red curves) and many-standards protocols (green) respectively. Data are adapted from [8]. (C) Average biphasic population activity in cortical L4 (top row), L6 (middle) and thalamus (bottom) in response to oddball standard (blue traces, recorded in D2), oddball deviant (red, recorded in C2) and many-standards deviant (green, recorded in D2) stimuli. Similar to (B), the left and right columns of panels show the early and late activities respectively across cortical layers and the thalamus. SSA and true deviance detection are initiated in the early L6 responses and subsequently enhanced in the late thalamic oscillation, which finally induces the secondary deviant-selective cortical response in L4. The third cycle of thalamic activity is too weak to produce a response in L4 since this thalamic input is below the rheobase (current threshold of gain function) of L4 neurons. Stimulation periods are marked in grey shading. Data in each protocol are averages of 120 stimulus presentations. In both stimulation paradigms, peripheral stimulus duration is 10 ms and inter-stimulus interval (ISI) is 1 s. Considerably larger early onset responses to the many-standards deviant compared to the oddball deviant are generated in the L6 population (see dashed boxes), because higher activity is propagated from the L6 standard column(s) into that of the deviant column in the many-standards compared to the oddball condition, and accordingly less synaptic resource is left and a weaker response is evoked upon arrival of the deviant. (D) Synaptic resources in L6 of D2 in response to the standard stimulus applied to D2 whisker embedded in oddball (solid traces) and many-standards (dashed) sequences. The standard stimulus to D2 in the many-standards condition, with its lower probability of appearance than the same stimulus in the oddball condition, allows the resource in L6 of D2 to have more time to recover fully before the next standard stimulus arrives and consequently to generate an adequate PS capable of spreading into the same layer of the deviant column. (E) Thalamocortical (TC) population activity in C2 barreloid evoked by the deviant in the oddball (red) and many-standards (green) protocols. The corticothalamic excitation of the deviant in the many-standards condition is normally strong enough to elicit late intermittent bursting, whereas that of the deviant in the oddball condition is not (dashed boxes). The time axes are aligned for panels A-E. In both stimulation paradigms, peripheral stimulus duration is 10 ms and inter-stimulus interval (ISI) is 1 s (onset-to-onset). In two sequences, the delivery of the deviant stimulus to the C1 whisker is the same, but three standard stimuli are employed either periodically (sequenced condition, top) or randomly (randomized condition, down) to C2, D1 and B1 whiskers. In the sequenced stimulus train, the deviant occasionally breaks the regularity by repeatedly taking the place of one of the standards in the order of C2, D1 and B1. In both paradigms, stimulus duration = 10 ms and inter-stimulus interval (ISI) = 250 ms (stimulus offset to onset). Deviant stimuli constitute 3% of overall stimuli (30 out of 1000 stimuli). (B) Average population activity across cortical L4 (left), L6 (middle) and thalamus (right) in response to the deviant in sequenced (red) and randomized (green) paradigms. A weakly differential response to sequenced and randomized deviant exclusively appears in the L6 population at short latency.
Owing to the thalamocortical loop, this subtle difference is successively manifested in late thalamic and cortical L4 bursting activities, which are magnified in insets for clarity. Stimulus duration is highlighted in grey. CSI values are 0.0115, 0.0254 and -0.0061 for early L6 (computed over [0, 40] ms) and late L4 and thalamic responses (computed over [40,250] ms), respectively.