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.

Figure 1:

Network of LIF neurons considered in the simulations. It is composed of 800 pyramidal cells and 200 interneurons connected with dynamic synapses of type E1 (Wang et al., 2006) among pyramidal cells and F2 (Gupta et al., 2000) among interneurons. Other synapses are not dynamic. Oscillatory input representing long-range connections is applied to both pyramidal cells and interneurons in section 3.6.

Figure 1:

Network of LIF neurons considered in the simulations. It is composed of 800 pyramidal cells and 200 interneurons connected with dynamic synapses of type E1 (Wang et al., 2006) among pyramidal cells and F2 (Gupta et al., 2000) among interneurons. Other synapses are not dynamic. Oscillatory input representing long-range connections is applied to both pyramidal cells and interneurons in section 3.6.

2.2.  Neuronal Dynamics.

The membrane potential Vi of a neuron i, where , is described by the following equations (Barbieri & Brunel, 2007):
formula
2.1
formula
2.2
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; INi, NMDA receptor-mediated currents; IAi, AMPA receptor-mediated currents; IGi, GABAA receptor-mediated currents; and Iexti, the external input current. We do not consider metabotropic GABAB receptor-mediated currents in this letter. When the membrane potential Vi reaches the threshold , a spike is emitted and the membrane potential is set to for a duration , which is the absolute refractory period.
In order to model external noise resulting notably from miniature postsynaptic potentials, the external input current Iexti is assumed to be gaussian and -correlated as follows:
formula
2.3
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.
The rise and decay of slow NMDA receptor-mediated synaptic currents are described by the following equations (Barbieri & Brunel, 2007):
formula
2.4
formula
2.5
where and are the decay and rise time, respectively, of NMDA receptor-mediated currents; , the proportion of NMDA receptor-mediated currents; Jij(t), the synaptic weight from neuron to i; tkj, the time of a spike k emitted from neuron j; and d, the synaptic delay. The synaptic delay is distributed uniformely in [0, dmax] where dmax is the maximal synaptic delay. I and U denote decaying and rising currents, respectively.
The rise and decay of rapid AMPA receptor-mediated synaptic currents are described by the following equations (Barbieri & Brunel, 2007):
formula
2.6
formula
2.7
where and are the decay and rise time, respectively, of AMPA receptor-mediated currents.
Finally, the rise and decay of GABAA receptor-mediated synaptic currents are described by
formula
2.8
formula
2.9
where and are the decay and rise time, respectively, of GABAA receptor-mediated currents and Jij(t) is the synaptic weight from neuron to i.
The dynamic synapses Jij(t) of type E1 and F2 are described by the following equation (Tsodyks, Pawelzik, & Markram, 1998; Mongillo et al., 2008):
formula
2.10
where Jij is of type E1 if or F2 if ; xij is the amount of available resources (or ratio of readily releasable vesicles); uij, the utilization factor (proportional to the calcium concentration in the intracellular media); and J0ij, the maximal weight value. For or , xij=uij=1. The dynamics of xij and uij are described by the following equations (Tsodyks et al., 1998; Mongillo et al., 2008):
formula
2.11
formula
2.12
where is the timescale of synaptic depression; , the timescale of synaptic facilitation ( if and if ); and , the parameter of synaptic facilitation. In numerical simulations, uij(t) and xij(t) are calculated just before a spike.
In section 3.6, the mean input oscillates to exhibit the effect of a low-frequency oscillatory input. Low-frequency oscillations may originate from recurrent loops through distal brain areas and intercolumnar interactions or via another local population of neurons. Accordingly, the mean input is described by the following equations:
formula
2.13
formula
2.14
where 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 (rAE=1); , the instantaneous phase; , the initial phase; , the angular velocity; and , the period of low-frequency oscillatory input.

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

