## Abstract

The transition from wakefulness to sleep represents the most conspicuous change in behavior and the level of consciousness occurring in the healthy brain. It is accompanied by similarly conspicuous changes in neural dynamics, traditionally exemplified by the change from “desynchronized” electroencephalogram activity in wake to globally synchronized slow wave activity of early sleep. However, unit and local field recordings indicate that the transition is more gradual than it might appear: On one hand, local slow waves already appear during wake; on the other hand, slow sleep waves are only rarely global. Studies with functional magnetic resonance imaging also reveal changes in resting-state functional connectivity (FC) between wake and slow wave sleep. However, it remains unclear how resting-state networks may change during this transition period. Here, we employ large-scale modeling of the human cortico-cortical anatomical connectivity to evaluate changes in resting-state FC when the model “falls asleep” due to the progressive decrease in arousal-promoting neuromodulation. When cholinergic neuromodulation is parametrically decreased, local slow waves appear, while the overall organization of resting-state networks does not change. Furthermore, we show that these local slow waves are structured macroscopically in networks that resemble the resting-state networks. In contrast, when the neuromodulator decrease further to very low levels, slow waves become global and resting-state networks merge into a single undifferentiated, broadly synchronized network.

## Introduction

Falling asleep is a gradual phenomenon that includes major changes in behavior, a progressive disconnection from the environment, and marked shifts in the level and content of consciousness, which can vary from the vivid awareness of full-fledged wakefulness to the near-annihilation of experience in early nonrapid eye movement (NREM) sleep (Nir and Tononi 2010). Corresponding changes in the electroencephalogram (EEG) are well known. During wakefulness, the EEG is characterized by waves of low amplitude and high frequency. This kind of EEG pattern is also known as desynchronized, which should be taken to indicate merely the absence of massive, low-frequency synchronization.

Recently, unit and local field recordings have shown that the transition from wakefulness to slow wave sleep has a marked local component. Using high-density EEG in humans, it has been shown that individual slow waves have specific sites of origin and propagation, and that they rarely involve the entire cortical mantle but remain more localized (Massimini et al. 2004; Murphy et al. 2009). Moreover, multiunit and local field potential recordings from up to 10 regions simultaneously have shown that slow sleep waves, corresponding to synchronous off-periods in neural activity, tend to be local, occurring in just a few areas at a time, and are only rarely truly global, primarily during deep sleep at the beginning of the night (Nir et al. 2011). Finally, multiarray recordings in rodents have shown that off-periods, associated with a slowing of the local field potential, can appear in individual cortical regions even during wakefulness. These brief episodes of “local sleep” occur while the rest of cortex remains active in a typical waking mode, and the animal's behavior is one of active waking with eyes open (Vyazovskiy et al. 2011).

Studies with functional magnetic resonance imaging (fMRI) also reveal changes in resting-state functional connectivity (FC) between the awake state and different sleep stages measured as correlated blood oxygen level-dependent (BOLD) signals among different brain regions. During wakefulness, spontaneous activity (unrelated to specific stimuli or tasks) is highly structured into specific spatio-temporal patterns, known as resting-state networks (RSNs; Biswal et al. 1995; Greicius et al. 2003; Fox et al. 2005; Fransson 2005; Raichle and Mintun 2006; Fox and Raichle 2007; Rogers et al. 2007; Vincent et al. 2007). More recently, RSNs have also been found and studied with magnetoencephalography (MEG) in the seminal papers of de Pasquale et al. (2010, 2012), Mantini et al. (2011), and Spadone et al. (2012). Initial studies showed that certain RSNs are largely preserved during light NREM sleep (Horovitz et al. 2008; Larson-Prior et al. 2009). Furthermore, Larson-Prior et al. (2009) did not detect changes in the default mode network (DMN) even during deep sleep. However, more recent studies using simultaneous EEG/fMRI techniques point to substantial changes in the coupling of cortical areas when subjects enter deeper stages of sleep (N3; Sämann et al. 2011; Boly et al. 2012; Tagliazucchi et al. 2012). In particular, the seminal work of Sämann et al. (2011) reported changes in the DMN connectivity and between the DMN and its anticorrelated network already in sleep stage 1. Specifically, they found that the contributions of the posterior cingulate (PC) to the DMN decrease as sleep becomes deeper.

Here, we employ a large-scale brain model to evaluate changes in resting-state FC at the transition between the awake and sleeping mode, which is achieved by modeling a progressive decrease in arousal-promoting neuromodulators. The main aim of this paper was to study how the resting-state FC of the cortex changes with the emergence of slow sleep waves, rather than to provide a detailed realistic biophysical model of the process of falling asleep, including its cortical and subcortical cellular components. Hence, we capture the overall effects of a progressive decrease in arousal-promoting neuromodulation by a single parameter, which we take to represent the level of cholinergic modulation due to the fact that acetylcholine (Ach) levels are selectively reduced during slow wave sleep when compared with both wakefulness and REM sleep (Jasper and Tessier 1971; Vanini et al. 2012). The results show that, when neuromodulators are parametrically decreased, local slow waves occur, but the overall organization of resting-state networks does not change. In contrast, when the neuromodulator decrease to very low levels, slow waves become global and resting-state networks merge into a single undifferentiated, synchronized network.

## Materials and Methods

Experimentally, assessing changes in resting-state FC (fMRI or MEG/EEG based) during the process of falling asleep is made difficult by the speed and nonstationarity of the transitions in brain dynamics. To do so systematically, large-scale models that can derive resting-state FC from anatomical connectivity data and systematically vary physiological parameters can be useful (Honey et al. 2007; Ghosh et al. 2008; Deco et al. 2009, 2011; Cabral et al. 2011; Deco and Jirsa 2012). These models take advantage of realistic neuroanatomical information from the macaque cortex provided by the CoCoMac neuroinformatics (Kötter 2004) and from the human provided by diffusion MRI (Hagmann et al. 2008). In particular, Deco and Jirsa (2012) considered a global attractor network of multiple cortical areas based on spiking neurons and realistic 2-amino-3-(3-hydroxy-5-methyl-isoxazol-4-yl) propanoic acid (AMPA), N-Methyl-d-aspartate (NMDA), and gamma-aminobutyric acid (GABA) synapses. In this model, the best fit of empirical resting-state fMRI correlations was obtained when the brain network is critical, that is, at the border of a dynamical bifurcation point. Here, we extend and use this model for investigating the process of falling asleep. In the next subsections, we describe the architecture and dynamics of this model.

### Brain Model

In this paper, only the cortex was modeled, because the resting-state experimental data (neuroanatomical- and fMRI BOLD-based FC matrices) used for fitting the model during the wake state included only cortical areas. We used here the experimental data of Hagmann et al. (2008), who analyzed the cortex of 5 healthy right-handed male human subjects using diffusion spectrum imaging (DSI; Wedeen et al. 2005) and white matter tractography (Hagmann et al. 2008; Honey et al. 2009). Imaging was performed on a 3-T scanner using a diffusion-weighted single-shot echo planar imaging (EPI) sequence with a time repetition (TR) of 4200 ms and a time echo (TE) of 89 ms and a maximal b-value of 9000 s/mm2. Q-space was sampled over 129 points located inside a hemisphere. The axial field-of-view was set to 224 mm with an in-plane resolution of 2 mm, and 36 slices of 3-mm thickness were acquired in an acquisition time of 18 min.

The gray matter was first parcellated in 66 areas (33 areas per hemisphere, see Fig. 1). The structure of the underlying neuroanatomical connectivity matrix was extracted by expressing as coupling weights the density with which 2 different brain areas are connected through white matter fiber tracts. In other words, the coupling weights were proportional to a normalized number of detected tracts linking the 2 brain areas. The neuroanatomical matrix was then averaged across the 5 human subjects. Figure 1 shows graphically the structure of the connectivity matrix by encoding the strengths of the different connections as a color map. Figure 1 also reveals the small-world structure of the cortex through the presentation of clusters of varying size, obtained by reordering the different brain regions according to modules that have substantially denser connectivity within the module than with the complementary part of the network, consistent with previous findings (Bullmore and Sporns 2009).

Figure 1.

Anatomical connectome derived by Hagmann et al. (2007) using DSI averaged over 5 healthy subjects. (Left) The parcellation scheme dividing the cortex into 33 anatomically segregated regions in each hemisphere [adapted from Hagmann et al. (2008)]. (Middle-left) White matter tracts detected using DSI and tractography [adapted from Honey et al. (2007)]. (Middle-right) Schematic representation of the anatomical network, where regions are represented by red spheres placed at their center of gravity and the link's thickness is proportional to the number of fiber tracts detected in each connection. (Right-top) The coupling weights are proportional to the number of tracts detected. White color means that no fiber connecting the 2 corresponding regions was detected. Weights were normalized so that they are between 0 and 1. (Right-bottom) Distance between regions (in mm) given as the average length of the fibers connecting a pair of regions. The list of brain regions is reported in Table 2.

Figure 1.

