Login using

You can login by using one of your existing accounts.

We will be provided with an authorization token (please note: passwords are not shared with us) and will sync your accounts for you. This means that you will not need to remember your user name and password in the future and you will be able to login with the account you choose to sync, with the click of a button.

Dynamics of networks of excitatory and inhibitory neurons in response to time-dependent inputs

We investigate the dynamics of recurrent networks of excitatory (E) and inhibitory (I) neurons in the presence of time-dependent inputs. The dynamics is characterized by the network dynamical transfer function, i.e., how the population firing rate is modulated by sinusoidal inputs at arbitrary frequencies. Two types of networks are studied and compared: (i) a Wilson–Cowan type firing rate model; and (ii) a fully connected network of leaky integrate-and-fire (LIF) neurons, in a strong noise regime. We first characterize the region of stability of the “asynchronous state” (a state in which population activity is constant in time when external inputs are constant) in the space of parameters characterizing the connectivity of the network. We then systematically characterize the qualitative behaviors of the dynamical transfer function, as a function of the connectivity. We find that the transfer function can be either low-pass, or with a single or double resonance, depending on the connection strengths and synaptic time constants. Resonances appear when the system is close to Hopf bifurcations, that can be induced by two separate mechanisms: the I–I connectivity and the E–I connectivity. Double resonances can appear when excitatory delays are larger than inhibitory delays, due to the fact that two distinct instabilities exist with a finite gap between the corresponding frequencies. In networks of LIF neurons, changes in external inputs and external noise are shown to be able to change qualitatively the network transfer function. Firing rate models are shown to exhibit the same diversity of transfer functions as the LIF network, provided delays are present. They can also exhibit input-dependent changes of the transfer function, provided a suitable static non-linearity is incorporated.

1 Introduction

Networks of neurons in the central nervous system are driven by stimuli that vary on a wide range of time scales, and need to encode these stimuli by the pattern of firing of their constituent neurons. To understand how this encoding is performed, one needs to understand the relationship between the input to the network (the set of spike trains of all neurons that are pre-synaptic to a given network) and its output (the set of spike trains of all neurons in the network). This is a question that is extremely hard to address experimentally, since one has typically access only to a very small fraction of neurons in a network (with the retina being a notable exception Pillow et al., 2008), or access to variables that only reflect population activity, such as the multi-unit activity (MUA) or local field potential (LFP). It is also very hard to address in general in theoretical models, since inputs and outputs live in a high-dimensional space.

A simpler question, and a useful first step, is to understand the relationship between the mean inputs to a network and its population firing rate (i.e., average instantaneous probability that a neuron in the network emits an action potential). However, even this question can be difficult to answer in the space of all possible time-dependent inputs. A natural strategy to study this question is to start with time-independent inputs. One then computes how the mean population rate depends on the mean input (in other words, computing a network “f–I curve”). The next step is to add to the time-independent input a sinusoidal input with weak amplitude, so that the population firing rate will depend linearly on the input, and therefore oscillate with the same frequency as the input. One can then characterize the dynamics of the network with only two functions: the amplitude and the phase of the population firing rate modulation as a function of frequency. These two quantities together characterize the transfer function of the network, and can be used to reconstruct the dynamics of the population in response to arbitrary time-dependent inputs, provided the amplitude of the time-dependent variations in the input is sufficiently small.

In the present paper, we investigate systematically the properties of this transfer function as a function of the connectivity of the network, for networks composed of two-populations, one excitatory (E), one inhibitory (I). Such networks have been extremely popular as simplified models of local networks in neocortex (Wilson and Cowan, 1972; van Vreeswijk and Sompolinsky, 1996; Amit and Brunel, 1997), hippocampus (Tsodyks et al., 1997), as well as other structures. Surprisingly however, extensive characterization of this transfer function is lacking, even for simple models. In fact, the goal of the present paper is to study this question systematically in two classes of models (1) “firing rate” models in which individual neurons are not described explicitly, rather all dynamical variables describe population quantities like average firing rates or average synaptic inputs (Wilson and Cowan, 1972; Ermentrout, 1998; Dayan and Abbott, 2001) (2) fully connected networks of leaky integrate-and-fire (LIF) neurons, in the presence of white noise (Brunel and Hansel, 2006). These two types of models have the advantage that the transfer function can be computed analytically, as a function of network parameters, provided the network is in a so-called “asynchronous state,” i.e., a state in which the population firing rate goes to a stable fixed point for time-independent inputs.

For each network model, we first characterize the region of stability of asynchronous states in the space of the parameters characterizing the coupling. Instabilities can arise through four qualitatively different mechanisms, which have been characterized extensively: (i) the “rate” instability (when the firing rate of the network can no longer be held at a given value, because of the excitatory feedback), present in all models; (ii) an oscillatory instability due to the E–I loop (Brunel and Wang, 2003), also present in both models; (iii) an oscillatory instability due to the delays in the I–I connections (Brunel and Hakim, 1999, 2008; Brunel and Wang, 2003; Brunel and Hansel, 2006), present in both models, provided delays are present in I–I connections, and finally; (iv) instability due to single neuron resonances at integer multiples of their firing rates, when the noise level is sufficiently small (Abbott and van Vreeswijk, 1993; Hansel and Mato, 2003; Brunel and Hansel, 2006). The last type of instability is only present in networks of spiking neurons.

