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.

1. 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 low-amplitude 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; Mölle et al., 2011). Sleep spindles are also observed in the deeper, N3 sleep stage, albeit with lower power in the fast spindle frequency range (Cox et al., 2017) 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; Mölle et al., 2011; 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 (Mölle and Born, 2011; Mölle et al., 2011).

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-wave-generating 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., 2016, 2017; Helfrich 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.

2. Materials and Methods

2.1. Model Design

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.

FIGURE 1

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.

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

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

$\begin{array}{ll}\tau {\stackrel{.}{V}}_{\alpha }=-(V_{\alpha }-E_{\text{leak}})-w_e{s}_e(V_{\alpha }-E_e)-w_i{s}_i(V_{\alpha }-E_i)-C_m^{-1}\tau I^{\text{intrinsic}},& (1)\end{array}$

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,

$\begin{array}{ll}r(V)=\frac{r_{max}}{1+exp\left(-(V-\theta )/\sigma \right)},& (2)\end{array}$

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,

$\begin{array}{ll}{I}_{LK}=g_{LK}(V-E_K),& (3)\end{array}$

and a T-type calcium current,

$\begin{array}{ll}{I}_T=g_T{m}_{\infty }^{2}h(V-E_{Ca}),& (4)\end{array}$

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_{\infty }^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:

$\begin{array}{ll}{I}_h=g_h(m_{h1}+g_{inc}{m}_{h2})(V-E_h),& (5)\end{array}$

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

FIGURE 2

Figure 2. Schematic of the spindle generation mechanisms. Refer to the text for a detailed explanation. Thalamic connectivity is shown in the upper left panel. The figure was adapted from Mayer et al. (2006).

Synaptic transmission in the thalamic model is conveyed by conductance-based synapses, where the spike rate $r_{k^{\prime }}$ of a presynaptic population k′ elicits a postsynaptic response s_lk in population k given by

$\begin{array}{ll}{s}_{lk}(t)={\displaystyle \sum _{k^{\prime }}}{\alpha }_l(t)\otimes N_{k{k}^{\prime }}{r}_{k^{\prime }}(t)+{\varphi }^{\prime }(t).& (6)\end{array}$

$N_{k{k}^{\prime }}$ 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

$\begin{array}{ll}{\alpha }_l(t)={\gamma }_l^2·t·\text{exp}(-{\gamma }_{l}t)& (7)\end{array}$

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)

$\begin{array}{ll}{\ddot{s}}_{lk}={\gamma }_l^2\left({\displaystyle \sum _{k^{\prime }}}{N}_{k{k}^{\prime }}{r}_{k^{\prime }}(t)+{\varphi }^{\prime }(t)-s_{lk}(t)\right)-2{\gamma }_l{ṡ}_{lk}(t).& (8)\end{array}$

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

$\begin{array}{l}{\tau }_t{\stackrel{.}{V}}_t=-(V_t-E^{\text{leak}})-w_e{s}_{et}(V_t-E^e)-w_i{s}_{it}(V_t-E^i)\\ -C_m^{-1}{\tau }_t\left(I_t^{LK}+I_t^T+I_h\right)& (9)\end{array}$

$\begin{array}{l}{\tau }_r{\stackrel{.}{V}}_r=-(V_r-E^{\text{leak}})-w_e{s}_{er}(V_r-E^e)-w_i{s}_{ir}(V_r-E^i)\\ -C_m^{-1}{\tau }_r\left(I_r^{LK}+I_r^T\right),& (10)\end{array}$

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.

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

$\begin{array}{ll}{C}_{m,c}{\stackrel{.}{V}}_j=I_{\text{ion}}(V_j)+I_{\text{syn}}(t)+I_{j,\text{ext}}(t)& (11)\end{array}$

$\begin{array}{ll}{I}_{\text{ion}}(V)=g_L(E_L-V)+g_L{\Delta }_T\text{exp}\left(\frac{V-V_T}{{\Delta }_T}\right)-I_{A,j}(t)& (12)\end{array}$

$\begin{array}{ll}{\tau }_A{İ}_{A,j}=a(V_j-E_A)-I_{A,j}.& (13)\end{array}$

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

$\begin{array}{ll}{I}_{\text{syn}}(t)=C_{m,c}\left(J_{\alpha E}{s}_{i,\alpha E}(t)+J_{\alpha I}{s}_{i,\alpha I}(t)\right),& (14)\end{array}$

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

