Next Article in Journal
Associations of Hyperactivity and Inattention Scores with Theta and Beta Oscillatory Dynamics of EEG in Stop-Signal Task in Healthy Children 7–10 Years Old
Next Article in Special Issue
Editorial to the Special Issue “Information Processing in Neuronal Circuits and Systems”
Previous Article in Journal
Metabolomics in Retinal Diseases: An Update
Previous Article in Special Issue
A Model of the Early Visual System Based on Parallel Spike-Sequence Detection, Showing Orientation Selectivity
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Modeling of Information Propagation during the Sleep–Waking Cycle

by
Farhad Razi
1,2,
Rubén Moreno-Bote
2 and
Belén Sancristóbal
1,*
1
Barcelona School of Design and Engineering, Universitat de Vic—Universitat Central de Catalunya, 08002 Barcelona, Spain
2
Center for Brain and Cognition, Department of Information and Communications Technologies, Universitat Pompeu Fabra, 08018 Barcelona, Spain
*
Author to whom correspondence should be addressed.
Submission received: 25 July 2021 / Revised: 9 September 2021 / Accepted: 16 September 2021 / Published: 22 September 2021
(This article belongs to the Special Issue Information Processing in Neuronal Circuits and Systems)

Abstract

:

Simple Summary

During the deep phases of sleep we do not normally wake up by a thunder, but we nevertheless notice it when awake. The exact same sound gets to our ears and cortex through the thalamus and still, it triggers two very different responses. There is growing experimental evidence that these two states of the brain—sleep and wakefulness—distribute sensory information in different ways across the cortex. In particular, during sleep, neuronal responses remain local and do not spread out across distant synaptically connected regions. On the contrary, during wakefulness, stimuli are able to elicit a wider spatial response. We have used a computational model of coupled cortical columns to study how these two propagation modes arise. Moreover, the transition from sleep-like to waking-like dynamics occurs in agreement with the synaptic homeostasis hypothesis and only requires the increase of excitatory conductances. We have found that, in order to reproduce the aforementioned observations, this parameter change has to be selectively applied: synaptic conductances between distinct columns have to be potentiated over local ones.

Abstract

Non-threatening familiar sounds can go unnoticed during sleep despite the fact that they enter our brain by exciting the auditory nerves. Extracellular cortical recordings in the primary auditory cortex of rodents show that an increase in firing rate in response to pure tones during deep phases of sleep is comparable to those evoked during wakefulness. This result challenges the hypothesis that during sleep cortical responses are weakened through thalamic gating. An alternative explanation comes from the observation that the spatiotemporal spread of the evoked activity by transcranial magnetic stimulation in humans is reduced during non-rapid eye movement (NREM) sleep as compared to the wider propagation to other cortical regions during wakefulness. Thus, cortical responses during NREM sleep remain local and the stimulus only reaches nearby neuronal populations. We aim at understanding how this behavior emerges in the brain as it spontaneously shifts between NREM sleep and wakefulness. To do so, we have used a computational neural-mass model to reproduce the dynamics of the sensory auditory cortex and corresponding local field potentials in these two brain states. Following the synaptic homeostasis hypothesis, an increase in a single parameter, namely the excitatory conductance g ¯ AMPA , allows us to place the model from NREM sleep into wakefulness. In agreement with the experimental results, the endogenous dynamics during NREM sleep produces a comparable, even higher, response to excitatory inputs to the ones during wakefulness. We have extended the model to two bidirectionally connected cortical columns and have quantified the propagation of an excitatory input as a function of their coupling. We have found that the general increase in all conductances of the cortical excitatory synapses that drive the system from NREM sleep to wakefulness does not boost the effective connectivity between cortical columns. Instead, it is the inter-/intra-conductance ratio of cortical excitatory synapses that should raise to facilitate information propagation across the brain.

1. Introduction

Sensory information processing is not fully blocked during sleep since animals need to identify relevant stimuli for their survival. However, when compared to wakefulness, focused attention to the environment is lost while we sleep. Such reduced sensitivity in the detection of external signals leads to a decrease in sensory awareness. Still, sleep is a vital recurrent state that reduces brain energy demands by slowing down the metabolic rate [1]. Behaviorally, this is clearly observed in the cessation of actions since sleep holds back reproduction, exploration, protection and nurture of the offspring. Electrophysiologically, the deep phases of sleep (NREM sleep) are distinguished from wakefulness based on the electrical activity emerging from neuronal populations. For instance, extracellular cortical recordings obtained during wakefulness show low-amplitude–high-frequency voltage fluctuations whereas during NREM sleep, and also during deep anesthesia [2,3,4], these fluctuations exhibit high-amplitude–low-frequency components and are bimodal [3,4,5,6,7,8]. The bimodality appears due to the alternation between periods of persistent firing activity (Up states) and silent periods (Down states), which provides the idiosyncratic dynamical features of NREM sleep, also known as slow wave activity (SWA). This leads to a reduction of the overall firing rate of both excitatory and inhibitory neuronal cell types in the time course of NREM sleep [4,5,6,7,8].
Despite these differences, neurons in the primary auditory cortex show comparable responses during sleep and wakefulness when pure tones are delivered to rats [7] and primates [9]. These experimental results challenge the hypothesis that during sleep cortical responses are weakened through thalamic gating, i.e., the mechanism by which thalamic spindles could prevent a reliable transmission of external signals to sensory cortices during sleep [10]. Interestingly, in higher-order cortical areas during sleep—in the perirhinal cortex in rats [8] and the human cingulate and prefrontal cortex [11] and during anesthesia—in the association cortex in humans [12], auditory responses are attenuated. These studies imply that information gating must occur within brain structures higher up the stream of sensory processing and seem to be more consistent with a breakdown of functional cortical connectivity rather than of thalamocortical connectivity. Indeed, fMRI data collected in humans indicate a reduced global effective connectivity during NREM sleep [13]. Even in the presence of cycling alternating patterns (CAP), a sign of sleep instability, EEG-based functional connectivity shows a small-world network type of structure [14]. Moreover, cortical responses to transcranial magnetic stimulation (TMS) have been shown to be local during NREM sleep and do not propagate to connected brain regions [15,16].
One possible neural correlate of decreased sensory responsiveness in the non-primary cortex during NREM sleep could be the downscaling of the conductance of excitatory synapses, a hypothesis known as synaptic homeostasis hypothesis (SHY) [17,18,19,20,21]. Experimental evidence [22,23,24,25] shows that potentiation of synapses during wakefulness caused by synaptic plasticity is compensated during NREM sleep with a general reduction of synaptic strength. We implement this synaptic conductance downscaling in our computational study by tuning a single parameter that is biophysically equivalent to an excitatory synaptic conductance. Moreover, decreasing its value shifts the model from wakefulness to NREM sleep. To our knowledge, this is the first mathematical model, among all the different proposals [26,27,28,29,30], that reconciles these two important findings of sleep research: synaptic downscaling and the simultaneous emergence of slow rhythms as the brain enters the deep phases of sleep from wakefulness.
We aim at understanding, within the neural mass model formalism used in this study and others [27,31], whether SHY not only explains the aforementioned change in brain dynamics, but also in the propagation span of evoked activities occurring throughout the sleep–waking cycle (SWC). To do so, we have simulated the local field potential (LFP) signals of two cortical columns operating in either NREM sleep or wakefulness. One cortical column—the perturbed one—responds to an externally applied excitatory input, reproducing the activation of a primary auditory area by a sound or a superficial region of the cerebral cortex stimulated by TMS. A second cortical column—the unperturbed one—is bidirectionally coupled to the perturbed column and does not directly receive the external input, as occurs with non-primary auditory areas and neighboring regions. We quantify the propagation of this external input as a function of the strength of the inter-column connectivity in each brain state. We show that the upscaling of the conductance of excitatory synapses shifts the LFP dynamics from NREM sleep to wakefulness in both columns, as they are intrinsically identical. However, a general build-up of this conductance does not lead to the propagation of the external input from the perturbed cortical column to the unperturbed cortical column. Interestingly, we have found that the unperturbed cortical column responds only significantly when the increase in the excitatory conductance selectively takes place in the inter-column synapses with respect to the intra-column connectivity. This finding points to a novel and different interpretation of SHY, namely that the tuning of synaptic conductances could be distance-dependent.
Our results are presented as follows. First, we introduce the model for a single cortical column in order to validate the parameter space where the system exhibits NREM and wakefulness dynamics. Second, we extend the network to two cortical columns coupled bidirectionally with excitatory synapses to generalize the results to a more complex connectivity diagram. Note that now the change in synaptic conductance that allows the system to shift dynamics, also affects the inter-column connectivity and thus directly influences signal transmission between columns. Finally, we show how efficient input propagation observed during wakefulness, and reduced during NREM sleep, can be obtained.

2. Materials and Methods

In the following section, we first describe the neural mass formalism [27,31] that is able to reproduce the extracellular dynamics of a single cortical column. Increasing the strength of the intra-cortical excitatory synapses, in agreement with SHY, allows us to move from NREM-like activity to wakefulness. Finally, we expand the model to two mutually connected cortical columns.

2.1. Neural Mass Model

Recordings of the extracellular neuronal field capture the superposition of all transmembrane currents, such as the synaptic and ionic currents [32]. In order to reproduce the dynamical features of these signals, we use a neural mass model [27,31] representing the average membrane potential of a neural cortical column that leads to the LFP characteristics of the SWC. A single cortical column consists of pyramidal, p, and inhibitory, i, population mutually connected (see Figure 1A for schematic).
The dynamics of the average membrane potential, V p / i , is governed by the mathematical formalism of the classical conductance-based model [33] with one leak, two synaptic currents and one activity-dependent potassium current as follows:
τ p V ˙ p = I L p I AMPA p I GABA p τ p C m 1 I KNa ,
τ i V i ˙ = I L i I AMPA i I GABA i ,
where τ p / i , I L p / i , I A M P A p / i and I G A B A p / i are, respectively, the membrane time constant and the leak, AMPAergic and GABAergic current of population p or i. I KNa is the activity-dependent potassium current of the pyramidal population (see Section 2.3). C m is the membrane capacitance in the Hodgkin–Huxley model.
The average firing rate of each population is modeled as an instantaneous function of the average membrane potential V p / i according to a sigmoid function [34,35]:
Q p / i ( V p / i ) = Q p / i m a x ( 1 + tanh ( C ( V p / i θ p / i ) / σ p / i ) ) 2 ,
where Q p / i m a x , θ p / i and σ p / i represent the maximal firing rate, firing rate threshold and inverse neural gain of population p or i, respectively (see Table 1 and Table 2 for symbol description and parameter value, respectively). Here, C = π 2 3 is a constant linking the neural gain to the slope of the sigmoid function [27].

2.2. Model Synapses