We then characterize regions with qualitatively different types of transfer functions, in the space of coupling parameters, in both models. Transfer functions are characterized both in terms of the behavior of the gain as a function of frequency (i.e., monotonically decaying as a function of frequency, or with a single or multiple resonant peaks), and in terms of the phase relationship between E and I cells, again as a function of frequency. We relate qualitatively different types of transfer function with nearby instabilities. Finally, we compare the behaviors of rate model and network of integrate-and-fire neurons, and point out similarities and differences between both models.

2 Models

We investigate the dynamics of a network composed of two-populations of neurons, one excitatory (E), the other inhibitory (I; see Figure 1). The network received average, time-dependent inputs μE,ext and μI,ext to E and I populations, respectively. All types of connections are potentially present in the model, with strengths JEE (E → E), JIE (E → I), JEI (I → E), JII (I → I). We consider in this paper two types of models: a rate model and a network of LIF neurons.

FIGURE 1

Figure 1. Schematic diagram of the network model. An open circle means that the connection is excitatory and a bar means that it is inhibitory.

where τE,I characterize the time constant of firing rate dynamics, while FE,I(x) are the static transfer functions that describe the relationship between a mean synaptic input and an average firing rate in stationary conditions. Various types of response functions FE,I(x) have been used in the literature. Here, in order to simplify the analysis, we take this function to be threshold-linear (see, e.g., Ben-Yishai et al., 1995; Hansel and Sompolinsky, 1998; Roxin et al., 2005)

In the simplest rate model [Eqs (1,2)], synaptic time constants are ignored. Another class of rate models assumes instantaneous firing rate dynamics, and dynamical variables are assumed to represent average synaptic currents (Ermentrout, 1998). Since synaptic time constants are known to play an important role in shaping the dynamics of the activity of neurons, we investigate a rate model in which both neuronal and synaptic time constants are present. Real synaptic currents are characterized by three qualitatively distinct time constants: a latency (interval between pre-synaptic spike and the beginning of the post-synaptic response), a rise time, and a decay time (see, e.g., Destexhe et al., 1998). The simplest way to implement such synaptic currents in a rate model is through two synaptic variables sab(t) and xab(t) (where a, b = E, I denote the post/pre-synaptic population, respectively) whose dynamics are described by the following equations

where Dab, τrab, τsab are the latencies, rise times, and decay times of synaptic currents from population b to population a (see, e.g., Brunel and Wang, 2003). In this model, Eqs (1,2) are replaced by

We also consider a version of the model where only the latency is described (Roxin et al., 2005). In this case, Eqs (4,5,6) are replaced by

This model corresponds to a special case of Eqs (4–6) with τrab = τsab = 0 for all a, b.

2.2 Network of Leaky Integrate-and-Fire Neurons

In this model E and I populations are composed of NE, NI LIF neurons (see, e.g., Tuckwell, 1988; Dayan and Abbott, 2001), with NE + NI = N. Each neuron i of population a is described by its membrane potential Vi(t), which obeys

where τma is the membrane time constant of neurons in population a = E, I, Vrest = −70 mV is the resting membrane potential, σηi(t) is a noise term modeling the variability in external inputs received by neurons in addition to the “signal” μa,ext(t), where σ is the standard deviation of the noise and ηi(t) is a white noise of unit variance, uncorrelated from neuron to neuron. Finally, sij(t) are individual synaptic currents from neuron j ∈ b = E, I to neuron i ∈ a = E, I, modeled as in the rate model with the help of an additional synaptic variable xij:

where again Dab, τrab, τsab are the latencies, rise times, and decay times of synaptic currents from population b to population a. The sum over k represents a sum over all spikes of pre-synaptic neuron j, that occur at times Note the factor τma in front of this sum, which is chosen such that (i) both x and s variables are dimensionless; (ii) the total charge of individual PSCs is independent of synaptic time constants. In this model, action potentials occur when the voltage crosses Vthr = −50 mV. At the time of the action potential, the voltage is reset at Vreset = −60 mV. For simplicity, no refractory period is included.

The instantaneous firing rate of each population rE,I(t) is given by the sum of the number of spikes of the corresponding population in a time window dt

In numerical simulations, rE,I(t) are computed from simulating a network of NE = 2000 and NI = 500 neurons, except in Figure 9 where NE = 8000 and NI = 2000, using a Runge–Kutta second order method with a time step dt = 0.01 ms. In the mean-field analysis, one takes the limit N → ∞, dt → 0 (see below).

2.3 External Inputs and Response of the Network

In this paper, we consider external inputs of the form

where μa0 is a constant external drive, and μa1 is the amplitude of a sinusoidal modulation in the average input at a frequency ω = 2πf.

We assume μa1 is sufficiently small so that the response of the network firing rate is linear in μa1. If this assumption is satisfied, then the firing rate response is given by

