The architecture of the thalamocortical model is shown in Figure 1. It consists of two thalamic neural populations, namely 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 spike-frequency 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).

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 cation-nonselective 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 rk′ of a presynaptic population k′ elicits a postsynaptic response s_lk in population k given by

Nkk′ 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.

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

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 voltage-dependent 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 by

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 mean V-α of the membrane potentials, the adaptive timescale τ_α, and the population firing rate r_α in the steady-state 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 synapses s-αβ obeys

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 σs,αβ2 are given by

The parameters are summarized in Supplementary Table S2.

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 Nkk′=Nctx→thal and rk′=rE(t-dctx,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.

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., 2021). neurolib is a computational framework for whole-brain 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) satisfying

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).

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.

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.

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 mechanisms—one 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 spindle-promoting 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², the TCR exhibit down state with a constant rate of 10 Hz, while for g_h = g_LK = 0 mS/cm² 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 T-type 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.

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 (Cakan et al., 2022).

We first investigate the state space of the thalamic model by studying the effects of external (cortical) inputs without 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 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 only. In this setting, this motif well reproduces the SO-correlated spindling observed during NREM sleep (Mölle et al., 2011; Ladenbauer et al., 2016, 2017; Helfrich 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.

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², g_h = 0.062 mS/cm². 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, and, for large enough coupling strength, to spindle activity from the thalamus being imprinted onto the UP state of the cortex.

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², g_h = 0.062 mS/cm².

Figure 9A summarizes cortical spindle activity dependent on thalamus → cortex and cortex → thalamus connection strengths 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).

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 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.

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 SO–spindle 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 (Ladenbauer et al., 2017; Helfrich et al., 2018). Mechanistically, a sustained UP state in the cortex 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 (Ladenbauer et al., 2017; 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 wave–spindle nesting (Sirota et al., 2003; Mölle et al., 2011; 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.