Anatomical connectome derived by Hagmann et al. (2007) using DSI averaged over 5 healthy subjects. (Left) The parcellation scheme dividing the cortex into 33 anatomically segregated regions in each hemisphere [adapted from Hagmann et al. (2008)]. (Middle-left) White matter tracts detected using DSI and tractography [adapted from Honey et al. (2007)]. (Middle-right) Schematic representation of the anatomical network, where regions are represented by red spheres placed at their center of gravity and the link's thickness is proportional to the number of fiber tracts detected in each connection. (Right-top) The coupling weights are proportional to the number of tracts detected. White color means that no fiber connecting the 2 corresponding regions was detected. Weights were normalized so that they are between 0 and 1. (Right-bottom) Distance between regions (in mm) given as the average length of the fibers connecting a pair of regions. The list of brain regions is reported in Table 2.

The brain model consists of coupled local attractor networks as specified in the next subsection. The coupling between the different local attractor networks at each node is specified by the neuroanatomical human matrix described above. We consider here that the white matter tract connections between 2 distinct brain areas describe synaptic connections between pyramidal neurons in those areas. We weight those interareal connections by the strength specified in the neuroanatomical matrix (normalized numbers of fibers connecting those regions) and by a common scaling factor denoted by W that multiplies all interareal connection strengths. We estimated this scaling factor by fitting the simulated BOLD FC to the empirical BOLD FC matrix acquired from the same human subjects from whom we had obtained the neuroanatomical matrices (see Deco and Jirsa 2012, for details). Resting-state activity was obtained on a 3-T MRI scanner with a 32-channel head coil using a gradient echo sequence with EPI read-out. Thirty-two slices of 3.6 mm with in-plane resolution of 3.3 mm (field-of-view 212 mm) were sampled at a TR of 1920 ms using a TE of 3 m ms. The fMRI BOLD signal was sampled during 20 min in the absence of any stimulation or task while the subject was asked to remain awake and look at a fix point in front of him. Both fMRI and DSI data were acquired from the same 5 subjects in order to be absolutely consistent and increase so the confidence in the results. After regressing out the global signal (Fox et al. 2005) and averaging across subjects, we obtained an empirical FC matrix (see Honey et al. 2009 for details). We used averaged values over the 5 subjects in order to minimize individual differences and artifacts that could appear because of the sporadic lack of vigilance. This empirical FC matrix reflects the correlation of the BOLD activity between different brain areas at rest. Next, after simulating in the model the firing activity of all 66 brain areas, we derived the BOLD fMRI signal by using the Balloon–Windkessel model described below. The simulated BOLD signal was also down sampled at 2 s to achieve the same resolution as in Honey et al. (2009), and the global signal was regressed out of the BOLD time series using the same procedure as in the experiments. Finally, we computed the simulated FC by calculating the correlation matrix of the BOLD activity between all brain areas. To identify the region of the parameter W for which the model best reproduced the empirical FC, we computed the Pearson correlation between the empirical and the simulated BOLD FC matrix. The best fit was obtained for W = 1.6.

### Single Cortical Area Model

Each local cortical area is modeled in a biophysically realistic framework of spiking attractor networks (Rolls and Deco 2010). In our case, the spiking attractor network consists of integrate-and-fire (IF) spiking neurons with excitatory (AMPA and NMDA) and inhibitory (GABA-A) synaptic receptor types (Brunel and Wang 2001). This attractor network constitutes a dynamical system specified by a set of coupled equations describing each neuron and synapse. This type of dynamical system has stationary fixed points called “attractors,” which are commonly expressed as stable patterns of firing activity. Furthermore, transitions between different stable attractors can be driven by noise that appears in the form of finite-size effects that amplify the effect of external stochastic background noise (as described below in eq. 8). Finite-size effects are due to the fact that the populations are described by a finite number, N, of neurons (Mattia and Del Giudice 2004).

The architecture of each local cortical area attractor network is given by a fully connected recurrent network of a population of, NE = 100, excitatory pyramidal neurons and a population of, NI = 100, inhibitory neurons. The recurrent self-excitation within each excitatory population is given by the weight w+, and within each inhibitory population is given by the weight w= 1. The connections between excitatory and inhibitory neurons have the weight w= 1.

Spiking activity of the neurons is described by the classical IF model. This model expresses the dynamical evolution of the membrane potential V(t) driven by incoming input currents coming from connected neurons or external inputs. An IF neuron consists of a basic resistor–capacitor circuit. Specifically, the membrane potential of each neuron in the network is given by the following equation:

$CmdV(t)dt=−gm[V(t)−VL]−IM−gAMPA,ext[V(t)−VE]∑j=1NextsjAMPA,ext(t)−gAMPA,rec[V(t)−VE]∑j=1NEwjsjAMPA,rec(t)−gNMDA[V(t)−VE]1+γe−βV(t)∑j=1NEwjsjNMDA(t)−gGABA[V(t)−VI]∑j=1NIwjsjGABA(t).$
If the membrane potential is below a given threshold Vthr (subthreshold dynamics). When the voltage across the membrane reaches the threshold Vthr, the neuron generates a spike. The spike is transmitted to other neurons and the membrane potential is instantaneously reset to Vreset and maintained there for a refractory time τref during which the neuron is unable to produce further spikes. In equation (1), gm is the membrane leak conductance, Cm, the capacity of the membrane, and VL, the resting potential. The membrane time constant is defined by τm= Cm/gm. The synaptic input current is given by the last 4 terms on the right-hand side of equation (1). The spikes arriving at the synapse produce a postsynaptic excitatory or inhibitory potential given by a conductance-based model specified by the synaptic receptors. The 4 terms correspond to: Glutamatergic AMPA (IAMPA,ext) external excitatory currents, AMPA (IAMPA,rec), and NMDA (INMDA) recurrent excitatory currents, and GABAergic recurrent inhibitory currents (IGABA). The respective synaptic conductances are gAMPA,ext, gAMPA,rec, gNMDA, and gGABA, and VE and VI are the excitatory and inhibitory reversal potentials, respectively. The dimensionless parameters wj of the connections are the synaptic weights described above. The gating variables $sji(t)$are the fractions of open channels of neurons and are given by:
$dsjI(t)dt=−sjI(t)τI+∑kδ(t−tjk),forI=AMPA or GABA,$

$dsjNMDA(t)dt=−sjNMDA(t)τNMDA,decay+αxjNMDA(t)[1−sjNMDA(t)],$

$dxjNMDA(t)dt=−xjNMDA(t)τNMDA,rise+∑kδ(t−tjk).$

The sums over the index k represent all the spikes emitted by the presynaptic neuron j (at times $tjk$). In equations (7–11), τNMDA,rise and τNMDA,decay are the rise and decay times for the NMDA synapses, and τAMPA and τGABA the decay times for AMPA and GABA synapses. The rise times of both AMPA and GABA synaptic currents are neglected because they are short (<1 ms). The values of the constant parameters and default values of the free parameters used in the simulations are displayed in Table 1.

Table 1

Neural and synaptic parameters for each local brain area

Excitatory neurons

Inhibitory neurons

Synapses

NE 100 neurons NI 100 neurons VE 0 mV
Cm 0.5 nF Cm 0.2 nF VI −70 mV
gm 25 nS gm 20 nS τAMPA 2 ms
VL −70 mV VL −70 mV τNMDA,rise 2 ms
Vthr −50 mV Vthr −50 mV τNMDA,decay 100 ms
Vreset −55 mV Vreset −55 mV τGABA 10 ms
τref 2 ms τref 1 ms α 0.5 kHz
gAMPA,ext 2.496 nS gAMPA,ext 1.944 nS β 0.062
gAMPA,rec 0.104 nS gAMPA,rec 0.081 nS γ 0.28
gNMDA,rec 0.327 nS gNMDA,rec 0.258 nS
gGABA 4.375 nS gGABA 3.4055 nS
Excitatory neurons

Inhibitory neurons

Synapses

NE 100 neurons NI 100 neurons VE 0 mV
Cm 0.5 nF Cm 0.2 nF VI −70 mV
gm 25 nS gm 20 nS τAMPA 2 ms
VL −70 mV VL −70 mV τNMDA,rise 2 ms
Vthr −50 mV Vthr −50 mV τNMDA,decay 100 ms
Vreset −55 mV Vreset −55 mV τGABA 10 ms
τref 2 ms τref 1 ms α 0.5 kHz
gAMPA,ext 2.496 nS gAMPA,ext 1.944 nS β 0.062
gAMPA,rec 0.104 nS gAMPA,rec 0.081 nS γ 0.28
gNMDA,rec 0.327 nS gNMDA,rec 0.258 nS
gGABA 4.375 nS gGABA 3.4055 nS
Table 2

Names and abbreviations of the brain regions considered in the human connectome from Hagmann et al. (2008) (in alphabetical order)

