http://orcid.org/0000-0002-9870-4924
Figures
Abstract
Few models exist that accurately reproduce the complex rhythms of the thalamocortical system that are apparent in measured scalp EEG and at the same time, are suitable for large-scale simulations of brain activity. Here, we present a neural mass model of the thalamocortical system during natural non-REM sleep, which is able to generate fast sleep spindles (12–15 Hz), slow oscillations (<1 Hz) and K-complexes, as well as their distinct temporal relations, and response to auditory stimuli. We show that with the inclusion of detailed calcium currents, the thalamic neural mass model is able to generate different firing modes, and validate the model with EEG-data from a recent sleep study in humans, where closed-loop auditory stimulation was applied. The model output relates directly to the EEG, which makes it a useful basis to develop new stimulation protocols.
Author Summary
Sleep plays a pivotal role for the consolidation of memory. While REM sleep had originally been the focus of research due to its similarity with wakefulness, more recent studies suggest that different sleep stages are responsible for the consolidation of different types of memory. To better understand the changes in neuronal dynamics between the different sleep stages, neural mass models are a valuable tool, as they relate directly to the large-scale dynamics measured by an EEG. Here, we present a model of the sleeping thalamocortical system. We first show that the isolated thalamic submodule is able to generate different oscillatory behavior found in vivo. Furthermore, the full thalamocortical model reproduces the EEG of sleep stages N2 and N3 and preserves the temporal relationship between cortical K-complexes/slow oscillations and thalamic sleep spindles. A comparison with event related potential data from a recent sleep study in humans demonstrates its possible application in predicting the outcome of external stimulation on EEG rhythms. Our study shows, that a neural mass model incorporating few key mechanisms is sufficient to reproduce the complex brain dynamics observed during sleep.
Citation: Schellenberger Costa M, Weigenand A, Ngo H-VV, 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(9): e1005022. https://doi.org/10.1371/journal.pcbi.1005022
Editor: Jean Daunizeau, Brain and Spine Institute (ICM), FRANCE
Received: May 18, 2015; Accepted: June 17, 2016; Published: September 1, 2016
Copyright: © 2016 Schellenberger Costa et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This work was supported by the grants "Plasticity and Sleep" (SFB 654) and the Graduate School "Computing in Medicine and Life Sciences" (GS-635) from the German Research Foundation and US-German Collaboration in Neuroscience (grant 01GQ1008) from the German Federal Ministry of Education and Research and the European Union Human Brain Project SP3 - Cognitive Architectures. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
In the human electroencephalogram (EEG) the most prominent features of non-rapid eye movement (NREM) sleep are neocortical slow oscillations (SO) at ∼ 0.8 Hz [1] and fast thalamocortical spindles characterized by a waxing and waning waveform, with a frequency of 12 to 15 Hz [2]. The cellular basis of the slow oscillation is a widespread alternation of activity in neocortical networks between an active (up) state and a hyperpolarized silent (down) state [2–4]. Fast sleep spindles are generated by the interaction of inhibitory reticular thalamic (RE) and excitatory thalamocortical (TC) neurons [5] and preferentially occur during the up state of SOs [6].
Many studies indicate a functional role of slow wave sleep (SWS) in the formation of memory [7, 8]. The synchronization of fast spindle activity to the depolarized up state is mediated by the thalamocortical loop. It appears to be critical for consolidation as it provides a window of opportunity and favorable conditions for plastic changes [9–11]. It has been shown that memory consolidation can be improved by transcranial electric and auditory stimulation [12, 13]. Auditory stimulation in synchrony with the brain’s own rhythm turned out to be particularly effective.
Detailed knowledge of how different stimulation modalities effect critical brain rhythms as well as the ramifications of potential stimulation parameters on their efficacy would enable an optimization of stimulation protocols and consequently an advantage for experiments in basic research and clinical applications. Mathematical models and computational approaches can yield meaningful insights into the underlying dynamics as well as provide predictions for further experiments. Multiple models based on networks of single cells have addressed brain activity during sleep in the cortex [14], thalamus [15–17], as well as the coupled thalamocortical system [18–21]. However, for investigations of macroscopic phenomena these models have disadvantages due to their complexity and computational load.
Neural mass models have shown great success in elucidating the generation of brain rhythms and evoked responses of the awake brain [22–25]. They describe the dynamics of a large number of cells by the evolution of a single population average and provide an output which directly relates to EEG signals [26, 27]. By still allowing an integration of physiologically motivated cell dynamics via intrinsic currents, neural mass models represent a compromise between a very detailed and an abstract modeling approach and provide insights into the dynamic repertoire of the neural populations via bifurcation analysis [28–30].
In previous work, we have shown that a cortical neural mass model equipped with an additive activity-dependent feedback current can generate a time series that closely resembles the EEG signal of sleep stages N2 and N3, without spindles [30]. Here, we extend the cortical model by adding a thalamic module to incorporate spindle activity and investigate the underlying dynamics of the coupled system. We test the evocability of SOs and spindles by auditory stimulation during NREM sleep and validate our results with scalp EEG data from a recent sleep study in humans [13]. Our findings add further support to the dynamic mechanisms proposed in [30].
Methods
In the following sections we detail the thalamocortical model. First we give a brief overview of the concept of neural mass models. We then introduce the cortical and thalamic modules and motivate the extensions with respect to sleep. The full equations of the thalamocortical model are given in S1 Text.
Neural mass framework
The conductance based neural mass model employed here is derived from [31–34]. Instead of considering the evolution of high-dimensional states in a large ensemble of single neurons, the population activity can be approximated by the evolution of the mean membrane voltage of that population. The complex spiking dynamics is replaced by an empirical firing rate function
(1)
with maximal firing rate
, firing threshold θ and neural gain σk. It converts the average membrane voltage Vk of the population k to an output spike rate. Here, k ∈ {p, i, t, r} stands for cortical pyramidal, cortical interneuron, thalamic relay and thalamic reticular populations, respectively. The firing rate function often has a sigmoidal shape and can be interpreted as stemming from the fluctuations of neuronal states or a distribution of thresholds in the population [35].
The spike rate Qk′ of the presynaptic population k′ then elicits a postsynaptic response smk within the receiving population k. The strength of this response can be calculated by the convolution
(2)
of the incoming spike rate Qk′, scaled with the averaged connectivity constant Nkk′ between the presynaptic population k′ and the postsynaptic population k and an alpha function
(3)
representing the average synaptic response to a single spike. Here, γm depicts the rate constant of the synaptic response, whereas m ∈ {e, g, r} denotes the type of synapse with e standing for excitatory AMPA and g, r for inhibitory GABA synapses in the cortex and thalamus, respectively.
The evolution of the population membrane voltage Vk obeys
(4)
where w denotes a synaptic input rate that scales the input smk and E the corresponding Nernst potential. While we lean on the naming convention of Hodgkin-Huxley models to highlight the structural similarity, please note that the above quantities J and w have different units. For the sake of simplicity we have normalized the synaptic rates w to 1, absorbing their numerical values into the connectivities Nkk′.
Intrinsic currents may additionally be included in the equation of the mean membrane voltage, given their time constant is large compared to the time constant of neuronal spiking [35, 36]. To emphasize the connection to physiology and to better differentiate between the core neural mass model and the additional mechanisms we will denote them with I and introduce a capacitive coupling via to the neural mass via the membrane capacity Cm.
The signal measured in the EEG stems mainly from the activity of pyramidal neurons [37]. We use the pyramidal membrane voltage as our model output, which is similar to other studies [38, 39]. Populations comprising multiple clusters have been considered in [40] and lead to interesting effects. In order to keep the complexity of the model low we consider a single point source. Therefore, filtering effects by the skull/scalp can be approximated by a linear scaling and do not affect the interaction between thalamus and cortex.
Cortex
We use the neural mass model described in [30], which captures essential features of the EEG during NREM sleep. In brief, it consists of an excitatory and an inhibitory neural mass, representing populations of pyramidal neurons (p) and interneurons (i), respectively. The pyramidal population contains a slow, additive and activity-dependent firing rate adaptation,
(5)
which is believed to be the main driver for the transition to the silent (down) state, while the active (up) state is maintained by synaptic transmission [3, 14, 41].
The membrane potentials evolve according to
(6)
and are linked by the AMPA and GABAergic synaptic inputs.
The transition from N2 to N3 is achieved by changing the neural gain, σp, and the strength of the adaptation, given by . Both parameters are known to be influenced by acetylcholine, serotonin, norepinephrine and dopamine [42–52], whose levels change throughout the sleep-wake cycle [53]. For consistency, we maintained the original model, adjusting only the excitatory connection strengths to compensate for additional input from the thalamic module.
Thalamus
The thalamic module comprises similarly an excitatory and an inhibitory neural mass, representing a thalamocortical (t) and the reticular (r) nucleus. As in the cortical module module they are coupled via AMPA and GABA synapses but have different synaptic time constants and only the RE population possesses a self-connection.
Both populations are equipped with additional currents. The inclusion of those currents within the thalamic submodule is necessary because spindle oscillations require rebound bursts. In classical neural mass models, this kind of bursting is not possible due to the monotonic firing rate function and demands the inclusion of additional mechanisms. The same argument was used in a previous neural mass model of spindle activity [54] and a thalamocortical neural mass model of epileptic activity [55]. Finally, in [40] the authors arrive at a HH-type extension of their population model of thalamic burst activity, which has been derived from integrate-and-fire-or-burst neurons. The potassium leak current is given by
(7)
as well as T-type calcium currents
(8)
which deinactivate upon hyperpolarization. They are essential for the generation of low-threshold spikes (LTSs) and rebound bursts. We use the description of IT given in [56] for the RE and the one in [16] for the TC population.
The TC population further includes the anomalous rectifier current
(9)
responsible for the waxing and waning structure of spindle oscillations in the isolated thalamus [15]. Other currents, such as the calcium-dependent potassium currents IKCa and ICAN, are also known to play a role in spindle oscillations, but are omitted for simplicity. The thalamic module is summarized by
(10)
Parameter settings for the currents are identical to [57], with the exception of the deactivation function
of the thalamic relay population, that is shifted towards more depolarized membrane voltages.
Full model
The model consists of one thalamic and one cortical module. We assume the long range afferents from the cortical pyramidal population project to both populations of the thalamic nuclei, and the long range afferents of the thalamic relay population project to both populations of the cortex, as depicted in Fig 1. The delays introduced by these long range afferents might play a crucial role in cortical dynamics [58, 59]. However, as the axonal conduction delay between thalamus and cortex is rather small [60–63], we approximate it by a convolution with an alpha function [64], which can be written as
(11)
where ϕk is the resulting delayed firing rate and ν depicts the axonal rate constant of that connection. We have provided a justification for that approximation in the supporting information S2 Text. In the case of short range connections ϕk can be replaced with Qk. The parameters of the full model are given in S1 Table.
Excitatory synapses are depicted by filled circles, inhibitory synapses by bars. Independent background noise entering the different populations is denoted by ϕn, and
, respectively. Stimulation is applied as an elevation in the mean of the background noise
of the thalamic relay population.
Auditory stimulation
We model an auditory stimulus as an elevation in background noise (square pulse) being gated through the thalamus. For all stimuli, we use a duration of 80 ms and 70 spikes per second.
Event triggered averages
KCs and SOs were detected similar to [10]. The model output was bandpass filtered between 0.25 and 4 Hz. Zero-crossings were detected to extract the negative half-waves. Negative half-waves with peaks below -69 mV were considered to be KCs/SOs.
Experimental data
Experimental data has been described in [65]. 11 Subjects were measured during two experimental nights in balanced order in a stimulation and control condition. For averages of the endogenous activity, data was taken from the control condition. See S1 Dataset.
Results
In the following section we perform a bifurcation analysis and demonstrate the ability of the thalamic module to generate spindle oscillations and reproduce different experimental observations. Afterwards we investigate the interplay between thalamus and cortex to reproduce the characteristics of different sleep stages. Finally, we examine the effect of auditory stimulation in the model and compare different stimulation protocols with experimental findings.
Thalamic spindle oscillations
In the isolated thalamic module, incorporation of the intrinsic currents may lead to oscillations in the spindle band. We follow closely the mechanisms established in the models by [15, 19, 20]. Physiologically, these oscillations are generated, through reciprocal interaction of the RE and TC populations. A LTS in the RE population causes hyperpolarization in the TC population, that deinactivates its T-type calcium current. Upon release from inhibition a rebound of activity occurs, that in turn drives the RE module to produce another LTS. Additionally, the deinactivation of the T-type calcium currents requires a strong tonic hyperpolarization by a potassium leak current [15, 19].
As previously shown in [5, 15, 20, 69], the rhythmicity of spindle occurrence and the waxing and waning of the spindle amplitude is caused by an anomalous rectifier channel Ih. A sequence of LTS leads to the build-up of calcium, which increases the effective conductivity of Ih. The ensuing depolarization of the TC population increasingly counteracts its ability to produce a LTS and terminates the spindle oscillation [70, 71].
Therefore, we chose and
as bifurcation parameters. A two-dimensional bifurcation analysis of the thalamic module reveals the existence of a Hopf bifurcation, as depicted in Fig 2, which generates continuous oscillations in the spindle band due to hyperpolarization induced rebound bursts, see Fig 3-CI and 3-CII for representative time series.
Here, we illustrate the bifurcation diagram of the isolated thalamus with respect to the two key parameters and
. The interaction between the currents incorporated into the thalamic module results in the emergence of two torus bifurcations via a blue sky catastrophe. They lead to spindle oscillations in the orange shaded regions. The left spindle regime (SI) is encased by a Hopf and a torus bifurcation, whereas the right spindle regime (SII) is constrained by two global bifurcations that are indicated by the dashed gray lines. The vertical line marks the emergence of the torus bifurcation, whereas the horizontal gray line marks the cusp bifurcation where the two saddle-nodes that accompany the left torus bifurcation vanish. The torus bifurcation on the right marks the transition from spindle oscillations to delta oscillations. The labeled points mark the parameter settings utilized in Fig 3, which are given in Table 1.
Here, we illustrate the different dynamic modes the isolated thalamic module exhibits. The left panels depict the thalamic relay membrane voltage, whereas the right panels illustrate that of the thalamic reticular population. The parameter values are depicted in Fig 2 and given in Table 1. SI and SII: The isolated thalamus generates rhythmic spindle oscillations via a balanced interplay between IT and Ih. The length and the average time between spindles is governed by . CI and CII: Outside of the spindle regime fast oscillations generated by the T-type calcium currents dominate and Ih is unable to sufficiently depolarize the thalamic relay population to cease them. DI and DII: For strong hyperpolarization through ILK the thalamic module switches into low frequency delta oscillations.
The torus bifurcations emerge from a blue sky catastrophe that is generated by the slow-fast interaction between the fast T-type channels and their slow modulation via Ih, which is similar to other models that exhibit switching between tonic spiking and structured bursting activity [72, 73].
As depicted in Fig 3-SI and 3-SII this leads to spindle like oscillations in the orange shaded regions in Fig 2. The spindles exhibit an oscillation frequency of around 13 Hz. The spindle frequency depends on the strength of the T-type calcium current . Importantly, spindle oscillations are initiated intrinsically. The thalamic module does not require modulatory input from external sources to initiate/terminate them.
Activation of Ih is responsible for a refractory period that follows a spindle. As long as Ih activation persists, LTS generation is impeded and stronger perturbations are necessary to trigger spindle oscillations. Consequently, an increase in results in a larger inter-spindle interval. The left spindle regime (SI) is encased by the Hopf and the torus bifurcation, whereas the right spindle regime (SII) is constrained by two global bifurcations that are indicated by the dashed gray lines. The vertical line marks the emergence of the torus bifurcation, whereas the horizontal gray line marks the cusp bifurcation where the two saddle-nodes that accompany the left torus bifurcation vanish.
Furthermore, for larger values of the model transitions from high frequency spindle oscillations to low frequency delta oscillations, e.g. Fig 3-DI and 3-DII.
K-complexes and spindles during sleep stage N2
In the coupled system, the cortex provides excitatory drive to the thalamic module, since it is predominantly in the active state. In order for the thalamic module to exhibit rhythmically occurring spindle oscillations we had to adjust and
slightly to account for that additional depolarization (see Table 2). We choose the right spindle regime, as it was suitable for reproduction of both sleep stage N2 and N3.
As can be seen in Fig 4 spindles may be triggered by KCs in the full model, but may also occur independent of KCs. During a KC the sudden drop of excitatory drive hyperpolarizes the RE and TC population, leading to deinactivation of IT. The ensuing depolarization upon the transition back to the active state triggers a LTS and a spindle sequence in turn. The spindle then projects back into the depolarizing phase of the KC. This is in good agreement with the grouping of spindles and KCs observed experimentally [6, 74].
Shown are membrane voltages of the cortical pyramidal (top) and the thalamic relay population (bottom). The spindle oscillations induced within the thalamic module project into the cortical module. While the spindle oscillations are generally induced by fluctuations in background noise, there is also a grouping between cortical KCs and thalamic spindles (see 7s-9s and 19s-21s). The grouping stems from the lack of depolarizing input during a cortical KC. Parameters are as in Table 2.
Although less likely the model can also give rise to KCs triggered by a spindle. This can be achieved by increasing the connection strength from the thalamic to the cortical module (model output not shown). During N2, KCs occur at a low rate. Hence, spindle initiation and termination are closely linked to the time course of Ih (Fig 3A), similar to the isolated thalamic module. The parameters for the output in Fig 4 are given in Table 2.
Given the parameter setting in Table 2, the cortical module is within a stable focus, close to a Hopf bifurcation accompanied by a canard explosion. This leads to noise driven medium amplitude background oscillations around the stable focus, that are interrupted by large amplitude deflections (KCs). In good agreement with experimental findings, KCs also appear within the isolated cortex, although they may be initiated through thalamic input.
SOs and spindles during sleep stage N3
On the transition to sleep stage N3 the canard phenomenon vanishes in a cusp bifurcation and only a high amplitude limit cycle remains. SOs are noise driven oscillations around a stable focus, close to a Hopf bifurcation [30].
In contrast to sleep stage N2 spindle initiation and termination are now dominated by the modulatory input from the cortical module, that overrules the Ih rhythm. Rather than occurring rhythmically spindles are time-locked to the depolarized phase of a SO. In Fig 5 an example time series is shown. Importantly, not every SO is able to trigger a spindle, as can be seen in Fig 5 (9–12 s, 13–15 s). We observed that in a sequence of SOs the first triggers a spindle, which leads to an activation of Ih. This reduces spindle amplitude or even inhibits spindle initiation by the following SO.
Shown are membrane voltages of the cortical pyramidal (top) and the thalamic relay population (bottom). During N3 the model shows ongoing slow oscillatory activity. In contrast to sleep stage N2, SOs cannot be identified as isolated events. Furthermore, there are no isolated spindle oscillations and spindle activity is time-locked to SOs. Parameters are given in Table 2.
Endogenous event triggered averages
To further validate the model, we determined averages of the generated EEG signal and fast spindle power time-locked to the negative peak of the endogenous KCs/SOs during N2 and N3. This method is often used to illustrate the grouping of spindles by SOs and morphological features of SOs, e.g. in [6, 13, 65]. Model output and data for N2 and N3 is depicted in Fig 6.
Averaged EEG signal (top) and fast spindle band power (bottom) time-locked to the negative peaks (t = 0 s) of all detected events from electrode Cz (black, left axis) and model output (red, right axis). (A) Detected KCs from data scored as sleep stage N2 (Experiment: 227,45 ± 19,22, Model: 238 events). (B) SO average from data scored as sleep stage N3 (Experiment: 983,64 ± 106,1, Model: 654 events). Each simulation was run for 3600 s with parameters set according to Table 2.
As can be seen in Fig 6, the grouping of spindles by SOs is present in the model. Spindle power is highest during the positive half-wave following the negative peak. However, there are some notable differences. Compared to the experimental data the initial depolarization preceding the transition to the down state is less prominent, leading to a shallower slope of the transition to the down state. In the thalamocortical model the transition to the depolarized up state occurs considerably earlier with a time to peak of 300 ms, compared to 440 ms in the data. This stems from strong depolarizing input by thalamic spindle bursts, which start directly after the negative peak of a KC/SO and push the cortex further into the depolarized state. However, this is still in line with other experimental studies, that find different timings of spindles for the supplementary motor area of the cortex [75].
Closed-loop auditory stimulation
In the following we show the ability of the model to reproduce data from a recent experiment in humans performing auditory closed-loop stimulation during NREM sleep [13]. The stimulation protocol is as follows: After the negative peak of a SO was detected, two auditory stimuli were applied phase-locked to the following positive peak of the depolarized up phase of the detected and the subsequent SO.
In the experimental study the delay time between the negative peak and the ensuing positive half-wave peak was determined for every subject independently. The second stimulus followed after a fixed interval of 1075 ms. Detection was then paused for 2.5 s. We accordingly determined the delay time from the model output, resulting in a delay of 450 ms for the N3 parameter setting. The second stimulus was chosen to occur 1075 ms after the first one and we also paused detection for 2.5 s. Stimuli are given as elevations in mean background noise of the thalamic relay population for a duration of 80 ms.
Fig 7 shows the averaged EEG signal and model output time-locked to the first stimulus (t = 0). There is a good agreement between model output and the experimental data. Especially the large amplitude, late components of the ERP are very close to the original waveform. The early component of the evoked potential, the P200, can be seen in the experimental data after each stimulus, but it is more pronounced in the model output.
The upper panel depicts in black the mean (± SEM) evoked potentials of human EEG data from electrode Cz during closed loop stimulation, time locked to the first stimulus (11 subjects, 245.6 ± 38.1 stimuli). In red the reproduction of the stimulation protocol with the model is shown (mean ± SD, 88 stimuli). The dashed line marks the stimulus onset. The lower panel shows the corresponding fast spindle power. Parameters used for model simulation are given in Table 2.
In addition, the evoked spindle responses of model and data also have similar time courses. In both cases spindle power is systematically increased during the depolarized up phases induced by the stimuli. However, the strong increase in spindle power seen in the data after the first stimulus is not visible in the model. We hypothesize this to stem from a recruitment effect, where the stimulus activates a larger fraction of the thalamus than the endogenous slow oscillation would. As our thalamic module is a point model without any spatial extent, these effects are excluded by construction.
Interestingly, in the experimental data there is a drop in spindle power after the second stimulus is applied. This seems to be a refractoriness of the thalamus after the second slow oscillation, which has also been observed in [76]. Despite the model showing such a refractory period in the isolated thalamus (Fig 3A), as well as during trains of endogenous SOs in the full model (Fig 8A), it lacks it upon stimulation (Fig 8B).
The upper two panels depict the membrane voltages of the pyramidal and thalamic relay populations, respectively. In the third panel the conductivity of the Ih current is shown. (A) Example time series of an unperturbed train of SOs during sleep stage N3. The first two SOs lead to an activation of Ih, that slowly declines back to baseline levels. As Ih activation is still well above baseline, the third SO is unable to trigger a spindle response. During the fourth SO Ih activation is sufficiently low so that a spindle occurs. (B) Shown is an example of closed loop stimulation during sleep stage N3, with the dashed lines indicating stimulus onset. In contrast to the endogenous case, the depolarization of the thalamic relay population induced by the stimulation leads to a rapid decrease in Ih activation, so that the following SO triggers a spindle. Parameters as in Table 2.
This happens because stimulation disturbs the Ih mediated spindle termination mechanism. As the stimulation depolarizes the TC population, the calcium concentration drops, because calcium influx through the IT current stops and calcium leaks out with a time constant of 10ms. Without the elevated calcium concentration, Ih deactivates back to baseline levels and immediately allows for a new full fledged spindle.
Discussion
We developed a neural mass model of the thalamocortical system that produces realistic time courses of EEG signals during sleep stages N2 and N3 and correctly replicates the timing of KCs and spindles. We validated our model with SO triggered averages of the EEG signal and spindle power. Finally, we used the model to reproduce evoked responses from closed-loop auditory stimulation during human NREM sleep.
Mechanisms of spindle generation
The model emphasizes the role of IT and Ih currents in the generation of thalamocortical rhythms as they were sufficient to reproduce the investigated EEG phenomena. It reproduces the grouping of spindles and KCs/SOs, observed in human EEG [6], that is thought to play a crucial role in the consolidation of memory [77, 78]. Additionally, it exhibits refractoriness of spindle oscillations, i.e. not every SO in a train of endogenous SOs triggers a spindle. Although adding extra currents increases dimensionality and parameter space, the model still preserves the overall simplicity and computational efficacy common to neural mass models.
Spindle timing
Relative to the negative deflection of a KC, spindles consistently start earlier than in the data. Consistently, the depolarizing up phase of endogenous KCs and SOs arrives earlier in the model than in the data. A comparison with the results from the isolated cortical module shows, that this is mostly due to strong depolarizing input from the thalamus. Yet, there is no clear explanation for the difference between model and experiment. It might be due to the simplification of the intrinsic mechanisms, e.g. firing rate adaptation in cortex and spindle dynamics in thalamus. On the other hand it could also be that finer details, e.g. spatial extension or the layered structure of the cortex are important for its temporal dynamics. Also the way conduction delays between cortex and thalamus were implemented, namely via an extra convolution with an alpha function, might play a role.
Auditory stimulation
A recent experimental study suggests that the refractoriness of thalamic spindles is a limiting factor for the impact of auditory stimulation upon memory consolidation [76]. They found, that longer trains of stimuli do not provide any benefit in memory consolidation compared to the two stimulus protocol. Remarkably, the first stimulus triggers a strong spindle, whereas the following stimuli show a diminished spindle response. This clearly indicates the importance of the grouping of spindles and SOs for the consolidation of memory. In contrast to these experimental findings, auditory stimulation in the model alleviates the refractoriness of the thalamic module, leading to spindle oscillations with similar amplitude following every stimulus. This is because strong depolarization of the thalamic populations by the stimulus interrupts the thalamic Ih rhythm. We see this as a challenge for the understanding of how auditory stimulation is processed during sleep and how it interacts with spindle generation.
Relation to other work
Recently, Cona et al. also developed a neural mass model to describe the sleeping thalamocortical system [79]. They combined two distinct firing modes via the activation of the T-type calcium current, showing that this multiplicative change in firing rate can lead to periodic spindle-like oscillations. However, in this study we include the currents directly into the equation of the membrane voltage, similar to [30, 54]. Our model relates directly to scalp EEG signals during natural sleep and auditory stimulation.
Effect of neuromodulators and sleep regulation
In our model, we induce the transition between the different sleep stages by changes of the three key parameters (gKNa and σp in the cortex and in the thalamus), that are directly linked to the action of neuromodulators [30, 80–82]. These parameters are known to be affected by neuromodulators, such as noradrenalin, serotonin and acetylcholine [44, 46, 50, 51, 83], whose concentrations vary over the night. Regulation of neuromodulator concentrations arises through complex interactions within different sleep regulatory networks [53, 84]. Recently there has been progress in the mathematical description of sleep regulatory networks [85–89]. However, as we focus on the different dynamical modes the thalamocortical system can exhibit and how thalamus and cortex interact, we do not include sleep regulation in this manuscript.
Are KCs biphasic or triphasic?
The waveform of a KC has been described as being biphasic, consisting of a large negative deflection (down state) followed by a pronounced depolarization (up state)—or triphasic, comprising an initial positive bump followed by a down state and an up state. Menicucci et al. [90] analyzed the shapes of KCs in N2 and N3 and found that on average a triphasic pattern, up-down-up, is present in both sleep stages. Our model does not show this sequence for sleep stage N2. In vivo, sleep stage N2 is rarely stationary and spans varying depths of sleep, as well as transitions to other sleep stages. In contrast, our model depicts idealized N2 at a single point in time, to separate it from wakefulness and N3. Choosing a parameter setting closer to N3 will naturally give rise to a depolarization preceding the down state. We predict that biphasic KCs should be found mostly in early N2 or very late N2, as in the second half of the night after the major SWS episodes.
Is there true bistability during natural sleep?
The model is consistent with the observation that during N2 and N3 of natural sleep the cortex is mostly in the active state [91]. We adopt the view of [30], where KCs were characterized as transient events—reversed spikes—initiated by a canard explosion. Consequently the down state is never stable in our model. This may seem counterintuitive as many intracellular recordings support the notion of bistability. However, neural mass models represent population averages, whereas intracellular recordings only sample individual members of a population, leaving open this alternative interpretation derived from stereotypical graphoelements in the EEG.
Supporting Information
S1 Text. Model equations.
This section provides the full mathematical description of the model presented here.
https://doi.org/10.1371/journal.pcbi.1005022.s001
(PDF)
S2 Text. Approximation of long range connection delay.
This section provides a justification of the approximation of the thalamocortical transmission delay by a convolution with an alpha function.
https://doi.org/10.1371/journal.pcbi.1005022.s002
(PDF)
S1 Table. Parameter values.
This table defines all constants used throughout this paper.
https://doi.org/10.1371/journal.pcbi.1005022.s003
(PDF)
S1 Dataset. Experimental data.
This dataset includes the experimental data used for Figs 6 and 7. Given are the time-locked averages for every subject. In Fig 6 only the sham condition was used to represent endogenous KCs/SOs, whereas Fig 7 displays the stimulation condition. For detailed information on data acquisition and processing please see the original study.
https://doi.org/10.1371/journal.pcbi.1005022.s004
(XLS)
Author Contributions
- Conceived and designed the experiments: MSC AW HVVN LM JB TM JCC.
- Performed the experiments: MSC AW HVVN.
- Analyzed the data: MSC AW HVVN LM JB TM JCC.
- Contributed reagents/materials/analysis tools: MSC AW HVVN LM JB TM JCC.
- Wrote the paper: MSC AW HVVN LM JB TM JCC.
References
- 1. Achermann P, Borbély AA. Low-frequency (<1 Hz) oscillations in the human sleep electroencephalogram. Neuroscience. 1997;81(1):213–22. pmid:9300413
- 2. Steriade M. The corticothalamic system in sleep. Frontiers in Bioscience. 2003;8(4):878–899.
- 3. Sanchez-Vives MV, McCormick DA. Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nature Neuroscience. 2000;3(10):1027–34. pmid:11017176
- 4. Peyrache A, Dehghani N, Eskandar EN, Madsen JR, Anderson WS, Donoghue Ja, et al. Spatiotemporal dynamics of neocortical excitation and inhibition during human sleep. Proceedings of the National Academy of Sciences of the United States of America. 2012;109(5):1731–6. pmid:22307639
- 5.
Timofeev I, Bazhenov M. Mechanisms and biological role of thalamocortical oscillations. Trends in Chronobiology Research. 2005;p. 1–47.
- 6. Mölle M, Marshall L, Gais S, Born J. Grouping of spindle activity during slow oscillations in human non-rapid eye movement sleep. The Journal of Neuroscience. 2002;22(24):10941–7. pmid:12486189
- 7. Walker MP, Stickgold R. Sleep-dependent learning and memory consolidation. Neuron. 2004;44(1):121–33. pmid:15450165
- 8. Rasch B, Born J. About sleep’s role in memory. Physiological Reviews. 2013;93(2):681–766. pmid:23589831
- 9. Rosanova M, Ulrich D. Pattern-specific associative long-term potentiation induced by a sleep spindle-related spike train. The Journal of Neuroscience. 2005;25(41):9398–9405. pmid:16221848
- 10. Mölle M, Bergmann TO, Marshall L, Born J. Fast and slow spindles during the sleep slow oscillation: disparate coalescence and engagement in memory processing. Sleep. 2011;34(10):1411–1421. pmid:21966073
- 11. Cox R, Hofman WF, Talamini LM. Involvement of spindles in memory consolidation is slow wave sleep-specific. Learning & Memory. 2012;19(7):264–267.
- 12. Marshall L, Helgadóttir H, Mölle M, Born J. Boosting slow oscillations during sleep potentiates memory. Nature. 2006;444(7119):610–3. pmid:17086200
- 13. Ngo HVV, Martinetz T, Born J, Mölle M. Auditory closed-loop stimulation of the sleep slow oscillation enhances memory. Neuron. 2013;78(3):545–53. pmid:23583623
- 14. Compte A, Sanchez-Vives MV, McCormick DA, Wang XJ. Cellular and network mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical network model. Journal of Neurophysiology. 2003;89(5):2707–25. pmid:12612051
- 15. Destexhe A, Bal T, McCormick DA, Sejnowski TJ. Ionic mechanisms underlying synchronized oscillations and propagating waves in a model of ferret thalamic slices. Journal of Neurophysiology. 1996;76(3):2049–70. pmid:8890314
- 16. Destexhe A, Neubig M, Ulrich D, Huguenard JR. Dendritic low-threshold calcium currents in thalamic relay cells. The Journal of Neuroscience. 1998;18(10):3574–88. pmid:9570789
- 17. Bazhenov M, Timofeev I. Cellular and network models for intrathalamic augmenting responses during 10-Hz stimulation. Journal of Neurophysiology. 1998;79:2730–2748. pmid:9582241
- 18. Destexhe A, Contreras D, Steriade M. Mechanisms underlying the synchronizing action of corticothalamic feedback through inhibition of thalamic relay cells. Journal of Neurophysiology. 1998;79(2):999–1016. pmid:9463458
- 19. Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ. Model of thalamocortical slow-wave sleep oscillations and transitions to activated States. The Journal of Neuroscience. 2002;22(19):8691–704. pmid:12351744
- 20. Destexhe A, Sejnowski TJ. Interactions between membrane conductances underlying thalamocortical slow-wave oscillations. Physiological Reviews. 2003;83(4):1401–53. pmid:14506309
- 21. Bonjean M, Baker T, Lemieux M, Timofeev I, Sejnowski TJ, Bazhenov M. Corticothalamic feedback controls sleep spindle duration in vivo. The Journal of Neuroscience. 2011;31(25):9124–34. pmid:21697364
- 22. Wilson HR, Cowan JD. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik. 1973;13(2):55–80. pmid:4767470
- 23. Lopes Da Silva FH, Hoeks A, Smits H, Zetterberg LH. Model of brain rhythmic activity. The alpha-rhythm of the thalamus. Kybernetik. 1974;15(1):27–37. pmid:4853232
- 24. Jansen BH, Zouridakis G, Brandt ME. A neurophysiologically-based mathematical model of flash visual evoked potentials. Electrical Engineering. 1993;283:275–283.
- 25. Kerr CC, Rennie CJ, Robinson PA. Physiology-based modeling of cortical auditory evoked potentials. Biological Cybernetics. 2008;98(2):171–84. pmid:18057953
- 26. Coombes S. Waves, bumps, and patterns in neural field theories. Biological Cybernetics. 2005;93(2):91–108. pmid:16059785
- 27. Deco G, Jirsa VK, Robinson PA, Breakspear M, Friston KJ. The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Computational Biology. 2008;4(8):e1000092. pmid:18769680
- 28. Spiegler A, Kiebel SJ, Atay FM, Knösche TR. Bifurcation analysis of neural mass models: Impact of extrinsic inputs and dendritic time constants. Neuroimage. 2010;52(3):1041–1058. pmid:20045068
- 29. Touboul J, Wendling F, Chauvel P, Faugeras O. Neural mass activity, bifurcations, and epilepsy. Neural Computation. 2011;23(12):3232–86. pmid:21919787
- 30. Weigenand A, Schellenberger Costa M, Ngo HVV, Claussen JC, Martinetz T. Characterization of K-Complexes and Slow Wave Activity in a Neural Mass Model. PLoS Computational Biology. 2014;10:e1003923. pmid:25392991
- 31. Robinson PA, Rennie CJ, Wright JJ. Propagation and stability of waves of electrical activity in the cerebral cortex. Physical Review E. 1997;56(1):826–840.
- 32. Liley DTJ, Cadusch PJ, Wright JJ. A continuum theory of electro-cortical activity. Neurocomputing. 1999;26–27(May 2000):795–800.
- 33. Liley DTJ, Cadusch PJ, Dafilis MP. A spatially continuous mean field theory of electrocortical activity. Network (Bristol, England). 2002;13(1):67–113.
- 34. Wilson MT, Sleigh JW, Steyn-Ross DA, Steyn-Ross ML. General anesthetic-induced seizures can be explained by a mean-field model of cortical dynamics. Anesthesiology. 2006;104(3):588–93. pmid:16508406
- 35. Marreiros AC, Daunizeau J, Kiebel SJ, Friston KJ. Population dynamics: variance and the sigmoid activation function. NeuroImage. 2008;42(1):147–157. pmid:18547818
- 36. Zandt BJ, Visser S, van Putten MJ, ten Haken B. A neural mass model based on single cell dynamics to model pathophysiology. Journal of Computational Neuroscience. 2014;37(3):549–568. pmid:25131270
- 37. Buzsáki G, Anastassiou CA, Koch C. The origin of extracellular fields and currents–EEG, ECoG, LFP and spikes. Nature reviews Neuroscience. 2012;13(6):407–20. pmid:22595786
- 38. Steyn-Ross DA, Steyn-Ross ML, Wilcocks L, Sleigh JW. Toward a theory of the general-anesthetic-induced phase transition of the cerebral cortex. II. Numerical simulations, spectral entropy, and correlation times. Physical Review E. 2001;64(1):1–12.
- 39. Sotero RC, Trujillo-Barreto NJ. Biophysical model for integrating neuronal activity, EEG, fMRI and metabolism. NeuroImage. 2008;39(1):290–309. pmid:17919931
- 40. Langdon AJ, Breakspear M, Coombes S. Phase-locked cluster oscillations in periodically forced integrate-and-fire-or-burst neuronal populations. Physical Review E. 2012;86(6):061903.
- 41. Benita JM, Guillamon A, Deco G, Sanchez-Vives MV. Synaptic depression and slow oscillatory activity in a biophysical network model of the cerebral cortex. Frontiers in Computational Neuroscience. 2012;6:64. pmid:22973221
- 42. McCormick DA. Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. Progress in Neurobiology. 1992;39(4):337–388. pmid:1354387
- 43. Mittmann T, Alzheimer C. Muscarinic inhibition of persistent Na+ current in rat neocortical pyramidal neurons. Journal of Neurophysiology. 1998;79:1579–1582. pmid:9497434
- 44. Timmons D, Geisert E, Stewart E, Lorenzon M, Foehring C. Alpha 2 -Adrenergic receptor-mediated modulation of calcium current in neocortical pyramidal neurons. Brain Research. 2004;1014:184–196. pmid:15213003
- 45. Mehaffey W, Doiron B. Deterministic multiplicative gain control with active dendrites. The Journal of Neuroscience. 2005;25(43):9968–9977. pmid:16251445
- 46. Zhang ZW, Arsenault D. Gain modulation by serotonin in pyramidal neurones of the rat prefrontal cortex. The Journal of Physiology. 2005;2:379–394.
- 47. Hasselmo ME, Giocomo LM. Cholinergic modulation of cortical function. Journal of Molecular Neuroscience. 2006;30(1):133–135. pmid:17192659
- 48. Disney A, Aoki C, Hawken M. Gain modulation by nicotine in macaque v1. Neuron. 2007;56:701–713. pmid:18031686
- 49. Thurley K, Senn W, Lüscher H. Dopamine increases the gain of the input-output response of rat prefrontal pyramidal neurons. Journal of Neurophysiology. 2008;99:2985. pmid:18400958
- 50. Gulledge A, Bucci D. M1 receptors mediate cholinergic modulation of excitability in neocortical pyramidal neurons. The Journal of Neuroscience. 2009;29(31):9888–9902. pmid:19657040
- 51. Soma S, Shimegi S. Cholinergic modulation of response gain in the primary visual cortex of the macaque. Journal of Neurophysiology. 2012;107:283–291. pmid:21994270
- 52. Polack PO, Friedman J, Golshani P. Cellular mechanisms of brain state-dependent gain modulation in visual cortex. Nature Neuroscience. 2013;16(9):1331–1339. pmid:23872595
- 53. Léna I, Parrot S, Deschaux O, Muffat-Joly S, Sauvinet V, Renaud B, et al. Variations in extracellular levels of dopamine, noradrenaline, glutamate, and aspartate across the sleep–wake cycle in the medial prefrontal cortex and nucleus accumbens of freely moving rats. Journal of Neuroscience Research. 2005;81(6):891–9. pmid:16041801
- 54. Żygierewicz J, Suffczyński P, Blinowska K. A model of sleep spindles generation. Neurocomputing. 2001;38–40:1619–1625.
- 55. Suffczyński P, Kalitzin S, Lopes Da Silva FH. Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience. 2004;126(2):467–84. pmid:15207365
- 56. Destexhe A, Contreras D. In vivo, in vitro, and computational analysis of dendritic calcium currents in thalamic reticular neurons. The Journal of Neuroscience. 1996;16(1):169–185. pmid:8613783
- 57. Chen JY, Chauvette S, Skorheim S, Timofeev I, Bazhenov M. Interneuron-mediated inhibition synchronizes neuronal activity during slow oscillation. The Journal of Physiology. 2012;590(Pt 16):3987–4010. pmid:22641778
- 58. Jirsa VK. Neural field dynamics with local and global connectivity and time delay. Philosophical Transactions Series A, Mathematical, Physical, and Engineering Sciences. 2009;367(February):1131–1143. pmid:19218155
- 59. Nakagawa TT, Woolrich M, Luckhoo H, Joensson M, Mohseni H, Kringelbach ML, et al. How delays matter in an oscillatory whole-brain spiking-neuron network model for MEG alpha-rhythms at rest. NeuroImage. 2014;87:383–394. pmid:24246492
- 60. Agmon A, Connors BW. Correlation between intrinsic firing patterns and thalamocortical synaptic responses of neurons in mouse barrel cortex [Journal Article]. The Journal of Neuroscience. 1992;12(1):319–329. pmid:1729440
- 61. Salami M, Itami C, Tsumoto T, Kimura F. Change of conduction velocity by regional myelination yields constant latency irrespective of distance between thalamus and cortex. Proceedings of the National Academy of Sciences of the United States of America. 2003;100(10):6174–9. pmid:12719546
- 62. Gentet LJ, Ulrich D. Electrophysiological characterization of synaptic connections between layer VI cortical cells and neurons of the nucleus reticularis thalami in juvenile rats [Journal Article]. European Journal of Neuroscience. 2004;19(3):625–633. pmid:14984412
- 63. Traub RD, Contreras D, Cunningham MO, Murray H, LeBeau FEN, Roopun A, et al. Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. Journal of Neurophysiology. 2005;93(4):2194–232. pmid:15525801
- 64. Biggio M, Storace M, Mattia M. Non-instantaneous synaptic transmission in spiking neuron networks and equivalence with delay distribution. BMC Neuroscience. 2013;14(Suppl 1):P267.
- 65. Ngo HVV, Claussen JC, Born J, Mölle M. Induction of slow oscillations by rhythmic acoustic stimulation. Journal of Sleep Research. 2013;22(1):22–31. pmid:22913273
- 66. Rößler A. Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations. SIAM Journal on Numerical Analysis. 2010;48(3):922–952.
- 67.
Schellenberger Costa M. Neural mass model of the isolated thalamus during sleep.; 2015. Available from: https://github.com/miscco/NM_Thalamus.
- 68.
Schellenberger Costa M. Neural mass model of the thalamocortical system during sleep.; 2015. Available from: https://github.com/miscco/NM_TC.
- 69. Lüthi A, McCormick DA. Modulation of a pacemaker current through Ca(2+)-induced stimulation of cAMP production. Nature Neuroscience. 1999;2(7):634–641. pmid:10404196
- 70. Lüthi A, McCormick DA. Periodicity of thalamic synchronized oscillations: the role of Ca2+- mediated upregulation of Ih. Neuron. 1998;20(3):553–563. pmid:9539128
- 71. Contreras D, Destexhe a, Sejnowski TJ, Steriade M. Spatiotemporal patterns of spindle oscillations in cortex and thalamus. The Journal of Neuroscience. 1997;17(3):1179–1196. pmid:8994070
- 72. Shilnikov A, Cymbalyuk G. Transition between tonic spiking and bursting in a neuron model via the blue-sky catastrophe. Physical Review Letters. 2005;94(4):2–5.
- 73. Mayer J, Schuster HG, Claussen JC. Role of inhibitory feedback for information processing in thalamocortical circuits. Physical Review E. 2006;73(3).
- 74. Contreras D, Steriade M. Cellular basis of EEG slow rhythms: a study of dynamic corticothalamic relationships. The Journal of Neuroscience. 1995;15(1):604–622. pmid:7823167
- 75. Andrillon T, Nir Y, Staba RJ, Ferrarelli F, Cirelli C, Tononi G, et al. Sleep spindles in humans: insights from intracranial EEG and unit recordings. The Journal of Neuroscience. 2011;31(49):17821–17834. pmid:22159098
- 76. Ngo HV, Miedema A, Faude I, Martinetz T, Mölle M, Born J. Driving sleep slow oscillations by auditory closed-loop stimulation—a self-limiting process. Accepted for publication in Journal of Neuroscience. 2015;.
- 77. Mölle M, Yeshenko O, Marshall L, Sara SJ, Born J. Hippocampal sharp wave-ripples linked to slow oscillations in rat slow-wave sleep. Journal of Neurophysiology. 2006;96:62–70. pmid:16611848
- 78. Diekelmann S, Born J. The memory function of sleep. Nature Reviews Neuroscience. 2010;11(2):114–26. pmid:20046194
- 79. Cona F, Lacanna M, Ursino M. A thalamo-cortical neural mass model for the simulation of brain rhythms during sleep. Journal of Computational Neuroscience. 2014;37:125–148. pmid:24402459
- 80. McCormick DA, von Krosigk M. Corticothalamic activation modulates thalamic firing through glutamate “metabotropic” receptors. Proceedings of the National Academy of Sciences of the United States of America. 1992;89:2774–2778. pmid:1313567
- 81. Hughes SW, Cope DW, Blethyn KL, Crunelli V. Cellular mechanisms of the slow (<1 Hz) oscillation in thalamocortical neurons in vitro. Neuron. 2002;33(6):947–58. pmid:11906700
- 82. Steriade M. Acetylcholine systems and rhythmic activities during the waking-sleep cycle. Progress in Brain Research. 2004;145:179–196. pmid:14650916
- 83. McCormick DA. Cholinergic and noradrenergic modulation of thalamocortical processing. Trends in Neurosciences. 1989;12(6):215–221. pmid:2473557
- 84. Lydic R, Baghdoyan HA. Sleep, anesthesiology, and the neurobiology of arousal state control. Anesthesiology. 2005;103(6):1268–1295. pmid:16306742
- 85. Tamakawa Y, Karashima A, Koyama Y, Katayama N, Nakao M. A quartet neural system model orchestrating sleep and wakefulness mechanisms. Journal of Neurophysiology. 2006;95(4):2055–2069. pmid:16282204
- 86. Diniz Behn CG, Booth V. Simulating microinjection experiments in a novel model of the rat sleep-wake regulatory network. Journal of Neurophysiology. 2010;103(4):1937–1953. pmid:20107121
- 87. Phillips AJK, Robinson PA. A quantitative model of sleep-wake dynamics based on the physiology of the brainstem ascending arousal system. Journal of Biological Rhythms. 2007;22(2):167–179. pmid:17440218
- 88. Rempe MJ, Best J, Terman D. A mathematical model of the sleep/wake cycle. Journal of Mathematical Biology. 2009;60(5):615–644. pmid:19557415
- 89. Kumar R, Bose A, Mallick BN. A mathematical model towards understanding the mechanism of neuronal regulation of wake-NREMS-REMS states. PLoS ONE. 2012;7(8):1–17.
- 90. Menicucci D, Piarulli A, Allegrini P, Laurino M, Mastorci F, Sebastiani L, et al. Fragments of wake-like activity frame down-states of sleep slow oscillations in humans: New vistas for studying homeostatic processes during sleep. International Journal of Psychophysiology. 2013;89(2):151–157. pmid:23384886
- 91. Chauvette S, Crochet S, Volgushev M, Timofeev I. Properties of Slow Oscillation during Slow-Wave Sleep and Anesthesia in Cats. The Journal of Neuroscience. 2011;31(42):14998–15008. pmid:22016533