Emergent effects of synaptic connectivity on the dynamics of global and local slow waves in a large-scale thalamocortical network model of the human brain

Slow-wave sleep (SWS), characterized by slow oscillations (SOs, <1Hz) of alternating active and silent states in the thalamocortical network, is a primary brain state during Non-Rapid Eye Movement (NREM) sleep. In the last two decades, the traditional view of SWS as a global and uniform whole-brain state has been challenged by a growing body of evidence indicating that SO can be local and can coexist with wake-like activity. However, the mechanisms by which global and local SOs arise from micro-scale neuronal dynamics and network connectivity remain poorly understood. We developed a multi-scale, biophysically realistic human whole-brain thalamocortical network model capable of transitioning between the awake state and SWS, and we investigated the role of connectivity in the spatio-temporal dynamics of sleep SO. We found that the overall strength and a relative balance between long and short-range synaptic connections determined the network state. Importantly, for a range of synaptic strengths, the model demonstrated complex mixed SO states, where periods of synchronized global slow-wave activity were intermittent with the periods of asynchronous local slow-waves. An increase in the overall synaptic strength led to synchronized global SO, while a decrease in synaptic connectivity produced only local slow-waves that would not propagate beyond local areas. These results were compared to human data to validate probable models of biophysically realistic SO. The model producing mixed states provided the best match to the spatial coherence profile and the functional connectivity estimated from human subjects. These findings shed light on how the spatio-temporal properties of SO emerge from local and global cortical connectivity and provide a framework for further exploring the mechanisms and functions of SWS in health and disease.


Introduction
Internal states of the brain rapidly and profoundly influence sensation, cognition, emotion, and action [1][2][3].A dramatic change of internal brain state occurs at the transition between wakefulness and sleep [3][4][5].The main electrophysiological hallmarks of slow wave sleep (SWS) and rapid eye movement (REM) sleep are shared across mammals.During SWS, brain dynamics are dominated by oscillations between an active (Up) and a silent (Down) states, the activity called slow oscillation (SO, <1Hz) [6,7].Alterations of these sleep dynamics have significant consequences for brain function and are at the core of many neuropsychiatric disorders [8].
Until recently, sleep/wake transitions or transitions between different stages of sleep were thought to occur uniformly throughout the cortex.Recent evidence in rodents, primates and humans points to the existence of electrographic events resembling local sleep during behavioral wakefulness [9][10][11][12][13].During wakefulness, local groups of cortical neurons tend to fall silent for brief periods, as they do during SWS.Similarly, during behavioral sleep, transient local wake-like activity can intrude SWS [14,15].SWS and REM sleep can also coexist in different cortical areas [16][17][18].Homeostatic sleep pressure can affect cortical activity locally [19].Spindles and slow waves may increase locally after learning [20,21], suggesting experience dependent regulation.Together these results suggest an intriguing idea that local slow-waves may have a functional significance enabling brain to process and learn in a parallel way.However, the underlying mechanisms of local vs global slow-wave dynamics as well as factors controlling transitions between them are unknown.
We previously developed computer models of sleep SO [22,23], including transitions between awake and sleep states [24].However, the highly detailed and computationally intensive nature of these models have prevented up-scaling to investigate slow wave dynamics at the whole-brain scale.
Recently, several whole-brain models of slow-wave sleep incorporating biologically grounded connectivity from probabilistic diffusion MRI tractography and providing a good fit to human sleep recordings were proposed [25,26].These models investigated the role of longrange connections in SO brain-wide propagation.However, the mean-field neural-mass representation of cortical regions in these models provides only limited insights on how local (within a few milimeters) and long-range connectivity and activity interact to affect initiation and propagation of sleep slow waves.
To bridge this gap between micro and macro-scale mechanisms, we have developed a multi-scale, whole-brain, thalamocortical network model with biologically grounded cortical connectivity that exhibits the essential activity states of wake and slow wave sleep.The model consists of 10,242 cortical columns per hemishpere, each containing spiking pyramidal (PY) and inhibitory (IN) neurons arrayed in 6 layers, and 642 thalamic cells of each of four types: thalamo-cortical (TC) and reticular (RE), belonging to either core or matrix.
To investigate the role of connectivity in emerging SO, we systematically modified the density, range, and relative strength of cortical connections, as well as distance-dependent synaptic delays.The model dynamics were compared to the spatial coherence profiles obtained from human subjects.Our study revealed that a characteristic balance of long and short-range connectivity is necessary to match the data and enable the emergence of complex mixed states of local and global SO [18,27,28].It predicts that changes in synaptic weights over the course of sleep may enable transitions between global synchronized and local SO.