Abbreviations Brain region
BSTS Bank of the superior temporal sulcus
CAC Caudal anterior cingulate cortex
CMF Caudal middle frontal cortex
CUN Cuneus
ENT Entorhinal cortex
FP Frontal pole
FUS Fusiform gyrus
IP Inferior parietal cortex
ISTC Isthmus of the cingulate cortex
IT Inferior temporal cortex
LING Lingual gyrus
LOCC Lateral occipital cortex
LOF Lateral orbitofrontal cortex
MOF Medial orbitofrontal cortex
MT Middle temporal cortex
PARC Paracentral lobule
PARH Parahippocampal cortex
PC Posterior cingulate cortex
PCAL Pericalcarine cortex
PCUN Precuneus
POPE Pars opercularis
PORB Pars orbitalis
PREC Precentral gyrus
PSTC Postcentral gyrus
PTRI Pars triangularis
RAC Rostral anterior cingulate cortex
RMF Rostral middle frontal cortex
SF Superior frontal cortex
SMAR Supramarginal gyrus
SP Superior parietal cortex
ST Superior temporal cortex
TP Temporal pole
TT Transverse temporal cortex.
Abbreviations Brain region
BSTS Bank of the superior temporal sulcus
CAC Caudal anterior cingulate cortex
CMF Caudal middle frontal cortex
CUN Cuneus
ENT Entorhinal cortex
FP Frontal pole
FUS Fusiform gyrus
IP Inferior parietal cortex
ISTC Isthmus of the cingulate cortex
IT Inferior temporal cortex
LING Lingual gyrus
LOCC Lateral occipital cortex
LOF Lateral orbitofrontal cortex
MOF Medial orbitofrontal cortex
MT Middle temporal cortex
PARC Paracentral lobule
PARH Parahippocampal cortex
PC Posterior cingulate cortex
PCAL Pericalcarine cortex
PCUN Precuneus
POPE Pars opercularis
PORB Pars orbitalis
PREC Precentral gyrus
PSTC Postcentral gyrus
PTRI Pars triangularis
RAC Rostral anterior cingulate cortex
RMF Rostral middle frontal cortex
SF Superior frontal cortex
SMAR Supramarginal gyrus
SP Superior parietal cortex
ST Superior temporal cortex
TP Temporal pole
TT Transverse temporal cortex.

We included in the IF-model spike frequency adapting mechanisms, including slow voltage-dependent potassium currents M-currents IM (Brown and Adams 1980) given by:

$IM=gMa(t)[V(t)−VK],$

In the last equation, VK= −80 mV is the reversal potential of the potassium channel. M-currents are mainly activated by supra-threshold membrane potential, and therefore they are implemented by incrementing the adaptation variable a(t) after each action potential by a small amount (α= 0.1), so that IM is incremented accordingly (Benda and Herz 2003). Between spikes, the a(t) dynamics is modeled as a leaky integrator with a decay constant τM= 500 ms. Thus, the M-current's dynamics can be described by the following system of equations:

$da(t)dt=−a(t)τM,$

$ifV(t)=Vthrthena(t)=a(t)+αandV=Vreset.$

Furthermore, all neurons in the network receive an external background input of uncorrelated Poisson spike trains with a time-varying rate $vextp(t)$, governed by:

$τndvextp(t)dt=−[vextp(t)−v0]+σv2τnnp(t),$
where τn= 30 ms, ν0= 2.4 kHz, σν is the standard deviation of $vextp(t)$, and np(t), normalized Gaussian white noise. Negative values of $vextp(t)$ that could arise due to the noise term are rectified to zero.

The emergence of high-structured spatio-temporal patterns across the brain at rest in neuroimaging experiments is typically characterized by the resting FC matrix. The FC reflects, in fact, the segregation of the activity of the idle brain in resting-state networks. In our model, the mutual interaction between the different local brain area dynamics is shaped by the underlying neuroanatomical structure. As illustrated in the next section, in order to adjust the global coupling between different brain areas, we used the empirical FC BOLD fMRI activity obtained in healthy humans during rest. It has been shown that, for this detailed spiking model, the best fit of empirical resting fMRI data is obtained when the brain network is critical, that is, the global scaling factor W is at the border of a dynamical bifurcation point, so that at that operating working point the system defines a meaningful dynamic repertoire that is inherent in the neuroanatomical connectivity (Deco and Jirsa 2012). Indeed, the resting FC emerges from the noisy fluctuations of the spiking activity of each brain area around a low activity spontaneous state that get correlated (between areas) through the underlying anatomical structural connectivity matrix. As described above, to identify the critical region of the parameter W for which the model best reproduced the empirical FC, we computed the Pearson correlation between the empirical and simulated BOLD FC matrix. Hence, it is necessary to convert spiking activity as simulated in the attractor network model to a BOLD fMRI signal. For this purpose, we used the Balloon–Windkessel hemodynamic model of Friston et al. (2003). The Ballon–Windkessel model describes the coupling of perfusion to the BOLD signal, with a dynamical model of the transduction of neural activity into perfusion changes. The model assumes that the BOLD signal is a static nonlinear function of the normalized total deoxyhemoglobin voxel content, normalized venous volume, resting net oxygen extraction fraction by the capillary bed, and resting blood volume fraction. The BOLD signal estimation for each brain area is computed by the level of neuronal activity summed over all neurons in both populations (excitatory and inhibitory populations) in that particular area. We also regressed out the global signal from the simulated BOLD time series in all calculations.

### Induction of Sleep: Acetylcholine Level

In the model, sleep was induced by simulating a progressive decrease in the level of arousal-promoting neuromodulatory inputs from subcortical systems. We capture the overall effects of a progressive decrease in arousal-promoting neuromodulation by a single parameter ζ, which we take to represent the level of cholinergic modulation. Although Ach affects cortical processing in a multitude of manners, we restricted our analysis to 2 of its main established actions, namely: (1) High ACh levels reduce the magnitude of cortico-cortical excitatory recurrent interactions through presynaptic inhibition of glutamate release via muscarinic receptors (Hasselmo and Bower 1992), and (2) high Ach level reduces spike frequency adaptation due to its effect on muscarinic M-currents (Hasselmo and Giocomo 2006). Both variables induce slow sleep oscillations in the system. By grading continuously the level of Ach, we can study carefully the transition from awake resting state to slow wave sleep. More specifically, we simulated 10 min of spontaneous dynamics for 35 different levels of Ach grading the level of sleep from the sleep mode (ζ= 0) to the awake mode (ζ= 1). Therefore, we scale the values of the recurrent excitatory weights w+ and of the intercortical weights W, and the value of the adaptation M-currents, as follows:

$w+=3−1.5ζW=3.2−1.6ζgM=9−9ζ$

Note that the values for ζ= 1 (wake mode) correspond to the working point that better fits the resting-state BOLD FC. In particular, we choose the parameters such that for the transition from the wake mode to the deepest sleep mode (ζ= 0) increases the level of adaptation and the recurrent intra- and intercortical excitatory waves to maximum values able to sustain slow sleep waves in most of the cortical areas. More specifically, for the wake mode (ζ= 1), the parameters in equation (9) are such that the wake resting state is modeled following the specific procedure of Deco and Jirsa (2012), that is, there is no adaptation (gM = 0) and for w+ = 1.5, the optimal W for fitting the empirical resting BOLD FC is obtained at the edge of the bifurcation (W = 1.6). For the sleep mode (ζ= 0), the adaptation gM = 9 was choosen such that slow sleep oscillations are properly created for the increased values of w+ = 3 and W = 3.2.

### Analysis of the Slow Sleep Waves

The analysis of the slow sleep waves is based on the simulated firing rates. In all calculations, we first simulated 10 min of the firing rate. The slow sleep waves are based on the excitatory populations. For each brain area, we filtered each excitatory population rate by a band pass filter between 0.5 and 4 Hz in order to focus on the slow sleep wave band. We also calculated the Hilbert transform to obtain the envelopes and phases of the power fluctuations of slow sleep waves. In concrete, we used the so-called “analytic signal” of the firing rate data. The analytical signal of the firing rate r(t) is: R(t) = r(t) + i H(t), that is, it is the complex number whose real part is the original data, and the imaginary part is the Hilbert transform H(t) of the data. The imaginary part is a version of the original real sequence with a 90° phase shift. The analytic signal can be expressed in polar coordinates as: R(t) = A(t) exp[(t)], where A(t) is the “envelope” or amplitude of the analytic signal, and ψ the “phase” of the analytic signal. Here, for characterizing the spatial structure of the slow sleep waves, we calculated the correlation matrix between the envelopes of the firing rates of the excitatory populations of all areas. We also calculated the correlation matrix between the corresponding phases in order to study their coherence. We also calculated the power spectrum of each time series using the Fast Fourier Transform and evaluated the power in the slow wave frequency range for each area and each condition.

### Estimation of the State Repertoire

