Cross-Frequency Slow Oscillation–Spindle Coupling in a Biophysically Realistic Thalamocortical Neural Mass Model

Sleep manifests itself by the spontaneous emergence of characteristic oscillatory rhythms, which often time-lock and are implicated in memory formation. Here, we analyze a neural mass model of the thalamocortical loop in which the cortical node can generate slow oscillations (approximately 1 Hz) while its thalamic component can generate fast sleep spindles of σ-band activity (12–15 Hz). We study the dynamics for different coupling strengths between the thalamic and cortical nodes, for different conductance values of the thalamic node's potassium leak and hyperpolarization-activated cation-nonselective currents, and for different parameter regimes of the cortical node. The latter are listed as follows: (1) a low activity (DOWN) state with noise-induced, transient excursions into a high activity (UP) state, (2) an adaptation induced slow oscillation limit cycle with alternating UP and DOWN states, and (3) a high activity (UP) state with noise-induced, transient excursions into the low activity (DOWN) state. During UP states, thalamic spindling is abolished or reduced. During DOWN states, the thalamic node generates sleep spindles, which in turn can cause DOWN to UP transitions in the cortical node. Consequently, this leads to spindle-induced UP state transitions in parameter regime (1), thalamic spindles induced in some but not all DOWN states in regime (2), and thalamic spindles following UP to DOWN transitions in regime (3). The spindle-induced σ-band activity in the cortical node, however, is typically the strongest during the UP state, which follows a DOWN state “window of opportunity” for spindling. When the cortical node is parametrized in regime (3), the model well explains the interactions between slow oscillations and sleep spindles observed experimentally during Non-Rapid Eye Movement sleep. The model is computationally efficient and can be integrated into large-scale modeling frameworks to study spatial aspects like sleep wave propagation.


