Excitable dynamics of NREM sleep: a unifying model for neocortex and hippocampus

During non-rapid eye movement (NREM) sleep, the neocortex and hippocampus alternate between periods of neuronal spiking and inactivity. By directly comparing experimental observations with a mean field model of an adapting, recurrent neuronal population, we find that the neocortical alternations reflect a dynamical regime in which a stable active state is interrupted by transient inactive states (slow waves) while the hippocampal alternations reflect a stable inactive state interrupted by transient active states (sharp waves). We propose that during NREM sleep, hippocampal and neocortical populations are excitable: each in a stable state from which internal fluctuations or external perturbation can evoke the stereotyped population events that mediate NREM functions.

Sleep function relies on internally-generated dynamics in neuronal populations. In the neocortex, non-rapid eye movement (NREM) sleep is dominated by a "slow oscillation" 1 : alternations between periods of spiking (UP states) and periods of hyperpolarization (DOWN states) that correspond to large "slow waves" (or "delta waves") in the local field potential (LFP) 2,3 ( Figure 1A,B, Supplemental Figure 1). In the hippocampus, NREM sleep is dominated by sharp wave-ripple dynamics: periods of spiking (SWRs) separated by periods of relative inactivity (inter-SWRs) 4 (Figure 1E,F). Slow waves and SWRs are bidirectionally and weakly coupled, in that each is more likely following the other [5][6][7][8] . The functional importance of these dynamics is well established: both slow waves and SWRs have been observed to perform homeostatic maintenance of the local synaptic network in the two regions [9][10][11] , and their temporal coupling has been found to support the consolidation of recently learned memories [12][13][14] .
However, it's unclear how the state of neuronal populations in the two regions promotes the generation of their respective dynamics, or how population state supports the propagation of neural activity between structures.
To study the state of hippocampal and neocortical populations during NREM sleep, we used an idealized model of an adapting recurrent neuronal population ( Figure 1C,G). Models with recurrence and adaptation have been directly matched to neocortical UP/DOWN alternations during anesthesia and in slice preparations [15][16][17] . These studies found that the UP/DOWN alternations in slice are adaptation-mediated oscillations 16 , while those in the anesthetized animal reflect noise-induced switches between bistable states 15 . However neuronal dynamics during NREM sleep in naturally sleeping rats 18 are distinct from those seen in anesthesia/slice 19 . With our adapting recurrent population model we are able to describe how effective physiological parameters determine the properties of alternation dynamics in a neuronal population. This framework allowed us to identify parameter domains that match the NREM data and, further, enabled description and understanding of both neocortical and hippocampal alternation dynamics with the same model.
We report that neocortical and hippocampal populations are neither endogenously oscillatory nor bistable during NREM sleep, but are excitable: each population rests in a stable state from which suprathreshold fluctuations can induce a transient population event that is terminated by the influence of adaptation. Specifically, the neocortex maintains a stable UP state with fluctuation-induced transitions to a transient DOWN state (slow waves), while the hippocampus rests in a stable DOWN state with fluctuation-induced transitions to a transient UP state (SWRs). Under the influence of noise, each region can generate its respective population event spontaneously (due to internally-generated fluctuations) or in response to an external perturbation (such as input from another brain structure). The result is alternations between active and inactive states in both structures with the asymmetric duration distributions observed during NREM sleep ( Figure 1D,H). We further observe that variation in the depth of NREM sleep corresponds to variation in the stability of the neocortical UP state. Our findings reveal a unifying picture of the state of hippocampal and neocortical populations during NREM sleep, which suggests that NREM function relies on excitable dynamics in the two regions.

UP/DOWN dynamics in an adapting excitatory population model
UP/DOWN alternations are produced in models of neural populations with recurrent excitation and slow adaptive feedback 15-17,20-24 25 . In our model, neuronal population activity is described in terms of the mean firing rate, ! ! , subject to a slow negative feedback (i.e. adaptation), ! ! (Figure 2A, see Supplemental Info for details).
Equations 1-2 describe how ! and ! evolve in time as a function of the net input to the population: the sum of the recurrent excitation with weight ! and a background level of drive with a tonic parameter !, and noisy fluctuations ! ! , minus adaptation weighted by gain parameter ! (See Supplemental Info for parameter interpretation). ! ! !"#$% is the "cellularlevel" input-output function, which defines the rate of the population given constant net input.
Similarly, ! ! ! defines the level of adaptation given a fixed population rate. To enable the analytical treatment of model dynamics in the following section, both ! ! !"#$% and ! ! ! are taken to be sigmoidal functions.
Model dynamics can be represented as a trajectory in the !-! phase plane 26  In the oscillatory regime (Figure 2Ci)  and modeled dynamics. The domain of high similarity was degenerate and remained in the Excitable UP regime with variation in the "fixed" parameters, ! ! , !, and the amplitude of the noise (Supplemental Figure 7). We thus found that NREM sleep in the rodent neocortex is characterized by an Excitable UP regime: a stable UP state with noise-induced transitions to a transient DOWN state.

Hippocampus is in an Excitable DOWN regime during NREM sleep
Since the burst-like dynamics of SWR is reminiscent of the Excitable DOWN regime of our model, we asked whether these patterns could also be explained by the same principles.
InterSWR durations are much longer (mean = 2.0±0.22s) compared to SWR events (mean = 0.06±0.005s ),and more variable (CV InterSWR = 1.3±0.10; CV SWR = 0.33±0.04) ( Figure 4E) suggesting a stable DOWN and transient UP state (SWR). We applied the duration distribution matching procedure to the SWR/inter-SWR duration distributions and confirmed that the !-! model can also mimic SWR dynamics, with a band of high data-model similarity in the Excitable DOWN regime ( Figure 4G). Interestingly, our idealized model is not able to capture the short-interval inter-SWR periods associated with occasional SWR "bursts" (Supplemental Figure   7), which suggest the presence of separate SWR-burst promoting mechanisms, possibly arising from interactions with the entorhinal cortex or spatially traveling patterns of SWRs in the hippocampus 27,28 . Accordingly, while the mean ratio and CV SWR of the best fitting model regime were within 2.5 standard deviations of those observed in vivo, the CV of inter-SWR periods was larger than expected from the model (i.e. CV>1). This finding suggests that during NREM sleep the hippocampus is in a stable DOWN-like state, from which internal 'noise' or an external perturbation can induce population-wide spiking events.

Changes in neocortical state correspond to changes in UP state stability
For our initial analysis of the neocortical NREM data we assumed that model parameters

Inhibition stabilizes a low-rate UP state and allows perturbation-evoked slow waves
Our previous analyses considered a constant (stationary) source of noise that produced "spontaneous" transitions out of stable states in our model. We now consider a brief input that evokes a transient event. For the hippocampal-like Excitable DOWN regime, a brief increase in drive will evoke a transient UP state, (i.e. a SWR, Supplemental Figure 9). In the absence of noise, perturbations must be of sufficient magnitude (i.e. suprathreshold). With noise, the probability to evoke a SWR increases with magnitude of the perturbation (Supplemental Figure   9). A converse situation is apparent for the neocortical-like Excitable UP regime --a brief decrease in drive is able to evoke a transient DOWN state (i.e. a slow wave, Supplemental Figure 9). However, as long-range projections tend to be excitatory, we wondered how an excitatory perturbation might evoke a neocortical UP->DOWN transition.
Neuronal spike rates during the UP state are generally low 18 with balanced excitatory and inhibitory synaptic inputs 30 . Previous work has shown that models with fast inhibition and slow adaptation can give UP/DOWN alternations in the same four regimes described above 15 with a low-rate UP state that is stabilized by feedback inhibition 31,32 . We hypothesized that the effect of inhibitory cells may support excitation-induced UP->DOWN transitions, and included an inhibitory population (! ! ≈ ! ! ) into our model ( Figure 6A): where adaptation acts on the excitatory population and ! !,! !!"#$ and ! !,! !"#$% are threshold power law I/O relations, as seen in the in vivo-like fluctuation-driven regime 33 (Supplemental Info).
Given that adaptation is slow we can treat ! as frozen and visualize model dynamics in the ! ! -! ! phase plane ( Figure 6B). The fixed point value of ! ! as a function of drive describes the effective I/O curve of the network (! !! , Figure 6C). Like the excitation-only model, strong recurrent excitation induces bistability at low levels of drive (Supplemental Figure 10). In the bistable condition, the ! ! -! ! phase plane shows stable UP and DOWN state fixed points, separated by a saddle point ( Figure 6B,C). With ! dynamic, the model can have steady state fixed points on either the UP or the DOWN branch of the I/O curve, resulting in the same regimes as the two-variable model described above 15 (Supplemental Figure 10).
We investigated Excitable UP dynamics in the adapting, inhibition-stabilized model 33 ( Figure 6D,E,F). Consider a transition from the UP to the DOWN state ( Figure 6F). As further implication is that a sufficiently strong excitatory input to the excitatory population ( Figure   6E) can recruit sufficient inhibition to force the entire network into a DOWN state. This threshold effect is seen as a trajectory in the phase plane that separates the basins of attraction of the UP and DOWN state (i.e. a separatrix, Figure 6F). The separatrix emerges (in reverse time) from the saddle and curves around the UP state fixed point. From this visualization we see that a brief excitatory input to either population can push the trajectory out of the UP state basin of attraction ( Figure 6F). Thus, a transient DOWN state (i.e. a slow wave) can be evoked by an excitatory perturbation to either population, as well as due to drops in the excitatory population rate.