In the last section of the data analysis, the repertoire of brain states was estimated from the variance of FC of parcellated brain regions. The variance of connectivity is calculated for all pair-wise correlation coefficients of between the BOLD signals of the parcellated brain regions. The variance of the strength of regional connectivity is presumed to reflect the repertoire of states accessed by the brain over time. All pair-wise correlation analysis was performed, followed by the calculation of variance of the obtained correlation coefficients for each level of Ach, that is, for each level of sleep. Voxel correlations were calculated for the whole time period of 10 min. First, the 66 regions from our parcellated functionally connected network were selected. All pair-wise correlation coefficients among the simulated BOLD signal of the selected regions were calculated for each Ach level. The variance of the so-obtained correlation coefficients characterizes the heterogeneity of resting-state brain connectivity, which presumably relates to the repertoire of brain states at rest. As previously argued (Tononi 2008), highly similar stereotypic activity across many brain regions or voxels (low variance) would imply that the repertoire of discriminable brain states was relatively small. Conversely, an increased heterogeneity of correlations (high variance) would entail a large repertoire of brain configurations or brain sates. A large repertoire translates to large information capacity and is thought to be a prerequisite of the wakeful conscious state (Tononi 2008).

## Results

Our main goal was to examine how local slow sleep waves emerge and evolve in the transition and how they are spatially structured in the cortex. Figure 2 (top) shows for each cortical area (in half hemisphere) the envelopes of slow sleep waves (based on the underlying firing rate) as a function of time for a snapshot of 10 s (after the first 10 min) and for 5 different Ach levels (ζ = 1, 0.71, 0.57, 0.42, 0, from left to right, respectively). Starting from the left side, the first panel corresponds to the wake mode and the subsequent panels reflect progressive “falling asleep” until a full-fledged sleep mode is reached. From the panels, it is easy to realize the presence of strong slow waves during sleep, which are globally synchronized for the deep sleep state. Reducing the level of Ach (i.e. moving toward the sleep mode) results at first in the emergence of local slow waves (in space and time), in line with the experimental observations of Vyazovskiy et al. (2011). Further reducing the level of Ach increases the incidence of slow waves, until in the full-fledged sleep mode (ζ= 0) the local slow waves become fully synchronized across cortical areas (rightmost panel).

Figure 2.

Top subpanels: Envelopes of the slow sleep oscillations of each cortical area (based on the firing rates) as a function of time, for 5 snapshot of 10 s (after the first 10 min). Only half hemisphere is shown. The 5 snapshots correspond to different Ach levels, namely: ζ = 1, 0.71, 0.57, 0.42, 0, and 1, from left to right, respectively. Increasing the level of Ach (i.e. getting awake) reduces the presence of slow sleep waves, so that only local (in space and time) slow sleep waves emerge according to the experimental observations of Vyazovskiy et al. (2011). Bottom panel (left): Total power (summed over all cortical areas) of slow sleep waves as a function of the Ach level. Bottom panel (right): Power of slow sleep waves for each cortical area as a function of the Ach level. Awakening decreases the power of the slow sleep waves.

Figure 2.

Top subpanels: Envelopes of the slow sleep oscillations of each cortical area (based on the firing rates) as a function of time, for 5 snapshot of 10 s (after the first 10 min). Only half hemisphere is shown. The 5 snapshots correspond to different Ach levels, namely: ζ = 1, 0.71, 0.57, 0.42, 0, and 1, from left to right, respectively. Increasing the level of Ach (i.e. getting awake) reduces the presence of slow sleep waves, so that only local (in space and time) slow sleep waves emerge according to the experimental observations of Vyazovskiy et al. (2011). Bottom panel (left): Total power (summed over all cortical areas) of slow sleep waves as a function of the Ach level. Bottom panel (right): Power of slow sleep waves for each cortical area as a function of the Ach level. Awakening decreases the power of the slow sleep waves.

The bottom subpanels of Figure 2 show the power of the slow sleep waves as a function of the Ach level (left total power summed over all areas and right for each area). In addition, slow waves are absent in the wake mode, then local slow waves begin to appear, and finally in deep sleep (low Ach level) high power in the slow sleep wave range is evident all over the cortex. Note that the slow sleep waves are also present at the spiking level even without any particular filtering. Figure 3 shows for a high (left) and low (right) Ach level the spiking activity of the model in a particular area (lPC) during 8 s. The top panel shows the spiking activity of 20 randomly chosen neurons of the lPC excitatory population. The middle panel shows the corresponding firing rate activity calculated using all neurons on the excitatory population of the lPC. These panels also show the results of the filtering and Hilbert transform, which confirms that the procedure properly captures the slow wave events (amplitude and phases). The bottom panels show the corresponding power spectrum.

Figure 3.

Firing rate activity of one particular area (lPC). Left (right) subpanels correspond to high (low) Ach levels (i.e. to awake, deep sleep phases, respectively). Top panels show the spiking activity of 20 randomly chosen neurons of the lPC excitatory population. The middle panels show the corresponding firing rate activity (in black) calculated using all neurons on the excitatory population of the lPC (and no filtering at all). These panels also show the results of the filtering and Hilbert transform (in red), which confirms that the procedure properly captures the slow waves events (amplitude and phases). The bottom panels show the corresponding power spectrum.

Figure 3.

Firing rate activity of one particular area (lPC). Left (right) subpanels correspond to high (low) Ach levels (i.e. to awake, deep sleep phases, respectively). Top panels show the spiking activity of 20 randomly chosen neurons of the lPC excitatory population. The middle panels show the corresponding firing rate activity (in black) calculated using all neurons on the excitatory population of the lPC (and no filtering at all). These panels also show the results of the filtering and Hilbert transform (in red), which confirms that the procedure properly captures the slow waves events (amplitude and phases). The bottom panels show the corresponding power spectrum.

Figure 4 plots the temporal correlation of the slow wave envelopes (based on the firing rates) between different areas and of the BOLD signals. The top panel shows the correlation of the envelopes of the slow sleep waves between a specific seed lPC and all cortical areas (in half hemisphere), for the same 5 snapshots corresponding to different levels of Ach (ζ = 1, 0.71, 0.57, 0.42, 0, from left to right, respectively). In the early sleep stages, at high levels of Ach, the spatial structure of the correlation between the envelopes of the slow sleep waves reflects the RSNs (i.e. the BOLD FC) and thus, the underlying neuroanatomical connectivity (This can be seen below in Fig. 6, more specifically). For low level of Ach (deep sleep), all envelopes are highly correlated due to the underlying global synchronization, that is, no particular neuroanatomical structure is reflected at this stage. The second row panel shows the correlation of the underlying fMRI BOLD signal, that is, the BOLD FC, for the same seed and for the same 5 snapshots. In addition, in the wake mode (leftmost subpanel), one observes the specific resting-state connectivity that is observed experimentally (see Deco and Jirsa 2012), whereas during deep sleep (rightmost subpanel) the resting-state connectivity is obliterated (but not completely lost), consistently with the observations of Sämann et al. (2011). Note that lPC is the main seed corresponding to the DMN and according to Sämman et al. (2011) it is one of the areas that get mainly decorrelated with the rest of the DMN during sleep. Figure 5 plots, for the wake mode (Ach level ζ= 1), the experimentally measured BOLD FC (left) and the simulated one (right) for the same seed. The model at the optimal critical point W = 1.6 (see Materials and Methods for details) is able to reproduce the resting-state FC with reasonable accuracy.

Figure 4.

Temporal correlation of the slow sleep envelopes (based on the firing rates) between different nodes and of the BOLD signal. First row top panels: Correlation of the envelopes of the slow sleep waves between a specific seed (left posterior cingulate, lPC) and all cortical areas, for 5 snapshots corresponding to different levels of Ach (ζ = 0, 0.42, 0.57, 0.71, and 1, from left to right, respectively). Second row panels: Correlation of the underlying fMRI BOLD signal, for the same seed and for the same 5 snapshots. Only half hemisphere is shown. Bottom panels (left): Total correlation (summed over all areas) between the envelopes of the slow sleep waves for the same seed as a function of the Ach level. Bottom panel (right): Correlation between the envelopes of the slow sleep waves for each cortical area and for the same seed as a function of the Ach level.

Figure 4.

Temporal correlation of the slow sleep envelopes (based on the firing rates) between different nodes and of the BOLD signal. First row top panels: Correlation of the envelopes of the slow sleep waves between a specific seed (left posterior cingulate, lPC) and all cortical areas, for 5 snapshots corresponding to different levels of Ach (ζ = 0, 0.42, 0.57, 0.71, and 1, from left to right, respectively). Second row panels: Correlation of the underlying fMRI BOLD signal, for the same seed and for the same 5 snapshots. Only half hemisphere is shown. Bottom panels (left): Total correlation (summed over all areas) between the envelopes of the slow sleep waves for the same seed as a function of the Ach level. Bottom panel (right): Correlation between the envelopes of the slow sleep waves for each cortical area and for the same seed as a function of the Ach level.

Figure 5.

Detailed comparison between the empirical and the simulated resting-state BOLD FC for the wake state (Ach level ζ = 1) and for an specific seed (left posterior cingulate, lPC). Left: Empirical resting-state BOLD FC. Right: Simulated resting-state BOLD FC. Both hemispheres are shown.

Figure 5.