Connectivity
We start with a brief overview of the model (Fig 1, see Methods for details).The model includes 10,242 cortical columns positioned uniformly across the surface of one hemisphere, corresponding to the vertices in the ico5 cortical reconstruction reported in [29].The medial wall includes 870 of these vertices, so all analyses for this one-hemisphere model were done on the remaining 9372 columns.Each column has 6 cortical layers (L2, L3, L4, L5a, L5b and L6), with 1 pyramidal cell and 1 inhibitory cell per layer.Cortical neurons are modeled using computationally efficient map-based [30] spiking PY and IN neurons.Intra-columnar connections follow the canonical cortical circuit (Fig 1A ), while long-range cortical connectivity (Fig 1B) is based on diffusion MRI (dMRI) tractography [31] of the Human Connectome Project (HPC) young adult dataset [32], organized into the 180 parcels of the HCP multimodal atlas [33]).The originating and terminating layers of long-range connections are assigned according to the parcel's relative myelination-based [34,35] hierarchical positions (Fig 1D and 1E).Column-wise connectivity is obtained by applying the parcel-wise relation between connection length and probability to synthesized intercolumnar distances (Fig 1C, see Diffusion MRI guided connection probability for details).The resulting connectivity (Fig 1G) retains the parcel-wise structure from the data [31] (Fig 1F, Pearson's correlation r = 0.81, p < 0.0001), with an approximately exponential distribution between connection frequency and length (see Fig 1C).The connection delays, which are proportional to the geodesic connection lengths, follow a similar distribution.The thalamus is simulated with layers of thalamocortical (TC) and reticular thalamic (RE) neurons [23,36,37], and consists of the core and the matrix [38,39], which selectively project to the neurons from different cortical layers satisfying the outlined corticothalamo-cortical loop [40].Following [24] the model included the effects of neuromodulators to induce transitions between Wake state and SWS.
In the following sections, we first discuss "baseline" awake and sleep state dynamics of the model.Next we discuss how sleep dynamics are affected by altering network connectivity and reveal conditions for coexistence of local and global SO.We finally present analysis of human recordings and compare different model dynamics to human data using coherence analysis across frequency and distance and functional connectivity analysis.In addition to the basic characteristics of SO: frequency, amplitude, and propagation speed, network effects are Arrows indicate directed connections between layers.Neurons are never connected to themselves.B) Connectivity diagram between different cortical columns and thalamocortical cells in the core (TC) or matrix (TCa).C) Exponential decay of connection probability with distance between columns.Previously reported inter-parcel dMRI connectivity [31] shows an exponential relation with fiber distance.Column-wise connection probability is obtained by applying the dMRI exponential function to inter-column fiber distances (synthesized from their corresponding geodesic distances) and then adding the residual inter-parcel dMRI connectivity, estimated by regressing out the exponential trend (see Diffusion MRI guided connection probability for details).D) Estimated myelination index [34,35] throughout the cortex.A hierarchical index inversely proportional to this myelination index is assigned to each of the 180 cortical parcels designated in the HCP-MMP atlas [33].E) Excitatory corticocortical connections belong to one of 6 classes based on the myelination-derived hierarchical index of the pre-and post-synaptic neurons: within the same column, within the same parcel, weakly or strongly feedforward (from lower to higher hierarchical index) and weakly or strongly feedback (from higher to lower hierarchical index).Based on experimental reports on connectivity, weights are scaled according to the connection type by the factor shown in the matrices (see Hierarchically guided laminar connectivity for details).F) Parcel-wise connectivity as previously reported [31].Fpt stands for "fraction of probabilistic tractography" and represents the log 10 probability of connection resulting from fractionally scaling raw streamline counts (detailed in [31]).G) Generated model connectivity, which retains the parcel-wise structure from data in panel F (Pearson's correlation r = 0.81, p < 0.0001).Fpt in the model is computed as the log 10 of the ratio of the primarily described in terms of 3 key metrics: synchrony, spread, and participation.Synchrony refers to the timing of the Up state onsets and offsets across all neurons participating in a slow wave.Spread refers to how far the slow wave travels across the cortex.Finally, participation refers to the percent of neurons that participate in a slow wave within the area of spread.

Wake state characterization
Wake state (Fig 2A) was achieved by applying biologically relevant changes in key parameter values (Wake to sleep transition), based on the general trends established in our previous work [24].The resultant network shows a consistent average firing rate of individual neurons at about 17 Hz.Given the broad distribution of firing rates reported in awake state [41,42], we believe this model is appropriate.There is indeed experimental evidence of average firing rates in awake state being 15 Hz [28], 20-50 Hz [42], or 5-16 Hz [43].However, it should be noted that different mean firing rates can be achieved with minimal adjustments to key parameters, namely AMPA and GABA strength.The power spectral density (Fig 2a .5)revealed a distinct 1/f phenomenon, as is common for in vivo recordings in awake state [44].Individual neurons showed irregular tonic firing (Fig 2a .1),with baseline membrane voltage around -60 mV (Fig 2a .4).The firing activity was not synchronized (Fig 2a .3),which was reflected in the low voltage LFP (Fig 2a .2). Overall, the network revealed properties representative of biological awake state activity.

Slow wave characterization
To attain a network displaying SWS-like activity from the wake state baseline, several cellular and neurochemical parameters were adjusted according to the general methodology presented in [24].Namely, we increased cortical AMPA and GABA conductances as well as the leak current and slow hyperpolarizing current in all cortical neurons (see Methods, In general, the network spent more time in the Down states, and the Up state onsets and offsets were well synchronized (see narrow onset and offset histograms in Fig 2b .5).This pattern of longer Down states is typical for recordings under urethane anesthesia [6], while Down states are shorter in non-anesthetized animals [45].Nevertheless, the frequency of SOs was within the range of human data [46].The latency maps in Fig 2D demonstrate how waves of activity travel from a single initiation site, which varied between Up states with no detectable dependency on the previous ones.Each Up state involved >97% of all excitatory neurons.In this baseline model, slow waves revealed both high spread and high participation.
Cortical inhibitory, thalamocortical, and reticular thalamic neurons all showed similar patterns of oscillations (Fig B in S1 Text), as previously found in [23].While thalamic cells followed cortical oscillations, they were not directly involved in the generation of SO.In fact, number of connections from parcel A into parcel B to the total number of connections that either originate at A or terminate at B, excluding within-parcel connections.https://doi.org/10.1371/journal.pcbi.1012245.g001 .This is in line with experimental results showing that SO survives extensive thalamic lesions [6], and is present in neocortical slices [47] and isolated cortical slabs [22].Nevertheless, removing the thalamus reduced synchrony of cortical slow waves (see in S1 Text: Effect of thalamus on SO synchronization), in agreement with in vivo data [48].

Connection density
In the next two sections, we test effects of alternating the intracortical connection density and the maximal radius of connections on SO properties.Importantly, in these simulations we scaled the strengths of remaining synaptic connections to compensate for the lost connections and to ensure an equivalent amount of net activation per neuron (see Modifying network structure for details).
In the baseline model, structural connectivity was set to match biological data (see Methods).To test the effects of connection density between regions, we gradually changed the density of connections, by applying a probability P to retain each possible connection (Fig 3A and  3B) regardless of the distance between neurons, thus keeping the connection distance distribution unchanged (Fig 3C).Remaining synapses were scaled to keep total synaptic input per neuron unchanged.
Setting the probability from 100% down to 50% (i.e., ablating a random half of all connections) had minimal effect on the slow waves, indicating the model's robustness to even significant loss of individual connections as long as the total strength of excitatory synaptic input per neuron remains unchanged.Slow waves exhibited a similar propagation pattern with only a slightly reduced amplitude (approx.92% of baseline SOs) and speed (Fig 3D and 3F).The synchrony of the Up state offsets, however, was impaired; the distribution of offset timings across neurons became much wider (Fig 3d.Video), there was continued loss of synchrony in Up state offsets and onsets, slow waves became less frequent (down to 0.17 Hz), with lower amplitude (50% of baseline), and lower participation (only 70% of neurons involved per Up state) (Fig 3F).Nevertheless, slow waves still propagated through the all brain regions on the macro scale, showing that long-range synchrony is still possible (if slightly impaired) with only a fraction of long-range connections in place as long as those connections are sufficiently strong.Below this level of remaining connections, slow waves no longer occurred and the network was largely silent.

Connection range
To determine the specific contribution of long-range connections, we varied the maximum connection radius; any connections longer than a specified distance R were set to 0 (Fig 4A).Again, the strengths of remaining connections were scaled up to compensate for reduction in the total number of connections and to keep total synaptic input per cell unchanged (see Modifying network structure for details).In contrast to the density reduction (see Fig 3C ), this alteration did shift the connection distance distribution (Fig 4C).details).The percent of active neurons during each Up state is shown below the corresponding latency map.Up states involve nearly every neuron in the cortex within 300ms from its initiation time.
https://doi.org/10.1371/journal.pcbi.1012245.g002The maximum connection length in baseline model was 226.1 mm, but most connections (approximately 80%) were maintained within 100 mm (Fig 4B ).A connection radius as small as 5 mm, which includes second-order neighbors of a given neuron, resulted in only approximately 10% of connections still intact.The alterations of the maximum connection length led to significantly different activity patterns than in the network with a spatially uniform drop to the same percentage of connections (Connection density).The network with a maximum connection radius of 5 mm (Fig 4D ) showed no reduction in the frequency of the slow waves, but both the onsets and offsets of the slow waves revealed dramatic reduction in synchrony.Additionally, there was a substantial drop in slow wave propagation speed, while changes in connection density had minimal effect on speed (compare Fig 4F and Fig 3F).Despite all this, characteristic properties of global slow waves (high spread and high participation) found in the baseline model were maintained.
However, reducing radius below 5mm (i.e., not including 2nd degree neighbors) caused a dramatic drop in global oscillation amplitude, akin to a phase transition from global to local oscillations.With a connection radius of 2.5 mm (Fig 4E), including only first order neighbors of each cell, global slow waves were lost in favor of local slow waves (Fig 5, low spread, high participation).Note that this radius retains only 3.7% of connections, where a spatially uniform drop to less than 10% of connections removed slow waves entirely.The local slow waves were generated independently in few selected areas, and occasionally showed regional propagation but never travelled across the full hemisphere (Fig 5F).This resulted in reduced global SO amplitude as long-range synchrony between locally oscillating neural populations was impaired (Fig 5C and 5D).This model can further be seen in S3 Video.While a large fraction of the cortex was not participating in SO, the active regions (such as regions 3 or 6, Fig 5 ) were consistent across Up states.A connectivity analysis revealed that these regions correspond to the network areas with the highest node degree (see Fig E in S1 Text).The same regions exhibited local slow waves for different simulations with the same connectivity parameters but different random seeds for network generation, implying that these regions are not random but emerge from the underlying structure imposed by the dMRI connectivity data.
As a whole, this analysis revealed that long-range connections (longer than 5 mm) are essential for the large-scale synchrony of the slow waves, but having only local connections (less than 5 mm) is still sufficient to make slow waves propagating over entire cortex (as long as remaining connections strength is scaled up to maintain total synaptic input per cell).Even with only first order neighbors, the network still maintained slow wave initiation, synchrony, and participation with local propagation.

Connection delay
The model implements synaptic delays scaled by geodesic distance between neurons with a maximum delay of 40 ms.The delay for the majority of synapses is under 2 ms because the vast majority of connections are relatively short-range.To test the effects of delay times, we changed the delays in two ways ( Setting the maximum delay to be smaller (all the way down to 0.1 ms) had little to no effect on slow wave amplitude, frequency, or synchrony.Setting the maximum delay to be larger also had minimal effect on network activity until biologically implausible delays were implemented Summary of the effect of reducing connection density on the global SO frequency and amplitude as well as the standard deviation of the onset/offset delays (i.e. the width of the onset/offset histograms in d.5 / e.5).https://doi.org/10.1371/journal.pcbi.1012245.g003(�1 second).Even at this extreme the network retained strong global slow waves, albeit with slightly lower amplitude, longer frequency, and de-synchronized onsets and offsets.Setting a uniform delay of 10 ms (where most delays were previously under 2 ms) further left the network behavior largely unchanged.Only when the uniform delay was increased to 30 ms did we observe a loss of amplitude and synchrony (Fig G in S1 Text).Since the delay times here shown to have an effect on network behavior are clearly beyond the range of biological plausibility, these results indicate that synaptic delays are not a primary determinant of slow wave spatio-temporal properties.

Connection strength
To probe the effects of absolute connection strength on the network activity, we varied the base connection strength between all cortical pyramidal neurons (Fig 6A and 6C).In contrast to experiments with changing connectivity density or radius, where remaining synapses were scaled to keep total synaptic input constant, this changed the total synaptic input per cell.Weight reductions or increments by a factor of 4 or more caused the slow oscillations to be barely detectable or undetectable (note the low amplitude at the extremes of Fig 6c .1).This was a result of low activity in the case of weights reduction, and of sustained Up state (Down state absence) in the case of weights increase.In general, within the detectable range, the SO amplitude, frequency, and participation decreased when synaptic strength was decreased (Fig 6C).Note that the model revealed an overall robustness to weight variation: the slow-wave propagation was minimally altered by strength increases or decreases by a factor of 2.
To further investigate the role of long range connectivity, we next performed N times synaptic weights reductions on only "long-range connections", i.e., connections over a defined minimum radius R th (Fig 6B and 6D).We varied the scaling factor N from 2 to 7 and radius threshold R th from 2.5 mm to 50 mm.Fig 6B illustrates the case: N=5, R th =10mm.As mentioned previously, we did not scale the total synaptic input per cell, so these alterations affected balance of excitation and inhibition in the network.This alternation can be interpreted as modeling the effect of stronger synaptic inputs nearest to the cell body [49,50].
When connections over R th =10mm were significantly weakened (N=5) (Figs 6B and 7 In this model, the local populations of neurons would regularly synchronize producing local slow waves but the global synchrony was lost in many cases except for a subset of all events when global transitions between Up and Down states were observed (Fig 7B and 7E).This was further confirmed by running longer 150 sec simulation (see S4 Video) that revealed multiple transitions between epochs of local and global slow-waves.With reduced spread and high participation, the local slow wave behavior resembles the results from simulations with exclusively local connectivity (see Connection range), but stems from the more biologically plausible assumption of weaker (rather than nonexistent) long range connectivity.Indeed, the propagation patterns in

Local vs global slow-waves
Given above analysis, we next sought to quantify the global v.s.local slow waves in the model as a function of synaptic strength.For each model, we calculated the regularity in each of the 180 parcels, where regularity is defined as the inverse of approximate entropy.The mean and variance of this distribution are then plotted against each other, as shown in Fig 8A .We define global models as models with high regularity, local models as those with low regularity, and mixed states falling in between.The Intact Score of each model, defined as the percentage of remaining synaptic strength compared to the base model, is further predictive of the regularity and network behavior (Fig 8B).Two groups of models in the middle of distribution with Intact Score 30-50% and regularity 4-12, were identified as mixed states models.In Fig 8C, ten random parcel LFPs from representative global, mixed, and local models show demonstrably different behavior-the global model parcels are all in synchrony, while the local model parcels can be seen to transit between Down and Up states largely independently.The mixed model shows a combination of these behaviors, i.e., transitions to Up state sometimes occurred independently and sometimes in sync across parcels.This group predominately includes models with long-range connectivity scaled down beyond a radius of 5-10 mm.

Coherence analysis
To validate our model, we conducted a comparative analysis with experimental data.Specifically, we analyzed coherence across 180 cortical parcels.We calculated coherence in successive narrow-band signals across all parcel pairs to generate a snapshot of coherence across both frequency and distance.We then took the average of the coherence matrices in the SO range (<2 Hz), scaled the values by parcel-to-parcel distances, and masked the model data to match the sparsity of the experimental data.These modified 2D coherence matrices were then used to calculate 2 metrics of functional connectivity: Percent connected and number of communities.Percent connected reports the percentage of parcel pairs with a coherence greater than a given threshold, while the number of communities is determined using Louvain community detection.
Subsequently, we fit the full coherence matrices (across all distances and frequencies) with an exponential model.From this model, we estimated the spatial length constant (λ) and the scaling term (α).These four metrics were used to compare the functional connectivity and coherence profiles of different biophysical models against in vivo data.Full details of coherence analysis can be found in the Methods under Coherence analysis.
2.9.1 Empirical data.The experimental data collection is fully described in Experimental data.Briefly, continuous stereo-electroencephalography telemetry data were collected for 236 patients with focal epilepsy.After significant data cleaning to ensure removal of epileptic activity and subsequent sleep scoring, NREM sleep data were selected for analysis.In Fig 9, the meaningfully quantified.D) Summary of changes in SO characteristics when only long-range connections are reduced as shown in (B).Across all plots, 5 mm is seen to be an inflection point (or "elbow") where network activity changes.https://doi.org/10.1371/journal.pcbi.1012245.g006

Model analysis.
To perform a comparable model analysis, we first averaged all the cell voltages per parcel (all layers, excitatory and inhibitory neurons) and added noise with signal-to-noise ratio (SNR) = 5 to mimic the signal that would be recorded by a depth electrode, as in experimental data [54,55].We then performed coherence analysis in an identical manner to the experimental data.After the full coherence matrices were computed, they were further masked to fit the sparsity of the experimental data before metrics of functional connectivity were computed.
The functional connectivity metrics reported the total percent of connected parcels, and the number of communities as determined by Louvain community detection (Fig 9A).Global models were highly connected, showing a low number of distinct communities.Conversely, local models showed very sparse connectivity, with multiple distinct communities detected.Mixed state models showed moderate values on both scales.The experimental data point can be seen to fall in the middle of the distribution, surrounded by mixed state models.
We then compared the exponential fits of the full coherence landscapes in the SO range (less than 2 Hz) using the mean α and the first λ value (at 0.5 Hz); these two describe the peak coherence and rate of decay.When plotted against each other (Fig 9B), all the mixed state models revealed intermediate α and low λ values.The experimental data point can be seen to fall in close proximity to the mixed state models.In particular, models in which connections longer than 10 mm were scaled down by a factor of 5, 6, or 7 were consistently very close to the experimental data point; this is also true for the model where connections longer than 5 mm are scaled down by a factor of 4 (Fig 9A .2 and 9B.2).Importantly, we found that the global slow wave model, that was initially set as a baseline model, was very far off.Note that different metrics identified slightly different variations of the models as the best match to the data (compare identified regions in Figs 8 and 9); nevertheless, we found high overlap between different measurements, all indicating that the models generating mixed states are in good agreement with in vivo data.
Fig 9C -9E shows coherence profiles obtained for in vivo data, selected mixed states model (model with a greater than 10 mm reduction, reduced 5-fold) and baseline global state model.Both the experimental data and the model displaying mixed local/global slow waves exhibit a sharp peak of coherence at low frequencies, which quickly tapers off with increased distances between parcels (Fig 9C and 9D).Notably, this pattern is not observed in the global slow wave baseline model (Fig 9E ), where only a slight drop in coherence is seen as the distance between parcels increases.This is expected, given the unrealistic levels of slow wave synchrony across the entire network in this model.Additionally, the global slow wave model shows very high levels of coherence (nearing 1, indicating perfect coherence) at low frequencies.
In summary, these results suggest that mixed states models, particularly those with a 5-10 mm radius of stronger local connections, are the best match for experimental data.

Discussion
Traditionally, brain states have been categorized into wakefulness, slow-wave sleep (SWS), and rapid eye movement (REM) sleep, assuming their global occurrence across brain regions.However, recent large-scale neural recordings have revealed the rich heterogeneity of neural dynamics on a brain-wide scale, demonstrating that sleep and wake signatures can coexist in different brain regions, influencing behavior with spatial and temporal specificity.While local circuit mechanisms of slow waves-alternations of positive/negative EEG/LFP waves, also called slow oscillation (SO)-have been well described, the factors leading to complex heterogeneous brain states with local slow waves remain unknown.Here, we present a wholebrain thalamocortical network model based on human dMRI connectivity data capable of generating awake-like state and sleep states with biologically realistic slow waves.The model was applied to make specific predictions about the effects of long-and short-range connectivity on the spatio-temporal dynamics of the slow waves.We identified intracortical synaptic connectivity properties that allow the coexistence of local and global slow waves and proposed hypothetical mechanisms for transitions between local sleep and global uniform sleep states.The model predictions were compared to human data to validate and fine-tune the model.
A prominent example of brain states heterogeneity is local sleep [9][10][11][12][13].During wakefulness, local groups of cortical neurons tend to fall silent for brief periods, as they do during SWS.Furthermore, during behavioral sleep, transient local wake-like activity can intrude SWS [14,15].It is worth noting that SWS and REM sleep can also coexist in different cortical areas [16][17][18].Additionally, the proportion of time allocated to wakefulness, SWS, and REM states varies significantly across cortical areas [16,56], highlighting the diverse predisposition of different regions to support these states.
From a mechanistic perspective, the basic need for local sleep can possibly be explained by homeostatic sleep pressure, that accumulates throughout the day, shaped by the duration and quality of preceding alertness.This pressure is evident in the intensity of SWS in the cortex, indicating the brain's need for recuperation [57].Homeostatic sleep pressure can influence cortical functions in specific areas [19].For example, spindles and slow waves may intensify in regions engaged in learning [20,21], suggesting that sleep pressure is regulated based on recent cognitive activities.Notably, after motor skill training, there can be a localized increase in SWS within the motor cortex during rest [58][59][60], and a visual perception learning task can lead to a rise in the initiation of slow waves in the lateral occipital cortex [61].
Why brain areas that are more actively involved in sensory processing during the day require more SWS is a challenging question.We can speculate that it may be related to the well-characterized function of SWS in memory consolidation [62].The slow oscillation (SO) repeatedly resets the thalamocortical network during the Down phase and temporally groups sleep spindles, which are associated with long-term potentiation (LTP)-like processes and learning, during the Up phase [63,64].These observations have led to the hypothesis that SO provides a global temporal frame within the cortex and between brain regions (e.g., hippocampus, thalamus, striatum) for offline memory processing and reactivation through the strengthening of neuronal circuits [2,[65][66][67][68].Local sleep may thus enable cortical networks to augment memory consolidation in specific brain regions that need it the most.
While specific mechanisms of slow-wave sleep (SWS) generation have been explored in many local circuit models [22,23,[69][70][71][72][73], only a few very recent studies [25,26] have proposed simulating realistic long-range cortical connectivity to test its effect on the spatio-temporal properties of SO.While these models take an important step in utilizing global brain connectivity data to study SWS, they apply a mean-field neural-mass representation of cortical regions, which limits their ability to study how the interaction of local (within a few millimeters) and long-range connectivity affects spatio-temporal SWS dynamics.
Here, we developed a large-scale thalamocortical network model based on spiking excitatory and inhibitory neurons in thalamus and cortex, organized in 6 cortical layers and thalamic core and matrix subsystems, with connectivity based on Human Connectome Project (HCP) diffusion-MRI (dMRI) tractography-based connection probabilities and HCP-derived degree of grey-matter myelination [31,32].To test the contribution of specific elements of connectivity to SO properties, we then systematically modified selected elements of the network connectivity: connection density and radius, synaptic delays and strength (Table 1).
While we inferred connectivity profiles from dMRI data, predicting the strength of synaptic connections is much harder problem, especially considering that our model is still a significant scaling down of the biological brain.Therefore, a substantial effort of this project was to explore the strength of connections to identify the SWS regimes that match those known from the literature and characterized in our data.We found that altering the balance of excitation and inhibition had a significant impact on the synchronization properties of slow waves.While our baseline model produced unrealistically synchronized global slow-waves, reducing the strength of "long-range" excitatory connectivity resulted in models with mixed global and local slow waves, similar to results reported in vivo [17,27].This occurred in models where connections longer than 5-10 mm were reduced four-to five-fold compared to the baseline model, to mimic the high density and strength of local connections in vivo [49,50].Interestingly, this aligns with the results presented by [74], which show an approximately 7mm radius within which slow waves are highly synchronous.Furthermore, [53] and [52] have reported that the occurrence of multi-peak slow waves, where each peak has a different origin, increases over the course of sleep, which was proposed to result from a gradual net decrease in the strength of cortico-cortical connections [28] (however, see [75]).Human recordings [18,51] reported two types of slow waves: Type I-global and more frequently occurring earlier in sleep, and Type II-local slow waves found later in sleep.A recent study suggests that global slow waves (referred to as SO) are responsible for the consolidation of new memories, while local events (referred to as delta oscillations) lead to forgetting [76].
To validate the model, we compared models with different synaptic connectivity profiles to human recordings during NREM sleep.This was done via a coherence analysis to determine correspondence of signals across distance and frequency and functional connectivity analysis to evaluate percentage of parcel pairs with a coherence greater than a given threshold and the number of communities.We found the experimental data to closely match the mixed state models, i.e., models showing coexistence of local and global slow-waves.Importantly, we also found that relatively modest changes of synaptic connectivity in these models lead to the models with only local or global slow-wave activity.
Taking these results together, our model predicts that changes of the cortical synapses strength over night, may explain more frequent occurrence of global slow waves earlier in sleep.However, a more extensive analysis of synaptic connectivity in the model and analysis of in vivo data are needed to confirm this prediction.
Several studies have revealed a link between long-range intracortical connectivity and the properties of sleep slow waves.In [77], it was shown that sleep slow oscillations propagate over longer distances with increasing age, and cortical SO propagation was positively correlated with intrahemispheric myelin content.Additionally, it was found that human subjects with a steeper rising slope of the slow wave exhibited higher axial diffusivity in the temporal fascicle and frontally located white matter tracts [78].These findings suggest that the profiles of sleep oscillations reflect not only the synaptic-level dynamics of the neuronal network but also the microstructural properties of its structural foundation, the white matter tracts.Factors such as normal aging [79] and Alzheimer's Disease (AD) [80] impact long-range connectivity.Specifically, in patients with severe Alzheimer's Disease, functional connectivity was notably reduced between regions separated by greater distances, and this loss of long-distance connectivity was associated with a less efficient global and nodal network topology [80].Alongside well-documented observations of reduced slow wave density in aging [81] and Alzheimer's Disease [82], these findings align with our model predictions about the critical role of balancing local and long-range connectivity for maintaining characteristic slow-wave dynamics.
The sleep SO is hypothesized to be a cortical rhythm [6,22,47,83].In line with these results, cortical SO were maintained in the model even when thalamic input was removed.However, other studies have shown that the normal pattern of the neocortical SO requires thalamic inputs to synchronize activity across large cortical areas [48].Here, we found that removing the thalamus reduced synchrony of slow-waves in the mixed state (but not the global) SO model, suggesting its role in long-range synchronization of cortical activity.A deeper study into thalamus involvement in SO synchrony will be a highly relevant avenue for future experimental and modeling studies.Moreover, despite not driving the SO rhythm that constitutes the object of this study, the thalamus allows this model to capture a range of human brain dynamics beyond SO.Indeed, natural extensions of our model include investigating sleep spindles, in which thalamus plays the major role [7,84], or impact of thalamic stimulation, which constitute an active field of neuroscience research [85,86].
Although our model design achieves high computational efficiency considering the size and complexity of the network, the number of the model parameters and simulation times do not easily allow a systematic fitting of the model output to electrophysiological brain signals, as it has been done with smaller computational models [25].It is possible that optimization techniques such as surrogate model optimization or simulation based inference may be used on a subset of parameters in the model with the addition of biological constraints on possible values.Future studies may wish to explore this, but care will have to be taken to maintain biological plausibility.Nevertheless, the use of realistic parameters for individual neuron dynamics in wake and sleep stages [23,24] as well as experimentally grounded cortical connectivity [31] allowed us to reproduce the essential characteristics of SWS activity in human brain.
Importantly, we have now established a biologically relevant model of network dynamics that can be used for a multitude of other applications.For example, this model could be used to study possible mechanisms of sleep disorders by testing the fit of various hypothesis-driven perturbations of the model to relevant experimental data.In that context, one particularly important application of the model is using it to reveal mechanisms of sleep alterations in aging and Alzheimer's Disease.Beyond the realm of sleep, this model could be used to study the response of a healthy brain to various types of traumatic brain injury or neurodegenerative diseases that may alter neural connectivity.Finally, this model could further provide preliminary insight on the network effects of drugs that act on specific neurotransmitters, either globally or in a targeted location.
In summary, we present a whole-hemisphere thalamocortical network model constrained by realistic local and long range connectivity based on human dMRI data and implementing effects of neuromodulators to account for sleep-wake transitions.The model displays both a biologically realistic awake state with characteristic random neuron firing as well as Type I and Type II slow wave sleep.The critical results of our study lie in revealing the separate contributions of structural connectivity (i.e., connectivity matrix) and functional connectivity (i.e., synaptic connectivity strength): while structural connectivity can alter the characteristics of individual slow waves (synchrony, spread, participation), it is primarily the proper balance of connectivity strength between local and global connections that leads to the mixed cortical states, where periods of synchronized global SO are intermittent with the periods of asynchronous local slow-waves.Importantly, only the models generating such intermittent states show spatio-temporal SO characteristics that match human recordings.Comparing the models to empirical data further supported the fundamental role of connectivity strength in generating biophysically realistic network behavior via mixed-state models.Furthermore, the model provides both the scale and the level of detail required to investigate how large-scale sleep dynamics arise from cellular and circuit level activity.

Ethics statement
Clinical data were recorded as standard of care and written informed consent, approved by each institution's review board (The Institutional Review Board of the University of California, San Diego, The Oregon Health and Sciences University Institutional Review Board, The University of Pennsylvania Institutional Review Board, The University of Alabama, Birmingham Institutional Review Board for Human Use, The Institutional Review Board at the Cleveland Clinic Foundation, The Mass General Brigham Institutional Review Boards), was obtained from each patient before the data were used for research purposes.

Map-based cortical neurons.
The model has 10,242 cortical columns uniformly positioned across the surface of one hemisphere, one for each vertex in the ico5 cortical reconstruction previously reported [29].The medial wall includes 870 of these vertices, so all analyses for this one-hemisphere model were done on the remaining 9372 columns.Each column has 6 cortical layers (L2, L3, L4, L5a, L5b and L6).For each layer in each column, one excitatory (PY) and one inhibitory (IN) neuron were simulated.To allow for large-scale network simulations we modeled cortical neurons with a model based on difference equations (MAP) [87,88] which has a number of distinct numerical advantages: the common problem of selecting the proper integration scheme is avoided since the model is already written in the form needed for computer simulations.The model parameters can be adjusted to match experimental data [88][89][90].Map-based models sample membrane potential using a large discrete time step compared to the widely used conductance-based models described by ordinary differential equations, and still capture neural responses from these models.Importantly, map-based models replicate spiking activity of cortical fast-spiking, regular-spiking and intrinsically bursting neurons [88].Variations of the map model have enabled the simulation of large-scale brain networks with emergent oscillatory activity [90], including realistic sleep spindles [29] and cortical slow oscillations [30].In our study, we used the previously proposed map model modification [30], which implemented nonlinear dynamical bias, activity dependent depolarizing mechanisms, and slow hyperpolarizing mechanisms capturing the effects of leak currents, the Ca 2+ dependent K + currents and the persistent Na + currents, which were found to be critical for simulating up/down state transitions during SO.Parameter values were initially set to the previously reported reference values [30], which were shown to produce biologically realistic neural behavior during SO, and then tuned with small variations to maintain adequate neural responses in our comparatively larger-scale network.The map model equations are described below, and all parameter values specific to this study are reported in Table 2.
The activity of individual pyramidal PY cells was described in terms of four continuous variables sampled at discrete moments of time n : x n , dictating the trans-membrane potential V n = 50x n − 15; y n , representing slow ion channel dynamics; u n , describing slow hyperpolarizing currents; and k n , determining the neuron sensitivity to inputs.The activity of IN cells, which biologically have faster spiking dynamics, was modeled using only variable x n and fixing the other variables to constant values.Dynamical PY variables evolve according to the following system of difference equations, where n = 1, 2, . . .indexes time steps of 0.5ms in size (as suggested in [90]): spike generation Table 2. Map-based model parameter definitions and values.Parameters γ and δ may take one of the two values shown (see text).Parameters were initially set to the values in [30] and then fine tuned to maintain biologically reasonable spiking behavior in the current larger-scale network.Sub-index X stands for "AMPA" or "GABA".For all parameters regarding synapses, the PY and IN columns denote the post-synaptic neuron.where function H(w) denotes the Heaviside step function, which takes the value 0 for w � 0 and 1 otherwise.The reduced model for IN cells is given by

Name
where y* is a fixed value of y n .We describe each of the system's components below.The model's parameters and their values are summarized in Table 2.
Spike-generating Systems.(Eqs 1 and 4): Variables x and y model the fast neuron dynamics, representing the effects of the fast Na + and K + currents responsible for spike generation [91].The nonlinear function f PY , modified by [30] from [87], determines the spike waveform as where α and w 0 are constant parameters that control the non-linearity of the spike-generating mechanism.Parameters μ 1 , μ 2 < 1 modulate the slower updating of variable y with respect to variable x, and parameter σ dictates the neuron resting potential as (1 − σ).The scaling constants β x , β y modulate the input current I total for variables x and y, respectively.The neuron's membrane potential is calculated as V n = 50x n − 15 and a spike is generated if and only if V n > 0.1.In simulations, the quantity β x I total (x n , u n ) was manually prevented from reaching 0 by imposing a lower bound of 10 −4 .For inhibitory neurons IN, function f IN in Eq 4 is given by A spike is generated if and only if 0 Adaptation Eq 2. The slow variable u n represents the effects of the slow hyperpolarizing potassium currents (I D , explained below) reducing neural excitability over the course of the Up state and involved in Up state termination.The value of γ is taken as 0.99 for x n < −1 and 0.995 otherwise.Parameter δ = 0.24 has an effect only if −0.5 < x n < 1; we assume δ = 0 for other values of x n .
Sensitivity Eq 3. The variable k n reduces neuron sensitivity to inputs during Up states (i.e. during high spiking activity) [45] by regulating the spike nonlinearity function f PY .This gives a net phenomenological representation of the combined effects of high-level synaptic activity, increase in voltage-gated hyperpolarizing currents, and/or adaptation of fast Na + channels.It adopts the value K 0 (low sensitivity) when voltage x n crosses threshold -0.5 (i.e.Up state) and switches to K 1 > K 0 (high sensitivity) when x n < −1 (i.e.cell leaves Up state).This prevents over-excitation and provides a physiological firing frequency [30].
Computation of currents.The total current I total is a sum of external applied DC current (I DC ), synaptic inputs (I syn ), persistent Na + current (I Nap ), slow hyperpolarizing currents (I D ), and leak currents (I leak ): with the last three dictated primarily by p Nap (maximal conductance of persistent Na + current), p D and p L (strength of leak current), respectively: The persistent Na + current contributes to initiating and maintaining Up states [22].It is modeled in I Nap through a sigmoidal activation function depolarizing the neuron with strength ρ Nap at high voltage values.The effect of slow hyperpolarizing currents contributing to Up state termination is represented by I D , which is regulated by the dynamic variable u n (described above).Lastly, ρ L represents the effect of potassium leak current, and here it is used to model the hyperpolarization of cortical neurons during sleep [23].For each neuron i, the synaptic current I i syn is the sum of individual AMPA, GABA A and NMDA synaptic currents.The individual synaptic current of type X from neuron j to neuron i has the general form for AMPA and GABA synapses, and for NMDA synapses, with an additional modulator φ n , Here g ij X ðnÞ represents conductance at time n for directed connections from neuron j to neuron i. Dynamic synaptic conductances were described by the first-order activation schemes: |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } for AMPA and GABA A and for NMDA type synapses.Here, g g X determines the decay rate of synaptic conductance and G syn X regulates the overall synaptic strength.To keep total synaptic input per cell constant, coefficients G syn X and G mini X are scaled by the number m i X of synapses of type X incoming to neuron i.The depression variable d X is dependent on whether the pre-synaptic neuron j spikes at time n or not, that is, if the voltage of neuron j surpasses a threshold V θ .Lastly, miniature post-synaptic potentials are only present in AMPA and GABA synapses.Variable Z mini n is a spontaneous event from a random Poisson process with mean frequency μ X , and G mini X is the constant by which conductance increases in an event (Z mini X ¼ 1), which produces miniature post-synaptic potentials ("minis") in neuron i [92].Following [23], the frequency μ X is modulated in time to account for the reduction in mini frequency after Up states [22], so that where n j 0 is the time of the most recent pre-synaptic spike.This modulation reduces the mini rate to 20% of its maximal value � m X after a pre-synaptic spike, and then allows recovery to � m X with time.

Hodgkin-Huxley thalamic neurons.
The thalamus was modeled using a network of core (specific) and matrix (non-specific) nuclei, each consisting of thalamic relay (TC) and reticular (RE) neurons.We simulated 642 thalamic cells of each of these four types.TC and RE cells were modeled as single compartments with membrane potential V governed by the Hodgkin-Huxley kinetics, so that the total membrane current per unit area is given by where C m is the membrane capacitance, g L is the non-specific (mixed Na + and Cl − ) leakage conductance, E L is the reversal potential.These parameters (Table 3) were unchanged from [36], where they were determined based on previous modeling studies grounded in experimental observations [93,94].Specifically, thalamic reticular and relay neurons included intrinsic properties necessary for generating rebound responses which were found to be critical for spindle generation [95,96].Furthermore, I int is a sum of active intrinsic currents, and I syn is a sum of synaptic currents.Intrinsic currents for both RE and TC cells included a fast sodium current, I Na , a fast potassium current, I K , a low-threshold Ca 2+ current, I T , and a potassium leak current, I KL = g KL (V − E KL ).For TC cells, an additional hyperpolarization-activated cation current, I h , was also included.The expressions and parameter values for each of these currents are given in [23,29,36,97].Synaptic currents (GABA, NMDA and AMPA) were modeled by first-order activation schemes [94], using the general form where g syn is the maximal conductance, [O](t) is the time-dependent fraction of open channels, and E syn is the reversal potential.Our previous studies show specific equations for all synaptic currents [23,98], as well as a detailed description of O(t) based on first-order activation schemes [22,37].

Wake to sleep transition.
The transition from wake to slow wave sleep stage was accomplished via tuning of several key parameters.The primary tunable parameters of the map model cortical neurons are p D and p L , which control the strength of the slow hyperpolarizing currents I D (similar to a calcium current) and the potassium leak currents I leak as described in Eqs 10 and 11.Both of these parameters significantly increased from wake to sleep.We further increased the strength of AMPA and GABA synaptic currents, as well as the strength of TC potassium leak currents g KL , similar to the changes described in [24] (see Table 4).

Network structure
Our whole hemisphere computational model has realistic local and long-range cortical connectivity.Connections within the column follow those described in the canonical cortical circuit (Fig 1A).Additionally, connections within 0.1 mm are strengthened by a factor of 5 to mimic in-vivo increased strength of proximal synapses [49,50].Inhibition is local: IN cells project only to PY cells in the same layer, with connections only in their own column, 1st and 2nd degree neighbors.The strength of inhibitory synapses is 2 times larger within the local column than outside.Finally, miniature postsynaptic potentials on intracortical AMPA synapses are modeled as a Poisson process (as in [97]).
Thalamocortical (TC) and reticular (RE) neurons are modeled using Hodgkin-Huxley dynamics [91], with 2 subtypes of each to represent matrix and core neurons.The matrix TC neurons have a 45 mm (and 80 mm) fanout radius to cortical excitatory (and inhibitory) neurons while core TC neurons have 12 mm and 13 mm radii, respectively.Within-thalamic connections have 11 mm radii.There are 642 neurons of each of these 4 types (thalamic matrix, thalamic core, reticular matrix, reticular core).

Hierarchically guided laminar connectivity.
Long-range connectivity in the model is primarily based on the cortical parcellation into 180 areas proposed by the Human Connectome Project (HCP)-MMP1.0 atlas [32,33].Each parcel was assigned a hierarchical index inversely proportional to its estimated bulk myelination [34,35].Excitatory cortical connections (PY-PY) were then split into 6 classes according to the pre-and post-synaptic hierarchical index: within the local column, within the local parcel (or between contralateral homologs), strongly (>50th percentile) or weakly (�50th percentile) feedforward (from lower to higher hierarchical index), and strongly or weakly (>, � 50th percentile respectively) feedback (from higher to lower hierarchical index).Based on previous connectivity reports [99,100], synaptic weights for these connections were scaled by the factors shown in Fig 1E .Note that, while these factor impose constraints on the relative inter-layer connection strengths, the absolute value of synaptic strength is a free parameter.4.3.2Diffusion MRI guided connection probability.The connection probability was further scaled based on previous diffusion MRI (dMRI) connectivity studies [31].These observations reveal an exponential decay relationship between inter-parcel dMRI connectivity and fiber distance (length constant λ = 23.4mm,scaling parameter β = 0.17).In order to obtain a probability of connection at the scale of columns, approximate intercolumnar fiber distances were estimated from geodesic distances.For parcel centroids, intra-hemispheric dMRI streamline distances and their corresponding geodesic distances were found to be related to a first approximation by a linear rational function with p1 = 295.6,p2 = −2256, and q1 = 69.4,where x denotes the geodesic distance.Columnwise fiber distances were estimated by applying this function to intercolumnar geodesic distances.Thus, a distance-dependent probability of connection between columns was obtained by applying the dMRI exponential relation to column-wise fiber distances (Fig 1C).Moreover, the residual parcel-wise dMRI connectivity (not distance related, obtained by regressing out the exponential trend) was added back to intercolumnar connections based on the columns' parcels.This distance-independent connectivity accounts for the functional specialization of each parcel.Lastly, conduction delays were set to be proportional to fiber tract distances, with the longest connection (*226mm) having an assigned delay of 40ms.

Analysis
4.4.1 Modifying network structure.We here attempted to create the most biologically realistic model and connectivity possible.However, there is no ground truth to dictate the absolute strength of connection between two individual neurons.While rules exist in our architecture to modify relative strength based on biological connectivity, the base connection strength to be modified is a completely free (and incredibly important) parameter.We first set the base strength level to be high enough to generate Type I slow waves which quickly synchronize across the entire brain.This makes slow wave property analysis more straightforward and less prone to spatial bias when modifying further network parameters, as is the goal of this paper.Hence, we use this global Type I slow wave model for all subsequent analysis aside from analyzing the effect of connection strength itself.
In the case of modulating connection density and range, the synaptic strengths were scaled up proportionally to the loss of connections to retain sufficient activity levels in the network.The scaling mechanism works on an individual neuron basis.For each neuron, the weights of all synapses of each type (AMPA, GABA, etc.) are scaled by the total number of synapses of that type.Thus, after removing connections, the weights are divided by a smaller quantity, which results in stronger connections.This approach allows us to test the relative effects of connectivity, e.g., local vs global, without dramatic changes in the network dynamics.This allowed us to isolate the effects of patterns of connection apart from changes in total synaptic input.The ablation of connections in these trials was accomplished by setting synaptic connection strengths to 0, removing all functional connectivity.Layer LFPs were used to detect cortical Up states.A peak detection algorithm (MATLAB's islocalmax function) was implemented for this purpose, requiring a minimum peak separation of 500ms and a minimum prominence of half the difference between the maximum and minimum values of the signal.The amplitude of each Up state was defined as the difference between the corresponding peak and the lowest point before the next peak.The inverse of each inter-peak interval was used to approximate the oscillation frequency.Average amplitudes and frequencies across all Up states of a simulation are reported in the summary plots of Figs 3, 4, 6 and G in S1 Text.
4.4.3Onset/Offset detection.In order to characterize how SO spread through the cortex, we identified the points in time when each neuron transitions from Down state to Up state (onset) and viceversa (offset).Our Up state detection algorithm follows [30], based on [101]: we set two thresholds, V + = −65mV and V − = −68mV (see dashed vertical lines in voltage histogram of Fig 2B ), and label any activity above V + as an Up state and any activity below V − as a Down state.This initial detection was further refined by merging any two Up states (or Down states) that were less than 100ms apart, and then removing any remaining Up state (or Down state) lasting less than 100ms, in a similar approach to [101].
To identify onsets and offsets for a particular global Up state, we considered the time window between consecutive minima of the average voltage signal around that Up state (e.g.see black upward triangles in average voltage signal of Fig 2B).For each neuron, we found the first time (if any) in this window when it reached an Up state (onset), and then the first time after the onset when it reached a Down state (offset).The minimum onset was then subtracted from all onsets so that the earliest neuron would have an assigned time of 0 and all other would have positive delays, or latencies.Similarly, the earliest offset was subtracted from all offsets.
In order to evaluate the synchrony of active and silent states, we constructed aggregated histograms of the onset and offset times across all Up states of each simulation.We used the standard deviation of the histograms as an measure for synchrony, with narrow histograms indicating highly synchronized global Up states.The percent of participating columns in each Up state was defined as the number of columns with onset delays greater than or equal to 0 over the total number of columns not belonging to the medial wall.Non-participating columns are shown in white in the latency maps for each Up state.The average percent of active columns per Up state is reported in the summary plots in Figs 3F, 3F, 6C, 6D and G in S1 Text.

Speed of propagation.
To obtain speed of propagation of each Up state, we first identified columns with onset delays under 10ms to get an approximate area of origin.We then selected the column with minimal sum of (geodesic) distances to all columns with onset delays under 10ms, excluding the medial wall.This column was identified as the origin of the Up state.We considered only an area of 150mm radius around the origin to compute propagation speed, in order to avoid interference from waves originating later in other cortical regions.A linear regression was performed between the onset delays of columns in the area and their distances to the origin.For regressions with an R 2 > .25, the speed of propagation was defined as the slope of the regression line.The average propagation speed over all Up states is reported in summary plots.Speed is not reported for simulations in which there is not a linear relation between onset delay and distance from the origin, that is, simulations in which regressions for all Up states yield R 2 � .25.

Experimental data
4.5.1 Patients and intracranial recordings.Continuous stereo-electroencephalography telemetry and structural imaging were obtained for 298 patients undergoing pre-resection intracranial monitoring for phamaco-resistant focal epilepsy at University of California, San Diego Medical Center, La Jolla CA, the Cleveland Clinic, Cleveland OH, Oregon Health and Science University Hospital, Portland OR, Massachusetts General Hospital, Boston MA, Brigham and Women's Hospital, Boston MA, the Hospital of the University of Pennsylvania, Philadelphia PA, and The University of Alabama Birmingham Hospital, Birmingham AL.Recordings were collected with a Nihon Kohen Neurofax (the Cleveland Clinic), Cadwell Arc (Oregon Health and Science University Hospital), or Natus Quantum amplifier (all other clinical sites) and acquired with a sampling frequency of at least 500 Hz.After preliminary inspection, 27 patients with prior resections, highly abnormal sleep patterns, or whose data were overcontaminated with excessive epileptiform activity or technical artifacts were excluded from analysis.A further 35 patients were excluded from further analysis after channel selection and quality control procedures, described below, resulted in them having with fewer than two clean bipolar channel pairs.After quality control, 236 patients (126 female, 101 male, 9 unknown, 36.5 ± 13.3 (mean ± stdev.) years old, 22 unknown age) were included in the final analysis.

Channel localization and inter-channel distance.
As part of standard of care, a pre-operative T1-weighted high-resolution structural MRI and inter-operative CT scan were obtained for each patient.SEEG contacts were localized as previouslt described [102].Briefly, the post-implant CT volume was co-registered to the pre-implant MR volume in standardized 1mm isotropic FreeSurfer space with the general registration module [103] in 3D Slicer [104] and each contact's position was determined by manually marking the contact centroid as visualized in the co-registered CT volume.Cortical surfaces were reconstructed from the MR volume with the FreeSurfer recon-all pipeline [105].Each transcortical bipolar channel was assigned to the white-gray surface vertex nearest to the midpoint of the adjacent contacts.Spherical folding-based surface registration [106] was used to map the individual subjects' bipolar channel vertices to the standardized fsaverage ico7 [105] and 32k_FS_LR grayordinate [107] atlases as well as to parcels of the HCP-MMP parcellation [33].The HCP-MMP parcel assignments of each channel were used to determine the approximate white matter fiber length of the connection between them, using the mean values for a cohort of 1,065 healthy subjects we have previously reported [31].

Telemetry preprocessing and channel selection.
For each patient, an average of 185.3 ± 112.0 hours of raw telemetry were analyzed.Data were anti-aliased filtered at 250 Hz and resampled to 500 Hz.To remove line noise, a notch filter was applied at 60, 120, 180, and 240 Hz.On each electrode shaft, data were re-referenced to a bipolar scheme by taking the difference of potentials in adjacent contacts, yielding 30,769 total channels.To ensure that bipolar channels were independent, contact pairs were selected such that no two pairs shared a contact, eliminating just under half of all channels from further consideration.Peaks in the delta band (0.5-2 Hz) signal during high-delta periods were identified and the peri-delta-peak high-gamma band (70-190 Hz) envelope obtained.Transcortical contact pairs were identified by (1) having a midpoint close to the gray-matter-white-matter interface (<5 mm), (2) having a high anticorrelation between the delta-band amplitude and the high gamma band  envelope during delta-band peaks, indicative of non-pathological Down states with a surfacenegative LFP deflection that quieted spiking, and (3) having high delta-band amplitude with and average delta peak of at least 40 μV, indicative of robust inversion in adjacent contacts.Only contacts with 5 mm pitch were analyzed to standardize absolute bipolar amplitudes and contact pairs located deep in the white matter or in subcortical gray matter were excluded from analysis.All channels were visually inspected to ensure they contained broadly normal non-pathological waveforms such as pathological delta activity.This procedure yielded 4,903 candidate channels which underwent NREM scoring, described below, yielding 9.2 ± 7.8 hours of putatively clean scored NREM sleep.After scoring, which excluded continuous periods of artifactual or pathological signals, remaining epileptogenic activity was detected with spike template convolution and gross artifacts, with amplitudes > 300 μV.Channels which contained persistent automatically detected interictal spikes or gross artifacts were removed from further analysis, yielding 3,665 channels.Periods containing interictal spikes or gross artifacts in any channel were excluded from analysis in all channels.Finally, after cross-spectral matrices were obtained for each subject as described below, channels with autospectral power more than 3 standard deviations above the pooled channelwise mean in any examined frequency were excluded, leaving 3,588 channels and 37,123 within-patient bivariate channel pairs in the final analysis.

NREM selection.
Clean periods of sleep were first manually classified based on ultradian delta power.Within manually identified nightly periods of sleep, NREM sleep was automatically identified by the presences of slow oscillations (SOs) in a manner derived from standard clinical procedures [108].The density of SOs was ascertained in each 30-second frame and NREM marked for the frame if at least 35% of the frame contained SOs in at least one bipolar channel.The reduction in the number of channels and increase in the percent of frame threshold relative to the clinical standard of 20% of the frame across scalp channels was to accommodate the increased spatial heterogeneity of sleep graphoelements observed intracranially, relative to the scalp [109].For this purpose, slow waves were automatically detected with the method described in [109], briefly, consecutive delta-band zero-crossings occurring within 0.23-3 s were detected and the top 20% of peaks were retained as SOs.Only periods without detected interictal spikes or gross artifacts in any bipolar channel were included in subsequent analysis.

Coherence analysis
4.6.1 Mock experimental data.Coherence analysis was performed on both experimental and model data in an identical manner following preprocessing.The model data was preprocessed by averaging the voltages from all cells belonging to each of the 180 parcels (all 6 layers, both excitatory and inhibitory) to mimic what would be recorded by an SEEG electrode.Additionally, white gaussian noise was added to each of the 180 parcel averaged signals.This was done to mimic experimental noise picked up from both other non-neural bodily signals (i.e., muscle movements) and noise from the recording electrode itself.[54] demonstrated that ECOG arrays show an SNR ranging from 2.22 to 4.37 in relation to eye blink movement artifacts.[55] further calculated the SNR of slow wave sleep Up states (signal) v.s.Down states (noise); they found low frequencies (less than 30 Hz) to have a steady SNR of 4-9, depending on the type of electrode used.At higher frequencies, the SNR dropped severely.With these studies in mind, we implemented an SNR of 5 to account for both body and electrode noise that our model data lacks.

Signal processing and epoch sampling.
A frequency resolved map of each patient's (or simulation's) neural covariability was achieved by estimating the complex cross-spectral matrix among bipolar SEEG channels (for empirical data) or model parcel averages.Narrowband signals were extracted with a series of second-order resonator filters (matlab's iirpeak) with peaks and 3 dB attenuation bandwidths shown in Table A in S1 Text.For each frequency the complex analytic signal was extracted via the Hilbert transform.Calculating the variancecovariance matrix of analytic signals yields the cross-spectral matrix.Cross-spectral matrices were normalized into magnitude-squared coherence by squaring the complex modulus of each cross-pairwise cross-spectrum and dividing by the product of the constituent autospectra.
For experimental data, an adaptive procedure was used to obtain stable representative estimates of spontaneous covariability for each patient.Continuous minutes of telemetry were iteratively sampled at random, without replacement, from artifact-free NREM sleep.Sampled minutes of telemetry were appended and prepended with 1.5 seconds of data for filtering and analytic signal extraction.These data were removed prior to covariance estimation to avoid edge artifacts.For each frequency, the coherence matrices were obtained as described above for the incoming minute, nþ1 P nþ1 , the cumulatively sampled telemetry, 1 P nþ1 , and, for the second sampled minute onward, the prior cumulatively sampled telemetry, 1 P n .The correlation matrix distance (CMD) [110] P nþ1 for all previously sampled minutes.The coherence matrix was considered to have converged when the CMD between 1 P n and nþ1 P nþ1 was less than 1 × 10 −5 and the cumulative cross-spectral matrix was then retained.To dilute the influence of the first sampled minute, this procedure was repeated ten times, and the mean of the resulting cross-spectral matrices was used for subsequent analyses.

Statistical analyses and model fitting.
The effect of inter-parcel fiber distance on coherence was assessed by fitting a linear mixed-effects model: log(coherence) distance + 1 + (1-patient), which accounted for the fixed effects of the predictor variable distance and an intercept term, as well the random effect intercept term grouped by individual (this term was excluded for simulated data).Because of the natural log transformation of the response variable coherence, the linear model is equivalent to an exponential model for the untransformed distance predictor.Therefore, we report the modeled slopes and intercepts as length constants (λ) and scaling terms (α, maximum coherence), respectively.Separate fits were made for each narrowband frequency.Because the parameters of the exponential trend appeared to change at distances greater that 100 mm in experimental, especially for higher frequencies, separate models were fit to electrode pairs 0-100 mm distant, >100 mm, as well as to all pairs.Despite having twice as many parameters, the two-domain model was favored over the one-domain model by both Bayesian and Akaike information criteria for all frequencies and behavioral states.Only the 0-100 mm model parameters are reported here for both experimental and simulation data.All mixed effects linear models were fitted with matlab's fitlme.To avoid the inherent skewness coherence coefficient distributions, coherence values were transformed into z-statistics [111] before model fitting or averaging for display then reverted into coherence coefficients after model fitting or averaging.
Finally, we take the 0.5 Hz lambda and the mean <2 Hz alpha values to describe the key rate of decrease and max coherence in the slow wave range as the primary characteristics of each coherence landscape to be compared across experimental and model data (Fig 9).

Fig 1 .
Fig 1. Model architecture.A) Connectivity between neurons in the same column follows the canonical circuit.Arrows indicate directed connections between layers.Neurons are never connected to themselves.B) Connectivity diagram between different cortical columns and thalamocortical cells in the core (TC) or matrix (TCa).C) Exponential decay of connection probability with distance between columns.Previously reported inter-parcel dMRI connectivity[31] shows an exponential relation with fiber distance.Column-wise connection probability is obtained by applying the dMRI exponential function to inter-column fiber distances (synthesized from their corresponding geodesic distances) and then adding the residual inter-parcel dMRI connectivity, estimated by regressing out the exponential trend (see Diffusion MRI guided connection probability for details).D) Estimated myelination index[34,35] throughout the cortex.A hierarchical index inversely proportional to this myelination index is assigned to each of the 180 cortical parcels designated in the HCP-MMP atlas[33].E) Excitatory corticocortical connections belong to one of 6 classes based on the myelination-derived hierarchical index of the pre-and post-synaptic neurons: within the same column, within the same parcel, weakly or strongly feedforward (from lower to higher hierarchical index) and weakly or strongly feedback (from higher to lower hierarchical index).Based on experimental reports on connectivity, weights are scaled according to the connection type by the factor shown in the matrices (see Hierarchically guided laminar connectivity for details).F) Parcel-wise connectivity as previously reported[31].Fpt stands for "fraction of probabilistic tractography" and represents the log 10 probability of connection resulting from fractionally scaling raw streamline counts (detailed in[31]).G) Generated model connectivity, which retains the parcel-wise structure from data in panel F (Pearson's correlation r = 0.81, p < 0.0001).Fpt in the model is computed as the log 10 of the ratio of the subsubsection 4.2.3).The resulting baseline model of SWS (Fig 2B) revealed synchronized and regular oscillations (*0.27Hz) between silent (Down) and active (Up) states (Fig 2b.2) in all cortical excitatory and inhibitory neurons (Fig 2b.1).This network dynamics was organized as the waves of spiking activity propagating across the entire cortex-slow-waves.(This can be seen in S1 Video).The slow-waves were further highly synchronized across layers (Fig A in S1 Text).Differences between layers include a relatively less active L4.This is probably due to L4 receiving within-column connections from only L6, while other layers have inputs from two or more other layers (see Fig 1E).

Fig 2 .
Fig 2. Baseline model activity in awake and sleep states.A) Wake state.a.1) Membrane voltages of all neurons from layer 2 over 5 seconds of simulation; a.2) Layer average voltage over time, simulating LFP; a.3) Activity of two representative neurons from layer 2, both showing irregular tonic firing; a.4) Voltage histogram of all neurons over the whole simulation time, note approximately -60mV peak representing baseline membrane voltage; a.5) Power spectral density (PSD) of the average voltage revealed a distinct 1/f phenomenon typical for in-vivo recordings.B) Slow-wave dynamics.b.1) Membrane voltages of all neurons (excluding medial wall) from layer 2 over 30 seconds of simulation, revealed synchronized bands of activity during Up states; b.2) Average voltage of all cortical layers (dashed black line) and layer 2 neurons (solid blue line, excluding medial wall) over time, with nearly identical behavior.Red triangles above and below the trace mark global Up and Down states, respectively, from layer 2 (coincident with global Up and Down states from the average of layers); a.3) Activity of two representative neurons from layer 2, both showing synchronized transitions between Up and Down states; a.4) Voltage histogram of all neurons over the whole simulation time, revealed the characteristic bi-modal distribution caused by Up and Down state alternations during SWS.Dashed vertical lines labeled V − and V + indicate the voltage thresholds used to detect Down states and Up states, respectively; b.5) Distribution of the Up state onsets and offsets for all neurons over the whole simulation time.Narrow histograms indicate highly synchronous initiation and termination of the Up states.C) Zoom into the first Up state from panel b.1, with neurons sorted from earliest to latest onset time, and a single cell voltage from panel b.2, showing the transition from Down to Up state, steady firing during the Up state, and transition to Down state.D) Latency map, calculated as the onset delay of each neuron with respect to the earliest onset, for each Up state in panel b.2 (see Onset/offset detection in Methods for 5 and 3f.5) compared to baseline model (Fig 2b.5).Further decreasing density below 30% caused a considerable drop in amplitude and participation (Fig 3f.1 and 3f.4).In the extreme case of only 10% of connections remaining (P = 0.1, Fig 3E) (see in S1 Text: Effects of connections density, range and delays, Fig F in S1 Text and S2

Fig 3 .
Fig 3. Effect of connection density on SO dynamics.A) Example target cortical column (yellow) with all columns connected to it through any layer (green) for different connection densities.Density is reduced by decreasing P, which denotes the probability of preserving a connection from the original dMRI-based connectivity.With P = 1 all connections are preserved, while for P = 0.5, P = 0.3 and P = 0.1 each connection is preserved with a 50%, 30% or 10% probability, respectively.B) Percent of the original number of connections retained for different values of P. C) Distribution of connection distances, or lengths, for different values of P. Red horizontal lines indicate medians, bottom and top box edges indicate the 25th and 75th percentiles, respectively, whiskers extend to the most extreme data points not considered outliers, and outliers are plotted individually using the '+' marker symbol.Note, connection density is reduced uniformly across all lengths.D-E) Individual analysis of SO dynamics for different values of P. Subpanels as in Fig 2. F) Fig G in S1 Text): (a) by changing the maximum delay and scaling all others proportionally, and (b) by setting one uniform delay to all connections.