where ra0 is a constant firing rate in absence of time-dependent modulation in the input (in other words, the background activity of the network), while ra1(ω) describes how the instantaneous firing rate is modulated by a sinusoidal input at frequency ω. It is a complex function which describes the gain and the phase of the modulation at frequency ω. In the following, we call ra0 the static transfer function of the network, while ra1(ω) is the dynamical transfer function of the network. The motivation for computing the response to sinusoidal inputs is twofold: (1) one can then reconstruct the response of the network to arbitrary-time dependent inputs; (2) oscillations are observed in many regions of the central nervous system, and are therefore a common type of input for post-synaptic structures.

3 Analytical Methods

The analysis of the network response is performed in three steps. (1) we compute the background activity as a function of the mean drive (static transfer function); (2) we check whether background activity with a constant firing rate is stable; (3) we compute the dynamical response of the network. The outcome of the analysis is therefore both the region of stability of states with time-independent firing rate (i.e., asynchronous states in networks of spiking neurons) in the space of parameters characterizing the network, and the dynamical transfer function as a function of such parameters. We now briefly outline how these three steps are performed in practice.

3.1 Static Transfer Function

3.1.1 Rate models

In the rate model, with constant external inputs (μa1 = 0 for a = E, I), the stationary firing rates are given by

where Ia0 is the synaptic input received by population a. In this paper, we only consider the situation where both input currents are suprathreshold (Ia0 > 0 for both a = E, I).

When Ia0 > 0 for both a = E, I, we have

It is sometimes useful to rewrite both equations as follows

where X and Y measure the strength of feed-forward and feedback inhibition, respectively:

and α0 is the ratio of the external inputs to the I population, divided by the external inputs to the E population, α0 = μI0/μE0. Note that in this model X should be smaller than one in order for the E rate to be positive.

3.1.2 LIF network

In a network of LIF neurons, the steady-state firing rates are given by

The basic idea is to add a small perturbation to the steady-state firing rate,

and to compute the eigenvalues λ of the system by solving the resulting equations at first order in δras. Instabilities arise when at least one eigenvalue has positive real part.

In the simplest rate model [Eqs (1,2)], the eigenvalues can be computed explicitly as a solution to the following quadratic equation:

As soon as delays are included, the eigenvalues no longer have simple analytical expressions. In order to compute the eigenvalues, we use the following strategy. The first step is, for a given perturbation in the rate δrb(t)eλt, to compute the resulting perturbation in the average synaptic variables

Figure 2. Neuronal (RLIF) and synaptic (S) transfer functions. (A) Modulus (amplitude) of RLIF as a function of frequency, for several values of the mean firing rate r0 and the amplitude of the noise σ. RLIF| is proportional to the slope of the static transfer function at low frequencies, and decays as at high frequency. It shows resonances at integer multiples of r0 at high firing rates and low noise. (B) Argument (phase shift) of the RLIF as a function of frequency, for several values of the mean firing rate r0 and the amplitude of the noise σ. (C) Modulus of S, for τs = 2 ms, τr = 0 ms (full line) and τs = 2 ms, τr = 1 ms (dashed line). Note that the modulus is independent of the delay. It decreases at high frequency as 1/(fτs) when the rise time is zero, and 1/(f2τsτr) when both rise and decay times are non-zero. (D) Argument of S, for three parameter sets: D = 0 ms, τs = 2 ms, τr = 0 ms (full line); D = 0 ms, τs = 2 ms, τr = 1 ms (dashed line); and D = 1 ms, τs = 2 ms, τr = 1 ms (dotted line).

From the perturbation in the individual synaptic variables, we can easily compute the perturbation in total average currents to neurons in a population

The second step is, from a perturbation in total average currents, to compute the perturbation in the instantaneous firing rate,

where M is a confluent hypergeometric function (Abramowitz and Stegun, 1970). The function RLIF is shown in Figure 2 for λ = iω. This figure shows some of the main qualitative features of RLIF: (i) when noise is large, it behaves as a low-pass filter, with a decay at high frequency in while resonances are present when noise is small; (ii) the gain depends markedly on both noise level, and mean firing rate.

The last step is to insert Eqs (27,29) in Eq. (30), and to divide both sides by eλt. Eliminating δras, we find that the eigenvalues λ should satisfy

where

describes the net effect of the firing rate of population b on population a. Note that the three terms in the r.h.s. of Eq. (34) potentially lead to three distinct instabilities, mediated by E → E, I → I, or the E → I → E feedback loop, respectively.

Note that in the simplest rate model, Sab(λ) = 1, Ra(λ) is given by Eq. (31), and Eq. (34) reduces to Eq. (26).

Equation (34) is solved numerically, using a globally convergent Newton routine which tracks the solutions λ = α + iβ to this equation, starting from many different initial conditions. Models are considered stable only if the routine converges toward negative αs for any initial condition.

3.3 Dynamical Transfer Function

Once linear stability of the background state has been checked, we can calculate the response of the system when the input has an additional oscillatory component. The oscillatory response of population firing rates is given by

where in the above equations Aab is given by Eq. (35) in which λ = iω, and Ra is given by Eq. (31) for the rate models, Eq. (32) for the network of integrate-and-fire neurons, in which λ = iω.

4 Results

4.1 One-Population Model

It is useful to start by summarizing the dynamics of a network of a single population of neurons, before moving to the two-population network. In a rate model, the stationary firing rate r0 is given by

where μ0 is the external input, and J is the recurrent coupling (excitatory if J > 0, inhibitory if J < 0). In a network of LIF neurons, the stationary firing rate r0 is given by