Detailed comparison between the empirical and the simulated resting-state BOLD FC for the wake state (Ach level ζ = 1) and for an specific seed (left posterior cingulate, lPC). Left: Empirical resting-state BOLD FC. Right: Simulated resting-state BOLD FC. Both hemispheres are shown.

The bottom subpanels of Figure 4 show the correlation between the envelopes of slow sleep waves for the same seed as a function of the Ach level (left total correlation summed over all areas and right for each area). Note that at higher levels of Ach (wake), slow sleep waves are indeed few (Fig. 2, bottom right, which also shows that the few extant slow waves are also localized). In addition, Figure 4 shows that, for high levels of Ach (wake), slow sleep waves are essentially uncorrelated (in Fig 4, bottom subpanels, at high levels of Ach, the correlation between the envelopes of slow sleep waves is indeed very small).

In the wake mode (higher Ach level), the spatial correlation between slow sleep waves resembles the pattern seen in the resting-state BOLD FC, as shown quantitatively in Figure 6. When transitioning to deep sleep (low Ach level), one observes a higher correlation between the seed and the rest of the cortex resulting from the underlying global synchronization.

Figure 6.

Top left panel: Correlation between the experimentally obtained resting-state BOLD FC (wakefulness only) and the envelopes of the slow sleep waves as a function of sleep stage. Top right panel: Correlation between the phases of the slow sleep waves and the same empirical resting-state BOLD FC. Bottom left panel: Correlation between the empirical and the simulated BOLD FC. Bottom right panel: Correlation between the simulated BOLD signal and the envelopes of the slow sleep waves. During awakening (increasing the Ach level), the correlation between slow sleep waves gets more shaped according to the resting-state FC.

Figure 6.

Top left panel: Correlation between the experimentally obtained resting-state BOLD FC (wakefulness only) and the envelopes of the slow sleep waves as a function of sleep stage. Top right panel: Correlation between the phases of the slow sleep waves and the same empirical resting-state BOLD FC. Bottom left panel: Correlation between the empirical and the simulated BOLD FC. Bottom right panel: Correlation between the simulated BOLD signal and the envelopes of the slow sleep waves. During awakening (increasing the Ach level), the correlation between slow sleep waves gets more shaped according to the resting-state FC.

During the process of falling asleep, slow waves that appear locally—in some cortical areas and not others—in line with what reported by Vyazovskiy et al. (2011)—are structured in agreement with the observed resting-state FC and thus, reflect the underlying anatomical structure. In other words, the correlation structures between the envelopes of the slow sleep waves resemble that of the RSNs, that is, of the BOLD FC. This can be seen explicitly in the results shown in Figure 6. The top left panel of Figure 6 plots the correlation between the experimentally obtained resting-state FC correlation matrix (with fMRI BOLD) and the correlation matrix between the envelopes of the slow sleep waves (based on the firing rates). This panel shows that at first there is an increase in the correlation between the BOLD awake FC, which reflects the underlying anatomical connectivity, and the envelopes of slow sleep waves, which remain localized. Then, when the envelopes of slow sleep waves become globally synchronized, the FC structure of the envelopes of the slow sleep waves is progressively obliterated and finally destroyed. This figure demonstrates that local slow sleep waves first become correlated and synchronized locally according to the neuroanatomical connectivity. Thus, early in the falling asleep process, overall brain dynamics are preserved in spite of local populations of neurons going “off-line.” The top right panel shows the evolution of the correlation between the correlation matrix between the phases of the slow sleep waves (extracted with the Hilbert transform) and the empirical resting-state BOLD FC correlation matrix. The bottom left figure shows that the empirical and the simulated FCs are even more correlated in the wake mode. In fact, the simulations show that the resting-state structure of the BOLD signal is progressively obliterated when falling asleep, in line with the experimental results of Sämann et al. (2011). The bottom right panel plots the correlation between the simulated resting-state BOLD FC and the correlation matrix between the envelopes of the slow sleep waves. In addition, BOLD signal and slow waves go hand in hand during the different sleep stages. Note that while FC is globally reduced, the correlation between the envelopes of slow sleep waves is globally increased (Fig. 4, bottom subpanels).

Finally, we examined how the variance of connectivity is altered by a graded reduction in cholinergic neuromodulation. The variance of the strength of regional connectivity reflects the heterogeneity of temporal interactions of brain regions and indirectly, the repertoire of brain states. The repertoire of brain states has been proposed to be a critical determinant of the level of consciousness (Tononi 2008). From the simulated BOLD time courses, the repertoire of brain states was estimated as the average variance of all pair-wise correlation coefficients of the BOLD signals of parcellated brain regions, repeated at various simulated Ach levels. The computational results in Figure 7 clearly show that the variance of connectivity was reduced as a function of the Ach level in a graded manner. A particularly sharp drop was observed at an intermediate Ach level (0.7–0.5), nominally associated with a transition to deep sleep. This transition appears to correspond to a global synchronization of neuronal activity, which can thus be taken as an indication of a diminished repertoire of available brain states. In other words, the analysis of the variability of RSNs during the falling asleep process reveals a parametric decrease in the repertoire of transient networks that parallels the decrease in neuromodulators.

Figure 7.

Variance of connectivity as a function of sleep stage. The variance of connectivity is calculated for all pair-wise correlation coefficients of between the BOLD signals of the parcellated brain regions. The variance of the strength of regional connectivity is presumed to reflect the repertoire of states accessed by the brain over time.

Figure 7.

Variance of connectivity as a function of sleep stage. The variance of connectivity is calculated for all pair-wise correlation coefficients of between the BOLD signals of the parcellated brain regions. The variance of the strength of regional connectivity is presumed to reflect the repertoire of states accessed by the brain over time.

In summary, we have shown using a large-scale model of cortico-cortical connectivity that mimicking the process of falling asleep results in the progressive appearance of slow sleep waves that are at first local, that is, confined to a few cortical areas at a given time, and then become more and more globally synchronous, in line with experimental results obtained in humans (Nobili et al. 2012). The FC expressed by the cortical model, which reflects the underlying anatomical connectivity derived from the DSI data, is in substantial agreement with the resting-state FC obtained experimentally in the wake mode. This specific FC is preserved when slow waves are still predominantly local, but is obliterated when slow waves become globally synchronized in deep sleep. In parallel with these gradual changes, the variability of functionally connected regions is reduced, suggesting a diminution of the repertoire of brain states as the brain transitions to deep sleep.

## Discussion

The characteristic EEG activity patterns that distinguish wakefulness from slow wave sleep are generated by the cortico-thalamic system (Steriade et al. 1997; Robinson et al. 2002, 2011; Nunez and Srinivasan 2006; Roberts and Robinson 2012). Furthermore, discrete populations of cells in the brainstem and hypothalamus provide ongoing neuromodulatory input to the cortex that promotes arousal (Moruzzi and Magoun 1949, Steriade and McCarley 1990; Robinson et al. 2002, 2011; Roberts and Robinson 2012). When the level of neuromodulators such as Ach, noerepinephrine, histamine, and orexin/hypocretin decreases, maintaining arousal becomes difficult, and the cortico-thalamic system tends to settle into a bistable pattern of activity, alternating between depolarized UP states and hyperpolarized DOWN states that constitute the sleep slow oscillation (Steriade et al. 2001). As more and more neurons become bistable, synchronized UP and DOWN states begin to appear, which are reflected in the EEG as slow sleep waves. As in a previous, detailed model of wake and sleep in the cortico-thalamic system (Hill and Tononi 2005; Esser et al. 2007), we here simulated the decrease in the level of arousal-promoting neuromodulators by progressively modifying some parameters governing the activation of cortical neurons. Our main goal has been to study how the resting-state FC of the cortex changes with the emergence of slow sleep waves. Hence, we captured the overall effects of a progressive decrease in arousal-promoting neuromodulation by a single parameter, which we took to represent the level of cholinergic modulation.

In the large-scale model of cortico-cortical connectivity presented in this paper, when the levels of arousal-promoting neuromodulators (Ach) are reduced, cortical neurons enter a bistable, on–off mode of firing, first locally and then globally, consistent with experimental observations (Steriade et al. 2001) and previous modeling work (Bazhenov et al. 2002; Hill and Tononi 2005; Esser et al. 2007). At high levels of ACh, neurons fire in a tonic manner, without off-periods and with no sign of low-frequency synchronization. Under these conditions, the resting-state FC obtained by convolving simulated neuronal activity with BOLD hemodynamics resembles empirical results obtained in awake humans. When Ach are slightly reduced, off-periods begin to appear in individual cortical areas at the spiking level, but no long-distance low-frequency synchronization is observed. As long as off-periods remain local, resembling experimentally observed local sleep (Vyazovskiy et al. 2011), BOLD resting-state FC patterns do not change appreciably, indicating that, overall, cortico-cortical interactions are not impaired. In contrast, when Ach levels are very low, most cortical neurons undergo off-periods synchronously, producing a near-global low-frequency synchronization of neural activity. Under these conditions, the slow sleep waves merge into a single undifferentiated, synchronized network. Further analysis indicates that, during the falling asleep process, the variability of BOLD RSNs decreases, indicating a reduction in the repertoire of available transient brain states.