The presynaptic activities have an intrinsic and an extrinsic origin. The intrinsic or recurrent activity is the input each population receives from local pyramidal and inhibitory populations. The extrinsic input comes from non-explicitly modeled distant cortical columns that impinges on population p and i only through excitatory synapses. This is in agreement with the morphology of cortical neurons, whereby pyramidal neurons have larger projecting axons with respect to the more local and constrained range of inhibitory axons [36,37,38]. This has also been reflected in numerous computational studies [30,39,40,41,42].
Given a presynaptic population k and a postsynaptic population k, the time course of synaptic activity, s k k , at time t, depends on the pyramidal or inhibitory nature of the presynaptic population k . Synapses are considered AMPAergic for excitation, i.e., when population k is the pyramidal one, and GABAergic for inhibition, i.e., when population k is the inhibitory one. The mathematical formalism adopted is the convolution of the presynaptic firing, Q k , with the average synaptic response to a single spike α k that has an exponential decay time course [43]:
s k k = N k k Q k ( V k ) α k ,
α k = γ k 2 t exp ( γ k t ) ,
where N k k is the mean number of synaptic connections from presynaptic population k to postsynaptic population k. γ k describes the time constant of the dynamics of a synapse activated by the presynaptic firing of population k . The set of Equations (4) and (5) lead to the second-order differential equation (derived in detail in Appendix A.1):
s ¨ k k = γ k 2 ( N k k Q k ( V k ) + ϕ s k k ) 2 γ k s ˙ k k ,
where Equation (6) describes the AMPAergic and GABAergic synaptic activity due to the firing of the pyramidal and inhibitory population k , respectively, onto postsynaptic population k. Moreover, the noise ϕ is simulated independently for each cortical population as a Gaussian process with zero autocorrelation time constant and zero mean. We assume that the standard deviation of ϕ is greater during NREM, i.e., 1.8 ms 1 , than during wakefulness, i.e., 1 ms 1 , because in the former dynamics, the firing rate fluctuates largely between the Up and Down states. The AMPAergic and GABAergic synaptic currents are defined as:
I AMPA k = g ¯ AMPA k s k p ( V k E AMPA ) ,
I GABA k = g ¯ GABA k s k i ( V k E GABA ) ,
where g ¯ AMPA k and g ¯ GABA k are, respectively, the average AMPAergic and GABAergic conductance on population k (also referred to as the recurrent or intra-conductances). E AMPA and E GABA are the reversal potential of AMPAergic and GABAergic currents.

2.3. Activity-Dependent Potassium Current

SWA is the hallmark of neuronal cortical activity during NREM sleep and anesthesia. It consists of silent (Down) and persistent (Up) firing patterns of activity that alternate. Activation of the slow activity-dependent K + current has been suggested as one possible mechanism triggering Down state initiation [3,44,45,46]. The sodium-dependent potassium current I KNa and sodium concentration [Na] are implemented in the model as:
I KNa = g ¯ KNa 0.37 1 + ( 38.7 [ Na ] ) 3.5 ( V p E K ) ,
τ Na [ Na ] ˙ = α Na Q p ( V p ) Na pump ( [ Na ] ) ,
where g ¯ KNa , E K , τ Na and α Na correspond to the average conductance of I KNa , the Nernst reversal potential of the I KNa current, the time constant of extrusion of sodium concentration and the average influx of sodium concentration per firing. Na pump is a function representing sodium pumps that determines the extrusion of sodium concentration through sodium pumps (see Appendix A for details).

2.4. Upscaling the AMPAergic Conductance

The neural mass model proposed by [27,31] is able to generate the physiological characteristics of the sleeping cortex for particular values of the parameters. There, the authors showed that by reducing σ p and g ¯ KNa , the model generates the low-amplitude–high-frequency fluctuations in the average membrane potential characteristic of wakefulness. However, with their proposed set of parameter values, the firing rate of the pyramidal population is saturated at its maximum, given by Q p max of the sigmoid function Q p ( V p ) . Therefore, we propose a different set of bifurcation parameters to place the model within the wakefulness regime.
We have chosen to upscale the conductance of excitatory synapses by increasing the average AMPAergic conductance, g ¯ AMPA k , to go from NREM sleep to wakefulness, which is in agreement with SHY. The average GABAergic conductance, g ¯ GABA k , has been modified on both pyramidal and inhibitory populations in order to maintain the steady state value of the average membrane potential of pyramidal V p and inhibitory V i populations equal to their values in the Up state during NREM sleep (see Table 2 to compare parameter values between brain states). These V p and V i values are determined by simulating 500 trials (10 s long in duration) of NREM dynamics. Up and Down events are determined from the LFP signals (see Section 2.8.1 below for a definition of LFP and Section 2.8.2 for the Up and Down separation method). Then, the peak of the distributions of V p and V i during the Up state events is used as the steady state value of the average membrane potential of pyramidal and inhibitory populations during wakefulness. The average GABAergic conductance on pyramidal and inhibitory populations is increased until these V p and V i values are obtained for the g ¯ AMPA k values responsible of wakefulness (see Appendix B.1 for dynamical constraints on the average GABAergic conductance on pyramidal and inhibitory populations in the one-cortical-column scenario).

2.5. External Input

The pyramidal population is the only one subject to an external input ξ (see Appendix A.2 for full model equations). Hence, it will be termed perturbed cortical population from here on and its cortical column will be termed perturbed cortical column. The external input is considered to be a transient 100 ms increase in the mean of the white noise ϕ (see Table 2 for parameter values). It is modeled as a square function (see Figure 1) of 1 ms 1 amplitude.

2.6. Model Extension to Two Cortical Columns

We have extend the model to two bidirectionally connected cortical columns (see Figure 1B). Besides the perturbed cortical column, now we model a second cortical column that does not directly receive the external input ξ and will be from here on termed unperturbed cortical column and its pyramidal population will be termed unperturbed cortical population. The rest of parameters are the same for the perturbed and unperturbed cortical columns. We consider that the connectivity between the two cortical columns is symmetric and excitatory due to the longer axons of pyramidal neurons as compared to the short-range axons of inhibitory neurons [30,36,37,38,39,40,41]. The average AMPAergic conductance of excitatory synapses between the two cortical columns will be referred here as inter-conductance, as opposed to the intra-conductance of each single cortical column. The average GABAergic conductance on pyramidal and inhibitory populations are here increased to counterbalance this inter-cortical column coupling that introduces more excitation to the pyramidal and inhibitory populations. By doing so, we keep the steady state value of V p and V i in the two-cortical-column scenario equal to the values in the one-cortical-column scenario both during the Up states of NREM sleep and wakefulness (see Table 2 for parameter values and Appendix B.2 for the dynamical constraints on the increased average GABAergic conductance during NREM sleep and wakefulness in the two-cortical-column model). Therefore, the second-order differential equations for AMPAergic postsynaptic responses in each cortical column due to the inter-cortical column coupling are as follows:
s ¨ p p inter = γ p 2 ( N p p Q p ( V p ) s p p inter ) 2 γ p s ˙ p p inter ,
s ¨ i p inter = γ p 2 ( N i p Q p ( V p ) s i p inter ) 2 γ p s ˙ i p inter ,
where N k p is the mean number of synaptic connections from a presynaptic pyramidal population p in one cortical column to postsynaptic populations k = p , i in the other cortical column (see Appendix A.3 for full model equations).
We introduce the parameter β as the ratio of the average AMPAergic inter-conductance to the average AMPAergic intra-conductance. Therefore, the AMPAergic current to the postsynaptic population k in each cortical column is as follows:
I AMPA k = g ¯ AMPA k s k p ( V k E AMPA ) + β · g ¯ AMPA k s k p inter ( V k E AMPA ) ,
where k is either the pyramidal or inhibitory postsynaptic population of either cortical column.
We quantify information propagation to the unperturbed population as a function of β and g ¯ AMPA k . First, we increase β from 2 to 5 in steps of one unit while g ¯ AMPA k = 2 is kept constant. To counterbalance the increase in excitation to pyramidal and inhibitory populations due to these changes in β , we increase the average GABAergic conductance on pyramidal and inhibitory populations (see Table 4 for parameter values).
Next, we repeat the same approach for two other values of g ¯ AMPA k = 6 and 10. In each case, β is varied within the aforementioned values. Again, to counterbalance the increase in excitation to the pyramidal and inhibitory populations, we increase the average GABAergic conductance on pyramidal and inhibitory populations of both cortical columns (seeTable 5 and Table 6 for g ¯ GABA k values at g ¯ AMPA k = 6 and 10, respectively).

2.7. Numerical Integration

The model was implemented and run in Python, using a stochastic Heun method [47] with a step size of 0.1 ms. The code is available at github [48]. Each simulation was run independently, choosing the initial conditions of all variables at random from a uniform distribution. Averages are obtained across 500 independent trials of length 10 s after discarding the initial 10 s from the simulation onset to eliminate transient dynamics towards the stable solution.

2.8. Data Analysis

All data analyses were carried out offline with Python. We simulated six signals from each cortical column: LFP, firing rate and average membrane potential signals from both pyramidal and inhibitory populations.

2.8.1. LFP Definition

The LFP signal is defined as the sum of the absolute value of AMPAergic and GABAergic synaptic currents to cortical population k of either cortical column [49,50]. The action potentials are fast events of less than 10 ms that are attenuated by the cortical tissue and, therefore, do not contribute to the LFP signals.
LFP = | I AMPA k | + | I GABA k | .

2.8.2. Analysis of Amplitude-Frequency Content of LFP Signal

In each trial, we removed the temporal mean of the LFP of prestimulus intervals of either cortical column (from the beginning of simulation up to five seconds later). We used the prestimulus intervals of either cortical column to analyze the amplitude-frequency content during NREM and wakefulness states and the poststimulus intervals (from external input onset up to five seconds later) to analyze the evoked responses in either cortical column due to the external input.
To asses the spread of the spontaneous LFP values, we computed the distribution of the LFP during the prestimulus interval of both cortical columns for each brain state separately with a bin size of 0.5 mV. To locate the peaks of the bimodal LFP distribution during NREM sleep, the distribution for pyramidal (inhibitory) population was first smoothed with a Gaussian kernel, exp ( x 2 2 c 2 ) , where c is equal to 5 5 mV (5 mV). The half distance between the Up and Down peaks in the LFP distribution during NREM sleep was used as a threshold to classify these events in NREM sleep. Events with the LFP signal higher than the threshold were classified as Up states, otherwise they were labeled as Down states. Then, the distribution of the firing rate, in 0.5 Hz bins, and the average membrane potential, in 0.5 mV bins, were computed separately for the Up and Down states. Distributions were normalized by the total number of events during NREM sleep. Besides, the distribution of the firing rate and the average membrane potential during wakefulness were computed and normalized in the same way.
In order to compute the frequency content of the LFP signal during NREM sleep and wakefulness, we computed its spectrogram. To do so, we first windowed the LFP signals in 2 s length intervals with a 90% overlap within which the LFP is considered quasi-stationary [51]. Then, each window was tapered with a Hann function to control for the spectral leakage, and we computed the discrete Short-Time Fourier Transform (STFT). Next, the power spectral density (PSD) on each time window was obtained from the amplitude of the STFT. Because NREM sleep and wakefulness were independently simulated and were not time-locked events, we averaged the PSD over all time windows [52]. To better observe the wide dynamic range of the LFP, we plotted the spectrogram in units of logarithmic decibels (dB), defined as:
dB = log 10 P P r ,
where P is the averaged PSD and P r the reference PSD, set at 1 ( mV ) 2 / Hz . Finally, the spectrogram of the LFP signals were averaged over 500 simulations for each brain state.
The distribution of the log ratio between high (>30 Hz) and low (<4 Hz) frequency bands of the prestimulus intervals was computed to further verify the different dynamics between NREM and wakefulness.

