Thalamic control of sensory processing and spindles in a biophysical somatosensory thalamoreticular circuit model of wakefulness and sleep

Summary Thalamoreticular circuitry plays a key role in arousal, attention, cognition, and sleep spindles, and is linked to several brain disorders. A detailed computational model of mouse somatosensory thalamus and thalamic reticular nucleus has been developed to capture the properties of over 14,000 neurons connected by 6 million synapses. The model recreates the biological connectivity of these neurons, and simulations of the model reproduce multiple experimental findings in different brain states. The model shows that inhibitory rebound produces frequency-selective enhancement of thalamic responses during wakefulness. We find that thalamic interactions are responsible for the characteristic waxing and waning of spindle oscillations. In addition, we find that changes in thalamic excitability control spindle frequency and their incidence. The model is made openly available to provide a new tool for studying the function and dysfunction of the thalamoreticular circuitry in various brain states.

In brief Iavarone et al. present a computational model of thalamoreticular microcircuitry providing insight into thalamic interactions during wakefulness and sleep. They find common mechanisms that regulate the frequency-dependent enhancement of sensory responses and fluctuations of spindle oscillations, and discover that thalamic excitability controls spindle frequency and incidence.

INTRODUCTION
The thalamus and thalamic reticular nucleus (Rt) lie at the heart of the thalamocortical (TC) system in the mammalian brain. 2 Thalamic relay cells send projections to the cortex and form excitatory collaterals with thalamic reticular neurons. These neurons then send inhibitory projections back to the thalamus, creating the thalamoreticular circuit. [3][4][5] This circuit plays a key role in several brain functions, such as transmitting sensory information to the cortex and regulating brain states such as sleep and wakefulness. It has also been linked to attentional processes and is responsible for the generation of spindle oscillations during sleep. [6][7][8][9][10] Changes in thalamic neuron firing and connectivity have been associated with abnormal brain rhythms, such as those seen in absence epilepsy. [11][12][13][14][15] Alterations in the incidence and density of spindle oscillations during sleep have been observed in various disorders, including schizophrenia, 16-20 neurodevelopmental disor-ders, 21 attention-deficit/hyperactivity disorder (ADHD), 22 and Alzheimer's disease. 23 In this work, we follow and extend the reconstruction pipeline presented in Markram et al. 24 to develop and validate a digital model of a thalamoreticular microcircuit of a portion of first-order somatosensory (ventral posterolateral nucleus [VPL]) thalamus in the adult mouse.
This study uses a computational model to examine the dynamics of the thalamoreticular circuit during wakefulness and sleep-like states. The model is able to replicate multiple experimental findings and provide new insights and predictions. In wakefulness, the model shows that the thalamoreticular circuitry can generate frequency-selective enhancement of thalamic responses. In sleep, the model reveals that the waxing and waning of spindle oscillations is generated intrathalamically, and that changes in cell excitability can alter spindle incidence and frequencies. We provide the experimental data and computational models as an open resource for further research. Electrical types Over 100 TC neurons, INs, and reticular neurons were characterized through patch-clamp recordings (see STAR Methods) in mouse brain slices. TC neurons were classified as adapting (cAD_ltb) and non-adapting (cNAD_ltb) e-types ( Figure 1B) based on their tonic firing responses, similar to those in rat TC cells. 26 Rt neurons also showed adapting and non-adapting responses. INs were classified as a single-burst-accommodating (bAC) e-type. 27 Morpho-electrical models A multi-objective optimization pipeline was used to build electrical models (e-models) for the e-types in Figure 1B, resulting in five e-models (see STAR Methods). These e-models were combined with 92,970 morphologies to generate 142,678 unique morpho-electrical models, which were assessed with a battery of tests to reject models with electrical features and firing behavior that differed significantly from the experiments.

Reconstructing thalamoreticular microcircuit structure
We defined the thalamoreticular microcircuit and determined its neuron numbers, composition, and positions ( Figures 1C-1G). The microcircuit spans parts of the Rt and VPL of the thalamus. We chose the VPL nucleus because it receives information from the hindlimb 28 and is representative of a more generic non-barreloid thalamus. The extent of the Rt_RC dendrites informed the horizontal dimension of the microcircuit (see STAR Methods), while the height of the circuit comes from the anatomical parcellation in the Allen Brain Atlas. 29 The resulting volume is $0.22 mm 4 with a length of 323 mm and a height of 800 mm. Using semi-auto-mated cell counting (see STAR Methods), we found an average cell density of 68,750 ± 1,976 cells/mm 3 in Rt and 57,467 ± 5,201 cells/mm 3 in VPL. 100% of Rt neurons are inhibitory, while 0.5% of VPL neurons are inhibitory (mean ± std, n = 37 slices).
The resulting model microcircuit contains a total of $14,000 neurons (averaged over seven microcircuit instances) composed of 4,909 ± 283 Rt_RCs, 8,952 ± 517 VPL_TCs, and 47 ± 2 VPL_INs. The Rt has 100% m-type Rt_RC neurons with 57% cAD_ltb and 43% cNAD_ltb e-types. VPL has VPL_IN m-type with bAC e-type, and VPL_TC m-type with 64% cAD_ltb and 36% cNAD_ltb e-types. Soma positions were determined using a pseudorandom algorithm and taking into account morphological constraints. 4,[28][29][30][31][32] Reconstructing and validating synaptic connectivity The connectivity between neurons in a microcircuit was established using an adapted algorithm based on Markram et al. 24 and Reimann et al. 33 A linear relation between the dendritic surface in the thalamus and bouton numbers on reticular axons 34 was found, suggesting predictability of functional synapse locations. 35 We used neuron morphologies and bouton densities as constraints. Our experimental dataset showed that TCs had 0.102 ± 0.021 boutons/mm (n = 9 axons) and Rt neurons 0.124 ± 0.002 boutons/mm (n = 2 axons). Synaptic connections were established by presynaptic axons and postsynaptic dendrites/ somata ( Figure 2). INs largely established connections through presynaptic dendrites, as observed in the visual thalamus. 36,37 Medial lemniscal (ML) and corticothalamic (CT) synapses were included using volumetric bouton densities as constraints (Figure 2D). Data on lemniscal synapses in the mouse VPM 38 was used and the relative proportions of CT to lemniscal synapses onto TCs and the ratio of CT to TCs synapses in the Rt. 39,40 The model resulted in 4.77 million intrathalamic synapses and 17,998 lemniscal and 40,905 corticothalamic synapses. We compared the synapse convergence onto reticular neurons 41 and the distribution of number of synapses per connection, 36 and both gave results comparable with the experimental counterpart ( Figures 2E and 2F).
A detailed map of synaptic physiology in a thalamic microcircuit We modeled the physiology of synapses in the thalamic microcircuit using data from experiments and literature on shortterm plasticity, postsynaptic potential amplitudes, time constant  Resource ll OPEN ACCESS of synaptic currents, and reversal potentials. We identified three types of short-term plasticity-inhibitory depressing (I2), excitatory depressing (E2), and excitatory facilitating (E1)-and used synapse models featuring stochastic transmission and short-term plasticity. 24,42 We constrained the parameters of the Tsodyks-Markram model with available thalamic data (See STAR Methods; Table S1; Figures 2G and 2H). The model provides a comprehensive map of synapse types with the main external afferents ( Figure 2H). The model was able to reproduce the stochastic nature of synaptic release and in silico pairedrecording data, closely matching the experimental results (Tables S2 and S3; Figures 2I-L).
Gap junction connectivity between reticular neurons is predicted by dendrodendritic appositions Rt neurons are connected through electrical synapses, contributing to synchronization in thalamic networks. [43][44][45][46] Gap junction (GJ)-coupled neurons are estimated to be 30%-50% of the total population in Rt 47 and play a role in synchronization and desynchronization balance in thalamic networks. 11 In the model, we defined GJ locations by identifying dendritesoma appositions and found that assigning GJs to 30% of appositions matched experimental results ( Figure 3A). The model predicts that each neuron is directly coupled with 2-20 other neurons at distances of 40-120 mm, with rare coupling at 300-400 mm ( Figure 3B), in line with experiment findings. Coupling strength was validated with a GJ conductance of 0.2 nS, 48 and calculated coupling coefficients between simulated paired recordings of 0.023 ± 0.008 (mean ± standard deviation, Figure 3C) were within the variability of mouse Rt neurons. 44 These results indicate that many aspects of GJ connectivity, in particular its distance dependence, can be predicted by the morphological properties of Rt_RC dendrites.
Spontaneous and evoked activity during wakefulness and sleep in the model thalamoreticular microcircuit Spontaneous and evoked sensory activity (wakefulnesslike) We first explored spontaneous and evoked activity in the resulting microcircuit in a simulated in vivo wakefulness-like condition ( Figure 4). In this condition, the model reproduces the distribu-tion of firing rates in first-order thalamic and Rt neurons during quiet wakefulness, characterized by low firing rates (<10 Hz). 49-51 The simulation exhibits irregular firing activity in all m-types, with higher firing rates in VPL_INs due to their lower spiking threshold in the model ( Figure 4B). Single-cell activities are dominated by single tonic spikes ( Figure 4B), consistent with their predominance over low-threshold bursts in wakefulness-like states. 52 We then compared spontaneous activity with sensory-evoked responses with brief activation of a subset of ML fibers, simulating a whisker flick or a brief stimulation of the hindlimb. 53 We found increased population firing rates in VPL_TC neurons, with a peak in firing in a short time window (5-10 ms) following the stimulus (Figures 4A and 4B; Video S1). This activity is followed by a period of silence in some cells, lasting 100-200 ms ( Figure 4B). Rt_RCs exhibit an increase in firing rate compared with their spontaneous activity as a result of the excitation they receive from VPL_TCs ( Figure 4B). The increased activity in Rt_RCs lasts for 50-100 ms after the stimulus and gradually decreases to baseline levels, resulting in a clear hyperpolarization of VPL_TCs ( Figure 4A). The longer activation of Rt_RCs is in line with experimental findings and suggests an important role of the Rt in limiting sensory responses and focusing them on rapid perturbations. 50 Simulating the activation of an increasing number of afferent ML fibers revealed increased recruitment of reticular inhibition ( Figure S2A) and decreased response latency and variability with increasing stimulus size, increasing synchronous responses ( Figure S2B). As expected, increased recruitment of inhibitory cells results in increased hyperpolarization of the VPL ( Figures S2A and S2C).
At the topographical level ( Figure S2D; Video S2), VPL cells located at the center responded to the stimulus, while the activation in the Rt was broader, with a degree of VPL-mediated depolarization in the periphery of the Rt. This result suggests that receptive fields in the Rt are larger than in the VPL, yet they provide topographically aligned inhibition to the thalamus in a focal manner, as shown for visual reticular neurons. 49,54 This result demonstrates that Rt can contribute to precise responses in the VPL not only by rapidly inhibiting directly responding neurons but also by limiting the response in the surrounding area.  Table S2). (J) Comparison of PSP amplitudes in the model for seven characterized pathways in house or in the literature (see Table S3).  [55][56][57] To examine the model response in such conditions, we simulated a transition from a wakefulness-like to nonrapid eye movement (NREM)-like state. We simulated the up and down states that characterize NREM by periodically removing the firing background from the cortical afferents for 500 ms and reactivating it for another 500 ms to the same level as the one used for our standard in vivo wakefulness-like condition ( Figures 4C and 4D; Video S3).
During the up state, CT input drives Rt activity and initiates spindle-like oscillations through post-inhibitory rebound in the VPL ( Figure 4C). Interestingly, during the down state, Rt_RCs become highly hyperpolarized (<₋70 mV), while VPL_TCs are less affected ( Figure 4C). During the simulated cortical up state, Rt_RCs fire robustly, predominantly low-threshold bursts, causing deep inhibitory postsynaptic potentials (IPSPs) in the VPL_TCs that, in turn, often respond with post-inhibitory rebound bursts, the hallmark of spindle-like activity in the thalamus ( Figure 4D).