INTRODUCTION
Sleep marks a pronounced change of the brain state as one of the vital means of persisting mental and physical health (Laureys et al., 2007). It manifests itself by the spontaneous emergence of characteristic oscillatory rhythms, most visible in the electroencephalogram (EEG) but also noticeable in intracellular recordings, electrooculography (EOG), and electromyography (EMG) (Steriade et al., 1993;Schomer and Da Silva, 2012). Distinct oscillatory features form the basis for sleep classification into several stages: rapid eye movement sleep (REM) and three stages of non-REM (NREM) sleep (N1 through N3) (Silber et al., 2007;Berry et al., 2012). The NREM sleep stages exhibit characteristic patterns: the N1 stage consists of slow eye movements and lowamplitude low-frequency δ-band EEG activity (Krishnan et al., 2016); sleep spindle oscillations dominate the N2 stage with a waxing and waning envelope and an underlying oscillation in the σ -band (∼12-15 Hz) (Steriade, 2003), in the case of "fast" spindles, or ∼9-12 Hz for "slow" spindles (De Gennaro and Ferrara, 2003;. Sleep spindles are also observed in the deeper, N3 sleep stage, albeit with lower power in the fast spindle frequency range  and lower density (Fernandez and Lüthi, 2020). Finally, slow oscillations (SOs), which is an alteration of active (UP) and silent (DOWN) cortical states at ∼1 Hz frequency, govern the deepest N3 sleep stage but are also present in the N1 and N2 sleep stages, albeit with a lower spectral power (Achermann and Borbely, 1997;Amzica and Steriade, 1997). Ripple oscillations are the second hallmark of the N3 sleep stage. Ripples are fast oscillations (80-140 Hz) that occur in hippocampal networks, often accompanied by a sharp wave, and signify reactivations (memory replay) of neural ensembles in these networks (Axmacher et al., 2008;Klinzing et al., 2019). The precise coordination of the slow oscillations, spindles, and ripples was shown to be vital for memory consolidation, of which the main manifestation is the reactivation of specific activity patterns, i.e., memory replay, during sleep (Walker and Stickgold, 2004;Popa et al., 2010;Bendor and Wilson, 2012;Rasch and Born, 2013). The hierarchical nesting of spindle waxing periods to the depolarized cortical UP states, mediated by the thalamocortical circuitry, is essential for consolidation by providing a "window of opportunity" and favorable conditions for plasticity for transferring episodic memories from short-term hippocampal storage to longer-term neocortical storage (Rosanova and Ulrich, 2005;Cox et al., 2012). Notably, it was recently shown that the interplay of the aforementioned rhythms is also vital for non-hippocampus-dependent consolidation (Klinzing et al., 2019) in rodents (Sawangjit et al., 2018), and humans (King et al., 2017).
The basis of slow oscillations consists of a widespread alternation of hyperpolarization and depolarization activity in neocortical networks (Sanchez-Vives and McCormick, 2000;Steriade, 2003;Peyrache et al., 2012). In contrast, spindles are generated by the interaction of inhibitory reticular thalamic and excitatory thalamocortical neurons (Timofeev and Bazhenov, 2005). Spindles occur in the isolated thalamus both in vivo and in vitro (Kim et al., 1995;Timofeev and Steriade, 1996). However, the cortex can also become actively involved in their initiation and termination (Bonjean et al., 2011) and their long-range synchronization (Contreras et al., 1997;Bonjean et al., 2012). Therefore, in the thalamocortical modeling, the assumption is that the cortical part of the model generates slow oscillations and the thalamic part of the model, in turn, generates spindles (Robinson et al., 2002;Suffczynski et al., 2004;Schellenberger Costa et al., 2016). Moreover, slow oscillations induce thalamic spindles, which then reflect back to the cortex (Oyanedel et al., 2020). Furthermore, it has been shown that the phase of SOs modulated the spindle power such that it exhibited UP and DOWN states similar to SOs themselves, with positive peaks during the depolarizing SO UP state, close to (or slightly before) the SO peak .
Due to the computational advantage, the ability to elucidate the dynamical repertoire, and the ability to describe macroscopic phenomena comparable with neuroimaging datasets (Deco et al., 2008;Touboul et al., 2011), herein we focus on neural mass models (NMMs hereafter). In particular, Suffczynski et al. (2004) proposed a thalamocortical NMM to explain the relationship between spindle-generating and spike-wavegenerating activities. More recently, Cona et al. (2014) proposed a new NMM to describe a "sleeping" thalamocortical system with tonic and bursting firing modes within thalamic neurons while Schellenberger Costa et al. (2016) probed a thalamocortical model that generates spindles and K-complexes. To this date, however, no study focused on the cross-frequency slow oscillation-spindle interaction in the model setting that would mimic the empirical findings (Ladenbauer et al., 2016Helfrich et al., 2018). We aim to close this gap in the present study.
This study effectively extends the results in Schellenberger Costa et al. (2016) by using a biophysically realistic model for the cortical node capable of generating in vivo-like slow oscillations (Cakan and Obermayer, 2020). By merging two modeling approaches, we also highlight the applicability of hybrid modeling approaches and their advantages, mainly the extensibility with other mass models. In our case, using an already explored cortical node facilitates the exploration of biophysically realistic stimulation protocols of the thalamocortical model. We present a thorough dynamical investigation of the thalamocortical model in the NREM sleep setting. Both cortical and thalamic modules are biophysically realistic while adhering to the notion of mass models and keeping their computational advantages. In addition to investigating the dynamical landscape, we study the nature of cross-frequency coupling between slow oscillations and sleep spindles, which are believed to set the stage for successful memory consolidation and transfer from the hippocampus to the neocortex (Sirota et al., 2003;Ji and Wilson, 2007).
The following sections first introduce the thalamic and cortical models, their basic dynamical properties, and their respective state spaces. Next, we investigate the effects of perturbations on both isolated models to understand their response to external stimuli. Finally, we study the full thalamocortical loop, including both feedforward and feedback connections. In the full model investigation, we focus on the interaction between cortical slow oscillations and thalamic spindles, the cause-effect mechanical understanding of the spindle and SO generation and timing, and, finally, their cross-frequency coupling.

Model Design
The architecture of the thalamocortical model is shown in Figure 1. It consists of two thalamic neural populations, namely FIGURE 1 | Schematic of the thalamocortical motif. The cortical node (top row) consists of one excitatory (EXC) and one inhibitory (INH) population. It is coupled to the thalamic node with its thalamocortical relay (TCR) and thalamic reticular nuclei (TRN) populations. Excitatory (inhibitory) populations are shown in orange (blue). Excitatory synapses are depicted with arrows, inhibitory synapses are depicted with filled circles. Squares depict noisy background inputs.
thalamocortical neurons (TCR) and thalamic reticular nucleus (TRN), which act as excitatory and inhibitory populations, respectively. The cortical module consists of mean-field approximation of excitatory and inhibitory exponential integrate-and-fire neurons grouped into two populations, of which the excitatory subpopulation exhibits a somatic spikefrequency adaptation mechanism. Both thalamic and cortical nodes rely on the notion of an empirical firing rate which effectively replaces the complex individual spiking dynamics of neural populations.
The model connectivity, both feedback and feedforward, relies on the fast ionotropic excitatory and inhibitory synapses conveyed by the AMPA and GABA A receptors. The cortical node has feedback and feedforward connections, whereas in the thalamic module, only the TRN possesses feedback connections (as shown in Figure 1) since thalamic relay cells generally do not form local connections within the population (Jones, 2001).
For the connection between the thalamus and the cortex, we assume that the long-range afferents from the cortical excitatory population project to both thalamic populations and that the TCR population projects to the cortical excitatory population, as depicted in Figure 1. We set delays of these long-range connections to a physiologically realistic value of 13 ms (Roux et al., 2013;Schellenberger Costa et al., 2016).

Thalamic Model
The thalamic model follows the approach detailed in Schellenberger Costa et al. (2016). The evolution of the mean membrane potential V α of the thalamic populations α ∈ {t, r} is described by with membrane time constant τ , synaptic input rate w e (w i ) that scales synaptic inputs s e (s i ) for excitatory (e) and inhibitory (i) synapses, the corresponding Nernst reversal potential E e (E i ), and the membrane capacitance C m . The mean membrane potential V α is then converted to a firing rate r(V) by a sigmoidal transfer function, with maximum firing rate r max , firing threshold θ , and gain coefficient σ . Both thalamic populations contain additional intrinsic currents [I intrinsic in Equation (1)] because spindle oscillations require rebound burst activity. Rebound bursting is impossible with a monotonic firing rate function and demands additional mechanisms. Following Schellenberger Costa et al. (2016), we employ the Hodgkin-Huxley-type extension, which has been derived from integrate-and-fire-or-burst neurons (Langdon et al., 2012).
The intrinsic currents in both thalamic populations include a potassium leak current, and a T-type calcium current, which de-inactivates upon depolarization. In both definitions, g LK and g T denote the conductance of the respective intrinsic current, E K and E Ca their respective Nernst reversal potential, and m 2 ∞ and h the gating functions of T-type calcium current. Both currents are essential for the generation of low-threshold spikes and rebound bursts (as shown in Figure 2). The definition of I T for TRN follows Destexhe et al. (1996b), while I T within the TCR population is given in Destexhe et al. (1998). In addition, the TCR population contains a hyperpolarization-activated cationnonselective current I h : with g h being its conductance, E h its Nernst reversal potential, m h1 and m h2 the gating functions, and g inc the conductivity scaling. This hyperpolarization-activated current is responsible for the waxing and waning of spindle oscillations in the isolated thalamus (Destexhe et al., 1996a). Synaptic transmission in the thalamic model is conveyed by conductance-based synapses, where the spike rate r k ′ of a presynaptic population k ′ elicits a postsynaptic response s lk in population k given by N kk ′ is the connection strength between the presynaptic, k ′ , and postsynaptic, k, populations. φ ′ (t) represents a background noise input, ⊗ denotes a convolution, and α l (t) is an alpha function given by representing the synaptic response to a single spike. γ l is the decay constant of the synaptic response, and l ∈ {e, i} denotes the type of synapse, i.e., excitatory AMPA or inhibitory GABA A . The convolution in Equation (6) is replaced in the numerical simulations by the second-order ordinary differential equation (ODE) The background noise input φ ′ (t) is modeled as an Ornstein-Uhlenback process with zero drift but finite variance and represents unresolved processes in our model, e.g., afferents from other brain regions that are not explicitly modeled here. In summary, the thalamic node is described by the set of equations where subscripts t (r) represent the TCR (TRN) population. Note that the background noise is included on the excitatory synaptic input at the TCR population, s et . The full set of equations is given in Supplementary Material, with all parameters summarized in Supplementary Table S1.

Cortical Model
The adaptive exponential integrate-and-fire (AdEx) neuron model (Brette and Gerstner, 2005) forms the basis for the derivation of the cortical mass model. Each population α ∈ {E, I} possesses N α neurons and the membrane voltage of neuron j in the population α is governed by C m,cVj = I ion (V j ) + I syn (t) + I j,ext (t) (11) The first equation describes the temporal evolution of neuron's j membrane voltage V j as a function of its internal current dynamics conveyed by I ion (V j ), its synaptic current I syn (t), and a background external current I j,ext (t) received from neural populations not specified by the computational model. C m,c denotes the cortical neuron's membrane capacitance. Note that inhibitory population α = I does not include the adaptation current. The first term in Equation (12) expresses the voltagedependent leak current with leak conductance g L and leak reversal potential E L ; the second term describes the exponential spike initiation mechanism with slope factor T and exponential threshold V T . Finally, the last term describes the somatic adaptation current, I A,j (t) (Equation 13), with subthreshold adaptation a, adaptation reversal potential E A , and adaptation time scale τ A . When the membrane voltage crosses the spiking threshold, V j ≥ V s , the voltage is reset with V j ← V r , clamped for the refractory time T ref , and the spike-triggered adaptation increment, b, is added to the adaptation current, I A,j ← I A,j + b.
The synaptic currents are described by a sum of excitatory and inhibitory contributions as with the coupling strength J αβ from population β to α, and the fraction of active synapses s i,αβ ∈ [0, 1]. The synaptic dynamics is then given bẏ where G ij is a random binary connectivity matrix. The first sum is over all afferent neurons j, while the second sum is over all incoming spikes k from neuron j emitted at time t k after a delay d αβ . The (1 − s i,αβ ) term in Equation (15) acts as a saturation term and integrates all incoming spikes only if s i,αβ < 1, i.e., only if there is a synaptic capacity available.
In the mean-field approximation, the distribution p(V) of the membrane potentials and the mean population firing rate r can be calculated using the Fokker-Planck equation in the thermodynamic limit, with the number of neurons N → ∞ (Brunel, 2000;Cakan and Obermayer, 2020). However, using the model reduction scheme described in Augustin et al. (2017) and Cakan and Obermayer (2020), we exploit the low-dimensional linear-nonlinear cascade model, where for a given mean membrane current µ α with standard deviation σ α , the meanV α of the membrane potentials, the adaptive timescale τ α , and the population firing rate r α in the steadystate can be captured by a set of nonlinear transfer functions (µ α , σ α ) (Richardson, 2007), i.e., for α ∈ {E, I}. In the case of the excitatory population (α = E), the mean adaptation currentĪ A is subtracted from the mean membrane current µ E in the computation of population firing rate, mean membrane voltage, and adaptive timescale via transfer functions, i.e., µ E → µ E −Ī A /C m,c . These transfer functions can be precomputed for a specific set of single AdEx neuron parameters (Augustin et al., 2017) and are shown in Supplementary Figure S1. The average population currents in the mean-field approximation are given by where the dynamics of mean membrane current µ α depends on the synaptic current µ syn α and an external noisy current µ ext α , which enters the system in the form of Ornstein-Uhlenbeck process with mean drift µ α , standard deviation σ α , and time scale τ OU . ξ α is drawn from a random Gaussian white noise process with zero mean and unit variance. In the subsequent text, we refer to mean drifts as µ E (µ I ) as an input to the excitatory (inhibitory) population with the original physical units of mV/ms. After multiplying µ α with membrane conductance C, we obtain input currents in physical units of A, but we will omit the multiplication by C in subsequent text and treat all µ α with units of A.
The mean of the fraction of active synapsess αβ obeyṡ where the mean r αβ and the variance ρ αβ of the effective input rate from population β to population α for a spike transmission delay d αβ are given by Finally, the current variance σ 2 α and the variance of the fraction of active synapses σ 2 s,αβ are given by The parameters are summarized in Supplementary Table S2.

Connecting the Thalamic and Cortical Models
The thalamic and cortical models introduced before are coupled into a thalamocortical model with feedback and feedforward connections via their firing rate. More concretely, the excitatory firing rate from the cortical node enters the dynamics of both TCR and TRN populations in the thalamus with a connection strength N ctx→thal (cf. Figure 1). In particular, the excitatory cortical firing rate comes in Equation 8 for the excitatory AMPA synapse with N kk ′ = N ctx→thal and r k ′ = r E (t − d ctx,thal ), with d ctx,thal being the thalamocortical delay of 13 ms. For the thalamocortical connection, we connect the excitatory firing rate from the TCR population onto excitatory population of the cortical model (cf. Figure 1). In particular, the TCR firing rate enters the cortical dynamics in Equation 24 with r β = N thal→ctx · r TCR , and d αβ = d ctx,thal . The transmission delay in the thalamocortical direction is set to the same value as in the corticothalamic direction, i.e., 13 ms.

Numerical Simulations
The whole delayed dynamical equations system [Equations (1)-(27)] was integrated with the forward Euler method. If not mentioned otherwise, simulated time was t = 30 s with an integration timestep of dt = 0.01 ms. After integration, time series were subsampled at dt samp = 10 ms. The thalamocortical model was simulated using the neurolib library (Cakan et al., Color-coded number of spindles per second as a function of conductance parameters g h and g LK for 4 different noise levels σ TCR . The model was simulated for 65 s. The first 5 s were not included in the statistics. (B) 10 s time series from the thalamic node model (red traces for TCR firing rate, blue traces for TRN firing rate) from the three regions marked in (A) with different parameters: g LK = 0.018 mS/cm 2 (I), g LK = 0.031 mS/cm 2 (II), g LK = 0.024 mS/cm 2 (III), and g h = 0.062 mS/cm 2 (all three regions). For other model parameters, see Supplementary Table S1. 2021). neurolib is a computational framework for wholebrain modeling written in Python. It provides a set of neural mass models and is designed to be extendable and allows for easy implementation of custom mass models. Moreover, it supports heterogeneous brain modeling by coupling more than one type of mass model together. It offers a custom parameter exploration and optimization module for fitting models to multimodal experimental data using evolutionary algorithms. All subsequent analyses were also done in Python, and the repository with the model and analysis code is available at https://github.com/ jajcayn/thalamocortical_model_study.
Background noise inputs are implemented as an independent Ornstein-Uhlenbeck processes (Bibbona et al., 2008) with the mean drift µ α and SD σ α , for α ∈ {E, I, TCR}, integration timestep dt. ξ is drawn from a random Gaussian white noise process with zero mean and unit variance. The Ornstein-Uhlenbeck processes, x, were pre-integrated and then inserted into the thalamic state equations as φ ′ into Equation (8), and to the cortical equations as µ ext α into Equation (22).

Spindle Detection From the Model Output
For automated spindle detection from model output, we used a modified version of the A7 spindle detection algorithm described in Lacourse et al. (2019) and implemented in the Yet Another Spindle Algorithm (YASA) package (Vallat and Jajcay, 2020) for Python.
Since the A7 spindle detection algorithm is designed to analyze empirical data, e.g., EEG, MEG, ECoG, or LFP (cf. Warby et al., 2014;Lacourse et al., 2020), we made adjustments to the algorithm parameters. For detecting spindles on the cortical model output, we lowered the threshold for the duration from 0.5 to 0.3 s and for the relative power in the fast spindle band from 0.2 to 0.15 s.

Cross-Frequency Coupling (CFC) Measures
In order to quantify the coupling between the phase of the slow oscillation and the spindle amplitude, we compute the Kullback-Leibler modulation index (KL-MI, cf. Tort et al., 2010) and the mean vector length (MVL, cf. Canolty et al., 2006). For quantifying the phase-phase CFC, we use the phase-locking value (PLV, cf. Cohen, 2008) and the mutual information (MI, cf. Paluš, 1997;Jajcay et al., 2018) between time series of instantaneous phases. We first filter the modeled outputs using a one-pass, zero-phase, non-causal finite impulse response (FIR) bandpass filter, implemented in the mne Python package (Gramfort et al., 2013) (either adapted to the SO, 0.1-3.0 Hz, or the fast spindle frequency, 12-15 Hz, ranges). The filtered signal is then passed to the Hilbert transform, which provides us with the estimate of the complex analytic signal s a (t) = s(t) + iŝ(t), where s(t) is the real-valued input signal (model output), andŝ(t) is the Hilbert-transformed signal. The instantaneous phase, φ(t), and amplitude, A(t), are then given by The KL-MI estimates the phase-amplitude coupling by computing the Kullback-Leibler divergence (Kullback and Leibler, 1951) between the distribution of spindle amplitudes A spindle (t) over the slow oscillation phase bins φ SO (t) and a uniform distribution (the null hypothesis of no phase-amplitude coupling). The KL-MI is defined by where N is the number of phase bins and H(P) is the Shannon entropy of the amplitude distribution. The MVL is computed by averaging the complex time series constructed by multiplying the spindle amplitude, A spindle (t), with the term containing the phase, φ SO , of the slow oscillations, thus The length (real part) of the MVL vector quantifies the strength of phase-amplitude coupling. Its phase denotes the mean phase of the slow oscillation at which the spindle amplitude is the strongest.
The PLV, as a measure of phase locking between two oscillations, is computed by temporally averaging the phase differences in unit circle between the phase of the slow oscillation and the phase of the spindle oscillation: As with the MVL, the length of the PLV vector indicates the strength of phase-locking, while the phase represents the phase shift.
Finally, the MI of two discrete variables (in our case, the time series of estimated phases) can be computed using where p(x, y) signifies the joint probability mass function of {X, Y}, and p(x) and p(y) represent the marginal probability mass functions of X and Y, respectively. All three probability mass functions are estimated using an equiquantal binning algorithm (Paluš, 1995) with 16 bins.
To test for statistical significance, we computed 1,000 Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogates from a randomization procedure that preserves the power spectrum and the amplitude distribution of the original time series (Schreiber and Schmitz, 2000). The surrogates were constructed by iterative replacements of Fourier amplitudes with the values from the original time series and by rescaling the distribution match between the distribution and the power spectrum of the original data. We opted to use IAAFT surrogates since our modeled data are firing rates with non-Gaussian distribution and do not meet the criteria for using the basic Fourier Transform surrogates. After constructing the surrogates, we computed a given CFC measure for each surrogate time series, yielding an empirical null distribution. Finally, we compared a result on the modeled output to this null distribution to obtain the empirical p-value.

Dynamical Repertoire of the Thalamic Model-Spindle Oscillations
The isolated thalamic node model generates spindle oscillations (as shown in Figure 2). As described previously (Bazhenov et al., 2002;Destexhe and Sejnowski, 2003;Schellenberger Costa et al., 2016), spindle oscillations emerge through the reciprocal interaction of the TRN, which acts as a pacemaker, and the TCR, which mediates spindle propagation to the cortex (Rasch and Born, 2013;Fernandez and Lüthi, 2020). A low-threshold burst discharge in the TRN population causes synchronous and robust inhibition of the TCR, which activates its T-type calcium current. The subsequent activity rebound drives the TRN population to elicit additional low-threshold bursts. Additionally, activating the T-type calcium current requires a strong tonic hyperpolarization caused by the potassium leak current (Destexhe et al., 1996a;Bazhenov et al., 2002).
The waxing and waning structure of the spindle oscillations is caused by the after depolarization in the TCR, mediated by the hyperpolarization-activated cation-nonselective channels, represented by the I h current (Fernandez and Lüthi, 2020). A sequence of low-threshold spikes leads to the build-up of calcium in the TCR cells, which increases the effective conductivity g h of I h . The depolarization of the TRN additionally counteracts its ability to produce a low-threshold spike, which conclusively ceases the spindle oscillation (Contreras et al., 1997;Lüthi and McCormick, 1998). Sleep spindle termination also involves cortical and brain stem mechanisms (Fernandez and Lüthi, 2020), which we omit for brevity here. The two mechanismsone creating fast spindle oscillation between 12 and 15 Hz, the other responsible for the waxing and waning structure with frequency < 1 Hz -are visualized in Figure 2.
We probed the thalamic model for a set of parameters that convey the generation of spontaneous spindles. Figure 3 summarizes spindle activity in the thalamic model as a function of three parameters: the conductances of the potassium leak, g LK , and rectifying, g h , currents, and the variance, σ TCR , of the background noise input. In the noise-free case (Figure 3A, σ TCR = 0.0 mV/ms 3/2 ), we observe two regions where spindle oscillations emerge spontaneously (marked as region I and region II in the figure). In these two regions, the interaction between the fast T-type current and its slow modulation via the rectifying current leads to spindle oscillations. The regions differ in the value of g LK (cf. Schellenberger Costa et al., 2016). Spindle oscillations in the region I are longer, more symmetric, and have a longer inter-spindle interval (Figure 3B rows I vs. II).
When we introduce noise to the thalamic model (second to the fourth column of Figure 3), we observe qualitative changes in the spindling behavior. The most noticeable is the emergence of the waxing and waning cycle as a function of noise strength in areas neighboring the original spindle-promoting regions. This is clearly seen in the third row of Figure 3B (region III), as the spindles emerge for sufficiently strong background noise. In the spindle-promoting regions I and II, noise randomizes spindle timing: it can push the TCR population into a spindle event or, reversely, cease spindling faster than the slow I h -driven negative feedback loop alone. The background noise in our model can act as a depolarizing or hyperpolarizing force, thus can either speed up the activation of I h and cease spindling, or, reversely, slow it down and prolong the spindling periods. Sufficiently strong noise may also change the number of spindles in the spindlepromoting regions I and II as seen for noise strengths σ TCR = 0.01 mV/ms 3/2 and larger.
To the best of our knowledge, there is no experimental evidence of how spindle properties change as a function of g LK . However, we know that spindle density ranges between 2-10 spindles/min in EEG data (Purcell et al., 2017;Fernandez and Lüthi, 2020), and spindles can have various levels of symmetry. Although, there seems to be a continuous range of spindle density and symmetry rather than two clusters in the conductance phase space. We consider the two separated spindling regions as a pure model feature (which is confirmed by the bifurcation analysis of the thalamic model in Schellenberger Costa et al., 2016, Figure 2), where the model can generate more symmetrical and less frequent spindles (region I, lower g LK , hence lower hyperpolarization), and less symmetrical and more frequent spindles (region II, higher g LK ).
Outside of the spindle-promoting regions I and II, the thalamic model displays continuous oscillation in the fast spindle band due to hyperpolarization induced rebound bursts (as in region III Figure 3B, but also for larger and smaller values of g h , hence above and below regions I and II). This behavior does not exhibit the waxing and waning structure since the T-type current dominates and I h is not strong enough to sufficiently depolarize the TRN population to cease the oscillation. Furthermore, for larger values of g LK (to the right of region II in Figure 3), the I LK current dominates, and the thalamic model switches to slow, δlike rate oscillations. In the corners of the state space spanned by the conductances g LK and g h , the thalamic model exhibits a stable fixed point behavior with constant firing rates (e.g., for maximal conductances g h = g LK = 0.08 mS/cm 2 , the TCR exhibit down state with a constant rate of 10 Hz, while for g h = g LK = 0 mS/cm 2 the TCR exhibit up state with a constant rate of 116 Hz).
In its default parametrization, our model exhibits a spindle frequency of approximately 13 Hz. The spindle frequency in the thalamic model depends on the conductance of the Ttype calcium current, g T (Schellenberger Costa et al., 2016). By changing the conductance value in the TCR population, the model is able to reproduce spindle frequencies in the whole range of fast spindle oscillations (approximately 12 -15 Hz). On the other hand, changing the value of g T in the TRN population has only a minor effect on the spindle frequency. This effect does not significantly change qualitatively with the introduction of noise into the thalamic model. Other spindle parameters, such as its duration, amplitude, or symmetry, are approximately invariant with respect to changes in g T in both populations.

Dynamical Repertoire of the Cortical Model-Slow Oscillations
The cortical node as a motif of delay-coupled excitatory and inhibitory populations with somatic spike-frequency adaptation can be parametrized into four different dynamical regimes dependent on the mean external input currents µ α to both populations α ∈ {E, I} (Cakan and Obermayer, 2020). The 2D slices of bifurcation diagrams of interest are shown in Figure 4A. When the excitatory input is weak, the system is in its DOWN state: a fixed point solution with a low firing rate of the excitatory population. As we increase the external input to the E population, the system undergoes a supercritical Hopf bifurcation, and a limit cycle emerges. For a weak background current to the I population, the interaction between inhibition and excitation creates an E-I oscillation with a frequency of ∼25 Hz. For a stronger background current to the I population, the interplay between adaptation and excitation gives rise to a slow limit cycle with frequencies of around 2 Hz and lower. The system is in a stable fixed point for stronger background current to the E population again, albeit with a higher firing rate, the so-called UP state (Cakan and Obermayer, 2020).
When noise is added to a system parametrized close to the border of the slow limit cycle and the DOWN state, we observe a DOWN state with occasional irregular UP state excursions (i.e., noise pushes the system intermittently into the limit cycle - Figure 4B middle). Reversely, along the border between the limit cycle and the UP state, the system exhibits a UP state with irregular DOWN state excursions caused by noise (Figure 4B  right). From the perspective of modeling slow oscillation during deep sleep, the optimal operating point of the cortical model is close to the border between the UP state and the slow limit cycle, where the cortical node undergoes irregular DOWN state excursions .

Thalamic Model Driven by Cortical Input
We first investigate the state space of the thalamic model by studying the effects of external (cortical) inputs without FIGURE 4 | Emergence of slow oscillations in the cortical model. (A) Slices of state space spanned by the mean external input currents µ α , α ∈ {E, I}, to the excitatory and inhibitory population with no (σ E = σ I = 0.0 mV/ms 3/2 , top panel) and finite (σ E = σ I = 0.05 mV/ms 3/2 , bottom panel) noise. The panels show (left to right) the following: the maximum firing rate of the excitatory population, the dominant frequency of its rate oscillations (computed as a frequency with maximum power in its Welch spectrum; white color denotes no oscillations, i.e., fixed point dynamics), and the difference between its maximum and minimum firing rates. White lines indicate boundaries between fixed points (UP and DOWN) and limit cycles (LC EI and LC aE ) computed on the noise-free case. (B) Simulated time series (from E in red, from I in blue) from selected states for no (σ E = σ I = 0.0 mV/ms 3/2 , top panel) and finite (σ E = σ I = 0.05 mV/ms 3/2 , bottom panel) noise. Panels show (left to right) traces obtained from inside the limit cycle [µ E = 0.56 nA, µ I = 0.4 nA, marked with green square in (A)], from the left [µ E = 0.466 nA, µ I = 0.4 nA, marked with green star in (A)] and the right [µ E = 0.7 nA, µ I = 0.4 nA, marked with green circle in (A)] border of the limit cycle. For other parameters, see Supplementary Table S2. closing the feedback loop back to the cortex. Figure 5 shows the number of thalamic spindles as a function of the potassium leak conductance g LK and the conductance g h of TCR's rectifying current for constant external rate inputs to both thalamic populations. For larger external inputs, the regions of thalamic spindling translate to lower values of g h independent of FIGURE 5 | Effects of constant external input on thalamic spindling. The panels show the color-coded rate of thalamic spindles as a function of the conductance parameters g h and g LK for three different noise levels σ TCR and three different strengths r stim of the external input. External coupling strengths are equal for the TCR and TRN populations and are set to one. The other model parameters of the thalamic node are given in Supplementary Table S1. Squares (g h = 0.05 mS/cm 2 , g LK = 0.033 mS/cm 2 ) and triangles (g h = 0.062 mS/cm 2 , g LK = 0.033 mS/cm 2 ) mark parameter values used for Figure 6. Symbol color refers to the different thalamic noise levels. Note that the typical excitatory firing rate of the cortical node's UP state is up to 50-60 Hz.
TCR noise levels. Increased external input leads to a higher concentration of intracellular calcium in the TCR, which activates the rectifying I h current more strongly and ultimately abolishes the waxing and waning cycle of spindle oscillations (cf. Figure 2). For smaller values of g h , I h is reduced, and the waxing and waning of spindle oscillations reappear.
We then stimulate the thalamic node with rectangular current pulses mimicking an idealized sequence of cortical UP and DOWN states. Figure 6A shows the resulting rates of the TCR and TRN populations as a function of time. For higher values of g h , thalamic spindles are only induced in the DOWN state, as already expected from Figure 5. For lower values of g h , spindling resumes during UP states, albeit with lower spindle event frequency. Finally, we stimulate the thalamic node with the excitatory rate of a cortical node undergoing sustained adaptation-induced slow oscillations ( Figure 6B) and noise-induced DOWN-state transitions ( Figure 6C). With few exceptions, thalamic spindles are induced right after cortical UP to DOWN transitions, but-given the short duration of the cortical DOWN state-every SO induces one thalamic spindle FIGURE 6 | Thalamic response to alternating cortical UP and DOWN states. The (A-C) panels show the firing rates of the TCR (red), the TRN (blue), and the cortical input (black) as a function of time. Note that values for TCR and TRN are marked on the primary y-axis (left), while the values for cortical input are marked on the secondary y-axis (right) for clarity. Corticothalamic coupling strength (N ctx→thal = 1.0) was equal for the TCR and TRN populations. (A) Thalamic response to square wave stimulation of 0.05 Hz frequency with the amplitude representing the firing rate of an external population of 60 Hz. (B) Thalamic response to excitatory rate input, which is generated by a cortical node undergoing limit cycle oscillations (cf. Figure 4B, upper left panel). (C) Thalamic response to excitatory rate input, which is generated by a cortical node undergoing noise-induced DOWN state transitions (cf. Figure 4B, lower right panel). (D) Distribution of delays between midpoints of cortical DOWN states and spindle-band peaks in the noise-induced slow oscillations from panel (C). The parameters for the thalamic node were as follows: g h = 0.05 mS/cm 2 (left panels, see also green squares in Figure 5), g h = 0.062 mS/cm 2 (right panels, see also green triangles in Figure 5), g LK = 0.033 mS/cm 2 , for all other parameters see Supplementary Table S1. The parameters for the cortical node were as follows: µ E = 0.56 nA, σ E = σ I = 0.0 mV/ms 3/2 for panels (B) and µ E = 0.7 nA, σ E = σ I = 0.05 mV/ms 3/2 for panels (C), for all other parameters see Supplementary Table S2. The parameters of connected thalamocortical model are given in Supplementary Table S3. only. In this setting, this motif well reproduces the SO-correlated spindling observed during NREM sleep Ladenbauer et al., 2016Ladenbauer et al., , 2017Helfrich et al., 2018). Changing the conductance g h of TCR's rectifying current, the ratio between the numbers of "free" and SO-correlated spindle events can be adapted.