Fig 4 .
Fig 4. Effect of connection range on SO dynamics.A) Example target column (yellow) with all its connected columns (green).The blue area indicates the connection range for the corresponding radius R. No connections are made in the black region outside the radius.R = 226mm encompasses all the original dMRI-derived connections.B) Percent of the original number of connections preserved for different values of R. C) Distance distribution for different values of R, with inset zooming into radii below 10mm.R imposes a maximum connection length and truncates the distance distribution accordingly.Red horizontal lines indicate medians, bottom and top box edges indicate the 25th and 75th percentiles, respectively, whiskers extend to the most extreme data points not ), the model displayed a number of clearly identifiable local slow waves originating and reliably traveling around the hemisphere in different directions (Fig 7F).Occasionally, these local slow waves would coincide to produce a global slow wave involving the entire cortex (see Fig 6b.2, where red downward triangles mark the four detectable global UP states analyzed in Fig 7F).

Fig 5 .
Fig 5. Local SO for small connection radius (R = 2.5mm).A) Ten cortical areas with a 5mm radius, that were used to calculate local dynamics.B) Average membrane voltage of layer II neurons, as in Fig 4e.2.C-E) For each region in (A), subpanels show (C) the single-cell voltage for two neurons in the area, (D) the local field potential (LFP) for the 5mm area, and (E) heatmap of individual voltages of all neurons in the area.Note that even within these very small regions, there are still subgroups of independently synchronized neurons.F) Latency map, calculated as the onset delay of each neuron with respect to the earliest onset, for each Up state identified in Fig 4e.2.Note that even the most global of slow waves has very low participation, around 35%, and can be seen to consist of many extremely small initiation sites that generally do not spread.https://doi.org/10.1371/journal.pcbi.1012245.g005

