## Abstract

We study a realistic model of a cortical column taking into account short-term plasticity between pyramidal cells and interneurons. The simulation of leaky integrate-and-fire neurons shows that low-frequency oscillations emerge spontaneously as a result of intrinsic network properties. These oscillations are composed of prolonged phases of high and low activity reminiscent of cortical up and down states, respectively. We simplify the description of the network activity by using a mean field approximation and reduce the system to two slow variables exhibiting some relaxation oscillations. We identify two types of slow oscillations. When the combination of dynamic synapses between pyramidal cells and those between interneurons accounts for the generation of these slow oscillations, the end of the up phase is characterized by asynchronous fluctuations of the membrane potentials. When the slow oscillations are mainly driven by the dynamic synapses between interneurons, the network exhibits fluctuations of membrane potentials, which are more synchronous at the end than at the beginning of the up phase. Additionally, finite size effect and slow synaptic currents can modify the irregularity and frequency, respectively, of these oscillations. Finally, we consider possible roles of a slow oscillatory input modeling long-range interactions in the brain. Spontaneous slow oscillations of local networks are modulated by the oscillatory input, which induces, notably, synchronization, subharmonic synchronization, and chaotic relaxation oscillations in the mean field approximation. In the case of forced oscillations, the slow population-averaged activity of leaky integrate-and-fire neurons can have both deterministic and stochastic temporal features. We discuss the possibility that long-range connectivity controls the emergence of slow sequential patterns in local populations due to the tendency of a cortical column to oscillate at low frequency.

## 1. Introduction

### 1.1. Ongoing Cortical Activity.

Realistic attractor neural networks composed of leaky integrate-and-fire (LIF) neurons store memories as stable equilibrium points of population-averaged activity (Amit & Brunel, 1997a, 1997b). In accordance with these models, it should be possible to observe a static state in the local field potential or electroencephalogram/magnetoencephalogram (EEG/MEG) of ongoing activity recorded from the cortex. However, it has been observed in many parts of the brain, such as the visual (Arieli, Sterkin, Grinvald, & Aertsen, 1996), auditory (Luczak, Bartho, & Harris, 2009) and somatosensory (Luczak et al., 2009) cortex, that the ongoing activity varies either with external stimuli or without them and that this variability cannot be averaged over trials. It is very likely that this variability in neural activity is not necessarily due to the noise such as synaptic noise (Holcman & Tsodyks, 2006; Parga & Abbott, 2007), quenched noise (Compte, Sanchez-Vives, McCormick, & Wang, 2003), and the finite size of the network (Destexhe, 2009) and that the population-averaged ongoing activity of the cortex may encode some sequential information.

Neural recordings are characterized by oscillations and irregular activity. Fast oscillations (30−80 Hz) have been associated with the binding of information (Singer & Gray, 1995; Buzsaki & Wang, 2012), whereas slower oscillations (<1 Hz) have been related to sleep and memory consolidation (Tononi & Cirelli, 2006; Diekelmann & Born, 2010). Slow wave sleep may cause the consolidation of memories (Marshall, Helgadottir, Molle, & Born, 2006) and the transfer of information from the hippocampus to the cortex (Squire, Stark, & Clark, 2004). This transfer could result from the replay of recent memory sequences (Wilson & McNaughton, 1994; Euston, Tatsuno, & McNaughton, 2007) that simultaneously occur in the cortex and hippocampus (Ji & Wilson, 2007). The process of consolidation during slow wave sleep might include long-term potentiation, notably during spindles nested in the up states (Rosanova & Ulrich, 2005), and long-term depression (Czarnecki, Birtoli, & Ulrich, 2007) inducing a reorganization of memories. Recent observations suggest, however, that low-frequency oscillations might also be related to cognition (Fries, Reynolds, Rorie, & Desimone, 2001; Petersen, Hahn, Mehta, Grinvald, & Sakmann, 2003; Gail, Brinksmeyer, & Eckhorn, 2004; Harris & Thiele, 2011).

### 1.2. Up and Down States.

Neurons in the cortex exhibit some characteristic up and down patterns, which are states of elevated and lowered activity (Steriade, Nunez, & Amzica, 1993; Contreras & Steriade, 1995; Metherate & Ashe, 1993), respectively, corresponding to phases of depolarization and hyperpolarization of the membrane potential or high and low concentrations of calcium (Ikegaya et al., 2004) in their intracellular media, respectively. Phases of high and low activity can also be detected based on the slow variation of multiunit activity (Sanchez-Vives et al., 2010) or spectral content of the local field potential (Mukovski, Chauvette, Timofeev, & Volgushev, 2007).

It has been shown that temporal patterns of up and down states can be observed at the scale of a local population of neurons (Sanchez-Vives & McCormick, 2000; Chauvette, Volgushev, & Timofeev, 2010), in addition to their presence in recordings of individual cells and electroencephalograms (Steriade et al., 1993). The slow oscillations can originate from the cortex only (Steriade et al., 1993; Sirota & Buzsaki, 2005). The onsets of up and down states in a local population are almost synchronous (Haider, Duque, Hasenstaub, & McCormick, 2006; Shu, Hasenstaub, & McCormick, 2003; Luczak et al., 2009). It has been measured that these states possibly result from a balance of excitation and inhibition (Shu et al., 2003), although it is unclear if this balance is dominated by excitation (Haider et al., 2006) or inhibition (Rudolph, Pospischil, Timofeev, & Destexhe, 2007; Shu et al., 2003). Some experiments have shown that regular spiking neurons maintain an elevated firing rate during the up state, whereas the activity of fast spiking neurons decreases more rapidly (Haider et al., 2006). On the other hand, other experiments have reported longer-lasting activity of inhibitory neurons rather than excitatory ones (Shu et al., 2003). In the latter case, the activity of inhibitory neurons could trigger the end of the up state (Shu et al., 2003; Brunel & Wang, 2001). Moreover, it has been shown that the synchrony of membrane potentials at the transition to the down state increases with the activity of interneurons (Chen, Chauvette, Skorheim, Timofeev, & Bazhenov, 2012), which may account for the greater synchrony observed at offset rather than onset of the up state (Chen et al., 2012; Volgushev, Chauvette, Mukovski, & Timofeev, 2006). Up states may thus share some similarities with persistent activity recorded during delayed-response tasks (Compte, Constantinidis, et al., 2003). These observations suggest that the mechanisms that generate these states result from intrinsic properties of local populations and the coordinated activity of their constituting neurons.

Previous studies argued that up and down states can be modeled by oscillations between two stable states of local populations. The mechanisms of transition from up to down state might include activity-dependent K^{+} currents (Compte, Constantinidis, et al., 2003), synaptic depression (Bazhenov, Timofeev, Steriade, & Sejnowski, 2002) and noise (Holcman & Tsodyks, 2006; Parga & Abbott, 2007). The increase in input resistance (Rigas & Castro-Alamancos, 2009) and reduction of release probability of excitatory postsynaptic potentials (Crochet, Chauvette, Boucetta, & Timofeev, 2005) during the up state argue that synaptic depression occurs between pyramidal cells before the termination of the up state. The transition from down to up state might result from miniature postsynaptic potentials (Timofeev, Grenier, Bazhenov, Sejnowski, & Steriade, 2000; Bazhenov et al., 2002) and the activity of large pyramidal cells in cortical layer V (Sanchez-Vives & McCormick, 2000; Chauvette et al., 2010).

The durations of up and down states depend on the experimental conditions. Usually down states are longer, up states shorter, and low-frequency oscillations more irregular in vitro (Sanchez-Vives & McCormick, 2000; Timofeev et al., 2000) than in vivo (Volgushev et al., 2006). The up and down states also tend to be more irregular during natural sleep (Okun, Naim, & Lampl, 2010; Luczak et al., 2009) than anesthesia (Volgushev et al., 2006; Crochet et al., 2005; Lampl, Reichova, & Ferster, 1999). This irregularity has been modeled by random transitions between the up and down states, between behaviorally activated states and the baseline activity (Kenet, Bibitchkov, Tsodyks, Grinvald, & Arieli, 2003; Sirota & Buzsaki, 2005).

### 1.3. Low-Frequency Oscillations.

Temporal patterns of up and down states are correlated with different types of oscillations, such as slow oscillations (Steriade et al., 1993; Sirota & Buzsaki, 2005) and sharp wave oscillations (Battaglia, Sutherland, & McNaughton, 2004). The classical view is that the slow oscillations depend mainly on the sleep cycle (Steriade & McCarley, 2005) and is a state of global synchrony (Volgushev et al., 2006).

Interestingly, ongoing activity correlated with slow wave signals depends on the state of arousal (Kisley & Gerstein, 1999) and the amount of sleep deprivation (Born, Rasch, & Gais, 2006; Vyazovskiy, Cirelli, Pfister-Genskow, Faraguna, & Tononi, 2008). The frequency of slow cortical and hippocampal oscillations can indeed be regulated by neuromodulatory substances related to attention (Salinas & Sejnowski, 2001) such as acetylcholine. Moreover, the amount of low-frequency oscillation recorded during sleep at the proximity of a local network depends on the prior use of this network (Born et al., 2006; Peyrache, Khamassi, Benchenane, Wiener, & Battaglia, 2009). The correlation between up and down states and oscillations may indeed be an important aspect of the offline replay of temporal patterns in the hippocampus (Battaglia et al. 2004, Ji & Wilson, 2007).

Moreover, it has been reported that the low-frequency power of local field potential spike-triggered average (LFP STA) decreases when attention is directed outside the receptive field of neurons near the electrode (Harris & Thiele, 2011; Fries et al., 2001).^{1} These observations suggest that the modulation of low-frequency oscillations during behavior might be local (Harris & Thiele, 2011; Fries et al., 2001; Mohajerani, McVea, Fingas, & Murphy, 2010; Lampl et al., 1999). It has been shown that low-frequency oscillations are modulated during activity (Gail et al., 2004; Petersen et al., 2003; Poulet & Petersen, 2008; Mohajerani et al., 2010) and representational switching (Ozaki et al., 2012). Moreover, down states similar to the ones recorded during sleep can be observed locally after sleep deprivation (Vyazovskiy et al., 2011). These experiments suggest that low-frequency oscillations might modulate neural activity and behavior.

It is known that up states propagate via local synaptic connections in vitro (Sanchez-Vives & McCormick, 2000; Compte, Sanchez-Vives, et al., 2003). If low-frequency oscillations are local, they might interact with each other through intracortical far-reaching axons (Buzas, Eysel, Adorjan, & Kisvarday, 2001), thalamocortical (Crunelli & Hughes, 2010), and hippocampal-cortical interactions (Battaglia et al., 2004; Ji & Wilson, 2007). It has indeed been reported that the modulation of low-frequency oscillations may modify the spatial distribution of local and long-range correlations (Smith & Kohn, 2008; Nauhaus, Busse, Carandini, & Ringach, 2009).

### 1.4. Dynamic Synapses.

Dynamic synapses (short-term plasticity; Markram, Wang, & Tsodyks, 1998), which may have an important role in the mechanism that generates up and down states (Bazhenov et al., 2002; Holcman & Tsodyks, 2006), may also participate in the encoding of memories in the working memory (Mongillo, Barak, & Tsodyks, 2008).

It has been shown that without dynamic synapses, networks of LIF neurons act as a low-pass filter for temporal variations of their external input (Fourcaud & Brunel, 2002; Fourcaud-Trocme & Brunel, 2005). On the other hand, with dynamic synapses, information encoded by the weights of instantaneously depressing synaptic connections is maximal for slow input waves (Fuhrmann, Segev, Markram, & Tsodyks, 2002), suggesting that encoding with dynamic synapses is optimal when operating at these frequencies.

Although inhibitory connections might also play an important role in shaping the up and down states (Sanchez-Vives et al., 2010; Compte, Reig, & Sanchez-Vives, 2009) and possibly synchronize neural activity during slow oscillations (Chen et al., 2012), the effect of short-term plasticity between interneurons (Gupta, Wang, & Markram, 2000) has not been as thoroughly studied as dynamic synapses between pyramidal cells (Wang et al., 2006).

### 1.5. Synchrony and Temporal Coding.

Previous models assume that the durations of down states are given by the time to recover from the up states, through some stochastic mechanism (Holcman & Tsodyks, 2006; Compte, Sanchez-Vives, et al., 2003; Bazhenov et al., 2002; Parga & Abbott, 2007). In such cases, no information is encoded by the time of transition between down and up state.

In the following, we study networks composed of leaky integrate-and-fire (LIF) neurons including recurrent dynamic synapses between pyramidal cells and between interneurons. We show that up and down phases, reminiscent of up and down states, respectively, emerge from intrinsic properties of the network. We develop a mean field model with the aim of explaining the slow dynamics of the network. In the mean field approximation, we demonstrate that the durations of up and down phases are deterministic and controlled by slow mechanisms of short-term plasticity. We distinguish two types of slow oscillations. The first type involves a competition between the slow dynamics of recurrent excitatory connections and that of recurrent inhibitory connections. It is characterized by asynchronous fluctuations of the membrane potentials at the end of the up phase. The second type involves mainly the slow short-term plasticity between interneurons. In this case, the activity of inhibitory neurons decreases more slowly than the activity of pyramidal cells, which induces a rapid change in the balance between excitation and inhibition in the network. Moreover, the fluctuations of membrane potentials are more synchronous at the end than at the beginning of the up phase. The simulation of the network of LIF neurons shows that the durations of up and down phases have both deterministic and stochastic contributions resulting notably from finite size effect. We study the role of finite size effect and slow currents on the duration of the up and down phases.

Moreover, we simulate the modulation of up and down phases by some low-frequency oscillatory input representing long-range intracortical, thalamocortical, or hippocampal-cortical interactions. Slow deterministic dynamics of neuronal population-averaged activity that takes the form of chaotic trajectories in the mean field limit can emerge from the interaction between low-frequency oscillations and dynamic synapses. Dynamic synapses induce relaxation oscillations, which, when modulated by the low-frequency oscillations, can take the form of chaotic relaxation oscillations. These oscillations exhibit some deterministic temporal patterns that encode sequential information by using the timing of transitions between up and down phases.

## 2. Model and Methods

### 2.1. Network Model.

We consider a local population of neurons (or a cortical column) composed of 800 pyramidal cells and 200 interneurons. The activity of neurons is simulated using the LIF model (see, e.g., Ricciardi, 1977; Brette et al., 2007). Note that the network is fully connected. Facilitation-dominant excitatory synaptic connections of type E1 (Wang et al., 2006), the most common type of excitatory connections in the prefrontal cortex (Wang et al., 2006), connect pyramidal cells. Depression-dominant inhibitory synaptic connections of type F2 (Gupta et al., 2000), the most common type of inhibitory connections in the neocortex (Gupta et al., 2000), connect interneurons (see Figure 1). Connections from excitatory to inhibitory neurons and vice versa are assumed to be static for simplicity even though they also exhibit short-term plasticity in the cortex (Gupta et al., 2000; Markram et al., 2004; Wang et al., 2006). A low-frequency oscillatory input is applied to both pyramidal cells and interneurons in section 3.6.

### 2.2. Neuronal Dynamics.

*V*of a neuron

_{i}*i*, where , is described by the following equations (Barbieri & Brunel, 2007): where depends on the type of neuron

*i*: if

*i*is a pyramidal cell (given as ) or if

*i*is an interneuron (). is the time constant of the membrane potential;

*I*, NMDA receptor-mediated currents;

^{N}_{i}*I*, AMPA receptor-mediated currents;

^{A}_{i}*I*, GABA

^{G}_{i}_{A}receptor-mediated currents; and

*I*, the external input current. We do not consider metabotropic GABA

^{ext}_{i}_{B}receptor-mediated currents in this letter. When the membrane potential

*V*reaches the threshold , a spike is emitted and the membrane potential is set to for a duration , which is the absolute refractory period.

_{i}*I*is assumed to be gaussian and -correlated as follows: where is the mean input; , the standard deviation of the gaussian distribution; and , uncorrelated white noise with a mean over their distribution and correlations . Finally, is the Kronecker delta.

^{ext}_{i}*J*(

_{ij}*t*), the synaptic weight from neuron to

*i*;

*t*, the time of a spike

^{k}_{j}*k*emitted from neuron

*j*; and

*d*, the synaptic delay. The synaptic delay is distributed uniformely in [0,

*d*

_{max}] where

*d*

_{max}is the maximal synaptic delay.

*I*and

*U*denote decaying and rising currents, respectively.

*J*(

_{ij}*t*) of type E1 and F2 are described by the following equation (Tsodyks, Pawelzik, & Markram, 1998; Mongillo et al., 2008): where

*J*is of type E1 if or F2 if ;

_{ij}*x*is the amount of available resources (or ratio of readily releasable vesicles);

_{ij}*u*, the utilization factor (proportional to the calcium concentration in the intracellular media); and

_{ij}*J*

^{0}

_{ij}, the maximal weight value. For or ,

*x*=

_{ij}*u*=1. The dynamics of

_{ij}*x*and

_{ij}*u*are described by the following equations (Tsodyks et al., 1998; Mongillo et al., 2008): where is the timescale of synaptic depression; , the timescale of synaptic facilitation ( if and if ); and , the parameter of synaptic facilitation. In numerical simulations,

_{ij}*u*(

_{ij}*t*) and

*x*(

_{ij}*t*) are calculated just before a spike.

*A*is the amplitude of low-frequency oscillations; , the ratio of the amplitude of low-frequency oscillations applied to inhibitory neurons and that applied to excitatory neurons (

*r*=1); , the instantaneous phase; , the initial phase; , the angular velocity; and , the period of low-frequency oscillatory input.

^{A}_{E}To simulate the network of LIF neurons, we use a basic clock-driven algorithm and store spike events in a circular matrix in order to consider the delay of synaptic transmission (Brette et al., 2007).

### 2.3. Mean Field Approximation.

In the following, we consider population-averaged quantities . We do not consider the transmission delay *d* for simplicity.^{2} Population-averaged quantities are sometimes simply given as , where is the index of a population, when . We first consider the case under which there is no low-frequency oscillatory input: .

To analyze the population-averaged dynamics of the network, we consider a mean field approximation (Amari, 1971; Amit & Tsodyks, 1991; Amit & Brunel, 1997a, 1997b; Tsodyks et al., 1998; Romani, Amit, & Mongillo, 2006). As in Romani et al. (2006), we assume that the correlations due to short-term plasticity in the current can be omitted. Synaptic weights are assumed to scale as such that the fluctuations in the recurrent currents are of the order and can be neglected when *N* is large enough (Barbieri & Brunel, 2007).

_{A}receptors can be described by the following equations:

^{4}we can write the following equations to describe the average steady-state membrane potential of the population (Amit & Tsodyks, 1991; Amit & Brunel, 1997b):

*p*(

*t*). Then, and , where is the Laplace transform of

*p*(

*t*), are given as follows (Alili, Patie, & Pedersen, 2005; Romani et al., 2006): where

*H*(

_{l}*z*) are the Hermite functions. For a numerical simulation, we use hypergeometric functions to evaluate the Hermite functions as follows (Lebedev & Silverman, 1972): where is the gamma function and

_{1}

*F*

_{1}is the hypergeometric function.

### 2.4. Critical Manifold.

The mean field approximation described in section 2.3 does not provide an instantaneous description of the population-averaged firing rate. To overcome this limitation, we exploit the fact that the variables of the system evolve with different timescales in order to apply fast-slow decomposition.

Networks with facilitation-dominant excitatory synapses are known to exhibit a bistable regime that can model the persistent activity during delayed-response tasks (Romani et al., 2006; Holcman & Tsodyks, 2006; Barbieri & Brunel, 2007; Mongillo et al., 2008). In fact, networks without dynamic synapses can also have a bistable regime (Amit & Brunel, 1997b). However, fast depression of recurrent excitatory synapses allows an increase in the parameter space region with bistability by limiting the increase in the firing rate of pyramidal cells. Because up states are very similar to persistent activity, we propose that the bistable regime of networks with fast depression in excitatory synapses can model up phases (such as in Bazhenov et al., 2002; Holcman & Tsodyks, 2006). Therefore, we group the variables representing the evolution of the population-averaged membrane potential , , with the relatively fast variables of the dynamic synapses *x _{EE}* and

*u*. These four variables represent the fast subsystem. We consider the case under which there is bistability in the fast system. The two other variables,

_{II}*u*and

_{EE}*x*, constitute the slow subsystem.

_{II}*I*, in which it is possible to use the implicit function theorem, we can solve the system

*f*(

*a*,

*b*)=0 and consider the following function: The dynamics of the slow system can then be studied independently of the fast one in each domain

*I*by considering the following system:

When , the system can jump from one sheet of the manifold to another. In the regions of the phase space located far from the jumps, the theory developed by Fenichel (1972, 1979) and others (Hirsch, Pugh, & Shub, 1977) ensures that if the critical manifold *M*_{0} is compact and normally hyperbolic, then there exists at a distance a “slow” manifold that has identical properties.^{5} We assume in the following that the critical manifold is normally hyperbolic.^{6} The fast dynamics converges rapidly to the stable regions of the slow manifold. The trajectory of the system then remains at proximity of the slow manifold until a jump occurs. We study the location of the critical manifold in the phase space, although the slow dynamics takes place on the -close slow manifold.

In the case of the system studied, the separation into fast and slow subsystems is motivated by the fact that and are much smaller than , , , and .

*T*is defined in equation 2.25 and the instantaneous variables and described by equations 2.28 and 2.29, respectively. These equations do not describe exactly the dynamics of the network (Romani et al., 2006). A more rigorous approach would require studying the nonstationary distribution of membrane potentials and variables of short-term plasticity. However, the calculation of such distribution is currently an unsolved problem (Lindner, 2009) and is out of the scope of this letter. The results of simulations using equations 2.41 and 2.42 show that this approximation allows for the description of slow oscillations (see section 3). Note that coupled networks of LIF neurons can typically exhibit fast oscillations due to delayed inhibition (Brunel & Hakim, 1999; Brunel & Wang, 2003; Geisler, Brunel, & Wang, 2005).

^{7}In summary, this approximation, referred to as the simple mean field approximation in the following, is described by the set of variables , and , and the set of equations 2.41, 2.42, 2.28, and 2.29.

*u*and

_{EE}*x*. By using equations 2.30 and 2.31 or 2.32 and 2.33 that describe the steady states of fast variables for the dynamic synapses, we can write the steady state of the fast variables and as a function of the steady-state membrane potential , where , and slow variable

_{II}*u*as follows:

_{EE}^{8}

*f*and