Cortical Model Driven by Thalamic Input
Secondly, we study the effects of external thalamic inputs to the cortical node with a one-way thalamus → cortex connection. Thalamic input (noise-free spindles) was simulated using our thalamic model parametrized in the spindling region I (Figure 3), i.e., with g LK = 0.018 mS/cm 2 , g h = 0.062 mS/cm 2 . Figure 7 summarizes cortical responses as a function of cortical parameters and thalamus → cortex connection strengths. We observe two distinct effects of thalamic stimulation. Firstly, it increases the slow oscillation frequency because additional thalamic input leads to an increase in the excitation, an effect similar to increasing µ E in the isolated cortical node.
The second effect is the imprinting of spindle oscillations into the cortical activity for stronger coupling strengths (Figure 7, lower two rows). Typically, a spontaneous thalamic spindle induces a transition to a cortical UP state, and spindle activity is superimposed on the UP state of the cortical activity. After the thalamic spindle ceases, the thalamic firing rate significantly decreases, and the cortex returns to the DOWN state due to the adaptation mechanism of excitatory neurons.
The addition of noise to the cortical node (σ E = σ I = 0.05 mV/ms 3/2 ) does not qualitatively change our previous observations. However, background noise can induce a DOWN swing (leading to irregular slow oscillation) in the cortex parametrized at the border between the limit cycle and the UP state (Supplementary Figure S2). Moreover, increasing the thalamus → cortex connection strength leads to prolonged UP states, to less frequent DOWN states, FIGURE 7 | Cortical response to thalamic spindle stimulation. The panels show the firing rates of the excitatory cortical population (red), the inhibitory cortical population (blue), and the thalamic input (black) as a function of time. Note that values for the excitatory and inhibitory cortical populations are marked on the primary y-axis (left), while the values for thalamic input are marked on the secondary y-axis (right) for clarity. Columns show different parametrizations of the cortical node (left to right): inside the slow limit cycle (µ E = 0.56 nA, cf. Figure 4B left), the border of the DOWN state and the limit cycle (µ E = 0.466 nA, cf. Figure 4B middle), and the border of the limit cycle and the UP state (µ E = 0.7 nA, cf. Figure 4B right). Different rows show different values of thalamus → cortex connection strength ranging from N thal→ctx = 0.0 to N thal→ctx = 0.1. The parameters for the cortical node were as follows: µ I = 0.4 nA, σ E = σ I = 0.0 mV/ms 3/2 , for all other parameters see Supplementary Table S2. The parameters for thalamic node were as follows: g LK = 0.018 mS/cm 2 , g h = 0.062 mS/cm 2 , σ TCR = 0.0 mV/ms 3/2 , for all other parameters see Supplementary Table S1. The parameters of connected thalamocortical model are given in Supplementary Table S3. and, for large enough coupling strength, to spindle activity from the thalamus being imprinted onto the UP state of the cortex.