Fig 6 .
Fig 6.Effect of cortical excitatory synaptic strength on the network dynamics.A) Reduction of all cortical connections by 3x.B) Reduction of only connections longer than 10 mm by 5x.All subpanels as in Fig 2. C) Summary of changes in frequency, amplitude, speed, participation, and onset/offset distribution spread when all connections are reduced as shown in (A).Increasing all weights shows little effect while decreasing all weights shows decreased amplitude and participation, with increases in onset/offset distributions (decreases in synchrony).The frequency of slow waves also becomes more variable.Note, Increasing all weights 5x results in the ablation of slow waves due to constant Up state, i.e., SO characteristics cannot be

Fig 7 .Fig 8 .
Fig 7. Mixed global (Type I) and local (Type II) Slow Waves: Connections greater than 10mm reduced 5-fold.A) Ten cortical areas with a 5mm radius, that were used to calculate local dynamics.B) Average membrane voltage of layer II neurons.C-E) For each region in (A), subpanels show (C) the single-cell voltage for two neurons in the area, (D) the local field potential (LFP) for the 5mm area, and (E) heatmap of individual voltages of all neurons in the area.This demonstrates that while there is some alignment, up states are not strongly global.F) Latency map, calculated as the onset delay of each neuron with respect to the earliest onset, for each Up state.Note, that while some Up states

