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

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.

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

with membrane time constant τ, synaptic input rate w_e (w_i) that scales synaptic inputs s_e (s_i) for excitatory (e) and inhibitory (i) synapses, the corresponding Nernst reversal potential E_e (E_i), and the membrane capacitance C_m. The mean membrane potential V_α is then converted to a firing rate r(V) by a sigmoidal transfer function,

with maximum firing rate r_max, firing threshold θ, and gain coefficient σ.

Both thalamic populations contain additional intrinsic currents [I^intrinsic in Equation (1)] because spindle oscillations require rebound burst activity. Rebound bursting is impossible with a monotonic firing rate function and demands additional mechanisms. Following Schellenberger Costa et al. (2016), we employ the Hodgkin-Huxley-type extension, which has been derived from integrate-and-fire-or-burst neurons (Langdon et al., 2012).

The intrinsic currents in both thalamic populations include a potassium leak current,

and a T-type calcium current,

which de-inactivates upon depolarization. In both definitions, g_LK and g_T denote the conductance of the respective intrinsic current, E_K and E_Ca their respective Nernst reversal potential, and $m_{\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:

with g_h being its conductance, E_h its Nernst reversal potential, m_h1 and m_h2 the gating functions, and g_inc the conductivity scaling. This hyperpolarization-activated current is responsible for the waxing and waning of spindle oscillations in the isolated thalamus (Destexhe et al., 1996a).

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

$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

representing the synaptic response to a single spike. γ_l is the decay constant of the synaptic response, and l ∈ {e, i} denotes the type of synapse, i.e., excitatory AMPA or inhibitory GABA_A. The convolution in Equation (6) is replaced in the numerical simulations by the second-order ordinary differential equation (ODE)

The background noise input ϕ′(t) is modeled as an Ornstein-Uhlenback process with zero drift but finite variance and represents unresolved processes in our model, e.g., afferents from other brain regions that are not explicitly modeled here.

In summary, the thalamic node is described by the set of equations

where subscripts t (r) represent the TCR (TRN) population. Note that the background noise is included on the excitatory synaptic input at the TCR population, s_et. The full set of equations is given in Supplementary Material, with all parameters summarized in Supplementary Table S1.

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

The first equation describes the temporal evolution of neuron's j membrane voltage V_j as a function of its internal current dynamics conveyed by I_ion(V_j), its synaptic current I_syn(t), and a background external current I_j,ext(t) received from neural populations not specified by the computational model. C_m,c denotes the cortical neuron's membrane capacitance. Note that inhibitory population α = I does not include the adaptation current. The first term in Equation (12) expresses the voltage-dependent leak current with leak conductance g_L and leak reversal potential E_L; the second term describes the exponential spike initiation mechanism with slope factor Δ_T and exponential threshold V_T. Finally, the last term describes the somatic adaptation current, I_A,j(t) (Equation 13), with subthreshold adaptation a, adaptation reversal potential E_A, and adaptation time scale τ_A. When the membrane voltage crosses the spiking threshold, V_j ≥ V_s, the voltage is reset with V_j ← V_r, clamped for the refractory time T_ref, and the spike-triggered adaptation increment, b, is added to the adaptation current, I_A,j ← I_A,j + b.

The synaptic currents are described by a sum of excitatory and inhibitory contributions as

with the coupling strength J_αβ from population β to α, and the fraction of active synapses s_i,αβ ∈ [0, 1]. The synaptic dynamics is then given by

where G_ij is a random binary connectivity matrix. The first sum is over all afferent neurons j, while the second sum is over all incoming spikes k from neuron j emitted at time t_k after a delay d_αβ. The (1 − s_i,αβ) term in Equation (15) acts as a saturation term and integrates all incoming spikes only if s_i,αβ < 1, i.e., only if there is a synaptic capacity available.

In the mean-field approximation, the distribution p(V) of the membrane potentials and the mean population firing rate r can be calculated using the Fokker-Planck equation in the thermodynamic limit, with the number of neurons N → ∞ (Brunel, 2000; Cakan and Obermayer, 2020). However, using the model reduction scheme described in Augustin et al. (2017) and Cakan and Obermayer (2020), we exploit the low-dimensional linear-nonlinear cascade model, where for a given mean membrane current μ_α with standard deviation σ_α, the mean ${\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.,

for α ∈ {E, I}. In the case of the excitatory population (α = E), the mean adaptation current Ī_A is subtracted from the mean membrane current μ_E in the computation of population firing rate, mean membrane voltage, and adaptive timescale via transfer functions, i.e., μ_E → μ_E − Ī_A/C_m,c. These transfer functions can be precomputed for a specific set of single AdEx neuron parameters (Augustin et al., 2017) and are shown in Supplementary Figure S1.

The average population currents in the mean-field approximation are given by

where the dynamics of mean membrane current μ_α depends on the synaptic current ${\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

where the mean r_αβ and the variance ρ_αβ of the effective input rate from population β to population α for a spike transmission delay d_αβ are given by

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

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

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

The KL-MI estimates the phase-amplitude coupling by computing the Kullback-Leibler divergence (Kullback and Leibler, 1951) between the distribution of spindle amplitudes A_spindle(t) over the slow oscillation phase bins ϕ_SO(t) and a uniform distribution (the null hypothesis of no phase-amplitude coupling). The KL-MI is defined by

where N is the number of phase bins and H(P) is the Shannon entropy of the amplitude distribution.

The MVL is computed by averaging the complex time series constructed by multiplying the spindle amplitude, A_spindle(t), with the term containing the phase, ϕ_SO, of the slow oscillations, thus

The length (real part) of the MVL vector quantifies the strength of phase-amplitude coupling. Its phase denotes the mean phase of the slow oscillation at which the spindle amplitude is the strongest.

The PLV, as a measure of phase locking between two oscillations, is computed by temporally averaging the phase differences in unit circle between the phase of the slow oscillation and the phase of the spindle oscillation:

As with the MVL, the length of the PLV vector indicates the strength of phase-locking, while the phase represents the phase shift.

Finally, the MI of two discrete variables (in our case, the time series of estimated phases) can be computed using

where p(x, y) signifies the joint probability mass function of {X, Y}, and p(x) and p(y) represent the marginal probability mass functions of X and Y, respectively. All three probability mass functions are estimated using an equiquantal binning algorithm (Paluš, 1995) with 16 bins.

To test for statistical significance, we computed 1,000 Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogates from a randomization procedure that preserves the power spectrum and the amplitude distribution of the original time series (Schreiber and Schmitz, 2000). The surrogates were constructed by iterative replacements of Fourier amplitudes with the values from the original time series and by rescaling the distribution match between the distribution and the power spectrum of the original data. We opted to use IAAFT surrogates since our modeled data are firing rates with non-Gaussian distribution and do not meet the criteria for using the basic Fourier Transform surrogates. After constructing the surrogates, we computed a given CFC measure for each surrogate time series, yielding an empirical null distribution. Finally, we compared a result on the modeled output to this null distribution to obtain the empirical p-value.

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

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

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.

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.

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.

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

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

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.

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.

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