### Resting Functional Connectivity and Sleep

Although for practical purposes most investigations of sleep are performed using EEG, our explicit purpose here was to study how resting-state FC changes when the behavior of nodes in the cortical network is varied parametrically, simulating changes in neuromodulation during the falling asleep process. The slow dynamics of resting-state FC are traditionally studies using fMRI, and while some progress is being made in evaluating resting-state FC based on EEG, there are still many issues related to source modeling, to the nonstationarity of the signal, and to the presence of various nested rhythms, which make EEG analysis of FC less straightforward than fMRI analysis. Fortunately, recent work comparing fMRI FC with EEG and intracranial local field potentials shows that the envelope of the power spectrum of the signal in different frequency bands shows a very similar structure of slow correlations as fMRI FC (e.g. He et al. 2008). In addition, only by considering the full spatial structure of fMRI FC could we investigate the relationship between the synchronization patterns of slow sleep waves at different levels of neuromodulation with the anatomical connectivity of the cerebral cortex. Nevertheless, despite the focus on fMRI resting-state FC, the results and predictions derived from the present paper are in principle applicable also to high-density EEG data, as long as one can extract correlations in the EEG envelopes that show a slow dynamics.

A second point that bears emphasizing is that this paper does not attempt at modeling the sleep EEG. A realistic model of the resting EEG is no easy undertaking, due to complex coupling of frequency bands, the bistability of neurons, and the nonstationarity of events such as spindles, unless one resorts to detailed, phenomenological models as in our previous work (e.g. Hill and Tononi 2005). However, such detailed models cannot currently be used for addressing the questions we intended to address at the scale of a large number of connected cortical areas. As we showed before, the simplest way to capture the resting slow functional correlations evidenced by fMRI, which is what we have attempted to model, is by assuming an underlying asynchronous dynamics (Deco and Jirsa 2012). This is a minimal model, which assumes at each node a very simple asynchronous dynamics. Since this is not at all a trivial manner, we eschewed completely any attempt to model different frequency bands, and merely implemented ad hoc, as it were, the occurrence of slow sleep waves.

### The Role of Neuromodulators in the Transition From Wake to Sleep

Since the beginning of unit recordings in vivo, many studies have investigated changes in neuronal firing across behavioral states. Early reports found that, throughout the cerebral cortex, neurons often show a burst-pause pattern in NREM sleep, compared with tonic firing in wakefulness and REM sleep. Later, intracellular recordings in vivo showed that, during NREM sleep, virtually all cortical neurons alternate between a depolarized up state, during which they are spontaneously firing, and a hyperpolarized, silent down state (Steriade et al. 2001). This so-called slow oscillation (<1 Hz) is reflected in an alternation between periods of multiunit activity (on-periods) and silence (off-periods) in extracellular recordings and to negative peaks (slow waves) in the sleep EEG. In contrast, wakefulness and REM sleep are characterized by a sustained depolarization and tonic firing and are associated with an absence of EEG slow waves. The main factor controlling the transition from the tonic, irregular firing of wakefulness and REM sleep to the synchronous, on–off firing of deep NREM sleep is the level of arousal-promoting neuromodulators, chief among them Ach (Steriade et al. 2001). Ach does so by affecting both intrinsic and synaptic conductances (Gil et al. 1997; Marder and Thirumalai 2002) that influence both membrane polarization and the tendency of neurons to synchronize. In the present simulations, in the wake mode, high Ach levels dampen the magnitude of cortico-cortical excitatory recurrent interactions through presynaptic inhibition of glutamate release via muscarinic receptors (Hasselmo and Bower 1992). Moreover, high Ach levels dampen the tendency of neurons to synchronize by reducing spike frequency adaptation due to its effect on muscarinic M-currents (Hasselmo and Giocomo 2006). When Ach was progressively decreased, to mimick the falling asleep process, neuronal activity in the model changed gradually from the irregular, tonic mode of wakefulness to the synchronous, on–off pattern of NREM sleep. While the complex effects of many different neuromodulators (including histamine, norepinephrine, serotonin, and hypocretin) on cortical circuits were not captured in this model, the changes in firing patterns and dynamics were qualitatively not unlike those observed in multiunit recordings in both humans (Nir et al. 2011) and animals (Steriade et al. 2001), suggesting that they constitute a sufficient basis for inferring overall changes in FC. Nevertheless, extending the present approach to consider the effects of specific neuromodulators would be valuable. For instance, classical pharmacological studies suggest that noradrenergic agents are less effective than cholinergic drugs in chronically affecting EEG activation and the occurrence of slow waves (Domino et al. 1968; Berridge and Espana 2000). On the other hand, recent stimulation studies show that the noradrenergic system may be particularly effective in producing a rapid and powerful EEG activation (Carter et al. 2010; Constantinople and Bruno 2011). Serotonin's role in regulating slow sleep waves is even more complex, because it may stem indirectly from its ability to affect the release of other neurotransmitters including Ach, norepinephrine, and GABA. The end result may be biphasic, with an early enhancement of wake followed by an increase in NREM sleep (Ursin 2002; Watson et al. 2010).

### Local Sleep and Preserved Functional Connectivity

It was recently shown that, if animals are kept awake beyond their normal sleep time, populations of neurons in different cortical areas can suddenly go off-line in a way that resembles the off-periods of NREM sleep (Vyazovskiy et al. 2011). Strikingly, subsets of neurons could enter an off-period in one cortical area, for example motor cortex, but not in another. The occurrence of local off-periods in the motor cortex was associated with errors in a motor task. On the other hand, the overall EEG remained typical of an awake state and the animal appeared behaviorally awake with eyes open and was responsive to stimuli. When the animal fell asleep, in contrast, the EEG displayed typical slow sleep waves and the animal became behaviorally immobile and unresponsive, with eyes closed. Correspondingly, cortical neurons in multiple brain areas began to show frequent on–off oscillations in the slow wave frequency range. During the falling asleep process, evidence for local sleep and wake has also been obtained in humans. For example, the thalamus and the hippocampal formation may show signs of sleep before the cortex, and within the cortex different regions may begin to show clear-cut slow waves at different times (Rey et al. 2007). Even during full-fledged NREM sleep, multiunit and local field potential recordings from up to 10 regions have shown that off-periods and associated slow sleep waves are often local, occurring in just a few areas at a time. Only rarely, especially during deep sleep at the beginning of the night, are slow sleep waves truly global (Nir et al. 2011). The present simulations show that, consistent with experimental evidence, when Ach levels were mildly reduced, one or another cortical areas begun to show brief off-periods, while all the other areas continued the irregular, tonic activity typical of wake. Importantly, as long as the off-periods remained local, the overall organization of FC remained similar to that observed during wake (cf. Horovitz et al. 2008; Larson-Prior et al. 2009). This result suggests if a patch or area of the cortex goes off-line transiently and intermittently, it should not interfere substantially with long-range cortico-cortical interactions among the remaining areas. This is likely the consequence of the large number of alternative paths that can mediate cortico-cortical interactions, not to mention cortico-subcortico-cortical pathways not examined here. This finding is in line with modeling studies of the resilience of cortico-cortical FC to localized lesions (Alstott et al. 2009). It remains to be seen if off-periods localized preferentially to cortical areas acting as hubs (van den Heuvel and Sporns 2011) may have a greater effect on FC, although it should be kept in mind that off-periods, unlike lesions, are extremely brief and intermittent.

### Global Synchronization and the Reduction in the Repertoire of Resting-State Networks