2.8.3. Statistical Test for Significance of the Evoked Response

We tested the significance of the evoked responses in the poststimulus intervals in the one-cortical-column and the two-cortical-column scenario. First, we tested whether the postistimulus signals differed significantly from the spontaneous activity captured in the prestimulus intervals during NREM sleep and wakefulness, separately. This test allowed us to quantify the emergence of a response to the external input. Second, we tested whether the postistimulus signals differed significantly between brain states. This test allowed us to quantify the sensitivity of the response to the external input to the spontaneous dynamics.
In the case of the two cortical columns scenario, we verified whether the postistimulus signals differed significantly from the spontaneous activity in the perturbed population to maintain the same conditions of the single cortical column. For the unperturbed population, the same test allowed us to quantify the propagation of the response.
These significance tests were run independently for all signals in the one- and two-cortical-column scenario. The statistical procedure was carried out using a temporal clustering nonparametric test [53] to control for family-wise error in multiple comparisons of temporal t-tests. Since the number of simulations was large (500 simulations) for the Central Limit Theorem to hold, we computed an independent t-test value (t-value) at every time point as t = ( x ¯ 1 x ¯ 2 ) / ( s 2 / n ) , where x ¯ 1 and x ¯ 2 are the averages across simulations of the signals that are being compared against the null hypothesis of no statistical difference between them. s is the sample standard deviation, whose variance is s 2 = ( ( n 1 ) s 1 2 + ( n 1 ) s 2 2 ) / ( 2 n 1 ) . s 1 / 2 2 is the variance of variable x 1 / 2 across simulations. Statistical significance is set at p-values below 0.01, which corresponds to a two-tailed critical t-value = ± 2.58 , for n = 500 . We computed the area of all clusters made of consecutive time points (referred as t-cluster statistic) where x ¯ 1 x ¯ 2 2.58 and sorted them in descending order. Clusters with an area below the largest t-cluster statistic computed from the prestimulus activity were disregarded. Note that when x 1 / 2 represent the prestimulus and poststimulus interval from the same population and brain state, this largest t-cluster is obtained by splitting in half the prestimulus activity x 2 .
Then, we built a null distribution for each of the remaining ordered t-cluster statistics as follows. We shuffled the simulations of x ¯ 2 with respect to signal x ¯ 1 and extracted again the t-cluster statistics, sorting them in descending order. We shuffled the simulations until we reached 1000 permuted clusters for the smallest t-cluster statistic. Each null distribution was assigned to the observed t-cluster statistic with the same ordinal position. Finally, for each observed t-cluster statistic, we computed its significance from the proportion of clusters that were larger than the observed one in the corresponding null distribution.

2.8.4. Quantification of Information Propagation

We quantified the propagation of the external input from the perturbed to the unperturbed population as a function of β and g ¯ AMPA k . When x ¯ 1 , 2 corresponds to the firing rate signals of the unperturbed population during the poststimulus/prestimulus interval, the magnitude and length of the first t-cluster statistic after stimulus onset were used to quantify the amount of information propagation. A bigger and longer t-cluster statistic was indicative of a higher transmission of information to the unperturbed population. The magnitude (in s) and length (in ms) of the first t-cluster statistic were plotted as a function of β and g ¯ AMPA k .

3. Results

In the following, we first present the results for the one-cortical-column model. We show that the model reproduces the dynamical features of NREM sleep and wakefulness, given the parameters shown in Table 2. Then, we show that in response to a transient increase in the input firing rate, the model reproduces the experimental results obtained by Nir and colleagues [7]. In particular, the amplitude of the response is larger and followed by a deeper negative deflection, corresponding to a poststimulus suppression of neuronal activity, in NREM sleep as compared to wakefulness.
Second, we extend the results to two coupled cortical columns and we show that the dynamics of each brain state are preserved in this scenario. Then, we study the propagation of the external input from the perturbed cortical column to the unperturbed cortical column. Finally, we show that the upscaling of the average AMPAergic conductance during wakefulness is not sufficient for propagating the response to the external input. Instead, a higher increase in the AMPAergic conductance of the inter-cortical excitatory connectivity with respect to the intra-cortical connectivity is necessary to transmit information.
All the figures shown in the main text refer to the pyramidal population. The corresponding figures for the inhibitory population are provided in the supplementary material (Figures S1, S2, S4–S7 and S9).

3.1. One-Cortical-Column Model

3.1.1. Spontaneous Activity across Brain States

Spontaneous activity differs between NREM sleep and wakefulness. In the model, it was obtained by potentiating g ¯ AMPA (see Section 2.4). All signal types—LFP, firing rate and membrane potential—showed larger fluctuations during NREM sleep than during wakefulness (see one example of a simulated signal for LFP, firing rate and V p in NREM sleep and wakefulness in Figure 2A). The distribution of the LFP signals during NREM sleep is bimodal, whereas for wakefulness is unimodal (see Figure 2B) in agreement with the occurrence of Up and Down states in the former case [3,4,5,6,7,8].
The firing rate and V p signals during NREM sleep were split into Up and Down states (see Section 2.8.2). The distribution of the firing rate during Down states peaks at zero, confirming that it is a quasi-silent activity mode of the NREM dynamics. On the contrary, the firing rate during Up states peaks at a higher value, corresponding to persistent activity (see the top row in Figure 2C). This bimodality was also present in the V p signal, where its mean value during Up states had a more depolarized value than the Down state (see the bottom row in Figure 2C). Note that the distribution of firing rates and V p during wakefulness is closer to the Up state distribution than to the Down state distribution.
LFP signals had more spectral power within the delta (0.5–2 Hz) and SWA frequency band (below 4 Hz) during NREM sleep than during wakefulness; whereas gamma power (>30 Hz) was highest in wakefulness (see left panel in Figure 2D). The distribution of high-/low-frequency (where high is above 30 Hz and low is below 4 Hz) power ratio in the LFP confirmed the spectral separation between NREM sleep and wakefulness (see right panel in Figure 2D). These results were confirmed in the inhibitory population as well (see Figure S1). Therefore, the simulated spontaneous activity in the NREM sleep and wakefulness parameter space exhibited the experimental dynamical features characteristic of each brain state.

3.1.2. Evoked Responses in the One-Cortical-Column Model

An external input consisting of an increase of 1 ms 1 in the mean of the noise term ϕ during 100 ms is applied to the pyramidal population. The evoked responses in poststimulus intervals were significantly different from the spontaneous activity in prestimulus intervals in all signal types during NREM sleep and wakefulness (see Figure 3A). Input responses showed significant deviations from spontaneous activity in two time intervals. The first period corresponded to a positive deflection of the signals above prestimulus values (corresponding to a poststimulus neuronal activity in response to the external input), while the second interval corresponded to a negative deflection well beyond stimulus offset (corresponding to a poststimulus suppression of neuronal activity) (see Figure 3A).
We also compared the evoked responses between brain states (see Figure 3B). Since the spontaneous value in wakefulness is higher than in NREM sleep, we removed its temporal average to compare the amplitude of the positive and negative deflections of the firing rate and V p signals. We observed that those deflections were significantly larger in NREM sleep than during wakefulness. Therefore, both the response to the external input and later suppression of neuronal activity were higher in NREM sleep than in wakefulness. However, the positive and negative deflections in the evoked response captured in the LFP were higher in wakefulness than in NREM sleep. These results regarding the firing rate and the negative deflections of the LFP are in agreement with the reported experimental behavior in [7]. These results were confirmed in the inhibitory population as well (see Figure S2).

3.2. Two-Cortical-Column Model

The dynamics of each brain state were preserved after coupling two cortical columns (see Figure S3) although the excitatory inter-connectivity caused a decrease in the occurrence of Down states and a depolarization of V p during the Up states of NREM sleep compared to the case of a single cortical model (see Figure S3B,C).
LFP signals in NREM sleep carried more spectral power within the delta and SWA frequency band. Gamma power was highest in wakefulness. Further analysis of high-/low-frequency power ratio in the LFP signals confirmed the spectral separation between NREM sleep and wakefulness in the two cortical model scenario (see Figure S3D). Thus, the coupling between the two cortical columns did not change qualitatively the intrinsic dynamics of each brain state. Results were consistent in the inhibitory population as well (see Figure S4).

3.2.1. Evoked Responses in the Two-Cortical-Column Model

The evoked responses of the perturbed population in the perturbed cortical column in the two-cortical-column scenario were maintained, i.e., they were significantly different from the spontaneous activity during the prestimulus interval in all signal types during NREM sleep and wakefulness (see Figure 4A). Similarly to Figure 3A, two significant temporal clusters appeared following stimulus onset.
On the contrary, the evoked responses in the unperturbed population captured in the firing rate and V p signals were not significantly different from the spontaneous activity neither during NREM sleep nor during wakefulness (Figure 4B). In the latter case, that was so despite the increase in the average AMPAergic conductance, g ¯ AMPA , affecting the excitatory inter- and intra-connections. This upscaling of the conductance of excitatory synapses that allowed us to move the system from NREM to wakefulness dynamics, was not able to produce a significant LFP response, neither in the firing rate nor in the membrane potential, of the unperturbed population in response to the perturbed population (see right column in Figure 4B). These results were confirmed in the inhibitory populations as well (see Figure S5).

3.2.2. Information Propagation as a Function of β and g ¯ AMPA

Figure 4B showed that the response of the perturbed population to the external input was not propagated to the unperturbed population in neither NREM sleep nor wakefulness. During NREM sleep, we aim at reproducing such local responses in agreement with the experiments [7,9]. However, this is at odds with a wider spread of electrical activity during wakefulness. Even increasing g ¯ AMPA from 2 to 6 for all inter- and intra-excitatory synapses (corresponding to an inter-/intra-conductance ratio of cortical excitatory synapses β = 1 ) was not able to propagate the input to the unperturbed column. Interestingly, a selective increase of only the inter-excitatory synapses ( β = 2 ) showed a significant positive deflection at the firing rate signal in the unperturbed cortical column. This corresponds to an increase of neural activity after the onset of the external input with respect to the spontaneous activity for both g ¯ AMPA = 2 and g ¯ AMPA = 6 (see Figure 5). These results were confirmed in the inhibitory populations as well (see Figure S6). The size and duration of such a positive deflection improved as a function of β in all signal types (see Figure 6B,C for the firing rate signal).
Our results show that keeping β = 1 does not trigger a significant response in the unperturbed population even when potentiating g ¯ AMPA to 10, at which value no significant t-cluster statistic still appears (see black squares at β = 1 in Figure 6B,C). This non-selective increase of g ¯ AMPA only gives rise to a larger amplitude fluctuation in the LFP but not in the firing rate nor in V p (results not shown here). Importantly, a selective increase of g ¯ AMPA for the inter-excitatory connections over the intra-excitatory connections, quantified by β , is required to enlarge the size and duration of the positive deflection. Results were consistent in the inhibitory population as well (see Figure S7). This dependence of information propagation on β still held when performing a Bonferroni corrected t-test rather than a temporal clustering nonparametric test (see Figures S8 and S9 for pyramidal and inhibitory population, respectively).
The observed dependence of information propagation on β is due to the fact that the input–output function of any system depends on the signal to noise ratio (SNR). In the case of neural populations, the local synaptic currents would play the role of a stochastic signal and the currents activated by the external input would form the signal. Therefore, the non-selective increase of synaptic excitatory afferents within and between cortical columns from NREM sleep to wakefulness would not change the SNR. In order to boost signal transmission, synaptic potentiation needs to be larger for the synapses between cortical columns than for the synapses within cortical columns.

