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 yellow star indicates the functional connectivity metrics (Fig 9A) and the values of λ and α obtained by fitting an exponential model (Fig 9B), based on in vivo data.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

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

4.2.1 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–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²⁺ 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

https://doi.org/10.1371/journal.pcbi.1012245.t002

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]): (1) (2) (3) 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 (4) (5) 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 (6) where α and w₀ are constant parameters that control the non-linearity of the spike-generating mechanism. Parameters μ₁, μ₂ < 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 β_xI_total(x_n, u_n) was manually prevented from reaching 0 by imposing a lower bound of 10⁻⁴.

For inhibitory neurons IN, function f_IN in Eq 4 is given by (7)

A spike is generated if and only if 0 < x_n ≤ α + y* + ψ(β_xI_syn) ∧ x_n−1 ≤ 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₀ (low sensitivity) when voltage x_n crosses threshold -0.5 (i.e. Up state) and switches to K₁ > K₀ (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): (8) 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: (9) (10) (11)

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 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 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: (12) with d(n) governed by (13) for AMPA and GABA_A and (14) for NMDA type synapses. Here, determines the decay rate of synaptic conductance and regulates the overall synaptic strength. To keep total synaptic input per cell constant, coefficients and are scaled by the number 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 is a spontaneous event from a random Poisson process with mean frequency μ_X, and is the constant by which conductance increases in an event (), 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 is the time of the most recent pre-synaptic spike. This modulation reduces the mini rate to 20% of its maximal value after a pre-synaptic spike, and then allows recovery to with time.

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Parameter definitions and values for conductance-based models are based on previous modeling of single-cell dynamics [23, 36, 93, 94].

https://doi.org/10.1371/journal.pcbi.1012245.t003

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²⁺ 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].

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

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Parameter changes between sleep and wake stages.

https://doi.org/10.1371/journal.pcbi.1012245.t004

4.3.1 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.2 Diffusion 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. Column-wise 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.

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.

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.

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.3 Onset/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.

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.

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.

4.4.5 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² > .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² ≤ .25.

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 over-contaminated 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.

4.5.2 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].

4.5.3 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 (70–190 Hz) envelope during delta-band peaks, indicative of non-pathological Down states with a surface-negative 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.

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

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.

4.6.2 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. Narrow-band 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 variance-covariance 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, , the cumulatively sampled telemetry, , and, for the second sampled minute onward, the prior cumulatively sampled telemetry, . The correlation matrix distance (CMD) [110], was calculated between and , and between and . Starting with the third sampled minute, incoming minutes were rejected if the CMD between and , was more than 3 standard deviations greater than average CMD between and for all previously sampled minutes. The coherence matrix was considered to have converged when the CMD between and was less than 1 × 10⁻⁵ 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.

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