When Ach levels were reduced substantially, most cortical areas in the model underwent repeated, frequent transitions between on- and off-periods, with a modal frequency around 1 Hz. These local changes were accompanied by a progressively more global synchronization of on- and off-periods among distant areas. The highly synchronous neuronal activity observed in the deepest sleep mode was associated with marked changes in FC, such that most cortical areas began to merge into a single undifferentiated, broadly synchronized network. To better characterize the changes in FC as a function of neuromodulation, we examined the variance of connectivity with Ach levels. The variance of connectivity characterizes the heterogeneity of the distribution of temporally coherent activity that the parcellated brain regions engage in the resting state. When the variance is low, most brain regions function in stereotypic synchrony, but when the variance is high, the wide range of correlations imply a large repertoire of possible interactions (states) among the brain regions. The results show that, in wake and early sleep, many different RSNs can be observed, whereas with the deepening of sleep the variance shrinks, with a steep decrease in variance at the transition to global BOLD synchronization. This finding is relevant in view of the possibility that the level of consciousness may be closely associated with the repertoire of neural states available to cortical networks (Tononi 2008). To the extent that the variance of RSNs evaluated with BOLD can be taken as reflecting the repertoire of neural states, the results of these simulations are consistent with the fact that consciousness is preserved at the beginning of the falling asleep process (stage N1), reports of mental content become more fleeting during light NREM sleep (stage N2), and may disappear altogether especially during deep NREM (N3) early in the night (Nir and Tononi 2010). The prediction of a decreased variance of RSNs with the deepening of NREM sleep could be evaluated experimentally in both sleep and anesthesia as well as in pathological conditions associated with the reduction in the level of consciousness such as vegetative and minimally conscious states. It should be noted that previous work in the area of nonlinear dynamical analysis of EEG data has also often found a reduction in various measures of EEG complexity with sleep, as well as with anesthesia, coma, and various brain disorders, using measures such as entropy, correlation dimension, Lyapunov exponents, generalized synchronization, and others (reviewed in Stam 2005). The present findings extend this approach to characterizing the variance of the repertoire of network states in the context of the analysis of resting-state networks modeled after fMRI data. Finally, it should be emphasized that studies using simultaneous EEG/fMRI techniques point to substantial changes in the coupling cortical areas when subjects enter deeper stages of sleep (N3) (Sämann et al. 2011; Boly et al. 2012). These studies suggest that, irrespective of changes in the repertoire of RSNs, sleep can also be accompanied by a loss of integration among specific cortical areas. Studies of effective connectivity using transcranial magnetic stimulation together with high-density EEG (Massimini et al. 2005; Ferrarelli et al. 2010) have also demonstrated that, in NREM sleep and anesthesia, there is both a reduction in the repertoire of brain states, as evidenced by the occurrence of global, stereotypic responses to stimuli, as well as in cortical integration, as shown by a breakdown of effective connectivity. In future work, the current model will be employed to assess changes in effective connectivity as a function of neuromodulation (cf. Esser et al. 2007), in addition to those in FC reported here.

In conclusion, we were able to demonstrate and predict: (1) that local slow sleep waves are structured macroscopically in networks that resemble the resting-state networks when they emerge during falling asleep; (2) in contrast, when the sleep level increases (neuromodulators decrease further to very low levels), slow waves become global and merge into a single undifferentiated, broadly synchronized network; (3) during the falling asleep process, the variability of RSNs decreases, indicating a reduction in the repertoire of available transient brain states and an obliteration of the RSNs.

## Funding

G.D. was supported by the ERC Advanced Grant: DYSTRUCTURE (no. 295129), by the Spanish Research Project SAF2010-16085, and by the CONSOLIDER-INGENIO 2010 Program CSD2007-00012 and the FP7-ICT BrainScales. G.T. was supported by the Paul Allen Family Foundation and by the McDonnell Foundation. P.H. was supported by Leenaards Foundation. The research reported herein was supported by the Brain Network Recovery Group through the James S. McDonnell Foundation.

## Notes

Conflict of Interest: None declared.

## References

Alstott
J
Breakspear
M
Hagmann
P
Cammoun
L
Sporns
O
Modeling the impact of lesions in the human brain
PLoS Comput Biol
,
2009
, vol.
5

6
pg.
e1000408

Bazhenov
M
Timofeev
I
M
Sejnowski
TJ
Model of thalamocortical slow-wave-sleep oscillations and transitions to activated states
J Neurosci
,
2002
, vol.
22
(pg.
8691
-
8704
)
Benda
J
Herz
A
A universal model for spike-frequency adaptation
Neural Comput
,
2003
, vol.
15
(pg.
2523
-
2564
)
Berridge
CW
Espana
RA
Synergistic sedative effects of noradrenergic alpha(1)- and beta-receptor blockade on forebrain electroencephalographic and behavioral indices
Neuroscience
,
2000
, vol.
99
(pg.
495
-
505
)
Biswal
B
Yetkin
F
Haughton
V
Hyde
J
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI
Magn Reson Med
,
1995
, vol.
34
(pg.
537
-
541
)
Boly
M
Perlbarg
V
Marrelec
G
Schabus
M
Laureys
S
Doyon
J
Pélégrini-Issac
M
Maquet
P
Benali
H
Hierarchical clustering of brain activity during human nonrapid eye movement sleep
,
2012
, vol.
109
(pg.
5856
-
5861
)
Brown
D
P
Muscarinic suppression of a novel voltage sensitive k+ current in a vertebrate neuron
Nature
,
1980
, vol.
183
(pg.
673
-
676
)
Brunel
N
Wang
XJ
Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition
J Comput Neurosci
,
2001
, vol.
11
(pg.
63
-
85
)
Bullmore
E
Sporns
O
Complex brain networks: graph theoretical analysis of structural and functional systems
Nat Rev Neurosci
,
2009
, vol.
10
(pg.
186
-
198
)
Cabral
J
Hugues
E
Sporns
O
Deco
G
Role of network oscillations in resting-state functional connectivity
Neuroimage
,
2011
, vol.
57
(pg.
130
-
139
)
Carter
ME
Yizhar
O
Chikahisa
S
Nguyen
H
A
Nishino
S
Deisseroth
K
de Lecea
L
Tuning arousal with optogenetic modulation of locus coeruleus neurons
Nat Neurosci
,
2010
, vol.
13
(pg.
1526
-
1533
)
Constantinople
CM
Bruno
RM
Effects and mechanisms of wakefulness on local cortical networks
Neuron
,
2011
, vol.
69
(pg.
1061
-
1068
)
de Pasquale
F
Della Penna
S
Snyder
AZ
Lewis
C
Mantini
D
Marzetti
L
Belardinelli
P
Ciancetta
L
Pizzella
V
Romani
GL
, et al.  .
Temporal dynamics of spontaneous MEG activity in brain networks
,
2010
, vol.
107

13
(pg.
6040
-
6045
)
de Pasquale
F
Della Penna
S
Snyder
AZ
Marzetti
L
Pizzella
V
Romani
GL
Corbetta
M
A cortical core for dynamic integration of functional networks in the resting human brain
Neuron
,
2012
, vol.
74

4
(pg.
753
-
764
)
Deco
G
Jirsa
V
Ongoing cortical activity at rest: criticality, multistability and ghost attractors
J Neurosci
,
2012
, vol.
32

10
(pg.
3366
-
3375
)
Deco
G
Jirsa
V
McIntosh
A
Emerging concepts for the dynamical organization of resting-state activity in the brain
Nat Rev Neurosci
,
2011
, vol.
12
(pg.
43
-
56
)
Deco
G
Jirsa
V
McIntosh
A
Sporns
O
Kötter
R
Key role of coupling, delay, and noise in resting brain fluctuations
,
2009
, vol.
106
(pg.
10302
-
10307
)
Domino
EF
Yamamoto
K
Dren
AT
Role of cholinergic mechanisms in states of wakefulness and sleep
Progr Brain Res
,
1968
, vol.
28
(pg.
113
-
133
)
Esser
SK
Hill
SL
Tononi
G
Sleep homeostasis and cortical synchronization: I. Modeling the effects of synaptic strength on sleep slow waves
Sleep
,
2007
, vol.
30
(pg.
1617
-
1630
)
Ferrarelli
F
Massimini
M
Sarasso
S
Casali
A
Riedner
BA
Angelini
G
Tononi
G
Pearce
R
Breakdown in cortical effective connectivity during midazolam-induced loss of consciousness
,
2010
, vol.
107
(pg.
2681
-
2686
)
Fox
M
Raichle
M
Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging
Nat Rev Neurosci
,
2007
, vol.
8
(pg.
700
-
711
)
Fox
M
Snyder
A
Vincent
J
Corbetta
M
Van Essen
D
Raichle
M
The human brain is intrinsically organized into dynamic, anticorrelated functional networks
,
2005
, vol.
102
(pg.
9673
-
9678
)
Fransson
P
Spontaneous low-frequency BOLD signal fluctuations: an fMRI investigation of the resting-state default mode of brain function hypothesis
Hum Brain Mapp
,
2005
, vol.
26
(pg.
15
-
29
)
Friston
K
Harrison
L
Penny
W
Dynamic causal modelling
NeuroImage
,
2003
, vol.
19
(pg.
1273
-
1302
)
Ghosh
A
Rho
Y
McIntosh
A
Kotter
R
Jirsa
V
Noise during rest enables the exploration of the brain's dynamic repertoire
PLoS Comput Biol
,
2008
, vol.
4
pg.
e1000196

Gil
Z
Connors
BW
Amitai
Y
Differential regulation of neocortical synapses by neuromodulators and activity
Neuron
,
1997
, vol.
19
(pg.
679
-
686
)
Greicius
M
Krasnow
B
Reiss
A
Menon
V
Functional connectivity in the resting brain: a network analysis of the default mode hypothesis
,
2003
, vol.
100
(pg.
253
-
258
)
Hagmann
P
Kurant
M
Gigandet
X
Thiran
P
Wedeen
VJ
Meuli
R
Thiran
JP
Mapping human whole-brain structural networks with diffusion MRI
PLoS One
,
2007
, vol.
2
pg.
e597

Hagmann
P
Cammoun
L
Gigandet
X
Meuli
R
Honey
C
Wedeen
VJ
Sporns
O
Mapping the structural core of human cerebral cortex
PLoS Biol
,
2008
, vol.
6
pg.
e159