4. Discussion

There is compelling experimental evidence that auditory stimuli elicit firing responses in the primary auditory cortex that are similar in NREM sleep and wakefulness [7,9]. This poses a challenge to the understanding of auditory perception across the wake–sleep cycle because sounds seem to bypass the thalamus in both brain states. However, the observation that neuronal responses to either auditory stimuli [7,8,9,11,12] and non-sensory stimuli [15,16] are elicited in areas located at distances from the stimulation site that are larger in wakefulness than in NREM sleep opens the door to a new hypothesis. In particular, researchers have suggested that it is the propagation of sounds throughout the cortex that is compromised during sleep, given that a neuronal response is detected in primary brain areas at the initial stages of auditory processing. This hypothesis is supported by the fact that effective connectivity during sleep is reduced with respect to wakefulness, a feature that has also been explored in large-scale models [13,39].
Here, we have used a computational approach to understand how input propagation is disrupted, beyond primary sensory encoding, by tuning the conductance of synaptic connections. By means of a neuronal mass model operating in two distinct brain dynamics, NREM sleep and wakefulness, we investigate the response of a neural system composed of two cortical columns to an external perturbation. Our simulated signals—LFP, firing rate and membrane potential—have a bimodal distribution in the parametric space belonging to NREM dynamics, which is indicative of the alternation between Up and Down states and the presence of slow oscillations. The higher power at low frequencies in NREM versus the higher power at higher frequencies in wakefulness of the LFP also validates the existence of these two different dynamical regimes.
Few models exhibiting both regimes have been published up to date [26,27,28,29,54]. Despite the fact that any computational model with multiple parameters exhibits a rich variety of bifurcations, we show that the sleep–waking transition can be obtained with the minimum set of possible changes. Naturally, the simplest solution is to vary a single parameter, in contrast with the proposal of [27,31]. Here, an increase in the AMPAergic conductance g ¯ AMPA k suffices to move the system from NREM-like to waking-like dynamics. This is so because in the model, the membrane potential relaxes back to its steady state value producing damped oscillations in response to a perturbation for the parameter space of the NREM state (see Figure S10) due to the presence of a stable focus [55]. This behavior generates the slow oscillation in the presence of noise. Moreover, this change is consistent with the waxing and waning of synaptic strengths across the wake–sleep transition that SHY proposes [17,18,19,20,21,22,23,24,25].
In this study, the connectivity between the two cortical columns is excitatory and symmetric, but only one of the pyramidal populations receives a transient external input. We investigate how this perturbation is transmitted to the other column in such a way that, in agreement with the aforementioned experiments, its response is enhanced in the awake state when compared to NREM. Previous computational work [30] also aimed at reproducing the distinct propagation patterns that arise across the SWC, obtaining similar results to the ones described here. In particular, a wider propagation of the external perturbation occurs during wakefulness than during NREM sleep. In our case, this is quantified by the response of the unperturbed population, whereas [30] uses the perturbational complexity index (PCI) [56]. However, the set of parameters’ changes used by [30] focus on an adaptation current that weakens during wakefulness, which in our opinion is less experimentally justified. Moreover, the authors showed that during wakefulness the ratio between poststimulus and prestimulus firing rate can be greater in an area that is not directly stimulated than in the stimulation site, which requires further experimental verification.
We show that when the conductance of the excitatory synapses, g ¯ AMPA k , within a cortical column and between neighboring cortical columns are equal, there is no propagation of the input to the unperturbed cortical column. Propagation only occurs when the synaptic conductance between cortical columns is selectively boosted with respect to the synaptic conductance within cortical columns, a ratio that we have quantified with the parameter β . It is important to note that this difference between the inter- and intra-AMPAergic conductance does not directly result in a higher AMPAergic external than local synaptic afferent current. In fact, in the model, local synaptic currents surpass external ones both during NREM sleep and wakefulness. However, during wakefulness, this synaptic imbalance increases. Therefore, our results suggest that the synaptic homeostasis hypothesis should be heterogeneously applied, that is, that synaptic potentiation has to be spatially organized and distant synaptic conductances have to be scaled with respect to the local ones. Although this result needs experimental validation, it has been reported that, in the cortex, perforated synapses enlarge their axon–spine interface after waking relative to sleep, in contrast to the lack of change found in non-perforated synapses [25]. This morphological change specific to perforated synapses also speaks for their major role in synaptic plasticity and the selective implementation of synaptic homeostasis that we propose.
Despite the fact that our model cannot account for whole brain interactions and is restricted to the dialogue between two networks, our results are derived from physiological assumptions and reproduce verified experimental observations. A more complex connectivity diagram comprising other cortical and subcortical structures could, nonetheless, provide spatial information about stimulus propagation. Subcortical structures, such as the brainstem and basal forebrain, are the substrate for the changes in the concentration of neuromodulators that occur throughout the sleep–weak cycle [57]. These neuromodulators target vesicular release at the presynaptic site or transmitter receptors at the postsynaptic site and alter synaptic strength [58,59]. We have phenomenologically implemented their role by tuning synaptic conductances but the precise effect of these substances (serotonin, acetylcholine, histamine, etc.) on target brain areas is out of the scope of this study, partly because the exact microcircuitry of the ascending arousal system is still unknown.
It is important to note that we have excluded REM sleep from our work, which is a short but relevant stage of sleep that follows NREM epochs. The LFP and EEG dynamics recorded during REM sleep show similar features to wakefulness. Since firing responses in the primary auditory cortex are preserved across the SWC, including REM sleep [7,8], computationally, we could use the same set of parameters to model REM than wakefulness at a single cortical column level. However, tracking these responses up to the perirhinal cortex, shows that their amplitude is attenuated from wakefulness to NREM sleep, and has intermediate values in REM sleep [8]. In our model, this could be replicated by using a β value for REM sleep smaller than for wakefulness (see Figure 3B,C). However, experimentally, SHY has only proven to be true when comparing wakefulness and NREM sleep [25]. Therefore, we need to be cautious before extending our conclusions to REM sleep.
Despite the limitations of our model, our work opens the door to explore which rules need to be implemented to reproduce experimental data and guide computational principles. This is particularly important while we still lack a clear understanding of how inputs scattered along the wide dendritic tree of a pyramidal neuron are added to combine synaptic sources at short and long distances with respect to the soma.

5. Conclusions

Sleep and wakefulness are distinguished by distinct behavioural, electrophysiological and molecular features. It is still unclear how these biological layers are causally connected. Here, we have explored one particular relationship between two experimentally observed changes that occur across the SWC: the potentiation of excitatory synapses and the broad distribution of incoming information throughout the cortex that takes place during wakefulness. The former happens during learning and strengthens synaptic connections, a process that is compensated during NREM sleep to enable synaptic homeostasis. The later evidence may speak for a cortical gating mechanism by which signals coming from the external world only reach higher processing brain areas during wakefulness and could possibly explain the emergence of consciousness. Our computational model links both phenomena by showing that the potentiation of excitatory conductances not only triggers a dynamic change of the electrical activity of the neuronal networks, but also of the propagation pattern across them. However, we predict that, in order for a spatially wider response to occur during wakefulness, it is necessary that such potentiation occurs preferentially between distinct networks over local and recurrent connections.

Supplementary Materials

Author Contributions

Conceptualization, F.R. and B.S.; methodology, F.R. and B.S.; software, F.R.; validation, F.R.; formal analysis, F.R.; investigation, F.R. and B.S.; resources, F.R. and B.S.; data curation, F.R.; writing—original draft preparation, F.R. and B.S.; writing—review and editing, F.R., B.S. and R.M.-B.; visualization, F.R.; supervision, B.S.; project administration, B.S.; funding acquisition, B.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Postdoctoral Junior Leader Fellowship Programme from La Caixa Banking Foundation (LCF/BQ/PI18/11630004), the Howard Hughes Medical Institute (HHMI ID 55008742) and the Human Brain Project (HBP SGA3 Grant Agreement No. 945539).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

We do not report any data. The simulation routine of the model is available at github [48].

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
NREM sleepNon-rapid eye movement sleep
SWASlow wave activity
TMSTranscranial magnetic stimulation
SHYSynaptic homeostasis hypothesis
LFPLocal field potential
SWCSleep-waking cycle
STFTShort-Time Fourier Transform
PSDPower spectral density
SNRSignal to noise ratio

Appendix A. Neural Mass Model Equations

In this section, we provide the full model equations for the one- and two-cortical-column models. The cortical column model is described by the neural mass model of pyramidal and inhibitory populations mutually connected through AMPAergic and GABAergic synapses.

Appendix A.1. Synapse Dynamics

Assuming that the postsynaptic response of a synapse s k k —given a presynaptic population k and a postsynaptic population k—is obtained by the convolution of the input at time t at the synaptic site H ( t ) with the average synaptic response to a single spike α k , s k k ( t ) takes the following form:
s k k ( t ) = H ( t ) α k ( t ) , = 0 t H ( τ ) α k ( t τ ) d τ ,
where α k ( t ) has an exponential decay time course:
α k ( t ) = γ k 2 t exp ( γ k t ) ,
Therefore, the first derivative of the s k k ( t ) with respect to t taking into account the Leibniz integral rule is as follows:
s ˙ k k ( t ) = d d t 0 t H ( τ ) α k ( t τ ) d τ , = H ( t ) α k ( t t ) + 0 t H ( τ ) d d t α k ( t τ ) d τ , = 0 t H ( τ ) γ k 2 exp ( γ k ( t τ ) ) γ k 3 ( t τ ) exp ( γ k ( t τ ) ) d τ , = 0 t H ( τ ) γ k 2 exp ( γ k ( t τ ) ) d τ γ k 0 t H ( τ ) γ k 2 ( t τ ) exp ( γ k ( t τ ) ) d τ , = 0 t H ( τ ) γ k 2 exp ( γ k ( t τ ) ) d τ γ k 0 t H ( τ ) α k ( t τ ) d τ , = 0 t H ( τ ) γ k 2 exp ( γ k ( t τ ) ) d τ γ k s k k ( t ) .
In the same way, the second derivative of the s k k ( t ) with respect to t is:
s ¨ k k ( t ) = γ k 2 H ( t ) + γ k 0 t H ( τ ) γ k 2 ( t τ ) exp ( γ k ( t τ ) ) d τ γ k s ˙ k k ( t ) , = γ k 2 H ( t ) γ k ( s ˙ k k ( t ) + γ k s k k ( t ) ) γ k s ˙ k k ( t ) , = γ k 2 H ( t ) γ k 2 s k k ( t ) 2 γ k s ˙ k k ( t ) , = γ k 2 H ( t ) s k k ( t ) 2 γ k s ˙ k k ( t ) ,
Finally, by substituting the H ( t ) = N k k Q k ( V k ) + ϕ —which is the presynaptic drive at synapse of population k—the postsynaptic response s k k takes the form of:
s ¨ k k = γ k 2 ( N k k Q k ( V k ) + ϕ s k k ) 2 γ k , s ˙ k k
Please note that the noise term ϕ is zero for the inhibitory synapses and is a stochastic Gaussian process with zero autocorrelation time constant for the excitatory synapses.

