Emergence and fragmentation of the alpha-band driven by neuronal network dynamics

Rhythmic neuronal network activity underlies brain oscillations. To investigate how connected neuronal networks contribute to the emergence of the α-band and to the regulation of Up and Down states, we study a model based on synaptic short-term depression-facilitation with afterhyperpolarization (AHP). We found that the α-band is generated by the network behavior near the attractor of the Up-state. Coupling inhibitory and excitatory networks by reciprocal connections leads to the emergence of a stable α-band during the Up states, as reflected in the spectrogram. To better characterize the emergence and stability of thalamocortical oscillations containing α and δ rhythms during anesthesia, we model the interaction of two excitatory networks with one inhibitory network, showing that this minimal topology underlies the generation of a persistent α-band in the neuronal voltage characterized by dominant Up over Down states. Finally, we show that the emergence of the α-band appears when external inputs are suppressed, while fragmentation occurs at small synaptic noise or with increasing inhibitory inputs. To conclude, α-oscillations could result from the synaptic dynamics of interacting excitatory neuronal networks with and without AHP, a principle that could apply to other rhythms.


Introduction
Electroencephalogram (EEG) is used to monitor the brain activity in various conditions such as sleep [1,2], coma [3] or meditation [4] and to reveal and quantify the presence of multiple frequency oscillations [5] over time [6]. This analysis can be used to assess the level of consciousness or depth of unconsciousness of the brain during anesthesia [7,8]. For example, during general anesthesia under propofol, the α-band (8-12Hz) emerges as a dominant oscillation [7,9]. However, the precise mechanisms underlying the emergence or disappearance of this αband remain unknown. Interestingly, when the level of sedation becomes too high, the EEG shows that the α-band can get fragmented and even disappear replaced by burst-suppression, consisting of alternation between periods of high frequency activity and iso-electric suppression where the EEG is almost flat [7]. In general, large doses of hypnotic in prolonged anesthesia in rodents alters brain synaptic architecture [10], confirming the need to avoid over sedation. Burst-suppression is a motif associated with a too deep anesthesia and its presence could indicate possible post anesthetic complications, although it has been attributed to ATP depletion [11]. Recently, it was shown that the loss of the α-band announces the appearance of burst-suppressions [8], and could be used as a prediction marker for these suppressions. However, this causality between α-band suppression and burst-suppressions remains unexplained.
The α-band revealed by the EEG signal is associated with the local Up and Down states activity [12][13][14], which corresponds to depolarized and hyperpolarized neuronal membrane voltage [15]. Alternation between Up and Down states generates slow wave oscillations present in Non-REM sleep, as reported in slices electrophysiology [16] as well as using modeling approaches [17,18]. Similarly, the emergence of the α-band during anesthesia could result from network interactions, as proposed by models based on the Hodgkin-Huxley formalism [19][20][21].
Since Up and Down states dynamics reflect neuronal activity at the population level [15,22], we propose here to investigate the emergence and fragmentation of the α-band using a modeling approach based on synaptic short-term plasticity [23,24], which is often used to obtain estimations for burst or interburst durations [25][26][27]. These models, based on facilitation and depression have recently been used to evaluate the working memory capacity to remember a sequence of words [28].
Here, we use a mean-field neuronal model that accounts for both synaptic short-term dynamics and afterhyperpolarization (AHP) [29] resulting in a refractory period during which neurons stops firing after a burst. As a result, at a population level, AHP can modify the type of oscillations [30], from waxing and waning spindle oscillations to slow waves.
We first study a single, two and then three interacting neuronal networks, a minimal configuration revealing the coexistence of α-oscillations with switching between Up and Down states. As we shall see, only the neuronal population with AHP can trigger spontaneous switching between Up and Down states while the other one, without AHP, is at the origin of the αoscillations in the Up state. We also investigate the role of synaptic noise and model the effect of propofol with an excitatory current on inhibitory neurons.