Linear stability analysis of the system shows that the eigenvalues λ satisfy

where R(λ) is given by Eq. (31/32) for rate model/network of LIF neurons, respectively, and S(λ) = 1 for the simplest rate model, and is otherwise given by Eq. (28).

When the asynchronous state is stable, the dynamical transfer function is

In an excitatory network (J > 0), the asynchronous state is stable provided J < Jc where Jc = 1 (rate model), Jc = 1/Φ′(Φ−1(r0)) (LIF network). When J = Jc, one of the eigenvalues becomes equal to zero, signaling the onset of a saddle-node bifurcation. This is the much-studied rate instability: because of the excitatory feedback, the network cannot maintain its rate at a value r0. Rather, the rate either decreases to a smaller value (typically close to zero), or increases to a larger value (typically close to maximal firing rate). When J < Jc, the effect of excitatory coupling is to amplify the response at low frequencies (the much-studied amplification effect, see Figure 3, red curves, and Seung, 2003), while at the same time slowing down the dynamics (through a reduction of the cut-off frequency, see, e.g., Seung, 2003). This effect is present in all models, independently of synaptic characteristics.

In the simplest inhibitory network (no synaptic time constants), the asynchronous state is always stable. Introducing a delay in synaptic transmission leads to an oscillatory instability (through a Hopf bifurcation) at a frequency ω. In the rate model in which no rise and decay times are present, the frequency ω satisfies the equation tan(ωD) = −τω, whose solution is ω ∼ π/2D when D ≪ τ (Roxin et al., 2005; Brunel and Hakim, 2008). This instability leads to an oscillation due to delayed inhibitory feedback when |J| > |Jc|, which has been characterized in many different inhibitory network models (see Brunel and Hakim, 2008 for a review). When |J| < |Jc|, the effect of inhibitory coupling is to suppress the response at low frequencies, which leads to an increase in the cut-off frequency (see, e.g., Seung, 2003). In the absence of delays, the response of the network is basically the response of a low-pass filter. In the presence of a delay, the response shows a resonance at a frequency close to the frequency of the oscillation at the Hopf bifurcation (see Figures 3B,D, blue curves).

4.2 Two-Population Rate Models

4.2.1 Stability of the asynchronous state in the simplest rate model

In the simplest rate model, the eigenvalues of the linear stability matrix are given by Eq. (26). From this equation it is easy to compute the conditions the coupling parameters must satisfy on the two possible types of instability lines.

A saddle-node bifurcation occurs when λ = 0; this happens whenever

or equivalently 1 + Y = 0. This instability corresponds to the rate instability induced by excitatory feedback onto the excitatory population, which is already present in the one-population model. The presence of feedback inhibition (as measured by the parameter Y) stabilizes the network with respect to this instability.

A Hopf bifurcation occurs when λ = iω this happens whenever

This bifurcation is a new feature of the two-population model; it depends both on E–E feedback and the E–I feedback loop. The two bifurcation lines are shown in Figure 4.

FIGURE 4

Figure 4. Qualitative behaviors of the network transfer function in the simplest two-population rate model. (A) Behaviors of the amplitude of the E network in the JEIJIE − JII plane, for three different values of JEE (indicated at the top of the corresponding panels). Red/blue areas show regions in which the asynchronous state is unstable due to the rate/Hopf instability, respectively. Above the dashed blue line, eigenvalues are real; below the line, they have a non-zero imaginary part. Green dashed/dotted lines indicate the boundary of the region in which the transfer function is of the low-pass type (above the lines) or resonant (below the lines), for different values of X (indicated on left panel). (B) Amplitude vs frequency curves, for two representative examples [+: JEE = 1.5, JEIJIE = 6, JII = 7, X = 0.5, *: JEE = 1.5, JEIJIE = 7, JII = 1, X = 0.5, both indicated by the corresponding symbols in the middle panel of (A)]. Both panels show the amplitude of the modulation of the E (red) and I (blue) populations. (C) Behaviors of the phase shift between E and I populations. Same conventions as in A except that black dashed/dotted lines indicate the boundary of the region in which the I population leads the E population at all frequencies (above the lines), or the opposite (below the lines), again for different values of X (indicated on left panel). (D) Phase vs frequency curves for two representative examples [+: JEE = 1.5, JEIJIE = 6, JII = 8, X = 0.5. *: JEE = 1.5, JEIJIE = 8, JII = 2, X = 0.5, again indicated by the corresponding symbols in the middle panel of (C)]. Both panels show the phase shift of the E (red) and I (blue) populations with respect to the input, as well as the phase shift between E and I populations (black). In all cases, τE = τI = 10 ms.

4.2.2 Dynamical transfer function in the simplest rate model

For a two-population rate model, the dynamical transfer function is given by Eqs (36,37), which can be rewritten for the simplest rate model as

Where is the effective time constant of the population a,

In writing Eqs (40,41), we have assumed that α1 ≡ μI1/μE1 is equal to α0 = μI0/μE0. This would be the case if both E and I networks receive inputs from the same structures. If E and I networks receive different types of inputs, then we could in principle have α1 ≠ α0. In this case, we would need to define two effective feed-forward parameters X0 (strength of “static” feed-forward inhibition) and X1 (strength of “dynamic” feed-forward inhibition), proportional to α0 and α1, respectively.