Appendix A.2. One-Cortical-Column Model

The one-cortical-column model consists of one pyramidal and one inhibitory population (see panel A in Figure 1). The model utilizes the mathematical formalism of the Hodgkin–Huxley model to describe the membrane voltage activity in terms of driving synaptic activity. The synaptic activity contains one leak, two synaptic (AMPAergic and GABAergic) currents and one activity-dependent potassium current. In the absence of extrinsic input (the Gaussian process), the system settles down on the steady state solution. At every time point, the extrinsic input pushes the system out of the steady state solution, however, the aforementioned synaptic currents drive the system back to the steady state solution. The full model equations for the one-cortical-column model are described by:
τ p V ˙ p = I L p I AMPA p I GABA p τ p C m 1 I KNa , τ i V i ˙ = I L i I AMPA i I GABA i , s ¨ p p = γ p 2 ( N p p Q p ( V p ) + ϕ s p p ) 2 γ p s ˙ p p , s ¨ p i = γ i 2 ( N p i Q i ( V i ) s p i ) 2 γ i s ˙ p i , s ¨ i p = γ p 2 ( N i p Q p ( V p ) + ϕ s i p ) 2 γ p s ˙ i p , s ¨ i i = γ i 2 ( N i i Q i ( V i ) s i i ) 2 γ i s ˙ i i , τ Na [ Na ] ˙ = α Na Q p ( V p ) Na pump ( [ Na ] ) ,
with the currents defined by:
I L k = g ¯ L ( V k E L k ) , I AMPA k = g ¯ AMPA k s k p ( V k E AMPA ) , I GABA k = g ¯ GABA k s k i ( V k E GABA ) , I KNa = g ¯ KNa 0.37 1 + ( 38.7 [ Na ] ) 3.5 ( V p E K ) ,
the sodium pump and the cortical firing rate functions are given by:
Na pump ( [ Na ] ) = R pump [ Na ] 3 [ Na ] 3 + 3375 [ Na ] eq 3 [ Na ] eq 3 + 3375 , Q k ( V k ) = Q k m a x ( 1 + tanh ( C ( V k θ k ) / σ k ) ) 2 , C = π 2 3
Symbol descriptions and the parameter values are provided in Table 1 and Table 2, respectively.

Appendix A.3. Two-Cortical-Column Model

The two-cortical-column model consists of two cortical columns each containing one pyramidal and one inhibitory population (see Appendix A.2) that are mutually coupled through AMPAergic connections (see panel B in Figure 1). The full model equations for the two-cortical-column model are described by:
τ p V ˙ p = I L p I AMPA p I GABA p τ p C m 1 I KNa p , τ i V i ˙ = I L i I AMPA i I GABA i , τ p V ˙ p = I L p I AMPA p I GABA p τ p C m 1 I KNa p , τ i V i ˙ = I L i I AMPA i I GABA i , s ¨ p p = γ p 2 ( N p p Q p ( V p ) + ϕ s p p ) 2 γ p s ˙ p p , s ¨ p i = γ i 2 ( N p i Q i ( V i ) s p i ) 2 γ i s ˙ p i , s ¨ i p = γ p 2 ( N i p Q p ( V p ) + ϕ s i p ) 2 γ p s ˙ i p , s ¨ i i = γ i 2 ( N i i Q i ( V i ) s i i ) 2 γ i s ˙ i i , s ¨ p p = γ p 2 ( N p p Q p ( V p ) + ϕ s p p ) 2 γ p s ˙ p p , s ¨ p i = γ i 2 ( N p i Q i ( V i ) s p i ) 2 γ i s ˙ p i , s ¨ i p = γ p 2 ( N i p Q p ( V p ) + ϕ s i p ) 2 γ p s ˙ i p , s ¨ i i = γ i 2 ( N i i Q i ( V i ) s i i ) 2 γ i s ˙ i i , s ¨ p p inter = γ p 2 ( N p p Q p ( V p ) s p p inter ) 2 γ p s ˙ p p inter , s ¨ i p inter = γ p 2 ( N i p Q p ( V p ) s i p inter ) 2 γ p s ˙ i p inter , s ¨ p p inter = γ p 2 ( N p p Q p ( V p ) s p p inter ) 2 γ p s ˙ p p inter , s ¨ i p inter = γ p 2 ( N i p Q p ( V p ) s i p inter ) 2 γ p s ˙ i p inter , τ Na [ Na ] ˙ p = α Na Q p ( V p ) Na pump ( [ Na ] p ) , τ Na [ Na ] ˙ p = α Na Q p ( V p ) Na pump ( [ Na ] p ) ,
where I KNa h is the sodium-dependent potassium current of either the perturbed pyramidal population ( h = p ) or the unperturbed pyramidal population ( h = p ). The currents are defined by:
I L k = g ¯ L ( V k E L k ) , I AMPA k = g ¯ AMPA k s k ( p / p ) ( V k E AMPA ) + β · g ¯ AMPA k s k ( p / p ) inter ( V k E AMPA ) , I GABA k = g ¯ GABA k s k ( i / i ) ( V k E GABA ) , I KNa h = g ¯ KNa 0.37 1 + ( 38.7 [ Na ] h ) 3.5 ( V h E K ) ,
the sodium pump and the cortical firing rate functions are given by:
Na pump ( [ Na ] h ) = R pump [ Na ] h 3 [ Na ] h 3 + 3375 [ Na ] eq 3 [ Na ] eq 3 + 3375 , Q k ( V k ) = Q k m a x ( 1 + tanh ( C ( V k θ k ) / σ k ) ) 2 , C = π 2 3
Symbol descriptions and the parameter values are provided in Table 1 and Table 3, respectively.

Appendix B. Dynamical Constrains on g ¯ GABA k

Appendix B.1. One-Cortical-Column Model

Upscaling the conductance of excitatory synapses by increasing the average AMPAergic conductance, g ¯ AMPA k , which is in agreement with SHY allows us to place the model from NREM sleep dynamics into wakefulness. The average GABAergic conductance on pyramidal and inhibitory population are increased to counterbalance the more excitation each of the populations receives due to upscaling the average AMPAergic conductance. To do so, we use the peaks of V p and V i values during Up states (see Section 2.8.2 for Up and Down definition) as the steady state value of the average membrane potential of pyramidal and inhibitory populations during wakefulness. The average GABAergic conductance on pyramidal and inhibitory populations are increased until these V p and V i values are obtained for the g ¯ AMPA k = 2 responsible of wakefulness.
The steady state solution of the model equations for the one-cortical-column (see Appendix A.2) are obtained by setting derivative of all variables to zero.
0 = I L p I AMPA p I GABA p τ p C m 1 I KNa p , 0 = I L i I AMPA i I GABA i , s k k = N k k Q k ( V k ) , [ Na ] = A · 3375 1 A 3 , A = α Na R pump Q p ( V p ) + [ Na ] eq 3 [ Na ] eq 3 + 3375 ,
therefore, the average GABAergic conductance on pyramidal and inhibitory populations, which is required to keep the average membrane potential of pyramidal and inhibitory populations during wakefulness equal to the peaks of V p and V i during Up states, is as follows:
g ¯ GABA p = g ¯ L ( V p E L p ) g ¯ AMPA p s p p ( V p E AMPA ) g ¯ KNa 0.37 1 + ( 38.7 [ Na ] ) 3.5 ( V p E K ) s p i ( V p E GABA ) , g ¯ GABA i = g ¯ L ( V p E L i ) g ¯ AMPA i s i p ( V p E AMPA ) s i i ( V i E GABA ) ,
The average GABAergic conductance of the pyramidal population and that of the inhibitory population during wakefulness are obtained by setting g ¯ AMPA p / i = 2 and the peaks of V p and V i values during Up states in the above equations.

Appendix B.2. Two-Cortical-Column Model

Inter-cortical column coupling introduces more excitation to the pyramidal and inhibitory populations. To keep the steady state value of V p and V i in the two-cortical-column scenario equal to the ones in the one-cortical-column model both during NREM sleep and wakefulness, the average GABAergic conductances on pyramidal and inhibitory populations are increased to counterbalance this inter-cortical column coupling.
The steady state solution of the model equations for the two-cortical-column (see Appendix A.3) is obtained by setting the derivatives of all variables to zero:
0 = I L p / p I AMPA p / p I GABA p / p τ p / p C m 1 I KNa p / p , 0 = I L i / i I AMPA i / i I GABA i / i , s k k = N k k Q k ( V k ) , s k k inter = N k k Q k ( V k ) [ Na ] p / p = A p / p · 3375 1 A p / p 3 , A p / p = α Na R pump Q p / p ( V p / p ) + [ Na ] eq 3 [ Na ] eq 3 + 3375 ,
by taking into account the symmetry between the two cortical columns, we set the steady state value of V p and V i equal to V p and V i ( V p = V p and V i = V i ), respectively. Therefore, to keep the average membrane potential of pyramidal and inhibitory populations in the two-cortical-column scenario equal to the ones in the one-cortical-column model both during NREM sleep and wakefulness, the average GABAergic conductances on pyramidal and inhibitory populations change as follows:
g ¯ GABA p = g ¯ L ( V p E L p ) g ¯ AMPA p s p p ( V p E AMPA ) β g ¯ AMPA p s p p ( V p E AMPA ) g ¯ KNa 0.37 1 + ( 38.7 [ Na ] ) 3.5 ( V p E K ) s p i ( V p E GABA ) , g ¯ GABA i = g ¯ L ( V p E L i ) g ¯ AMPA i s i p ( V p E AMPA ) β g ¯ AMPA i s i p ( V p E AMPA ) s i i ( V i E GABA ) ,
The average GABAergic conductances on pyramidal and inhibitory populations during NREM sleep in the two-cortical-column scenario are obtained by setting g ¯ AMPA p / i = 1 and the steady state values of V p and V i during NREM sleep in the one-cortical-column scenario. Therefore, the one cortical column in NREM sleep regime is run for 15 seconds using a stochastic Heun method [47] with a step size of 0.1 ms in the absence of noise to obtain the steady state values of V p and V i .
The average GABAergic conductances on pyramidal and inhibitory populations during wakefulness in the two-cortical-column scenario are obtained by setting g ¯ AMPA p / i = 2 , β = 1 and the peaks of V p and V i values during Up states in the one-cortical-column scenario in the aforementioned equations.
Moreover, the average GABAergic conductances on pyramidal and inhibitory populations during wakefulness in the two-cortical-column scenario for variations of β and g ¯ AMPA k are obtained by setting the corresponding β , g ¯ AMPA p / i and the peaks of V p and V i values during Up states in the one-cortical-column scenario in the aforementioned equations.