EEG reveals the dynamics of the α-band during general anesthesia
General anesthesia can be monitored using EEG (Fig 1A) that often reveals a stable α-band which persists over time (Fig 1B). When anesthesia involves the propofol agent, it activates the inhibitory neurons, resulting in the emergence of a stable α-band. Increasing propofol leads to a deeper anesthesia characterized by transient removal of the α-band ( Fig 1C) leaving a dominant δ-band in the spectrum. α-band disappearance can be quantified by two parameters, defining a fragmentation level which accounts for the percentage of presence (persistence level) P α of the α-power over a period of time and the number of disruptions D α in the power band (see Methods). In the remaining part of this manuscript we shall develop a mean-field model based on synaptic properties to study this emergence and fragmentation.

Single neuronal population can exhibit α-oscillations or slow switching between Up and Down states
Short-term depression and facilitation generates locked α-oscillations. To analyze changes between a persistent α and a δ-band, we develop mean-field models of neuronal networks based on short term synaptic plasticity (Fig 2A). The first model consists of one excitatory neuronal network in which neurons are sufficiently connected so that they exhibit a synchronized activity. In this context, the average activity of the population represents each neuron's activity. We describe this neuronal network by three variables: the mean voltage h, the synaptic facilitation x and the depression y, resulting in a stochastic dynamical system (see Methods, Eq 3) showing bi-stability. The dynamics presents a first attractor corresponding to the Down state (hyperpolarized, low frequency oscillations) and a second attractor associated with the Up state (depolarized, high frequency oscillations). One fundamental parameter of the model (first equation in 3) is the connectivity level J, which quantifies the average number of connections (synapses) between neurons in the population, that we shall vary (Fig 2B). We found that such a system can generate a dominant oscillatory band where the peak value is an increasing function of the connectivity J (Fig 2C). In the present scenario the network dynamics can be locked into an Up state and then the dominant oscillation frequency is generated by the imaginary part of the eigenvalues at the Up state attractor (see paragraph Origin of oscillations in the Up state in Methods for details). Sustaining these oscillations requires an additive noise, otherwise the dynamics would converge to the Up state attractor. This result shows that the persistent oscillations are the consequence of the noise and of the synaptic properties as well as the biophysical parameters (S1 Table). Changing the synaptic properties can lead to a The EEG signal is composed of multiple bands as shown in the spectrogram (lower panel) composed of two major bands: the δ-band (0-4Hz) and the α-band (8-12Hz) tracked by its maximum (black curve) revealing the persistence of the αband during anesthesia. The color indicates power in dB/Hz (see colorbar). C: Same as A for a case of fragmented α-band. Data from the database VitalDB [31].
https://doi.org/10.1371/journal.pcbi.1009639.g001 fragmentation and disappearance of the band so that most of the remaining spectral energy is located in the δ-band, as quantified by the spectral edge frequency at 95% (SEF95, S1 Fig). In addition, varying the noise amplitude allows to either fragment the band (Fig 2D, σ = 7) or to increase its power and persistence (Fig 2D, σ = 15) but it does not affect the value of the peak of the dominant oscillation frequency (Fig 2E). Indeed, similarly to the EEG analysis (Fig 1), we quantified the fragmentation level for different noise amplitudes and found that it varies from (P α , D α ) = (52%, 34/min) for a noise amplitude σ = 5 to (86%, 16/min) for σ = 7 and to (100%, 0/min) for σ = 15, similar to the magnitude observed in the human EEG data. This fragmentation of the α-band results only from the changes in the noise amplitude and can occur even though the neuronal population is locked into the Up state, suggesting that the loss the αband could be independent from the switch between Up and Down states.
AHP with short-term depression and facilitation generates Up and Down states but no α-oscillations. The short-term plasticity equations proposed in the previous section did not generate switching between Up and Down states, thus we decided to add the consequences of AHP to the model (see Methods). With this additional effect, it is possible to explore a larger range of neuronal dynamics. In brief, we added an AHP component to the mean-field Eq (3) presented above (see Methods). The novel dynamics exhibits bi-stability characterized by random transitions between Up and Down states External stimulation in Up states cannot destabilize the oscillation rhythm. Adding an input current on the mean voltage h during the Up states simulates a situation where the observed network projects an excitatory input on a second network that would send a positive feedback when activated. The second network would only get activated by such stimuli when the first (observed) network is in the Up state. In previous studies with a 2D model (modeling only the firing rate and depression) we showed that such stimulus stabilizes the Up state [24]. Here we ran simulations for the cases with and without AHP where we added a constant input current only when the system was in the Up state (

Effect of inhibition on an excitatory neuronal network with and without AHP
To explore the range of oscillatory behaviors, we connected an inhibitory neuronal network to an excitatory one with or without AHP (see Methods, Eq 5). We also added a constant input current I i on the inhibitory population ( Fig 4A).
In the case where the excitatory network has AHP the dynamics exhibit switching between Up and Down states ( Fig 4B). Interestingly, by increasing the current I i we can modulate the fraction of time spent in Up state by the excitatory system (Fig 4C and 4D). Furthermore, when the value of the current I i injected to the inhibitory population reaches the threshold I i � 60, the excitatory population is locked in the Down state, therefore silencing the network ( Fig  4C and 4D right panel). However, independently of the value of I i , the Up state does not show any persistent α-oscillations ( Fig 4E). Finally, when connected to an excitatory population without AHP, the inhibitory network cannot induce switching between Up and Down states, and thus the dominant α-band generated by the excitatory population in the Up state remains stable ( Fig 4F).

Coexistence of Up and Down states and α-oscillations in two excitatory and one inhibitory networks
To define the conditions under which the Up and Down states can coexist with an α-band, we coupled two excitatory networks with an inhibitory network (see Methods, Eq 6). With such topological organization, we explored whether α-oscillations could be generated as described for the thalamo-cortical loop [7,13]. In our model, the thalamo-cortical excitatory subsystem is composed of two populations: 1) α and 2) U/D networks connected by reciprocal connections ( Fig 5A). They receive an inhibitory input from the inhibitory subsystem NR. The NR component sends reciprocal connections to the U/D component and can also be activated by an external stimulation I i . We study the mean of the three voltage components as a proxy for a recorded EEG. During general anesthesia with propofol, increasing the dose leads to a fragmentation and transient disappearance of the α-band. To assess under which conditions this phenomenon could be reproduced in the model, we followed the same protocol: we first investigated the effect of switching off all external stimuli, and then by increasing an injected current to the inhibitory neuronal component, we simulated an increase of the propofol concentration.

Spontaneous emergence of α-oscillations after the suppression of external stimuli.
To mimic the loss of consciousness at the start of a general anesthesia, we switched off external sensory inputs in the model. To account for this transition, we first added external stimuli modeled as excitatory input currents I ext = 300 + 20ξ (resp. I ext,2 = 100 + 20ξ) where ξ is a Gaussian white noise of mean 0 and variance 1. We applied I ext and I ext,2 to the three components of the model (Fig 5A, blue and green) for the first 40 seconds of the simulation. To model the loss of consciousness at the beginning of anesthesia, we set the external stimuli I ext and I ext,2 to zero for the rest of the simulation (Fig 5B). We found no dominant oscillatory band in the spectrogram during wakefulness (I ext > 0 from 0 to 40s). However, after external inputs suppression (I ext = 0) during the interval 40 to 120s, the network stabilizes in the Up state and a dominant α-band emerges (Fig 5B). Here external stimuli are modeled by direct inputs including a noise term ξ that accounts for the variations in the sensory inputs. Note that we could have obtained a similar result with constant values of I ext = 300 (resp. I ext,2 = 100). However, the noise term s x _ o x in equation for h x in the three populations x 2 {U/D, α, R} is essential for the emergence of the α-band. Without this noise term, the dynamics would stay trapped at the attractor A Up with no oscillations. To conclude, the present simulations suggest that the α-band emerges as a stable state once the external stimuli are suppressed. This behavior results from the interactions between the three networks.
External inputs on the inhibitory network modulate switching between Up and Down states. To better characterize switching between Up and Down states, we decided to investigate the role of inhibitory inputs. To this end, we increased the inhibitory input current I i (Fig 6A) and estimated the fraction of time the networks spend in Up and Down states. For I i = 0, we found that the dynamics is characterized by a large proportion of time spent in the Up state (99%), showing persistent α-oscillations (Fig 6B). If the transition from Up to Down is characterized by a disappearance of the α-band, the transition from Down to Up comes with a burst which can either lead to the emergence of an α-band or a return to the Down state. By increasing I i from 0 to 50 and 150, we also found that the fraction of time spent in Up states decreases from 99% to 89% and then to 4.5% We found that a small connectivity J UD = J α = 5 is associated with a large number of Down states (88%) and that transient bursts rarely lead to a stable α-band ( Fig 6F). By increasing the connectivity J UD = J α to 5.6 and 6.5, the fraction of Up states increases to 43% and 99% respectively ( Fig 6I) leading to stable Up states associated with a persistent α-band. To conclude, the amplitude of inhibitory inputs as well as the connectivity level are fundamental parameters that modify the proportion of Up states and thus the persistence of the α-oscillations.
Transient responses of the thalamo-cortical model to step and stairs inputs. To study, the possible thalamo-cortical response to a bolus of the hypnotic agent propofol, we simulated two protocols: 1) a step input (protocol 1) and 2) a stairs input (protocol 2) as shown in Fig 7A. To analyze the response to protocol 1, we ran simulations with N = 2500 iterations lasting T = 2min, where a strong injection is delivered by simulating a positive input current I i = 1000 on the inhibitory network (NR) lasting t i = 20s (Fig 7B and 7C). To quantify the response we focused on two parameters: 1) the duration t C after which the α-band disappears after the step function begins. 2) The duration t U after which the α-band reappears after the end of the step function. Interestingly, for some realizations the α-band does not disappear (Fig 7C), we thus characterized this effect by the collapse probability p C . We found a probability value of p C = 53% and a duration t C = 9.42 ± 5.36s while t U = 4.39 ± 2.58s (Fig 7D). The histogram for t C is characterized by an abrupt decay at 20s confirming that α-band suppression can only occur during the stimulation period. However, the time t U is dominated by an exponential decay, a classical feature of dynamical systems driven by noise.
Moreover, each network has a different contribution to the EEG: the U/D excitatory network with AHP shows an α-band while the inhibitory network NR only exhibits weak power in the slow δ-wave region for frequencies � 1Hz.  measuring its re-emergence, as well as the probability to collapse p C are given by p C = 40%, t C = 53.39 ± 4.44s and t U = 4.20 ± 2.43s (Fig 7G). Similarly, the histogram of t C is characterized by an abrupt decay at 60s confirming that the suppression of the α-band can only occur during the stimulation period, while the histogram of t U is dominated by an exponential decay. These results quantify the impact of the intensity and duration of a propofol bolus on the α-band.

Discussion
We presented here a parsimonious computational approach based on coarse-grained neuronal network models in order to generate α-oscillations. Interestingly, a single neuronal population driven by synaptic short-term plasticity can elicit oscillations at a defined frequency, which directly depends on the value of the network connectivity: a higher connectivity generates faster oscillations (Fig 2B and 2C). Interestingly, we show here that stable α-oscillations result from the combination of network connectivity, the overall synaptic strength and the intrinsic biophysical properties of the network, leading to the emergence of a focus attractor, around which the mean population voltage oscillates in the phase-space S7(A) and S7(B) Fig. Moreover, we showed that spontaneous switching between Up and Down states in a single neuronal population is modulated by AHP and also that the network connectivity controls the proportion of Up vs Down states: a higher connectivity J results in a dominant percentage of time spent in Up states (Fig 3C and 3D). The stability of the oscillations during Up states for a population without AHP could result from the intrinsic network regulation: indeed, interactions between inhibitory interneurons and hippocampal pyramidal excitatory neurons can redistribute the firing load to maintain the oscillation frequency even when up to 25% of the synapses are deactivated [32].
When we added an excitatory input current on an inhibitory population coupled to an excitatory population, the proportion of time spent in Up states decreased and, after reaching a threshold value I i = 60, the network became completely silenced, characterized by Down states only (Fig 4C and 4D). When coupling two excitatory populations and one inhibitory population (Fig 5A), α-oscillations, generated by the excitatory component without AHP, coexisted with spontaneous switching between Up and Down states induced by the excitatory population with AHP, as summarized in Fig 8. Stimulating the inhibitory population modulates directly the proportion of Up vs Down states and thus induces a fragmentation of the αband (Fig 6B-6E).
Finally, we suggest that synaptic noise has two main roles on the network properties: 1) increasing the noise intensity stabilizes the α-band ( Fig 2D) and 2) when an external stimulation is applied to the inhibitory system in a step or stairs input, the network can react with two opposite behaviors: either the network activity collapses, leading to a suppression of all oscillatory bands in the EEG or a stable persistent α-band emerges during the entire stimulation (Fig 7). Finally, we propose that three connected neuronal populations are sufficient to generate an α-oscillation that could be fragmented by increasing the inhibitory pathway, as suggested during general anesthesia [7]. The present model could be generalized to other brain rhythms to study the emergence and disappearance of other oscillations such as the θ-oscillations occurring during REM sleep [2,33,34].