$\begin{array}{ll}{ṡ}_{i,\alpha \beta }=-\frac{s_{i,\alpha \beta }}{{\tau }_s}+\frac{c_{\alpha \beta }}{J_{\alpha \beta }}\left(1-s_{i,\alpha \beta }\right){\displaystyle \sum _j}{G}_{ij}{\displaystyle \sum _k}\delta (t-t_k^j-d_{\alpha \beta }),& (15)\end{array}$

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 ${\stackrel{-}{V}}_{\alpha }$ 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.,

$\begin{array}{ll}{r}_{\alpha }={\Phi }_r({\mu }_{\alpha },{\sigma }_{\alpha }),& (16)\end{array}$

$\begin{array}{ll}{\stackrel{-}{V}}_{\alpha }={\Phi }_V({\mu }_{\alpha },{\sigma }_{\alpha }),& (17)\end{array}$

$\begin{array}{ll}{\tau }_{\alpha }={\Phi }_{\tau }({\mu }_{\alpha },{\sigma }_{\alpha }),& (18)\end{array}$

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

$\begin{array}{ll}{\tau }_{\alpha }{\stackrel{.}{\mu }}_{\alpha }={\mu }_{\alpha }^{\text{syn}}+{\mu }_{\alpha }^{\text{ext}}-{\mu }_{\alpha },& (19)\end{array}$

$\begin{array}{ll}{\mu }_{\alpha }^{\text{syn}}=J_{\alpha E}{\stackrel{-}{s}}_{\alpha E}+J_{\alpha I}{\stackrel{-}{s}}_{\alpha I},& (20)\end{array}$

$\begin{array}{ll}{\stackrel{.}{Ī}}_A={\tau }_A^{-1}\left(a({\stackrel{-}{V}}_E-E_A)-{Ī}_A\right)+b·r_E,& (21)\end{array}$

$\begin{array}{ll}{\mu }_{\alpha }^{\text{ext}}=\frac{{\mu }_{ctx,\alpha }-{\mu }_{\alpha }^{\text{ext}}}{{\tau }_{OU}}+{\sigma }_{\alpha }{\xi }_{\alpha },& (22)\end{array}$

