Oscillations and Irregular Emission in Networks of Linear Spiking ...
Journal of Computational Neuroscience 11, 249–261, 2001 c 2002 Kluwer Academic Publishers. Manufactured in The Netherlands.
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons GIANLUIGI MONGILLO AND DANIEL J. AMIT Department of Physics, University of Rome La Sapienza, Ple. Aldo Moro, 00185, Rome, Italy [email protected]
Received February 12, 2001; Revised July 25, 2001; Accepted October 12, 2001 Action Editor: Roger Traub
Abstract. The dynamics of a network of randomly connected inhibitory linear integrate and fire (LIF) neurons (with a floor for the depolarization), in the presence of stochastic external afferent input, is considered in various parameter regimes of the neurons and of the network. Applying a technique recently introduced by Brunel and Hakim, we classify the regimes in which such a network has stable stationary states and in which spike emission rates oscillate. In the vicinity of the bifurcation line, the oscillation frequency and its amplitude are computed and compared with simulations. As for leaky IF neurons, the space of parameters can be compactified into two. Yet despite significant technical differences between the two models, related to both the different dynamics of the depolarization as well as to the different boundary conditions, the qualitative behavior is rather similar. The significance of LIF neurons and of the differences with leaky IF neurons is discussed. Keywords: spiking neuron, recurrent networks, collective oscillations.
1.
Introduction
Spiking neurons that integrate linearly (rather than exponentially) and have a reflective floor for their depolarization (LIF), present several attractive features: they are natural candidates for Very Large Scale Integration (VLSI) implementations (Mead, 1989; Fusi et al., 2001); the single-cell spike dynamics is relatively simple (Mattia, 1997; Fusi and Mattia, 1999); and due to their depolarization floor, networks of such neurons have both low (negative drift) and high spike rates (positive drift) (Mattia, 1997; Fusi and Mattia, 1999) and hence coexistence of spontaneous and selective activity states (Mattia, 1997). Recently (Rauch et al., 2001), it has been observed that the emission rate of a single cell to a gaussian injected current is rather well described, with few free parameters, by the transduction function of the LIF neuron, in a rather wide range of low rates. Both facts naturally raise the question as to whether networks of spiking LIF neurons would exhibit collec-
tive behavior that is qualitatively similar to that of the more common RC neurons. Historically, such neurons have been studied in the absence of a depolarization floor (Gerstein and Mandelbrot, 1964; Knight, 1972). In that case they emit spikes regularly only for positive drift; their characteristics are rather transparent and not very rich in phenomenology. In a recent study Brunel and Hakim (1999) have paved the way for a detailed investigation of the dynamic characteristics of rather complex networks of spiking neurons. Following on the empirical observation (Amit and Brunel, 1997a) that spike dynamics in a network is dominated mostly by the evolution of the global activity of the network, Brunel and Hakim have set up a perturbative study of the network’s spike dynamics around its mean-field stationary states. The approach constitutes a powerful tool that allows an effective classification, in network parameter space, of activity states of different types, including stationary
250
Mongillo and Amit
states of irregular spiking, oscillatory states for the global population with irregularly firing neurons, and synchronous states. The initial study considered networks constituted exclusively of inhibitory neurons driven by external currents. It was then extended to networks with both excitatory and inhibitory IF neurons, which can sustain spontaneous activity and selective delay activity following learning (Brunel, 2001a). Beyond the clarification of the dependence of different network operating states on critical network parameters, the technique developed allows also the investigation of the reaction of a network of spiking neurons, in a stationary state, to an external perturbation, such as the change in external afferents or the presentation of a stimulus. Here we apply this tool to a network of inhibitory spiking LIF neurons. The neurons in the network are connected randomly, and external currents are stochastic. This study opens the way for the analysis of the domains of different collective behavior of networks with LIF neurons, which is relevant to the correspondence between VLSI networks and networks of more biologically plausible RC neurons. Given the complexity of the system it is quite rewarding that the theoretical technique, while somewhat complicated, allows for detailed quantitative predictions of parameter regions for states of different dynamical characteristics (stationary, oscillatory), mean emission rates, and oscillation periods and amplitudes of the rates, where the network transits into oscillating behavior. All these are compared against simulations. 2. 2.1.
current tends to hyperpolarize the neuron. Emission in this case is due to fluctuations. We define the drift as µ ≡ −β + I , where I is the mean afferent current. When it is positive, firing is essentially regular, and the rate of emission is rather insensitive to fluctuations. The presence of the two firing regimes is a sufficient condition for networks of such neurons to have both spontaneous and working memory activity states (Brunel, 2001b). The RC neuron also has such two regimes of operation, depending on whether the mean afferent current is below or above threshold, respectively. We therefore compare (Fig. 1) the behavior of the two models by superposing the corresponding dynamical regimes. In other words, we map the negative drift region of the LIF neuron on the subthreshold regime of the RC neuron (ND regimes) and the complementary (SD regimes) on each other. To render the variables describing the two neurons of equal dimensions, we define µ∗ as µ for the LIF neuron and ( I − θ )/τ for the RC neuron.1 The correspondence of the noise 2 2 variable is σLIF to σ RC /τ , where τ is the membrane time constant of RC neuron. Note that if the reflecting barrier is absent, V can become arbitrarily negative,
Methods LIF Neuron and Stationary States of Networks
The evolution of the depolarization V (t) of the LIF neuron (Mead, 1989) is V˙ = −β + I (t),
(1)
where β is the linear loss, and I (t) is the afferent current charging the membrane. This equation is accompanied by the standard spike emission condition at an upper threshold θ, when V (t k ) = θ a spike is emitted and V is reset to → Vr . The VLSI-LIF neuron has a floor for its depolarization; V cannot be hyperpolarized below V = 0. When the depolarization is pushed to negative values, it remains at 0, a reflecting barrier. This allows spike emission also when the mean afferent
Figure 1. Comparison of current-to-rate transduction functions of LIF (full curves) and exponential (RC) neuron (dotted), for two values of the corresponding variances of the afferent current. The correspondence used: matching the subthreshold regime of the RC neuron (µRC < θ ) with the ND regime of the LIF neuron (see text). The plots exhibit the emission rates of the two neurons vs µ∗ , where µ∗ = (µRC − θ )/τ , for the RC neuron and is µ for the LIF neuron. The variances put in correspondence are σ 2 for the LIF neuron and σ 2 /τ for the RC neuron. Note the higher rates of the LIF neuron at corresponding negative values of the µ∗ and the consequent smaller gain of this neuron at µ∗ = 0. θ = 20 mV for both types of neurons, τRC = 20 ms.
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
and the LIF neuron operates exclusively in the positive drift regime (see, e.g., Fusi and Mattia, 1999). This is of little consequence for the exponential neuron. If the afferent current to the neuron is stationary, Gaussian, with mean µ, with variance σ 2 per unit time, and independent at different times, the spike emission rate is calculated to be (Fusi and Mattia, 1999) (µ, σ ) = 2
−1 σ 2 2µθ − 2µθ + e σ2 − 1 . 2µ2 σ 2
(2)
The significance of this result is that in a network of interacting spiking neurons with low, random connectivity and low rates, the spike trains afferent along different presynaptic channels often are independent to a good approximation. Moreover, the synaptic efficacies are rather small relative to the spike-emission threshold. Hence, the recurrent currents can be described as Gaussian, with a mean and width determined by the characteristics of functional subpopulations of neurons and by the average spike rates in these populations. In this way one arrives at a mean-field description of the average spike dynamics in the entire network (Amit and Brunel, 1997b; Mascaro and Amit, 1999). In particular, if the network is a homogeneous population of inhibitory neurons connected by synapses of connectivity C, of efficacy with average J and variance 2 J 2 , then the mean and variance of the recurrent afferent currents is µl (t) = −C J ν(t); σl2 (t) = C J 2 (1 + 2 )ν(t), (3) where ν(t) is the spike rate in the population. When external currents are added, with mean and variance 2 µext , σext , respectively, one has µ(t) = −β + µl (t) + µext (t); 2 σ 2 (t) = σl2 (t) + σext (t).
(4)
Given the hypotheses of mean-field, with mean and variance linear functions of the instantaneous rate, we can write the equation for the stationary states of the networks, those activity states that reproduce in output the rate appearing on input: ν = (µ(ν), σ 2 (ν)).
(5)
The solutions, ν0 , of this equation must be accompanied by a test of stability. A stability analysis, for the stationary states of a network of leaky IF neurons, has
251
been developed (Amit and Brunel, 1997a), based on dynamic equations for µ and σ in conditions of very slow variations for the underlying integrator. Both simulations (Van Vreeswijk et al., 1994; Gerstner, 1995) and experiment (Buzsaki et al., 1992; Traub et al., 1996) (for many additional references, see, e.g, Brunel and Hakim, 1999) indicate that networks of the type considered may oscillate at high frequency unrelated to the spike emission rate, even in ranges of parameters where the test predicts stability. Moreover, such an asymptotic regime, in the dynamics of the depolarization, does not exist for the LIF neuron. Hence the requirement for a finer tool of control over the dynamics of the network and of the stability of its stationary states. Following Brunel and Hakim, one finds that not only can one obtain a proper classification of the parameter regimes of stable and unstable (oscillatory) states, but the frequency and amplitude of the spike emission rates in the oscillatory states can be fully computed near the line of instability.
2.2.
Perturbation of Stationary States
To go beyond the study of rates in stationary states, one observes (Brunel and Hakim, 1999; Amit and Brunel, 1997b) that in randomly and sparsely connected networks, firing patterns of different neurons are essentially independent. The network dynamics can be therefore characterized by the instantaneous distribution of depolarizations in the neural population together with the spike emission probability. In networks in which there is more than one functionally distinct population one would have, of course, as many distribution functions and spike emission probabilities as distinct populations (see, e.g., Brunel, 2000a). The temporal evolution of the probability distribution of the depolarization (the Fokker-Plank equation) corresponding to a neuron whose dynamics is governed by Eq. (1) is σ 2 (t) ∂ 2 P ∂P ∂P = , − µ(t) 2 ∂t 2 ∂V ∂V
(6)
where µ and σ are given in Eqs. (3) and (4). Note that the rate ν entering the mean and variance at time t is the rate at t − δ, the rate a synaptic delay constant (δ) earlier. Like for the exponential IF neuron, the boundary condition at the spike emission threshold is P(θ, t) = 0. Moreover, the probability flux at θ gives
252
Mongillo and Amit
the emission probability at time t, ν(t), that is, 1 ∂P ν(t) = − σ 2 (t) (θ, t). 2 ∂V
where
At the reset potential (Vr = 0, for simplicity), the boundary condition is different from the one for the exponential IF neuron. For the linear neuron, there is a reflecting barrier (Fusi and Mattia, 1999.) No process can cross V = 0 from below. The only flux across this barrier is due to spike emission. Hence, 1 2 ∂P 1 ∂P σ (t) (0, t) − µ(t)P(0, t) = σ 2 (t) (θ, t). 2 ∂V 2 ∂V (8) In a stationary state, Eq. (6) for the stationary pdf, P0 (V ), becomes σ02 ∂ 2 P0 ∂ P0 − µ0 = 0, 2 ∂V 2 ∂V whose solution is (Mattia, 1997; Fusi and Mattia, 1999) P0 (V )
µ0 ν0 1 − exp −2 2 (θ − V ) = in V ∈ [0, θ ]. µ0 σ0 (9)
Imposing the normalization condition on P0 one obtains the transduction Eq. (2), for the stationary emission rate. Next we rescale the variables in Eq. (6) as follows: P=
σl2 CJ 2 (1 + 2 )ν0 = , 2 σ02 CJ (1 + 2 )ν0 + σext CJν0 µl B= =− . −β − CJν0 + µext µ0
(7)
2ν0 σ2 Q; V = 0 V¯ ; ν(t) = ν0 (1 + n(t)). µ0 µ0 (10)
A=
2
Moreover, τ ≡ σ02 /µ20 will play below the role of an intrinsic time constant, which is not present in the evolution of the depolarization: it determines the temporal scale on which there is a cross-over from small to large synaptic delays. The differences remain in the form of the operators L and in the boundary conditions. Note that A and B (in distinction to those of the IF neuron) are, respectively, the ratio of the variance and the mean of the local currents to the total variance and currents. A remains in the range [0, 1]. But B, positive definite for the IF neurons (= −µl /σl ), is no longer so in our case: µ0 will be indeed negative for negative drift. The boundary conditions for Q become ¯ t) = 0; Q(θ, (13) ∂Q 1 + n(t) ¯ t) = − (θ, 1 + An(t − δ) ∂ V¯ 1 ∂Q A ∂Q × − Q + n(t − δ) + BQ 2 ∂ V¯ 2 ∂ V¯ V¯ =0
1 = − [1 + n(t)], 2 where θ¯ = 2.3.
(14)
µ0 θ. σ02
Linear Stability Analysis
We write the deviations from the stationary state as The new variable are all dimensionless. This transformation leads Eq. (6) to an equation rather similar to the one of the exponential neuron, with only two parameters: τ
∂Q ∂t
1 ∂2 Q ∂Q = − + n(t − δ) 2 ¯ 2 ∂V ∂ V¯
Q(V¯ , t) = Q 0 (V¯ ) + Q 1 (V¯ , t) + O( 2 ) (15) n(t) = n 1 (t) + O( 2 ).
(16)
At first order, one has
∂Q A ∂2 Q +B 2 ¯ 2 ∂V ∂ V¯
(11)
≡ L[Q](V¯ , t) + n(t − δ)L AB [Q](V¯ , t),
(12)
τ
∂ Q1 = L[Q 1 ](V¯ , t) + n 1 (t − δ)L AB [Q 0 ](V¯ ), (17) ∂t
in which L and L AB are linear operators defined in Eq. (12). This leads (see, e.g., Appendix A) to an
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
253
eigenvalue equation of the form 1 −z 1 θ¯ z 1 (αz 2 − β)e −1 z2 − z1 2 1 ¯ z2 − (αz 1 − β)e−z2 θ − 1 + He−λδ/τ = − , 2 2 (18) where −λδ/τ α(λ) = (A + B) e λ β(λ) = −1 + Ae−λδ/τ + 2α (19) 1 ¯ ¯ H = B(1 − e−2θ ) − Ae−2θ 2 √ and z 1,2 = 1 ∓ 1 + 2λ with the positive (negative) sign corresponding to positive (negative) drift. 2.4.
Limit of Short Delays
For the linear neuron, for a wide range of parameters, σ02 is relatively large, ≈ 100 ms, while δ is a few ms. µ20 So δ/τ 1. In this limit √ the eigenvalues become very large, and z 1,2 = ∓ 2λ. In Eq. (18) one of the two exponentials in square brackets dominates, the one with z i θ¯ > 0, and the equation reduces to (Appendix A2) αz = β.
(20)
Note that θ¯ is positive (negative) for positive (negative) drift. Hence, dealing with the LIF neuron, we use Eq. (20), after having checked that the solutions for the eigenvalues agree very closely with those of the full Eq. (18). 3. 3.1.
Results Stability, Instability, Bifurcation
Using Eq. (20) one identifies the regions, in the space of parameters, in which the stationary state of the network is stable: all eigenvalues have a negative real part and is unstable: some eigenvalues have a positive real part. There are also bifurcation lines: those at which the eigenvalue with the largest real part is purely imaginary. To start with, the number of parameters is quite high: β, C, , νext , δ, and so on. The exposition is simplified
Figure 2. Regions in the AB plane of (theoretically expected) different dynamical behavior, at low, random connectivity and low rates. SS: stationary state stable; OS: oscillations. The curve has parts for B > 0, positive drift as well B < 0, negative drift. To obtain a continuous curve, between positive and negative B, we used B −1 as abscissa. Circles: points at which simulations have confirmed stability of stationary state; squares: simulations confirmed oscillations. Points marked “Sim” correspond to simulations discussed below; the point Pc on the transition line is the critical point (Ac = 0.468; Bc = 51.9) used in the expansion for the oscillation amplitude. Parameters of Network 1, Table 1.
by the fact that the dependence on all parameters has been reduced to two combinations, A and B. Once the different regions are delineated in the plane of A and B, one can obtain the classification by other parameters via the dependence of A and B on the desired parameter. In Fig. 2, corresponding to Network 1 in Table 1, we present in the AB plane a curve that separates the region in which the stationary state is stable (SS) from where oscillations appear (OS). To obtain a continuous curve that covers both positive and negative drift conditions, we have used the parameter B −1 . The bifurcation line is drawn fixing B and varying A until the eigenvalue with largest real part becomes purely imaginary; Table 1. Parameters for two networks used in the analysis and the simulations. In square brackets in the first column are the units in which the parameter is given. Network 1
Network 2
β [θ/ms]
.006
J [θ ]
.027
.03
[θ ]
.25
.5
C δ [ms]
.008
500
500
2
2
254
Mongillo and Amit
Figure 3. Regions in the AB plane of (theoretically expected) different dynamical behavior, at low, random connectivity and low rates. SS: stationary state stable; OS: oscillations. Circles: points at which simulations have confirmed SS state; squares: simulations confirmed oscillations. Parameters of Network 2, Table 1. Note that this network crosses from SS to OS dynamics operating always at negative drift. Notations as in Fig. 2: (Ac = 0.473; Bc = −96.8).
in this way we obtain the points (+) on the bifurcation line in Fig. 2. The various points in the plane indicate where simulations have been carried out to confirm the predicted dynamical state. In Fig. 3 we show the phase diagram of a network corresponding to Network 2 in Table 1. For these parameters the network operates always in the negative drift regime. If all parameters but one are fixed, one traces a curve in the AB plane, crossing from a regime of stable
stationary dynamics to an oscillatory regime. The transition occurs when the real part of the largest eigenvalue changes sign. We exhibit this in two situations: first varying νext and then varying δ. Figure 4A shows the variation of the real part of the maximal eigenvalue vs νext , at a fixed synaptic delay of 2 ms. The point at which the curve intersects the horizontal axis corresponds to crossing the line separating the SS from the OS regimes, in a figure like Fig. 2, along a specific trajectory. The bifurcation diagrams provide us with extra information: the real part of the leading eigenvalue is an estimate of the response time (|Re(λ)|−1 ) of the network to perturbations in a given dynamical state. In the SS states the response time represents the time scale of approach to the fixed point; in the OS states the time scale of the approach to the limit cycle. Finally, in Fig. 4B we draw a bifurcation curve as function of the synaptic delay δ emphasizing the crucial role of this parameter in the stability of the network’s dynamics. The network is Network 1 of Table 1 again, but this time we keep fixed the external rate νext (= 5 Hz) and vary δ. We find that the stationary state of this network becomes unstable for synaptic delays longer than ≈ 1.7 ms and remains stable for shorter delays.
3.2.
Instability: The Appearance of Oscillations
Near the bifurcation line, into the region of instability, the deviation, n 1 (t), of the rate begins to increase in
Figure 4. Bifurcation diagram—variation (across zero) of the real part of the bifurcating eigenvalue. Left: varying νext (δ = 2 ms); right: (c) = 4.32 Hz and δ (c) = 1.667 ms, respectively. Parameters varying δ (νext = 5 Hz). λ on the ordinate is in units of ms−1 . Critical parameters: νext as in Fig. 2.
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
magnitude and to oscillate. It assumes the form n 1 (t) = eλt/τ nˆ 1 (λ), where on the critical line λ = iωc is the bifurcating eigenvalue. The amplitude nˆ is controlled by an equation of the form (Brunel and Hakim, 1999; Bender and Orszag, 1987) τc
d nˆ 1 = W1 nˆ 1 − W3 |nˆ 1 |2 nˆ 1 . dt
(21)
W1 and W3 are complex numbers calculated in an expansion of the perturbation to the stationary state in powers of the deviation from the bifurcation line. Given that one is dealing with a perturbation in presence of resonant terms, secular terms are produced, which require an analysis at different time scales (see, e.g., Brunel and Hakim, 1999; Bender and Orszag, 1987). Details are reported in Appendix B. W1 , as expected, turns out to be a linear function of the deviation of A and B from their critical values. Its real part changes sign across the critical line. The particular critical points used below are marked as Pc on the separation line in Figs. 2 and 3. The results, though, are invariant to the particular choice of Pc . If one moves into a region of stability, Re(W1 ) < 0 and nˆ → 0 with t. Moving in the direction of instability, Re(W1 ) > 0 and the deviation begins to grow. Re(W3 ) > 0; hence the growth is saturated by the thirdorder term. The stable limit cycle of this supercritical Hopf bifurcation is nˆ 1 (t) = Rexp(iωt/τc ), where R =
Re(W1 ) , Re(W3 )
ω = I m(W1 ) − I m(W3 )
Re(W1 ) . Re(W3 )
Since n 1 (t) = nˆ 1 (t) exp(iωc t/τc ) + c.c., one has ν(t) = ν0 [1 + 2R cos(ωc + ω)t/τc ]. 3.3.
(22)
Rate Oscillations: Theory Versus Simulations
To test the theory, we have carried out simulations of networks in the various regions of expected dynamic behavior. The networks were composed of 104
255
inhibitory neurons with the dynamics of Eq. (1). The parameters of networks are reported in Table 1. We use a simulation program written by Massimo Mascaro (1998). The synaptic efficacies are drawn randomly and independently from a Gaussian distribution, with mean J and variance 2 J 2 ; the connection probability is C/N . The synaptic delays are uniform and fixed; the main effects of variability in the synaptic delays are a more stable stationary state (that is, an increased dimension of the SS region in the AB plane) and a slight reduction of the oscillation frequency (BH). The simulation is initiated by preparing the network in a state of global activity corresponding to mean-field theory prediction. In other words, the MF Eq. (5) is solved for the stationary rate, ν0 , of a network with the chosen parameters. For each neuron a spike is emitted with probability ν0 dt (dt is the integration time step) and Vi (0) = 0; otherwise Vi (0) is chosen in [0, θ ] with uniform probability. Each neuron is receiving an external stochastic afferent current, drawn from a Poisson distribution of mean νext dt. The network reaches its dynamical attractor in a very short time (≈10 ms). After the network has reached an attractor—a fixed point in the SS states or a stable limit cycle in the OS states—it is possible to estimate various observables: autocorrelation of global activity (AC), mean emission rate in both types of states, frequency and amplitude of oscillation in OS states, and so on. SS states would have a flat AC, while the OS states would have a sinusoidal AC. Alternatively, one expects a narrow distribution of the temporal intervals (T ) between successive maxima of the global activity in the OS states and a wide one for the SS states. The estimate of the amplitude of the oscillations is obtained averaging the maxima of the instantaneous rate (global spike activity); the period is estimated averaging over temporal intervals between successive maxima of the instantaneous rate. We moved from region to region by varying either the rate of the external current νext or the magnitude of the synaptic delays δ. The fixed, common parameters characterizing the networks we simulated are summarized in Table 1. The only exception is one case in which δ is varied. A first example is shown in Fig. 5A. We plot the total number of spikes emitted in the network (total activity) in bins of 0.5 ms, normalized to give the instantaneous rate (see, e.g., Amit and Brunel, 1997a). The network operates in a regime of positive drift, corresponding to the point (A = .456; B = 52.5, marked “Sim” in region SS) in Fig. 2.
256
Mongillo and Amit
Figure 5. Simulation vs theory: Left: spike activity (instantaneous rate in 0.5 ms bins), simulation in regime of stability. Parameters of Network 1, (c) Table 1, positive drift; νext = 3 Hz < νext . Horizontal line is mean-field prediction: νth = 2.51 Hz; observed ν = 2.53 ± 0.98 Hz. Right: AC function and normalized histogram of T . In SS states AC ∼ 1 and distribution of T is wide: T ∗ = 2.9 ± 1.4 ms (in contrast to OS states, see also Fig. 6).
The external rate is relatively low (νext = 3 Hz), which places the network in an SS regime, as theoretically expected. The agreement between the observed rate (2.53 ± 0.98 Hz) and the rate predicted by stationary mean field theory (2.51 Hz) is excellent. On the right of Fig. 5 we present the measured (simulation) AC and the normalized histogram of time intervals (T ) between successive maxima of the instantaneous rate. In an SS state, one expects the AC to be flat and the distribution of T to be wide and flat. In fact, both expectations are confirmed in this state as is underlined by the contrast with the AC and T distribution in Fig. 6 for the oscillating state of the same network. The ratio
Var(T )/T¯ 2 = 0.25 (compared to 0.0042 for the oscillating state). In Fig. 6 the simulation is carried out with an external rate above critical (νext = 6.5 Hz (A = 0.477; B = 46.5)). Spike emission becomes oscillatory, as predicted by theory and supported by the AC and T distribution on the right of the figure. The mean emission rate agrees again with the stationary mean-field prediction (solution of Eq. (5)), which corresponds also to Eq. (22) averaged over many oscillation periods. Theory: 5.93 Hz; simulation: 6.3 ± 0.36 Hz.2 The amplitude and the period of the oscillation oscillation are estimated to be νmax = 14.78 ± 1.25 Hz and
Figure 6. Left: Simulation vs theory in positive drift. Emission rate sampling is as in Fig. 5. Horizontal line represents the theoretical prediction for the oscillation amplitude. Oscillation period, theory: 5.81 ms; simulation: 6.0 ± .39; rate amplitude, theory: 15.02 Hz simulation: 14.78 ± 1.25 Hz; νext = 6.5 Hz. Right: AC and distribution of T in OS state. AC is sinusoidal and the distribution of T is peaked at T ∗ = 6 ms. (see also Fig. 5).
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
Figure 7. Simulation for Network 2 (Table 1), which operates only in negative drift. Top: SS state (νext = 3.5 Hz), horizontal line represents the theoretical prediction for the mean emission rate νth = 2.99 Hz νsim = 3.05 ± 0.21 Hz. Bottom: OS state (νext = 6 Hz), horizontal line: theoretical prediction for the rate oscillation amplitude. Oscillation period, theory: 6.03 ms; simulation: 6.2 ± 0.44 ms; rate oscillation amplitude, theory: 11.17 Hz simulation: 11.66 ± 1.4 Hz.
Tsim = 6.0 ± 0.39 ms, respectively. The theoretical reth sults are νmax = 15.03 Hz and Tth = 5.87. Network 2 operates only in the negative drift regime. Figure 7 exhibits the existence of both SS and OS dynamical regimes for this network. Figure 7A presents the instantaneous rate of emission vs t in the SS state (νext = 3.5 Hz; A = .462; B = −90.3) (Fig. 3). The average emission rate is measured to be νsim = 3.05 ± 0.21 Hz and predicted νth = 2.9 Hz, In this state, Var(T )/T¯ 2 = 0.88. Figure 7A presents the instantaneous rate in the OS state (νext = 6.0 Hz; (A = .479; B = −120.3), (in Fig. 3 the point marked “Sim” in the region OS). The average rate in this state is νsim = 5.87 ± 0.23 Hz, while theory gives νth = 5.51 Hz. For this state Var(T )/T¯ 2 = 0.0048, which should be contrasted with the corresponding value in the SS state. For the oscillation amplitude we find theory: 11.17 Hz, simulation: 11.66 ± 1.4 Hz. 4.
257
have certain attractive features, such as implementability in VLSI and some potential for modeling the effects of slow NMDA currents. In the latter we refer not to the saturation effects of NMDA, but rather to the fact that NMDA introduces a long time constant into the dynamics, which adds stability and reduces fluctuations and a weaker dependence on the finite size of the network (Brunel and Wang, 2001). Networks of LIF neurons have a similar characteristic. Several items are on the immediate agenda for networks of this type. Given that the neurons have apparently long time constants, what would be the response time of such a network to different perturbations from a stationary state? Simulations lead us to believe that such times will be relatively short and the technique can be applied to compute these times explicitly. In addition, the dynamical characteristics of working memory states, along the lines of Brunel (2000a), should be clarified and compared with simulations. Finally the behavior of such networks should be compared more carefully with that of exponential neurons networks to clarify if this model captures some aspects of neuromodulation. Appendix A A.1. Linear Stability The eigenmodes of Eq. (17) would be Q 1 (V¯ , t) = eλt/τ nˆ 1 (λ) Qˆ 1 (V¯ ), n 1 (t) = nˆ 1 (λ)eλt/τ , where Qˆ 1 (V¯ ) obeys the ordinary differential equation (λ − L) Qˆ 1 (V¯ ) = e−λδ/τ L AB [Q 0 ](V¯ )
(23)
together with the boundary conditions ∂ Qˆ 1 ¯ = 0, ¯ = −1 + Ae−λδ/τ , Qˆ 1 (θ) (θ) ¯ ∂ V 1 ∂ Qˆ 1 A ∂ Q 0 −λδ/τ − Qˆ 1 + e + BQ0 2 ∂ V¯ 2 ∂ V¯
(24)
V¯ =0
1 =− . 2 (25)
Discussion
The solution of the homogeneous part of Eq. (23) is By studying the various dynamical regimes of a purely inhibitory network of LIF neurons, we have set the ground for a more detailed investigation of networks of neurons of this type in more generic and structured situations, as was done for networks of exponential IF neurons (Brunel, 2000a). Networks of LIF neurons
) Qˆ (H = C1 φ1 (V¯ ) + C2 φ2 (V¯ ), 1
with ¯
φi = e zi V
and
zi = 1 ∓
√
1 + 2λ.
(26)
258
Mongillo and Amit
A particular solution of Eq. (23) is e−λδ/τ Qˆ 1 p = L AB [Q 0 ]. λ
At the onset of an oscillatory instability, the eigenvalue of the largest real part is purely imaginary, λ = iωc , where ωc is the frequency of the emerging oscillation:
The boundary conditions in (24) give a linear system for C1 and C2
e−iωc δ/τc ∓(1 − i)B √ = −1 + Ae−iωc δ/τc . ωc
e−λδ/τ L AB [Q 0 ](θ¯ ) = α(λ) λ (27) −λδ/τ ∂ e C1 φ1 (θ¯ ) + C2 φ2 (θ¯ ) = −1 + Ae−λδ/τ − λ ∂ V¯ × L AB [Q 0 ](θ¯ ) ≡ β(λ), (28)
C1 φ1 (θ¯ ) + C2 φ2 (θ¯ ) = −
αz 2 − β −1 αz 1 − β −1 ¯ φ1 (θ¯ ), C2 = − φ (θ). z2 − z1 z2 − z1 2 (29)
The boundary condition (25), ensuring equal flux at the threshold and at the reflecting lower barrier, gives Eq. (18) for the eigenfrequencies A.2. Limit δ/τ → 0 In this limit, the eigenvalues |λ| → ∞. That this is the case can be qualitatively understood as follows: the oscillations are due to the fact that an increase of activity in the network provokes an increase in the feedback inhibitory input after about one synaptic (delay) time. The effect of increased recurrent inhibitory input is a decrease in global activity. This argument predicts √ ω ∼ τ/δ. When |λ| → ∞, z 1,2 ∓ 2λ (Eq. 26). In Eq. (18) one of the exponentials dominates, and it reduces to αz = β,
(30)
with α(λ) and β(λ) given in Eq. (19). z is positive (negative) for positive(negative) drift. We then note √ that for Eq. (30) 2to have such a root, we need B ∼ |λ|. Since σl A = σ 2 +σ 2 < 1, one can neglect the term proportional ext l to A in the left-hand side and the term 2α in the righthand side. The equation becomes B
√ Bc = ∓ ωc sin(ωc δ/τc ) , Ac = sin(ωc δ/τc ) + cos(ωc δ/τc )
for negative and positive drift regimes, respectively. Appendix B
which leads to C1 =
Separating real and imaginary parts, one has
e−λδ/τ z = −1 + Ae−λδ/τ . λ
To determine the amplitude of the collective oscillation near (Ac , Bc ), we expand in powers of the deviation from the bifurcation line. The following discussion has a large overlap with that of Brunel and Hakim (1999). It is reproduced here for the sake of completeness in presence of significant technical differences. Hopf’s theorem ensures that the development parameter is of order of the square root of this displacement; therefore, A(= A − Ac ), B(= B − Bc ) enter only at the third-order. First we determine the second-order term. At the second order Eq. (11) becomes τc
∂ Q2 = L[Q 2 ] + n 2 (t − δ)L AB [Q 0 ] ∂t + n 1 (t − δ)L AB [Q 1 ],
(31)
with the boundary conditions ¯ t) = 0, Q 2 (θ, ∂ Q2 ¯ t) = − n 2 (t) + Ac n 2 (t − δ) (θ, ∂ V¯ − A2c n 21 (t − δ) + Ac n 1 (t)n 1 (t − δ), 1 ∂ Q2 − Q 2 + n 2 (t − δ)K AB [Q 0 ] 2 ∂ V¯ + n 1 (t − δ)K AB [Q 1 ] V¯ =0
1 = − n 2 (t), 2
(32)
where K AB = A2 ∂∂V¯ + B. The forcing terms on the righthand side of Eq. (31) contains terms at frequencies 2
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons ωc and 0; therefore,
together with the boundary conditions
Q 2 (V¯ , t) = |nˆ 1 |2 Qˆ 20 + nˆ 21 e2iωc t/τc Qˆ 22 + c.c.
(33)
n 2 (t) = ρ20 |nˆ 1 |2 + ρ22 e2iωc t/τc nˆ 21 + c.c.
(34)
Substituting into (31) one obtain an ordinary differential equation for Qˆ 22 — (2iωc − L) Qˆ 22 = ρ22 e−2iωc δ/τc L AB [Q 0 ] + e−iωc δ/τc L AB [ Qˆ 1 ]
(35)
—with the boundary conditions
A2c e−2iωc δ/τc
+ Ac e
¯ = 0, (41) Qˆ 20 (θ) ˆ ∂ Q 20 ¯ = ρ20 (−1 + Ac ) + 2Ac cos(ωc δ/τc ) − A2c (θ) ∂ V¯ and 1 ∂ Qˆ 20 ∗ − Qˆ 20 + ρ20 K AB [Q 0 ] + e−iωc δ/τc K AB 2 ∂ V¯ (hom) 1 (hom) iω δ/τ + e c c K AB Qˆ 1 × Qˆ 1 = − ρ20 . 2 ¯ V =0
Qˆ 22 (θ¯ ) = 0 ∂ Qˆ 22 (θ¯ ) = ρ22 −1 + Ac e−2iωc δ/τc ∂ V¯ −
259
(42) We look for a solution of the form
−iωc δ/τc
ˆ (lo) Qˆ 20 = Qˆ (hom) + ρ20 Qˆ (so) 20 20 + Q 20 ,
1 ∂ Qˆ 22 − Qˆ 22 + ρ22 e−2iωc δ/τc × 2 ∂ V¯ × K AB [Q 0 ] + e−iωc δ/τc K AB [ Qˆ 1 ]
where ¯ Qˆ (hom) = D1 + D2 e 2 V 20
V¯ =0
1 = − ρ22 . 2
¯ 2(V¯ −θ¯ ) Qˆ (so) 20 = (Ac + Bc ) V e e−iωc δ/τc ∗ ˆ (hom) Qˆ (lo) L AB Q 1 20 = iωc eiωc δ/τc − L AB Qˆ (hom) . 1 iωc
(36)
The general solution of Eq. (35) is Qˆ 22 = Qˆ (hom) + 22 (lo) ˆ + Q , where ρ22 Qˆ (so) 22 22 Qˆ (hom) (V¯ ) = γ1 ψ1 (V¯ ) + γ2 ψ2 (V¯ ) e−2iωc δ/τc Qˆ (so) = L AB [Q 0 ] 22 2iωc
(37) (38)
(43) (44)
(45)
As above, one determines the three unknowns D1 , D2 and ρ20 using Eqs. (41) and (42). Third Order
and ¯
ψi (V¯ ) = eξi V ; ξ1,2 = 1 ∓
1 + 4iωc .
Since [L , L AB ] = 0 and L[Q 0 ] = 0, (hom) e−2iωc δ/τc 2 e−iωc δ/τc Qˆ (lo) = L L [Q 0 ]. Qˆ 1 + AB 22 iωc 2(iωc )2 AB The boundary conditions in (36) determine γ1 and γ2 , in terms of ρ22 and the previously determined functions, as well as a linear algebraic equation for ρ22 . The component Qˆ 20 obeys −L[ Qˆ 20 ] = ρ20 L AB [Q 0 ] + e−iωc δ/τc L ∗AB [ Qˆ 1 ] (39) + eiωc δ/τc L AB [ Qˆ 1 ] (hom) −iωc δ/τc ∗ = ρ20 L AB [Q 0 ] + e L AB Qˆ 1 + eiωc δ/τc L AB Qˆ (hom) , (40) 1
At the third order, the naive expansion fails due to the appearance of a resonant term at frequency ωc , leading to secular terms that are an artifact of the truncation procedure (see, e.g., Bender and Orszag, 1987). The third-order terms obey τc
∂ Q3 = L[Q 3 ] + n 3 (t − δ)L AB [Q 0 ] ∂t + n 2 (t − δ)L AB [Q 1 ] + n 1 (t − δ)L AB [Q 2 ] 1 ∂ 2 Q0 + n 1 (t − δ) (A − Ac ) 2 ∂ V¯ 2 ∂ Q0 d nˆ 1 iωc τ + (B − Bc ) Q1 − − ¯ τc dt ∂V iωc t/τc iωc (t−δ)/τc ˆ + δe L AB [Q 0 ] . × τc Q 1 e (46)
260
Mongillo and Amit
The boundary conditions are Q 3 (θ¯ , t) = 0 ∂ Q3 (θ¯ , t) = −n 3 (t) + Ac n 3 (t − δ) ∂ V¯
(47)
together with boundary conditions:
− 2A2c n 1 (t − δ)n 2 (t − δ)
+ Ac (n 1 (t)n 2 (t − δ) + n 1 (t − δ)n 2 (t)) + A3c n 31 (t − δ) − A2c n 1 (t)n 21 (t − δ)
+ (A − Ac )n 1 (t − δ) d nˆ 1 iωc (t − δ)/τc − Ac δ e dt
(48)
K [Q 3 ] + n 3 (t − δ)K AB [Q 0 ] + n 2 (t − δ)K AB [Q 1 ] + n 1 (t − δ)K AB [Q 2 ] + n 1 (t − δ) 1 ∂ Q0 × (A − Ac) + (B − Bc )Q 0 2 ∂ V¯ d nˆ 1 iωc (t−δ)/τc −δ e K AB [Q 0 ] dt V¯ =0 1 = − n 3 (t). (49) 2 The terms in ddtnˆ 1 are due to the slow time scale, T ≡ t, needed to cancel secular terms. The first term arises from explicit time differentiation; the second comes from a first order expansion in of n(t − δ) = n(t − δ, T − δ). This term appears also in the boundary condition. Variation of A and B from Ac , Bc imply an additional variation of τ ; formally, τ (A, B) = τc (Ac , Bc ) + τ (A, B), where τ is linear in A and B. Indeed, τ depends on µ0 and σ02 , which, in turn, can be expressed in terms of A and B. The forcing terms on the right-hand side of Eq. (46) oscillate at frequencies 3ωc and ωc . Therefore, Q 3 (V¯ , t) = e3iωc t/τc Qˆ 33 + eiωc t/τc Qˆ 31 + c.c. n 3 (t) = nˆ 33 e3iωc t/τc + nˆ 31 eiωc t/τc + c.c. The terms at frequency ωc are resonant with the firstorder terms: (iωc − L) Qˆ 31 = nˆ 31 e−iωc δ/τc L AB [Q 0 ] + nˆ 1 |nˆ 1 |2 × ρ22 e−2iωc δ/τc L ∗AB [ Qˆ 1 ] + ρ20 L AB [ Qˆ 1 ] + e−iωc δ/τc L AB [ Qˆ 20 ] + eiωc δ/τc L AB [ Qˆ 22 ]
iωc τ ˆ Q1 + nˆ 1 e−iωc δ/τc L A−Ac ,B−Bc [Q 0 ] − τc d nˆ 1 ˆ − τc (50) ( Q 1 + δ/τc e−iωc δ/τc L AB [Q 0 ]), dt
¯ =0 Qˆ 31 (θ) ∂ Qˆ 31 ¯ = nˆ 31 (−1 + Ac e−iωc δ/τc ) + nˆ 1 1 (θ) ∂ V¯ d nˆ 1 + nˆ 1 |nˆ 1 |2 3 + (51) d dt 1 ∂ Qˆ 31 − Qˆ 31 + nˆ 31 e−iωc δ/τc K AB [Q 0 ] 2 ∂ V¯ V¯ =0 1 d nˆ 1 = − nˆ 31 + F nˆ 1 + G nˆ 1 |nˆ 1 |2 + H , 2 dt
(52)
in which we have denoted by 1 , 3 , d the coefficients of the linear, third-order, and derivative terms, respectively. The general solution of Eq. (50) is Qˆ 31 (V¯ ) ¯ = K¯ 1 φ1 (V¯ ) + K¯ 2 φ2 (V¯ ) + nˆ 31 Qˆ 1 p (V¯ ) + Qˆ (lo) 31 ( V ), where ψ1 , ψ2 , and Qˆ 1 p are the functions that appear at first order, and d nˆ 1 ¯ Qˆ (lo) Q d (V¯ ) + nˆ 1 Q l (V¯ ) + nˆ 1 |nˆ 1 |2 Q c (V¯ ). 31 ( V ) = dt Using (51), one obtains ¯ + K¯ 2 φ2 (θ) ¯ = −nˆ 31 Qˆ 1 p (θ) ¯ − Qˆ (lo) ¯ K¯ 1 φ1 (θ) 31 (θ) (53) ¯ + K¯ 2 φ2 (θ) ¯ K¯ 1 φ1 (θ) ∂ ˆ ¯ = nˆ 31 −1 + Ac e−iωc δ/τc − Q 1 p (θ) ∂ V¯ ∂ ˆ (lo) ¯ + , − (54) Q 31 (θ) ∂ V¯ where stands for the sum of the last three terms in Eq. (52). The coefficients of nˆ 31 in the right-hand side of (53) are respectively α(iωc ) and β(iωc ), defined in Eq. (19). Therefore, K¯ 1 = nˆ 31 C1 + K 1 with C1 is defined in Eq. (29) and K 1 is the remainder ˆ 31 Qˆ 1 + Qˆ 3 . Qˆ 31 = nˆ 31 Qˆ 1 + K 1 φ1 + K 2 φ2 + Qˆ (lo) 31 = n (55)
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
Using (52), one finally obtains 1 ∂ Qˆ 1 −iωc δ/τc ˆ nˆ 31 K AB [Q 0 ] − Q1 + e 2 ∂ V¯ 1 ∂ Qˆ 3 1 + − Qˆ 3 = − nˆ 31 + F nˆ 1 + G nˆ 1 |nˆ 1 |2 2 ∂ V¯ 2 V¯ =0
d nˆ 1 . dt Noting that the term multiplying nˆ 31 is the eigenvalue Eq. (18), the equation reduces to 1 ∂ Qˆ 3 ˆ d nˆ 1 − Q3 . (56) = F nˆ 1 +G nˆ 1 |nˆ 1 |2 + H 2 ∂ V¯ dt ¯ V =0 +H
When Qˆ 3 is expressed, from Eq. (55) in terms of the first order, third order, and the derivative of nˆ 1 , one arrives at the (Landau) equation of motion for nˆ 1 : d nˆ 1 = W1 nˆ 1 − W3 nˆ 1 |nˆ 1 |2 (57) dt In the limit δ¯ → 0 the expressions for W1 and W3 simplify: Eq. (56) becomes τc
∂ ˆ (lo) ¯ Qˆ (lo) Q 31 (θ¯ ) + = 0, 31 (θ )z − ∂ V¯ √ where z = ± 2iωc . For W1 and W3 , one obtains Q l (θ¯ )z − W1 = − Q d (θ¯ )z − Q c (θ¯ )z − W3 = Q d (θ¯ )z −
∂ Ql ¯ (θ ) + 1 ∂ V¯ ∂ Qd ¯ (θ ) + d ∂ V¯ ∂ Qc ¯ (θ ) + 3 ∂ V¯ . ∂ Qd ¯ (θ ) + d ∂ V¯
, (58)
Acknowledgments We are indebted to Dr. Nicolas Brunel for very helpful discussions and to Dr. Massimo Mascaro for providing us with his simulation program. This study was supported by a Center of Excellence grant “Statistical Mechanics and Complexity” of the INFM, Roma, Italia and Center of Excellence grant “changing your mind” of the Israel Science Foundation. Notes 1. We follow the convention of taking the afferent current for the RC neuron in units of potential. 2. To eliminate the effect of the oscillations on the estimate of the variance in the oscillating states, we have averaged the rate in entire cycles and calculated the variance of the rate from cycle to cycle.
261
References Amit DJ, Brunel N (1997) Dynamics of a recurrent network of spiking neurons before and following learning. Network Comput. Neural Sys. 8: 373–404. Amit DJ, Brunel N (1997) Model of global spontaneous activity and local structured (learned) delay activity during delay periods in cerebral cortex. Cerebral Cortex 7: 237. Bender CM, Orszag SA (1987) Advanced Mathematical Methods for Scientists and Engineers. McGraw-Hill, Singapore. Brunel N (2000) Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neuroscience 8: 183. Brunel N (2000) Persistent activity and the single cell frequencycurrent curve in a cortical network model. Network 11: 261. Brunel N, Hakim V (1999) Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11: 1621. Brunel N, Wang XJ (2001) A cortical network model of object working memory: Resistance to distractors and neuromodulation. J. Comput. Neurosci. 11: 63. Buzsaki G, Hovrath Z, Urioste R, Hetke J, Wise K (1992) Highfrequency network oscillations in hippocampus. Science 256: 1025. Compte A, Brunel N, Goldman-Rakic PS, Wang XJ (2000) Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cerebral Cortex 10: 910. Fusi S, Mattia M (1999) Collective behavior of networks with linear (VLSI) integrate and fire neurons. Neural Comput. 11: 633. Fusi S, Annunziato M, Badoni D, Salamon A, Amit DJ (2000) Spikedriven synaptic plasticity: Theory, simulation, VLSI implementation. Neural Comput. 12: 2227. Gerstein GL, Mandelbrot B (1964). Random walk models for the spike activity of a single neuron. Biophysical J. 4: 41. Gerstner W (1995) Time structure of the activity in neural network models. Phys. Rev. E 51: 738. Kahn PB, Zarmi Y (1998) Nonlinear Dynamics: Exploration Through Normal Forms. Wiley, New York. Knight BW (1972) Dynamics of encoding in a population of neurons. J. Gen. Physiol. 59: 734. Rauch A, La Camera G, Fusi S, Senn W, Luscher HR (2001) Discharge rate of neocortical pyramidal cells in response to noisy input currents. Submitted for oral presentation to CNS2001. Mascaro M (1998) Tesi di laurea, University of Roma La Sapienza. (in Italian). Funzione Di Risposta Efficace Di Una Sotto Popolazione Di Neuroni In Una Rete Multimodulare. Mascaro M, Amit DJ (1999) Effective neural response function for collective population states. Network 10: 351. Mattia M (1997), Tesi di laurea, University Roma La Sapienza. (in Italian). Dinamica Di Una Rete Di Neuroni Impucsivi Con Depolarizzazione Lineave. Mead C (1989) Analog VLSI and Neural System. Addison-Wesley, Reading, MA. Traub RD, Whittington MA, Collins SB, Buzsaki G, Jefferys JE (1996) Analysis of gamma rhythms in the rat hippocampus in vitro and in vivo. J. Physiol. 493: 471. Van Vreeswijk C, Abbott L, Ermentrout GB (1994) When inhibitio not excitation synchronizes firing. J. Comput. Neurosci. 1: 313.
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons GIANLUIGI MONGILLO AND DANIEL J. AMIT Department of Physics, University of Rome La Sapienza, Ple. Aldo Moro, 00185, Rome, Italy [email protected]
Received February 12, 2001; Revised July 25, 2001; Accepted October 12, 2001 Action Editor: Roger Traub
Abstract. The dynamics of a network of randomly connected inhibitory linear integrate and fire (LIF) neurons (with a floor for the depolarization), in the presence of stochastic external afferent input, is considered in various parameter regimes of the neurons and of the network. Applying a technique recently introduced by Brunel and Hakim, we classify the regimes in which such a network has stable stationary states and in which spike emission rates oscillate. In the vicinity of the bifurcation line, the oscillation frequency and its amplitude are computed and compared with simulations. As for leaky IF neurons, the space of parameters can be compactified into two. Yet despite significant technical differences between the two models, related to both the different dynamics of the depolarization as well as to the different boundary conditions, the qualitative behavior is rather similar. The significance of LIF neurons and of the differences with leaky IF neurons is discussed. Keywords: spiking neuron, recurrent networks, collective oscillations.
1.
Introduction
Spiking neurons that integrate linearly (rather than exponentially) and have a reflective floor for their depolarization (LIF), present several attractive features: they are natural candidates for Very Large Scale Integration (VLSI) implementations (Mead, 1989; Fusi et al., 2001); the single-cell spike dynamics is relatively simple (Mattia, 1997; Fusi and Mattia, 1999); and due to their depolarization floor, networks of such neurons have both low (negative drift) and high spike rates (positive drift) (Mattia, 1997; Fusi and Mattia, 1999) and hence coexistence of spontaneous and selective activity states (Mattia, 1997). Recently (Rauch et al., 2001), it has been observed that the emission rate of a single cell to a gaussian injected current is rather well described, with few free parameters, by the transduction function of the LIF neuron, in a rather wide range of low rates. Both facts naturally raise the question as to whether networks of spiking LIF neurons would exhibit collec-
tive behavior that is qualitatively similar to that of the more common RC neurons. Historically, such neurons have been studied in the absence of a depolarization floor (Gerstein and Mandelbrot, 1964; Knight, 1972). In that case they emit spikes regularly only for positive drift; their characteristics are rather transparent and not very rich in phenomenology. In a recent study Brunel and Hakim (1999) have paved the way for a detailed investigation of the dynamic characteristics of rather complex networks of spiking neurons. Following on the empirical observation (Amit and Brunel, 1997a) that spike dynamics in a network is dominated mostly by the evolution of the global activity of the network, Brunel and Hakim have set up a perturbative study of the network’s spike dynamics around its mean-field stationary states. The approach constitutes a powerful tool that allows an effective classification, in network parameter space, of activity states of different types, including stationary
250
Mongillo and Amit
states of irregular spiking, oscillatory states for the global population with irregularly firing neurons, and synchronous states. The initial study considered networks constituted exclusively of inhibitory neurons driven by external currents. It was then extended to networks with both excitatory and inhibitory IF neurons, which can sustain spontaneous activity and selective delay activity following learning (Brunel, 2001a). Beyond the clarification of the dependence of different network operating states on critical network parameters, the technique developed allows also the investigation of the reaction of a network of spiking neurons, in a stationary state, to an external perturbation, such as the change in external afferents or the presentation of a stimulus. Here we apply this tool to a network of inhibitory spiking LIF neurons. The neurons in the network are connected randomly, and external currents are stochastic. This study opens the way for the analysis of the domains of different collective behavior of networks with LIF neurons, which is relevant to the correspondence between VLSI networks and networks of more biologically plausible RC neurons. Given the complexity of the system it is quite rewarding that the theoretical technique, while somewhat complicated, allows for detailed quantitative predictions of parameter regions for states of different dynamical characteristics (stationary, oscillatory), mean emission rates, and oscillation periods and amplitudes of the rates, where the network transits into oscillating behavior. All these are compared against simulations. 2. 2.1.
current tends to hyperpolarize the neuron. Emission in this case is due to fluctuations. We define the drift as µ ≡ −β + I , where I is the mean afferent current. When it is positive, firing is essentially regular, and the rate of emission is rather insensitive to fluctuations. The presence of the two firing regimes is a sufficient condition for networks of such neurons to have both spontaneous and working memory activity states (Brunel, 2001b). The RC neuron also has such two regimes of operation, depending on whether the mean afferent current is below or above threshold, respectively. We therefore compare (Fig. 1) the behavior of the two models by superposing the corresponding dynamical regimes. In other words, we map the negative drift region of the LIF neuron on the subthreshold regime of the RC neuron (ND regimes) and the complementary (SD regimes) on each other. To render the variables describing the two neurons of equal dimensions, we define µ∗ as µ for the LIF neuron and ( I − θ )/τ for the RC neuron.1 The correspondence of the noise 2 2 variable is σLIF to σ RC /τ , where τ is the membrane time constant of RC neuron. Note that if the reflecting barrier is absent, V can become arbitrarily negative,
Methods LIF Neuron and Stationary States of Networks
The evolution of the depolarization V (t) of the LIF neuron (Mead, 1989) is V˙ = −β + I (t),
(1)
where β is the linear loss, and I (t) is the afferent current charging the membrane. This equation is accompanied by the standard spike emission condition at an upper threshold θ, when V (t k ) = θ a spike is emitted and V is reset to → Vr . The VLSI-LIF neuron has a floor for its depolarization; V cannot be hyperpolarized below V = 0. When the depolarization is pushed to negative values, it remains at 0, a reflecting barrier. This allows spike emission also when the mean afferent
Figure 1. Comparison of current-to-rate transduction functions of LIF (full curves) and exponential (RC) neuron (dotted), for two values of the corresponding variances of the afferent current. The correspondence used: matching the subthreshold regime of the RC neuron (µRC < θ ) with the ND regime of the LIF neuron (see text). The plots exhibit the emission rates of the two neurons vs µ∗ , where µ∗ = (µRC − θ )/τ , for the RC neuron and is µ for the LIF neuron. The variances put in correspondence are σ 2 for the LIF neuron and σ 2 /τ for the RC neuron. Note the higher rates of the LIF neuron at corresponding negative values of the µ∗ and the consequent smaller gain of this neuron at µ∗ = 0. θ = 20 mV for both types of neurons, τRC = 20 ms.
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
and the LIF neuron operates exclusively in the positive drift regime (see, e.g., Fusi and Mattia, 1999). This is of little consequence for the exponential neuron. If the afferent current to the neuron is stationary, Gaussian, with mean µ, with variance σ 2 per unit time, and independent at different times, the spike emission rate is calculated to be (Fusi and Mattia, 1999) (µ, σ ) = 2
−1 σ 2 2µθ − 2µθ + e σ2 − 1 . 2µ2 σ 2
(2)
The significance of this result is that in a network of interacting spiking neurons with low, random connectivity and low rates, the spike trains afferent along different presynaptic channels often are independent to a good approximation. Moreover, the synaptic efficacies are rather small relative to the spike-emission threshold. Hence, the recurrent currents can be described as Gaussian, with a mean and width determined by the characteristics of functional subpopulations of neurons and by the average spike rates in these populations. In this way one arrives at a mean-field description of the average spike dynamics in the entire network (Amit and Brunel, 1997b; Mascaro and Amit, 1999). In particular, if the network is a homogeneous population of inhibitory neurons connected by synapses of connectivity C, of efficacy with average J and variance 2 J 2 , then the mean and variance of the recurrent afferent currents is µl (t) = −C J ν(t); σl2 (t) = C J 2 (1 + 2 )ν(t), (3) where ν(t) is the spike rate in the population. When external currents are added, with mean and variance 2 µext , σext , respectively, one has µ(t) = −β + µl (t) + µext (t); 2 σ 2 (t) = σl2 (t) + σext (t).
(4)
Given the hypotheses of mean-field, with mean and variance linear functions of the instantaneous rate, we can write the equation for the stationary states of the networks, those activity states that reproduce in output the rate appearing on input: ν = (µ(ν), σ 2 (ν)).
(5)
The solutions, ν0 , of this equation must be accompanied by a test of stability. A stability analysis, for the stationary states of a network of leaky IF neurons, has
251
been developed (Amit and Brunel, 1997a), based on dynamic equations for µ and σ in conditions of very slow variations for the underlying integrator. Both simulations (Van Vreeswijk et al., 1994; Gerstner, 1995) and experiment (Buzsaki et al., 1992; Traub et al., 1996) (for many additional references, see, e.g, Brunel and Hakim, 1999) indicate that networks of the type considered may oscillate at high frequency unrelated to the spike emission rate, even in ranges of parameters where the test predicts stability. Moreover, such an asymptotic regime, in the dynamics of the depolarization, does not exist for the LIF neuron. Hence the requirement for a finer tool of control over the dynamics of the network and of the stability of its stationary states. Following Brunel and Hakim, one finds that not only can one obtain a proper classification of the parameter regimes of stable and unstable (oscillatory) states, but the frequency and amplitude of the spike emission rates in the oscillatory states can be fully computed near the line of instability.
2.2.
Perturbation of Stationary States
To go beyond the study of rates in stationary states, one observes (Brunel and Hakim, 1999; Amit and Brunel, 1997b) that in randomly and sparsely connected networks, firing patterns of different neurons are essentially independent. The network dynamics can be therefore characterized by the instantaneous distribution of depolarizations in the neural population together with the spike emission probability. In networks in which there is more than one functionally distinct population one would have, of course, as many distribution functions and spike emission probabilities as distinct populations (see, e.g., Brunel, 2000a). The temporal evolution of the probability distribution of the depolarization (the Fokker-Plank equation) corresponding to a neuron whose dynamics is governed by Eq. (1) is σ 2 (t) ∂ 2 P ∂P ∂P = , − µ(t) 2 ∂t 2 ∂V ∂V
(6)
where µ and σ are given in Eqs. (3) and (4). Note that the rate ν entering the mean and variance at time t is the rate at t − δ, the rate a synaptic delay constant (δ) earlier. Like for the exponential IF neuron, the boundary condition at the spike emission threshold is P(θ, t) = 0. Moreover, the probability flux at θ gives
252
Mongillo and Amit
the emission probability at time t, ν(t), that is, 1 ∂P ν(t) = − σ 2 (t) (θ, t). 2 ∂V
where
At the reset potential (Vr = 0, for simplicity), the boundary condition is different from the one for the exponential IF neuron. For the linear neuron, there is a reflecting barrier (Fusi and Mattia, 1999.) No process can cross V = 0 from below. The only flux across this barrier is due to spike emission. Hence, 1 2 ∂P 1 ∂P σ (t) (0, t) − µ(t)P(0, t) = σ 2 (t) (θ, t). 2 ∂V 2 ∂V (8) In a stationary state, Eq. (6) for the stationary pdf, P0 (V ), becomes σ02 ∂ 2 P0 ∂ P0 − µ0 = 0, 2 ∂V 2 ∂V whose solution is (Mattia, 1997; Fusi and Mattia, 1999) P0 (V )
µ0 ν0 1 − exp −2 2 (θ − V ) = in V ∈ [0, θ ]. µ0 σ0 (9)
Imposing the normalization condition on P0 one obtains the transduction Eq. (2), for the stationary emission rate. Next we rescale the variables in Eq. (6) as follows: P=
σl2 CJ 2 (1 + 2 )ν0 = , 2 σ02 CJ (1 + 2 )ν0 + σext CJν0 µl B= =− . −β − CJν0 + µext µ0
(7)
2ν0 σ2 Q; V = 0 V¯ ; ν(t) = ν0 (1 + n(t)). µ0 µ0 (10)
A=
2
Moreover, τ ≡ σ02 /µ20 will play below the role of an intrinsic time constant, which is not present in the evolution of the depolarization: it determines the temporal scale on which there is a cross-over from small to large synaptic delays. The differences remain in the form of the operators L and in the boundary conditions. Note that A and B (in distinction to those of the IF neuron) are, respectively, the ratio of the variance and the mean of the local currents to the total variance and currents. A remains in the range [0, 1]. But B, positive definite for the IF neurons (= −µl /σl ), is no longer so in our case: µ0 will be indeed negative for negative drift. The boundary conditions for Q become ¯ t) = 0; Q(θ, (13) ∂Q 1 + n(t) ¯ t) = − (θ, 1 + An(t − δ) ∂ V¯ 1 ∂Q A ∂Q × − Q + n(t − δ) + BQ 2 ∂ V¯ 2 ∂ V¯ V¯ =0
1 = − [1 + n(t)], 2 where θ¯ = 2.3.
(14)
µ0 θ. σ02
Linear Stability Analysis
We write the deviations from the stationary state as The new variable are all dimensionless. This transformation leads Eq. (6) to an equation rather similar to the one of the exponential neuron, with only two parameters: τ
∂Q ∂t
1 ∂2 Q ∂Q = − + n(t − δ) 2 ¯ 2 ∂V ∂ V¯
Q(V¯ , t) = Q 0 (V¯ ) + Q 1 (V¯ , t) + O( 2 ) (15) n(t) = n 1 (t) + O( 2 ).
(16)
At first order, one has
∂Q A ∂2 Q +B 2 ¯ 2 ∂V ∂ V¯
(11)
≡ L[Q](V¯ , t) + n(t − δ)L AB [Q](V¯ , t),
(12)
τ
∂ Q1 = L[Q 1 ](V¯ , t) + n 1 (t − δ)L AB [Q 0 ](V¯ ), (17) ∂t
in which L and L AB are linear operators defined in Eq. (12). This leads (see, e.g., Appendix A) to an
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
253
eigenvalue equation of the form 1 −z 1 θ¯ z 1 (αz 2 − β)e −1 z2 − z1 2 1 ¯ z2 − (αz 1 − β)e−z2 θ − 1 + He−λδ/τ = − , 2 2 (18) where −λδ/τ α(λ) = (A + B) e λ β(λ) = −1 + Ae−λδ/τ + 2α (19) 1 ¯ ¯ H = B(1 − e−2θ ) − Ae−2θ 2 √ and z 1,2 = 1 ∓ 1 + 2λ with the positive (negative) sign corresponding to positive (negative) drift. 2.4.
Limit of Short Delays
For the linear neuron, for a wide range of parameters, σ02 is relatively large, ≈ 100 ms, while δ is a few ms. µ20 So δ/τ 1. In this limit √ the eigenvalues become very large, and z 1,2 = ∓ 2λ. In Eq. (18) one of the two exponentials in square brackets dominates, the one with z i θ¯ > 0, and the equation reduces to (Appendix A2) αz = β.
(20)
Note that θ¯ is positive (negative) for positive (negative) drift. Hence, dealing with the LIF neuron, we use Eq. (20), after having checked that the solutions for the eigenvalues agree very closely with those of the full Eq. (18). 3. 3.1.
Results Stability, Instability, Bifurcation
Using Eq. (20) one identifies the regions, in the space of parameters, in which the stationary state of the network is stable: all eigenvalues have a negative real part and is unstable: some eigenvalues have a positive real part. There are also bifurcation lines: those at which the eigenvalue with the largest real part is purely imaginary. To start with, the number of parameters is quite high: β, C, , νext , δ, and so on. The exposition is simplified
Figure 2. Regions in the AB plane of (theoretically expected) different dynamical behavior, at low, random connectivity and low rates. SS: stationary state stable; OS: oscillations. The curve has parts for B > 0, positive drift as well B < 0, negative drift. To obtain a continuous curve, between positive and negative B, we used B −1 as abscissa. Circles: points at which simulations have confirmed stability of stationary state; squares: simulations confirmed oscillations. Points marked “Sim” correspond to simulations discussed below; the point Pc on the transition line is the critical point (Ac = 0.468; Bc = 51.9) used in the expansion for the oscillation amplitude. Parameters of Network 1, Table 1.
by the fact that the dependence on all parameters has been reduced to two combinations, A and B. Once the different regions are delineated in the plane of A and B, one can obtain the classification by other parameters via the dependence of A and B on the desired parameter. In Fig. 2, corresponding to Network 1 in Table 1, we present in the AB plane a curve that separates the region in which the stationary state is stable (SS) from where oscillations appear (OS). To obtain a continuous curve that covers both positive and negative drift conditions, we have used the parameter B −1 . The bifurcation line is drawn fixing B and varying A until the eigenvalue with largest real part becomes purely imaginary; Table 1. Parameters for two networks used in the analysis and the simulations. In square brackets in the first column are the units in which the parameter is given. Network 1
Network 2
β [θ/ms]
.006
J [θ ]
.027
.03
[θ ]
.25
.5
C δ [ms]
.008
500
500
2
2
254
Mongillo and Amit
Figure 3. Regions in the AB plane of (theoretically expected) different dynamical behavior, at low, random connectivity and low rates. SS: stationary state stable; OS: oscillations. Circles: points at which simulations have confirmed SS state; squares: simulations confirmed oscillations. Parameters of Network 2, Table 1. Note that this network crosses from SS to OS dynamics operating always at negative drift. Notations as in Fig. 2: (Ac = 0.473; Bc = −96.8).
in this way we obtain the points (+) on the bifurcation line in Fig. 2. The various points in the plane indicate where simulations have been carried out to confirm the predicted dynamical state. In Fig. 3 we show the phase diagram of a network corresponding to Network 2 in Table 1. For these parameters the network operates always in the negative drift regime. If all parameters but one are fixed, one traces a curve in the AB plane, crossing from a regime of stable
stationary dynamics to an oscillatory regime. The transition occurs when the real part of the largest eigenvalue changes sign. We exhibit this in two situations: first varying νext and then varying δ. Figure 4A shows the variation of the real part of the maximal eigenvalue vs νext , at a fixed synaptic delay of 2 ms. The point at which the curve intersects the horizontal axis corresponds to crossing the line separating the SS from the OS regimes, in a figure like Fig. 2, along a specific trajectory. The bifurcation diagrams provide us with extra information: the real part of the leading eigenvalue is an estimate of the response time (|Re(λ)|−1 ) of the network to perturbations in a given dynamical state. In the SS states the response time represents the time scale of approach to the fixed point; in the OS states the time scale of the approach to the limit cycle. Finally, in Fig. 4B we draw a bifurcation curve as function of the synaptic delay δ emphasizing the crucial role of this parameter in the stability of the network’s dynamics. The network is Network 1 of Table 1 again, but this time we keep fixed the external rate νext (= 5 Hz) and vary δ. We find that the stationary state of this network becomes unstable for synaptic delays longer than ≈ 1.7 ms and remains stable for shorter delays.
3.2.
Instability: The Appearance of Oscillations
Near the bifurcation line, into the region of instability, the deviation, n 1 (t), of the rate begins to increase in
Figure 4. Bifurcation diagram—variation (across zero) of the real part of the bifurcating eigenvalue. Left: varying νext (δ = 2 ms); right: (c) = 4.32 Hz and δ (c) = 1.667 ms, respectively. Parameters varying δ (νext = 5 Hz). λ on the ordinate is in units of ms−1 . Critical parameters: νext as in Fig. 2.
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
magnitude and to oscillate. It assumes the form n 1 (t) = eλt/τ nˆ 1 (λ), where on the critical line λ = iωc is the bifurcating eigenvalue. The amplitude nˆ is controlled by an equation of the form (Brunel and Hakim, 1999; Bender and Orszag, 1987) τc
d nˆ 1 = W1 nˆ 1 − W3 |nˆ 1 |2 nˆ 1 . dt
(21)
W1 and W3 are complex numbers calculated in an expansion of the perturbation to the stationary state in powers of the deviation from the bifurcation line. Given that one is dealing with a perturbation in presence of resonant terms, secular terms are produced, which require an analysis at different time scales (see, e.g., Brunel and Hakim, 1999; Bender and Orszag, 1987). Details are reported in Appendix B. W1 , as expected, turns out to be a linear function of the deviation of A and B from their critical values. Its real part changes sign across the critical line. The particular critical points used below are marked as Pc on the separation line in Figs. 2 and 3. The results, though, are invariant to the particular choice of Pc . If one moves into a region of stability, Re(W1 ) < 0 and nˆ → 0 with t. Moving in the direction of instability, Re(W1 ) > 0 and the deviation begins to grow. Re(W3 ) > 0; hence the growth is saturated by the thirdorder term. The stable limit cycle of this supercritical Hopf bifurcation is nˆ 1 (t) = Rexp(iωt/τc ), where R =
Re(W1 ) , Re(W3 )
ω = I m(W1 ) − I m(W3 )
Re(W1 ) . Re(W3 )
Since n 1 (t) = nˆ 1 (t) exp(iωc t/τc ) + c.c., one has ν(t) = ν0 [1 + 2R cos(ωc + ω)t/τc ]. 3.3.
(22)
Rate Oscillations: Theory Versus Simulations
To test the theory, we have carried out simulations of networks in the various regions of expected dynamic behavior. The networks were composed of 104
255
inhibitory neurons with the dynamics of Eq. (1). The parameters of networks are reported in Table 1. We use a simulation program written by Massimo Mascaro (1998). The synaptic efficacies are drawn randomly and independently from a Gaussian distribution, with mean J and variance 2 J 2 ; the connection probability is C/N . The synaptic delays are uniform and fixed; the main effects of variability in the synaptic delays are a more stable stationary state (that is, an increased dimension of the SS region in the AB plane) and a slight reduction of the oscillation frequency (BH). The simulation is initiated by preparing the network in a state of global activity corresponding to mean-field theory prediction. In other words, the MF Eq. (5) is solved for the stationary rate, ν0 , of a network with the chosen parameters. For each neuron a spike is emitted with probability ν0 dt (dt is the integration time step) and Vi (0) = 0; otherwise Vi (0) is chosen in [0, θ ] with uniform probability. Each neuron is receiving an external stochastic afferent current, drawn from a Poisson distribution of mean νext dt. The network reaches its dynamical attractor in a very short time (≈10 ms). After the network has reached an attractor—a fixed point in the SS states or a stable limit cycle in the OS states—it is possible to estimate various observables: autocorrelation of global activity (AC), mean emission rate in both types of states, frequency and amplitude of oscillation in OS states, and so on. SS states would have a flat AC, while the OS states would have a sinusoidal AC. Alternatively, one expects a narrow distribution of the temporal intervals (T ) between successive maxima of the global activity in the OS states and a wide one for the SS states. The estimate of the amplitude of the oscillations is obtained averaging the maxima of the instantaneous rate (global spike activity); the period is estimated averaging over temporal intervals between successive maxima of the instantaneous rate. We moved from region to region by varying either the rate of the external current νext or the magnitude of the synaptic delays δ. The fixed, common parameters characterizing the networks we simulated are summarized in Table 1. The only exception is one case in which δ is varied. A first example is shown in Fig. 5A. We plot the total number of spikes emitted in the network (total activity) in bins of 0.5 ms, normalized to give the instantaneous rate (see, e.g., Amit and Brunel, 1997a). The network operates in a regime of positive drift, corresponding to the point (A = .456; B = 52.5, marked “Sim” in region SS) in Fig. 2.
256
Mongillo and Amit
Figure 5. Simulation vs theory: Left: spike activity (instantaneous rate in 0.5 ms bins), simulation in regime of stability. Parameters of Network 1, (c) Table 1, positive drift; νext = 3 Hz < νext . Horizontal line is mean-field prediction: νth = 2.51 Hz; observed ν = 2.53 ± 0.98 Hz. Right: AC function and normalized histogram of T . In SS states AC ∼ 1 and distribution of T is wide: T ∗ = 2.9 ± 1.4 ms (in contrast to OS states, see also Fig. 6).
The external rate is relatively low (νext = 3 Hz), which places the network in an SS regime, as theoretically expected. The agreement between the observed rate (2.53 ± 0.98 Hz) and the rate predicted by stationary mean field theory (2.51 Hz) is excellent. On the right of Fig. 5 we present the measured (simulation) AC and the normalized histogram of time intervals (T ) between successive maxima of the instantaneous rate. In an SS state, one expects the AC to be flat and the distribution of T to be wide and flat. In fact, both expectations are confirmed in this state as is underlined by the contrast with the AC and T distribution in Fig. 6 for the oscillating state of the same network. The ratio
Var(T )/T¯ 2 = 0.25 (compared to 0.0042 for the oscillating state). In Fig. 6 the simulation is carried out with an external rate above critical (νext = 6.5 Hz (A = 0.477; B = 46.5)). Spike emission becomes oscillatory, as predicted by theory and supported by the AC and T distribution on the right of the figure. The mean emission rate agrees again with the stationary mean-field prediction (solution of Eq. (5)), which corresponds also to Eq. (22) averaged over many oscillation periods. Theory: 5.93 Hz; simulation: 6.3 ± 0.36 Hz.2 The amplitude and the period of the oscillation oscillation are estimated to be νmax = 14.78 ± 1.25 Hz and
Figure 6. Left: Simulation vs theory in positive drift. Emission rate sampling is as in Fig. 5. Horizontal line represents the theoretical prediction for the oscillation amplitude. Oscillation period, theory: 5.81 ms; simulation: 6.0 ± .39; rate amplitude, theory: 15.02 Hz simulation: 14.78 ± 1.25 Hz; νext = 6.5 Hz. Right: AC and distribution of T in OS state. AC is sinusoidal and the distribution of T is peaked at T ∗ = 6 ms. (see also Fig. 5).
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
Figure 7. Simulation for Network 2 (Table 1), which operates only in negative drift. Top: SS state (νext = 3.5 Hz), horizontal line represents the theoretical prediction for the mean emission rate νth = 2.99 Hz νsim = 3.05 ± 0.21 Hz. Bottom: OS state (νext = 6 Hz), horizontal line: theoretical prediction for the rate oscillation amplitude. Oscillation period, theory: 6.03 ms; simulation: 6.2 ± 0.44 ms; rate oscillation amplitude, theory: 11.17 Hz simulation: 11.66 ± 1.4 Hz.
Tsim = 6.0 ± 0.39 ms, respectively. The theoretical reth sults are νmax = 15.03 Hz and Tth = 5.87. Network 2 operates only in the negative drift regime. Figure 7 exhibits the existence of both SS and OS dynamical regimes for this network. Figure 7A presents the instantaneous rate of emission vs t in the SS state (νext = 3.5 Hz; A = .462; B = −90.3) (Fig. 3). The average emission rate is measured to be νsim = 3.05 ± 0.21 Hz and predicted νth = 2.9 Hz, In this state, Var(T )/T¯ 2 = 0.88. Figure 7A presents the instantaneous rate in the OS state (νext = 6.0 Hz; (A = .479; B = −120.3), (in Fig. 3 the point marked “Sim” in the region OS). The average rate in this state is νsim = 5.87 ± 0.23 Hz, while theory gives νth = 5.51 Hz. For this state Var(T )/T¯ 2 = 0.0048, which should be contrasted with the corresponding value in the SS state. For the oscillation amplitude we find theory: 11.17 Hz, simulation: 11.66 ± 1.4 Hz. 4.
257
have certain attractive features, such as implementability in VLSI and some potential for modeling the effects of slow NMDA currents. In the latter we refer not to the saturation effects of NMDA, but rather to the fact that NMDA introduces a long time constant into the dynamics, which adds stability and reduces fluctuations and a weaker dependence on the finite size of the network (Brunel and Wang, 2001). Networks of LIF neurons have a similar characteristic. Several items are on the immediate agenda for networks of this type. Given that the neurons have apparently long time constants, what would be the response time of such a network to different perturbations from a stationary state? Simulations lead us to believe that such times will be relatively short and the technique can be applied to compute these times explicitly. In addition, the dynamical characteristics of working memory states, along the lines of Brunel (2000a), should be clarified and compared with simulations. Finally the behavior of such networks should be compared more carefully with that of exponential neurons networks to clarify if this model captures some aspects of neuromodulation. Appendix A A.1. Linear Stability The eigenmodes of Eq. (17) would be Q 1 (V¯ , t) = eλt/τ nˆ 1 (λ) Qˆ 1 (V¯ ), n 1 (t) = nˆ 1 (λ)eλt/τ , where Qˆ 1 (V¯ ) obeys the ordinary differential equation (λ − L) Qˆ 1 (V¯ ) = e−λδ/τ L AB [Q 0 ](V¯ )
(23)
together with the boundary conditions ∂ Qˆ 1 ¯ = 0, ¯ = −1 + Ae−λδ/τ , Qˆ 1 (θ) (θ) ¯ ∂ V 1 ∂ Qˆ 1 A ∂ Q 0 −λδ/τ − Qˆ 1 + e + BQ0 2 ∂ V¯ 2 ∂ V¯
(24)
V¯ =0
1 =− . 2 (25)
Discussion
The solution of the homogeneous part of Eq. (23) is By studying the various dynamical regimes of a purely inhibitory network of LIF neurons, we have set the ground for a more detailed investigation of networks of neurons of this type in more generic and structured situations, as was done for networks of exponential IF neurons (Brunel, 2000a). Networks of LIF neurons
) Qˆ (H = C1 φ1 (V¯ ) + C2 φ2 (V¯ ), 1
with ¯
φi = e zi V
and
zi = 1 ∓
√
1 + 2λ.
(26)
258
Mongillo and Amit
A particular solution of Eq. (23) is e−λδ/τ Qˆ 1 p = L AB [Q 0 ]. λ
At the onset of an oscillatory instability, the eigenvalue of the largest real part is purely imaginary, λ = iωc , where ωc is the frequency of the emerging oscillation:
The boundary conditions in (24) give a linear system for C1 and C2
e−iωc δ/τc ∓(1 − i)B √ = −1 + Ae−iωc δ/τc . ωc
e−λδ/τ L AB [Q 0 ](θ¯ ) = α(λ) λ (27) −λδ/τ ∂ e C1 φ1 (θ¯ ) + C2 φ2 (θ¯ ) = −1 + Ae−λδ/τ − λ ∂ V¯ × L AB [Q 0 ](θ¯ ) ≡ β(λ), (28)
C1 φ1 (θ¯ ) + C2 φ2 (θ¯ ) = −
αz 2 − β −1 αz 1 − β −1 ¯ φ1 (θ¯ ), C2 = − φ (θ). z2 − z1 z2 − z1 2 (29)
The boundary condition (25), ensuring equal flux at the threshold and at the reflecting lower barrier, gives Eq. (18) for the eigenfrequencies A.2. Limit δ/τ → 0 In this limit, the eigenvalues |λ| → ∞. That this is the case can be qualitatively understood as follows: the oscillations are due to the fact that an increase of activity in the network provokes an increase in the feedback inhibitory input after about one synaptic (delay) time. The effect of increased recurrent inhibitory input is a decrease in global activity. This argument predicts √ ω ∼ τ/δ. When |λ| → ∞, z 1,2 ∓ 2λ (Eq. 26). In Eq. (18) one of the exponentials dominates, and it reduces to αz = β,
(30)
with α(λ) and β(λ) given in Eq. (19). z is positive (negative) for positive(negative) drift. We then note √ that for Eq. (30) 2to have such a root, we need B ∼ |λ|. Since σl A = σ 2 +σ 2 < 1, one can neglect the term proportional ext l to A in the left-hand side and the term 2α in the righthand side. The equation becomes B
√ Bc = ∓ ωc sin(ωc δ/τc ) , Ac = sin(ωc δ/τc ) + cos(ωc δ/τc )
for negative and positive drift regimes, respectively. Appendix B
which leads to C1 =
Separating real and imaginary parts, one has
e−λδ/τ z = −1 + Ae−λδ/τ . λ
To determine the amplitude of the collective oscillation near (Ac , Bc ), we expand in powers of the deviation from the bifurcation line. The following discussion has a large overlap with that of Brunel and Hakim (1999). It is reproduced here for the sake of completeness in presence of significant technical differences. Hopf’s theorem ensures that the development parameter is of order of the square root of this displacement; therefore, A(= A − Ac ), B(= B − Bc ) enter only at the third-order. First we determine the second-order term. At the second order Eq. (11) becomes τc
∂ Q2 = L[Q 2 ] + n 2 (t − δ)L AB [Q 0 ] ∂t + n 1 (t − δ)L AB [Q 1 ],
(31)
with the boundary conditions ¯ t) = 0, Q 2 (θ, ∂ Q2 ¯ t) = − n 2 (t) + Ac n 2 (t − δ) (θ, ∂ V¯ − A2c n 21 (t − δ) + Ac n 1 (t)n 1 (t − δ), 1 ∂ Q2 − Q 2 + n 2 (t − δ)K AB [Q 0 ] 2 ∂ V¯ + n 1 (t − δ)K AB [Q 1 ] V¯ =0
1 = − n 2 (t), 2
(32)
where K AB = A2 ∂∂V¯ + B. The forcing terms on the righthand side of Eq. (31) contains terms at frequencies 2
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons ωc and 0; therefore,
together with the boundary conditions
Q 2 (V¯ , t) = |nˆ 1 |2 Qˆ 20 + nˆ 21 e2iωc t/τc Qˆ 22 + c.c.
(33)
n 2 (t) = ρ20 |nˆ 1 |2 + ρ22 e2iωc t/τc nˆ 21 + c.c.
(34)
Substituting into (31) one obtain an ordinary differential equation for Qˆ 22 — (2iωc − L) Qˆ 22 = ρ22 e−2iωc δ/τc L AB [Q 0 ] + e−iωc δ/τc L AB [ Qˆ 1 ]
(35)
—with the boundary conditions
A2c e−2iωc δ/τc
+ Ac e
¯ = 0, (41) Qˆ 20 (θ) ˆ ∂ Q 20 ¯ = ρ20 (−1 + Ac ) + 2Ac cos(ωc δ/τc ) − A2c (θ) ∂ V¯ and 1 ∂ Qˆ 20 ∗ − Qˆ 20 + ρ20 K AB [Q 0 ] + e−iωc δ/τc K AB 2 ∂ V¯ (hom) 1 (hom) iω δ/τ + e c c K AB Qˆ 1 × Qˆ 1 = − ρ20 . 2 ¯ V =0
Qˆ 22 (θ¯ ) = 0 ∂ Qˆ 22 (θ¯ ) = ρ22 −1 + Ac e−2iωc δ/τc ∂ V¯ −
259
(42) We look for a solution of the form
−iωc δ/τc
ˆ (lo) Qˆ 20 = Qˆ (hom) + ρ20 Qˆ (so) 20 20 + Q 20 ,
1 ∂ Qˆ 22 − Qˆ 22 + ρ22 e−2iωc δ/τc × 2 ∂ V¯ × K AB [Q 0 ] + e−iωc δ/τc K AB [ Qˆ 1 ]
where ¯ Qˆ (hom) = D1 + D2 e 2 V 20
V¯ =0
1 = − ρ22 . 2
¯ 2(V¯ −θ¯ ) Qˆ (so) 20 = (Ac + Bc ) V e e−iωc δ/τc ∗ ˆ (hom) Qˆ (lo) L AB Q 1 20 = iωc eiωc δ/τc − L AB Qˆ (hom) . 1 iωc
(36)
The general solution of Eq. (35) is Qˆ 22 = Qˆ (hom) + 22 (lo) ˆ + Q , where ρ22 Qˆ (so) 22 22 Qˆ (hom) (V¯ ) = γ1 ψ1 (V¯ ) + γ2 ψ2 (V¯ ) e−2iωc δ/τc Qˆ (so) = L AB [Q 0 ] 22 2iωc
(37) (38)
(43) (44)
(45)
As above, one determines the three unknowns D1 , D2 and ρ20 using Eqs. (41) and (42). Third Order
and ¯
ψi (V¯ ) = eξi V ; ξ1,2 = 1 ∓
1 + 4iωc .
Since [L , L AB ] = 0 and L[Q 0 ] = 0, (hom) e−2iωc δ/τc 2 e−iωc δ/τc Qˆ (lo) = L L [Q 0 ]. Qˆ 1 + AB 22 iωc 2(iωc )2 AB The boundary conditions in (36) determine γ1 and γ2 , in terms of ρ22 and the previously determined functions, as well as a linear algebraic equation for ρ22 . The component Qˆ 20 obeys −L[ Qˆ 20 ] = ρ20 L AB [Q 0 ] + e−iωc δ/τc L ∗AB [ Qˆ 1 ] (39) + eiωc δ/τc L AB [ Qˆ 1 ] (hom) −iωc δ/τc ∗ = ρ20 L AB [Q 0 ] + e L AB Qˆ 1 + eiωc δ/τc L AB Qˆ (hom) , (40) 1
At the third order, the naive expansion fails due to the appearance of a resonant term at frequency ωc , leading to secular terms that are an artifact of the truncation procedure (see, e.g., Bender and Orszag, 1987). The third-order terms obey τc
∂ Q3 = L[Q 3 ] + n 3 (t − δ)L AB [Q 0 ] ∂t + n 2 (t − δ)L AB [Q 1 ] + n 1 (t − δ)L AB [Q 2 ] 1 ∂ 2 Q0 + n 1 (t − δ) (A − Ac ) 2 ∂ V¯ 2 ∂ Q0 d nˆ 1 iωc τ + (B − Bc ) Q1 − − ¯ τc dt ∂V iωc t/τc iωc (t−δ)/τc ˆ + δe L AB [Q 0 ] . × τc Q 1 e (46)
260
Mongillo and Amit
The boundary conditions are Q 3 (θ¯ , t) = 0 ∂ Q3 (θ¯ , t) = −n 3 (t) + Ac n 3 (t − δ) ∂ V¯
(47)
together with boundary conditions:
− 2A2c n 1 (t − δ)n 2 (t − δ)
+ Ac (n 1 (t)n 2 (t − δ) + n 1 (t − δ)n 2 (t)) + A3c n 31 (t − δ) − A2c n 1 (t)n 21 (t − δ)
+ (A − Ac )n 1 (t − δ) d nˆ 1 iωc (t − δ)/τc − Ac δ e dt
(48)
K [Q 3 ] + n 3 (t − δ)K AB [Q 0 ] + n 2 (t − δ)K AB [Q 1 ] + n 1 (t − δ)K AB [Q 2 ] + n 1 (t − δ) 1 ∂ Q0 × (A − Ac) + (B − Bc )Q 0 2 ∂ V¯ d nˆ 1 iωc (t−δ)/τc −δ e K AB [Q 0 ] dt V¯ =0 1 = − n 3 (t). (49) 2 The terms in ddtnˆ 1 are due to the slow time scale, T ≡ t, needed to cancel secular terms. The first term arises from explicit time differentiation; the second comes from a first order expansion in of n(t − δ) = n(t − δ, T − δ). This term appears also in the boundary condition. Variation of A and B from Ac , Bc imply an additional variation of τ ; formally, τ (A, B) = τc (Ac , Bc ) + τ (A, B), where τ is linear in A and B. Indeed, τ depends on µ0 and σ02 , which, in turn, can be expressed in terms of A and B. The forcing terms on the right-hand side of Eq. (46) oscillate at frequencies 3ωc and ωc . Therefore, Q 3 (V¯ , t) = e3iωc t/τc Qˆ 33 + eiωc t/τc Qˆ 31 + c.c. n 3 (t) = nˆ 33 e3iωc t/τc + nˆ 31 eiωc t/τc + c.c. The terms at frequency ωc are resonant with the firstorder terms: (iωc − L) Qˆ 31 = nˆ 31 e−iωc δ/τc L AB [Q 0 ] + nˆ 1 |nˆ 1 |2 × ρ22 e−2iωc δ/τc L ∗AB [ Qˆ 1 ] + ρ20 L AB [ Qˆ 1 ] + e−iωc δ/τc L AB [ Qˆ 20 ] + eiωc δ/τc L AB [ Qˆ 22 ]
iωc τ ˆ Q1 + nˆ 1 e−iωc δ/τc L A−Ac ,B−Bc [Q 0 ] − τc d nˆ 1 ˆ − τc (50) ( Q 1 + δ/τc e−iωc δ/τc L AB [Q 0 ]), dt
¯ =0 Qˆ 31 (θ) ∂ Qˆ 31 ¯ = nˆ 31 (−1 + Ac e−iωc δ/τc ) + nˆ 1 1 (θ) ∂ V¯ d nˆ 1 + nˆ 1 |nˆ 1 |2 3 + (51) d dt 1 ∂ Qˆ 31 − Qˆ 31 + nˆ 31 e−iωc δ/τc K AB [Q 0 ] 2 ∂ V¯ V¯ =0 1 d nˆ 1 = − nˆ 31 + F nˆ 1 + G nˆ 1 |nˆ 1 |2 + H , 2 dt
(52)
in which we have denoted by 1 , 3 , d the coefficients of the linear, third-order, and derivative terms, respectively. The general solution of Eq. (50) is Qˆ 31 (V¯ ) ¯ = K¯ 1 φ1 (V¯ ) + K¯ 2 φ2 (V¯ ) + nˆ 31 Qˆ 1 p (V¯ ) + Qˆ (lo) 31 ( V ), where ψ1 , ψ2 , and Qˆ 1 p are the functions that appear at first order, and d nˆ 1 ¯ Qˆ (lo) Q d (V¯ ) + nˆ 1 Q l (V¯ ) + nˆ 1 |nˆ 1 |2 Q c (V¯ ). 31 ( V ) = dt Using (51), one obtains ¯ + K¯ 2 φ2 (θ) ¯ = −nˆ 31 Qˆ 1 p (θ) ¯ − Qˆ (lo) ¯ K¯ 1 φ1 (θ) 31 (θ) (53) ¯ + K¯ 2 φ2 (θ) ¯ K¯ 1 φ1 (θ) ∂ ˆ ¯ = nˆ 31 −1 + Ac e−iωc δ/τc − Q 1 p (θ) ∂ V¯ ∂ ˆ (lo) ¯ + , − (54) Q 31 (θ) ∂ V¯ where stands for the sum of the last three terms in Eq. (52). The coefficients of nˆ 31 in the right-hand side of (53) are respectively α(iωc ) and β(iωc ), defined in Eq. (19). Therefore, K¯ 1 = nˆ 31 C1 + K 1 with C1 is defined in Eq. (29) and K 1 is the remainder ˆ 31 Qˆ 1 + Qˆ 3 . Qˆ 31 = nˆ 31 Qˆ 1 + K 1 φ1 + K 2 φ2 + Qˆ (lo) 31 = n (55)
Oscillations and Irregular Emission in Networks of Linear Spiking Neurons
Using (52), one finally obtains 1 ∂ Qˆ 1 −iωc δ/τc ˆ nˆ 31 K AB [Q 0 ] − Q1 + e 2 ∂ V¯ 1 ∂ Qˆ 3 1 + − Qˆ 3 = − nˆ 31 + F nˆ 1 + G nˆ 1 |nˆ 1 |2 2 ∂ V¯ 2 V¯ =0
d nˆ 1 . dt Noting that the term multiplying nˆ 31 is the eigenvalue Eq. (18), the equation reduces to 1 ∂ Qˆ 3 ˆ d nˆ 1 − Q3 . (56) = F nˆ 1 +G nˆ 1 |nˆ 1 |2 + H 2 ∂ V¯ dt ¯ V =0 +H
When Qˆ 3 is expressed, from Eq. (55) in terms of the first order, third order, and the derivative of nˆ 1 , one arrives at the (Landau) equation of motion for nˆ 1 : d nˆ 1 = W1 nˆ 1 − W3 nˆ 1 |nˆ 1 |2 (57) dt In the limit δ¯ → 0 the expressions for W1 and W3 simplify: Eq. (56) becomes τc
∂ ˆ (lo) ¯ Qˆ (lo) Q 31 (θ¯ ) + = 0, 31 (θ )z − ∂ V¯ √ where z = ± 2iωc . For W1 and W3 , one obtains Q l (θ¯ )z − W1 = − Q d (θ¯ )z − Q c (θ¯ )z − W3 = Q d (θ¯ )z −
∂ Ql ¯ (θ ) + 1 ∂ V¯ ∂ Qd ¯ (θ ) + d ∂ V¯ ∂ Qc ¯ (θ ) + 3 ∂ V¯ . ∂ Qd ¯ (θ ) + d ∂ V¯
, (58)
Acknowledgments We are indebted to Dr. Nicolas Brunel for very helpful discussions and to Dr. Massimo Mascaro for providing us with his simulation program. This study was supported by a Center of Excellence grant “Statistical Mechanics and Complexity” of the INFM, Roma, Italia and Center of Excellence grant “changing your mind” of the Israel Science Foundation. Notes 1. We follow the convention of taking the afferent current for the RC neuron in units of potential. 2. To eliminate the effect of the oscillations on the estimate of the variance in the oscillating states, we have averaged the rate in entire cycles and calculated the variance of the rate from cycle to cycle.
261
References Amit DJ, Brunel N (1997) Dynamics of a recurrent network of spiking neurons before and following learning. Network Comput. Neural Sys. 8: 373–404. Amit DJ, Brunel N (1997) Model of global spontaneous activity and local structured (learned) delay activity during delay periods in cerebral cortex. Cerebral Cortex 7: 237. Bender CM, Orszag SA (1987) Advanced Mathematical Methods for Scientists and Engineers. McGraw-Hill, Singapore. Brunel N (2000) Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neuroscience 8: 183. Brunel N (2000) Persistent activity and the single cell frequencycurrent curve in a cortical network model. Network 11: 261. Brunel N, Hakim V (1999) Fast global oscillations in networks of integrate-and-fire neurons with low firing rates. Neural Comput. 11: 1621. Brunel N, Wang XJ (2001) A cortical network model of object working memory: Resistance to distractors and neuromodulation. J. Comput. Neurosci. 11: 63. Buzsaki G, Hovrath Z, Urioste R, Hetke J, Wise K (1992) Highfrequency network oscillations in hippocampus. Science 256: 1025. Compte A, Brunel N, Goldman-Rakic PS, Wang XJ (2000) Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model. Cerebral Cortex 10: 910. Fusi S, Mattia M (1999) Collective behavior of networks with linear (VLSI) integrate and fire neurons. Neural Comput. 11: 633. Fusi S, Annunziato M, Badoni D, Salamon A, Amit DJ (2000) Spikedriven synaptic plasticity: Theory, simulation, VLSI implementation. Neural Comput. 12: 2227. Gerstein GL, Mandelbrot B (1964). Random walk models for the spike activity of a single neuron. Biophysical J. 4: 41. Gerstner W (1995) Time structure of the activity in neural network models. Phys. Rev. E 51: 738. Kahn PB, Zarmi Y (1998) Nonlinear Dynamics: Exploration Through Normal Forms. Wiley, New York. Knight BW (1972) Dynamics of encoding in a population of neurons. J. Gen. Physiol. 59: 734. Rauch A, La Camera G, Fusi S, Senn W, Luscher HR (2001) Discharge rate of neocortical pyramidal cells in response to noisy input currents. Submitted for oral presentation to CNS2001. Mascaro M (1998) Tesi di laurea, University of Roma La Sapienza. (in Italian). Funzione Di Risposta Efficace Di Una Sotto Popolazione Di Neuroni In Una Rete Multimodulare. Mascaro M, Amit DJ (1999) Effective neural response function for collective population states. Network 10: 351. Mattia M (1997), Tesi di laurea, University Roma La Sapienza. (in Italian). Dinamica Di Una Rete Di Neuroni Impucsivi Con Depolarizzazione Lineave. Mead C (1989) Analog VLSI and Neural System. Addison-Wesley, Reading, MA. Traub RD, Whittington MA, Collins SB, Buzsaki G, Jefferys JE (1996) Analysis of gamma rhythms in the rat hippocampus in vitro and in vivo. J. Physiol. 493: 471. Van Vreeswijk C, Abbott L, Ermentrout GB (1994) When inhibitio not excitation synchronizes firing. J. Comput. Neurosci. 1: 313.