Modeling the dynamics of the α-band
For several decades, the origin of the α-band [14] had remained unclear. Early modeling efforts using the Hodgkin-Huxley framework [30] suggested a key role of ionic currents such as sodium, potassium currents, low threshold calcium, AHP and synaptic currents (GABAs, AMPA) that could reproduce various patterns of oscillations [19] as well as initiation, propagation and termination of oscillations (see also [35]). By varying the GABA conductances other models have reproduced the dominant α-oscillation observed in propofol anesthesia [20,21]. Indeed GABA conductances regulate the firing frequencies and synchronization of pyramidal neurons [36]. In contrast to conductance based models, we used here short-term synaptic plasticity driven by noise and found that the α-band can be generated only when the mean voltage reaches the Up state, suggesting that the ionic mechanisms are not necessary to generate the αband, but contributes to the termination of the Up states. Furthermore, adding an input current to the inhibitory population allows to modulate the transitions between Up and Down states (Figs 4C, 4D and 6B-6E). In addition, we found that the α-band can be stabilized by increasing the noise amplitude, while the peak frequency of the α-band was unchanged ( Fig  2D-2E). Thus, we propose that the synaptic noise could be responsible for the stabilization of the α-band, as quantified by the present fragmentation analysis (Fig 2D-2E). Finally, we reported a second mechanism responsible for an α-band fragmentation which is associated with the transition between Up and Down states (Fig 6).
Interestingly, in the medical context, the α-band is persistent in young subjects and becomes sparser with age [37]. Possibly, a higher neuronal activity (in younger subjects) leads to higher extracellular potassium which, in turn, increases the synaptic noise [38]. Another possible mechanism for fragmenting the α-band could involve the metabolic pathway, when the ATP coupled to the sodium concentration is decreased: during a burst, a high sodium concentration depletes ATP that deregulates the potassium current and thus leads to a phase of iso-electric suppression [11].

Up and Down states and the α-band
Neuronal networks exhibit collective transitions from Up to Down states [15,16,22]. This transitions from Up to Down states are responsible for the slow wave oscillations in the δ range, observed in deep sleep as well as anesthesia and probably originates from the thalamic network [16][17][18]. Previous modeling efforts have accounted for this phenomenon using single compartment neuron models with fast and slow currents responsible for neuronal spiking and slow modulation of the membrane voltage, respectively [39]. It was also shown that increasing hyperpolarization in a Hodgkin-Huxley type model, including calcium and hyperpolarizing currents [30] that was previously used to reproduce waxing and waning oscillations around 4Hz, could reduce the oscillation frequency [40]. Our approach here is thus in agreement with the previous theory even though we use a different modeling strategy. Indeed, we reported that the α-oscillation is only generated when the neuronal ensemble is in the Up state. We could not generate simultaneously in a single neuronal population both the α-oscillation and the Up-Down states transitions. Rather we needed a minimum of two coupled excitatory neuronal populations, one of them representing the thalamic neurons (Fig 5A U/D network) the function of which was to trigger the transitions between Up and Down states (Fig 6). Furthermore, we have shown that this neuronal population needs to account for hyperpolarization in order to trigger these Up-Down transitions (Fig 3) as also seen with classical Hodgkin-Huxley type models [30]. We reported here that the fraction of time spent in Up and Down states depends on the level of synaptic connectivity (Figs 3C, 3D and 6F-6I). However, by adding an inhibitory network, we were able to modulate the proportion of time spent in Up vs Down states by changing the input stimulation current (Figs 4C, 4D and 6B-6E) on this inhibitory population without varying the connectivity. This mechanism is feasible because the inhibitory input on both excitatory populations allows to destabilize the Up state and thus increase the transitions to Down state modulating the overall fraction of time spent in the Up state. To conclude, in the extreme case where the Down states are dominant, the overall voltage dynamics could resemble iso-electric suppressions, but this dynamics does necessarily rely on a metabolic stress [41,42].