where the dynamics of mean membrane current μ_α depends on the synaptic current ${\mu }_{\alpha }^{\text{syn}}$ and an external noisy current ${\mu }_{\alpha }^{\text{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 ${\stackrel{-}{s}}_{\alpha \beta }$ obeys

$\begin{array}{ll}{\stackrel{.}{\stackrel{-}{s}}}_{\alpha \beta }=-{\tau }_{s,\beta }^{-1}{\stackrel{-}{s}}_{\alpha \beta }+(1-{\stackrel{-}{s}}_{\alpha \beta })·r_{\alpha \beta },& (23)\end{array}$

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

$\begin{array}{ll}{r}_{\alpha \beta }=\frac{c_{\alpha \beta }}{J_{\alpha \beta }}{K}_{\beta }·r_{\beta }(t-d_{\alpha \beta }),& (24)\end{array}$

$\begin{array}{ll}{\rho }_{\alpha \beta }=\frac{c_{\alpha \beta }}{J_{\alpha \beta }}·r_{\alpha \beta }.& (25)\end{array}$

Finally, the current variance ${\sigma }_{\alpha }^2$ and the variance of the fraction of active synapses ${\sigma }_{s,\alpha \beta }^2$ are given by

$\begin{array}{ll}{\sigma }_{\alpha }^2={\displaystyle \sum _{\beta \in \left\{E,I\right\}}}\frac{2J_{\alpha \beta }^2{\sigma }_{s,\alpha \beta }^2{\tau }_{s,\beta }{\tau }_m}{(1+r_{\alpha \beta }){\tau }_m+{\tau }_{s,\beta }}+{\sigma }_{\text{ext},\alpha }^2,& (26)\end{array}$

$\begin{array}{l}{\stackrel{.}{\sigma }}_{s,\alpha \beta }^2={(1-{\stackrel{-}{s}}_{\alpha \beta })}^2·{\rho }_{\alpha \beta }\\ +{\tau }_{s,\beta }^{-2}\left({\rho }_{\alpha \beta }-2{\tau }_{s,\beta }({\rho }_{\alpha \beta }+1)\right)·{\sigma }_{s,\alpha \beta }^2.& (27)\end{array}$

The parameters are summarized in Supplementary Table S2.

2.4. 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_{k{k}^{\prime }}=N_{ctx\to thal}$ and $r_{k^{\prime }}=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.

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

$\begin{array}{ll}\text{d}x=\frac{{\mu }_{\alpha }-x}{\tau }\text{d}t+{\sigma }_{\alpha }\xi ,& (28)\end{array}$

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 ${\mu }_{\alpha }^{\text{ext}}$ into Equation (22).

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

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

$\begin{array}{ll}\varphi (t)=\text{arctan}\frac{ŝ(t)}{s(t)}& (29)\end{array}$

$\begin{array}{ll}A(t)=\sqrt{s^2(t)+{ŝ}^2(t)}.& (30)\end{array}$

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

$\begin{array}{ll}KL\text{–}MI=\frac{log(N)-H(P)}{log(N)},& (31)\end{array}$

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

$\begin{array}{ll}MVL={\displaystyle \sum _{t=0}^T}{A}_{spindle}(t)e^{i{\varphi }_{SO}(t)}.& (32)\end{array}$

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:

$\begin{array}{ll}PLV={\displaystyle \sum _{t=0}^T}{e}^{i({\varphi }_{SO}(t)-{\varphi }_{spindle}(t))}.& (33)\end{array}$

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

$\begin{array}{ll}MI={\displaystyle \sum _{y\in Y}}{\displaystyle \sum _{x\in X}}p(x,y)log\left(\frac{p(x,y)}{p(x)p(y)}\right),& (34)\end{array}$

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.

3. Results

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

FIGURE 3

Figure 3. The generation of thalamic sleep spindles depends on the conductances g_LK and g_h. (A) 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² (I), g_LK = 0.031 mS/cm² (II), g_LK = 0.024 mS/cm² (III), and g_h = 0.062 mS/cm² (all three regions). For other model parameters, see Supplementary Table S1.

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.

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

FIGURE 4

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.

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

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

FIGURE 5

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², g_LK = 0.033 mS/cm²) and triangles (g_h = 0.062 mS/cm², g_LK = 0.033 mS/cm²) 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.

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.

FIGURE 6

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² (left panels, see also green squares in Figure 5), g_h = 0.062 mS/cm² (right panels, see also green triangles in Figure 5), g_LK = 0.033 mS/cm², 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.

3.4. 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², 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.

FIGURE 7

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

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.

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

FIGURE 8

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.

3.6. Cortical Spindle Activity in the Thalamocortical Loop Model

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

FIGURE 9

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.

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.

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

FIGURE 10

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.

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.

FIGURE 11

Figure 11. Phase-amplitude CFC in the thalamocortical motif. The figure shows the distribution of the average amplitudes of the thalamic spindle dependent on the phase of the cortical slow oscillations in the simulation with short cortical UP states (time series from Figure 10 with μ_E = 0.61 nA). Dark gray bars denote the distribution in the simulated thalamocortical model, while light gray bars denote the mean and 95^th percentile distribution computed using 1000 IAAFT surrogates (for details, see Cross-Frequency Coupling (CFC) Measures).

4. 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 (Ladenbauer et al., 2017; Helfrich et al., 2018, 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 mean-field 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, 2011). 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 layer-resolving 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.

Data Availability Statement

The original contributions presented in the study are publicly available. This data can be found here: GitHub, https://github.com/jajcayn/thalamocortical_model_study.

Author Contributions

NJ and KO conceived and designed the study. NJ implemented the thalamocortical model, ran numerical simulations, and data analyses. NJ, CC, and KO discussed the results and wrote the article. All authors contributed to the article and approved the submitted version.

Funding

NJ was funded by the Operational Programme Research, Development and Education, Ministry of Education, Youth and Sport of the Czech Republic (co-funded by the EU)—Project Number CZ.02.2.69/0.0/0.0/19_074/0016209 (Modelling the sleeping brain: toward a neural mass model of sleep rhythms and their interactions) and by the Czech Science Foundation project number 21-32608S. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project number 327654276–SFB 1315 (CC and KO).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Some of the computational resources used in this study were supplied by the project “e-Infrastruktura CZ” (e-INFRA CZ LM2018140) supported by the Ministry of Education, Youth and Sports of the Czech Republic. This manuscript has been released as a pre-print at BioRxiv (Jajcay et al., 2022).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.769860/full#supplementary-material

References

Achermann, P., and Borbely, A. (1997). Low-frequency (< 1 Hz) oscillations in the human sleep electroencephalogram. Neuroscience 81, 213–222. doi: 10.1016/S0306-4522(97)00186-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Alfonsi, V., D'Atri, A., Gorgoni, M., Scarpelli, S., Mangiaruga, A., Ferrara, M., et al. (2019). Spatiotemporal dynamics of sleep spindle sources across NREM sleep cycles. Front. Neurosci. 13, 727. doi: 10.3389/fnins.2019.00727

PubMed Abstract | CrossRef Full Text | Google Scholar

Amzica, F., and Steriade, M. (1997). The K-complex: its slow (< 1-Hz) rhythmicity and relation to delta waves. Neurology 49, 952–959. doi: 10.1212/WNL.49.4.952

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, M., Ladenbauer, J., Baumann, F., and Obermayer, K. (2017). Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: comparison and implementation. PLoS Comput. Biol. 13, e1005545. doi: 10.1371/journal.pcbi.1005545

PubMed Abstract | CrossRef Full Text | Google Scholar

Axmacher, N., Elger, C. E., and Fell, J. (2008). Ripples in the medial temporal lobe are relevant for human memory consolidation. Brain 131, 1806–1817. doi: 10.1093/brain/awn103

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazhenov, M., Timofeev, I., Steriade, M., and Sejnowski, T. J. (2002). Model of thalamocortical slow-wave sleep oscillations and transitions to activated states. J. Neurosci. 22, 8691–8704. doi: 10.1523/JNEUROSCI.22-19-08691.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bendor, D., and Wilson, M. A. (2012). Biasing the content of hippocampal replay during sleep. Nat. Neurosci. 15, 1439–1444. doi: 10.1038/nn.3203

PubMed Abstract | CrossRef Full Text | Google Scholar

Berry, R. B., Brooks, R., Gamaldo, C. E., Harding, S. M., Marcus, C., Vaughn, B. V., et al. (2012). “The AASM manual for the scoring of sleep and associated events,” in Rules, Terminology and Technical Specifications, Vol. 176 (Darien, IL: American Academy of Sleep Medicine), 2012.

PubMed Abstract | Google Scholar

Bibbona, E., Panfilo, G., and Tavella, P. (2008). The ornstein-uhlenbeck process as a model of a low pass filtered white noise. Metrologia 45, S117. doi: 10.1088/0026-1394/45/6/S17

CrossRef Full Text | Google Scholar

Bonjean, M., Baker, T., Bazhenov, M., Cash, S., Halgren, E., and Sejnowski, T. (2012). Interactions between core and matrix thalamocortical projections in human sleep spindle synchronization. J. Neurosci. 32, 5250–5263. doi: 10.1523/JNEUROSCI.6141-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonjean, M., Baker, T., Lemieux, M., Timofeev, I., Sejnowski, T., and Bazhenov, M. (2011). Corticothalamic feedback controls sleep spindle duration in vivo. J. Neurosci. 31, 9124–9134. doi: 10.1523/JNEUROSCI.0077-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Brady, S., Siegel, G., Albers, R. W., and Price, D. (2011). Basic Neurochemistry: Principles of Molecular, Cellular, and Medical Neurobiology. Cambridge, MA: Academic Press.

Google Scholar

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neurosci. 8, 183–208. doi: 10.1023/A:1008925309027

PubMed Abstract | CrossRef Full Text | Google Scholar

Cakan, C., Dimulescu, C., Khakimova, L., Obst, D., Flöel, A., and Obermayer, K. (2022). Spatiotemporal patterns of adaptation-induced slow oscillations in a whole-brain model of slow-wave sleep. Front. Comput. Neurosci. 15, 129. doi: 10.3389/fncom.2021.800101

PubMed Abstract | CrossRef Full Text | Google Scholar

Cakan, C., Jajcay, N., and Obermayer, K. (2021). neurolib: A simulation framework for whole-brain neural mass modeling. Cogn. Comput. doi: 10.1007/s12559-021-09931-9

CrossRef Full Text | Google Scholar

Cakan, C., and Obermayer, K. (2020). Biophysically grounded mean-field models of neural populations under electrical stimulation. PLoS Comput. Biol. 16, e1007822. doi: 10.1371/journal.pcbi.1007822

PubMed Abstract | CrossRef Full Text | Google Scholar

Canolty, R. T., Edwards, E., Dalal, S. S., Soltani, M., Nagarajan, S. S., Kirsch, H. E., et al. (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science 313, 1626–1628. doi: 10.1126/science.1128115

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, M. X. (2008). Assessing transient cross-frequency coupling in eeg data. J. Neurosci. Methods 168, 494–499. doi: 10.1016/j.jneumeth.2007.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cona, F., Lacanna, M., and Ursino, M. (2014). A thalamo-cortical neural mass model for the simulation of brain rhythms during sleep. J. Comput. Neurosci. 37, 125–148. doi: 10.1007/s10827-013-0493-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Contreras, D., Destexhe, A., Sejnowski, T. J., and Steriade, M. (1997). Spatiotemporal patterns of spindle oscillations in cortex and thalamus. J. Neurosci. 17, 1179–1196. doi: 10.1523/JNEUROSCI.17-03-01179.1997

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, R., Hofman, W. F., and Talamini, L. M. (2012). Involvement of spindles in memory consolidation is slow wave sleep-specific. Learn. Mem. 19, 264–267. doi: 10.1101/lm.026252.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, R., Schapiro, A. C., Manoach, D. S., and Stickgold, R. (2017). Individual differences in frequency and topography of slow and fast sleep spindles. Front. Hum. Neurosci. 11, 433. doi: 10.3389/fnhum.2017.00433

PubMed Abstract | CrossRef Full Text | Google Scholar

De Gennaro, L., and Ferrara, M. (2003). Sleep spindles: an overview. Sleep Med. Rev. 7, 423–440. doi: 10.1053/smrv.2002.0252

PubMed Abstract | CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., Robinson, P. A., Breakspear, M., and Friston, K. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4, e1000092. doi: 10.1371/journal.pcbi.1000092

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Bal, T., McCormick, D. A., and Sejnowski, T. J. (1996a). Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070. doi: 10.1152/jn.1996.76.3.2049

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Contreras, D., Steriade, M., Sejnowski, T. J., and Huguenard, J. R. (1996b). In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. J. Neurosci. 16, 169–185. doi: 10.1523/JNEUROSCI.16-01-00169.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Neubig, M., Ulrich, D., and Huguenard, J. (1998). Dendritic low-threshold calcium currents in thalamic relay cells. J. Neurosci. 18, 3574–3588. doi: 10.1523/JNEUROSCI.18-10-03574.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., and Sejnowski, T. J. (2003). Interactions between membrane conductances underlying thalamocortical slow-wave oscillations. Physiol. Rev. 83, 1401–1453. doi: 10.1152/physrev.00012.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandez, L. M. J., and Lüthi, A. (2020). Sleep spindles: mechanisms and functions. Physiol. Rev. 100, 805–868. doi: 10.1152/physrev.00042.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., et al. (2013). MEG and EEG data analysis with MNE-Python. Front. Neurosci. 7, 267. doi: 10.3389/fnins.2013.00267

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagler, D. J., Ulbert, I., Wittner, L., Erőss, L., Madsen, J. R., Devinsky, O., et al. (2018). Heterogeneous origins of human sleep spindles in different cortical layers. J. Neurosci. 38, 3013–3025. doi: 10.1523/JNEUROSCI.2241-17.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Helfrich, R. F., Lendner, J. D., Mander, B. A., Guillen, H., Paff, M., Mnatsakanyan, L., et al. (2019). Bidirectional prefrontal-hippocampal dynamics organize information transfer during sleep in humans. Nat. Commun. 10, 1–16. doi: 10.1038/s41467-019-11444-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Helfrich, R. F., Mander, B. A., Jagust, W. J., Knight, R. T., and Walker, M. P. (2018). Old brains come uncoupled in sleep: slow wave-spindle synchrony, brain atrophy, and forgetting. Neuron 97, 221–230. doi: 10.1016/j.neuron.2017.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Jajcay, N., Cakan, C., and Obermayer, K. (2022). Cross-frequency slow oscillation-spindle coupling in a biophysically realistic thalamocortical neural mass model. bioRxiv [Preprint]. doi: 10.1101/2021.08.29.458101.

CrossRef Full Text | Google Scholar

Jajcay, N., Kravtsov, S., Sugihara, G., Tsonis, A. A., and Paluš, M. (2018). Synchronization and causality across time scales in El Niño Southern Oscillation. npj Clim. Atmos. Sci. 1, 1–8. doi: 10.1038/s41612-018-0043-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, D., and Wilson, M. A. (2007). Coordinated memory replay in the visual cortex and hippocampus during sleep. Nat. Neurosci. 10, 100–107. doi: 10.1038/nn1825

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, X., Gonzalez-Martinez, J., and Halgren, E. (2019). Posterior hippocampal spindle ripples co-occur with neocortical theta bursts and downstates-upstates, and phase-lock with parietal spindles during NREM sleep in humans. J. Neurosci. 39, 8949–8968. doi: 10.1523/JNEUROSCI.2858-18.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, E. G. (2001). The thalamic matrix and thalamocortical synchrony. Trends Neurosci. 24, 595–601. doi: 10.1016/S0166-2236(00)01922-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, U., Bal, T., and McCormick, D. A. (1995). Spindle waves are propagating synchronized oscillations in the ferret LGNd in vitro. J. Neurophysiol. 74, 1301–1323. doi: 10.1152/jn.1995.74.3.1301

PubMed Abstract | CrossRef Full Text | Google Scholar

King, B. R., Hoedlmoser, K., Hirschauer, F., Dolfen, N., and Albouy, G. (2017). Sleeping on the motor engram: the multifaceted nature of sleep-related motor memory consolidation. Neurosci. Biobehav. Rev. 80, 1–22. doi: 10.1016/j.neubiorev.2017.04.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Klinzing, J. G., Niethard, N., and Born, J. (2019). Mechanisms of systems memory consolidation during sleep. Nat. Neurosci. 22, 1598–1610. doi: 10.1038/s41593-019-0467-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnan, G. P., Chauvette, S., Shamie, I., Soltani, S., Timofeev, I., Cash, S. S., et al. (2016). Cellular and neurochemical basis of sleep stages in the thalamocortical network. Elife 5:e18607. doi: 10.7554/eLife.18607

PubMed Abstract | CrossRef Full Text | Google Scholar

Kullback, S., and Leibler, R. A. (1951). On information and sufficiency. Ann. Math. Stat. 22, 79–86. doi: 10.1214/aoms/1177729694

CrossRef Full Text | Google Scholar

Lacourse, K., Delfrate, J., Beaudry, J., Peppard, P., and Warby, S. C. (2019). A sleep spindle detection algorithm that emulates human expert spindle scoring. J. Neurosci. Methods 316, 3–11. doi: 10.1016/j.jneumeth.2018.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Lacourse, K., Yetton, B., Mednick, S., and Warby, S. C. (2020). Massive online data annotation, crowdsourcing to generate high quality sleep spindle annotations from EEG data. Sci. Data 7, 1–14. doi: 10.1038/s41597-020-0533-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ladenbauer, J., Külzow, N., Passmann, S., Antonenko, D., Grittner, U., Tamm, S., et al. (2016). Brain stimulation during an afternoon nap boosts slow oscillatory activity and memory consolidation in older adults. Neuroimage 142, 311–323. doi: 10.1016/j.neuroimage.2016.06.057

PubMed Abstract | CrossRef Full Text | Google Scholar

Ladenbauer, J., Ladenbauer, J., Külzow, N., de Boor, R., Avramova, E., Grittner, U., et al. (2017). Promoting sleep oscillations and their functional coupling by transcranial stimulation enhances memory consolidation in mild cognitive impairment. J. Neurosci. 37, 7111–7124. doi: 10.1523/JNEUROSCI.0260-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Langdon, A. J., Breakspear, M., and Coombes, S. (2012). Phase-locked cluster oscillations in periodically forced integrate-and-fire-or-burst neuronal populations. Physi. Rev. E 86, 061903. doi: 10.1103/PhysRevE.86.061903

PubMed Abstract | CrossRef Full Text | Google Scholar

Latchoumane, C.-F. V., Ngo, H.-V. V., Born, J., and Shin, H.-S. (2017). Thalamic spindles promote memory formation during sleep through triple phase-locking of cortical, thalamic, and hippocampal rhythms. Neuron 95, 424–435. doi: 10.1016/j.neuron.2017.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Laureys, S., Perrin, F., and Brédart, S. (2007). Self-consciousness in non-communicative patients. Conscious Cogn. 16, 722–741. doi: 10.1016/j.concog.2007.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lüthi, A., and McCormick, D. A. (1998). Periodicity of thalamic synchronized oscillations: the role of Ca2+-mediated upregulation of Ih. Neuron 20, 553–563. doi: 10.1016/S0896-6273(00)80994-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayer, J., Schuster, H. G., and Claussen, J. C. (2006). Role of inhibitory feedback for information processing in thalamocortical circuits. Phys. Rev. E 73, 031908. doi: 10.1103/PhysRevE.73.031908

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A. (1992). Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. Progr. Neurobiol. 39, 337–388. doi: 10.1016/0301-0082(92)90012-4

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A., and Williamson, A. (1991). Modulation of neuronal firing mode in cat and guinea pig lgnd by histamine: possible cellular mechanisms of histaminergic control of arousal. J. Neurosci. 11, 3188–3199. doi: 10.1523/JNEUROSCI.11-10-03188.1991

PubMed Abstract | CrossRef Full Text | Google Scholar

Mölle, M., Bergmann, T. O., Marshall, L., and Born, J. (2011). Fast and slow spindles during the sleep slow oscillation: disparate coalescence and engagement in memory processing. Sleep 34, 1411–1421. doi: 10.5665/SLEEP.1290

PubMed Abstract | CrossRef Full Text | Google Scholar

Mölle, M., and Born, J. (2011). Slow oscillations orchestrating fast oscillations and memory consolidation. Prog. Brain Res. 193, 93–110. doi: 10.1016/B978-0-444-53839-0.00007-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Mölle, M., Marshall, L., Gais, S., and Born, J. (2002). Grouping of spindle activity during slow oscillations in human non-rapid eye movement sleep. J. Neurosci. 22, 10941–10947. doi: 10.1523/JNEUROSCI.22-24-10941.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Nir, Y., Staba, R. J., Andrillon, T., Vyazovskiy, V. V., Cirelli, C., Fried, I., et al. (2011). Regional slow waves and spindles in human sleep. Neuron 70, 153–169. doi: 10.1016/j.neuron.2011.02.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Oyanedel, C. N., Durán, E., Niethard, N., Inostroza, M., and Born, J. (2020). Temporal associations between sleep slow oscillations, spindles and ripples. Eur. J. Neurosci. 52, 4762–4778. doi: 10.1111/ejn.14906

PubMed Abstract | CrossRef Full Text | Google Scholar

Paluš, M. (1995). Testing for nonlinearity using redundancies: Quantitative and qualitative aspects. Physica D 80, 186–205. doi: 10.1016/0167-2789(95)90079-9

CrossRef Full Text | Google Scholar

Paluš, M. (1997). Detecting phase synchronization in noisy systems. Phys. Lett. A 235, 341–351. doi: 10.1016/S0375-9601(97)00635-X

CrossRef Full Text | Google Scholar

Peyrache, A., Dehghani, N., Eskandar, E. N., Madsen, J. R., Anderson, W. S., Donoghue, J. A., et al. (2012). Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep. Proc. Natl. Acad. Sci. U.S.A. 109, 1731–1736. doi: 10.1073/pnas.1109895109

PubMed Abstract | CrossRef Full Text | Google Scholar

Piantoni, G., Halgren, E., and Cash, S. S. (2016). The contribution of thalamocortical core and matrix pathways to sleep spindles. Neural Plast. 2016, 3024342. doi: 10.1155/2016/3024342

PubMed Abstract | CrossRef Full Text | Google Scholar

Popa, D., Duvarci, S., Popescu, A. T., Léna, C., and Paré, D. (2010). Coherent amygdalocortical theta promotes fear memory consolidation during paradoxical sleep. Proc. Natl. Acad. Sci. U.S.A. 107, 6516–6519. doi: 10.1073/pnas.0913016107

PubMed Abstract | CrossRef Full Text | Google Scholar

Potjans, T. C., and Diesmann, M. (2014). The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cereb. Cortex 24, 785–806. doi: 10.1093/cercor/bhs358

PubMed Abstract | CrossRef Full Text | Google Scholar

Purcell, S., Manoach, D., Demanuele, C., Cade, B., Mariani, S., Cox, R., et al. (2017). Characterizing sleep spindles in 11,630 individuals from the national sleep research resource. Nat. Commun. 8, 1–16. doi: 10.1038/ncomms15930

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasch, B., and Born, J. (2013). About sleep's role in memory. Physiol. Rev. 93, 681–766. doi: 10.1152/physrev.00032.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Richardson, M. J. (2007). Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive. Phys. Rev. E 76, 021919. doi: 10.1103/PhysRevE.76.021919

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, P., Rennie, C., and Rowe, D. (2002). Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E 65, 041924. doi: 10.1103/PhysRevE.65.041924

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosanova, M., and Ulrich, D. (2005). Pattern-specific associative long-term potentiation induced by a sleep spindle-related spike train. J. Neurosci. 25, 9398–9405. doi: 10.1523/JNEUROSCI.2149-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, F., Wibral, M., Singer, W., Aru, J., and Uhlhaas, P. J. (2013). The phase of thalamic alpha activity modulates cortical gamma-band activity: evidence from resting-state MEG recordings. J. Neurosci. 33, 17827–17835. doi: 10.1523/JNEUROSCI.5778-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubio-Garrido, P., Pérez-de Manzo, F., Porrero, C., Galazo, M. J., and Clascá, F. (2009). Thalamic input to distal apical dendrites in neocortical layer 1 is massive and highly convergent. Cereb. Cortex 19, 2380–2395. doi: 10.1093/cercor/bhn259

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez-Vives, M. V., and McCormick, D. A. (2000). Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat. Neurosci. 3, 1027–1034. doi: 10.1038/79848

PubMed Abstract | CrossRef Full Text | Google Scholar

Sawangjit, A., Oyanedel, C. N., Niethard, N., Salazar, C., Born, J., and Inostroza, M. (2018). The hippocampus is crucial for forming non-hippocampal long-term memory during sleep. Nature 564, 109–113. doi: 10.1038/s41586-018-0716-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Schellenberger Costa, M., Weigenand, A., Ngo, H.-V. V., Marshall, L., Born, J., Martinetz, T., et al. (2016). A thalamocortical neural mass model of the EEG during NREM sleep and its response to auditory stimulation. PLoS Comput. Biol. 12, e1005022. doi: 10.1371/journal.pcbi.1005022

PubMed Abstract | CrossRef Full Text | Google Scholar

Schomer, D. L., and Da Silva, F. L. (2012). Niedermeyer's Electroencephalography: Basic Principles, Clinical Applications, and Related Fields. Oxford: Oxford University Press. doi: 10.1093/med/9780190228484.001.0001

CrossRef Full Text | Google Scholar

Schreiber, T., and Schmitz, A. (2000). Surrogate time series. Physica D 142, 346–382. doi: 10.1016/S0167-2789(00)00043-9

CrossRef Full Text | Google Scholar

Silber, M. H., Ancoli-Israel, S., Bonnet, M. H., Chokroverty, S., Grigg-Damberger, M. M., Hirshkowitz, M., et al. (2007). The visual scoring of sleep in adults. J. Clin. Sleep Med. 3, 121–131. doi: 10.5664/jcsm.26814

CrossRef Full Text | Google Scholar

Sirota, A., Csicsvari, J., Buhl, D., and Buzsáki, G. (2003). Communication between neocortex and hippocampus during sleep in rodents. Proc. Natl. Acad. Sci. U.S.A. 100, 2065–2069. doi: 10.1073/pnas.0437938100

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M. (2003). The corticothalamic system in sleep. Front. Biosci. 8, 1043. doi: 10.2741/1043

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., McCormick, D. A., and Sejnowski, T. J. (1993). Thalamocortical oscillations in the sleeping and aroused brain. Science 262, 679–685. doi: 10.1126/science.8235588

PubMed Abstract | CrossRef Full Text | Google Scholar

Suffczynski, P., Kalitzin, S., and Da Silva, F. L. (2004). Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience 126, 467–484. doi: 10.1016/j.neuroscience.2004.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Timofeev, I., and Bazhenov, M. (2005). “Mechanisms and biological role of thalamocortical oscillations,” in Trends in Chronobiology Research (Hauppauge, NY: Nova Science Publishers, Inc.), 1–47.

Google Scholar

Timofeev, I., and Steriade, M. (1996). Low-frequency rhythms in the thalamus of intact-cortex and decorticated cats. J. Neurophysiol. 76, 4152–4168. doi: 10.1152/jn.1996.76.6.4152

PubMed Abstract | CrossRef Full Text | Google Scholar

Tort, A. B., Komorowski, R., Eichenbaum, H., and Kopell, N. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J. Neurophysiol. 104, 1195–1210. doi: 10.1152/jn.00106.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Touboul, J., Wendling, F., Chauvel, P., and Faugeras, O. (2011). Neural mass activity, bifurcations, and epilepsy. Neural Comput. 23, 3232–3286. doi: 10.1162/NECO_a_00206

PubMed Abstract | CrossRef Full Text | Google Scholar

Ujma, P. P., Hajnal, B., Bódizs, R., Gombos, F., Erőss, L., Wittner, L., et al. (2021). The laminar profile of sleep spindles in humans. Neuroimage 226, 117587. doi: 10.1016/j.neuroimage.2020.117587

PubMed Abstract | CrossRef Full Text | Google Scholar

Vallat, R., and Jajcay, N. (2020). YASA. Available online at: https://github.com/raphaelvallat/yasa.

Vanini, G., Wathen, B. L., Lydic, R., and Baghdoyan, H. A. (2011). Endogenous GABA levels in the pontine reticular formation are greater during wakefulness than during rapid eye movement sleep. J. Neurosci. 31, 2649–2656. doi: 10.1523/JNEUROSCI.5674-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, M. P., and Stickgold, R. (2004). Sleep-dependent learning and memory consolidation. Neuron 44, 121–133. doi: 10.1016/j.neuron.2004.08.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Warby, S. C., Wendt, S. L., Welinder, P., Munk, E. G., Carrillo, O., Sorensen, H. B., et al. (2014). Sleep-spindle detection: crowdsourcing and evaluating performance of experts, non-experts and automated methods. Nat. Methods 11, 385. doi: 10.1038/nmeth.2855

PubMed Abstract | CrossRef Full Text | Google Scholar

Werth, E., Achermann, P., Dijk, D.-J., and Borbély, A. A. (1997). Spindle frequency activity in the sleep EEG: individual differences and topographical distribution. Electroencephalogr. Clin. Neurophysiol. 103, 535–542. doi: 10.1016/S0013-4694(97)00070-9

PubMed Abstract | CrossRef Full Text | Google Scholar