Thalamic responses to sensory input exhibit adaptation and cortical enhancement
Numerous studies have shown adaptation to sensory stimuli in the lemniscal pathway due to short-term depression of lemniscal excitatory postsynaptic potentials (EPSPs). [58][59][60][61][62] To explore this, we set out to use the model to recreate a recent study that showed that depressed responses in the somatosensory thalamus can be enhanced by cortical activation in anesthetized mice. 63 Consistent with experimental findings, activating lemniscal fibers, in the model, with trains of stimuli at 8 Hz results in high response probability to the first stimuli, while subsequent EPSPs exhibit decreased amplitudes ( Figure 5A). The activation of cortical afferents 200 ms before the sensory stimulus with noisy input at a mean firing rate of 4 Hz increased the firing probability to all the other stimuli in the train, thus counterbalancing the adaptation ( Figure 5B), as seen in the original study. 63 Cortical enhancement of sensory responses is frequency dependent and frequency selective The effects of the cortex on thalamic activity depend on the frequency and pattern of cortical activation. 64,65 In this study, when sensory stimuli of elevated frequencies (e.g., 10-20 Hz) were presented, thalamic responses were highly depressed and the probability of firing was greatly reduced ($0.1) (Figures 5C and  5D). When the mean firing rate of the cortex was increased, the efficacy of corticothalamic (CT) inputs in counterbalancing sensory adaptation also increased, as indicated by an increase in the firing probability of thalamic cells ( Figure 5D). The cortical enhancement of sensory responses was found to be frequency selective, with the greatest enhancement observed when sensory stimuli were at a frequency of around 10 Hz. This was due to the recruitment of Rt inhibition and the activation of intrinsic subthreshold currents in thalamic cells ( Figures 5E, S3F, and S3G). At higher frequencies of sensory stimuli, the effects of the cortex on thalamic activity were less pronounced due to synaptic short-term depression and limited activation of the low-threshold calcium current in thalamic cells ( Figure 5E). The selective enhancement of 10-Hz sensory stimuli was most prominent when the mean CT firing rate was R10 Hz.

Spindle-like oscillations emerge through thalamoreticular interactions
The model exhibited spindle-like activity during simulated cortical up states. It is known that sleep spindles are generated through interactions between cells in the Rt and thalamic cells. 10,66 Numerous studies have explored the mechanisms underlying the generation of spindle-like oscillations. 6,10,[67][68][69][70][71] This process involves the activation of Rt neurons, which hyperpolarize TC cells via inhibitory connections between the two cell types. 66,72 This hyperpolarization primes TC cells for rebound bursts, which perpetuate the spindle cycle through a "pingpong" interaction between Rt and TC cells. 73 The model was not specifically designed to replicate these known mechanisms, but it was found to do so. Activating the Rt increases thalamic bursting and initiates spindle-like oscillations A subset of cells in the Rt were briefly activated with a current pulse (20 ms) to simulate optogenetic activation 74 ( Figure S4). This caused brief ($250 ms) oscillatory responses in both Rt and VPL cells ( Figure S4A). The local field potential in the VPL showed oscillations at $10 Hz ( Figure S4B), which is consistent with the spindle frequency range in vivo. Single-cell responses Resource ll OPEN ACCESS also showed increased bursting in both Rt and VPL cells (Figure S4D). This synchronized burst firing in Rt cells was found to be a potent trigger of spindle oscillations ( Figures S4D-S4G). These results confirm the important role of intrinsic neuronal mechanisms in the generation of spindle-like oscillations 6,75-79 and validate that optogenetic-like activation of the Rt in the model can trigger them, as shown by Halassa et al. 74 Spindle-like oscillations are maintained by ping-pong interactions between the Rt and thalamus To investigate the mechanisms of spindle-like oscillations, the network of thalamoreticular cells was simulated in a silent, hyperpolarized state (only receiving spontaneous synaptic release). When a group of neighboring Rt_RCs was activated (20-ms pulse; see Video S4), the duration of the resulting pingpong interactions between VPL_TCs and Rt_RCs that mediate the spindle oscillation was influenced by the synaptic release probability (P rel ) ( Figure 6A). Short-term synaptic depression, controlled by the recovery time constant (t RecDep ), also affected the duration of the oscillations ( Figure 6B). The frequency of the oscillations ($5-6 Hz; see Figure 6C) was similar to the barrages of IPSPs recorded in spindle waves in ferrets in vitro. 6 The study found that synaptic mechanisms, including P rel and short-term synaptic depression, play a role in the duration of spindle-like oscillations, which was not previously known. GJs increase the duration of spindle-like oscillations by propagating low-threshold bursts across the reticular network GJs between reticular neurons contribute to the duration of spindle-like oscillations by synchronizing activity and transmitting low-threshold bursts, as shown in previous studies. 8,11,44,46 Our results show direct evidence of GJs' effect, as removal of GJs prevented spindle-like oscillations from occurring. Comparison of spiking activity and spatiotemporal patterns in the Rt and VPL revealed increased spiking and longer-lasting excitation in Resource ll OPEN ACCESS the Rt and higher inhibition in the VPL with the presence of GJs ( Figures 6D-6I). GJs transmitted slow signals efficiently, with decreased amplitude in peripheral cells ( Figure 6H; see Video S5 for comparison of spatiotemporal activity).
The termination of spindle-like oscillations is determined by short-term depression and buildup of intra-reticular inhibition Several intrathalamic mechanisms have been identified as contributing to the termination of spindle activity. These include the Ca 2+ -dependent upregulation of I h current in thalamic cells, 68,80 hyperpolarization of reticular neurons, with the activation of Na + -and Ca ++ -dependent K + currents. 81 External inputs, such as desynchronized cortical activity and noradrenergic input from the locus coeruleus, have also been hypothesized to play a role in spindle termination. 82,83 The model revealed additional cellular and synaptic mechanisms responsible for spindle-like oscillation termination, not previously investigated. Intrathalamic connections are governed by depressing synapses (Figure 2), which leads to decreased efficacy between Rt and VPL neurons. By scaling the time constant of recovery from short-term synaptic depression, it was found that reducing it to 85% of the control value in Rt_RC to VPL_TC connections caused the oscillation to not terminate (Figure 6B). Short-term depression in the VPL_TC to Rt_TC pathway had limited effect on termination and could even be removed if present in Rt_RCs to VPL_TCs. Rt_RC bursting activity is important for spindle-like oscillation initiation and duration, but these cells inhibit each other through GABAergic synapses. Inhibitory Rt-Rt synapses were studied and found to be crucial, as reducing their conductance led to non-terminating oscillations ( Figure 6C) with an increase in Rt_RC cells participating in each cycle. Waxing and waning of spindle-like oscillations emerges due to intrinsic cellular and synaptic mechanisms Waxing and waning are defining characteristics of spontaneous sleep spindles recorded at the cortical level in the electroencephalogram (EEG), local field potential, and in thalamic recordings. 6 Despite this characteristic pattern of activity, the mechanisms underlying the spindle-shaped oscillation have not been clear. 84 It has been thought that TC and corticothalamic interactions were responsible for waxing and waning. 83,[85][86][87] In the model, we found that the waxing and waning of spindles can be generated within the thalamoreticular microcircuit alone ( Figure 6L). Specifically, the waxing of a spindle oscillation is created by the rhythmic recruitment of neurons, first in the Rt, and second through the ping-pong interaction with the thalamus: each ping from the Rt successively recruits additional neurons and generate a stronger pong, via low-threshold bursts in VPL_TC cells. GJs augment this process through their ability to rapidly recruit Rt neurons. The waning of the spindle-like oscillation is a result of the progressive reduction in the probability of synaptic release (due to short-term synaptic depression), the subsequent decrease in postsynaptic potential (PSP) amplitudes, and the consequential reduction in recruitment of neuron firing during the ping-pong interaction between Rt and VPL neurons. At the same time, mutual inhibition between Rt_RCs progressively builds up, as additional Rt_RCs are recruited, and prevents the spread of the activity, acting as a self-limiting mechanism.
Depolarizing the Rt decreases spindle-like oscillation duration Spindle oscillations in naturally sleeping (i.e., non-anesthetized) rodents are more easily evoked when thalamic activity is mildly synchronized and thalamic neurons are less active. 74,88 Their features, such as frequency and duration, evolve during NREM episodes. 88,89 In the model, we show that spindle duration can vary as a result of the membrane potential in the thalamus and the Rt.
We studied how membrane potential levels in Rt and VPL affect spindle-like oscillation duration, frequency, and peak firing ( Figure 7). This was done by injecting noisy current into Rt_RCs and VPL_TCs to approximate the influence of neuromodulators on thalamic and reticular activity. A range of depolarization levels was explored, from hyperpolarized to near firing threshold.
Stimulating a subset of Rt_RCs with a 20-ms current pulse resulted in different spiking responses based on the depolarization level. When Rt_RCs were depolarized, spindle-like oscillation duration decreased ( Figure 7A) and the number of spikes per burst decreased ( Figure 7B), leading to reduced inhibition to VPL and decreased VPL_TC firing. As small as a 2-to 3-mV increase in membrane potential was sufficient to observe this effect. The waxing and waning in firing responses also tended to become predominantly waning when the Rt was depolarized ( Figure 7A). Depolarizing the thalamus increases spindle-like oscillation duration and frequency When the membrane potential in the Rt is held constant (around ₋76 mV) and the VPL is depolarized, the oscillation increases in (E) A schematic explaining why cortical enhancement is greater for sensory stimuli at around 10 Hz compared with higher frequencies (around 20 Hz). Sensory stimuli around 10 Hz are timed with post-inhibitory rebounds and activation of low-threshold calcium and produce larger EPSPs that can reach firing threshold with cortical activation. For higher stimulus frequencies, EPSPs decrease in amplitude due to synaptic depression, and cortical inputs are no longer sufficient to reach the firing threshold and counterbalance the adaptation. Resource ll OPEN ACCESS duration and frequency ( Figure 7C). VPL depolarization results in deeper and faster IPSPs ( Figure 7D). Deeper IPSPs result in stronger post-inhibitory rebound responses in VPL_TCs, which in turn excite more Rt_RCs causing a longer period of pingpong interactions between the two populations, increasing the oscillation duration. Faster repolarization after the IPSPs is associated with post-inhibitory rebound responses occurring at shorter intervals, driving the ping-pong of activity at higher frequencies.
Differential depolarization of the Rt and thalamus modulates spindle properties For each combination of membrane potentials in Rt_RCs and VPL_TCs, we calculated oscillation duration, frequency, and peak firing ( Figures 7E-7G). The oscillation duration map shows that clear spindle-like oscillations (duration R 500 ms) can only be evoked in a region where Rt_RC membrane potentials are below ₋75 mV and VPL_TCs are more depolarized than Rt_RCs. If we assume that VPL_TCs neurons are in general more depolarized during wakefulness than during sleep, and that Rt_RC neurons are more hyperpolarized, this result suggests that spindle-like oscillations are easier to evoke at the transition between wakefulness-like and sleep-like states, consistent with experimental observations. 8,90,91 When both VPL_TCs and Rt_RCs are at their baseline membrane potential, the frequency decreases to 5-6 Hz.
Simulating the effects of neuromodulation on thalamic and reticular neurons causes spindle-like oscillations to cease Sleep spindles are a defining characteristic of stage 2 NREM sleep, and appear less frequently in deeper stages of NREM sleep. 8,90,91 The change in the incidence of spindles are due, at least in part, to neuromodulatory changes. 92-94 Our previous simulations, in the in vitro-like condition, showed that membrane potentials have a strong effect on oscillation duration and frequency ( Figure S5). They also indicated that, when the VPL is hyperpolarized and the Rt is depolarized, spindles are less easy to evoke. We hypothesized that such differential polarization levels of the Rt and VPL would resemble the transition to deeper stages of NREM sleep and result in decreased incidence of spindles. We tested this hypothesis with the model in the in vivo-like condition using simulated cortical up and down states (as in Figure 4). We then progressively depolarized the Rt and hyperpolarized the VPL, approximating the differential effect of neuromodulators on thalamic and reticular populations. 95-97 These simulated neuro-modulatory changes cause the spindle-like oscillations to occur with reduced amplitude and then to cease, while the corticallygenerated up-down states continue ( Figure S5).

DISCUSSION
In this study, we developed a detailed model of the thalamoreticular microcircuit, a network of brain cells involved in sensory processing and rhythm generation. The model incorporates experimental data on the anatomy, physiology, and connectivity of individual neurons, as well as the three-dimensional (3D) organization of the reticular and VPL nuclei. It was validated across different simulated conditions, including wakefulness and sleep. The model showed that short-term synaptic plasticity and mutual inhibition within the Rt play key roles in the termination of spindlelike oscillations. It also found that GJs, corticothalamic feedback, and membrane potentials play a key role in modulating the duration and frequency of spindle-like oscillations. The model provides a comprehensive account of thalamoreticular network dynamics and can be used to interpret alterations in corticothalamic feedback and rhythm generation in health and disease. The model and accompanying data are openly available for future research (see section ''data and code availability'' below).

Comparison with prior work
The model differs from previous ones in several aspects, including scale (in terms of number of neurons), the level of biological detail, and scope. 68,69,83,98-102 While previous models were explicitly built and tuned to study network phenomena, such as sleep spindles, oscillations emerge in our model only by fitting cellular and synaptic properties, without tuning for network behavior. This data-driven approach enables novel insights on spindle generation and termination (see below). Prior models have largely used single-compartment neurons or simplified multicompartmental models, 103 while this model uses reconstructed 3D neuron morphologies to constrain the biophysical models and connectivity.

Morphological constraints on thalamoreticular connectivity
The model suggests that many properties of intrathalamic connectivity can be predicted by the axodendritic and dendrodendritic overlap of neuron morphologies. This includes the distribution of synapse locations and the number of synapses per connection Resource ll OPEN ACCESS established by thalamic INs, as well as the convergence of synapses from TC axon collaterals, other reticular neurons, and corticothalamic afferents onto reticular neurons. The model accurately recreates the distribution of single and multi-synaptic connections, 36 and it recapitulates distance-dependent GJ connectivity between reticular neurons, as observed in experimental studies. 45 The model also shows that dendrodendritic appositions provide a sufficient basis for determining GJ locations, and that the distancedependent distribution of GJs between reticular neurons is determined by the extent of their dendrodendritic overlap.

Common cellular and synaptic mechanisms between sensory enhancement and spindle generation
The model recreates the network dynamics of spontaneous and sensory-evoked activity in wakefulness-like states. It allows for simultaneous observation of direct and indirect sensory responses in the thalamus and Rt. Responses to stimuli are clearly visible in the Rt, as has been shown experimentally with sensory, auditory, or visual stimuli. 50,53, 104 The model also shows that robust inhibition from the Rt generates strong surround inhibition in the VPL, and that cortical activation sharpens the spatial prop- erties of sensory inputs by evoking stronger Rt-mediated surround inhibitions. This is in line with recent research on the visual thalamus. 49 As previously reported in anesthetized animals, the responses to trains of sensory stimuli were adapting for relatively low frequencies of 4-5 Hz, and the degree of adaptation increased with the increase of stimulus frequency. 59,60,62 The adaptation was reduced by simulated cortical activation, approximating cortical activity during activated states, allowing the transmission of high-frequency inputs. 63 The model showed that the thalamoreticular microcircuit exhibits frequencyselective cortical enhancement. When simulating different frequencies of input and corticothalamic activation, peak enhancement of thalamic responses was observed for lemniscal inputs around 10 Hz. This enhancement occurs because cortical activation recruits sufficient Rt inhibition to activate I h and low-threshold Ca 2+ currents, enhancing the gain of input to TC cells at $10 Hz. This is notable because rhythmic activity around this frequency emerges in TC networks during different behavioral states such as sleep spindles and alpha oscillations during attention. 105 Our simulations suggest that $10-Hz rhythms are intrinsic to thalamoreticular networks and could be responsible for enhanced gain of sensory inputs and TC activity around that frequency. These results are consistent with the observation of overlapping thalamic mechanisms between sensory processing, attention, and sleep. 105 A novel cellular and synaptic account of spindle generation The model provides novel insights into the mechanisms underlying the generation of sleep spindles in thalamic networks, including the role of TC-Rt synapse efficacy, short-term synaptic Cell Reports 42, 112200, March 28, 2023 13 Resource ll OPEN ACCESS depression, mutual inhibition between Rt neurons and GJs in the thalamic generation of waxing-and-waning spindle-like oscillations.
The model generates spindle-like oscillations without being explicitly designed for this purpose. The model is based on in vitro findings, but many aspects closely resemble thalamic activities during spindle oscillations in rodents in vivo. 88,89,106 This indicates that the model is robust and can generate spindle-like oscillations within experimentally plausible parameter ranges.
Consistent with previous experimental and modeling studies, spindle oscillations in the model are generated through a combination of intrinsic mechanisms, namely low-threshold calcium bursting in reticular neurons 75,77 and the synaptic interactions between reticular and TC neurons. 67,68,107 Dendrodendritic GJs between reticular neurons have been known to transmit low-threshold bursts between cells, promote spiking correlations and synchronized activity. 44,46 The model shows that GJs also influence spindle duration. Furthermore, the spatial organization of Rt dendrites and their connections through GJs also shows a clear functional role in propagating the stimuli along the horizontal dimension of the microcircuit and enhancing spindle-like oscillations.
The model showed that spindle termination can occur through synaptic mechanisms alone, without requiring specific ionic mechanisms such as Ca 2+ -dependent I h current upregulation in TC cells or cortical input desynchronization. 68,83,100,102 The model revealed that mechanisms such as short-term synaptic depression and mutual inhibition between reticular neurons limit the duration of the oscillation. While mutual inhibition has been previously suggested as a factor, 11,13,14,70,108 short-term synaptic depression has not been considered as a termination mechanism. Based on the role of Rt-Rt mutual inhibition, which decreases spindle duration, and GJs, which tend to prolong it, we propose that Rt interactions have a self-limiting contribution to spindle oscillations. GJ coupling promotes synchronized IPSPs and TC excitation, recruiting more Rt neurons. When a critical recruitment level is reached, the overall excitation is countered by reciprocal inhibition and short-term synaptic depression, and the oscillatory activity limits itself.
Thalamoreticular circuitry is responsible for the waxing and waning of spindles The gradual increase in amplitude (waxing) and the gradual decline (waning) of spindle oscillations has been hypothesized to be due to changes in cortical activity impinging on the Rt and thalamus. 83,85-87 However, the model demonstrates that cellular and synaptic mechanisms of thalamoreticular circuitry that underlie the increased recruitment of additional neurons during each cycle of a spindle can produce the waxing phenomena, while synaptic depression and the buildup of inhibition within the Rt are sufficient to explain the waning phenomena.
Relative differences in the excitability of thalamic and reticular nuclei can explain differences in spindle frequency, density, and amplitude We found that changes in depolarization level of thalamic and reticular neurons affect oscillation frequency, duration, strength, and incidence of spindle-like oscillations. Our model approximates the dynamic change of thalamic network states during transitions between wakefulness and sleep, based on known mechanisms of action of neuromodulators such as acetylcholine (ACh) and noradrenaline (NA). ACh levels decrease during the transition from wakefulness to NREM sleep, causing hyperpolarization of the Rt and depolarization of the thalamus. 95-97 NA levels fluctuate during NREM sleep in vivo, leading to depolarized membrane potentials that affect low-threshold bursting activity underlying sleep spindles. 93 Other alterations in thalamic and reticular excitability could be due to changes in synaptic activity and plasticity of corticothalamic and corticoreticular projections. 109 These findings have important implications when considering changes in sleep spindle frequency, density, and amplitude as potential biomarkers of disease. The model could be used to study the mechanisms that affect spindle oscillations and their properties during sleep in conditions such as schizophrenia, 16-19,21 ADHD, 22 and Alzheimer's disease. 23 For example, spindle densities have been found to decrease in patients with Parkinson's disease, and their density and duration are sensitive to Alzheimer's disease. 23, [110][111][112] Previous research has suggested that the cortex is required for the generation of slow spindles, 113 but the present model shows that the thalamus can modulate the frequency of spindle-like oscillations, indicating that thalamic mechanisms may be sufficient to generate both types of spindles.

Limitations of the study
The model presented here is a first-draft reconstruction of thalamoreticular microcircuitry and is not intended to be a complete representation of the full circuitry. The model does not account for reciprocal connectivity with the cortex and therefore is not sufficient to explore the role of the thalamus in higher-order functions that depend on the cortex, such as attention and cognition. We are making the model openly available to facilitate future extensions and refinements.
This model represents a primary somatosensory microcircuit and includes only the "core" cell characteristic of sensory thalamic nuclei that target cortical layers III and IV, rather than the ''matrix'' cells, found in higher-order thalamic nuclei that target supragranular layers. Different thalamic nuclei have different mixtures of excitatory cells that have unique projection characteristics to the cortex. 114,115 Future refinements could take into account these different cell populations.
In addition, the cellular e-types in the present model do not capture the full diversity of known ion channels implicated in bursting behavior in reticular neurons or dendritic properties of thalamic INs. Inclusion of additional and more specific ion channel mechanisms, (e.g., variants of low-threshold Ca 2+ Cav3.1, Cav3.2, and Cav3.3) in TC and Rt neurons could more accurately reproduce differences in bursting behavior in thalamic and reticular neurons, as well as dendritic properties of thalamic INs. 75,77,80,116,117 Further refinements of the model, such as including the distribution of different morphological types within the thalamic and reticular domains, could improve its accuracy. [118][119][120][121] Additionally, incorporating the laminar structure 120 Resource ll OPEN ACCESS could provide a more comprehensive understanding of the role of TC interactions in sensory processing, attention, and rhythm generation across wakefulness and sleep.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:     includes code for the ion channel models, single neuron models, synapse models, circuit files in SONATA format, 1 simulation output. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Experimental data were collected in conformity with the Swiss Welfare Act and the Swiss National Institutional Guidelines on Animal Experimentation for the ethical use of animals. The Swiss Cantonal Veterinary Office approved the project following an ethical review by the State Committee for Animal Experimentation.
In vivo-stained TC and Rt morphologies were obtained from C57Bl/6J 3 month-old male mice.
In vivo labeling of reticular neurons were employed in SSt-Cre; Ai139 P54 mice. Electrophysiological data for TC and Rt neurons and interneurons (INs) were characterized in vitro from brain slices of P12-35 GAD67-eGFP or C57Bl/6J mice.
For all mouse lines, if not otherwise specified, both male and female subjects were used.