Discussion
To account for cortical dynamics during NREM sleep, we used a firing rate model that

UP/DOWN dynamics of the neocortical NREM slow oscillation
Despite the widely used term slow "oscillation" 1 34 and temporal correlations that emerge with strong recurrent connections 35 . Similarly, the level of afferent activity from thalmo-or cortico-cortical projections would be expected to fluctuate. We also note that while the isolated cortex can produce UP/DOWN state alternations 36 , we should consider the thalamocortical system for an understanding of slow wave dynamics in vivo 37 . Because the cortex and corresponding thalamic nuclei are highly interconnected, cortex and thalamus may transition UP and DOWN together and reflect interacting (as opposed to independent) systems. However, it was recently found that cortex tends to lead the thalamus into the DOWN state 38 . Future work should expand the model to include a thalamic population, which would also allow a better understanding of the interaction of slow waves with thalamocortical spindle oscillations 8,39,40 .
We found that the depth of NREM sleep reflects an evolution of the stability of the UP state in a manner that resembles the stages of NREM/SWS sleep in humans 41 . In light NREM sleep (N1 stage, using the human clinical term), long UP states are occasionally punctuated by neuronal silence-associated positive delta or slow waves, which can be localized at one or few recording sites across the cortical mantle 8  As the relevant parameter is the relative strength of adaptation and recurrence, the different nature of recurrent connectivity in the two regions may also be responsible for their differing dynamics. Strongly recurrent pyramidal cell populations are found in neocortical layer 5 and the hippocampal CA2 and CA3a subregions 48 , the loci of UP state and sharp wave initiation, respectively 49,50 . However, crucial differences exist between connectivity of neocortical layer 5 and hippocampal CA2-3 regions. The neocortex is a modularly organized structure. In contrast, the hippocampus can be conceived as a single expanded cortical module 51 . Excitatory connectivity in layer 5 is local (200 µm), dense (up to 80% connection probability), and follows a 'Mexican hat' excitatory-inhibitory spatial structure with strong local excitatory connections and spatially extensive inhibition 52 . In contrast, excitatory connectivity in the hippocampus is sparse and spatially extensive 48 , with local inhibitory connections 53,54 . While layer 5 excitatory synapses are relatively strong, the transmitter release probability of synapses between hippocampal pyramidal neurons is very low, resulting in comparatively weak synapses 55 . Together, these factors indicate that the effective strength of recurrence in the hippocampus is lower than that in neocortex, which would result in DOWN-dominated as opposed to UP-dominated dynamics, as are observed. To further understand the physiological factors responsible for the distinct NREM dynamics in the two regions will require experimental manipulations that independently manipulate adaptation, recurrent excitation, and excitability.

NREM function through stochastic coordination of excitable dynamics
According to the two-stage model of memory consolidation 56,57 , the hippocampus acts as a fast, but unstable, learning system. In contrast, the neocortex acts as a slow learning system that forms long-lasting memories after many presentations of a stimulus. The two-stage model proposes that recently-learned patterns of activity are reactivated in the hippocampus during SWRs, which act as a "training" signal for the neocortex, and that the neocortical consolidation of those patterns relies on SWR-slow wave coupling 5,58 . Excitable dynamics provide a mechanism for coordination of slow waves and SWRs ( Figure 7B): the excitatory kick of a hippocampal SWR can induce a neocortical UP->DOWN transition by briefly disrupting the neocortical excitatory/inhibitory balance, while the population burst at the neocortical DOWN->UP transition can induce a hippocampal SWR.
Extensive experimental evidence points towards temporal coordination between slow waves and SWRs. Slow waves in higher-order neocortical regions are more likely following SWRs 5,6 , and SWR->slow wave coupling is associated with reactivation in the neocortex 5,58,59 .
As is observed in vivo, the ability of a transient input to evoke a slow wave in our model is probabilistic, due to an interaction between the magnitude of perturbation, local noise, and the stability of the UP state. The probability of SWR->slow wave induction likely varies by brain state, cortical region, and even SWR spiking content. Further work to investigate how these factors shape SWR->slow wave coupling will likely shed light on the brain-wide mechanisms of memory consolidation.
How then, does a SWR-induced neocortical slow wave induce changes in the neocortex? Recent work has found that SWR->slow wave coupling alters spiking dynamics at the subsequent neocortical DOWN->UP transition 58 (aka the k complex), which acts a window of opportunity for synaptic plasticity that supports NREM functions 10,60-62 . Interestingly, the interaction between excitation and inhibition produces a transient (gamma-like) oscillation at the DOWN->UP transition in our model. This brief oscillation is reminiscent of the gamma (~60-150Hz) activity following slow waves in vivo 18 and may act to coordinate and promote plasticity between cell assemblies 63 .
In turn, the burst of neocortical activity during the k complex could induce a SWR in the hippocampus. The functional role of slow wave->SWR coupling is less well understood, but hippocampal SWRs are more likely immediately following slow waves in some neocortical regions -including the entorhinal cortex 5,7 . Slow wave->SWR coupling could provide a mechanism by which neocortical activity is able to bias SWR content, or another mechanism by which the SWR could bias neocortical activity at the DOWN->UP transition. Further, a SWRslow wave-SWR loop could produce the occasional SWR bursts not captured by our model of hippocampal SWR activity in isolation. Future work on regional or state-dependent differences in the directionality of slow wave-SWR coupling could provide insight into the physiological mechanisms that support memory consolidation.

Conclusions
Our results reveal that NREM sleep is characterized by structure-specific excitable dynamics in the mammalian forebrain. We found that a model of an adapting recurrent neural population is sufficient to capture a variety of UP/DOWN alternation dynamics comparable to those observed in vivo. The neocortical "slow oscillation" is well-matched by the model in an

Datasets
The

NREM Detection
Sleep state was detected using an automated scoring algorithm as described previously (Watson et  In the hippocampal dataset, manual NREM scoring as reported in Grosmark and Buzsaki 2016 was used for this study.

Slow Wave Detection
Slow waves were detected using the coincidence of a two-stage threshold crossing in two signals (Supplemental Figure 1A,B): a drop in high gamma power (100-400Hz, representative of spiking (Watson et al 2017)) and a peak in the delta-band filtered signal (0.5-8Hz). The gamma power signal was smoothed using a sliding 80ms window, and locally normalized using a modified (non-parametric) Z-score in the surrounding 20s window, to account for nonstationaries in the data (for example due to changes in brain state and noise), that could result in local fluctuations in gamma power. The channel used for detection was determined as the channel for which delta was most negatively correlated with spiking activity, while gamma was most positively correlated with spiking activity.
Two thresholds were used for event detection in each LFP-derived signal, a "peak threshold" and a "window threshold". Time epochs in which the delta-filtered signal crossed the peak threshold were taken as putative slow wave events, with start and end times at the nearest crossing of the window threshold. Peak/window thresholds were determined for each recording individually to best give separation between spiking (UP states) and non-spiking (DOWN states) (Supplemental Figure 1C). To determine the delta thresholds, all peaks in the delta-filtered signal greater than 0.25 standard deviations were detected as candidate delta peaks and binned by peak magnitude. The peri-event time histogram (PETH) for spikes from all cells was calculated around delta peaks in each magnitude bin, and normalized by the mean rate in all bins. The smallest magnitude bin at which spiking (i.e. the PETH at time = 0) was lower than a set rate threshold (the "sensitivity" parameter, Supplemental Figure 1D) was taken to be the peak threshold. For example, a sensitivity of 0.5 means that the delta peak threshold is set to the smallest threshold for which spiking drops below 50% of mean spiking activity. The window threshold was set to the average delta value at which the rate crosses this threshold in all peak magnitude bins. The gamma thresholds were calculated similarly, but using drops below a gamma power magnitude instead of peaks above a delta magnitude.
Once the thresholds were calculated, candidate events were then detected in the delta and gamma power signals, and further limited to a minimum duration of 40ms. Slow wave events were then taken to be overlapping intervals of both the gamma and delta events. DOWN states with spiking above the sensitivity threshold were thrown out.
Detection quality was checked using a random sampling and visual inspection protocol.
LFP and spike rasters for random 10s windows of NREM sleep were presented to a manual scorer, who marked correct SW detections, false alarms, and missed SWs. This protocol was used to estimate the detection quality (miss %, FA %) for each recording (Supplemental Figure   1E), and to optimize the detection algorithm. 1,085 -21,147 slow waves (i.e. UP/DOWN states) were detected per recording and used for subsequent analysis.

Model Implementation
Phase plane and bifurcation analysis of the model in the absence of noise was implemented in XPP, and a similar code was implemented in MATLAB for simulations of the model with noisy input, for the analysis of UP/DOWN state durations. Noise was implemented using Ornstein-Uhlenbeck noise.
Simulations of equations [1][2] and [3][4][5] were performed in Matlab using the ode45 solver, with input noise ! ! pre-computed independently for each simulation using forward Euler method with time step dt=0. as threshold crossings between high and low rate states. To avoid spurious transition detection due to noise, a "sticky" threshold was used: the threshold for DOWN->UP transitions was taken to be the midpoint between positive crossings of a threshold between the high rate peak of the rate distribution and the inter-peak trough, while the threshold for UP->DOWN transitions was the midpoint between the low rate peak of the rate distribution and the inter-peak trough.

UP/DOWN State Duration Matching
In vivo and simulated UP/DOWN state durations were compared using a non-parametric distribution matching procedure (Supplemental Figure 6). Similarity was calculated as