How to interpret the α-band during general anesthesia?
During general anesthesia, the physiological mechanism leading to the emergence of the αband shortly after propofol injection remains unclear [7,43]. Possibly, during wakefulness, the amount of external stimuli suppresses the emergence of α-oscillations [2]. When the external stimuli ceases with propofol injection, the α-oscillation could become dominant (Fig 5B). The present model suggests that the initial state represents an already anesthetized brain where the neuronal networks do not receive any external stimuli, leading to the spontaneous emergence of the α-band.
General anesthesia needs to be sustained over the whole course of a surgery and thus controlling the optimal anesthetic injection to prevent cortical awareness or a too deep anesthesia remains a difficult problem [8,37]. Population models such as the one presented here could be used to test different activation pathways of anesthetic drugs. The present model accounts for the appearance an iso-electric suppression [44] in the EEG of the order of a few seconds ( Fig  7B and 7E) induced by increasing transiently the hypnotic. We predicted here that the fragmentation of the α-band results from a shift between Up and Down states dominance that could be tested with in vivo experiments. It would be interesting to further account for longer term consequences of an anesthetic input (several minutes). Indeed, the causality between αsuppressions and burst-suppressions [8] remains unexplained, suggesting that this relation could involve other mechanisms than the ones we modeled here based on synaptic plasticity, AHP and network connectivity.

Temporal fragmentation of an oscillatory band
We define here a measure for the persistence in time of the α-band. First we detect the peak spectral value S α (t) of the spectrogram as the highest power value in the extended range α min = 4 − α max = 16 Hz. When the condition S α (t) > T α is satisfied, we consider that the band is present and attribute a value x pr (t) = 1, otherwise x pr (t) = 0. When the time interval between 0 and T is divided into N bins at times t k , the persistence of the α band is defined by The persistence level P α measures the proportion of time where the α-band is present. To further quantify the fragmentation level, we introduce the disruption number D α that counts the number of times per minute where the peak spectral value S α (t) goes under the threshold T α .
x pr ðt k Þð1 À x pr ðt kþ1 ÞÞ: ð2Þ We call the fragmentation level the pair F α = (P α , D α ) (S6 Fig). For the human EEG data (Fig 1), we used a bin size w = 0.5s and a threshold value T α = 1.5dB. For the simulated data, we use the same bin size and T α = 10dB.