_{x}*f*are not detailed here for conciseness and can be obtained using simple algebraic operations. We can then write the definition of the critical manifold as follows. The critical manifold is the set of variables , such that the following condition holds: where and are defined as follows (see equations 2.27, 2.43 and 2.46, and 2.47):

_{u}Equations 2.48 and 2.49 are solved using a numerical solver (Levenberg Marquardt algorithm (Marquardt, 1963)) from initial guesses of the solutions taken on a grid. This method allows obtaining a set of discrete points on the critical manifold. In order to obtain approximate expressions of equations 2.52 and 2.53 in the entire domain of the definition of *u _{EE}* and

*x*, that is, [0, 1]

_{II}^{2}, the discrete points obtained are interpolated by locally weighted scatterplot smoothing (Cleveland, 1979). In order to increase the quality of the interpolation, the set of points on the fold curves is also included in the interpolation (see section 2.5).

In appendix A, we show the stable regions of the critical manifold. The stability is determined by calculating the linearization of the four equations describing the fast system and the eigenvalues of the corresponding Jacobian matrix (see equations 2.41, 2.42, 2.28, and 2.29).^{9} Moreover, we verify that *x _{EE}* and

*u*rapidly reach their steady state and , respectively. In the following, we simplify the analysis of stability by considering only the time evolution of the population-averaged membrane potential (Tuckwell, 1988; Amit & Brunel, 1997b) and find out the corresponding linearization (the Jacobian matrix) to determine local stability. The stability condition derived is thus related to only the fastest variable . The Jacobian matrix on which the stability of the fastest system depends is given in appendix D. The stability analysis achieved using all four fast variables and only these two fast variables gives identical results (see Figures 12 and 2, respectively).

_{II}### 2.5. Fold Curves.

*u*and

_{EE}*x*are varied. To do so, we use the Jacobian matrix obtained in appendix D and derive expressions for the eigenvalues and of the Jacobian matrix as functions of the variables , where (see appendix E). Using a numerical solver with initial guesses of taken on a grid, we solve the following systems in order to determine the set of points for which bifurcations occur (note that we do not consider here Bogdanov-Takens bifurcations):

_{II}The set of fold (saddle-node) bifurcation points is called a fold curve. Fold curves are approximated by interpolating between the set of points obtained from the numerical solver.

### 2.6. Slow System.

*I*for which equations 2.48 and 2.49 can be solved:

*I*, the system of equations 2.56 and 2.57 is analyzed using classical methods of nonlinear analysis to determine the steady states of the slow dynamics defined as A numerical solver is used to determine the steady states in each domain

*I*. Moreover, we calculate the Jacobian matrix corresponding to the slow dynamics in order to determine the stability of these steady states (not shown here).

## 3. Results

### 3.1. Critical Manifold.

The simulation parameters are given in appendix B if not stated otherwise. Our initial consideration is that there is no low-frequency oscillatory input. The network simulated with the parameters given in appendix B exhibits some slow oscillatory dynamics resulting from its intrinsic properties.

The critical manifold obtained by the method described in the previous section is shown in Figures 2A1, 2A2, and 2A3 using equations 2.30 and 2.31. The critical manifold obtained using equations 2.32 and 2.33 is shown in Figures 2B1, 2B2, and 2B3. We plot the values of the firing rate for pyramidal cells and interneurons as functions of the slow variables *u _{EE}* and

*x*in the domains to which the implicit function theorem is applicable. The manifold is composed of three sheets: two stable sheets and an unstable one. The first stable sheet is located at the lower left corner in the definition domain of

_{II}*u*and

_{EE}*x*([0, 1]

_{II}^{2}). In this part of the critical manifold, the firing rate of pyramidal cells varies from approximately 0 to 2 Hz. This part of the manifold can be interpreted as a phase of low activity. Moreover, the firing rate at the second stable sheet, located at the upper right corner, varies from 5 to 40 Hz. This part of the manifold can be interpreted as a phase of high activity. These two stable sheets are delimited and connected to the unstable manifold via two fold curves. The first fold curve given as

*FC*

_{1}in Figures 2A3 and 2B3 separates the stable upper sheet from the unstable one. The slow dynamics of the system may jump from the upper sheet (up phase) to the lower sheet (down phase) at this curve. The second fold curve given as

*FC*

_{2}separates the stable lower sheet from the unstable one. The slow dynamics of the system may jump from the lower sheet (down phase) to the upper sheet (up phase) at this curve. Moreover, a trajectory obtained by using the simple mean field approximation and simulating the networks of LIF neurons is plotted in Figures 2A1 to 2A3 and 2B1 to 2B3, respectively. The trajectory of the simple mean field approximation and simulation of LIF neurons are at the proximity of the critical manifold. The fast-slow decomposition described in section 2.4 is thus a good approximation. We will detail in section 3.3 the origin of differences between the critical manifold and the simulated trajectories.

Between the fold curves *FC*_{1} and *FC*_{2}, there is a region of bistability where the system can be in either the up or down phase (see Figures 2A3 and 2B3). The size of this bistable region depends on the parameters of the network, particularly the synaptic weights (see Figure 3).

For each domain delimited by the fold curves given as *I*_{1}, *I*_{2} and *I*_{3} in Figures 2A3 and 2B3, the implicit function theorem can be used to determine the surfaces , described in equations 2.52 and 2.53 for each sheet of the critical manifold. *I*_{1}, situated to the left of curve *FC*_{2}, is the domain of the lower sheet. Similarly, *I*_{3} situated to the right of curve *FC*_{1} is the domain of the higher sheet. There are no Hopf bifurcations for the parameter values considered here.

### 3.2. Simple Mean Field Approximation.

The slow dynamics of the simple mean field approximation is studied separately in domains *I*_{1} and *I*_{3}, corresponding to the slow dynamics of the down and up phase, respectively. Figures 3A1, 3B1, and 3C1 show the nullclines of the system described in equations 2.56 and 2.57, whose intersections give the steady states of the slow system (see equations 2.58 and 2.59). A pair of nullclines can exist per stable sheet of the critical manifold (two sets in total). These nullclines are defined only within the domain in which the corresponding sheets of the critical manifold are defined. We describe in Figure 3 three cases corresponding to three distinct set of parameters.

The simple mean field approximation is simulated by solving the differential equations in equations 2.41 and 2.42 for the membrane potential and in equations 2.28 and 2.29 for the dynamic synapses, by the Euler method.

The parameters used for the first case shown in Figures 3A1 to 3A3 are identical to the ones used in Figure 2. The slow system exhibits no stable equilibrium point (see Figure 3A1) on either the upper or lower sheet of the critical manifold. Because the slow variables are bounded within the domain [0, 1]^{2}, the slow motion is necessarily dynamical like periodic, quasi-periodic, or chaotic. We show in Figures 3A2 and 3A3 that the slow dynamics is periodic. The trajectory remains on the upper sheet of the critical manifold (up phase) for s and on the lower sheet (down phase) for s. The network thus exhibits oscillations with a frequency Hz due to its intrinsic properties.

The second case (see Figures 3B1 to 3B3) corresponds to an increase in strength of recurrent excitatory connections *J*^{0}_{EE}. Similar to the first case, there is no equilibrium for the slow system, and the trajectory is periodic. Stronger recurrent connections *J*^{0}_{EE} increase the duration of the up phase (+0.04 s) and reduce the duration of the down phase (−0.56 s).

### 3.3. Network of LIF Neurons.

The network of LIF neurons is simulated by solving the differential equations 2.1 to 2.12 by the Euler method.

In Figure 4, we simulate the network of LIF neurons with the same three set of parameters used for Figure 3 in order to compare the results with the simple mean field approximation. In order to determine the critical manifold, we consider here the more rigorous approach and use equations 2.32 and 2.33. First, the steady states of the slow system obtained by the more rigorous method are similar to the simple mean field approximation. In Figures 4A1 and 4B1, there is no stable equilibrium for the slow system. In Figure 4C1, there is a unique stable equilibrium on the upper sheet of the critical manifold.

As predicted by the mean field approximation, we observe low-frequency oscillations of the population-averaged firing rate in Figures 4A2 and 4B2. The up and down phases can be detected based on a threshold of the firing rate. This threshold can be defined using the bimodal distribution of the logarithm of the firing rate. This technique has been used to study the duration of up and down states recorded in vitro (Sanchez-Vives et al., 2010). In Figure 4A2, the mean durations of up and down phases are s and s, respectively, with a 95% confidence interval. These durations are similar to those of up and down states observed in vitro (Sanchez-Vives & McCormick, 2000; Sanchez-Vives et al., 2010). Figures 4A3 and 4B3 show that the first and second cases, respectively, exhibit a bimodality. Consequently, up and down phases can be detected. There is no clear bimodal distribution for the third case (see Figure 4C3).

There are important differences between the mean field approximation and the dynamics of LIF neurons. First, the oscillations shown in Figures 4A2 and 4B2 are irregular compared to the corresponding oscillation shown in Figures 3A3 and 3B3. Moreover, although the mean field approximation of Figure 4C1 predicts that the trajectory of the system settles in the upper sheet of the critical manifold, the simulated activity shows some occasional transitions to the down phase (see Figure 4C2). Finally, the mean durations of up and down phases in Figure 4A2 are 0.75 s and 4.75 s, respectively, which are longer than the durations predicted by the simple mean field approximation (see section 3.2). These differences may result from the simplifications of the mean field approximation or random fluctuations perturbing the system. In the following, we consider the role of finite size effect and slow NMDA currents.

In the simulation of LIF neurons, we observe that the population-averaged firing rates and fluctuate at proximity of the critical manifold in Figures 2B1 and 2B2, respectively. When the trajectory is between the fold curves *FC*_{1} and *FC*_{2}, the system is bistable, and a rapid change of the firing rate can induce a transition from up to down phase or vice versa. The finite size effect can result in fluctuations and induce some random transitions between up and down phases. In order to show this, we study the influence of the network size on the mean durations of up and down phases in Figures 5A1 and 5A2, respectively. The size of the network clearly affects the mean duration of the down state (see Figure 5A1). Moreover, we compute the standard deviation of the instantaneous firing rate , defined as where *N _{s}*(

*t*) is the number of spikes fired in the network between the discrete time steps

*t*and

*t*+

*dt*. Figures 5B1, 5B2, 5C1, and 5C2 show that the standard deviation of the firing rate decreases with the network size. Finite size effect induces some fluctuations in the network at the first order

*n*

_{1}(

*t*) of the firing rate , (Brunel & Sergi, 1998).

Second, we have considered in equation 2.21 that the NMDA currents have slower dynamics than membrane potentials, although . The slow currents might in fact participate in the slow dynamics (Brunel & Sergi, 1998). To study the role of NDMA currents, we compute the mean durations of up and down phases for different proportions of NMDA currents. Figure 6B shows that NMDA currents induce a significant increase in the duration of up phases.

In summary, finite size effect and NMDA currents account for a decrease of the frequency of slow oscillations and fluctuations in the population-averaged dynamics.

### 3.4. Mechanism of Up and Down Phases.

Previous models describing oscillations resulting from short-term plasticity considered dynamic synapses between pyramidal cells (Bazhenov et al., 2002; Holcman & Tsodyks, 2006; Mongillo et al., 2008; Barak & Tsodyks, 2007) or from pyramidal cells to interneurons (Melamed, Barak, Silberberg, Markram, & Tsodyks, 2008). In networks with dynamic synapses between pyramidal cells, self-sustained oscillations of population spikes have been reported (Barak & Tsodyks, 2007; Mongillo et al., 2008). In this case, the time between two population spikes depends on the recovery from synaptic depression, and its order of magnitude is thus similar to the time constant of synaptic depression (Mongillo et al., 2008). Moreover, the mechanism terminating the population spike is the rapid decrease of the number of neurotransmitters *x _{EE}* (Mongillo et al., 2008). As suggested by the experimental results of Wang et al. (2006), the order of magnitude of observed in vivo would correspond to population spikes with a frequency close to theta oscillations (Mongillo et al., 2008).

Similarly to these previous models, the network analyzed in this letter includes recurrent dynamic synapses between pyramidal cells. In addition, we consider here the role of dynamic synapses between interneurons. Consequently, the mechanism generating oscillations depend on the dynamic synapses between interneurons. Figure 7 shows that the shapes of oscillations are modified by the short-term plasticity between interneurons.

Initially, the system is in the down phase between s and s (see Figure 7). The strength of recurrent excitatory connections *J _{EE}* increases slowly because of synaptic facilitation, whereas the strength of recurrent inhibitory connections

*J*is approximately constant. At a certain time, the stable equilibrium of the fast system disappears, which induces a rapid increase of the firing rate and . This corresponds to the onset of the up phase at s in Figure 7.

_{II}During the up phase, the decrease in firing rate of pyramidal cells is concurrent with a decrease in the strength of recurrent inhibitory connections *J _{II}*. The interneurons are thus less inhibited, and the firing rate decreases more slowly compared to the case without dynamic synapses between interneurons (cf. Mongillo et al., 2008). Therefore, the higher firing rate of pyramidal cells is balanced by an increase of inhibition resulting in a stable equilibrium for the fast system on the upper sheet of the critical manifold (see Figures 2A1, 2A2, 2B1, and 2B2). The trajectory then moves slowly at proximity of the critical manifold. The firing rates and decrease almost linearly with time (see Figure 7A1), and the amount of neurotransmitter

*x*recovers gradually.

_{EE}Further decreases of *J _{II}* and of

*J*, resulting from a slow increase of

_{EE}*u*and decrease of

_{EE}*x*, respectively, lead the fast subsystem to the fold curve on the upper sheet of the critical manifold. The fast systems then return to the lower sheet.

_{II}The balance between excitation and inhibition in the network is thus preserved during the down and up phases, as reported in vivo and in vitro (Haider et al., 2006; Shu et al., 2003; Rudolph et al., 2007). This dynamic balance allows an increase in durations of the up and down phases compared to networks including solely dynamic synapses between pyramidal cells (Mongillo et al., 2008; Barak & Tsodyks, 2007). Although the order of magnitude of synaptic depression between pyramidal cells is approximately 0.03 s, the frequency of oscillations ( Hz) described in this letter corresponds to slow oscillations.

### 3.5. Two Types of Slow Oscillations.

We now consider the effect of changing the balance in favor of inhibition by reducing the time constant of short-term facilitation between pyramidal cells . In Figure 8A, we show that slow oscillations can also occur when s, that is, when *u _{EE}* is considered as a fast variable and

*x*is the only slow variable. In this case, Figure 8B2 shows that the end of the up phase is characterized by synchronous fluctuations of the membrane potentials (; see appendix F for a definition of the measure of synchrony used), whereas the beginning of the up phase is relatively less synchronous (). Note that the firing rate is greater at the beginning than at the end of the up phase; Hz and Hz, respectively. This characteristic is reminiscent of the in vivo observation that the transition to the down state can be more synchronous than the transition to the up state (Volgushev et al., 2006; Chen et al., 2012).

_{II}^{11}In the previous case ( s), the membrane potential fluctuations at the end of up phases were more asynchronous (see Figure 8B1).

Moreover, Figures 8C1 and 8C2 show that when s, the activity of pyramidal cells and that of interneurons decreases at the same moment, whereas when s, the inhibitory neurons remain active longer during the offset of the up phase. This induces a change in the balance between excitation and inhibition in favor of inhibition, which results in turn in a rapid termination of the up phase. Finally, Figure 8B3 shows that when all variables of short-term plasticity except *x _{II}* are approximated by their steady-state values (see equations 2.30 and 2.31), the simple mean field approximation exhibits some slow oscillations similar to Figure 8B2, which confirms that only the slow variable

*x*(

_{II}*t*) drives the slow oscillations in this case.

Note that some spikes are emitted during the down phase and that the distribution of membrane potentials does not display a clear bimodality (Bazhenov et al., 2002; Volgushev et al., 2006) although we were able to detect up and down phases with the technique described in Sanchez-Vives et al. (2010; see the discussion in section 4.2).

### 3.6. Effect of Low-Frequency Input.

We now consider the effect of a low-frequency oscillatory input. It is known that relaxation oscillations, which are forced by a periodic signal, such as the forced van der Pol oscillator (Guckenheimer & Holmes, 1983) or forced FitzHugh-Nagumo model (Fitzhugh, 1966; Izhikevich, 2007), can exhibit phase locking through the forcing signal as well as chaotic dynamics.

When the simple mean field approximation is simulated with an oscillatory input of frequency *f _{in}*, the resulting dynamics of the simple mean field approximation exhibits some complex deterministic trajectories (see Figure 9).

In order to study the frequency of the relaxation oscillations in this case, we consider a Poincaré section between the upper and lower stable sheet of the critical manifold and a return map on this section. The time interval between the moment the trajectory leaves the Poincaré section and returns to it from the down state (i.e., the time interval between two successive transitions to the up phase from the down phase) is given as . In Figure 9A, is plotted as a function of the frequency of the oscillatory input noted *f _{in}*. We observe that the slow oscillations produced by the network synchronize with the oscillatory input and induce subharmonic synchronization. Moreover, between the domains of subharmonic synchronization, the trajectory of the network can become chaotic. Figure 9B shows that bifurcates to chaos through a succession of period-doubling bifurcations with

*f*as a bifurcation parameter. Different modes of forced oscillations are shown in Figures 9C1 to 9C3 and 9D1 to 9D3.

_{in}When the network of LIF neurons is simulated with an oscillatory input of frequency *f _{in}*, the resulting dynamics exhibits some complex trajectories (see Figure 10) similar to the ones observed using the mean field approximation (see Figure 9).

In the following, we call the interburst interval the duration between the onsets of two consecutive up phases. Moreover, we represent the duration of the up phase . In Figure 10A, similar to the case described in Figure 9A, it appears that the slow oscillations of the network synchronize with the oscillatory input at different modes of synchronization. However, some fluctuations persist owing to noise, which causes the trajectory to fluctuate. Consequently, period-doubling bifurcations that lead to chaotic trajectories cannot be well detected by the simulation of LIF neurons in Figure 10A whereas it is visible in Figure 9B. Because the frequency range in which the period-doubling bifurcations occur is very narrow in Figure 9B, it is difficult to observe period-doubling bifurcations when simulating LIF neurons. Nonetheless, the distribution of interburst intervals depends on the frequency of oscillatory input in the simulation of LIF neurons. The bimodal distribution of the interburst intervals shown in Figure 10A is reminiscent of the chaotic attractor of Figure 9B, which exhibits two islands at Hz. Note that bifurcation diagrams obtained using the mean field approximation and the simulation of LIF neurons are different notably because of the finite size effect and NMDA currents as shown in section 3.3. Figure 10B shows that the distribution of the up phase duration becomes multimodal for oscillatory frequencies Hz. The network of LIF neurons exhibits various types of synchronization, as shown in Figures 10C1 to 10C3 and 10D1 to 10D3.

In Figure 10D2, the trajectory of population-averaged activity occupies the phase space and is not close to a periodic trajectory. This spreading of the dynamics to larger regions of the phase space is partly due to deterministic dynamics. The timing of transitions between up and down phases depends here on the contributions of deterministic and stochastic components.

In order to verify whether two successive interburst intervals, represented as and , are not purely independent random variables but reflect determinism to some extent, the mutual information (Peng, Long, & Ding, 2005) between interburst intervals and their first and second successors, and , respectively, is calculated to quantify their dependency. Higher mutual information of two consecutive interburst intervals indicates some deterministic contribution to the dynamics. Figure 11A shows a deterministic contribution to the interburst intervals, notably for Hz. Moreover, Figure 11B shows that there is also a deterministic contribution to the durations of up phases. These observations suggest that interburst intervals and up phase durations encode temporal information when the oscillatory input modulates the slow system.

The mutual information and trajectories show that temporal information can be encoded in the population-averaged quantities of the network, namely, the timing and duration of up phases, although the activities of single neurons are asynchronous and subject to large fluctuations. This network can therefore encode some temporal information in the imprecise but deterministic timing of population spikes.

## 4. Conclusion

### 4.1. Summary.

We have shown that short-term plasticity of type E1 between pyramidal cells and of type F2 between interneurons can result in slow (<1 Hz) relaxation oscillations at the population level. The mean approximation introduced in this letter allows simplifying the dynamics of the network to two slow variables *u _{EE}* and

*x*. We have shown that the dynamics of these slow variables is sufficient to explain the occurrence of the slow oscillations.

_{II}We have identified two types of slow oscillations according to the timescale of short-term facilitation between pyramidal cells. When both *u _{EE}* and

*x*participate in the slow dynamics, the balance between the activity of pyramidal cells and that of interneurons is maintained and the membrane potential fluctuations are asynchronous during the end of the up phase. This case is reminiscent of the experimental observations described in Haider et al. (2006). When

_{II}*x*is the only slow variable, slow oscillations are characterized by synchronous membrane potential fluctuations at the end of the up phase. Moreover, the decrease in inhibition between interneurons is slower than the decrease in excitation between pyramidal cells and, consequently, the activity of interneurons dominates over excitation. These characteristics have been described in vivo and in vitro (see Shu et al., 2003, and Volgushev et al., 2006), but were not fully explained by previous models (Bazhenov et al., 2002; Compte et al., 2003; Holcman & Tsodyks, 2006; Chen et al., 2012). Moreover, the two types of slow oscillations described in this letter may be related to the two types recently observed in Luczak and Bartho (2012). Further research is needed in order to verify if short-term plasticity between interneurons can account for the discrepancy between experimental observations described in Haider et al. (2006) and Shu et al. (2003). Note that slow GABA

_{II}_{B}currents may modify the shape of transitions to the down state (Mann, Kohl, & Paulsen, 2009).

We have studied the differences between the result of our mean field approximation and simulation of LIF neurons. Finite size effect induces some fluctuations in the firing rate. These fluctuations can cause some random transitions between up and down phases. In smaller networks, the durations of up and down phases become irregular, and down phases shorten. We have shown that slow NMDA currents can increase the duration of up phases.

Finally, low-frequency oscillations modulate these relaxation oscillations, resulting in synchronization and eventually chaos. In the simulation of a network of LIF neurons, for low-frequency oscillations around 0.5 Hz, the dynamics has both a stochastic and deterministic contribution, which can encode temporal features and sequential memory. The calculation of the mutual information between successive interburst intervals shows that the timing of up phases encodes some sequential information.

### 4.2. Bimodality of the Membrane Potential.

We have not observed a bimodal distribution of membrane potentials (Bazhenov et al., 2002; Volgushev et al., 2006) in our simulations. Note that it has been reported that slow oscillations can occur without two distinct states for the membrane potential (Lampl et al., 1999). The absence of bimodality may result from the choice of leaky integrate-and-fire neurons, which simplifies many aspects of spike generation. A simulation involving a more realistic conductance model might exhibit clear bimodality.