Full Thalamocortical Loop
In the following section, we study the dynamics of the full thalamocortical loop consisting of one cortical and one thalamic node, coupled according to Figure 1. To summarize our hypotheses regarding the thalamocortical motif, we expect slow oscillations to emerge in the cortex and spindles in the thalamus (Krishnan et al., 2016;Latchoumane et al., 2017). Additionally, we expect slow oscillation activity to affect the timing of thalamic spindles (Hagler et al., 2018;Jiang et al., 2019) and thalamic spindles to cause cortical spindles during cortical UP states (Mölle et al., 2002;Ladenbauer et al., 2017;Helfrich et al., 2018). Figure 8 shows the amplitude difference and the power in two spectral bands for the firing rate of the cortical excitatory population, depending on the connection strengths of the thalamus → cortex and cortex → thalamus. Increasing thalamus → cortex connection strength leads to a broader region in which the cortical node is oscillating in the slow oscillation band, albeit the amplitude of cortical oscillation decreases ( Figure 8A). Moreover, increasing the thalamus → cortex connection strength also causes a power decrease in the slow oscillation band (0.1-3.0 Hz), and the region of increased slow oscillation power shifts slightly toward lower background excitation values (to the left in the panels), as depicted in Figure 8B. Finally, the power in the fast spindle band (12-15 Hz) increases as a function of both the thalamus → cortex and cortex → thalamus connection strengths (Figure 8C). Adding background noise to both cortical (σ E = σ I = 0.05 mV/ms 3/2 ) and thalamic (σ TCR = 0.005 mV/ms 3/2 ) nodes does not qualitatively change these observations (not shown). Thalamic parameterization was chosen from the spindling region II (Figure 3), i.e., the thalamus was simulated with g LK = 0.033 mS/cm 2 , g h = 0.062 mS/cm 2 . Figure 9A summarizes cortical spindle activity dependent on thalamus → cortex and cortex → thalamus connection strengths FIGURE 8 | Impact of changes of thalamocortical connection strengths on cortical SO and spindle activity. The panels show a single slice of the state space diagram of the cortical model on par with Figure 4, zoomed into the slow limit cycle region, for the recurrently connected thalamocortical model of Figure 1. The full model was simulated with different strengths of connections in both directions: thalamus → cortex is shown on the x-axis, and cortex → thalamus is shown on the y-axis. Panels show (A) the color-coded difference between the maximum and minimum firing rate in the excitatory cortical population, (B) the color-coded mean spectral power in the slow oscillation range (0.1-3 Hz), and (C) the color-coded mean spectral power in the fast spindle range (12-15 Hz). The mean spectral power in both bands was computed by averaging the power spectral density computed using Welch's method within the slow oscillation band, 0.1-3.0 Hz, and fast spindle band, 12-15 Hz. Note that the bottom left panels (N ctx→thal = 0.0, N thal→ctx = 0.0) parallel the bifurcation diagram of isolated cortex in Figure 4. All simulations were conducted without noise, for other parameters see Supplementary Tables S1-S3. and the level of background excitation and inhibition (µ E and µ I , respectively). The density of cortical spindles is directly proportional to the thalamus → cortex connection strength. Cortical spindle activity emerges mainly in the region where cortical node exhibits oscillatory activity (cf. Figure 8A).

Cortical Spindle Activity in the Thalamocortical Loop Model
The time series of firing rates in both nodes in the thalamocortical model is presented in Figure 9B. For selected values of model parameters, at the border between the cortical DOWN state and the limit cycle, a spontaneous thalamic spindle causes the cortical node to go into the UP state, and the cortical FIGURE 9 | Interaction between cortical slow oscillation and thalamic spindles in the noise-free case. (A) Estimated number of cortical spindles per second as a function of the mean external input currents µ ext to the excitatory and inhibitory population (encoded by x and y-axis in each smaller panel, respectively) and the thalamus → cortex and cortex → thalamus connectivity strength (encoded by x and y-axis over panels, respectively). The thalamocortical model was simulated for 65 s without noise. (B) 15 s excerpts of time series from the excitatory population in the cortical node (black) and TCR in the thalamic node (gray) for two sets of connection strengths (N thal→ctx = 0.04 and N ctx→thal = 0.4, denoted by lime green symbols; and N thal→ctx = 0.12 and N ctx→thal = 1.2, denoted by aqua blue symbols) and three values of µ E as per three dynamical states of the cortex: DOWN state with µ E = 0.45 nA (star), slow limit cycle with µ E = 0.55 nA (square), and UP state with µ E = 0.8 nA for weaker connection strengths, and µ E = 0.6 nA for stronger connection strength (circle). The points of interest are also denoted in (A). For all points µ I = 0.7 nA, and other parameters were kept constant as per Supplementary Tables S1-S3. activity is then shaped by the incoming spindle, depending on the thalamus → cortex connection strength as seen in Figure 9B (aqua blue and lime green star). In the slow limit cycle, cortical parametrization for the higher connection strengths, the cortical UP states, and thalamic spindle waxing become coupled ( Figure 9B, aqua blue and lime green square). Typically, the thalamic node elicits a spindle after cortical UP to DOWN state transition, which creates a thalamic DOWN state. This allows the thalamic node to go into hyperpolarization and generates a spindle.
The cortical UP state activity is slightly modulated by the spindle oscillations of the thalamic node (Figure 9B, lime green and aqua blue circles). With the addition of noise, this regime becomes interesting for our investigation, as seen in the next section.

Slow Oscillation-Spindle Interaction in the UP State Regime
In the UP state-dominant regime, the cortical model is parametrized in the UP state, with irregular DOWN swings caused by the background noise pushing the state into the limit cycle. The average length of the UP state depends on the amount µ E of background excitation.
For µ E = 0.61 nA, i.e., close to the bifurcation, UP states are of shorter duration (as shown in Figure 10). In this case, the UP and DOWN states are relatively regular, with cortical DOWN states providing a window of opportunity for hyperpolarization of the TCR and subsequent spindle generation, which then embed spindle oscillatory activity in the cortical node typically just after the DOWN to UP state transition. This behavior can be observed by following vertical white dashed lines in the cortical time series (Figure 10A), which denotes the midpoints of cortical DOWN states which are followed by a spindle within a 1.5 s window. These lines extend to the time-frequency representation, and we can see an intermittent increase of the cortical power in the spindle band, thus demonstrating cortical spindles nested in the UP states. The transition from the DOWN to the UP state in the cortical node is accompanied by an overshoot in the cortical activity, which then decreases due to the adaptation current, as shown in Figure 10B. The mechanistic explanation of the SOspindle relationship is illustrated by cortical DOWN-state and thalamic spindle time-locked plots (Figures 10B,D): the cortical DOWN state is followed by a waxing period in the thalamic activity, the spindle activity is then imprinted onto the cortical activity, followed by an approximately 1 second long cortical spindle. Thus, spindles are typically observed during the peak of the cortical UP state or shortly after, as shown in circular histograms in Figures 10B,D. Our model results align well with the data-driven studies Helfrich et al., 2018). Mechanistically, a sustained UP state in the cortex FIGURE 10 | Spindles in the thalamocortical motif with short UP states. The figure shows various variables from 120 s of simulation of full thalamocortical model simulation with the cortical node being parametrized at the right border between the slow limit cycle and the UP state, relatively close to the bifurcation line (µ E = 0.61 nA). Individual panels show the following: (A) 20 s time series excerpt of the firing rate of the excitatory population in the cortical node (black) and its slow oscillation phase (blue) computed using the Hilbert transform of low-pass filtered cortical excitatory firing rate. The panel below shows the time-frequency representation of the cortical excitatory firing rates computed using the Short Time Fourier Transform with a 2-s time window. Dashed vertical lines (red in time series plot and white in time-frequency plot) denote the midpoints of cortical DOWN states which are followed by a cortical spindle in the 1.5-s window. (B) the mean ± the SEM of the cortical excitatory firing rates (red) and TCR firing rates (black), locked on the cortical DOWN states for the whole interval of 120 s. The panel below shows the distribution of cortical slow oscillation phases (shown in (A) with blue) for the maximum of the cortical fast spindle band peak. Shown are histogram bars in thin red and circular mean (± circular STD) in thick (dashed) red. (C) 20 s time series excerpt of the firing rate of the thalamocortical relay population in the thalamic node (black). The panel below shows the time-frequency representation of the TCR firing rates computed using the Short Time Fourier Transform with a 2-s time window. Dashed vertical lines (red in time series plot and white in time-frequency plot) denote the midpoints of cortical DOWN states which are followed by a cortical spindle in the 1.5-s window. (D) the mean ± the SEM of the TCR firing rates (red) and cortical excitatory firing rates (black), locked on the thalamic spindle peaks in the whole interval of 120 s. The panels below show as follows: the distribution of delays between midpoints of cortical DOWN states and spindle-band peaks and the distribution of cortical slow oscillation phases for the maximum of the thalamic fast spindle band peak. Shown are histogram bars in thin red and circular mean (± circular STD) in thick (dashed) red. The model was simulated with N thal→ctx = 0.12, N ctx→thal = 1.2, µ I = 0.4 nA, σ E = σ I = 0.05 mV/ms 3/2 , and σ TCR = 0.005 mV/ms 3/2 , while other parameters were kept constant as per Supplementary Tables S1-S3.
suppresses spindling via the corticothalamic connections, since TCR cannot reach necessary hyperpolarization. Therefore, for this particular parametrization (predominant cortical UP states with irregular DOWN swings), we conclude the causal direction goes in the sense of the cortex → thalamus.
For a parametrization further away from the bifurcation (µ E = 0.66 nA), prolonged UP states are observed (Supplementary Figure S3). The SO-spindle interactions are similar as to the case with shorter UP states (Figure 10), however, these prolonged UP states are more similar to what is seen during human NREM sleep Helfrich et al., 2018).
We also probed our connected thalamocortical model in the remaining two cortical parametrizations. In particular, when the cortex is parametrized in the DOWN state (with µ E = 0.36 nA, other parameters unchanged with respect to previous text), the causal pathway between cortex and thalamus reverses its direction. In this case, the cortex is predominantly in the DOWN state. Without a sustained cortical drive, the thalamus can generate free spindles, and these spindles, in turn, drive the cortex into the UP state as they are projected (Supplementary Figure S4). Finally, when the cortex is parametrized in the limit cycle (µ E = 0.42 nA), its UP states hinder the generation of free spindles in the thalamus, and the thalamic spindles are generated only in prolonged DOWN states. When the thalamus generates a spindle within the window of the cortical DOWN state, the projection of the spindle then forces the cortex to the DOWN state-UP state transition (Supplementary Figure S5).
Following previous data-driven results on the slow wavespindle nesting (Sirota et al., 2003;Ladenbauer et al., 2017;Helfrich et al., 2018), we finally quantified the phase-phase and phase-amplitude coupling between these two rhythms. Figure 11 shows the distribution of average amplitudes of thalamic spindles as a function of the phase of the cortical slow oscillation. Spindle amplitudes are highest between slow oscillation phases 0 and π/4. The phase-amplitude relationship was indeed significant as shown by both tested measures (KL − MI = 0.0109, p < 0.001 and MVL = 11.3136, p < 0.001), while the phase-phase relationship cannot be deemed significant (PLV = 0.0020, p = 0.372 and MI = 0.0049, p = 0.266). We conclude that the thalamocortical model possesses significant phase-amplitude coupling between the phase of slow oscillations and the amplitude of thalamic spindles, but no phase-phase coupling was found.