Modeling a single neuronal population based on short-term depressionfacilitation
For a sufficiently well connected ensemble of neurons, we use a mean-field system of equations to study bursting dynamics, AHP and the emergence of Up and Down states. This stochastic dynamical system consists of three equations [23,26,29] for the mean-field variable h, the depression y, and the synaptic facilitation x: where h + = max(h, 0) is the population mean firing rate [24]. The term Jxy reflects the combined synaptic short-term dynamics with the network activity. The second equation describes facilitation, and the third one depression. The parameter J accounts for the mean number of synaptic connections per neuron [23,45]. We previously distinguished [26] the parameters K and L which describe how the firing rate is transformed into synaptic events that are changing the duration and probability of vesicular release respectively. The time scales τ f and τ r define the recovery of a synapse from the network activity. We account for AHP with two features: 1) a new equilibrium state representing hyperpolarization after the peak response of the burst 2) two timescales for the medium and slow recovery to the resting membrane potential to describe the slow transient to the steady state. Finally, _ o is an additive Gaussian noise and σ its amplitude, representing fluctuations in the mean voltage.
In the case of a neuronal network that does not exhibit AHP the resting membrane potential is constant T • τ 0 = τ sAHP and T 0 = 0 in the subspace O sAHP ¼ fy < 1 1þLxðhÀ T 0 Þ and (Y AHP < y or h < H AHP )}, which represents the slow recovery to resting membrane potential.
The threshold parameters defining the three phases are Y h = 0.5, Y AHP = 0.85 and H AHP = −7.5. In this study, we varied the network connectivity parameter J 2 [5.6, 8.6] and all other parameters are described in S1 Table. To convert the mean-field variable h into a mean voltageh in mV, we use the following conversionh where V rest = −70 mV. We identified h min = −100 and h max = 1200 based on numerical simulations and chose A max = 200mV according to the classical amplitude of intracellular recordings.
defined piece-wise when it is present, as already discussed for model 3. The inhibitory population is always modeled without AHP and thus τ is constant and T = 0. All other parameters are described in the central columns called "2 populations" of S1 Table. deterministic dynamics oscillates around the point A Up at a frequency which corresponds to the dominant spectral band observed in Fig 2. The oscillation eigenfrequency ω Up further depends on the network connectivity J (Fig 2B and 2C), but not on the noise amplitude (Fig 2D and 2E). Note that the noise allows to generate persistent oscillations compared to the case of the pure deterministic system. Note that if we take τ = 0.025s, τ r = 0.5s and τ f = 0.3s, then the eigenvalues of A Up become Up state stability due to multiple re-entries in its basin of attraction. To explain the locking in the Up state, we recall that the stochastic trajectories starting inside the basin of attraction of the Up state can cross the separatrix Γ and fall into the Down state. However, because the deterministic vector field of system (3) is very shallow near Γ, the additive noise on the h variable can push the trajectories back into the Up state, where the field is stronger, and thus the trajectory is brought back into a neighborhood of A Up and continues to oscillate, as shown in S7(B) Fig (see inset).
To explain the other frequencies (than the eigenfrequency ω Up ) observed in the spectrum of h S7(C) Fig, we note that when a trajectory falls back in the Up state, it can produce a longer or shorter loop depending on its initial distance to the attractor A Up . These oscillations between the two basins of attraction define stochastic oscillations that contribute to the spectrogram of h.
Oscillations between Up and Down state in the phase-space. For a neuronal network with an AHP component, the Up state attractor has only real negative eigenvalues S8(C) Fig,  thus no oscillations are expected near the attractor. However, the presence of a slow AHP component S8 Fig (A-B pink) can push the dynamics into the Down state, as opposed to the case without AHP. Finally, in the Down state, the trajectories fluctuate with the noise until they escape. Once trajectories cross the separatrix Γ, they follow an almost deterministic path close to that of the unstable manifold of S S8 Fig (A-B grey) showing a long excursion in the phasespace before falling back near the attractor A Up . This dynamics explains the recurrent switches between Up and Down states.

S1 Table. Models (3) (1 population), (5) (2 populations) and (6) (3 populations) parameters.
For models (5) and (6), the inhibitory population is always without AHP and excitatory populations can be with or without AHP. For model (6) E 1 corresponds to the network with AHP (U/D), and E 2 to the network without AHP (α). The parameters values are chosen in agreement with [23,24,26,29,46]. The timescales τ f and τ r and the facilitation and depression rates K and L have been scaled in order to obtain oscillations in the range 5 − 15Hz around the Up state attractor A Up . The network connectivity values J ii (i = {I, E 1 , E 2 }) are adjusted (J = 4.21 in [29]) in order to reach a bifurcation that transforms the fixed point A Up from a saddle-point in [29] to an attractor. To obtain a bi-stable system, a minimal connectivity level is needed in the network. (PDF)

S1 Fig. Effect of network connectivity J and noise amplitude σ on model (3) without AHP.
A: Time-series and spectrograms of h (60s simulations) with SEF95 (blue curve) for J = 5.6 (upper), 6.6 (center) and 7.6 (lower). B: Mean value of the SEF95 for J 2 [3.8, 10]. C: Timeseries and spectrograms of h (60s simulations) with SEF95 (blue curve) for σ = 7 (upper), 11 (center) and 15 (lower). D: Mean value of the SEF95 for σ 2 [5,15]. Synaptic plasticity timescales: τ = 0.025s,τ r = 0.5s and τ f = 0.3s. A: 3D phase-space of the system with the two attractor points A Down (purple), A Up (red) and the saddle-point S (cyan) with its 2-dimensional stable manifold Γ (blue surface) which defines the separatrix. Stable trajectories (black curves) and unstable manifold of S (grey) and deterministic trajectories starting below (purple, resp. above light red) Γ falling to A Down (resp. A Up ). The phase-space is separated into 3 subspaces defining the different dynamics: fast O fast (above pink and orange meshes), medium O mAHP (below the orange mesh) and slow O sAHP (below the pink mesh).