Computational Modeling of Information Propagation during the Sleep–Waking Cycle

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.


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.

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.

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: where τ p/i , I p/i L , I p/i AMPA and I p/i GABA 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]: where Q max p/i , θ p/i and σ p/i represent the maximal firing rate, firing rate threshold and inverse neural gain of population p or i, respectively (see Tables 1 and 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].  Tables 1 and 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 Tables 1 and 3 for symbol description and parameter values). ξ represents the transient external input impinging only in one cortical column.

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 kk , 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]: where N kk 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): 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: whereḡ k AMPA andḡ k GABA 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. 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 k X 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

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: whereḡ 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).

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ḡ KNa , the model generates the low-amplitudehigh-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 max p 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. Table 3. Parameter values for the two-cortical-column model. We have chosen to upscale the conductance of excitatory synapses by increasing the average AMPAergic conductance,ḡ k AMPA , to go from NREM sleep to wakefulness, which is in agreement with SHY. The average GABAergic conductance,ḡ k GABA , 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ḡ k AMPA 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).

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.

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: where N kp 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: 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ḡ k AMPA . First, we increase β from 2 to 5 in steps of one unit whileḡ k AMPA = 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ḡ k AMPA = 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 (see Tables 5 and 6 forḡ k GABA values atḡ k AMPA = 6 and 10, respectively).

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.

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.

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.

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 2c 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: 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.

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 twocortical-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), wherex 1 andx 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 2 1 + (n − 1)s 2 2 )/(2n − 1). s 2 1/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 ofx 2 with respect to signalx 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.

Quantification of Information Propagation
We quantified the propagation of the external input from the perturbed to the unperturbed population as a function of β andḡ k AMPA . Whenx 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ḡ k AMPA .

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

Spontaneous Activity across Brain States
Spontaneous activity differs between NREM sleep and wakefulness. In the model, it was obtained by potentiatingḡ 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.

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

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

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,ḡ 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). 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ḡ 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ḡ AMPA = 2 andḡ 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). Pyramidal Population of the Unperturbed Column Our results show that keeping β = 1 does not trigger a significant response in the unperturbed population even when potentiatingḡ 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ḡ 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ḡ 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). Increasingḡ 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ḡ AMPA . In the case ofḡ AMPA = 2, the t-cluster statistics are significant (p < 0.001) for β = 2, 3, 4 and 5. In the case ofḡ AMPA = 6, the t-cluster statistics are significant (p < 0.003) for β > 1. In the case ofḡ AMPA = 10, the t-cluster statistics are significant (p < 0.001) for β > 2.

Information Propagation as a Function of β andḡ AMPA
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.

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ḡ k AMPA 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,ḡ k AMPA , 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 6B,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.

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.

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 conflicts 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:

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. 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: with the currents defined by: Symbol descriptions and the parameter values are provided in Tables 1 and 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 the sodium pump and the cortical firing rate functions are given by: Symbol descriptions and the parameter values are provided in Table 1  Upscaling the conductance of excitatory synapses by increasing the average AMPAergic conductance,ḡ k AMPA , 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.
[Na] 3 eq + 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: The average GABAergic conductance of the pyramidal population and that of the inhibitory population during wakefulness are obtained by settingḡ p/i AMPA = 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: [Na] p/p = 3 A p/p · 3375 1 − A p/p , [Na] 3 eq + 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: The average GABAergic conductances on pyramidal and inhibitory populations during NREM sleep in the two-cortical-column scenario are obtained by settingḡ p/i AMPA = 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ḡ p/i AMPA = 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ḡ k AMPA are obtained by setting the corresponding β,ḡ p/i AMPA and the peaks of V p and V i values during Up states in the one-cortical-column scenario in the aforementioned equations.