Hasselmo
M
Bower
J
Cholinergic suppression specific to intrinsic not afferent fiber synapses in rat piriform (olfactory) cortex
J Neurophysiol
,
1992
, vol.
67

5
(pg.
1222
-
1229
)
Hasselmo
M
Giocomo
L
Cholinergic modulation of cortical function
J Mol Neurosci
,
2006
, vol.
30
(pg.
133
-
135
)
He
BJ
Snyder
AZ
Zempel
JM
Smyth
MD
Raichle
ME
Electrophysiological correlates of the brain's intrinsic large-scale functional architecture
,
2008
, vol.
105

41
(pg.
16039
-
16044
)
Hill
S
Tononi
G
Modeling sleep and wakefulness in the thalamocortical system
J Neurophysiol
,
2005
, vol.
93
(pg.
1671
-
1698
)
Honey
C
Kötter
R
Breakspear
M
Sporns
O
Network structure of cerebral cortex shapes functional connectivity on multiple time scales
,
2007
, vol.
104
(pg.
10240
-
10245
)
Honey
C
Sporns
O
Cammoun
L
Gigandet
X
Thiran
J
Meuli
R
Hagmann
P
Predicting human resting-state functional connectivity from structural connectivity
,
2009
, vol.
106
(pg.
2035
-
2040
)
Horovitz
SG
Fukunaga
M
de Zwart
JA
van Gelderen
P
Fulton
SC
Balkin
TJ
Duyn
JH
Low frequency BOLD fluctuations during resting wakefulness and light sleep: a simultaneous EEG-fMRI study
Hum Brain Mapp
,
2008
, vol.
29
(pg.
671
-
682
)
Jasper
H
Tessier
J
Acetylcholine liberation from cerebral cortex during paradoxical (REM) sleep
Science
,
1971
, vol.
172

3983
(pg.
601
-
602
)
Kötter
R
Online retrieval, processing, and visualization of primate connectivity data from the CoCoMac database
Neuroinformatics
,
2004
, vol.
2
(pg.
127
-
144
)
Larson-Prior
LJ
Zempel
JM
Nolan
TS
Prior
FW
Snyder
AZ
Raichle
ME
Cortical network functional connectivity in the descent to sleep
,
2009
, vol.
106
(pg.
4489
-
4494
)
Mantini
D
Della Penna
S
Marzetti
L
de Pasquale
F
Pizzella
V
Corbetta
M
Romani
GL
A signal-processing pipeline for magnetoencephalography resting-state networks
Brain Connect
,
2011
, vol.
1

1
(pg.
49
-
59
)
Marder
E
Thirumalai
V
Cellular, synaptic and network effects of neuromodulation
Neural Netw
,
2002
, vol.
15
(pg.
479
-
493
)
Massimini
M
Ferrarelli
F
Huber
R
Esser
SK
Singh
H
Tononi
G
Breakdown of cortical effective connectivity during sleep
Science
,
2005
, vol.
309
(pg.
2228
-
2232
)
Massimini
M
Huber
R
Ferrarelli
F
Tononi
G
Sleep slow oscillations as traveling waves: origins and pathways of propagation in humans
Sleep
,
2004
, vol.
27
(pg.
71
-
79
)
Mattia
M
Del Giudice
P
Finite-size dynamics of inhibitory and excitatory interacting spiking neurons
Phys Rev E
,
2004
, vol.
70
pg.
052903

Moruzzi
G
Magoun
HW
Brain stem reticular formation and activation of the EEG
Clin Neurol
,
1949
, vol.
1
(pg.
455
-
473
)
Murphy
MJ
Riedner
BA
Huber
R
Massimini
M
Ferrarelli
F
Tononi
G
Source modeling sleep slow waves
,
2009
, vol.
106
(pg.
1608
-
1613
)
Nir
Y
Staba
RJ
Andrillon
T
Vyazovskiy
V
Cirelli
C
Fried
I
Tononi
G
Regional slow waves and spindles in human sleep
Neuron
,
2011
, vol.
70
(pg.
153
-
169
)
Nir
Y
Tononi
G
Dreaming and the brain: from phenomenology to neurophysiology
Trends Cogn Sci
,
2010
, vol.
14

2
(pg.
88
-
100
)
Nobili
L
De Gennaro
L
Proserpio
P
Moroni
F
Sarasso
S
Pigorini
A
De Carli
F
Ferrara
M
Local aspects of sleep: observations from intracerebral recordings in humans
Prog Brain Res
,
2012
, vol.
199
(pg.
219
-
232
)
Nunez
P
Srinivasan
R
Electric fields of the brain: the neurophysics of EEG
,
2006
2nd ed
New York
Oxford University Press
Raichle
M
Mintun
M
Brain work and brain imaging
Ann Rev Neurosci
,
2006
, vol.
29
(pg.
449
-
476
)
Roberts
J
Robinson
P
Cortico-thalamic dynamics: structure of parameter space, spectra, instabilities, and reduced model
Phys Rev E
,
2012
, vol.
85

1
pg.
011910

Robinson
P
Phillips
A
Fulcher
B
Puckeridge
M
Roberts
J
Quantitative modelling of sleep dynamics
Philos Trans Ser A Math Phys Eng Sci
,
2011
, vol.
13369

1952
(pg.
3840
-
3854
)
Robinson
P
Rennie
C
Rowe
D
Dynamics of large-scale brain activity in normal arousal states and epileptic seizures
Phys Rev E
,
2002
, vol.
65

4
pg.
041924

Rogers
BP
Morgan
VL
Newton
AT
Gore
JC
Assessing functional connectivity in the human brain by fMRI
Magn Reson Imaging
,
2007
, vol.
25
(pg.
1347
-
1357
)
Rolls
E
Deco
G
The noisy brain
,
2010
Oxford (UK)
Oxford University Press
Sämann
P
Wehrle
R
Hoehn
D
Spoormaker
V
Peters
H
Tully
C
Holsboer
F
Czisch
M
Development of the brain's default mode network from wakefulness to slow wave sleep
Cereb Cortex
,
2011
, vol.
21
(pg.
2082
-
2093
)
S
de Pasquale
F
Mantini
D
Della Penna
S
A K-means multivariate approach for clustering independent components from magnetoencephalographic data
Neuroimage
,
2012
, vol.
62

3
(pg.
1912
-
1923
)
Stam
CJ
Nonlinear dynamical analysis of EEG and MEG: review of an emerging field
Clin Neurophysiol
,
2005
, vol.
116
(pg.
2266
-
2301
)
M
Jones
E
McCormick
D
Thalamus
,
1997
, vol.
1–2

Amsterdam
Elsevier
M
McCarley
RW
Brainstem control of wakefulness and sleep
,
1990
New York
Plenum Press
M
Timofeev
I
Greiner
F
Natural waking and sleep states: a view from inside neocortical neurons
J Neurophysiol
,
2001
, vol.
85
(pg.
1969
-
1985
)
Tagliazucchi
E
von Wegner
F
Morzelewski
A
Borisov
S
Jahnke
K
Laufs
H
Automatic sleep staging using fMRI functional connectivity data
Neuroimage
,
2012
, vol.
63

1
(pg.
63
-
72
)
Tononi
G
Consciousness as integrated information: a provisional manifesto
Biol Bull
,
2008
, vol.
215
(pg.
216
-
242
)
Rey
M
Bastuji
H
Garcia-Larrea
L
Guillemant
P
Mauguiere
F
Magnin
M
Human thalamic and cortical activities assessed by dimension of activation and spectral edge frequency during sleep wake cycles
Sleep
,
2007
, vol.
30
(pg.
907
-
912
)
Ursin
R
Serotonin and sleep
Sleep Med Rev
,
2002
, vol.
6
(pg.
55
-
69
)
van den Heuvel
M
Sporns
O
Rich-club organization of the human connectome
J Neurosci
,
2011
, vol.
31

44
(pg.
15775
-
15786
)
Vanini
G
Lydic
R
Baghdoyan
HA
GABA-to-ACh ratio in basal forebrain and cerebral cortex varies significantly during sleep
Sleep
,
2012
, vol.
35
(pg.
1325
-
1334
)
Vincent
J
Patel
G
Fox
M
Snyder
A
Baker
J
Van Essen
D
Zempel
J
Snyder
L
Corbetta
M
Raichle
M
Intrinsic functional architecture in the anaesthetized monkey brain
Nature
,
2007
, vol.
447
(pg.
83
-
86
)
Vyazovskiy
V
Olcese
U
Hanlon
E
Nir
Y
Cirelli
C
Tononi
G
Local sleep in awake rats
Nature
,
2011
, vol.
472
(pg.
443
-
447
)
Watson
CJ
Baghdoyan
HA
Lydic
R
Neuropharmacology of sleep and wakefulness
Sleep Med Clin
,
2010
, vol.
5
(pg.
513
-
528
)
Wedeen
VJ
Hagmann
P
Tseng
WY
Reese
TG
Weisskoff
RM
Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging
Magn Reson Med
,
2005
, vol.
54
(pg.
1377
-
1386
)