We focus here on two quantities: the amplitude of the excitatory transfer function |rEI(ω)|, and the phase shift between the two-populations, ΔΦEI(ω). In this model, it is easy to show (see Appendix) that the amplitude of the transfer function can have only two qualitatively different behaviors: (i) a monotonically decaying function of the frequency ω (“low-pass” behavior – see, e.g., Figure 4B, panel with the cross); or (ii) a function that first increases with frequency, peaks at a finite frequency, and then decays toward zero at high frequencies (see Figure 4B, panel with the star). Note that in the rate model the amplitude of the response does not depend on ra0. Therefore the examples in the Figures 4, 5, and 6 are shown with amplitudes normalized by ra1(f = 0) and not ra0. The regions where both types of behaviors occur are shown in the plane of couplings involving the inhibitory network (self-inhibitory coupling JII vs strength of EI loop JEIJIE), for different values of JEE, in Figure 4A. We see that strong values of JII tend to favor low-pass behavior, while strong values of JEIJIE tend to favor resonant behavior. Close to the saddle-node bifurcation, the transfer function is low-pass, with a zero-frequency response that diverges on the bifurcation line. Close to the Hopf bifurcation, the transfer function is resonant, with a peak that diverges on the bifurcation line, at the frequency of the oscillation that appears on the line.

FIGURE 5

Figure 5. Qualitative behaviors of the rate model with delays, for τE = τI = 10 ms, JEE = 1. 5. (A,B)DE = 2 ms, DI = 1 ms, (C,D)DE = DI = 1.5 ms; (E,F)DE = 1 ms, DI = 2 ms. (A,C,E) frequency of the oscillatory instability on the bifurcation lines, as a function of JEIJIE. (B,D,F) behaviors of the gain of the E network in the JEIJIE − JII plane for three different values of X (indicated at the top of the corresponding panels). White areas show regions in which the asynchronous state is unstable due to the rate or Hopf instabilities, shown as black lines. The stable region is subdivided in regions of different colors, indicating qualitatively different behaviors of the E transfer function. Colors indicate the number of extrema at non-zero frequencies of the amplitude vs frequency curve (0: black, 1: green, 2: red, 3: orange, 4: yellow). (G) Representative transfer functions (top: amplitude; bottom: phase) for each of the regions. From left to right, JEIJIE = 5, 18, 23, 80, 19; JII = 5, 4, 15.5, 14, 7. Symbols on the top left of the amplitude plots refer to the location of the network parameters in (B,F). Black dots on top of the red curves indicate the location of the extrema.

The boundary separating the two regions depends on the strength of feed-forward inhibition X: the larger X, the more pronounced the resonant behavior. This is because feed-forward inhibition tends to suppress inputs at low frequencies. In fact, feed-forward inhibition can create a resonance even in the absence of feedback E–I coupling (see Figure 4A, case where JEE = 0). In the limit X → 1, the resonant behavior becomes present in the whole stable region. Note also that the presence of a resonance is not equivalent to the presence of complex eigenvalues of the linear stability analysis, as happens in a two-variable single neuron model (Richardson et al., 2003): the network can be resonant in the absence of complex eigenvalues, and conversely, the network can be low-pass in the presence of such complex eigenvalues.

Another quantity of interest is the phase shift between excitatory and inhibitory populations. In this model, one can show (see Appendix for details) that again only two qualitatively different behaviors can occur: (i) the excitatory population leads the inhibitory population at all frequencies f (see, e.g., Figure 4D, panel with the cross); or (ii) the opposite occurs and the inhibitory population leads the excitatory population at all f (see, e.g., Figure 4D, panel with the star). The regions where both types of behavior occur are shown in the plane of couplings involving the inhibitory network (self-inhibitory coupling JII vs strength of EI loop JEIJIE), for different values of JEE, in Figure 4C. Close to the rate instability, inhibition tends to lead excitation, while the opposite occurs close to the Hopf bifurcation. Note that for all parameters, the phase shift between the two-populations tends to zero in both f → 0 and f → ∞ limits, and that it reaches a maximum at some finite value. In the resonant region, this value is close to the location of the resonant peak (compare Figures 4B,D, panels with the star).

4.3 Rate Model with Non-Zero Synaptic Time Constants

How are the dynamics modified when the synaptic time scales are taken into account in the rate model? We show here that new behaviors emerge, due to the presence of delays in synaptic interactions.

4.3.1 Stability of the asynchronous state

The stability of the asynchronous state can be evaluated by computing the eigenvalues, i.e., the solutions of Eqs (34,35,28,31). We first discuss the stability analysis when synaptic interactions are simply described by a delay (i.e., rise and decay times are equal to zero), and both E and I populations have the same rate time constant. In fact, all the qualitative features of the general case (i.e., with strictly positive rise times and decay times, and/or with τE ≠ τI) are already present in this simplified model.

The main qualitative changes that occur when delays are introduced are as follows: (1) there is a new oscillatory instability induced by I–I connectivity, which is reached when JII becomes strong enough (Brunel and Hakim, 1999; Roxin et al., 2005); (2) the oscillatory instability induced by the E–I feedback loop is now present at all values of JEE. As a result of the new oscillatory instability, the space of parameters for which the asynchronous state is stable is now bounded for any value of JEE (see Figure 5). The number of distinct bifurcation lines that form the boundary of the region of stability can vary between one and three, depending on parameters. First, the rate instability line is only present when JEE > 1, as in the simplest rate model without synaptic time constant. Second, there can be either one or two distinct oscillatory instability lines, depending on the relative values of the excitatory and inhibitory delays, DE and DI.

When DE > DI, there exist two distinct instability lines, characterizing the I–I and E–I scenarios; these two lines cross at some point in the JII − JEIJIE plane (see Figures 5A,B). The frequencies of the corresponding instabilities remain well separated. The frequency of the I–I instability is generically higher than the frequency of the E–I instability, since the I–I induced-oscillations involve only two synapses from peak to peak (increase in I rate leads to a decrease in I rate due to I feedback, which in turn leads to the next increase in I rate), while E–I induced-oscillations involve four synapses (E rate increases, leading to an increase in I rate, leading to a decrease in E rate, then decrease in I rate, followed by the next increase in E rate). For example, for the parameters of Figure 5A, I–I instability leads to frequencies in the 150 to 200-Hz range, while E–I instability leads to frequencies which are in the 10 to 80-Hz frequency range. When DE = DI, the two instability lines meet at a cusp; at this point the frequencies of the two instabilities become equal (Figure 5C). Finally, when DE < DI there is a single oscillatory instability line, smoothly connecting the “pure” I–I and “pure” E–I scenarios in the JII − JEIJIE plane (Figure 5F). In this case the frequency decreases smoothly along the line as the JII/JEIJIE ratio is decreased. As a result, the region of stability of the asynchronous state exhibits either a triangular shape, when DE ≥ DI (see Figures 5B,F) or a D shape, when DE < DI (see Figure 5F).

4.3.2 Dynamical transfer function

The next step is to compute the dynamical transfer function, Eqs (36,37). We show in Figure 5 how the qualitative shape of this transfer function depends synaptic parameters and on the connectivity strengths. In the vicinity of the rate instability, the transfer function has either low-pass characteristics (black regions in Figure 5, see a representative transfer function in Figure 5G), or shows a trough followed by a resonant peak, when the system is close to the I–I instability (red region in Figure 5, see a representative transfer function in Figure 5G). In the vicinity of the Hopf bifurcation lines, the transfer function always exhibits resonant peaks, but in qualitatively different ways. There can be a single resonant peak – the typical situation when the stability region is D-shaped, i.e., there is a single bifurcation line (see green, red, and yellow regions in Figure 5, and the corresponding transfer functions in Figure 5G). Note that in the yellow region, there is a second peak in the transfer function at a frequency which is larger than the one of the resonant peak, but its magnitude is always small compared to the amplitude of the other peak. Alternatively, there can be two separate resonant peaks, typically when the stability region is triangular, i.e., with two separate bifurcation lines, with two separate frequencies (see orange regions in Figure 5B, and the transfer function in Figure 5G), and the network is close to the intersection between the two lines. The frequency of the resonant peaks is close to the frequency of the corresponding instabilities. For parameters of Figure 5, the resonant frequency is typically in the gamma range when the system is close to the E–I bifurcation line, while it is typically very high (more than 100 Hz) when the system is close to the I–I line. As in the simplest rate model, the boundary between the different regions also depends on the strength of feed-forward inhibition: in particular, the pure low-pass region tends to disappear as X increases, since a high value of X decreases the response at low frequencies.

The transfer functions can also be classified according to the behavior of the phase shift between excitatory and inhibitory populations. The main new behavior introduced by the presence of synaptic delays is that it is now possible for the phase shift to change sign when the frequency is increased. This leads to a larger number of qualitatively different behaviors for the phase shift. Some of the typical behaviors of the phase shift are represented in Figure 5G. Note in particular that resonances are associated with sharp changes in the phase shift. Close to the E–I instability line, the E population leads the I population at the corresponding resonant frequency.

The behavior of a model with τE ≠ τI is qualitatively very similar to the one of the model with τE = τI. The main difference is that the transition between the 1 Hopf and 2 Hopf bifurcation behaviors no longer occurs at DE = DI, but for slightly different values of the delays. For example, for τE = 20 ms, τI = 10 ms, and DE = 1.5 ms, the transition between both behaviors occurs at DI ∼ 1.55 ms.

4.3.3 Behavior of the model with finite rise and decay times

A similar analysis can be performed for the rate model with finite delay, rise, and decay times. The qualitative behaviors of the transfer function in the region of stability of the asynchronous state are shown in Figure 6. This figure shows that the same types of behaviors are observed as in Figure 5. At the quantitative level, the frequencies of the instabilities are reduced by the presence of finite rise/decay times (compare Figures 6A,C with Figures 5A,E). This decrease in frequency is stronger for the I–I instability. Consequently, the two possible resonant frequencies are now much closer than in the case in which only delays are present, which decreases the size of the region in which both resonances are observed (compare orange regions in Figures 5 and 6). Another quantitative difference is that phase differences between E and I populations tend to be smaller.

4.4 Network of LIF Neurons

The same analysis can be performed for a network of LIF neurons, as outlined in Section 3.2. Interestingly, the stability analysis of the asynchronous state revealed the same types of instabilities as in the rate model with synaptic time constants, i.e., a rate instability (mostly governed by E–E interactions), and two types of oscillatory instabilities (governed by E–I or I–I interactions). As in the rate model, the two types of oscillatory instabilities merge smoothly in the space of synaptic connection strengths if DE < DI, while they remain disjoint if DE = DI. There are however a few qualitative differences between the two models:

• The high frequency behavior of the LIF model is different from the rate model: the amplitude of the response decays as , and the phase shifts go to π/4, in the LIF network, as compared to a 1/f decay with a phase shift of π/2 in the rate model (see Section 5)

• In the LIF network, the single neuron transfer function depends both on its average inputs (and therefore its average firing rate) and on the level of noise. As a consequence, the network transfer function changes as a function of both mean and variance of the external inputs.

In this section we first show how the qualitative behaviors of the LIF network compare with those of the rate model (at fixed mean firing rate and level of noise), and then show how external inputs can qualitatively change the network transfer function in the LIF network. This property is not present in a rate model with threshold-linear static transfer function, but can be reproduced with a suitable non-linearity in the static transfer function.

4.4.1 Comparison between LIF network and rate model

Figure 7 shows the behaviors of the transfer function of the LIF network in strong noise regime (σ = 5 mV). To make an explicit comparison with the rate model, we choose to fix the firing rates of E and I populations to specific values (10 Hz and 30 Hz, respectively), independently of the matrix of coupling. This means that on any point of the phase diagrams showed in Figure 7 the external inputs must be computed such that the network has the required mean firing rates. We also need to generalize the definition of the feed-forward and feedback inhibition parameters to the LIF network. A natural generalization is

In Eqs (42,43), we see that these parameters are now governed by “effective coupling strengths,” where the coupling strength Jab is multiplied by the gain of the transfer function of the pre-synaptic population at the corresponding firing rate.

To compare the LIF network with the rate model, we then modulated independently the sinusoidal inputs to the E and I populations, but with a ratio between amplitudes taken in order to keep the feed-forward inhibition parameter X constant in each plot. Synaptic parameters was taken to be identical to the rate model. Figure 7 shows that, qualitatively, the LIF network behaves essentially as the rate model. There are of course quantitative differences due to the fact the transfer function of the LIF neuron is different from that of the rate model. Importantly, the high frequency behavior is different, with a decay of the gain, and a phase shift of 45° at high frequencies in the LIF model, as already mentioned above.

4.4.2 Dependence of the network transfer function on external inputs

In the LIF model, the single-cell transfer function RLIF depends on the inputs to the network. In particular, changing the overall firing rate alters the static gain of the transfer function Φ, which changes the effective coupling strengths and in turn the effective feed-forward and feedback inhibition parameters.

Figure 8 shows an example of a network where an increased input to both populations led to qualitative changes of the network transfer function. The coupling strengths were chosen such that when the external inputs to both populations increase, the firing rates of both populations first increase. Then, inhibition becomes strong enough to reduce excitatory rates even with the increased external inputs. Interestingly, this leads to a feed-forward parameter X larger than one. The network transfer function goes through a number of changes as the external inputs are varied: at low input strength/low rates, the transfer function is low-pass, due to the fact that the gain of single neuron f–I curves is small, and hence effective coupling strengths are small. When the inputs increase, the slopes of the transfer functions increase, and therefore the effective couplings increase. This leads to a resonant transfer function, due to feedback between E and I populations. Finally, when the inputs increase further, the E rate reaches a peak and then decreases with increased inputs. This leads to a major change of the response to external inputs: the E and I populations are now in antiphase for slowly varying inputs (see Figure 8D). This is due to the fact that X > 1. When X = 1, then the phase shift between the two-populations is 90°, as shown by Figure 8D. In both cases, the phase shift between the two-populations decreases with frequency, and goes to zero in the high frequency limit.

A qualitative change in network transfer function with increasing external input is a feature that is not present in a rate model with threshold-linear static transfer function, but it can be obtained using a suitable non-linear transfer function. In particular, Figure 8D shows when one uses a rate model with the static transfer function of the LIF neuron, Eq. (24), one obtains network transfer functions which are very similar to the LIF network.

4.4.3 Dependence of the network transfer function on external noise

The single-cell transfer function also depends on the magnitude of external noise, as shown by Figure 2. In particular, as was characterized in a number of previous studies (Knight, 1972; Brunel et al., 2001; Fourcaud and Brunel, 2002; Fourcaud-Trocmé et al., 2003), the single-cell transfer function exhibits pronounced resonances at integer multiples of the single-cell firing rate at low noise, while these resonances disappear at high noise and the transfer function becomes low-pass.

This qualitative change in the single-cell transfer function affects the network behavior, as seen in Figure 9, which shows the qualitative behaviors of the network transfer function as a function of a parameter g measuring the ratio between inhibitory and excitatory connection strengths (i.e., we choose JII = JEI = gJEE = gJIE), and the strength of external inputs. Parameters are chosen such that the network can potentially have a double resonance, due to E–I feedback coupling and I–I coupling, respectively (the orange region on the plots). The top graphs show the “phase diagram” in the high noise regime (σ = 5 mV) that has been explored so far. We see that the asynchronous state is stable in a finite range of external inputs. As the external inputs increase, the network transfer function develops a single or several resonances, depending on the value of g. Above a critical value of the external inputs the asynchronous state destabilizes due to a Hopf bifurcation, as in Brunel (2000).

The other two graphs show how the picture is modified when the level of external noise decreases. There are two main effects: first, the region of stability of the asynchronous state becomes smaller, as expected from previous studies in purely inhibitory networks (Brunel and Hakim, 1999; Brunel and Hansel, 2006); second, qualitatively new behaviors of the network transfer function emerge. In particular, the network can develop a new resonance which is associated with the single-cell resonance (see Figures 9D,F).

5 Discussion

In this paper, we have investigated how the transfer function of an E–I network is controlled by its connectivity and the temporal parameters of its synapses. We have shown that a two-population rate model, as well as a fully connected network of excitatory and inhibitory integrate-and-fire neurons, exhibit many qualitatively different behaviors, from low-pass to multiple resonances. Resonance phenomena appear generically in both models provided inhibition is sufficiently strong. With realistic synaptic currents, oscillatory instabilities can arise due both to the E–I feedback loop and to the I–I connectivity. Close to these instabilities, a single or double resonance can appear, depending on synaptic parameters. We have also shown that for a given level of external input, and in a high noise regime, the behaviors of the network transfer function of LIF networks is qualitatively similar to the rate models with identical synaptic dynamics. However, an important difference is that in the network of LIF neurons the transfer function depends on the level of external inputs, and external noise. In a rate model with a threshold-linear transfer function, the “single-cell” dynamical transfer function is independent of the input. If a sigmoidal static transfer function had been used instead of the threshold-linear one, then the gain at zero-frequency (i.e., the slope of the f–I curve) would depend on the input, which could also lead to qualitative changes of the network transfer function as a function of the input. However, the LIF neuron has the additional feature that the single-cell dynamical transfer function can exhibit qualitative changes as a function of its inputs, such as the appearance of a resonance at single-cell firing frequency when the noise is sufficiently weak.

Input-dependent changes in the network transfer function could underlie several prominent phenomena observed in vivo in different areas. In visual cortex, the power spectrum of the LFP typically decays monotonically at low contrast, while a plateau or even a peak in the gamma frequency emerges at high contrast (Henrie and Shapley, 2005). Input-dependent changes of the transfer function at gamma frequency (due to the emergence of a resonance) has been modeled in a framework similar to ours by Kang et al. (2010). Such input-dependent changes could underlie the high information content of the gamma range in the LFP (Belitski et al., 2008), as shown by Mazzoni et al. (2008). It could also underlie the phenomenon of cross-frequency coupling, through which the power of a fast oscillation is strongly modulated by the phase of a slower oscillation, as is the case in the hippocampus (Bragin et al., 1995) as well as in primary sensory cortices (Whittingstall and Logothetis, 2009; Mazzoni et al., 2010).

This study could be generalized in a number of ways. First, it could be generalized to more realistic single-cell models such as the exponential integrate-and-fire (EIF) neuron (Fourcaud-Trocmé et al., 2003), since an efficient method has been introduced to compute both the static and dynamic single-cell transfer function for this model (Richardson, 2007). Interestingly, the transfer function of the EIF model has a high frequency behavior which is qualitatively identical to the rate model (a 1/f decay at high frequency, associated with a phase shift of 90°). This led Ostojic and Brunel (2011) to propose a reduction of the EIF model to an (adaptive) rate model, whose parameters (time constant and static transfer function) can be computed analytically from the parameters of the underlying EIF model. From these results, we expect the dynamics of networks of EIF neurons to be even closer to rate models, compared to LIF networks, provided of course the noise is sufficiently strong.

A second generalization would be to study networks of more than two-populations. It would be of interest, for example, to understand how the presence of functionally different classes of interneurons affect network dynamics. Finally, it would also be of interest to study the transfer function of randomly connected networks (Brunel and Hakim, 1999; Brunel, 2000), or of networks with more structured connectivity (Ben-Yishai et al., 1995; Roxin et al., 2005).

In this paper, we have computed how the population firing rate reacts to time-dependent inputs. Experimentally, the firing rate of small populations of neurons can be observed through MUA. Another experimental observable is the LFP. In cortex and hippocampus, the LFP is assumed to be related to the average synaptic inputs received by pyramidal cells of the main network. Given a model of how the LFP depends on synaptic currents (e.g., Mazzoni et al., 2008), one could compute the dynamics of the LFP signal using the methods used in this paper.

Finally, this work offers insight into the relationship between network connectivity and experimental observables such as phase shifts of different populations in the network. It could be used to predict network connectivity from observed patterns of resonances and phase shifts, in areas where this information is available.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Boris Barbour and Vincent Hakim for useful discussions at the initial stage of the project, and for a critical reading of the manuscript. We also thank David Hansel for useful suggestions.

Hansel, D., and Mato, G. (2003). Asynchronous states and the emergence of synchrony in large networks of interacting excitatory and inhibitory neurons. Neural Comput. 15, 1–56. doi: 10.1371/journal.pcbi.1000239