However, several properties of the slow oscillations described in this letter argue in favor of a different interpretation. First, it has been reported that very few neurons fire during the down state, whereas the mechanism of upward transition detailed in section 3.4 relies on short-term plasticity and spiking activity during the down phase. Second, the local network described in this letter can generate very low-frequency (smaller than 0.1 Hz) oscillations (data not shown). The current model may thus also account for the generation of slow oscillations during the awake state (Gail et al., 2004; Petersen et al., 2003; Poulet & Petersen, 2008; Mohajerani et al., 2010) or the default mode network (Greicius, Krasnow, Reiss, & Menon, 2003). These oscillations might occur in local populations (Harris & Thiele, 2011; Fries et al., 2001; Mohajerani et al., 2010; Lampl et al., 1999). Consequently, these local oscillations might interact via long-range connections (Buzas et al., 2001; Smith & Kohn, 2008; Nauhaus et al., 2009). This study shows that local populations of neurons, or cell assemblies, can oscillate at low frequency.

### 4.3. Slow Oscillation Code Hypothesis.

It has been proposed that the activation of up and down states may be the result of a synchronous propagation of activity with millisecond precision through synfire chains (Ikegaya et al., 2004; Abeles et al., 1993; Kumar, Rotter, & Aertsen, 2010). However, the existence of precisely timed sequential patterns of up and down states, which are the signature of synfire chains, is controversial (Baker & Lemon, 2000; McLelland & Paulsen, 2007; Mokeichev et al., 2007; Roxin, Hakim, & Brunel, 2008). The onset of up state might nonetheless be characterized by some stereotypical sequence of activity (Luczak et al., 2009).

Our analysis of a realistic local cortical network suggests that slow patterns of up and down phases, when modulated by long-range oscillatory inputs corresponding to a global signal, might encode sequential information despite their imprecise timing. It is unclear how neural systems are able to deal with this imprecision. Note that the precision in a real neural network might be improved compared to the one described in this letter because the finite size effect would be reduced. The memory of a sequence may be the result of the coordinated dynamics of multiple encoding populations (as instrument groups in an orchestra). The precision of encoding in distributed representation among multiple local populations is an important issue to be investigated in the future. This encoding scheme is a priori not incompatible with the synfire chain hypothesis, and both schemes might be used concurrently by the brain to encode sequential information.

### 4.4. Influence of Low-Frequency Oscillations.

It has been observed that low-frequency oscillations are correlated with the state of arousal in vivo (Steriade, 1999; Destexhe & Contreras, 2006) owing to the modulating role of basal ganglia (Wichmann, Kliem, & Soares, 2002; Ruskin et al., 1999), dopamine and acethylcholine systems.

The state of arousal can be modeled by reducing the signal-to-noise ratio (Brunel & Wang, 2001). Here, we suggest that low-frequency oscillations constitute an important factor in shaping the attention level. Indeed, we have seen that the temporal patterns of up and down states are modulated by the low-frequency oscillations. In particular, they can synchronize to the oscillatory input for frequencies described in Figures 9A and 10A. At these frequencies, there is a low mutual information between the interburst intervals in Figure 11A. This case would correspond to a phase of restfulness during which regular slow oscillatory patterns have been observed in vivo.

### 4.5. Chaos.

We have shown that the population-averaged firing rate of a local network, when modulated by a global oscillatory input, can exhibit chaotic trajectories in the mean field approximation. It is known that large neural networks can generate chaotic trajectories of activity (Samuelides & Cessac, 2007; Sprott, 2008; Dauce, Quoy, Cessac, Doyon, & Samuelides, 1998; Dauce, Moynot, Pinaud, & Samuelides, 2001; Cessac, Doyon, Quoy, & Samuelides, 1994), especially when excitation and inhibition are balanced (van Vreeswijk & Sompolinsky, 1996; Sompolinsky, Crisanti, & Sommers, 1988; Zillmer, Brunel, & Hansel, 2009; Battaglia, Brunel, & Hansel, 2007). Chaotic trajectories exhibit interesting properties when applied to the encoding of sequences. Moreover, chaotic dynamics can naturally encode sequences into a hierarchy contained in the fractal attractor of chaotic trajectories (Tsuda & Kuroda, 2001). Chaotic patterns of activity in local field potentials (Hayashi & Ishizuka, 1995, 1992) or EEGs have been observed experimentally (Bressler & Freeman, 1980).

### 4.6. Future Directions.

A local population of excitatory and inhibitory neurons connected with two distinct types of dynamic synapses has been analyzed. Previous studies usually consider networks with a single type of dynamic synapses (Tsodyks et al., 1998; Holcman & Tsodyks, 2006; Barak & Tsodyks, 2007; Melamed et al., 2008; Mongillo et al., 2008; Barbieri & Brunel, 2007) or rely on simulations (Verduzco-Flores, Ermentrout, & Bodner, 2012). The difficulty in analyzing networks with multiple dynamic synapses is attributable to the large number of variables to be taken into account. In this letter, we have proposed a method to simplify this analysis by approximating the system to coupled fast-slow subsystems. We have considered here only the case of recurrent dynamic synapses in the excitatory and inhibitory population. However, similar arguments could be used to study other types of networks, particularly those composed of subnetworks with different values of the parameter (i.e., of types E1a, E1b, and so on; see Wang et al., 2006, and Barak & Tsodyks, 2007).

Moreover, it would be interesting to consider the population-averaged firing rates at least at the first order in the mean field approximation to include the effect of fluctuations in the mean field approximation (Brunel & Sergi, 1998). Indeed, the simple mean field approximation based on equations 2.41 and 2.42 does not explain the synchrony observed at the end of up phases when s. An increase in synchrony has also been observed at the onset of representational switching (Sakamoto et al., 2008) and of epileptic seizures (Mormann, Andrzejak, Elger, & Lehnertz, 2007). In the model described here, transitions from up to down phases occur when trajectories reach a fold bifurcation. It has been proposed that early signs of such critical transitions (Scheffer et al., 2009) can include an increase in spike synchrony. The two types of slow oscillations described in this letter (see Figures 8B1 and 8B2) emerge with or without synchrony at the end of the up phase. This illustrates that the increase in synchrony is not a sign of the offset transition in certain cases. We propose that the calculation of synchrony at the end of up phases in vivo provides a method that could be used to discriminate between the two types of slow oscillations we have described.

The simple mean field approximation developed in this letter accounts for the slow oscillations of the network but does not describe the fast oscillations that can occur within the up state. First, fast oscillations can occur due to delayed inhibition (Brunel & Hakim, 1999). Second, short-term plasticity between excitatory neurons can induce a transient response at the transition to the up state (Barak & Tsodyks, 2007), which may take the form of neural avalanches (Millman, Mihalas, Kirkwood, & Niebur, 2010). We are currently investigating the fast oscillations resulting from such transient dynamics (Leleu & Aihara, 2013).

Finally, we have seen that noise has an important influence on the relaxation oscillations. Fluctuations induce random transitions between up and down phases when the trajectory is at the proximity of a fold curve. This kind of stochastic transition is characteristic of systems operating close to a fold bifurcation and can be a common phenomenon in neural systems (Korn & Faure, 2003). It has been proposed that stochastic resonance may embed some deterministic signal in vivo (Faure & Korn, 2001). It would hence be informative to calculate by methods described in Korn and Faure (2003) the amount of determinism of the LIF model described in this letter, especially in the case of a chaotic motion of the population-averaged trajectory.

## Appendix A: Fast-Slow Decomposition

In order to define the critical manifold in section 2.4, we have considered that the variables *u _{II}*(

*t*) and