METHOD DETAILS
Constraining and validating the model with experimental data The model microcircuit was built by constraining and validating it at multiple levels, using experimental data and algorithmic approaches, based on methods published previously. 24 For validation we performed direct comparison of the model properties with experimental measurements that were not used during the model building steps. Before describing the details of the reconstruction, validation and simulations, we provide a list of data used for constraining the model, the validation data and further validations at the network level.
Experimental data used to constrain the model The following experimental data was used to constrain the model, further details on the experimental procedures and literature references are provided below. In vitro-stained neurons were mainly located in primary somatosensory nuclei (VPL, and ventral posteromedial nucleus -VPM) and the somatosensory sector of the reticular nucleus. 34,126,127 Although the main focus of our model was the VPL, we found that specifically targeting this nucleus was challenging. To increase the available experimental data we thus included recordings from the VPM. Reconstructions used the Neurolucida system (MicroBrightField) and were corrected for shrinkage along the thickness of the slice. Shrinkage along other dimensions was taken into account during the unraveling step (see below). Dendrites were reconstructed with a 1003 magnification (oil immersion objective) and axons at 603 (water immersion objective).
In vivo-stained TC and Rt morphologies were obtained through different experimental techniques. In some cases, neurons were labeled by injection of replication-defective Sindbis virus particles in the thalamus or Rt nucleus in C57Bl/6J adult mice 128 or electroporation of RNA of the same virus. 129 The virus labeled the membrane of the neurons thanks to a palmitoylation signal linked to a green fluorescent protein (GFP). Brains were cut in 50 mm serial sections and immunostained against GFP and enhanced with glucose oxidase-nickel staining. 130 Neurons were reconstructed from sequentially-ordered slices under bright-field optics using the Neurolucida system (MicroBrightField). The complete method is described elsewhere. 131 In vivo-labelled TC morphologies were obtained from the Janelia Mouselight project, from sparsely-labeled adult C57/BL6 mice brains; the method is described in detail elsewhere 25 and summarized here. Brains were then delipidated, fluorescence was enhanced by immunolabeling and imaged with a 403 oil-immersion objective. This procedure generated large datasets of high-resolution image stacks. The 3D reconstructions were conducted combining semi-automated segmentation of the neurites, human annotation and quality control. Janelia Mouselight reconstructions lacked diameter variations in their neurites, which is important for accurate electrical modeling of neurons. 132 For this reason, we only used their axons in order to increase the variability of our axonal reconstructions. We obtained 96 morphologies whose soma was located in the thalamus and we visually inspected their shape along with 3D meshes of the reticular nucleus of the thalamus using the Janelia MouseLight Project (RRID:SCR_016668; http://ml-neuronbrowser.janelia.org/). Since most thalamocortical neurons project to the Rt on their path to the cortex, 114,127 we selected the 41 morphologies which gave off collaterals in the reticular nucleus. We assumed that neurons without collaterals in the Rt were partially labeled and/or reconstructed, since those collaterals are often very thin. 133 Given the limited number of reconstructed morphologies of neurons in VPL and VPM in the Janelia MouseLight dataset, we included 27 axons (with collaterals to Rt) from other thalamic nuclei. To ensure that the connectivity would not be impacted, we analyzed the geometrical properties of the Rt collaterals and found that the difference within the same nucleus was as high as the difference between nuclei.
For in vivo labeling of reticular neurons virus injections for sparse labeling of whole brain neuron morphologies were employed in SSt-Cre; Ai139 adult mice. 134,135 Brains were imaged using fluorescence micro-optical sectioning tomography (fMOST). 136 Neurons were manually reconstructed from high resolution image stacks obtained after slicing. Further details of the method are available in related publications. 137 Morphology analysis, alignment and visualization Raw morphological data did not have a common orientation along a principal axis, which is necessary to place them in the microcircuit volume according to biologically-plausible constraints (see below). We thus computed a rotation matrix so that the principal axis of the morphology was parallel to the vertical axis of the microcircuit. The principal axis of TC morphologies was the one 24 Cell Reports 42, 112200, March 28, 2023 Resource ll OPEN ACCESS connecting the center of the soma and the center of mass of the axon collaterals in the Rt nucleus (see below). For Rt neurons, the principal axis connected the soma and the center of mass of the axonal arborization in the thalamus. After rotating the morphologies, we visually validated the results. Rotation of the INs was not performed, since no orientation information relative to known landmarks was available.
For morphology analysis we used the open-source library NeuroM (https://github.com/BlueBrain/NeuroM). To identify the TC axon collaterals projecting to the Rt we selected the morphological sections which had branch order R 1 and path distance from the soma <2,500 mm and visually validated the results. For some morphologies, we selected those having path distance % 2,000 mm, because some TC neurons have collaterals projecting to other subcortical regions (e.g striatum), see Clascá et al. 114 Raw morphological data were algorithmically corrected for slicing artifacts and processed to generate a large pool of unique morphologies for building the microcircuit and connectivity. Spurious sections, which were accidentally introduced during manual reconstruction, were identified as those having 0 mm diameter and removed. The details are described in Supplemental Experimental Procedure of the neocortical microcircuit model, and summarized below. 24 The morphology images in Figure 2 were created using NeuroMorphoVis. 138 Unraveling morphologies Since we found that 3D reconstructions from in vitro-stained neurons had increased tortuosity in their dendrites as a result of tissue shrinkage, we unraveled them using an existing algorithm. 24 This process resulted in an increase of the reach of the morphologies, while preserving the original length of the branches. Briefly, unraveling was performed by sections and for each section a sliding window composed of a given number of successive points was created. The number of points in the sliding window (N) was the only parameter of the algorithm and we found that N = 5 previously used performed well on thalamic morphologies. The general direction of the points in the window was computed using principal component analysis (PCA). The segment at the middle of the window was then aligned along this direction. It meant that its direction was set to the one of the sliding windows but it retained its original length. The sliding window was moved over all points of the section and the algorithm was applied to all sections.

Repairing morphologies
Most of the in vitro-stained morphologies were truncated at slice edges and in the case of some TC morphologies, which have very dense dendritic arborization, this resulted in a significant decrease in dendritic mass. We applied an existing algorithm 24,139 to repair missing dendritic branches. First, the algorithm detects cut points on the XY plane, i.e. the plane parallel to the slice, along the Z direction (parallel to the slice thickness). The 3D coordinate system was centered on the morphology soma. Although the algorithm was designed to detect cut points on two planes, we found that our morphologies were truncated on the top plane. We improved the algorithm by searching the cut points before unraveling the morphologies and updated their position during the unraveling step. Cut detection required a tolerance parameter to detect terminal points within a certain distance from maximum Z extents. We found that 15 mm gave the most accurate results by visual inspection of the morphology. Some terminal points were then tagged cut points and dendrites were repaired. The dendrite repair process created new dendritic sections starting at the identified cut points. Dendrite repair did not aim to recover the initial morphology, but rather recreated it in a statistical manner, under the assumption of statistical symmetry of the morphology. This method analyzed the behavior of intact branches as a function of branch order and euclidean distance from the soma. For each branch order, probability density clouds of branch continuation, bifurcation or termination were calculated in a series of concentric spheres. 140 At each cut point, the behavior of the branch was sampled according to the calculated probabilities. The factor governing the direction of the re-grown branches was adjusted to achieve final branches tortuosity comparable with our experimental data. To address neurite swelling artifacts at cut points, the diameters of the re-grown branches were set to the average diameter of the last section.

Morphology diversification
We increased the variability of the reconstructed and repaired morphologies to ensure robust and invariant connectivity patterns. 35, 141 We followed a previously published method 24 to generate a unique branching pattern for each morphology, while maintaining the general morphological and electrical structure for each m-type. In summary, branch lengths and rotations at each bifurcation point were varied according to random numbers drawn from Gaussian distributions with mean 0% and standard deviation 20% for branch lengths and mean 0 and standard deviation 20 for branch rotations. A sample of the resulting morphologies was visually validated, and we did not find significant alterations of their structure for any of the m-types ( Figure S1).
We then applied a mix-and-match procedure to maximize the utilization of good morphological reconstruction data. This procedure divided dendrites from axons and allowed us to combine good dendritic reconstructions of TC and Rt dendrites from in vitro and in vivo-stained neurons and good axonal reconstructions from in vivo-stained neurons. In vitro-stained neurons typically lacked reconstruction of the full axon due to the slicing procedure and/or poor labeling. For each morphology, we manually annotated which dendrites and axons were to be kept. The decision in most cases depended on the labeling method (in vitro vs. in vivo). To increase the probability that in vivo-stained morphologies and in particular the axons of TC and Rt morphologies were compatible with the microcircuit dimensions (see below) we duplicated and scaled the morphologies along their principal axis (Y axis) by ± 2.5%.
Reconstructing the electrical diversity of neurons Electrophysiological data The firing patterns of TC, Rt neurons and interneurons (INs) were characterized in vitro from brain slices of P12-35 GAD67-eGFP or C57Bl/6J mice and expert-classified into five electrical types (see results). The detailed electrophysiological protocol has been published elsewhere. 26 Neurons were sampled from the ventrobasal complex of the thalamus (VPL and VPM nuclei) and the somatosensory sector of the reticular nucleus. 126,127 We used responses to step-like currents to build electrical models, ramp and noise currents to validate them, 26 along with excitatory postsynaptic-like currents (EPSC) injected into the dendrites. All the recordings were corrected for liquid junction potential by subtracting 14 mV from the recorded voltage. Neuron models Multicompartmental conductance-based models employed 3D morphological reconstructions. Active ion currents and a simple intracellular calcium dynamics model were distributed in the somatic, dendritic and axonal compartments. Only the axonal initial segment (AIS) -not the complete axon -was modeled. 24 The axons were substituted by a 60 mm stub constituted by two sections, five segments each. For each segment, the diameter was extracted from the original axon in order to preserve its tapering. Morphologies were divided into compartments of 40 mm maximal length. Specific membrane capacitance was set to 1 mF/cm 2 and specific intracellular resistivity to 100 Ucm.

Ion channel models
We included ion current models whose kinetics were obtained from previously published ion current models or published experimental data. All ion channel models were corrected for liquid junction potential and for simulation at different temperatures whenever possible. Simulation temperature was always set to 34 C.
The details of the ion channel kinetics and calcium dynamics used for low-threshold bursting neurons (TC and Rt) have been described elsewhere 26 and are summarized here. The type of ionic currents present in TC and Rt were: transient sodium current, delayed potassium current and low-threshold calcium from a previous model of TC neurons 142 ; h-current model (I h ) was built from published data 26,143-145 ; persistent sodium, based on an existing models 146,147 and published data, 148 A-type transient potassium was taken from an existing model, 146 based on published data 144 ; high-threshold calcium was based on published models and data, 144,146 SK-type calcium-activated potassium, was taken from previous models and published data. 147,149 Intracellular calcium dynamics was modeled with an exponential decay mechanism that linked low-threshold and high-threshold calcium currents to the calcium-activated potassium.
Since interneurons had firing patterns similar to cortical ones, we used the same ion channel models of the cortical microcircuit model, 24 which were based on existing models or published data. The type of ionic currents were transient sodium, 150 low-threshold calcium, 151 h-current, 152 persistent sodium, 153 transient potassium, 154 high-threshold calcium 155 and potassium Kv3.1. 156 The reversal potential of sodium, potassium and h-current were set to 50 mV, À90 mV and À43 mV, respectively.
Ion channel models were distributed uniformly and with different peak conductance values for somatic, dendritic and axonal compartments, except for the h-current in the interneurons, whose distribution increased exponentially from the soma to the dendrites. 24 Optimization of neuron models Five electrical models (e-models), corresponding to each electrical-type (e-type), were fitted using a multiobjective optimization algorithm using the Python library BluePyOpt. 26, 157 The free parameters of the model were the peak conductances of the different mechanisms, parameters of the intracellular calcium dynamics (time constant of decay and percent of free calcium, gamma) and the reversal potential of the passive mechanism that contributes to the resting membrane potential. Each e-model was fitted with an exemplar morphology.
The optimization objectives were the electrical features extracted from the electrophysiological recordings. All modeled neurons included transient sodium, persistent sodium, A-type transient potassium, delayed potassium, low-threshold calcium, highthreshold calcium, calcium-activated potassium (SK-type), h-current ion channel currents. All cellular compartments (i.e. somata, dendrites and axon initial segments) contained active membrane mechanisms. Rt_RC e-models followed the same approach as for VPL_TC cells, as previously published, 26 with the additional electrical features to quantify the deeper post-burst afterhyperpolarization observed in Rt_RCs. 158,159 To validate and test the generalization of the neuron models, we used features from current stimuli not used during the optimization phase.
The detailed experimental protocol and the type of current stimuli and features are described elsewhere, 26 and summarized here. For all the e-types, two hyperpolarizing steps (À20/À40% and À120/À140% of the threshold current) were used to constrain passive properties (input resistance, resting membrane potential) and current activated by hyperpolarization, e.g. h-current (sag amplitude). Three levels of depolarizing steps (150%, 200%, 250% of the threshold current) were used to constrain firing pattern (adaptation index or inverse of the first and last interspike intervals, spike count, mean frequency) and spike shape-related features (action potential amplitude, depth of the after-hyperpolarization, action potential duration). All these protocols were applied in combination with a hyperpolarizing holding current (to reach stable membrane potential of À84 mV, after liquid junction potential correction). When low-threshold bursting cells are hyperpolarized compared to their resting membrane potential and then stimulated, they fire stereotypical low-threshold bursts. One step (200% of firing threshold) on top of a hyperpolarizing current was therefore used to constrain the bursting response, while three depolarizing steps on top of a depolarizing holding current (to reach À64 mV) were used to constrain the tonic firing responses, as explained above. For reticular neurons, a new feature (initburst_sahp) was added for the afterhyperpolarization after the burst. For Rt and TC cells, two additional protocols without any current injection or only holding currents were used to ensure that the e-models were not firing without stimulus or with the holding current only.
Electrical features from the experimental recordings and model traces were extracted using the open source library eFEL (https:// github.com/BlueBrain/eFEL).
We considered a model a good fit to the experimental data if all the feature errors (i.e. the Z-scores) were below 3. Quality assurance of morpho-electrical models After fitting the five e-models, they were combined with the 92,970 morphologies generated as output of the morphology diversification step. An automated pipeline tested the e-models in combination with the different morphologies (me-models) and filtered out those that deviated significantly from the experimental electrical features. To decide which me-model was to be accepted, we used the repaired exemplar morphology (i.e. the morphology used during the optimization, after being repaired) as a benchmark: a memodel passed if it had all the feature errors were below 5 standard deviations of the repaired exemplar. 24 To account for the input resistance given by the different morphologies, we devised an algorithm, based on binary search, to find the appropriate holding and threshold current for each me-model. 26 In addition, we first ran this pipeline on a small subset of the morphologies generated after morphology repair. In this way, we could visually inspect if the accepted me-models were generating biologically plausible firing behavior and the reasons why other memodels had high feature errors. In some cases, after inspecting the me-model voltage responses, we set less stringent criteria on some features, to ensure that we had enough different me-models for building the microcircuit. At the same time, we set more stringent criteria to reject me-models that were active without any input, since we did not find neurons that were spontaneously active in our experimental recordings.

Measuring neuron density
Immunohistochemistry of Rt and VPL for cell counting We complemented the neuron densities values from the Blue Brain Cell Atlas 32 (RRID:SCR_019266) by counting neurons in adult mouse brain slices. The brain was cryosliced at 50 mm on the sagittal plane and stained following standard immunohistological procedures with antibodies anti-GABA (for inhibitory neurons), anti-NeuN (for neurons) and DAPI (for all cells), using an existing protocol. 24 The slices were imaged with a confocal microscope (Zeiss, 710). The immunohistology and imaging of the region of interest (ROI) was completed for one P21 C57B1/6J mouse. Semi-automated cell counting and cell densities The images were aligned to the Allen Reference Atlas to create proper boundaries for the Rt and VPL. We used Imaris software (Bitmap) to create the ROI, for counting the neurons and to estimate the volume for density calculation. For a chosen ROI, the software detected the difference of signal intensity, created a 3D shape around the detected cells and extracted statistics (e.g. count, positions) following given parameters. These parameters were defined by running multiple trials so that the results from semi-automated cell counting were as close as possible to those from manual cell counting. The semi-automated counting method results in very low error rates compared to manual counting (2.25%) and is less time consuming. A 3D shape of the entire ROI was created in order to extract the volume for density calculation. Neuron densities were calculated as the ratio between neuron counts in a ROI and the volume as calculated in Imaris for each slice. For modeling we used the average cell densities for Rt and VPL neurons.
Reconstructing the dimensions and structure of a thalamoreticular microcircuit Since the thalamus does not have a clear laminar structure, we approximated a thalamic microcircuit as a cylindrical volume having its base parallel to a portion of the Rt and its vertical dimension (y axis) running through the VPL and Rt (see Figure 1C).
The horizontal dimensions of the microcircuit were calculated from the density of dendritic fibers at the center of the circuit, following an approach published previously. 24 For each m-type, we began by considering all the morphologies (after repairing them) that had their somata located within 25 mm from the circuit center on the horizontal plane (XZ). We then increased the maximal distance in steps of 25 mm which resulted in an increase of dendritic densities at the center. The microcircuit horizontal dimension (radius) resulting from this process was 294 mm, corresponding to the distance where 95% of the asymptotical maximal density of reticular neuron dendrites was reached. As a comparison, considering only thalamocortical cell morphologies would have resulted in a circuit with radius 125 mm, while considering only interneurons the radius would have been 279 mm.
We used hexagonal boundaries with the same area as the resulting circle to facilitate tiling of multiple microcircuits, while keeping asymmetrical edge effects minimal. The resulting side of the hexagon was 323 mm and the longest diagonal (vertex-to-vertex) measured 646 mm.
To calculate the vertical dimension of the microcircuit, we extracted a 3D subvolume within the VPL and the Rt. We started from the thalamus parcellation of the Allen Brain Atlas version 3 (25 mm resolution). 29 A spherical coordinate system was fitted to the volume of the Rt, which can be approximated by a spheroidal surface. We chose a ROI located approximately in the middle of the VPL nucleus and computed the probability distribution of widths in the ROI for the VPL and Rt. The widths were calculated along the radius of the Cell Reports 42, 112200, March 28, 2023 27 Resource ll OPEN ACCESS spherical coordinate system. The resulting thickness corresponds to the median of the distributions, which was 550 mm for the VPL and 250 mm for the Rt.
Soma positions and me-type model assignment The horizontal and vertical extents resulted in a microcircuit having the shape of a hexagonal prism, that was 646 mm wide (at the widest point) and 800 mm high; 69% of the volume was occupied by the VPL and 31% by the Rt. This volume was then populated by defining somata positions according to the experimentally measured neuron densities in the Rt and VPL. The positions were distributed according to an algorithm based on Poisson disc sampling. 160,161 This algorithm avoids clustering normally obtained with sampling according to uniform distributions, by using a parameter for the minimum distance between points. To calculate the minimum distance, we used the cell densities to calculate the expected number of cell positions per voxel. Each soma position was assigned an m-type according to the excitatory/inhibitory fractions and an electrical model in agreement with the me-types composition ( Figures 1E-1G). Moreover, each position was associated with a random rotation around the y axis to be applied to each morphology.

Morphology placement
Our pool of experimental morphologies and the ones derived from the morphology diversification process contained morphologies with different sizes and shapes. Moreover, it contained TC morphologies whose somata was not located in the VPL nucleus and Rt morphologies whose axons were not arborizing in the VPL nucleus. We adapted a placement scoring algorithm 24 to ensure that each position was assigned a suitable morphology considering its geometrical properties and the microcircuit vertical dimension.
We thus defined placement rules that took into account the known properties of Rt and TC neurons' arborizations relative to the anatomical boundaries of thalamic nuclei. 162,163 Each reconstruction of TC and Rt neuron morphologies was manually annotated according to the placement rules. For TC cells, we identified the axonal arborization projecting to the Rt (and that should be located in the Rt part of the model). For Rt cells, the densest part of axonal arborization was annotated, which should be located in the VPL. For IN, the only constraint is that the full morphology should be contained within the VPL and not crossing into the Rt. 36 Each annotation was automatically carried over during the unraveling, repairing and diversification steps. Moreover, we included a stricter rule to avoid that Rt morphologies were located outside the top of the circuit boundary, with a 30 mm tolerance.
Given the placement rules, each morphology was assigned a score based on the microcircuit position and the constraints set by the placement rules. 24

Generating different microcircuit instances
We created five different microcircuit instances to assess the model robustness to different input parameters. The experimentallymeasured cell densities were jittered by +/À 5%, resulting in microcircuits with different total number of neurons and number of neurons for each m-type (see Figure 1F).
Reconstructing the synaptic connectivity of a thalamoreticular microcircuit Connectivity based on morphological appositions After placing the morphologies in the 3D microcircuit volume we generate the first version of the connectivity by detecting zones of geometrical overlap (''touches'') using an existing touch detection algorithm. 24,164 Briefly, this algorithm sub-divided the circuit 3D space into sub-volumes ensuring that each sub-volume contained the same amount of data, i.e. the same number of morphological segments. Each sub-volume was processed in parallel on different cores and written in parallel to disk. All geometrical overlaps were considered as touches if their distance was smaller or equal to 1 mm (''touch distance'').
Touches were then filtered according to biological rules: touches were allowed between all m-types, except between interneurons and reticular cells, because interneurons are only located in the thalamus and are not expected to have neurites extending into the reticular nucleus. 36 Touches between VPL neurons were removed, in agreement with experimental findings showing that excitatory connections between TC neurons disappear during development. 165 Interneurons also form axonal and dendritic inhibitory synapses. 116,166,167 For all other m-type combinations, touches were formed between presynaptic axons, postsynaptic dendrites and somata.
The same algorithm was used to detect touches between Rt_RC dendrites, i.e. the locations of putative gap junctions. Since gap junctions are established with close appositions of cell membranes, we used a touch distance of 0 mm in this case.
At the end of this process the resulting contacts (or ''appositions'') are normally higher compared to experimental findings and are pruned further to arrive at the final functional synapses. 33 Determining functional synapse positions We employed an existing algorithm to decide which appositions were to be pruned according to biological constraints. 33 The main constraints were the experimental bouton densities (number of boutons/axonal length) from 3D neuron reconstructions (n = 9 TC axons and n = 2 Rt axons) and the coefficient of variation of number of synapses per connections (i.e. the number of functional synapses, between a pair of neurons) from presynaptic INs and post-synaptic INs and TCs. 36 In the first two steps, the algorithm tried to match the distribution of synapses per connection, using the coefficient of variation of appositions per connections and the coefficient of variation of synapses per connection. Then, in step 3, it compared the current 28 Cell Reports 42, 112200, March 28, 2023 Resource ll OPEN ACCESS bouton density to the target value and removed multi-synaptic connections until the target value was matched. The number of synapses per connections, N func , was predicted from the number of appositions per connections (N app ) resulting from the previous steps, similarly to uncharacterized pathways in cortical microcircuitry. 33 N func was predicted from N app according to a simple formula (N func = 1 , N app ) for each m-type to m-type connection. We used a generalized coefficient of variation for N func of 0.9 for all connections, extracted from published data (Morgan and Lichtman, 36 Figure 2E). The coefficient of variation was combined with the predicted N func to calculate its standard deviation, as detailed in Reimann et al. 33 At the end of this pruning process, we verified that the bouton densities in the model matched the experimental ones (see Figure 2A). The shape of a geometric distribution for N func was a prediction from our touch detection process.

Connections from lemniscal and corticothalamic afferents
We followed an approach similar to the generation of thalamic input to the cortical microcircuit model 24 to model afferent synapses in the thalamus from the sensory periphery (medial lemniscus) and from the cortex. The algorithm uses volumetric bouton densities and the morphologies already placed in a circuit to map synapses from afferent ''virtual'' fibers to postsynaptic morphologies.
We built medial lemniscus (ML) and corticothalamic (CT) afferents separately for one microcircuit. Since data for lemniscal innervation in the mouse VPL was not available we calculated volumetric bouton density from data of mouse VPM, 38 see Figure 2D for the exact values. Volumetric bouton densities for the CT pathway were derived from known proportions between CT synapses and other synapses onto TC and Rt neurons, as found in electron microscope investigations. 39,40 Each synapse was assigned a virtual ML or CT fiber. We estimated a number of 2,601 ML fibers; this number took into account the ratio between the putative number of neurons from the dorsal column nuclei projecting to the thalamus 168 and the number of neurons in the VPL (see Jones 2 for a similar calculation). The number of CT fibers was 75,325, about ten times the number of thalamocortical fibers in a microcircuit. 64, 169,170 To take into account the correlation between synaptic inputs onto postsynaptic neurons innervated from the same afferent fiber, the mapping between postsynaptic synapses and fibers took into account their respective positions, i.e. synapses that were closer together were more likely to be innervated by the same presynaptic fiber. As in the neocortical microcircuit model, 24 the probability (P) that a synapse was assigned to a fiber depended on the distance between the synapse and the fiber: where S pre represents the mapping of a synapse S to the presynaptic fiber i, T pre is its spatial location, f i the spatial location of fiber i and s denoted the degree of spatial mapping, that was set to 25 mm.

Modeling synapse physiology
Stochastic synaptic transmission and short-term plasticity We used available paired-recording data and generalization principles to assign synaptic conductance values (g syn ) to match experimentally recorded PSP amplitudes. We predicted that single g syn from inhibitory neurons are in general small (e.g. 0.9 ± 0.23 nS for VPL_IN to VPL_IN), while conductances from VPL_TCs and lemniscal afferents are larger (>2 nS), consistent with being ''driver'' synapses, 171,172 while corticothalamic synapses have small g syn (<0.5 nS), but are facilitating.
To model synapse kinetics, we used existing models of synaptic currents 24 and included literature findings on decay time constants, reversal potentials and the relative contribution of AMPA, NMDA, GABA A and GABA B currents, summarized in Table S1. 118,[173][174][175][176] These models consist of a 2-state Markov process, with recovered and unrecovered states. When a pre-synaptic event occurs (pre-synaptic spike or spontaneous release) the synapse will release if it is the recovered state. If there is release, the synapse will transition to the unrecovered state. The ensemble average response is equivalent to the phenomenological Tsodyks-Markram model. 42, 177 The underlying assumptions were derived from the classical model of quantal synaptic release, in which each synapse is assumed to have N independent release sites, each has a probability p of releasing a single quantum q. 178,179 The number of release sites was assumed to be equivalent to the number of synapses per connection. 24 The detailed implementation of the synapse models can be downloaded from the neuron model packages in the Neocortical Microcircuit Portal 180 (rrid:SCR_022032).
We modeled short-term synapse plasticity with depressing (E2 and I2) and facilitating synapses (E1), see Figure 2H. In our experimental recordings, in agreement with experimental findings, all existing intrathalamic (between TC, Rt neurons and INs) and lemniscal connections were depressing, 27,171,181-183 while corticothalamic ones were facilitating. 64,183-186 When sufficient experimental paired recordings data were available, the parameters of the Tsodyks-Markram model of short-term synaptic plasticity were fitted (see above). The data used for fitting were the excitatory postsynaptic potentials (EPSPs) or inhibitory postsynaptic potentials (IPSPs) peaks amplitudes (or EPSCs/IPSCs in the case of voltage-clamp recordings), evoked by stimulating the presynaptic cell with a train of eight pulses followed by a recovery pulse (see Figure 2G). The parameters were: U -release probability, D -time constant of recovery from depression, and F -time constant of recovery from facilitation. Postsynaptic data were filtered and deconvolved for easier automatic identification of the peaks. 187  Synapse dynamic parameters in the model were different for each synapse and drawn for truncated Gaussian distributions. Spontaneous miniature potentials were modeled as independent Poisson processes at each synapse that triggered release at low rates (0.01 Hz).

Synapse models
Excitatory synaptic transmission was modeled with AMPA and NMDA receptor kinetics, and GABA A receptors were used for inhibitory connections. The rise and decay phases of the currents were described using mono-exponential functions. We used time constants from thalamic experiments performed at 34-35 C, when available, or from cortical synapses models when thalamic-specific ones were missing. The rise time and decay time constants for AMPA receptors were 0.2 ms and 1.74 ms, respectively. 189 For TC to Rt connections the AMPA decay time constant was 1.58 ms and CT afferents to Rt was 2.74 ms. 190 The rise and decay time constants of the NMDA component were 0.29 and 43 ms. 191 The magnesium concentration was set to 1 mM 192 and the reversal potential of the AMPA and NMDA currents was 0 mV. Experimentally measured ratios of NMDA and AMPA conductances were gathered from the literature and are summarized in Table S1. 173,174,190 Inhibitory synaptic transmission was modeled with GABA A receptor kinetics. The rise and decay time constants were 0.2 ms and 8.3 ms, respectively. 24 The reversal potential of GABA A current was set to À82 mV for all inhibitory pathways, except for connections onto postsynaptic TC neurons, where it was À94 mV, consistent with lower chloride reversal potentials in TC compared to Rt neurons. 193 Constraining synapse conductance values Synaptic conductance values were optimized by performing in silico paired recordings to match the postsynaptic potential (PSP) amplitudes measured experimentally whenever data was available, similarly to other morphologically detailed models. 24,188 For each pathway, 50 neuron pairs were simulated, and each pair was recorded for 30 trials. Experimentally characterized values in rodents are summarized in Table S2. For all other pathways, we extrapolated the quantal synapse conductances from similar pathways, according to the same generalization principles applied for short-term plasticity parameters (see Table S1).

Modeling gap junctions
Along with excitatory and inhibitory chemical synapses, the microcircuit included detailed gap junction (GJ) connectivity established between the dendrites of Rt neurons. We used the same touch detection algorithm described above to find appositions between Rt neuron dendrites and somata. Since we did not have any experimental data on the number of GJs between connected neurons or the density of GJs (number of GJs per unit length of dendrite or volume), in this first draft we randomly removed a certain fraction of GJs until we matched data on neuron divergence ( Figure 3A). To analyze the number of coupled neurons and their spatial properties, we reproduced the experimental protocol, 45 by analyzing a sample of 33 Rt neurons in a 90 mm vertical slice located at the center of the microcircuit.
Functionally, GJ were modeled as conductances that coupled the membrane potential of the adjacent morphological compartments (simple resistors). We predicted the value of GJ conductance for all GJs and validated their functional properties by comparing coupling coefficient values with experiments. [43][44][45][46] Once the structural properties of GJs-coupled neurons were validated, we performed in silico paired recordings and measured the coupling coefficients for each pair of neurons. We found that the mean coupling coefficients in the model compared well with the experiments for GJ conductance values of 0.2 nS (nS).
After adding GJs to the circuit, the input resistance of the neurons changed. To guarantee that the electrical properties of the neurons did not change, thus changing the responses to synaptic inputs, we devised an algorithm to compensate for the change in input resistance. 48 The algorithm changed the conductance of the leak current (g pas ) to restore the input resistance of the neuron before adding gap junctions. This compensation resulted in a different g pas value for each neuron.

Simulation methods and conditions
Simulation software and high-performance computing resources The reconstructed microcircuit was simulated using software based on the NEURON simulation package 194 (RRID:SCR_005393). A collection of tools and templates were written in order to handle simulation configuration, in silico network experiments and to save the results. We used the CoreNEURON simulator engine, 195 which has been optimized for efficient large-scale simulations. A typical simulation run of a microcircuit for 3,500 ms of simulation time took $45 min on 16 Intel Xeon 6140 CPUs (288 cores, with HyperThreading enabled).

Simulating in vivo-like conditions
To simulate spontaneous activity in in vivo wakefulness-like states, we activated lemniscal and CT fibers with Poisson spike trains at 25 and 4 Hz, respectively. To find the optimal value of lemniscal and cortical spontaneous rates we started from literature values 196,197 and then explored different combinations until we reached a biologically-plausible level of spontaneous activity in the model. We considered it as biologically plausible uncorrelated (asynchronous) firing in all cell types. From the literature, we know that spontaneous activity from L6 corticothalamic afferents is low in spontaneous conditions and is lower than lemniscal activity. We lowered the extracellular calcium concentration from 2 mM (in vitro-like conditions) to 1.2 mM, with the effect of reducing synapse release probabilities and PSPs amplitudes. 24 Since the dependence of PSP amplitudes on the extracellular Ca 2+ concentration in the thalamus is not known, we assumed that PSPs were dependent on calcium concentration in the same way as uncharacterized pathways as previously published (Markram et al., 24 Figure S11 -Intermediate [Ca 2+ ] o dependence). This condition was used in all simulations of in vivo wakefulness-like activity, if not stated otherwise.
To simulate cortical UP and DOWN states, we removed the background activity from the CT afferents (Figures 4 and S5).

Simulating lightly-anesthetized in vivo-like conditions
To simulate lightly anesthetized in vivo-like states, we followed the same methodology as for in vivo-like conditions, however the spontaneous firing induced at the thalamus through the ML afferents was reduced from 25 Hz to 10 Hz to reflect the presumed hyperpolarizing influence of the anesthetic (Figure 4). After exploring different values for the ML background activity, we found that for the 10 Hz ML background TC cells responded to the stimulus with low-threshold spikes of bursts as shown in the related publication. 63 Simulating in vitro-like conditions To simulate in vitro-like states, all neurons were left at their resting potentials (which ranged between À75 and À70 mV as recorded experimentally) and the only source of input was the spontaneous synaptic release from intrathalamic, medial lemniscus and corticothalamic synapses (at a rate of 0.01 Hz). The spontaneous synaptic release rate was taken from other modeling studies. After verifying that changing this value had no significant effect on the microcircuit network activity, we did not modify it in any simulation condition. The extracellular calcium concentration was set to 2 mM, as used in our in-vitro experiments. Simulating depolarization levels As a first approximation of the action of neuromodulators in the VPL and Rt, we applied constant current injections to the soma of each neuron. All neurons in the VPL or Rt regions were depolarized to the same target baseline membrane potential. The amplitude of the current was different for each neuron, to take into account the different input resistance of each morpho-electrical model.

Simulation analysis
The spectrogram in Figure S4 was calculated using the function scipy.signal.spectrogram, with inputs the sampling frequency of simulated membrane potential (10 kHz), interval = 5000, overlap = 0.99, and the other parameters with the default values.
Burst probabilities ( Figure S4) were calculated as the ratio between the numbers of spikes belonging to a burst and the overall number of spikes. We considered a spike belonging to a burst when the interspike intervals were % 15 ms and the first spike in the burst was preceded by a pause R 50ms. For this analysis, we considered neurons that had a baseline activity between 1 and 20 Hz, as shown in corresponding publication. 74 To analyze the percentage of neurons firing for each m-type during each cycle of the oscillation ( Figure 6C) we started by finding the oscillation peaks. The peaks were extracted from the firing rate histograms as input, using the scipy.signal.find_peaks function. We then added a peak corresponding to the time of the stimulus injected into Rt neurons (cycle 0). Spikes for each m-type were then assigned to the different cycles if they occurred within 30 ms of the oscillation peak.
Oscillation strength (Figure 7) was calculated as the maximal value of the power spectral density (PSD). The PSD was obtained using the function scipy.signal.periodogram.
To calculate the oscillation duration (in ms) we used firing rate histograms for all the neurons and extracted their peaks, using the scipy.signal.find_peaks function (Figure 7). Peaks were counted only if they were significantly higher than baseline firing rates. Oscillation duration was then calculated as the time difference between the last and the first peak.
To calculate oscillation frequency (Figure 7), we computed the normalized autocorrelation of the firing rate histograms and extracted the time (oscillation period) corresponding to the first non-zero peak. 73 The inverse of the oscillation period corresponded to the oscillation frequency.

Supplemental information
Thalamic control of sensory processing and spindles in a biophysical somatosensory thalamoreticular circuit model of wakefulness and sleep Table S1. Synapse kinetics and short-term plasticity parameters, related to Figure 2 Synaptic parameters for all pathways in the model. Quantal synaptic conductance g syn (in nanosiemens nS), τ d is the decay time constant of AMPA and GABA A currents for excitatory and inhibitory connections. U SE (utilisation of synaptic efficacy, analogous to release probability), D (time constant of recovery from depression), F (time constant of recovery from facilitation) are the short-term plasticity parameters. Values are expressed as mean ± standard deviation. All the parameters were fitted to in-house paired-recordings or generalized from similar pathways.          (A) Simulated UP and DOWN states evoke reticular and thalamic depolarizations through afferent input, resulting in the "ping-pong" generation of spindle-like oscillations (as seen in Fig.  11). To approximate the differential effects of neuromodulators (e.g. acetylcholine) onto Rt_RC and VPL_TC we applied constant currents to depolarize Rt_RC and hyperpolarize VPL_TC cells. This resulted in spindle-like oscillations (left) being abolished (right). (B) Example single cell recordings from the simulation in A, note that while Rt_RC cells fire preferentially low threshold bursts during the cortical UP states (left), they transition to single spike modes when depolarized (right). The change in Rt_RC firing mode and hyperpolarization of VPL_TC cells resulted in a significant decrease of large amplitude IPSPs in VPL_TC cells and rebound bursts.