We can then derive the following equations to describe the population-averaged synaptic currents (Barbieri & Brunel, 2007) mediated by NMDA receptors:
formula
2.15
formula
2.16
where are the indices of the excitatory or inhibitory populations; , the number of neurons in the population ; , the population-averaged synaptic weight from to ; and , the population-averaged firing rate of . Similarly, the population-averaged synaptic currents mediated by AMPA receptors can be described by the following equations (Barbieri & Brunel, 2007):
formula
2.17
formula
2.18
Likewise, the population-averaged synaptic currents mediated by GABAA receptors can be described by the following equations:
formula
2.19
formula
2.20
We consider, moreover, that synaptic currents, including those mediated by NMDA receptors, vary much faster than membrane potentials such that the currents can be expressed as follows:3
formula
2.21
formula
2.22
formula
2.23
It is usual to assume that in the diffusion approximation (Ricciardi, 1977), the membrane potentials of LIF neurons are described by an Ornstein-Uhlenbeck process. Consequently, the probability density function of the membrane potential satisfies the Fokker-Planck equation. The diffusion approximation can be obtained when postsynaptic potentials (PSP) are considerably smaller than the rates of arrival of presynaptic spikes. The average firing rate of the network can be obtained by calculating the inverse of the mean first passage time using Siegert's recursion formula (Siegert, 1951). Then, by exploiting the condition of self-consistency for the rates in the network (Amit & Tsodyks, 1991),4 we can write the following equations to describe the average steady-state membrane potential of the population (Amit & Tsodyks, 1991; Amit & Brunel, 1997b):
formula
2.24
The average firing rate of population is described by the following equation (Tsodyks et al., 1998; Ricciardi, 1977):
formula
2.25
where ; , the Gauss error function; and , the steady state of when it exists.
For approximating the population-averaged dynamic synaptic weights, we assume that the amount of available resources and the utilization factor of a neuron are independent of each other (Tsodyks et al., 1998). The relative error of this approximation can be estimated using the Cauchy-Schwarz inequality of the probability theory and is bounded from above by the coefficient of variations of the amount of available resources and the utilization factor (Tsodyks et al., 1998). In Tsodyks et al. (1998), the relative error does not exceed 5% when the spike trains have Poisson statistics. We thus consider the following assumption:
formula
2.26
Consequently, the population-averaged weights are given as
formula
2.27
where .
The approximation of the population-averaged dynamic synaptic weights can be achieved by two methods. The first method assumes that for the calculation of , the arriving spikes have Poisson statistics (Tsodyks et al., 1998). In this case, the equations describing the population-averaged synaptic dynamics are given as follows ( or ):
formula
2.28
formula
2.29
From equations 2.28 and 2.29, the steady states of and , given as and , respectively, are derived as follows:
formula
2.30
formula
2.31
The second, more rigorous method is to consider the diffusion approximation of the membrane potential for the calculation of (Romani et al., 2006), similar to the calculation of population-averaged membrane potential. The stable population-averaged values of and , , and , respectively, can then be written as a function of the Laplace transform of the interspike interval distribution function, which in turn can be obtained from the Fokker-Plank equation (Romani et al., 2006). We extend the method described in Romani et al. (2006) to calculate the steady-state population-averaged value of the amount of available resources and the utilization factor . In this case, and are given as follows (see appendix  C):
formula
2.32
formula
2.33
where and . Further, we note the distribution of interspike intervals p(t). Then, and , where is the Laplace transform of p(t), are given as follows (Alili, Patie, & Pedersen, 2005; Romani et al., 2006):
formula
2.34
formula
2.35
where Hl(z) are the Hermite functions. For a numerical simulation, we use hypergeometric functions to evaluate the Hermite functions as follows (Lebedev & Silverman, 1972):
formula
2.36
where is the gamma function and 1F1 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 xEE and uII. These four variables represent the fast subsystem. We consider the case under which there is bistability in the fast system. The two other variables, uEE and xII, constitute the slow subsystem.

Fast-slow systems can be studied using classical methods of the singular perturbation theory (Verhulst, 2005; Guckenheimer & Holmes, 1983; Hoppensteadt & Izhikevich, 1997; Guckenheimer, Wechselberger, & Young, 2006). Consider the fast-slow system described as follows:
formula
2.37
formula
2.38
where is a very small constant. In the limit , the fast system converges to its critical manifold, defined by , while the slow dynamics is constrained on this manifold. In domain 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:
formula
2.39
The dynamics of the slow system can then be studied independently of the fast one in each domain I by considering the following system:
formula
2.40

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

We consider to a first approximation that the instantaneous population-averaged membrane potentials and , respectively, are described by the following equations (Romani et al., 2006; Amit & Brunel, 1997b):
formula
2.41
formula
2.42
where 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.
The critical manifold is defined as follows (see equations 2.24 to 2.27):
formula
2.43
formula
2.44
formula
2.45
The slow subsystem is then described exclusively by the dynamics of uEE(t) and xII(t) on the critical manifold.
To study the critical manifold, we need to solve equation 2.43 for in the slow variables uEE and xII. 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 uEE as follows:8
formula
2.46
formula
2.47
The expressions of fx and fu 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:
formula
2.48
formula
2.49
where and are defined as follows (see equations 2.27, 2.43 and 2.46, and 2.47):
formula
2.50
formula
2.51
By solving equations 2.48 and 2.49 for variables and , we can decouple the slow and fast variables, as shown in equations 2.39 and 2.40. In each domain I where the implicit function theorem is applicable,
formula
2.52
formula
2.53

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 uEE and xII, that is, [0, 1]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 xEE and uII 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).

2.5.  Fold Curves.

To investigate the jumps between sheets of the critical manifold, we compute the bifurcations of the fast system on the critical manifold when the slow variables uEE and xII 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):
formula
2.54
formula
2.55

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.