DISCUSSION
This study investigated the dynamical states of a biophysically realistic neural mass model of a thalamocortical motif. To assess the contribution of each part of the model to the dynamics, we started with an isolated cortical and thalamic node and examined their individual dynamical landscapes. We perturbed each node with an external stimulus resembling their counterpart, i.e., the cortical model was perturbed by stimulation with spindle-like oscillations, and the thalamic model was perturbed using a square pulse with a low frequency, mimicking the idealized UP and DOWN state sequence of cortical SO activity. Next, we connected the isolated nodes to a full thalamocortical network and focused on spindle imprinting onto the cortical UP state activity and the interactions between spindles and slow oscillations, both being hallmark activity in the human brain during slow-wave sleep.
The results of our modeling study are in line with previous neuroimaging studies focusing on various aspects of sleep spindles and their interactions with cortical slow oscillations. In particular, the frequency of the spindles in our model (as shown in Figures 10A,C, Supplementary Figures S3A,C) both in thalamic and cortical nodes matches the experimental values for fast spindles with an average of around 12 Hz (Nir et al., 2011;Purcell et al., 2017;Alfonsi et al., 2019;Ujma et al., 2021). Typical spindle duration between 0.5-1 s also agrees with experimental values (Purcell et al., 2017). Finally, recent neuroimaging studies showed that sleep spindles possess a significant phase-amplitude relationship with cortical slow oscillations. More concretely, waxing periods of spindles are nested in the cortical UP states and the spindle peak typically occurs right at (or just after) the slow oscillation UP state peak Helfrich et al., 2018Helfrich et al., , 2019. The presented thalamocortical model well reproduces this SO-spindle relationship (as shown in Figures 10B,D circular histograms and Figure 11).
By dissecting the thalamocortical model, we sought the mechanistic explanation of the observed SO-spindle relationship. From the modeling perspective, we needed to recognize the difference between two different cortical parametrizations and their mechanistic consequences. The causal direction follows the path cortex → thalamus in the cortical parametrization with predominantly UP state activity with irregular, noise-driven DOWN state excursions. A sustained UP state in the cortex suppresses spindling via the corticothalamic connections since the TCR population cannot reach the necessary hyperpolarization (cf. Figure 2). On irregular cortical DOWN swings caused by the background noise, the cortex is pushed into the DOWN state, which creates a window of opportunity for the thalamic DOWN state and a subsequent spindle waxing period. Finally, this spindle is then imprinted onto the cortical activity due to the thalamocortical projections. Reversely, with DOWN state-dominant cortical activity with irregular UP state excursions (cf. Figure 4B middle), the causal direction reverses. In this case, without sustained cortical drive, the thalamus can generate "free" spindles, which subsequently push cortical activity into the UP state.
We left some parameters untouched despite exploring the model dynamics with various parameter settings. In particular, the conductances of calcium T-type currents g T in both thalamic populations control the underlying frequency of spindles in the spindle band (cf. Figure 2). Allowing heterogeneous connections from the excitatory population in the cortical node to TCR and TRN populations (as shown in Figure 1) would control the balance between cortically driven excitation and inhibition, leading to an altered shape of the thalamic spindles. In particular, for a higher TRN/TCR input ratio, i.e., higher strength for cortex → TRN connection, the waxing periods are short-individual spindles contain only one or two oscillations. The change in shape is visible for higher ratios than 10 (not shown), and thus we have not included this parameter in our overall investigation. Note that the underlying frequency in the spindle band is not modulated by heterogeneous connection strengths between the cortex and thalamus. Finally, one of the caveats of our study is the unknown realistic connection strengths in a thalamocortical motif in humans, which is why we treated both connection strengths from cortex → thalamus and thalamus → cortex as free parameters of our investigation. Given the maximum firing rates of both populations (cortex ∼40 Hz, thalamus ∼400 Hz), we indeed expected the connectivity strengths in both directions to have approximately a 10-fold difference.
The dynamics of slow oscillations and spindles vary significantly between distinct sleep stages (Amzica and Steriade, 1997;Mölle et al., 2002;Steriade, 2003). Various neuromodulators, such as acetylcholine (ACh) or histamine (HA), are known to vary significantly during sleep and awake, as well as across sleep stages (Brady et al., 2011;Vanini et al., 2011). Specific effects of these neuromodulators can be implemented by changing the strength of the intrinsic and synaptic currents in the cortical and thalamic populations. In their biophysically realistic thalamocortical model, Krishnan et al. (2016) identified a minimal set sufficient to account for characteristic changes in the brain's electrical activity across the sleep-wake cycle. We have not included these effects and focused solely on the N3 sleep stage. Briefly, a reduction of ACh can be implemented by an increase in potassium leak conductance (g LK ) in the thalamic node and an increase in spike-frequency adaptation (b) in the cortical node (McCormick, 1992), while the effect of HA can be implemented as a shift in the activation curve of a hyperpolarization-activated current I h in the TCR (McCormick and Williamson, 1991).
This study merged two mass modeling approaches: a meanfield approximation of the Fokker-Planck equations for the cortical dynamics (Cakan and Obermayer, 2020) and a mass model based on conductance-based average membrane voltage dynamics for the thalamus (Schellenberger Costa et al., 2016), and as such, showcase the applicability of the hybrid modeling approaches. We are not aware of any caveats of this approach, as long as a natural link in the sense of coupling variable between the frameworks can be established. In our case, both frameworks operate with the notion of firing rate and use firing rate as the main model output and a coupling variable between their subpopulations. Using the already probed cortical node facilitates studying biophysically realistic stimulation protocols of the thalamocortical model in our future work. Our motivation was 2-fold: the cortical model introduced by Cakan and Obermayer (2020), which we used here, is biophysically realistic and very well studied with respect to its dynamical states. Moreover, a previous study also probed biophysically realistic stimulation protocols on par with transcranial direct current stimulation (tDCS) during sleep (Cakan and Obermayer, 2020). Finally, by merging different mass modeling approaches in our current work, we are also setting a stage for future extensibility of the thalamocortical model by including more cortical nodes, or, alternatively, nodes simulating different brain areas, e.g., the hippocampus.
One of the extension possibilities lies in topographic mass models. A well-known property of cortical spindles is their heterogeneity, i.e., mainly the difference between faster and slower spindles, where faster spindles are usually observed in parietal regions, while slower spindles are found in the frontal regions (Werth et al., 1997;Mölle et al., 2002. The main reason for this heterogeneity is the core and matrix thalamocortical pathways (Rubio-Garrido et al., 2009;Piantoni et al., 2016): core thalamocortical neurons are spatially selective and topographically organized, target a single cortical area, and project mainly to the granular layer. On the other hand, matrix neurons have diffuse, multiarea projections, characterized by multiple distant arbors, and reach mostly superficial layers of the cortex (Rubio-Garrido et al., 2009;Piantoni et al., 2016). The advantage of mass models over spiking models is the computational efficiency and the fast transition into layerresolving 2D models by modeling more cortical layers with the same base model and a different between-layer corticocortical and corticothalamic connectivity and fan-outs (Potjans and Diesmann, 2014), accounting for distinct core and matrix thalamocortical projections.
As already mentioned above, a realistic outlook and extension possibility is taking a step further in the direction of whole-brain dynamics and creating a whole-brain model consisting of many cortical nodes coupled using the structural connectome with the thalamus, able to undergo stimulation. Recently, Cakan et al. (2022) constructed a deep sleep whole-brain model and studied the dynamics of local and global slow oscillation events. Each node in their model consisted of an excitatory-inhibitory pair of the mean-field approximation model of the AdEx neurons. Hence, we might build on our current investigation of the thalamocortical motif.
We suggest that our current investigation of the computationally efficient thalamocortical microcircuit model allows us to dive deeper into the sleeping brain and shed light on the exact temporal structure and interaction of sleep rhythms involved in episodic memory consolidation.