References

  1. Brown, R.E.; Basheer, R.; McKenna, J.T.; Strecker, R.E.; McCarley, R.W. Control of Sleep and Wakefulness. Physiol. Rev. 2012, 92, 1087–1187. [Google Scholar] [CrossRef] [Green Version]
  2. Borbély, A.A. A two-process model of sleep regulation. Hum. Neurobiol. 1982, 1, 195–204. [Google Scholar] [PubMed]
  3. Sanchez-Vives, M.V.; McCormick, D.A. Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat. Neurosci. 2000, 3, 1027–1034. [Google Scholar] [CrossRef] [PubMed]
  4. Steriade, M.; Timofeev, I.; Grenier, F. Natural waking and sleep states: A view from inside neocortical neurons. J. Neurophysiol. 2001, 85, 1969–1985. [Google Scholar] [CrossRef]
  5. Vyazovskiy, V.V.; Olcese, U.; Lazimy, Y.M.; Faraguna, U.; Esser, S.K.; Williams, J.C.; Cirelli, C.; Tononi, G. Cortical Firing and Sleep Homeostasis. Neuron 2009, 63, 865–878. [Google Scholar] [CrossRef] [PubMed]
  6. Vyazovskiy, V.V.; Cui, N.; Rodriguez, A.V.; Funk, C.; Cirelli, C.; Tononi, G. The Dynamics of Cortical Neuronal Activity in the First Minutes after Spontaneous Awakening in Rats and Mice. Sleep 2014, 37, 1337–1347. [Google Scholar] [CrossRef] [Green Version]
  7. Nir, Y.; Vyazovskiy, V.V.; Cirelli, C.; Banks, M.I.; Tononi, G. Auditory Responses and Stimulus-Specific Adaptation in Rat Auditory Cortex are Preserved Across NREM and REM Sleep. Cereb. Cortex 2015, 25, 1362–1378. [Google Scholar] [CrossRef] [Green Version]
  8. Sela, Y.; Krom, A.J.; Bergman, L.; Regev, N.; Nir, Y. Sleep differentially affects early and late neuronal responses to sounds in auditory and perirhinal cortices. J. Neurosci. 2020, 40, 2895–2905. [Google Scholar] [CrossRef]
  9. Issa, E.B.; Wang, X. Sensory Responses during Sleep in Primate Primary and Secondary Auditory Cortex. J. Neurosci. 2008, 28, 14467–14480. [Google Scholar] [CrossRef]
  10. Bastuji, H.; Lamouroux, P.; Villalba, M.; Magnin, M.; Garcia-Larrea, L. Local sleep spindles in the human thalamus. J. Physiol. 2020, 598, 2109–2124. [Google Scholar] [CrossRef]
  11. Portas, C.M.; Krakow, K.; Allen, P.; Josephs, O.; Armony, J.L.; Frith, C.D. Auditory Processing across the Sleep-Wake Cycle. Neuron 2000, 28, 991–999. [Google Scholar] [CrossRef] [Green Version]
  12. Krom, A.J.; Marmelshtein, A.; Gelbard-Sagiv, H.; Tankus, A.; Hayat, H.; Hayat, D.; Matot, I.; Strauss, I.; Fahoum, F.; Soehle, M.; et al. Anesthesia-induced loss of consciousness disrupts auditory responses beyond primary cortex. Proc. Natl. Acad. Sci. USA 2020, 117, 11770–11780. [Google Scholar] [CrossRef]
  13. Jobst, B.M.; Hindriks, R.; Laufs, H.; Tagliazucchi, E.; Hahn, G.; Ponce-Alvarez, A.; Stevner, A.B.A.; Kringelbach, M.L.; Deco, G. Increased Stability and Breakdown of Brain Effective Connectivity During Slow-Wave Sleep: Mechanistic Insights from Whole-Brain Computational Modelling. Sci. Rep. 2017, 7, 4634. [Google Scholar] [CrossRef] [PubMed]
  14. Parrino, L.; Ferri, R.; Bruni, O.; Terzano, M.G. Cyclic alternating pattern (CAP): The marker of sleep instability. Sleep Med. Rev. 2012, 16, 27–45. [Google Scholar] [CrossRef] [PubMed]
  15. Massimini, M. Breakdown of Cortical Effective Connectivity During Sleep. Science 2005, 309, 2228–2232. [Google Scholar] [CrossRef] [PubMed]
  16. Massimini, M.; Ferrarelli, F.; Murphy, M.J.; Huber, R.; Riedner, B.A.; Casarotto, S.; Tononi, G. Cortical reactivity and effective connectivity during REM sleep in humans. Cogn. Neurosci. 2010, 1, 176–183. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Tononi, G.; Cirelli, C. Sleep and synaptic homeostasis: A hypothesis. Brain Res. Bull. 2003, 62, 143–150. [Google Scholar] [CrossRef]
  18. Tononi, G.; Cirelli, C. Sleep function and synaptic homeostasis. Sleep Med. Rev. 2006, 10, 49–62. [Google Scholar] [CrossRef]
  19. Tononi, G.; Cirelli, C. Sleep and the Price of Plasticity: From Synaptic and Cellular Homeostasis to Memory Consolidation and Integration. Neuron 2014, 81, 12–34. [Google Scholar] [CrossRef] [Green Version]
  20. Tononi, G.; Cirelli, C. Sleep and synaptic down-selection. Eur. J. Neurosci. 2020, 51, 413–421. [Google Scholar] [CrossRef] [Green Version]
  21. Cirelli, C.; Tononi, G. The why and how of sleep-dependent synaptic down-selection. Semin. Cell Dev. Biol. 2021, in press. [Google Scholar] [CrossRef] [PubMed]
  22. Cirelli, C.; Tononi, G. Differential expression of plasticity-related genes in waking and sleep and their regulation by the noradrenergic system. J. Neurosci. 2000, 20, 9187–9194. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Vyazovskiy, V.V.; Cirelli, C.; Pfister-Genskow, M.; Faraguna, U.; Tononi, G. Molecular and electrophysiological evidence for net synaptic potentiation in wake and depression in sleep. Nat. Neurosci. 2008, 11, 200–208. [Google Scholar] [CrossRef] [PubMed]
  24. Liu, Z.W.; Faraguna, U.; Cirelli, C.; Tononi, G.; Gao, X.B. Direct evidence for wake-related increases and sleep-related decreases in synaptic strength in rodent cortex. J. Neurosci. 2010, 30, 8671–8675. [Google Scholar] [CrossRef] [Green Version]
  25. Cirelli, C.; Tononi, G. Effects of sleep and waking on the synaptic ultrastructure. Philos. Trans. R. Soc. B Biol. Sci. 2020, 375. [Google Scholar] [CrossRef] [PubMed]
  26. Hill, S.; Tononi, G. Modeling sleep and wakefulness in the thalamocortical system. J. Neurophysiol. 2005, 93, 1671–1698. [Google Scholar] [CrossRef]
  27. Weigenand, A.; Schellenberger Costa, M.; Ngo, H.V.V.; Claussen, J.C.; Martinetz, T. Characterization of k-complexes and slow wave activity in a neural mass model. PLoS Comput. Biol. 2014, 10, e1003923. [Google Scholar] [CrossRef]
  28. Wei, Y.; Krishnan, G.P.; Komarov, M.; Bazhenov, M. Differential roles of sleep spindles and sleep slow oscillations in memory consolidation. PLoS Comput. Biol. 2018, 14, e1006322. [Google Scholar] [CrossRef] [Green Version]
  29. Bensaid, S.; Modolo, J.; Merlet, I.; Wendling, F.; Benquet, P. COALIA: A Computational Model of Human EEG for Consciousness Research. Front. Syst. Neurosci. 2019, 13, 1–18. [Google Scholar] [CrossRef]
  30. Goldman, J.S.; Kusch, L.; Yalcinkaya, B.H.; Depannemaecker, D.; Nghiem, T.A.E.; Jirsa, V.; Destexhe, A. Brain-scale emergence of slow-wave synchrony and highly responsive asynchronous states based on biologically realistic population models simulated in The Virtual Brain. bioRxiv 2020. [Google Scholar] [CrossRef]
  31. Costa, M.S.; Born, J.; Claussen, J.C.; Martinetz, T. Modeling the effect of sleep regulation on a neural mass model. J. Comput. Neurosci. 2016, 41, 15–28. [Google Scholar] [CrossRef] [PubMed]
  32. Buzsáki, G.; Anastassiou, C.A.; Koch, C. The origin of extracellular fields and currents-EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 2012, 13, 407–420. [Google Scholar] [CrossRef] [PubMed]
  33. Hodgkin, A.L.; Huxley, A.F.; Katz, B. Measurement of current-voltage relations in the membrane of the giant axon of Loligo. J. Physiol. 1952, 116, 424–448. [Google Scholar] [CrossRef] [PubMed]
  34. Wilson, H.R.; Cowan, J.D. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 1972, 12, 1–24. [Google Scholar] [CrossRef] [Green Version]
  35. David, O.; Friston, K.J. A neural mass model for MEG/EEG. Neuroimage 2003, 20, 1743–1755. [Google Scholar] [CrossRef]
  36. Budd, J.M.; Kisvárday, Z.F. Local lateral connectivity of inhibitory clutch cells in layer 4 of cat visual cortex (area 17). Exp. Brain Res. 2001, 140, 245–250. [Google Scholar] [CrossRef]
  37. Buzás, P.; Kovács, K.; Ferecskó, A.S.; Budd, J.M.; Eysel, U.T.; Kisvárday, Z.F. Model-based analysis of excitatory lateral connections in the visual cortex. J. Comp. Neurol. 2006, 499, 861–881. [Google Scholar] [CrossRef] [PubMed]
  38. Antolík, J. Rapid Long-Range Disynaptic Inhibition Explains the Formation of Cortical Orientation Maps. Front. Neural Circuits 2017, 11, 1–15. [Google Scholar] [CrossRef] [Green Version]
  39. Esser, S.K.; Hill, S.; Tononi, G. Breakdown of effective connectivity during slow wave sleep: Investigating the mechanism underlying a cortical gate using large-scale modeling. J. Neurophysiol. 2009, 102, 2096–2111. [Google Scholar] [CrossRef] [Green Version]
  40. Deco, G.; Ponce-Alvarez, A.; Hagmann, P.; Romani, G.L.; Mantini, D.; Corbetta, M. How Local Excitation-Inhibition Ratio Impacts the Whole Brain Dynamics. J. Neurosci. 2014, 34, 7886–7898. [Google Scholar] [CrossRef] [Green Version]
  41. Cohen, B.P.; Chow, C.C.; Vattikuti, S. Dynamical modeling of multi-scale variability in neuronal competition. Commun. Biol. 2019, 2, 319. [Google Scholar] [CrossRef]
  42. Barbero-Castillo, A.; Mateos-Aparicio, P.; Porta, L.D.; Camassa, A.; Perez-Mendez, L.; Sanchez-Vives, M.V. Impact of gabaa and gabab inhibition on cortical dynamics and perturbational complexity during synchronous and desynchronized states. J. Neurosci. 2021, 41, 5029–5044. [Google Scholar] [CrossRef]
  43. Gerstner, W.; Kistler, W.M.; Naud, R.; Paninski, L. Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition, 1st ed.; Cambridge University Press: Cambridge, UK, 2014; pp. 58–63. [Google Scholar]
  44. Compte, A.; Sanchez-Vives, M.V.; McCormick, D.A.; Wang, X.J. Cellular and network mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical network model. J. Neurophysiol. 2003, 89, 2707–2725. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Neske, G.T. The slow oscillation in cortical and thalamic networks: Mechanisms and functions. Front. Neural Circuits 2016, 9, 1–25. [Google Scholar] [CrossRef] [Green Version]
  46. Sanchez-Vives, M.V. Origin and dynamics of cortical slow oscillations. Curr. Opin. Physiol. 2020, 15, 217–223. [Google Scholar] [CrossRef]
  47. San Miguel, M.; Toral, R. Stochastic Effects in Physical Systems. In Instab. Nonequilibrium Struct. VI; Tirapegui, E., Martínez, J., Tiemann, R., Eds.; Springer: Dordrecht, The Netherlands, 2000; pp. 35–127. [Google Scholar] [CrossRef] [Green Version]
  48. Razi, F. Simulation Routine of the Model. Available online: https://github.com/fraziphy/cortex (accessed on 20 September 2021).
  49. Mazzoni, A.; Panzeri, S.; Logothetis, N.K.; Brunel, N. Encoding of naturalistic stimuli by local field potential spectra in networks of excitatory and inhibitory neurons. PLoS Comput. Biol. 2008, 4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  50. Sancristóbal, B.; Rebollo, B.; Boada, P.; Sanchez-Vives, M.V.; Garcia-Ojalvo, J. Collective stochastic coherence in recurrent neuronal networks. Nat. Phys. 2016, 12, 881–887. [Google Scholar] [CrossRef]
  51. Freeman, W.J.; Quiroga, R.Q. Imaging Brain Function with EEG; Springer: New York, NY, USA, 2013. [Google Scholar] [CrossRef]
  52. Welch, P. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electroacoust. 1967, 15, 70–73. [Google Scholar] [CrossRef] [Green Version]
  53. Maris, E.; Oostenveld, R. Nonparametric statistical testing of EEG- and MEG-data. J. Neurosci. Methods 2007, 164, 177–190. [Google Scholar] [CrossRef]
  54. di Volo, M.; Romagnoni, A.; Capone, C.; Destexhe, A. Biologically Realistic Mean-Field Models of Conductance-Based Networks of Spiking Neurons with Adaptation. Neural Comput. 2019, 31, 653–680. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Strogatz, S.H. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, 1st ed.; Perseus Books Publishing: New York, NY, USA, 1994; pp. 150–152. [Google Scholar]
  56. Casali, A.G.; Gosseries, O.; Rosanova, M.; Boly, M.; Sarasso, S.; Casali, K.R.; Casarotto, S.; Bruno, M.A.; Laureys, S.; Tononi, G.; et al. A Theoretically Based Index of Consciousness Independent of Sensory Processing and Behavior. Sci. Transl. Med. 2013, 5, 198ra105. [Google Scholar] [CrossRef] [PubMed]
  57. Lee, S.H.; Dan, Y. Neuromodulation of Brain States. Neuron 2012, 76, 209–222. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Nadim, F.; Bucher, D. Neuromodulation of neurons and synapses. Curr. Opin. Neurobiol. 2014, 29, 48–56. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Marder, E.; Thirumalai, V. Cellular, synaptic and network effects of neuromodulation. Neural Netw. 2002, 15, 479–493. [Google Scholar] [CrossRef]