We assume here an exponential distribution for the interspike interval distribution, although the position of the critical manifold can be obtained by assuming the Ornstein-Uhlenbeck process as shown in equations 2.32 and 2.33. We can express the slow dynamics of the system on the critical manifold as follows by using equations 2.52 and 2.53 in each domain I for which equations 2.48 and 2.49 can be solved:
formula
2.56
formula
2.57
For each domain 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
formula
2.58
formula
2.59
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 uEE and xII 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 uEE and xII ([0, 1]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 FC1 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 FC2 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.

Figure 2:

(A1, A2) Critical manifold representing the steady states of the fast system plotted in the space and , respectively. The dots are the solutions obtained by solving equations 2.48, 2.49, and 2.54 using a numerical solver. The dots are stable fixed points of the fast subsystem on the stable manifold and are unstable fixed points on the unstable manifold. The gray surfaces are the result of interpolation between the dots. (A3) Fold curves solution of equation 2.54 (dots) and interpolation (curves). The left-most fold curve is given as FC1 and the right-most one as FC2. Domain I1 is located to the left of the fold curve FC2; I3, to the right of the fold curve FC1; and I2, between FC1 and FC2. A trajectory obtained by using the simple mean field approximation and parameters given in appendix  B is superimposed (black curve) in panels A1, A2, and A3. (B1, B2, B3) The same as panels A1, A2, and A3, but equations 2.32 and 2.33 are used to calculate the critical manifold, and the superimposed trajectories are the population-averaged dynamics simulated using the network of LIF neurons smoothened with a moving average window of 0.025 s.

Figure 2:

(A1, A2) Critical manifold representing the steady states of the fast system plotted in the space and , respectively. The dots are the solutions obtained by solving equations 2.48, 2.49, and 2.54 using a numerical solver. The dots are stable fixed points of the fast subsystem on the stable manifold and are unstable fixed points on the unstable manifold. The gray surfaces are the result of interpolation between the dots. (A3) Fold curves solution of equation 2.54 (dots) and interpolation (curves). The left-most fold curve is given as FC1 and the right-most one as FC2. Domain I1 is located to the left of the fold curve FC2; I3, to the right of the fold curve FC1; and I2, between FC1 and FC2. A trajectory obtained by using the simple mean field approximation and parameters given in appendix  B is superimposed (black curve) in panels A1, A2, and A3. (B1, B2, B3) The same as panels A1, A2, and A3, but equations 2.32 and 2.33 are used to calculate the critical manifold, and the superimposed trajectories are the population-averaged dynamics simulated using the network of LIF neurons smoothened with a moving average window of 0.025 s.

Between the fold curves FC1 and FC2, 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).

Figure 3:

(A1, B1, C1) Nullclines and steady states of the slow system for domains I1 and I3 corresponding to the lower and upper sheets of the critical manifold, respectively. The nullclines are denoted by the set of dark gray points, which are the result of calculation using the simple mean field approximation, connected by dark gray lines, which are interpolated from the data. The steady states, which are the intersection of the nullclines, are denoted by black circles. The fold curves FC1 and FC2 delimiting the domains I3 and I1, respectively, are denoted by light gray curves. (A2, B2, C2) Fold curves (black dots resulting from calculation and gray curves interpolated) and a sample trajectory simulated using the simple mean field approximation (black curve directed by an arrow). (A3, B3, C3) Time series of the population-averaged firing rates of pyramidal cells and interneurons, and , respectively, and the slow variables uEE and xII obtained using the simple mean field approximation. The gray square wave denotes up and down phases; panels A1, A2, and A3 are simulated using the parameters given in appendix  B. (B1, B2, B3) V. (C1, C2, C3) V and UdEE=0.46.

Figure 3:

(A1, B1, C1) Nullclines and steady states of the slow system for domains I1 and I3 corresponding to the lower and upper sheets of the critical manifold, respectively. The nullclines are denoted by the set of dark gray points, which are the result of calculation using the simple mean field approximation, connected by dark gray lines, which are interpolated from the data. The steady states, which are the intersection of the nullclines, are denoted by black circles. The fold curves FC1 and FC2 delimiting the domains I3 and I1, respectively, are denoted by light gray curves. (A2, B2, C2) Fold curves (black dots resulting from calculation and gray curves interpolated) and a sample trajectory simulated using the simple mean field approximation (black curve directed by an arrow). (A3, B3, C3) Time series of the population-averaged firing rates of pyramidal cells and interneurons, and , respectively, and the slow variables uEE and xII obtained using the simple mean field approximation. The gray square wave denotes up and down phases; panels A1, A2, and A3 are simulated using the parameters given in appendix  B. (B1, B2, B3) V. (C1, C2, C3) V and UdEE=0.46.

For each domain delimited by the fold curves given as I1, I2 and I3 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. I1, situated to the left of curve FC2, is the domain of the lower sheet. Similarly, I3 situated to the right of curve FC1 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 I1 and I3, 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 J0EE. Similar to the first case, there is no equilibrium for the slow system, and the trajectory is periodic. Stronger recurrent connections J0EE increase the duration of the up phase (+0.04 s) and reduce the duration of the down phase (−0.56 s).

Finally, the third case corresponds to an increase in JEE and UdEE.10 In this case, the slow system has a unique equilibrium on the upper manifold (see Figure 3C1). Consequently, the trajectory of the slow system remains in the up phase (see Figures 3C2 and 3C3).

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.

Figure 4:

(A1, B1, C1) Nullclines and steady states of the slow system for domains I1 and I3 corresponding to the lower and upper sheets of the critical manifold, respectively, obtained using equations 2.32 and 2.33. (A2, B2, C2) Time series of the population-averaged firing rates of pyramidal cells and interneurons, and , respectively, and the slow variables uEE and xII obtained by simulating the network of LIF neurons. The up and down phases denoted by the gray square wave in panels A2, B2, and C2 are detected based on the bimodal distribution of the logarithm of the firing rate shown in panels A3, B3, and C3, respectively.

Figure 4:

(A1, B1, C1) Nullclines and steady states of the slow system for domains I1 and I3 corresponding to the lower and upper sheets of the critical manifold, respectively, obtained using equations 2.32 and 2.33. (A2, B2, C2) Time series of the population-averaged firing rates of pyramidal cells and interneurons, and , respectively, and the slow variables uEE and xII obtained by simulating the network of LIF neurons. The up and down phases denoted by the gray square wave in panels A2, B2, and C2 are detected based on the bimodal distribution of the logarithm of the firing rate shown in panels A3, B3, and C3, respectively.

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 FC1 and FC2, 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 Ns(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 n1(t) of the firing rate , (Brunel & Sergi, 1998).

Figure 5:

(A1, A2) Effect of the number of neurons N in the network of LIF neurons (normalized using the number of neurons given in appendix  B) on the mean duration of down and up phases, respectively. (B1, B2) Standard deviation of the instantaneous firing rates and , respectively, during the down phase. (C1, C2) The same as panels B1 and B2, but during the up phase.

Figure 5:

(A1, A2) Effect of the number of neurons N in the network of LIF neurons (normalized using the number of neurons given in appendix  B) on the mean duration of down and up phases, respectively. (B1, B2) Standard deviation of the instantaneous firing rates and , respectively, during the down phase. (C1, C2) The same as panels B1 and B2, but during the up phase.

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.

Figure 6:

(A, B) Effect of the proportion of NMDA currents on the mean duration of down and up phases, respectively.

Figure 6:

(A, B) Effect of the proportion of NMDA currents on the mean duration of down and up phases, respectively.

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 xEE (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.

Figure 7:

(A1, B1) Time series of the population-averaged firing rate of pyramidal cells and interneurons, and , respectively. (A2, B2) Time series of the population-averaged amount of available resource xEE, utilization factor uEE, and normalized instantaneous synaptic strength JEE=xEEuEE between pyramidal cells. (A3, B3) The same as panels A2 and B2, but for dynamic synapses between interneurons. (A1, A2, A3) Simple mean field approximation. (B1, B2, B3) Network of LIF neurons.

Figure 7:

(A1, B1) Time series of the population-averaged firing rate of pyramidal cells and interneurons, and , respectively. (A2, B2) Time series of the population-averaged amount of available resource xEE, utilization factor uEE, and normalized instantaneous synaptic strength JEE=xEEuEE between pyramidal cells. (A3, B3) The same as panels A2 and B2, but for dynamic synapses between interneurons. (A1, A2, A3) Simple mean field approximation. (B1, B2, B3) Network of LIF neurons.

Initially, the system is in the down phase between s and s (see Figure 7). The strength of recurrent excitatory connections JEE increases slowly because of synaptic facilitation, whereas the strength of recurrent inhibitory connections JII 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.

During the up phase, the decrease in firing rate of pyramidal cells is concurrent with a decrease in the strength of recurrent inhibitory connections JII. 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 xEE recovers gradually.

Further decreases of JII and of JEE, resulting from a slow increase of uEE and decrease of xII, 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.

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 uEE is considered as a fast variable and xII 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).11 In the previous case ( s), the membrane potential fluctuations at the end of up phases were more asynchronous (see Figure 8B1).

Figure 8:

Two types of slow oscillations. (A) Frequency fn of network oscillation in hertz represented by shades of gray in the parameter space {, JpEE} with V. In the lower left and upper right corners, the network settles to the down and up state, respectively. (B1, B2) Raster plot superimposed on the measure of synchrony of the membrane potentials (see appendix  F). Pyramidal cells and interneurons are evenly spaced along a line of length 8 mm and indexed by their distance di in millimeters to the center. Excitatory and inhibitory spikes are represented by dots and crosses, respectively. The right and left boxes are the values of synchrony measure calculated at the beginning and end of the up phase, respectively. The black lines crossing these boxes are the 95% confidence intervals. (B1) s, . (B2) s, . (B3) Firing rates of excitatory and inhibitory neurons simulated using the simple mean field approximation are denoted by dashed black and gray curves, respectively. The solid curves denote the same as the dashed curves with , , and . The parameters are the same as in panel B2. (C1, C2) Averaged firing rates of excitatory (black) and inhibitory (gray) neurons at the offset of the up phase rescaled to fit in [0, 1] and smoothened with a moving average window of 0.006 s. In panels C1 and C2, the parameters are the same as in panels B1 and B2, respectively.

Figure 8:

Two types of slow oscillations. (A) Frequency fn of network oscillation in hertz represented by shades of gray in the parameter space {, JpEE} with V. In the lower left and upper right corners, the network settles to the down and up state, respectively. (B1, B2) Raster plot superimposed on the measure of synchrony of the membrane potentials (see appendix  F). Pyramidal cells and interneurons are evenly spaced along a line of length 8 mm and indexed by their distance di in millimeters to the center. Excitatory and inhibitory spikes are represented by dots and crosses, respectively. The right and left boxes are the values of synchrony measure calculated at the beginning and end of the up phase, respectively. The black lines crossing these boxes are the 95% confidence intervals. (B1) s, . (B2) s, . (B3) Firing rates of excitatory and inhibitory neurons simulated using the simple mean field approximation are denoted by dashed black and gray curves, respectively. The solid curves denote the same as the dashed curves with , , and . The parameters are the same as in panel B2. (C1, C2) Averaged firing rates of excitatory (black) and inhibitory (gray) neurons at the offset of the up phase rescaled to fit in [0, 1] and smoothened with a moving average window of 0.006 s. In panels C1 and C2, the parameters are the same as in panels B1 and B2, respectively.

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 xII 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 xII(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 fin, the resulting dynamics of the simple mean field approximation exhibits some complex deterministic trajectories (see Figure 9).

Figure 9:

Effect of slow oscillatory input on the oscillations of the network using the simple mean field approximation. (A) Bifurcation diagram of , where is the duration between the onsets of two consecutive up phases with the frequency of oscillatory input Hz as a bifurcation parameter. We observe that at Hz, there is 1:1 synchronization between slow oscillations produced by the local network and the oscillatory input; at Hz, there is 1:2 subharmonic synchronization; at Hz, there is 1:3 subharmonic synchronization; and at Hz, there is 1:4 subharmonic synchronization. (B) The same as panel A, but with Hz. (C1, C2, C3) Time series of the population-averaged firing rate and . The gray square wave denotes up and down phases. The black dots above the square wave denote the time when the oscillatory input peaks. (C1) fin=0.6125 Hz. (C2) fin=0.992 Hz. (C3) fin=1.1 Hz. (D1, D2, D3) Critical manifold and superimposed trajectory simulated using the simple mean field approximation in the space corresponding to panels C1, C2, and C3, respectively.

Figure 9:

Effect of slow oscillatory input on the oscillations of the network using the simple mean field approximation. (A) Bifurcation diagram of , where is the duration between the onsets of two consecutive up phases with the frequency of oscillatory input Hz as a bifurcation parameter. We observe that at Hz, there is 1:1 synchronization between slow oscillations produced by the local network and the oscillatory input; at Hz, there is 1:2 subharmonic synchronization; at Hz, there is 1:3 subharmonic synchronization; and at Hz, there is 1:4 subharmonic synchronization. (B) The same as panel A, but with Hz. (C1, C2, C3) Time series of the population-averaged firing rate and . The gray square wave denotes up and down phases. The black dots above the square wave denote the time when the oscillatory input peaks. (C1) fin=0.6125 Hz. (C2) fin=0.992 Hz. (C3) fin=1.1 Hz. (D1, D2, D3) Critical manifold and superimposed trajectory simulated using the simple mean field approximation in the space corresponding to panels C1, C2, and C3, respectively.

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 fin. 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 fin as a bifurcation parameter. Different modes of forced oscillations are shown in Figures 9C1 to 9C3 and 9D1 to 9D3.

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

Figure 10:

The same as Figure 9, but by simulating the network of LIF neurons and using equations 2.32 and 2.33 for determining the critical manifold in (D1,D2,D3). Hz. NE=4800, NI=1200. (A) Note that synchronization modes similar to Figure 9A appear in the case of the simulation of LIF neurons. Nonetheless, fluctuate at the proximity of the harmonics and subharmonics. (B) The same as panel A, but for the inverse of up phase duration . (C1, D1) fin=0.32 Hz. (C2, D2) fin=0.53 Hz. (C3, D3) fin=1 Hz.

Figure 10:

The same as Figure 9, but by simulating the network of LIF neurons and using equations 2.32 and 2.33 for determining the critical manifold in (D1,D2,D3). Hz. NE=4800, NI=1200. (A) Note that synchronization modes similar to Figure 9A appear in the case of the simulation of LIF neurons. Nonetheless, fluctuate at the proximity of the harmonics and subharmonics. (B) The same as panel A, but for the inverse of up phase duration . (C1, D1) fin=0.32 Hz. (C2, D2) fin=0.53 Hz. (C3, D3) fin=1 Hz.

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.

Figure 11:

(A) Mutual information between the durations of interburst intervals and their first and second successors, and , respectively, as a function of the oscillatory input frequency fin. (B) Mutual information between the durations of up phases and their first and second successors, and , respectively, as a function of the oscillatory input frequency fin. (C1, C2, C3) Phase space of the durations of interburst intervals and successors . (D1, D2, D3) Phase space of up phase durations and successors . (C1, D1) fin=0.32 Hz. (C2, D2) fin=0.53 Hz. (C3, D3) fin=1 Hz. Hz. NE=4800, NI=1200.

Figure 11:

(A) Mutual information between the durations of interburst intervals and their first and second successors, and , respectively, as a function of the oscillatory input frequency fin. (B) Mutual information between the durations of up phases and their first and second successors, and , respectively, as a function of the oscillatory input frequency fin. (C1, C2, C3) Phase space of the durations of interburst intervals and successors . (D1, D2, D3) Phase space of up phase durations and successors . (C1, D1) fin=0.32 Hz. (C2, D2) fin=0.53 Hz. (C3, D3) fin=1 Hz. Hz. NE=4800, NI=1200.

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 uEE and xII. We have shown that the dynamics of these slow variables is sufficient to explain the occurrence of the slow oscillations.

We have identified two types of slow oscillations according to the timescale of short-term facilitation between pyramidal cells. When both uEE and xII 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 xII 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 GABAB 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 uII(t) and xEE(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 uII(t) and xEE(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 .

Figure 12:

Verification of the fast-slow decomposition. (A1, A2, A3, A4) Stability analysis of the critical manifold based on the linearization of the four equations describing the fast system. Dots and crosses denote stable and unstable fixed points of the fast subsystem representing the critical manifold. The values of , , , and are shown in panels A1, A2, A3, and A4, respectively. (B) The time series of the instantaneous population-averaged amount of available resources xEE(t) obtained by using the simple mean field approximation and the corresponding steady state (see equation 2.46). We see that xEE(t) and are superimposed. Similarly, uII(t) and are superimposed.

Figure 12:

Verification of the fast-slow decomposition. (A1, A2, A3, A4) Stability analysis of the critical manifold based on the linearization of the four equations describing the fast system. Dots and crosses denote stable and unstable fixed points of the fast subsystem representing the critical manifold. The values of , , , and are shown in panels A1, A2, A3, and A4, respectively. (B) The time series of the instantaneous population-averaged amount of available resources xEE(t) obtained by using the simple mean field approximation and the corresponding steady state (see equation 2.46). We see that xEE(t) and are superimposed. Similarly, uII(t) and are superimposed.

Appendix B: Simulation Parameters

Table 1:
Simulations Parameters (1/2).
NotationDetailValue
Timescales of simulation   
   dt Discrete time step for simple mean field approximation 
   dt Discrete time step for network of LIF neurons 
Neuron parameters   
    Membrane time constant of pyramidal cells 
    Membrane time constant of interneurons 
   ,  Absolute refractory period 
   VrE Reset value of membrane potential for pyramidal cells 
   VrI Reset value of membrane potential for interneurons 
   ,  Spike threshold of pyramidal cells 
Network parameters   
   NE Number of pyramidal cells 800 
   NI 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 
NotationDetailValue
Timescales of simulation   
   dt Discrete time step for simple mean field approximation 
   dt Discrete time step for network of LIF neurons 
Neuron parameters   
    Membrane time constant of pyramidal cells 
    Membrane time constant of interneurons 
   ,  Absolute refractory period 
   VrE Reset value of membrane potential for pyramidal cells 
   VrI Reset value of membrane potential for interneurons 
   ,  Spike threshold of pyramidal cells 
Network parameters   
   NE Number of pyramidal cells 800 
   NI 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 
Table 2:
Simulations Parameters (2/2).
NotationDetailValue
Synaptic transmission parameters 
   dmax Maximal synaptic delay 
   J0EE Synaptic weights from E (pyramidal cells) to E 
   J0EI Synaptic weights from I (interneurons) to E 
   J0IE Synaptic weights from E to I 
   J0II Synaptic weights from I to I 
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 GABAA 0.005 s 
    Rise time of GABAA 0.0005 s 
    Proportion of NMDA currents 0.25 
Dynamic synapse parameters   
   UdEE Ud Parameter for type E1 dynamic synapses 0.35 
   UdII Ud Parameter for type F2 dynamic synapses 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  
   rAE Amplitude ratio of pyramidal cells 
   rAI Amplitude ratio of interneurons 
NotationDetailValue
Synaptic transmission parameters 
   dmax Maximal synaptic delay 
   J0EE Synaptic weights from E (pyramidal cells) to E 
   J0EI Synaptic weights from I (interneurons) to E 
   J0IE Synaptic weights from E to I 
   J0II Synaptic weights from I to I 
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 GABAA 0.005 s 
    Rise time of GABAA 0.0005 s 
    Proportion of NMDA currents 0.25 
Dynamic synapse parameters   
   UdEE Ud Parameter for type E1 dynamic synapses 0.35 
   UdII Ud Parameter for type F2 dynamic synapses 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  
   rAE Amplitude ratio of pyramidal cells 
   rAI Amplitude ratio of interneurons 

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 xij and uij (given as simply x and u in the following for the sake of clarity), , are independent random variables. Moreover, the refractory period is not considered for simplicity.

From equations 2.12, we can write the following recurrence equation:
formula
C.1
formula
C.2
where un and are the values of u(t) immediately before and after a spike n and .
From equations C.1 and C.2, we can obtain the following expression:
formula
C.3
Then, by averaging a population of neurons in equations C.3, and taking , we obtain the following equation:12
formula
C.4
formula
C.5
formula
C.6
where p(t) is the distribution of interspike intervals and the Laplace transform of p(t). After some algebra, and by considering the refractory period, we obtain equations 2.33.
A recurrence relation similar to equation C.2 can be obtained for x (Romani et al., 2006):
formula
C.7
where xn is the value of x(t) immediately before the spike n. By averaging the population, we have the following relation ():
formula
C.8

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 , xEE, and uII. 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: uEE and xII are constant and , (see equations 2.46 and 2.47).

The time evolution of the variable can be written as follows (Tuckwell, 1988; Amit & Brunel, 1997b):
formula
D.1
formula
D.2
The Jacobian matrix corresponding to the linearization of the functions and of equations D.1 and D.2 is given as follows:
formula
D.3
The derivatives of the population-averaged firing rates according to the population-averaged membrane potential, , are given as
formula
D.4
From the approximations given in equations 2.30 and 2.31, the derivatives of the population-averaged amount of the available resource and the utilization factor according to the membrane potential and , respectively, are given as
formula
D.5
formula
D.6
From the approximations given in equations 2.32 and 2.33, the derivatives of the population-averaged amount of the available resource and the utilization factor according to the membrane potential and , respectively, are given as
formula
D.7
formula
D.8
where
formula
D.9
and
formula
D.10

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

Points on the critical manifold are given by solving equations 2.48 and 2.49. Using equations 2.50 and 2.51 and the approximation of equations 2.30 and 2.31 for fx and fu, respectively, we can write the slow variable uEE as a function of and as
formula
E.1
When equations 2.32 and 2.33 are used for fx and fu, respectively, the slow variable uEE can be expressed as
formula
E.2
Lastly, the slow variable xII can be written as function of and as (using equations 2.31 or 2.33 for ):
formula
E.3
Therefore, the Jacobian matrix in equations D.3 of appendix  D can be written as a function of and as
formula
E.4

Using the Jacobian matrix , the eigenvalue and of equations 2.54 and 2.55 can be obtained.

Appendix F:  Measure of Synchrony

In Figures 8B1 and 8B2, we have calculated the measure of synchrony of membrane potential fluctuations. Because we consider that the fast subsystem reaches a steady state during the up phase of the slow oscillation, we can apply the measure of synchrony of membrane potential fluctuations described in Hansel and Sompolinsky (1992) and Ginzburg and Sompolinsky (1994) to the fast subsystem. The synchrony is calculated using two time frames Wd and Wu, located at the beginning and end of up phases, respectively, defined as
formula
F.1
formula
F.2
where td and tu are the times of transition from down to up and up to down phases, respectively; w, the half duration of the time frame defined as ; and wt, the offset duration with wt=0.05 s. We use an offset duration to ensure that the fast subsystem has reached its steady state.
Then, is calculated as follows (Hansel & Sompolinsky, 1992; Ginzburg & Sompolinsky, 1994) within each time frame:
formula
F.3
where is the ensemble average, that is, and is the time average over W, or . Then can be expressed as follows (Hansel & Sompolinsky, 1992; Ginzburg & Sompolinsky, 1994):
formula
F.4

We consider the membrane potential fluctuations of the pyramidal cells only for the calculation of . Using a linear regression with K<NE, 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.

References

Abeles
,
M.
,
Bergman
,
H.
,
Margalit
,
E.
, &
Vaadia
,
E.
(
1993
).
Spatiotemporal firing patterns in the frontal-cortex of behaving monkeys
.
Journal of Neurophysiology
,
70
(
4
),
1629
1638
.
Alili
,
A.
,
Patie
,
P.
, &
Pedersen
,
J. L.
(
2005
).
Representations of the first hitting time density of an Ornstein-Uhlenbeck process
.
Stochastic Models
,
21
(
4
),
967
980
.
Amari
,
S. I.
(
1971
).
Characteristics of randomly connected threshold-element networks and network systems
.
Proceedings of the Institute of Electrical and Electronics Engineers
,
59
(
1
),
35
47
.
Amit
,
D. J.
, &
Brunel
,
N.
(
1997a
).
Dynamics of a recurrent network of spiking neurons before and following learning
.
Network-Computation in Neural Systems
,
8
(
4
),
373
404
.
Amit
,
D. J.
, &
Brunel
,
N.
(
1997b
).
Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex
.
Cerebral Cortex
,
7
(
3
),
237
252
.
Amit
,
D. J.
, &
Tsodyks
,
M. V.
(
1991
).
Quantitative study of attractor neural network retrieving at low spike rates.1. Substrate spikes, rates and neuronal gain
.
Network Computation in Neural Systems
,
2
(
3
),
259
273
.
Arieli
,
A.
,
Sterkin
,
A.
,
Grinvald
,
A.
, &
Aertsen
,
A.
(
1996
).
Dynamics of ongoing activity: Explanation of the large variability in evoked cortical responses
.
Science
,
273
(
5283
),
1868
1871
.
Baker
,
S. N.
, &
Lemon
,
R. N.
(
2000
).
Precise spatiotemporal repeating patterns in monkey primary and supplementary motor areas occur at chance levels
.
Journal of Neurophysiology
,
84
(
4
),
1770
1780
.
Barak
,
O.
, &
Tsodyks
,
M.
(
2007
).
Persistent activity in neural networks with dynamic synapses
.
Plos Computational Biology
,
3
(
2
),
e35
.
Barbieri
,
F.
, &
Brunel
,
N.
(
2007
).
Irregular persistent activity induced by synaptic excitatory feedback
.
Frontiers in Computational Neuroscience
,
1
,
5
.
Battaglia
,
D.
,
Brunel
,
N.
, &
Hansel
,
D.
(
2007
).
Temporal decorrelation of collective oscillations in neural networks with local inhibition and long-range excitation
.
Physical Review Letters
,
99
(
23
),
238106
.
Battaglia
,
F. P.
,
Sutherland
,
G. R.
, &
McNaughton
,
B. L.
(
2004
).
Hippocampal sharp wave bursts coincide with neocortical “up-state” transitions
.
Learning and Memory
,
11
(
6
),
697
704
.
Bazhenov
,
M.
,
Timofeev
,
I.
,
Steriade
,
M.
, &
Sejnowski
,
T. J.
(
2002
).
Model of thalamocortical slow-wave sleep oscillations and transitions to activated states
.
Journal of Neuroscience
,
22
(
19
),
8691
8704
.
Born
,
J.
,
Rasch
,
B.
, &
Gais
,
S.
(
2006
).
Sleep to remember
.
Neuroscientist
,
12
(
5
),
410
424
.
Bressler
,
S. L.
, &
Freeman
,
W. J.
(
1980
).
Frequency-analysis of olfactory system EEG in cat, rabbit, and rat
.
Electroencephalography and Clinical Neurophysiology
,
50
(
1–2
),
19
24
.
Brette
,
R.
,
Rudolph
,
M.
,
Carnevale
,
T.
,
Hines
,
M.
,
Beeman
,
D.
,
Bower
,
J. M.
, … &
Destexhe
,
A.
(
2007
).
Simulation of networks of spiking neurons: A review of tools and strategies
.
Journal of Computational Neuroscience
,
23
(
3
),
349
398
.
Brunel
,
N.
, &
Hakim
,
V.
(
1999
).
Fast global oscillations in networks of integrate-and-fire neurons with low firing rates
.
Neural Computation
,
11
(
7
),
1621
1671
.
Brunel
,
N.
, &
Sergi
,
S.
(
1998
).
Firing frequency of leaky integrate-and-fire neurons with synaptic current dynamics
.
Journal of Theoretical Biology
,
195
(
1
),
87
95
.
Brunel
,
N.
, &
Wang
,
X. J.
(
2001
).
Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition
.
Journal of Computational Neuroscience
,
11
(
1
),
63
85
.
Brunel
,
N.
, &
Wang
,
X. J.
(
2003
).
What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance
.
Journal of Neurophysiology
,
90
(
1
),
415
430
.
Buzas
,
P.
,
Eysel
,
U. T.
,
Adorjan
,
P.
, &
Kisvarday
,
Z. F.
(
2001
).
Axonal topography of cortical basket cells in relation to orientation, direction, and ocular dominance maps
.
Journal of Comparative Neurology
,
437
(
3
),
259
285
.
Buzsaki
,
G.
, &
Wang
,
X. J.
(
2012
).
Mechanisms of gamma oscillations
.
Annual Review of Neuroscience
,
35
,
203
225
.
Cessac
,
B.
,
Doyon
,
B.
,
Quoy
,
M.
, &
Samuelides
,
M.
(
1994
).
Mean-field equations, bifurcation map and route to chaos in discrete-time neural networks
.
Physica D
,
74
(
1–2
),
24
44
.
Chauvette
,
S.
,
Volgushev
,
M.
, &
Timofeev
,
I.
(
2010
).
Origin of active states in local neocortical networks during slow sleep oscillation
.
Cerebral Cortex
,
20
(
11
),
2660
2674
.
Chen
,
J. Y.
,
Chauvette
,
S.
,
Skorheim
,
S.
,
Timofeev
,
I.
, &
Bazhenov
,
M.
(
2012
).
Interneuron-mediated inhibition synchronizes neuronal activity during slow oscillation
.
Journal of Physiology–London
,
590
(
3
),
3987
4010
.
Cleveland
,
W. S.
(
1979
).
Robust locally weighted regression and smoothing scatterplots
.
Journal of the American Statistical Association
,
74
(
368
),
829
836
.
Compte
,
A.
,
Constantinidis
,
C.
,
Tegner
,
J.
,
Raghavachari
,
S.
,
Chafee
,
M. V.
,
Goldman-Rakic
,
P. S.
, &
Wang
,
W. J.
(
2003
).
Temporally irregular mnemonic persistent activity in prefrontal neurons of monkeys during a delayed response task
.
Journal of Neurophysiology
,
90
(
5
),
3441
3454
.
Compte
,
A.
,
Reig
,
R.
, &
Sanchez-Vives
,
M. V.
(
2009
).