*x*(

_{EE}*t*) are approximately equal to their steady-state value and . For confirmation of this hypothesis, we verify that the critical manifold is stable (see Figure 12A) and we compare the values of

*u*(

_{II}*t*) and

*x*(

_{EE}*t*) obtained by the simple mean field approximation (see equations 2.41, 2.42, 2.28, and 2.29) and the corresponding steady states and defined in equations 2.46 and 2.47, respectively. Figure 12B shows that and .

## Appendix B: Simulation Parameters

Notation . | Detail . | Value . |
---|---|---|

Timescales of simulation | ||

dt | Discrete time step for simple mean field approximation | s |

dt | Discrete time step for network of LIF neurons | s |

Neuron parameters | ||

Membrane time constant of pyramidal cells | s | |

Membrane time constant of interneurons | s | |

, | Absolute refractory period | s |

V ^{r}_{E} | Reset value of membrane potential for pyramidal cells | V |

V ^{r}_{I} | Reset value of membrane potential for interneurons | V |

, | Spike threshold of pyramidal cells | V |

Network parameters | ||

N _{E} | Number of pyramidal cells | 800 |

N _{I} | Number of interneurons | 200 |

Mean input to pyramidal cells | 0.78 Vs^{−1} | |

Mean input to interneurons | 1.65 Vs^{−1} | |

Standard deviation of input to pyramidal cells | 0.02 Vs^{−0.5} | |

Standard deviation of input to interneurons | 0.029 Vs^{−0.5} |

Notation . | Detail . | Value . |
---|---|---|

Timescales of simulation | ||

dt | Discrete time step for simple mean field approximation | s |

dt | Discrete time step for network of LIF neurons | s |

Neuron parameters | ||

Membrane time constant of pyramidal cells | s | |

Membrane time constant of interneurons | s | |

, | Absolute refractory period | s |

V ^{r}_{E} | Reset value of membrane potential for pyramidal cells | V |

V ^{r}_{I} | Reset value of membrane potential for interneurons | V |

, | Spike threshold of pyramidal cells | V |

Network parameters | ||

N _{E} | Number of pyramidal cells | 800 |

N _{I} | Number of interneurons | 200 |

Mean input to pyramidal cells | 0.78 Vs^{−1} | |

Mean input to interneurons | 1.65 Vs^{−1} | |

Standard deviation of input to pyramidal cells | 0.02 Vs^{−0.5} | |

Standard deviation of input to interneurons | 0.029 Vs^{−0.5} |

Notation . | Detail . | Value . |
---|---|---|

Synaptic transmission parameters | ||

d_{max} | Maximal synaptic delay | s |

J^{0}_{EE} | Synaptic weights from E (pyramidal cells) to E | V |

J^{0}_{EI} | Synaptic weights from I (interneurons) to E | V |

J^{0}_{IE} | Synaptic weights from E to I | V |

J^{0}_{II} | Synaptic weights from I to I | V |

Current parameters | ||

Decay time of NMDA | 0.1 s | |

Rise time of NMDA | 0.002 s | |

Decay time of AMPA | 0.005 s | |

Rise time of AMPA | 0.0005 s | |

Decay time of GABA_{A} | 0.005 s | |

Rise time of GABA_{A} | 0.0005 s | |

Proportion of NMDA currents | 0.25 | |

Dynamic synapse parameters | ||

U ^{d}_{EE} | U Parameter for type E1 dynamic synapses _{d} | 0.35 |

U ^{d}_{II} | U Parameter for type F2 dynamic synapses _{d} | 0.2 |

Type E1 depression time constant | 0.027 s | |

Type E1 facilitation time constant | 1 s | |

Type F2 depression time constant | 1 s | |

Type F2 facilitation time constant | 0.03 s | |

Parameters of the low-frequency oscillatory input | ||

A | Amplitude of low-frequency oscillations | 0.039 |

Phase of low-frequency oscillations | ||

r ^{A}_{E} | Amplitude ratio of pyramidal cells | 1 |

r ^{A}_{I} | Amplitude ratio of interneurons | 2 |

Notation . | Detail . | Value . |
---|---|---|

Synaptic transmission parameters | ||

d_{max} | Maximal synaptic delay | s |

J^{0}_{EE} | Synaptic weights from E (pyramidal cells) to E | V |

J^{0}_{EI} | Synaptic weights from I (interneurons) to E | V |

J^{0}_{IE} | Synaptic weights from E to I | V |

J^{0}_{II} | Synaptic weights from I to I | V |

Current parameters | ||

Decay time of NMDA | 0.1 s | |

Rise time of NMDA | 0.002 s | |

Decay time of AMPA | 0.005 s | |

Rise time of AMPA | 0.0005 s | |

Decay time of GABA_{A} | 0.005 s | |

Rise time of GABA_{A} | 0.0005 s | |

Proportion of NMDA currents | 0.25 | |

Dynamic synapse parameters | ||

U ^{d}_{EE} | U Parameter for type E1 dynamic synapses _{d} | 0.35 |

U ^{d}_{II} | U Parameter for type F2 dynamic synapses _{d} | 0.2 |

Type E1 depression time constant | 0.027 s | |

Type E1 facilitation time constant | 1 s | |

Type F2 depression time constant | 1 s | |

Type F2 facilitation time constant | 0.03 s | |

Parameters of the low-frequency oscillatory input | ||

A | Amplitude of low-frequency oscillations | 0.039 |

Phase of low-frequency oscillations | ||

r ^{A}_{E} | Amplitude ratio of pyramidal cells | 1 |

r ^{A}_{I} | Amplitude ratio of interneurons | 2 |

## Appendix C: Population-Averaged Amount of Available Resources x*_{αβ} and Utilization Factor u*_{αβ}

We extend the method desribed in Romani et al. (2006) to the calculation of and , , under the simplifying assumption that *x _{ij}* and

*u*(given as simply

_{ij}*x*and

*u*in the following for the sake of clarity), , are independent random variables. Moreover, the refractory period is not considered for simplicity.

*u*and are the values of

_{n}*u*(

*t*) immediately before and after a spike

*n*and .

When we assume *u* and *x* as independent and consider the refractory period, then we obtain equations 2.32 after some calculation.

## Appendix D: Stability of the Fast System

The fast subsystem is composed of the following variables: , with , *x _{EE}*, and

*u*. All of these variables should be considered in order to analyze the stability of this subsystem on the critical manifold. However, the population-averaged membrane potential has faster dynamics than the other variables. We consider the stability of the fastest variables , where , assuming that the other variables are not dynamic:

_{II}*u*and

_{EE}*x*are constant and , (see equations 2.46 and 2.47).

_{II}The expression becomes similar to the formula given in equation D.9 when *D* is replaced by *F*.

## Appendix E: Fold Curve on the Critical Manifold

## Appendix F: Measure of Synchrony

*W*and

_{d}*W*, located at the beginning and end of up phases, respectively, defined as where

_{u}*t*and

_{d}*t*are the times of transition from down to up and up to down phases, respectively;

_{u}*w*, the half duration of the time frame defined as ; and

*w*, the offset duration with

_{t}*w*=0.05 s. We use an offset duration to ensure that the fast subsystem has reached its steady state.

_{t}We consider the membrane potential fluctuations of the pyramidal cells only for the calculation of . Using a linear regression with *K*<*N _{E}*, can be determined.

## Acknowledgments

This research is supported by the Global Center of Excellence, “Secure-Life Electronics” sponsored by the Ministry of Education, Culture, Sports, Science and Technology, and by the the Aihara Innovative Mathematical Modelling Project; the Japan Society for the Promotion of Science through the “Funding Program for World-Leading Innovative R&D on Science and Technology (FIRST Program),” initiated by the Council for Science and Technology Policy.