Fig 9 .
Fig 9. Functional Connectivity and Coherence Analysis of Experimental v.s.Model data.Coherence analysis was performed on all data sets according to the procedure described in Coherence analysis.A) The parcel-to-parcel coherence matrices are first averaged in the slow wave frequency band and scaled by parcel-to-parcel distances to determine functional connectivity via percent of connected parcels and number of communities (determined by Louvain community detection).B) Full (unaveraged, unscaled) coherence matrices are then fitted with an exponential function to determine the full shape of the coherence landscape across distance and frequency.The resultant first (0.5 Hz) Lambda and mean Alpha in the slow wave frequency range (<2 Hz) are taken to describe the shape of the coherence landscape, and plotted in (B).(A.1 and B.1) show zoom-ins of each respective plot; the labels show the radius and factor of strength reduction-e.g., 10/5 indicates a model where connections longer than 10 mm are reduced by a factor of 5.This shows that the models with primarily mixed slow waves are the closest fit to experimental results across all 4 metrics.Full 3D coherence plots across distance and frequency are shown for the experimental data (C), the 10mm / 5 model (D), and the global slow wave base model (E), showing that the global slow wave model has a fundamentally different shape (lacking a dependence of coherence on distance) and much higher levels of coherence at low frequencies.https://doi.org/10.1371/journal.pcbi.1012245.g009 ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl } miniature postsynaptic potentials ð12Þ with d(n) governed by

4 . 4 . 2
Local field potentials, amplitude and frequency.In Figs 2, 3 and 4, the local field potential (LFP) was obtained as the average voltage of Layer 2 at each timestep scaled by a factor of 20.In Figs 5, 7 and F in S1 Text the LFP for each local region was further smoothed using a running mean with 20ms window.

4 . 4 . 4
Latency and participation.Latency maps are shown in Fig 2F (also Figs 5F and 7F and FF in S1 Text) for each global Up state, representing the delay for activation of each cell after the earliest onset in that Up state.The 3D representation of latencies allows for a visual identification of the source (or sources) of each global Up state, characterized by clusters of adjacent neurons with near-zero latencies.

Table 1 . Summary of results for each network manipulation.
Decreasing connection density resulted in decreased slow wave frequency, amplitude, synchrony, and spread, but slow waves remained global.Decreasing connection radius resulted in decreased slow wave amplitude, speed, participation, and spread.Increasing delays resulted in effects are seen only for manipulations far past biological plausibility, as indicated with a single down arrow for frequency, amplitude, and participation.Decreases in synaptic strength caused reduced amplitude, synchrony, frequency and participation.Finally, decreasing the strength of only long-range connections resulted in traveling local slow waves with decreased amplitude, speed, and synchrony.

Table 4 . Parameter changes between sleep and wake stages.
, was calculated between 1 P n and 1