Figure 1. Diagram of the neural mass model. (A) One cortical column containing one pyramidal and one inhibitory population. Coupling between pyramidal and inhibitory populations are mediated through AMPAergic and GABAergic connections (see Table 1 and Table 2 for symbol description and parameter values). (B) Neural mass model describing two mutually coupled cortical columns. Connections between two cortical columns are only AMPAergic, i.e., only the pyramidal populations target both pyramidal and inhibitory populations from the other cortical column (see Table 1 and Table 3 for symbol description and parameter values). ξ represents the transient external input impinging only in one cortical column.
Figure 1. Diagram of the neural mass model. (A) One cortical column containing one pyramidal and one inhibitory population. Coupling between pyramidal and inhibitory populations are mediated through AMPAergic and GABAergic connections (see Table 1 and Table 2 for symbol description and parameter values). (B) Neural mass model describing two mutually coupled cortical columns. Connections between two cortical columns are only AMPAergic, i.e., only the pyramidal populations target both pyramidal and inhibitory populations from the other cortical column (see Table 1 and Table 3 for symbol description and parameter values). ξ represents the transient external input impinging only in one cortical column.
Biology 10 00945 g001
Figure 2. Spontaneous activity ( ξ = 0 ) during NREM sleep and wakefulness in the one-cortical-column model. (A) Simulated signals from one random simulation showing large amplitude fluctuations during NREM sleep between Up and Down states (left column). Lower amplitude fluctuations appear during wakefulness (right column). The first, second and third row correspond to LFP, firing rate and V p signal, respectively. (B) The distribution of the LFP signals is bimodal during NREM sleep and unimodal during wakefulness. The half distance between the two peaks of the NREM distribution is used to extract Up and Down events from the LFP. (C) Top row: The distribution of the firing rate during Up states (dark green) shows the persistent activity mode, while the Down state distribution (light green) is a silent activity mode. The distribution of firing rate during wakefulness shows a unimodal distribution. The distribution during wakefulness and Up states have similar means. Bottom row: The distribution of the V p signal during Up (dark green) and Down (light green) events differ. It has a more depolarized mean value during Up states (close to the mean value during wakefulness) than Down states. The bigger area of the distributions of Up events indicates that they are more prevalent than Down events in the modeled NREM sleep. (D) Spectral content of LFP signals during NREM sleep and wakefulness. Left: Average power spectrum density of LFP signals during NREM sleep (green) and wakefulness (black). The shaded area corresponds to the standard deviation of the mean of PSD over 500 simulations during NREM sleep and wakefulness. Right: Distribution of high-/low-frequency (where high is above 30 Hz and low is below 4 Hz) power ratio in the LFP during NREM sleep (green) and wakefulness (black).
Figure 2. Spontaneous activity ( ξ = 0 ) during NREM sleep and wakefulness in the one-cortical-column model. (A) Simulated signals from one random simulation showing large amplitude fluctuations during NREM sleep between Up and Down states (left column). Lower amplitude fluctuations appear during wakefulness (right column). The first, second and third row correspond to LFP, firing rate and V p signal, respectively. (B) The distribution of the LFP signals is bimodal during NREM sleep and unimodal during wakefulness. The half distance between the two peaks of the NREM distribution is used to extract Up and Down events from the LFP. (C) Top row: The distribution of the firing rate during Up states (dark green) shows the persistent activity mode, while the Down state distribution (light green) is a silent activity mode. The distribution of firing rate during wakefulness shows a unimodal distribution. The distribution during wakefulness and Up states have similar means. Bottom row: The distribution of the V p signal during Up (dark green) and Down (light green) events differ. It has a more depolarized mean value during Up states (close to the mean value during wakefulness) than Down states. The bigger area of the distributions of Up events indicates that they are more prevalent than Down events in the modeled NREM sleep. (D) Spectral content of LFP signals during NREM sleep and wakefulness. Left: Average power spectrum density of LFP signals during NREM sleep (green) and wakefulness (black). The shaded area corresponds to the standard deviation of the mean of PSD over 500 simulations during NREM sleep and wakefulness. Right: Distribution of high-/low-frequency (where high is above 30 Hz and low is below 4 Hz) power ratio in the LFP during NREM sleep (green) and wakefulness (black).
Biology 10 00945 g002
Figure 3. Evoked response in the one-cortical-column model. (A) Average LFP, firing rate and V p across 500 simulations for NREM sleep (left column) and wakefulness (right column). Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus versus the prestimulus intervals. In all signals, responses to the external stimulus consist of one positive and one negative significant deflection. The positive ( p < 0.0006 ) and negative ( p < 0.001 ) deflections are significant in NREM and wakefulness. (B) Average difference Δ LFP, Δ firing rate and Δ V p between NREM sleep and wakefulness across 500 simulations. Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus intervals during NREM sleep with wakefulness. The positive ( p < 0.0005 ) and negative ( p < 0.001 ) deflections in all signals are significant. The amplitude of both deflections are higher in NREM sleep than in wakefulness in the firing rate and V p . The opposite happens in the LFP signal. In all cases, the dashed vertical line, aligned at zero, indicates stimulus ξ onset. Stimulus duration is represented by a black horizontal line and light gray area.
Figure 3. Evoked response in the one-cortical-column model. (A) Average LFP, firing rate and V p across 500 simulations for NREM sleep (left column) and wakefulness (right column). Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus versus the prestimulus intervals. In all signals, responses to the external stimulus consist of one positive and one negative significant deflection. The positive ( p < 0.0006 ) and negative ( p < 0.001 ) deflections are significant in NREM and wakefulness. (B) Average difference Δ LFP, Δ firing rate and Δ V p between NREM sleep and wakefulness across 500 simulations. Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus intervals during NREM sleep with wakefulness. The positive ( p < 0.0005 ) and negative ( p < 0.001 ) deflections in all signals are significant. The amplitude of both deflections are higher in NREM sleep than in wakefulness in the firing rate and V p . The opposite happens in the LFP signal. In all cases, the dashed vertical line, aligned at zero, indicates stimulus ξ onset. Stimulus duration is represented by a black horizontal line and light gray area.
Biology 10 00945 g003
Figure 4. Evoked responses in the two-cortical-column model where β = 1 and g ¯ AMPA = 2 . (A) Average LFP, firing rate and V p across 500 simulations for NREM sleep (left column) and wakefulness (right column) in the perturbed pyramidal population. Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus versus the prestimulus intervals. In all signals, as in Figure 3A, responses to the external stimulus consist of one positive and one negative significant deflection. The positive ( p < 0.0007 ) and negative ( p < 0.001 ) deflections are significant in NREM and wakefulness. Note that the coupling between the two cortical columns does not change the significance of positive and negative deflections in the evoked responses. (B) Results are as in panel (A) but for the case of the unperturbed pyramidal population. The evoked responses do not show significant deflections during NREM sleep and wakefulness in the LFP signal, neither in the firing rate nor in V p , where the increase in activity is not significantly higher than the spontaneous activity. In all cases, the dashed vertical line, aligned at zero, indicates stimulus ξ onset. Stimulus duration is represented by a black horizontal line and light gray area.
Figure 4. Evoked responses in the two-cortical-column model where β = 1 and g ¯ AMPA = 2 . (A) Average LFP, firing rate and V p across 500 simulations for NREM sleep (left column) and wakefulness (right column) in the perturbed pyramidal population. Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus versus the prestimulus intervals. In all signals, as in Figure 3A, responses to the external stimulus consist of one positive and one negative significant deflection. The positive ( p < 0.0007 ) and negative ( p < 0.001 ) deflections are significant in NREM and wakefulness. Note that the coupling between the two cortical columns does not change the significance of positive and negative deflections in the evoked responses. (B) Results are as in panel (A) but for the case of the unperturbed pyramidal population. The evoked responses do not show significant deflections during NREM sleep and wakefulness in the LFP signal, neither in the firing rate nor in V p , where the increase in activity is not significantly higher than the spontaneous activity. In all cases, the dashed vertical line, aligned at zero, indicates stimulus ξ onset. Stimulus duration is represented by a black horizontal line and light gray area.
Biology 10 00945 g004
Figure 5. Evoked responses in the two-cortical-column model for β = 1 , 2 and g ¯ AMPA = 2 , 6 . Average across 500 simulations of the firing rate of the unperturbed pyramidal population. Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus versus the prestimulus interval. Firing rate of the unperturbed population shows a positive deflection after stimulus onset in all cases. Increasing g ¯ AMPA from 2 to 6 evokes a significant response ( p < 0.003 ) in the unperturbed population only if β is also increased from 1 to 2. In all cases, the dashed vertical line, aligned at zero, indicates stimulus ξ onset. Stimulus duration is represented by a black horizontal line and light gray area.
Figure 5. Evoked responses in the two-cortical-column model for β = 1 , 2 and g ¯ AMPA = 2 , 6 . Average across 500 simulations of the firing rate of the unperturbed pyramidal population. Red horizontal bars show the location and length of significant t-cluster statistics when comparing the poststimulus versus the prestimulus interval. Firing rate of the unperturbed population shows a positive deflection after stimulus onset in all cases. Increasing g ¯ AMPA from 2 to 6 evokes a significant response ( p < 0.003 ) in the unperturbed population only if β is also increased from 1 to 2. In all cases, the dashed vertical line, aligned at zero, indicates stimulus ξ onset. Stimulus duration is represented by a black horizontal line and light gray area.
Biology 10 00945 g005
Figure 6. Information propagation as a function of β and g ¯ AMPA . (A) Example of extraction of t-cluster statistics when comparing the evoked responses in firing rate signals with the spontaneous activity in the prestimulus intervals at g ¯ AMPA = 2 and β = 5 . The yellow horizontal dashed line defines the critical t-value. The concatenation of consecutive t-values above and below this threshold that forms the t-cluster statistics is shown as rectangles in the bottom line. Significant and non significant t-cluster statistics are indicated in red and black, respectively. The area and duration of the first t-cluster statistic, shown in gray, are used to quantify response propagation to the unperturbed population. Stimulus duration is represented by a black horizontal line. (B) Changes in the area of the first significant t-cluster statistic in seconds (y-axis) with respect to β (x-axis). (C) Changes in the length of the first significant t-cluster statistic in milliseconds (y-axis) with respect to β (x-axis). Violet, dark violet and indigo correspond to g ¯ AMPA = 2 , g ¯ AMPA = 6 and g ¯ AMPA = 10 , respectively. Figures show that the propagated information to the unperturbed population depends on β . Increasing g ¯ AMPA from 2 to 10 and keeping β = 1 does not lead to a significant cluster in the unperturbed population (black squares). Increasing β from 1 to 3 shows a significant cluster (red squares) in all values of g ¯ AMPA . In the case of g ¯ AMPA = 2 , the t-cluster statistics are significant ( p < 0.001 ) for β = 2, 3, 4 and 5. In the case of g ¯ AMPA = 6 , the t-cluster statistics are significant ( p < 0.003 ) for β > 1 . In the case of g ¯ AMPA = 10 , the t-cluster statistics are significant ( p < 0.001 ) for β > 2 .
Figure 6. Information propagation as a function of β and g ¯ AMPA . (A) Example of extraction of t-cluster statistics when comparing the evoked responses in firing rate signals with the spontaneous activity in the prestimulus intervals at g ¯ AMPA = 2 and β = 5 . The yellow horizontal dashed line defines the critical t-value. The concatenation of consecutive t-values above and below this threshold that forms the t-cluster statistics is shown as rectangles in the bottom line. Significant and non significant t-cluster statistics are indicated in red and black, respectively. The area and duration of the first t-cluster statistic, shown in gray, are used to quantify response propagation to the unperturbed population. Stimulus duration is represented by a black horizontal line. (B) Changes in the area of the first significant t-cluster statistic in seconds (y-axis) with respect to β (x-axis). (C) Changes in the length of the first significant t-cluster statistic in milliseconds (y-axis) with respect to β (x-axis). Violet, dark violet and indigo correspond to g ¯ AMPA = 2 , g ¯ AMPA = 6 and g ¯ AMPA = 10 , respectively. Figures show that the propagated information to the unperturbed population depends on β . Increasing g ¯ AMPA from 2 to 10 and keeping β = 1 does not lead to a significant cluster in the unperturbed population (black squares). Increasing β from 1 to 3 shows a significant cluster (red squares) in all values of g ¯ AMPA . In the case of g ¯ AMPA = 2 , the t-cluster statistics are significant ( p < 0.001 ) for β = 2, 3, 4 and 5. In the case of g ¯ AMPA = 6 , the t-cluster statistics are significant ( p < 0.003 ) for β > 1 . In the case of g ¯ AMPA = 10 , the t-cluster statistics are significant ( p < 0.001 ) for β > 2 .
Biology 10 00945 g006
Table 1. Parameter description for the neural mass model.
Table 1. Parameter description for the neural mass model.
SymbolDescription
Q k max Maximal firing rate of population k
θ k Firing rate threshold of population k (sigmoid function half activation)
σ k Inverse neural gain of the sigmoid function of population k
τ k Membrane time constant of population k
C m Membrane capacitance in the HH model
ϕ Stochastic Gaussian process
N k k Mean number of synaptic connections from presynaptic
population k to postsynaptic population k
γ m Time constant of the postsynaptic response of synapse type m
g ¯ X k Average X-ergic conductance on population k
E X Reversal potential of the X-ergic current
g ¯ KNa Average conductance of sodium-dependent potassium channel
E K Nernst reversal potential of potassium channel
τ Na Time constant of sodium extrusion
α Na sodium influx through population firing rate
R pump Strength of the sodium pump
Na eq Resting state sodium concentration equilibrium
Table 2. Parameter values for the one-cortical-column model.
Table 2. Parameter values for the one-cortical-column model.
SymbolNREMWUnit
Q p max 0.030.03 ms 1
Q i max 0.060.06 ms 1
θ p , θ i −58.5−58.5mV
σ p 6.76.7mV
σ i 66mV
τ p , τ i 3030ms
C m 11 μ F / cm 2
ϕ 1.81 ms 1
N p p 160160-
N i p 4040-
N p i 160160-
N i i 4040-
γ p 70 × 10 3 70 × 10 3 ms 1
γ i 58.6 × 10 3 58.6 × 10 3 ms 1
g ¯ AMPA p 12ms
g ¯ AMPA i 12ms
g ¯ GABA p 12.294ms
g ¯ GABA i 12.313ms
E AMPA 00mV
E GABA −70−70mV
E L p −66−66mV
E L i −64−64mV
g ¯ KNa 1.91.9 mS / cm 2
E K −100−100mV
τ Na 1.71.7ms
α Na 22 mM · ms
R pump 0.090.09mM
Na eq 9.59.5mM
Table 3. Parameter values for the two-cortical-column model.
Table 3. Parameter values for the two-cortical-column model.
SymbolNREMWUnit
Q p max , Q p max 0.030.03 ms 1
Q i max , Q i max 0.060.06 ms 1
θ p , θ i , θ p , θ i −58.5−58.5mV
σ p , σ p 6.76.7mV
σ i , σ i 66mV
τ p , τ i , τ p , τ i 3030ms
C m 11 μ F / cm 2
ϕ 1.81 ms 1
N p p , N p p 160160-
N i p , N i p 4040-
N p i , N p i 160160-
N i i , N i i 4040-
N p p , N p p 88-
N i p , N i p 22-
γ p , γ p 70 × 10 3 70 × 10 3 ms 1
γ i , γ i 58.6 × 10 3 58.6 × 10 3 ms 1
g ¯ AMPA p , g ¯ AMPA p 12ms
g ¯ AMPA i , g ¯ AMPA i 12ms
g ¯ GABA p , g ¯ GABA p 1.0822.446ms
g ¯ GABA i , g ¯ GABA i 1.0662.445ms
E AMPA 00mV
E GABA −70−70mV
E L p , E L p −66−66mV
E L i , E L i −64−64mV
g ¯ KNa 1.91.9 mS / cm 2
E K −100−100mV
τ Na 1.71.7ms
α Na 22 mM · ms
R pump 0.090.09mM
Na eq 9.59.5mM
β 11-
Table 4. Average GABAergic conductance values for the two-cortical-column model where g ¯ AMPA p , i , p , i = 2 .
Table 4. Average GABAergic conductance values for the two-cortical-column model where g ¯ AMPA p , i , p , i = 2 .
Symbol β = 1 2345Unit
g ¯ GABA p , g ¯ GABA p 2.4462.5992.7522.9053.058ms
g ¯ GABA i , g ¯ GABA i 2.4452.5762.7082.842.971ms
Table 5. Average GABAergic conductance values for the two-cortical-column model where g ¯ AMPA p , i , p , i = 6 .
Table 5. Average GABAergic conductance values for the two-cortical-column model where g ¯ AMPA p , i , p , i = 6 .
Symbol β = 1 2345Unit
g ¯ GABA p , g ¯ GABA p 8.8699.3279.78610.24510.704ms
g ¯ GABA i , g ¯ GABA i 7.9728.3678.7619.1569.551ms
Table 6. Average GABAergic conductance values for the two-cortical-column model where g ¯ AMPA p , i , p , i = 10 .
Table 6. Average GABAergic conductance values for the two-cortical-column model where g ¯ AMPA p , i , p , i = 10 .
Symbol β = 1 2345Unit
g ¯ GABA p , g ¯ GABA p 15.29116.05516.8217.58518.349ms
g ¯ GABA i , g ¯ GABA i 13.49914.15714.81515.47316.131ms
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Razi, F.; Moreno-Bote, R.; Sancristóbal, B. Computational Modeling of Information Propagation during the Sleep–Waking Cycle. Biology 2021, 10, 945. https://0-doi-org.brum.beds.ac.uk/10.3390/biology10100945

AMA Style

Razi F, Moreno-Bote R, Sancristóbal B. Computational Modeling of Information Propagation during the Sleep–Waking Cycle. Biology. 2021; 10(10):945. https://0-doi-org.brum.beds.ac.uk/10.3390/biology10100945

Chicago/Turabian Style

Razi, Farhad, Rubén Moreno-Bote, and Belén Sancristóbal. 2021. "Computational Modeling of Information Propagation during the Sleep–Waking Cycle" Biology 10, no. 10: 945. https://0-doi-org.brum.beds.ac.uk/10.3390/biology10100945

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop