1

Network Dynamics Group, Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany

2

Bernstein Center for Computational Neuroscience, Göttingen, Germany

3

Center for Brain Science, Faculty of Arts and Sciences, Harvard University, Cambridge, MA, USA

Large sparse circuits of spiking neurons exhibit a balanced state of highly irregular activity under a wide range of conditions. It occurs likewise in sparsely connected random networks that receive excitatory external inputs and recurrent inhibition as well as in networks with mixed recurrent inhibition and excitation. Here we analytically investigate this irregular dynamics in finite networks keeping track of all individual spike times and the identities of individual neurons. For delayed, purely inhibitory interactions we show that the irregular dynamics is not chaotic but stable. Moreover, we demonstrate that after long transients the dynamics converges towards periodic orbits and that every generic periodic orbit of these dynamical systems is stable. We investigate the collective irregular dynamics upon increasing the time scale of synaptic responses and upon iteratively replacing inhibitory by excitatory interactions. Whereas for small and moderate time scales as well as for few excitatory interactions, the dynamics stays stable, there is a smooth transition to chaos if the synaptic response becomes sufficiently slow (even in purely inhibitory networks) or the number of excitatory interactions becomes too large. These results indicate that chaotic and stable dynamics are equally capable of generating the irregular neuronal activity. More generally, chaos apparently is not essential for generating the high irregularity of balanced activity, and we suggest that a mechanism different from chaos and stochasticity significantly contributes to irregular activity in cortical circuits.

Most neurons in the brain communicate by emitting and receiving electrical pulses, called action potentials or spikes, via chemically operating synaptic connections. Local cortical circuits often exhibit spiking dynamics that is highly irregular and appears as if it were random. Such irregular activity at low neuronal firing rate is thus considered a basic “ground state”. It is characterized by individual neurons that display largely fluctuating membrane potentials and highly variable inter-spike-intervals (ISIs) as well as by low correlations between the neurons (v.Vreeswijk and Sompolinsky, 1996
, 1998
; Brunel, 2000
; Vogels and Abbott, 2005
; Kumar et al., 2007
). Originally, this dynamical state seemed to be in contradiction to cortical anatomy, where each neuron receives a huge number of synapses, typically 10

^{3}–10^{4}(Braitenberg and Schüz, 1998 ): One might expect that a large number of uncorrelated, or weakly correlated synaptic inputs to one neuron, given the central limit theorem, sums up to a regular total input signal with only small relative fluctuations, therefore excluding the emergence of irregular dynamics. So the finding of highly irregular activity might be surprising.This issue was resolved by the idea of a “balanced state” (v.Vreeswijk and Sompolinsky, 1996
), in which excitatory (positive) and inhibitory (negative) input balances such that the average membrane potential is sub-threshold and strong fluctuations once in a while are sufficiently depolarizing to initiate a spike. The original work “Chaos in neuronal networks with balanced excitatory and inhibitory activity” (v.Vreeswijk and Sompolinsky, 1996
) was an analysis of a self-consistent, highly irregular “balanced” activity for sparse random networks of binary model neurons. It was shown that the balanced state naturally and robustly occurs in large networks if the inhibitory and excitatory coupling strengths and their respective numbers of synapses appropriately scale with each other. Moreover, flipping the state of only one of the binary neurons in a large network (i.e. applying the smallest possible non-zero perturbation in such a system) leads to a supra-exponential divergence between the perturbed and the unperturbed realizations of the network dynamics, exemplifying the extremely chaotic nature of the balanced activity (v.Vreeswijk and Sompolinsky, 1996
, 1998
).

Later analysis (Amit and Brunel, 1997
; Brunel, 2000
) of networks of integrate-and-fire neurons demonstrated that the mean field description of the balanced activity for these continuous-state neuron networks is very similar to that of the original binary neuron networks. These and related results (Amit and Brunel, 1997
; Brunel, 2000
; Timme et al., 2002
; Timme and Wolf, 2008
) indicate that statistically the same balanced activity persists both in networks with external excitatory inputs and recurrent inhibition only as well as in networks with equal total amounts of recurrent inhibition and recurrent excitation. Given these findings about the robustness of the balanced state, the original work (v.Vreeswijk and Sompolinsky, 1996
, 1998
) together with common intuition may suggest that highly irregular activity originates from chaotic network dynamics. This hypothesis, however, has not been systematically investigated so far. Recent research even points towards the contrary: it shows that in globally coupled networks without delay the dynamics tends to converge to stable periodic orbits if inhibition dominates (Jin, 2002
). Numerical investigations of weakly diluted inhibitorily coupled networks without delay show that although the dynamics may be irregular, its Lyapunov exponent is negative (Zillmer et al., 2006
). These numerical simulations also demonstrate by example that the dynamics converges to a periodic orbit after long quasi-stationary transients. Interestingly, a related article (Zillmer et al., 2009
) also gives numerical evidence that there can be chaos and long transients even in networks with only inhibitory connections.

In this article we show analytically in the limit of fast synaptic response, that in inhibitory networks with inhomogeneous delay distribution and arbitrary, strongly connected topology (A network is

*strongly connected*if there is a directed path of connections between any ordered pair of neurons.) any generic trajectory is asymptotically stable. After a (typically long) stable transient characterized by irregular activity the dynamics converges to a periodic orbit that is also stable, in agreement with the results presented by Memmesheimer and Timme (2006a) and Timme and Wolf (2008) . In particular the transients are not chaotic in contrast to the ones occurring in purely excitatorily coupled networks (Zumdieck et al., 2004 ). We show that this collective dynamics is robust upon increasing the synaptic response time from zero and upon replacing some inhibitory by excitatory interactions. Nevertheless, if the synaptic response becomes too slow or the number of excitatory interaction too large, the collective dynamics becomes chaotic via a transition where the Lyapunov exponent changes smoothly and the spiking activity stays highly irregular. Thus the irregularity equally prevails in networks with stable as well as in networks with chaotic dynamics (Figure 1 ), leaving no evidence that chaos generates the irregularity.**Figure 1. Highly irregular spiking activity equally emerges from chaotic and from stable circuit dynamics.**Irregular dynamics in purely inhibitorily

**(A-C)**and inhibitorily and excitatorily

**(D–F)**coupled random networks of identical leaky integrate-and-fire neurons (

*N*= 400, γ

*≡ 1,*

_{i}*V*

_{Θ,i}≡ 1, |Pre(

*i*)| ≡ 80).

**(A,D)**Spiking dynamics

**(A)**

*I*≡ 4,

_{i}*N*= 0;

_{E}**(D)**,

*N*= 10000, where the

_{E}*N*excitatory couplings are distributed such that each neuron has the same number of excitatory inputs. The upper panel displays the spiking times (blue lines) of the first 40 neurons. The lower panel displays the membrane potential trajectory of neuron

_{E}*i*= 1 (spikes of height Δ

*V*= 1 added at firing times).

**(B,E)**Histogram of mean firing rates ν

*.*

_{i}**(C,F)**Histogram of the coefficients of variation CV

*:= σ*

_{i}*/μ*

_{i}*; averaged over time.*

_{i}Some analytically accessible aspects of the stable irregular dynamics have been briefly reported before (Jahnke et al., 2008a
). Parts of this work have been presented by Memmesheimer (2007) and at a Bernstein Symposium (Jahnke et al., 2008b
).

We consider networks of

*N*neurons with directed couplings that interact by sending and receiving spikes. If such a directed connection exists from neuron*j*to neuron*i*, we call*i*postsynaptic to neuron*j*. We denote the set of all postsynaptic neurons of neuron*j*by Post(*j*). Neuron*j*is then presynaptic to neuron*i*, the set of all presynaptic neurons is denoted by Pre(*i*). The membrane potential*V*(_{i}*t*) of some neuron*i*evolves according to:where a smooth function

*f*specifies the internal dynamics, ε_{i}_{ij}are the coupling strengths from presynaptic neurons*j*to*i*, δ(.) is the Dirac delta-distribution, τ_{ij}> 0 are the delay times of the connection and denotes the time when neuron*j*sends the*k*th spikes. If not stated otherwise, in the following we consider inhibitory networks, i.e. ε_{ij}≤ 0. Spikes sent to postsynaptic neurons with different delay times from neuron*j*are considered as separate spikes, so if the delays are all different, e.g. if they are chosen randomly, neuron*j*sends |Post(*j*)| spikes at time When neuron*j*reaches the threshold*V*_{Θ,j}of the potential, i.e.*V*(_{j}*t*^{−}) =*V*_{Θ,j}, it generates spikes at for some*k*and is reset, The neuronal dynamics is therefore smooth except at times when events, namely sendings or receivings of spikes happen. Simultaneous sendings of spikes by one neuron are treated as one event as well as simultaneous receptions of spikes sent by the same neuron. We require that the*f*satisfy_{j}*f*(_{j}*V*) > 0 and for all_{j}*j*and*V*≤_{j}*V*_{Θ,j}, such that in isolation each neuron*j*exhibits oscillatory dynamics with a periodThe network dynamics can equivalently be described by a phase-like variable ϕ

_{j}(*t*) ∈ (−∞,ϕ_{Θ,j}] satisfyingat all non-event times (Mirollo and Strogatz, 1990
). When the phase threshold is reached, the phase is reset, and a spike is generated. This spike travels to the postsynaptic neurons, arrives after a delay time τ

_{ij}at neuron*i*and induces a phase change according towith the transfer function

where each

*U*(_{i}*t*) is the free (all ε_{ij}= 0) solution of (1) through the initial condition*U*(0) = 0, yielding_{i}*U*′ > 0 and_{i}*U*″ < 0 and (cf. Memmesheimer and Timme, 2006a ). Figure 2 illustrates the relation between phase dynamics and membrane potential._{i}**Figure 2. Relation between membrane potential and phase dynamics.**

**(A)**Membrane potential of neuron 1. At the membrane potential

*V*

_{1}(

*t*) crosses the threshold potential

*V*

_{Θ,1}which leads to a reset of the potential to and spikes are emitted.

**(B)**Membrane potential of neuron 2 which is postsynaptic to neuron 1, i.e. 2 ∈ Post(1). In this example the coupling is excitatory. The spike is received at and induces a jump of size

*ε*

_{21}in the potential.

**(C)**Phase dynamics of neuron 1. The phase increases linearly until the phase threshold ϕ

_{Θ,1}is reached, then it is set to zero.

**(D)**Phase dynamics of neuron 2. When the spike is received at it induces a phase jump: .

The analysis below is valid for general

*U*(ϕ); in the numerical simulations we employ leaky integrate-and-fire neurons,_{i}*f*(_{i}*V*) :=*I*− γ_{i}_{i}*V*with time scale and equilibrium potential the membrane potential has the functional dependence on ϕ and the oscillation period of a free neuron is given by We consider arbitrary generic spike sequences in which all neurons are active, i.e. there is a finite constant*T*> 0, such that in every time interval [*t*,*t*+*T*),*t*∈ ℝ, every neuron fires at least once. Further, we assume that the dynamics is sufficiently irregular such that two events occur at the same time with zero probability.Due to the delay, the state space is formally infinite dimensional. However, it becomes finite dimensional after some finite time

*t*′ (cf. Ashwin and Timme, 2005 ). At a given time*t*>*t*′ the network dynamics is completely determined by the phases ϕ_{i}(*t*) and by the spikes which have been sent but not yet received by the postsynaptic neurons at*t*. Their number is bounded by some constant*ND*′. Due to the inhibitory character of the network couplings, each neuron*i*needs at least the free oscillation period to generate a spike after the last reset. Consequently, at mostspikes per neuron are in transit and the state space stays finite, with dimensionality smaller than or equal to

*N*·(1 +*D*′) (cf. also Ashwin and Timme, 2005 ). Here ⌈*x*⌉ denotes the ceiling function, the smallest integer larger or equal to*x*.We now introduce variables to describe spikes, which are already sent at time

*t*by neuron*j*to the postsynaptic neuron*i*and not yet received. A single spike in transit is characterized by the state variable σ_{ijk}(*t*) ∈ [0, τ_{ij}]. The index*k*= 1, 2, 3…≤*D*′/*N*numbers the different spikes traveling from neuron*j*to*i*at time*t*in the order of arrival at the postsynaptic neuron*i*, starting with*k*= 1 for the next spike to arrive. When spikes are emitted at time for some*n*, σ_{ijk}(*t*) is set to The spike index*k*equals the number of spikes already in transit plus one. It thus depends on the actual network state at time Between two events σ_{ijk}(*t*) increases linearly with slope one, when the spike is received by the postsynaptic neuron*i*, where it induces a phase jump according to Eq. 3. After spike reception we cancel the spike arrived (which has index*k*= 1) and renumber the indices*k*> 1 as*k*→*k*− 1 such that specifies the spike sent by neuron*j*which arrives next at the postsynaptic neuron*i*(cf. Figures 3 A,B for illustration).**Figure 3. Phase dynamics of neuron j and definition of state variables.**

**(A)**At

*t*=

*t*and

_{n}*t*=

*t*

_{n + 1}the phase ϕ

_{j}(

*t*) reaches the threshold and is reset to 0.

**(B)**The spikes emitted travel to the postsynaptic neurons

*i*∈ Post(

*j*). They are described by σ

*(*

_{ijk}*t*). In this example we show two spikes traveling from neuron

*j*to one specific postsynaptic neuron

*i*, described by σ

*(*

_{ij1}*t*) (black) and σ

*(*

_{ij2}*t*) (red). At

*t*=

*t*, σ

_{n}*(*

_{ij1}*t*) is set to 0; here

_{n}*k*= 1 because there is no spike in transit at When neuron

*j*spikes again at

*t*=

*t*, σ

_{n + 1}*(*

_{ij2}*t*) is set to 0. At

_{n + 1}*t*=

*t*+ τ

_{n}*the spike emitted at*

_{ij}*t*=

*t*arrives at the postsynaptic neuron

_{n}*i*and induces a phase jump in ϕ

*(*

_{i}*t*) (not shown, cf. Figure 2 ). After spike reception, we renumber

*k*→

*k*− 1 such that σ

_{ij2}→ σ

_{ij1}.

**(C)**Definition of the phase shifts. The phase curve ϕ

_{i}(

*t*) of neuron

*i*(blue) before and after the reception of a spike at

*t*=

*t*is shown together with (black). is the difference of neuron

_{n}*i*’s phases in the unperturbed and the perturbed dynamics taken at corresponding event times

*t*and . denotes the difference of event times

_{n}*t*and i.e. the temporal offset between both sequences. Finally, is some phase shift of neuron

_{n}*i*with the temporal offset taken into account.

### Lyapunov Stability of Arbitrary Generic Spike Sequences

In this section we study the stability properties of the spike sequences. We compare the microscopic dynamics of two sequences, that slightly differ in the timing of the spikes, but have the same ordering. We show that the distance between these trajectories is bounded by the initial distance. Assuming that one sequence is generated by a perturbation of the other, this implies Lyapunov stability for the considered spike pattern. Distances and perturbation sizes are measured using the maximum norm.

Since the distance between trajectories only changes at event times, we can choose an event-based perspective. The time of the

*n*th event in the entire network is denoted by*t*in the first sequence, and by in the second one. Analogously, we denote the phases of a neuron_{n}*i*at a given time*t*by ϕ_{i}(*t*) and and the spikes in transit at time*t*by σ_{ijk}(*t*) and in the different sequences. Letdenote the difference between the phase difference, and the temporal offset, after the

*n*th and before the (*n*+ 1)th event (cf. Figure 3 C). Similarlylabels the shift of the

*k*th spike sent by neuron*j*and not yet arrived at neuron*i*after the*n*th and before the (*n*+ 1)th event. Between two consecutive events, both and increase linearly and only at event times the phases and spike variables are updated nonlinearly as described above. Therefore, to study the stability of the system it is sufficient to consider the phase shifts after events.In the following we investigate the evolution of shifts at the discrete event times. There are two different kinds of events: (i) sending and (ii) receiving of spikes. In the first case, the shifts and stay unchanged, but new spikes with new spike variables are generated. These variables inherit the perturbation of the sending neuron. In the second case, the phase shift of the neuron which receives the spike changes and the spikes in transit are reordered. The resulting phase shift of the neuron receiving a spike turns out to be a weighted sum of previous shifts. This can be shown by studying both cases in detail.

#### Transfer of perturbations without change of size

If as (

*n*+ 1)th event the phase of some neuron*j** reaches its threshold and a spike is emitted, the shifts of all neurons’ phases stay unchanged,Similarly, the shifts of the spikes in transit stay unchanged

Additionally, new spikes are generated σ

_{ij*k*}(*t*_{n + 1}) = 0 and where*k** =*k**(*i*,*j**,*n*+ 1) is the appropriate spike number, cf. Figure 3 . The shifts of the new spike variables depend on the phase shift of the sending neuron*j** according to#### Averaging of prior perturbations

If as (

*n*+ 1)th event some spike arrives, say σ_{i*j*1}(*t*_{n + 1}) = τ_{i*j*}, it induces a phase jump in the postsynaptic neuron*i**. According to Eq. 3, the phase shift can be computed as:where and are the phases “just before” spike reception. Using the definitions (6) we find the identity

Applying the mean value theorem in Eq. 11 and the relation

yields

where is given by the derivative

for for ϕ takes values ϕ ∈ If neuron

*i** is not connected to neuron*j**, ε_{i*j*}= 0, the function becomes the identity, such that the phase shift stays unchanged, indeed independent of ϕ.For ε

_{i*j*}< 0, is bounded bywhere and We have used that is monotonic increasing with ε

The shifts of traveling spikes stay unchanged on spike reception, cf. Eq. 9,

for all spikes with

*i*≠*i** ∨*j*≠*j**. For*i*=*i**,*j*=*j** the spike variables are renumbered andholds, except for

*k*= 1, because σ_{i*j*1}(*t*) is the variable describing the spike received and therefore canceled.We will now ascertain that the coefficients in Eq. 14 lie in a compact interval within (0, 1), such that a true averaging takes places when interactions happen. Formally, the phases of neurons can achieve values ϕ

_{i}∈ (−∞, ϕ_{Θ,i}]. Each neuron fires at least once within a time interval of length*T*, therefore the phases are certainly bounded to the compact interval ϕ_{i}∈ [−*T*+ ϕ_{Θ,i}, ϕ_{Θ,i}]. Further, in inhibitory networks the phase after an interaction is smaller than before,because together with

*U*also is strictly monotonic increasing and ε < 0. The strict concavity of_{i}*U*(ϕ) implies_{i}for any finite ϕ (cf. Figure 4
for illustration). The derivative is continuous in ϕ, therefore the image of [−T + ϕ

_{Θ,i},ϕ_{Θ,i}] under the map is compact. Together with Eq. 22 it follows that**Figure 4. (A)**The transfer function for a leaky integrate-and-fire neuron. For ε = 0 (black), is the identity (black), for ε < 0 (blue: ε = −0.5, green: ε = −1) the phase ϕ after receiving an input is smaller than before. For all inhibitory inputs we find

**(B)**The derivative of

*U*(ϕ) is monotonic decreasing, therefore for ε ≤ 0 (cf. Eq. 22).

_{i}Taken together, Eqs. 8–10, 14 and 23 imply that a true averaging between shifts already present in the system takes place when a spike is received. For other events the shifts stay unchanged. As a consequence, the maximal and minimal shift after the

*n*th event,are bounded by the initial shifts for all future events,

as long as the order of events in both sequences is the same. Here the minima and maxima are taken over

*i*,*j*∈ {1,…,*N*} and*k*numbers the spikes traveling from neuron*j*to*i*at time*t*. An initial perturbation cannot grow, thus the trajectory is Lyapunov stable. We note that we did not make any assumptions about the network connectivity, the results hold for any network structure and the described class of trajectories._{n}### Asymptotic Stability

In this section we prove that for strongly connected networks even asymptotic stability holds under the condition that the perturbed and the unperturbed sequences have the same order of events, i.e. the order of events is unchanged by small perturbations. The central idea is as follows: We study the dynamics and convergence of two neighboring trajectories. We will track the propagation of the perturbation of one specific neuron

*l*_{0}through the entire network. Since there is a directed connection between every pair of neurons in the network and any spike reception leads to an averaging of shifts, there is an averaging over all perturbations in the network. For large times all perturbations converge towards the same value, such that both sequences become equivalent, only shifted by a constant temporal offset. Further details and the derivation of the following Eqs. 26–28 is provided in the Appendix.We find the upper bound

for the maximal perturbation after

*K*:= 2*NM*events and analogously the lower boundfor the minimal perturbation. The averaging factor

*c** is determined by the network parameters and bounded toA bound for the difference of the maximal and minimal perturbation after

*K*events is therefore given byThe spread of any perturbation through the network has a contracting effect on the total perturbation, it leads to a decay of the difference between the extremal perturbations at least by a factor (1 −

*c**). Inequality (29) implies together with the bound (28) that for the considered trajectories in the long-time limit the maximal and minimal perturbation are the same,Thus, for

*t*→ ∞ the events are just shifted by a constant temporal offset(cf. Figure 5
for illustration), and both sequences become equivalent Thus all sequences considered are asymptotically stable for all strongly connected networks and perturbations decay exponentially fast with at least

where ⌊·⌋ is the floor function. The actual and numerically measured exponential decay is much faster than the estimation given by Eq. 32, because for deriving Eqs. 26 and 27 in the Appendix we assumed a worst-case scenario. The main reasons for the faster decay are that (i) the mean path length is more meaningful for estimating the number of events until all neurons have received an input from the starting one and (ii) it is impossible that the neurons receive the worst-case perturbation at each reception.

**Figure 5. Stable dynamics in the network of Figures 1 A–C.**

**(A)**Exponential decay of the maximal perturbation (blue dots) and minimal margin κ

^{(n)}(gray line) for one microscopic dynamics. The initial perturbation decays exponentially fast.

**(B)**Exponential convergence of the temporal offset δ

*t*

^{(n)}. For large

*n*the original and the perturbed sequences are equivalent, just shifted by a constant offset

### Margins

The stability results in the previous sections hold for the class of patterns, where a small perturbation does not change the order of spikes. In this section we show that typical spike patterns, generated by networks with a complex connectivity, belong to this class.

In heterogeneous networks with purely inhibitory interactions the occurrence of events at

*identical*times has a zero probability. There is no mechanism causing simultaneous spiking, like supra-threshold inputs in excitatorily coupled networks. As long as two events do not occur at the same event time, there is a non-zero perturbation size keeping the order unchanged in any*finite*time interval. However, the requirement of an unchanged event order yields more and more conditions over time such that the allowed size of a perturbation could decay more quickly with time than the actual perturbation. This is excluded if a temporal*margin*μ^{(n)}(cf. Jin, 2002 ) stays larger than the dynamical perturbation for*infinite*time. Formally, after time*t*denote the_{n}*k*th potential future event time (of the original trajectory) that would arise if there were no future interactions by θ_{n,k},*k*∈ ℕ, and the temporal*margin*by μ^{(n)}:= θ_{n,2}− θ_{n,1}. A sufficiently small perturbation, satisfying cannot change the order of the (*n*+ 1)th event.#### Stability of generic periodic orbits

This directly implies that almost all

*periodic orbits*(all those with non-degenerate event times*t*) consisting of a finite number of_{n}*P*events are stable because there is a minimal marginfor every non-degenerate periodic pattern (c.f. also Memmesheimer and Timme, 2009 under revision).

#### Stability of irregular non-periodic spiking activity

We study the stability properties of

*irregular non-periodic spike sequences*by considering the minimal margin κ^{(n)}over the first*n*events. For simplicity, we consider delay distributions where τ_{ij}is independent of the spike receiving neuron*i*. The irregular spiking dynamics of the entire network is well modeled by a Poisson point process with rate ν^{s}, where ν^{s}specifies the mean firing rate of the network (Tuckwell, 1988 ; Brunel and Hakim, 1999 ; Brunel, 2000 ; Burkitt, 2006 ). We assume that, along with the irregular dynamics, the temporal margins are also generated by a Poisson point process with the network event rate ν = 2ν^{s}, because any spike sending time generates one receiving event due to the independence of τ_{ij}from*i*and the definition of events. The distribution function of margins is given byTherefore, the probability that the minimal margin κ

^{(n)}after*n*events is smaller or equal to μ is determined by the probabilities that not all individual margins μ^{(n)}are larger than μ such thatwith density ρ

_{n}(μ) := d*P*(κ^{(n)}≤ μ)/dμ =*n*νexp(−*n*νμ). This implies an*algebraic*decay with the number*n*of events for the expected minimal marginthat depends only on the event rate and is independent of the specific network parameters. Numerical simulations show excellent agreement with this algebraic decay (36); a typical example is shown in Figure 6
.

**Figure 6. Margins in the network given in Figures 1 A–C.**

**(A,B)**Probability distribution

*P*(κ

^{(n)}≤ μ) after

**(A)**

*n*= 1 and

**(B)**

*n*= 10

^{5}events. The blue curve shows the distribution measured over 2500 samples, the red dotted line is the analytical prediction (no free fit parameter, rate ν is measured; cf. Eq. 35).

**(C)**Algebraic decay of the average minimal margin, (green dashed line, averaged over 250 random initial conditions) and its analytical prediction (no free fit parameter; black solid line). Additionally we show the minimal margin κ

^{(n)}for three exemplary initial conditions (gray lines), including that of Figure 5 A.

This already strongly suggests that a sufficiently small perturbation stays smaller than the minimal margin for all times. However, in each step, the exponential distribution of margins has a finite density for arbitrary small values of μ, i.e. in each step the margin can fall below the level of perturbation with positive probability. We will show that the probability that there is at least one step in which the margin falls below the perturbation size, goes to zero if the size of the initial perturbation goes to zero. Thus, also the probability that there is a change in the order of events goes to zero. Of course, we cannot expect to reach zero for nonzero perturbation.

We derive a lower bound for the probability that the margin stays larger than for infinite time. We show that it converges to one when the size of the perturbation goes to zero and thus prove that sufficiently small perturbations have arbitrarily high probability to stay smaller than the minimal margin for all times. We start from the upper bound for the evolution of the perturbation, Eq. 32. Using and Eq. 28 leads to

where we introduced

In particular,

*C*→ 0 if the initial perturbation goes to zero, i.e. while α is independent of the initial perturbation. The probability that all margins are larger than all perturbations is given bysince the margins are independent. Using Eq. 37 and

*P*(μ^{(n)}> μ) = exp(−νμ) yieldswhich goes to one if the initial perturbation (and thus

*C*) goes to zero.We note that the assumption of a constant lower bound of the minimal margin is not necessary in contrast to Jin (2002)
. Indeed this assumption would be highly problematic, because in an irregular dynamics arbitrarily small margins naturally occur with positive probability in every step. So, assuming some lower bound would exclude generic irregular dynamics. In contrast, the novel approach introduced above enabled us to show that the generic irregular dynamics is stable.

For infinitely large networks in a mean field approach, where one takes the limit of infinitely many neurons in the beginning (e.g. v.Vreeswijk and Sompolinsky, 1996
, 1998
; Brunel and Hakim, 1999
), our method is not applicable in a straightforward way. In this limit the average inter-spike interval goes to zero and a positive margin cannot be presumed. In our analysis we employ the fact that the minimal margin stays finite, i.e. a sufficiently small perturbation does not change the order of events. This assumption and therefore our results hold for arbitrarily large finite systems.

### Convergence to Periodic Orbits

Interestingly, the asymptotic stability together with the finite phase space imply that generic spike sequences converge to a periodic orbit. To show convergence analytically, we extend and explicate the ideas presented by Jin (2002) and Jahnke et al. (2008a)
. In the following we focus on a finite subsequence

*s** of a spike sequence generated by a given network. The number*E** of events in*s** is called the length*E** of*s**. As discussed towards the end of Section “Network Model”, the considered system is finite dimensional with dimension bounded from above by*N*+*ND*′. Thus, if*s** is sufficiently long it contains at least two disjoint subsequences*s*_{1}and*s*_{2}of length*E*, where the ordering of events is identical. The maximal*E*for which the existence of*s*_{1}and*s*_{2}is guaranteed, is given by the largest integer*E*satisfyingWhen increasing the observation length

*E** also the possible length*E*of the subsequences increases. Both the phases and the variables encoding the spikes in transit are bounded to a finite interval, ϕ_{i}(*t*) ∈ [−*T*+ ϕ_{Θ,i},ϕ_{Θ,i}] and σ_{ijk}(*t*) ∈ [0, max_{ij}{τ_{ij}}] at any given time*t*. Therefore, the maximal event-based distance between two trajectories is also bounded to a finite size,Thus, comparing sequences

*s*_{1}and*s*_{2}the initial event-based distance between their underlying trajectories at the beginning of*s*_{1}and at the beginning of*s*_{2}is bounded by Φ_{max}, i.e. By definition, the order of events in both sequences*s*_{1}and*s*_{2}is the same; therefore the distance between them shrinks according to Eq. 32. After*E*events the distance is bounded by (cf. Eq. 37)where we used the definition (39). We note that, if we increase our observation length

*E** and therefore the subsequence length*E*, the maximal possible distance between the trajectories underlying the sequences*s*_{1}and*s*_{2}decreases. If is sufficiently small, i.e. the distance between*s*_{1}and*s*_{2}after*E*events is smaller than the average margin, there is a high probability that also the order of events in the sequence following*s*_{1}and*s*_{2}are the same. Analogous to Eq. 41 the probability is given by:which goes to one when the observation length

*E** tends to infinity.This implies a periodic dynamics: Assume that

*s*_{1}occurs first in*s** at the*a*th event and after a certain event number*a*+*L*≤*E** −*E*the sequence*s*_{2}begins. We have seen that (with a certain probability) the ordering of events for (infinite) sequences starting at the*a*th and the (*a*+*L*)th event is the same for all future times. Therefore also the sequence starting at the (*a*+ 2*L*)th event is identical to the ones mentioned before, so the ordering of events is periodic. Together with the exponential convergence of equally ordered sequences this implies that also the spike timings converge towards a periodic orbit. For arbitrarily long observed sequences this happens with an arbitrary large probability that tends to one as*E** → ∞.In Figure 7
we show a typical example: The mean margin, decays algebraically on the transient and saturates after the periodic orbit is reached. Interestingly, the periodic attractor (shown in Figure 7
D) of the sparse random network resembles the “splay state” known from globally coupled networks (Strogatz and Mirollo, 1993
). In a splay state the firing pattern is characterized by equally spaced ISIs. It has been shown that it is possible to design networks, which exhibit more complex periodic spike patterns (Memmesheimer and Timme, 2006a
,b
). Indeed, in different parameter regimes we observe such complex periodic orbits, with a large periodicity, cf. the heterogeneous globally coupled network in Figure 8
, where the periodic orbit is reached after a small number of events compared to sparse networks. As shown previously (Timme et al., 2002
; Timme and Wolf, 2008
), highly irregular spiking activity may coexist with even the simplest (fully synchronous) periodic orbits that exhibits regular, maximally ordered activity.

**Figure 7. Convergence towards a periodic orbit in a random network.**(

*N*= 40, γ

_{i}≡ 1, I

_{i}≡ 3.0,

*V*

_{Θ,i}≡ 1.0, |Pre(

*i*)| ≡ 8, )

**(A)**Coupling matrix, each realized connection is indicated by a black square.

**(B)**Average minimal margin (averaged over 250 random initial conditions, cf. Figure 6 C) decays as a power-law (region

*A*) and saturates after about 10

^{7}events (region

*B*) when the periodic orbit is reached. Inset: Margin μ

^{(n)}(black) and minimal margin κ

^{(n)}(gray) for a trajectory started from one specific initial condition. The margin μ

^{(n)}fluctuates strongly on the transient and is comparatively large after the sequence becomes periodic; thus the minimal margin κ

^{(n)}does not decrease further for future events

*n*.

**(C,D)**: Snapshots of irregular spike sequences

**(C)**after

*n*≈ 15.000 events on the transient and

**(D)**after

*n*≈ 10

^{8}events on the periodic orbit.

**Figure 8. Long periodic orbit in a random network.**(

*N*= 50, γ

_{i}≡ 1,

*I*

_{i}≡ 1.2,

*V*

_{Θ,i}≡ 1.0, ) The coupling strengths are randomly independently drawn from the uniform distribution [−1, 0] and normalized afterwards The final spike pattern (shown after a transient of 9 × 10

^{6}events) repeats every 11012 events and is highly irregular.

Although the attractor is reached after a finite number of events, the transient becomes very important in systems with strong inhibition or large number of neurons. As formerly found in networks with excitatory coupling (Zumdieck et al., 2004
), and also in weakly diluted networks with purely inhibitory interactions (Zillmer et al., 2006
), the transient length grows rapidly with network size such that the dynamics is governed by the transient for large time scales. We study inhibitory random networks with an arbitrary network structure, typically far away from the weakly diluted topology. To perform numerical measurements of the transient length in dependence on the network size

*N*, we define the length of the transient,*t*by the number of events after which the order of events stays periodic. When increasing the network size_{r}*N*, we leave the sum of the external excitatiory current,*I*, and the internal inhibition, constant. Thus, on average each neuron receives a constant effective input independent of_{i}*N*and the mean firing rate of a single neuron 〈ν_{i}〉 is approximately conserved. Figure 9 shows the increase of transient lengths with network size for different sizes of internal inhibition and different scalings of the single connection strengths. We observe an exponential increase of the mean transient length with network size*N*independent of the scaling of the coupling strengths, ε_{ij}, and the number of incoming connections,*K*:= |Pre(_{i}*i*)|. This is qualitatively similar to the scaling of the transient lengths in weakly diluted networks (Zillmer et al., 2006 ). Assuming that the rapid increase continues for much larger networks, the transient will dominate the dynamics essentially forever in networks of biological relevant sizes. In this sense the transient becomes quasi-stationary (cf. also Zillmer et al., 2006 , 2009 ). If a larger network is in the balanced state (cf. Figure 1 ), the stable transients typically dominate the network dynamics.**Figure 9. Scaling of the transient lengths in sparse random networks.**

**(A)**The sum of the external excitation and the internal inhibition is fixed, such that the mean firing rate, 〈ν

_{i}〉

_{i}, of the single neurons stays approximately constant. The number of incoming connections to each neuron,

*K*:= |Pre(

_{i}*i*)|, is fixed to

*K*≡ 15 (green +),

_{i}*K*≡ 20 (black ×),

_{i}*K*≡ 25 (red *). The single connections are set to with

_{i}*e*

_{0}= −0.38. The solid lines show the best exponential fit. We observe an exponential increase of the trial averaged transient length 〈

*t*〉 with network size

_{r}*N*. The inset shows the dependence of the transient length on the in-degree for fixed network size (

*N*= 100).

**(B)**Networks with fixed fraction of connections, where

*K*≡ ⌊

_{i}*pN*⌋ and

*p*= 1/4. We scale the coupling strengths as such that and the external current

*I*is adjusted appropriately to fix the mean firing rate 〈ν

_{i}_{i}〉

_{i}. The mean transient lengths for different

*e*

_{0}are shown together with the best exponential fit (

*e*

_{0}= −0.34 (green +),

*e*

_{0}= −0.38 (black ×),

*e*

_{0}= −0.42 (red *)). Again, we observe an exponential increase of the transient lengths with

*N*. The inset shows the increase of the transient length with inhibitory coupling strength for fixed network size (

*N*= 100).

**(C)**Networks with fixed fraction of connections, where

*K*≡ ⌊

_{i}*pN*⌋ and

*p*= 1/4. The external currents

*I*≡ 3.0 are fixed and the internal coupling is normalized to (green +), −1.5 (blue ×), −1.75 (black *), −1.9 (red □), such that the mean firing rate, 〈ν

_{i}_{i}〉

_{i}, of the single neurons is approximately the same for each curve. The number of incoming connections per neuron,

*K*, increases linearly with network size, so the coupling strengths are scaled as ε

_{i}_{ij}∝ 1/

*K*. Here we also observe an exponential increase of 〈

_{i}*t*〉 with network size. The inset shows the fast increase of the transient length with inhibitory coupling strength for fixed network size and external current (

_{r}*N*= 100,

*I*≡ 3.0).

_{i}**(D)**The transient length is broadly distributed as shown in the histogram for 2500 trials started from random initial conditions, where the initial phases were randomly independently drawn from the uniform distribution on [0, 1] (

*N*= 100,

*I*≡ 3.0, ).

_{i}### Robustness of Stability and Smooth Transition to Chaos

In the following we will check the robustness of our results. The considerations above hold for networks with inhibitory interactions without temporal extent. We investigate the influence of excitatory interaction and pulses with a finite duration. For small deviation from the networks considered above the stability properties are similar, for large fractions of excitatory connections and larger temporal extent we observe a transition to a chaotic regime. For temporally extended couplings we assume that after a neuron is reset all previous input is lost. Therefore the state of a neuron is specified by its last spike time and all spikes it has received afterwards. The phase representation is thus not meaningful anymore and we track two trajectories by comparing the differences in the last spiking times of the neurons and in the spike arrival times since these last spikings. To keep the section consistent, we adopt this view when studying the Lyapunov exponents of the excitatory dynamics.

#### Excitatory interactions

We have shown that in networks with purely inhibitory interactions the dynamics is typically stable. If the connection from neuron

*j*to neuron*i*is excitatory, the phase shift of the postsynaptic neuron*i*after receiving a spike from neuron*j*as the (*n*+ 1)th event may exceed its shift before and the shift of the received spike. Figures 10 A,B gives an illustration: A spike is simultaneously received in the perturbed and in the unperturbed dynamics The phase shift before and after the application of the transfer function is shown for (A) an inhibitory input and (B) an excitatory input. For inhibitory input the phase shift is reduced, this leads to the stable dynamics as described in the sections above. For excitatory input, the phase shift increases when the spike is received.**Figure 10. Destabilizing effect of excitation.**

**(A)**Simultaneous receiving of an inhibitory input decreases the phase shift.

**(B)**Simultaneous receiving of an excitatory input increases the phase shift.

**(C)**The largest Lyapunov exponent is measured for an increasing fraction of excitatory connections starting with the network of Figures 1 A–C and ending with the network of Figures 1 D–F. Here we keep the firing rate constant. Inset shows exemplarily the convergence of a finite time Lyapunov exponent with

*n*.

**(D)**Average firing rate 〈ν

_{i}〉

_{i}versus number of excitatory connections in the network for different external currents

*I*(black crosses, values belonging to the same

_{i}*I*are connected by a dashed line). The neurons’ firing rate stays almost constant if we reduce the external current linearly with the number of excitatory connections (blue crosses). The inset displays the current strength employed to maintain firing rate of 〈ν

_{i}_{i}〉

_{i}≈ 0.23. (Further details see text).

Indeed, since the inverse of

*U*(ϕ) is monotonically increasing with ϕ, we find for a given ε > 0_{i}In contrast to Eq. 22, the derivative of the transfer function is bounded from below by

According to Eq. 14, this can lead to an increase of (in particular extremal) perturbations and to a destabilization of the trajectory. The upper bound of

*c*_{max}< 1 (cf. Eq. 23) does not hold anymore. However, in a network with a small fraction of excitatory connections, the trajectory is still stable. At an interaction the perturbation may increase, but the stabilizing effect of inhibitory inputs dominates the dynamics. We study the transition from the stable regime to chaotic dynamics (a discussion of the chaotic dynamics in networks with purely excitatory interactions can be found in Zumdieck et al., 2004 ). When increasing the number of excitatory couplings, we increase the mean effective input current to the neurons. Thus we additionally decrease the external input*I*to keep the network rate ν constant. Indeed, in good approximation, the current has to be decreased linearly with_{i}*N*, the number of excitatory connections,_{E}*I*≡_{i}*I*−*kN*where_{E}*i*is the original input current.To quantify the transition we estimate the largest Lyapunov exponent of the system: At the

*n*th event, we denote*n*−*W*(*n*) the earliest event which still influences the future dynamics of the system explicitly. We apply an initial perturbation of size ||Δ_{0}||_{1}to the event times*t*_{0},*t*_{−1},…,*t*_{−W(0)}, where Δ_{n}= (Δ*t*, Δ_{n}*t*_{n − 1},…,Δ*t*_{n − W(n)}) is the perturbation vector at the*n*th event time and Δ*t*is the perturbation of_{i}*t*. Here ||.||_{i}_{1}is the 1-norm, ||*x*||_{1 }:= Σ_{i}|*x*|. We evolve the system and rescale the perturbation vector Δ_{i}_{n}by*a*after each event, such that the rescaled perturbation vector is of the same size as the initial perturbation,_{n}The largest Lyapunov exponent, λ

_{max}, is then given byWe observe a transition from a stable to a chaotic regime, characterized by a positive Lyapunov exponent. For small fraction of excitatory neurons the dynamics is typically stable, the effect of the inhibitory pulses dominates the dynamics and, on average, a perturbation do not grow over time. With increasing

*N*the Lyapunov exponent increases until the dynamic becomes chaotic. Of course, in our simulations we can only study finite time Lyapunov exponents with very large_{E}*n*and estimate the value to which they converge. The chaotic dynamics may thus be transient. However, it dominates the dynamics at least over very long times.Estimating the maximal Lyapunov exponent in networks including excitatory interactions can be difficult (Hansel et al., 1998
; Zumdieck et al., 2004
; Brette et al., 2007
; Cessac and Vieville, 2008
; Kirst and Timme, 2009
). Suprathreshold excitation, together with the infinitely fast response of neurons receiving a spike and the sharp threshold may induce synchronous events. Thus even an infinitesimal small perturbation may change the order of spikes. Nonetheless generically the perturbation will stay infinitesimal small, in particular for a small fraction of excitatory connections, such that we estimate the largest Lyapunov exponent in the following way: We evaluate at each time step the resulting temporal perturbation on the actual event as a result of earlier perturbation under the assumption that the order of spikes stays the same. This gives us the new perturbation vector Δ

_{n}, which is rescaled according to (48). For long times, λ(*n*) then will give an estimate of the largest Lyapunov exponent and describe the generic behavior of the trajectory under the influence of sufficiently small perturbations. However, we cannot exclude the occurrence of macroscopic perturbations in general.Figures 10
C,D show some numerical results: In Figure 10
C the largest Lyapunov exponent is measured for an increasing fraction of excitatory connections starting with the network of Figures 1
A–C and ending with the network of Figures 1
D–F. The number of excitatory connections is increased by successively choosing one incoming inhibitory connection per neuron to be excitatory. The external current,

*I*, is reduced linearly to keep the network rate unchanged according to_{i}*I*≡_{i}*I*−*kN*, where_{E}*I*= 4.0 is the initial external current,*k*≈ 1.3 10^{−4}and*N*is the number of excitatory connections. For a large fraction of excitatory couplings we observe a transition to an unstable, chaotic regime. The inset demonstrate the convergence of (exemplarily shown for_{E}*N*= 10000). For a constant external current the rate increases with increasing fraction of excitatory connections (black crosses, for_{E}*I*∈ {2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 4.5}). The neurons’ firing rate stays almost constant, if we reduce the external current linearly with the number of excitatory connections (blue crosses). We determined the value of_{i}*I*at the intersection point of the 〈ν_{i}_{i}〉_{i}vs.*N*curves with the desired frequency by linear interpolation. The values (_{E}*I*,_{i}*N*) that give rise to the desired frequency lie in good approximation on a straight line with slope −_{E}*k*≈ −1.3 10^{−4}(cf. Inset to (D)).The result is particularly remarkable since in mean-field descriptions of balanced networks, as long as the mean input to each cell is the same, the regime where

*N*= 0 is comparable to the regime where_{E}*N*> 0 with appropriately reduced external excitatory current_{E}*I*._{i}#### Temporally extended interactions

Up to now, we considered δ-coupling, where the response to an action potential is instantaneous. However, in biological neuronal systems the postsynaptic current has finite temporal extent. In the following, we investigate the influence of such temporally extended interactions. The analysis gets more complicated, because neurons are permanently influenced by incoming signals. As mentioned above, in our model we assume that the neuron looses the information about previously received spikes when it reaches the threshold and is reset.

We modify Eq. 1 by introducing a temporally extended interaction kernel

*g*(*t*), such that the evolution of the membrane potential is given byIn the following analysis we consider single exponential couplings,

*g*(*t*) = Θ (*t*)·β*e*^{−βt}with time scales γ^{−1}> β^{−1}> 0, the time constant of the postsynaptic current is shorter than the membrane time constant. As an exemplary neuron model we study the leaky integrate-and-fire neuron,*f*(_{i}*V*) = −γ_{i}*V*(_{i}*t*) +*I*, but the analysis can easily be extended to more complex neuron models and interaction kernels._{i}Numerical simulations show that the stability of the dynamics is robust against introduction of synaptic currents with small temporal extent, but on increase of temporal extension a transition to chaos occurs. In Figure 11
A, the largest Lyapunov exponent, λ

_{max}, in a random network is estimated in dependence of the decay time constant β^{−1}of the synaptic current. For small time constant β^{−1}, the dynamics behaves similar to the dynamics with δ-pulse interactions, in particular it is stable, the largest Lyapunov exponent is negative. For increasing β^{−1}the temporal extension becomes more and more influential and there is a transition to an unstable, chaotic regime with positive largest Lyapunov exponent.**Figure 11. Robustness and transition to chaos with increasing temporal extent of the postsynaptic current.**

**(A)**Largest Lyapunov exponent λ

_{max}in a random network (

*N*= 50, |Pre(

*i*)| ≡ 10, γ = 1,

*I*≡ 4.0, ) versus time constant β

_{i}^{−1}of the synaptic current.

**(B)**Stable periodic pattern in a network with temporally extended interactions.

**(C)**Time compressed, thus unstable pattern in a network with temporally extended interactions. (For details see text).

We now study the linear stability properties analytically. We denote the last spiking time of neuron

*i*before*t*by_{n}at

*t*=*t*_{0}(*n*,*i*) the potential of neuron*i*was reset to zero. The solution of Eq. 50 together with the initial condition*V*(_{i}*t*_{0}(*n*,*i*)) = 0 is then given between the*n*th and (*n*+ 1)th network event bywhere is the reception time of the spike sent at by neuron

*j*at neuron*i*. The sum in Eq. 52 takes into account all spikes which are received by neuron*i*between*t*_{0}(*n*,*i*) ≤*t*≤*t*and therefore influence the potential_{n}*V*_{i,n}(*t*). In the limit of very short temporal extension of the postsynaptic current, β → ∞, Eq. 52 becomes a solution of Eq. 1. After the*n*th event, neuron*i*would reach the threshold at some time*t*′ under the assumption that there are no further inputs after*t*. According to Eq. 52,_{n}*t*′ is implicitly given bywhere

*A*is the vector of the original event times*t*,…,_{n}*t*_{n − W},where we introduced

*W*= max_{n}{*W*(*n*)}. We now estimate the effect of a small perturbation Δ*t*,…,Δ_{n}*t*_{n − W}of the event times*t*,…,_{n}*t*_{n − W}on the hypothetical event time*t*′. By Eq. 53, the Jacobian of*t*′,*Dt*′, with respect to former spike times,*t*,…,_{n}*t*_{n − W}, is given asThe linearized estimation of the displacement Δ

*t*′of*t*′ is then given byThe special structure of

*V*_{i,n}(*t*) (cf. Eq. 52), more precisely the fact that*V*_{i,n}(*t*) depends on*t*via*t*−*t*for_{k}*k*∈ {*n*−*W*,…,*n*}, yields the identityUnder the condition

we can combine Eqs. 56 and 57 and find bounds for the displacement

Condition (58) implies that if neuron

*i*sends or receives a spike earlier, also the threshold is crossed earlier. This always holds for δ-couplings, for interactions with temporal extend it restricts the class of patterns as we show below. Equation 59 is an analog to Eq. 25, sufficiently small perturbations stay bounded by the initial ones for finite times. This directly implies Lyapunov stability for periodic orbits. For general irregular dynamics and to prove asymptotic stability, the propagation of pulses through the network has to be studied as for the nonlinear stability analysis in the main part.We now want to specify a class of periodic patterns which are stable in a network with temporally extended synaptic currents. The influence of various events on

*V*_{i,n}(*A*,*t*′) is as follows: For an influential spike receiving*t*, Eq. 52 yields_{k}where ε* < 0 is the coupling strength from the sending neuron. For the last spike sending of neuron

*i*,*t*=_{k}*t*_{0}(*n*,*i*),For any other event

*t*a displacement of_{k}*t*has no influence on_{k}*V*_{i,n}(*t*), hereTherefore condition (58) reduces to a condition on the left and right hand side of Eq. 60 and can be reformulated as

where

*t*are the spikes arrival times at neuron_{k}*i*since the last reset*t*_{0}(*n*,*i*). This means that the class of patterns where each neuron*i*does not cross the threshold for a time period*T*after receiving a spike are stable. For β → ∞ the system tends to the δ – pulse coupled system and indeed_{d}*T*vanishes, such that any non-degenerated orbit is stable. However, for temporally extended interactions unstable periodic orbits exists and also chaotic dynamics is possible (cf. Figure 11 A)._{d}To illustrate our analytical findings, in Figures 11
B,C, we used a generalization of a recently introduced method (Memmesheimer and Timme, 2006a
,b
) to design two networks realizing predefined spike patterns in a network of five neurons (

*V*_{Θ,i}≡ 1.0,*I*≡ 2.4, β = 8, τ := τ_{i}_{ij}≡ 0.125) with temporally extended couplings. Both patterns are the same, but with different ISIs. In (B) all spikes are separated by Δ*T*=*T*+ τ, which ensures that a neuron never spikes within a time period_{d}*T*after receiving a spike; in (C) we choose the ISIs smaller Δ_{d}*T*= (*T*+ τ)/2. The lower panels illustrates the stability properties: The spike times of the different neurons are plotted relative to the spike time of neuron 1 in vertical direction. The horizontal direction is simulated time, different colors indicate spike times of the five different neurons. At certain points in time (blue arrows) the network dynamics is perturbed. The dynamics in (B) is stable: After perturbations of size ≈0.2 (maximum norm), the dynamics converge towards the periodic orbit. (C) The dynamics is unstable: a perturbation of size ≈10_{d}^{−12}leads to a divergence from the unstable periodic orbit.Irregular spiking activity that robustly arises in balanced state models constitutes a generic feature of cortical dynamics. Here, for a class of models including e.g. the common leaky integrate-and-fire neuron, we have shown that generic trajectories which give rise to the irregular balanced state can be exponentially stable. In particular, they are stable in purely inhibitory strongly connected networks of neurons with delayed couplings and with infinitesimal synaptic time course. We showed that the dynamics even converges to a periodic orbit. However, the length of the irregular yet stable transient grows rapidly with network size such that for larger networks, in particular for biological relevant sizes, transients may dominate the dynamics on all relevant time scales.

Furthermore, we found that the phenomenon of stable yet irregular dynamics is robust against introducing some excitatory interactions or against increasing the synaptic time scales from zero. If the synaptic responses become too slow or excitatory interactions too many, we revealed a smooth transition from stable irregular dynamics to chaotic, equally irregular dynamics. We emphasize that we kept the network rate during this transition (and thus keep the balance) and that the mean field descriptions (v.Vreeswijk and Sompolinsky, 1998
; Brunel, 2000
) of networks in both regimes are identical when the parameters are suitably chosen. Thus, highly irregular spiking dynamics occurs independent of the stability properties of the network.

Earlier studies on balanced neural activity considered a priori the limit of infinitely many neurons in sparse networks (v.Vreeswijk and Sompolinsky, 1996
, 1998
; Amit and Brunel, 1997
; Brunel and Hakim, 1999
; Brunel, 2000
). In this mean field limit the collective dynamics is well understood. In particular in infinitely large networks of binary neurons with balanced excitatory and inhibitory interactions the dynamics are chaotic (v.Vreeswijk and Sompolinsky, 1996
, 1998
). Further studies of finite networks found stable dynamics in weakly diluted networks of inhibitory coupled neurons (Zillmer et al., 2006
), as well as in globally coupled networks with dominating inhibition (Jin, 2002
). Recent analytical evidence confirmed the existence of stable dynamics in inhibitory coupled networks of integrate-and-fire neurons with a more complex structure (Jahnke et al., 2008a
). As the inter-event times that underly our analysis shrink inversely proportional to the network size (at a given individual neuron-spiking-rate), the methods applied here, however, are not applicable in a straightforward way in associated mean field models. Thus, one cannot make strict statements about stability in the limit of infinitely many neurons. Nevertheless, as shown above, generic transients and periodic trajectories in arbitrarily large inhibitorily coupled networks are stable.

Taken together, the results show that the microscopic dynamics in the considered purely inhibitorily coupled networks differs substantially from the dynamics of networks that explicitly include excitatory couplings. Whereas the latter may be chaotic, the former are generically stable – despite both showing the same irregularity features. In particular, chaotic as well as stable dynamics are equally well capable of generating highly irregular spiking activity. The smoothness of the transition to chaos, without essential change of the irregularity (e.g. of the large coefficients of variation) further suggests that chaos is not the main dynamical origin of the high irregularity. We thus suggest that a mechanism different from chaos contributes to the irregularity of cortical firing patterns in a substantial way.

Nevertheless, chaos as well as stochastic network properties such as unreliable synapses, may support the robust occurrence of irregular activity in cortical networks and also modify its computational features. It is thus an important future task to investigate which anatomical and dynamical features of cortical networks are indeed of crucial relevance for their spiking activity and their functions.

In the following we will derive the bounds Eqs. 26–28 stated in the Section “Asymptotic Stability”. We will track the propagation of the perturbation of one specific neuron

*l*_{0}through the entire network.All neurons spike at least once in a sufficiently large but finite time interval

*T*. Moreover, after τ_{max}= max_{i,j}(τ_{ij}) all spikes in transit have certainly arrived at the postsynaptic neurons. We label the maximal number of events possible in the time interval [*t*,*t*+ max{*T*,τ_{max}}] by*M*. For purely inhibitory networks,*M*< ∞ due to the bounded neural spike rate. We denote the set of postsynaptic neurons of*l*_{0}bythus a neuron

*l*∈ Post_{i}^{i}(*l*_{0}) is connected to*l*_{0}by a directed path of length*i*(cf. Figure 12 ). Further, we define Post^{0}(*l*_{0}) := {*l*_{0}}.We estimate the bounds of the perturbation following one specific path from a neuron

*l*_{0}to a neuron*l*. In a strongly connected network,*l*∈ Post^{j}(*l*_{0}) for some*j*≤*N*− 1, so there is a directed path between*l*_{0}and*l*=*l*via neurons_{j}*l*_{1},…,*l*_{j − 1}in the network. As the consideration holds for an arbitrary path, the result is an universal bound of the perturbation. Initially, at*n*= 0, the neurons are perturbed by The first spiking of neuron*l*_{0}after*n*= 0 is labeled by*s*_{0}≤*M*. After a delay time this spike is received by the postsynaptic neuron*l*_{1}∈ Post(*l*_{0}), we call the event*r*_{1}≤ 2*M*. After at most*M*further events, at*s*_{1}≤ 3*M*, the neuron*l*_{1}emits a spike. In general, we recursively define*s*as the first spiking event of neuron_{i}*l*∈ Post_{i}^{i}(*l*_{0}) after*r*and_{i}*r*as the event when the spike generated by_{i}*l*_{i − 1}∈ Post^{i − 1}(*l*_{0}) at s_{i − 1}is received (cf. Figure 12 ). Due to the definition of*M*, the relations*s*≤ (2_{i}*i*+1)*M*and*r*≤ 2_{i}*iM*hold.**Figure 12. Tracking the propagation of a pulse.**Here neuron

*l*

_{0}= 1 (black) is fixed as the initial neuron. The sets of postsynaptic neurons are: Post(1) = {2} (blue), Post

^{2}(1) = {1, 3} (green), Post

^{3}(1) = {1, 2} (red). Following one specific path through the network, we label the first spike event of neuron

*l*

_{i}∈ Post

^{i}(

*l*

_{0}) after receiving an input from

*l*

_{i − 1}∈ Post

^{i − 1}(

*l*

_{0})

*s*

_{i}. The event when the generated spike is received by

*l*

_{i + 1}is labeled

*r*

_{i + 1}. In a strongly connected network of size

*N*the union contains all the neurons, because any two of them are connected by a directed path of maximal length

*N*− 1.

First, we prove by induction that the perturbation of the neuron

*l*before sending of a spike at_{i},*s*is bounded from above by_{i},(1) Initially neuron

*l*_{0}is perturbed by Before*l*_{0}generates a spike it receives at most*s*_{0}− 1 inputs. According to Eq. 14, if the neuron*l*_{0}receives an input its perturbation may increase. To find an upper bound, we assume that at every event 0 <*n*<*s*_{0}neuron*l*_{0}receives an input with the maximal initial perturbation, and a minimal averaging constant*c*_{min}, which moves the average into the direction of the maximal possible perturbation. Repeated application of Eq. (14) yieldswhich is the inductive statement (65) for

*i*= 0.(2) We assume that the statement (65) holds for which is neuron

*l*’s perturbation as inherited by the spike sent at_{i}*s*(cf. Eq. 10). After at most_{i}*M*events the spike is received by the postsynaptic neuron*l*_{i + 1}at event*r*_{i + 1}. In our worst- (or worse than worst-) case estimation, we assume that neuron*l*_{i + 1}is maximally perturbed before it receives the spike, and that the interaction factor is maximal,*c*_{max}, such that again the average is moved into the direction of the maximal perturbation. Therefore the perturbation after the interaction is bounded byBefore

*s*_{i + 1}>*r*_{i + 1}>*s*, neuron_{i}*l*_{i + 1}receives at most (*s*_{i + 1}− 1 −*r*_{i + 1}) inputs. Analogously to Eq. 66, we assume that with each event*l*_{i + 1}receives a spike which is maximally perturbed (with ), and the averaging constant is minimal,*c*_{min}. This yieldsWe replace by in Eq. 68, thereby increasing the right-hand side, because

*s*− 1 −_{i}*r*_{i + 1}< 0. This directly yields the induction statement forBased on Eq. 65, we now derive an upper bound of the perturbation of all neurons after event

*s*_{N − 1}. After this event, every neuron has sent at least one spike which is influenced by the initial perturbation of neuron*l*_{0}, because in a strongly connected network the union contains all neurons of the network (cf. Figure 12 ). After the*s*th event, neuron_{i}*l*can still receive spikes. Before the_{i}*s*_{N − 1}th event, taken as reference, it receives in the worst case scenario (*s*_{N − 1}−*s*) inputs with maximal initial perturbation and minimal averaging factor_{i}*c*_{min}. Using Eq. 65 we repeatedly apply Eq. 14 (*s*_{N − 1}−*s*) times which leads to_{i}The right-hand side increases with

*i*, therefore the perturbation of an arbitrary neuron*j*∈ {1,…,*N*} after*s*_{N − 1}events is bounded from above byAt the

*s*_{N − 1}th event there can be*D*′ spikes per neuron in transit which are, in the worst-case scenario, assumed to have the maximal perturbation. Due to their arrival after the*s*_{N − 1}th event, the perturbations of neurons can still increase. However, after*s*_{N − 1}+*M*events all spikes generated before the*s*_{N − 1}th event have arrived at the corresponding postsynaptic neurons. Taking into account the arrival of these spikes using Eqs. 14 and 70, we find an upper bound for the perturbation after*s*_{N − 1}+*M*events,Due to the fact that generated spikes inherit a perturbation present in the phases at spike sending time, the bound (71) holds also for the perturbation of spikes generated after the

*s*_{N − 1}th and before the (*s*_{N − 1}+*M*)th event, because the bound (71) limits the maximal perturbation for all neurons between the*s*_{N − 1}th and the (*s*_{N − 1}+*M*)th event.We conclude that the perturbations of the neurons and the spikes in transit after

*K*:= 2*NM*≥*s*_{N − 1}+*M*events are bounded byTherefore, we find an upper bound for the maximal perturbation after

*K*events,with

Similarly, we find a lower bound for the minimal perturbation after

*K*eventsby an estimation analogous to the one above, where only has to be replaced and the relation “≤” has to be replaced by “≥”. We note, that we did not have to specify the perturbation to derive this result.

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.

We thank Fred Wolf for helpful comments, David Hansel for pointing out the article (Zillmer et al., 2009
) and the Federal Ministry of Education and Research (BMBF) Germany for partial support under Grant No. 01GQ0430. RMM thanks the Sloan Swartz Foundation for partial support.

Brette, R., Rudolph, M., Carnevale, T., Hines, M., Beeman, D., Bower, J. M., Diesmann, M., Morrison, A., Goodman, P. H., Harris, F. C. Jr, Zirpe, M., Natschlager, T., Pecevski, D., Ermentrout, B., Djurfeldt, M., Lansner, A., Rochel, O., Vieville, T., Muller, E., Davison, A. P., El Boustani, S., and Destexhe, A. (2007). Simulation of networks of spiking neurons: a review of tools and strategies.

*J. Comput. Neurosci.*23, 349–398.