How Chaotic is the Balanced State?

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.


INTRODUCTION
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 fi ring rate is thus considered a basic "ground state". It is characterized by individual neurons that display largely fl uctuating membrane potentials and highly variable inter-spike-intervals (ISIs) as well as by low correlations between the neurons (v.Vreeswijk andSompolinsky, 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 fl uctuations, therefore excluding the emergence of irregular dynamics. So the fi nding 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 fl uctuations once in a while are suffi ciently 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 How chaotic is the balanced state?

Sven Jahnke 1,2 *, Raoul-Martin Memmesheimer 3 and Marc Timme 1,2
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 fi nite 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 suffi ciently 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 signifi cantly contributes to irregular activity in cortical circuits. 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.
Some analytically accessible aspects of the stable irregular dynamics have been briefl y 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).

NETWORK MODEL
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 i specifi es the internal dynamics, ε 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 t jk s denotes the time when neuron j sends the kth 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 t jk s . When neuron j reaches the threshold V Θ,j of the potential, i.e. V j (t − ) = V Θ,j , it generates spikes at t t jk =: s for some k and is reset, V t j jk ( ) . s = 0 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 j satisfy ( ) 0 for all j and V j ≤ V Θ,j , such that in isolation each neuron j exhibits oscillatory dynamics with a period T j free . The network dynamics can equivalently be described by a phase- at all non-event times (Mirollo and Strogatz, 1990). When the phase threshold is reached, φ φ 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 fi nite constant T > 0, such that in every time interval [t, t + T ), t ∈ R, every neuron fi res at least once. Further, we assume that the dynamics is suffi ciently irregular such that two events occur at the same time with zero probability.
Due to the delay, the state space is formally infi nite dimensional. However, it becomes fi nite dimensional after some fi nite 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 φ = Θ,i i T free to generate a spike after the last reset. Consequently, at most spikes per neuron are in transit and the state space stays fi nite, 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 t jn s for some n, σ ijk (t) is set to σ ( ) = ijk jn t s 0. The spike index k equals the number of spikes already in transit plus one. It thus depends on the actual network state at time t jn s . Between two events σ ijk (t) increases linearly with slope one, when σ ( + τ ) = τ ijk jn ij ij t s 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 σ ( + τ ) ij jn ij t 1 s specifi es the spike sent by neuron j which arrives next at the postsynaptic neuron i (cf. Figures 3A,B for illustration).

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 nth event in the entire network is denoted by t n in the fi rst sequence, and by t n in the second one. Analogously, we denote the phases of a neuron i at a given time t by φ i (t) and φ ( ) i t and the spikes in transit at time t by σ ijk (t) and ijk t σ ( ) in the different sequences. Let denote the difference between the phase difference, δ := φ − , after the nth and before the (n + 1)th event (cf. Figure 3C). Similarly : labels the shift of the kth spike sent by neuron j and not yet arrived at neuron i after the nth and before the (n + 1)th event. Between two consecutive events, both φ φ i i t t ( ), ( ) and σ σ ijk ijk t t ( ), ( ) 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 suffi cient 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 fi rst case, the shifts Δ i n ( ) and Δσ ijk n ( ) 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, Additionally, new spikes are generated σ ij * k * (t n + 1 ) = 0 and ij k n t * * σ ( ) + = 1 0 where k * = k * (i, j * , n + 1) is the appropriate spike number, cf.

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 Δ * i n ( ) +1 can be computed as: where φ φ * are the phases "just before" spike reception. Using the defi nitions (6) we fi nd the identity

FIGURE 3 | Phase dynamics of neuron j and defi nition of state variables. (A)
At t = t n and 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 specifi c postsynaptic neuron i, described by σ ij1 (t) (black) and σ ij2 (t) (red). At t = t n , σ ij1 (t n ) is set to 0; here k = 1 because there is no spike in transit at t t n = − . When neuron j spikes again at t = t n + 1 , σ ij2 (t n + 1 ) is set to 0. At t = t n + τ ij the spike emitted at t = t n arrives at the postsynaptic neuron 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) Defi nition of the phase shifts. The phase curve φ i (t) of neuron i (blue) before and after the reception of a spike at t = t n is shown is the difference of neuron i's phases in the unperturbed and the perturbed dynamics taken at corresponding event times t n and t n . δ = t t t n n n ( ) − denotes the difference of event times t n and t n , i.e. the temporal offset between both sequences.
is some phase shift of neuron i with the temporal offset taken into account.

Frontiers in Computational
Applying the mean value theorem in Eq. 11 and the relation The shifts of traveling spikes stay unchanged on spike reception, cf. Eq. 9, holds, 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 coeffi cients c i j n * * ( ) +1 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 fi res 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 i also U i −1 is strictly monotonic increasing and ε < 0. The strict concavity of U i (φ) implies for any fi nite φ (cf. Figure 4 for illustration). The derivative ∂ 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 nth event, are bounded by the initial shifts for all future events, and (25) 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 n . 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.

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 specifi c 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 The transfer function H i ε φ ( ) ( ) for a leaky integrate-and-fi re neuron. For ε = 0 (black), H i 0 φ φ ( ) ( ) = 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 fi nd H i constant temporal offset. Further details and the derivation of the following Eqs. 26-28 is provided in the Appendix. We fi nd the upper bound for the maximal perturbation after K := 2NM events and analogously the lower bound for the minimal perturbation. The averaging factor c * is determined by the network parameters and bounded to A bound for the difference of the maximal and minimal perturbation after K events is therefore given by The 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)  0 Thus all sequences considered are asymptotically stable for all strongly connected networks and perturbations decay exponentially fast with at least where ⎣·⎦ is the fl oor 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.

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 suprathreshold 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 fi nite 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 infi nite time. Formally, after time t n denote the kth potential future event time (of the original trajectory) that would arise if there were no future interactions by θ n,k , k ∈ N, and the temporal margin by µ (n) := θ n,2 − θ n,1 . A suffi ciently small perturbation, satisfying Δ Δ <μ , max ( ) min ( ) ( ) n n n − 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 n ) consisting of a fi nite number of P events are stable because there is a minimal margin for 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 fi rst 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 specifi es the mean fi ring 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 defi nition of events. The distribution function of margins is given by Therefore, 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 that that depends only on the event rate and is independent of the specifi c network parameters. Numerical simulations show excellent agreement with this algebraic decay (36); a typical example is shown in Figure 6. This already strongly suggests that a suffi ciently small perturbation stays smaller than the minimal margin for all times. However, in each step, the exponential distribution of margins has a fi nite 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 P n n n n (∃ ∈ : μ ≤ Δ Δ , 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 Δ −Δ max ( ) min ( ) n n for infi nite time. We show that it converges to one when the size of the perturbation goes to zero and thus prove that suffi ciently 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 n K n K ⎢ ⎣ ⎥ ⎦ ≥ −1 and Eq. 28 leads to where we introduced In particular, C → 0 if the initial perturbation goes to zero, i.e. Δ −Δ → , since the margins are independent. Using Eq. 37 and P(µ (n) > µ) = exp(−νµ) yields  Algebraic decay of the average minimal margin, κ ( ) n (green dashed line, averaged over 250 random initial conditions) and its analytical prediction (no free fi t parameter; black solid line). Additionally we show the minimal margin κ (n) for three exemplary initial conditions (gray lines), including that of Figure 5A. which 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 infi nitely large networks in a mean fi eld approach, where one takes the limit of infi nitely many neurons in the beginning (e.g. v. Vreeswijk andSompolinsky, 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 fi nite, i.e. a suffi ciently small perturbation does not change the order of events. This assumption and therefore our results hold for arbitrarily large fi nite systems.

CONVERGENCE TO PERIODIC ORBITS
Interestingly, the asymptotic stability together with the fi nite 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 fi nite 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 fi nite dimensional with dimension bounded from above by N + ND′. Thus, if s * is suffi ciently 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 satisfying (42) When 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 fi nite 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 fi nite 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 defi nition, 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 defi nition (39). We note that, if we increase our observation length E * and therefore the subsequence length E, the maximal possible distance Φ * max between the trajectories underlying the sequences s 1 and s 2 decreases. If Φ * max is suffi ciently 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 infi nity. This implies a periodic dynamics: Assume that s 1 occurs fi rst in s * at the ath 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 (infi nite) sequences starting at the ath and the (a + L)th event is the same for all future times. Therefore also the sequence starting at the (a + 2L)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, κ ( ) , n decays algebraically on the transient and saturates after the periodic orbit is reached. Interestingly, the periodic attractor (shown in Figure 7D) of the sparse random network resembles the "splay state" known from globally coupled networks (Strogatz and Mirollo, 1993). In a splay state the fi ring 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.
Although the attractor is reached after a fi nite 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 defi ne the length of the transient, t r by the number of events after which the order of events stays periodic. When increasing the network size N, we leave the sum ˆ: of the external excitatiory current, I i , and the internal inhibition, ∑ ε j N ij =1 , constant. Thus, on average each neuron receives a constant effective input independent of N and the mean fi ring 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 i := ⏐Pre(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 quasistationary (cf. also Zillmer et al., 2006Zillmer et al., , 2009. If a larger network is in the balanced state (cf. Figure 1), the stable transients typically dominate the network dynamics.

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 infl uence of excitatory interaction and pulses with a fi nite 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 specifi ed 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 Δ i n ( ) +1 of the postsynaptic neuron i after receiving a spike from neuron j as the (n + 1)th event may exceed its shift Δ i n ( ) before and the shift Δσ ij n 1 ( ) of the received spike. Figures 10A,B gives an illustration: A spike is simultaneously received in the perturbed and in the unperturbed dynamics (i.e. Δσ Δ i n ( ) is reduced, this leads to the stable dynamics as described in the sections above. For excitatory input, the phase shift Δ i n ( ) increases when the spike is received.
Indeed, since the inverse of U i (φ) is monotonically increasing with φ, we fi nd for a given ε > 0 In contrast to Eq. 22, the derivative of the transfer function 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 ij n ( ) , 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 i to keep the network rate ν constant. Indeed, in good approximation, the current has to be decreased linearly with N E , the number of excitatory connections, I i ≡ I − kN E where I is the original input current.
To quantify the transition we estimate the largest Lyapunov exponent of the system: At the nth event, we denote n − W(n) the earliest event which still infl uences 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 nth event time and Δt i is the perturbation of t i . Here ||.|| 1 is the 1-norm, ||x|| 1 := Σ i |x i |. We evolve the system and rescale the perturbation vector Δ n by a n after each event, such that the rescaled perturbation vector is of the same size as the initial perturbation, The largest Lyapunov exponent, λ max , is then given by We 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 E the Lyapunov exponent increases until the dynamic becomes chaotic. Of course, in our simulations we can only study fi nite time Lyapunov exponents with very large 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 diffi cult (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 infi nitely fast response of neurons receiving a spike and the sharp threshold may induce synchronous events. Thus even an infi nitesimal small perturbation may change the order of spikes. Nonetheless generically the perturbation will stay infi nitesimal 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 infl uence of suffi ciently small perturbations. However, we cannot exclude the occurrence of macroscopic perturbations in general. The result is particularly remarkable since in mean-fi eld descriptions of balanced networks, as long as the mean input to each cell is the same, the regime where N E = 0 is comparable to the regime where N E > 0 with appropriately reduced external excitatory current 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 fi nite temporal extent. In the following, we investigate the infl uence of such temporally extended interactions. The analysis gets more complicated, because neurons are permanently infl uenced 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 by In 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-fi re neuron, f i (V i ) = −γV i (t) + I i , but the analysis can easily be extended to more complex neuron models and interaction kernels.
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 11A, 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 infl uential and there is a transition to an unstable, chaotic regime with positive largest Lyapunov exponent. We now study the linear stability properties analytically. We denote the last spiking time of neuron i before t n by 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 nth and (n + 1)th network event by where t t i j k j k i j r s = +τ is the reception time of the spike sent at t jk s 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 n and therefore infl uence the potential 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 nth event, neuron i would reach the threshold at some time t′ under the assumption that there are no further inputs after t n . According to Eq. 52, t′ is implicitly given by where 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 as The linearized estimation of the displacement Δt′of t′ is then given by The 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 k for k ∈ {n − W,…,n}, yields the identity 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, suffi ciently small perturbations stay bounded by the initial ones for fi nite 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 infl uence of various events on V i,n (A,t′) is as follows: For an infl uential spike receiving t k , Eq. 52 yields 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 k a displacement of t k has no infl uence on Therefore condition (58) reduces to a condition on the left and right hand side of Eq. 60 and can be reformulated as where t k are the spikes arrival times at neuron 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 d after receiving a spike are stable. For β → ∞ the system tends to the δ -pulse coupled system and indeed T d vanishes, lim , β→∞ T d = 0 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 11A).
To illustrate our analytical fi ndings, in Figures 11B,C, we used a generalization of a recently introduced method (Memmesheimer and Timme, 2006a,b) to design two networks realizing predefi ned spike patterns in a network of fi ve neurons (V Θ,i ≡ 1.0, I i ≡ 2.4, β = 8, τ := τ 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 d + τ, which ensures that a neuron never spikes within a time period T d after receiving a spike; in (C) we choose the ISIs smaller ΔT = (T d + τ)/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 fi ve 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 −12 leads to a divergence from the unstable periodic orbit.

DISCUSSION
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-fi re 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 infi nitesimal 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 fi eld 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 infi nitely many neurons in sparse networks (v.Vreeswijk andSompolinsky, 1996, 1998;Amit and Brunel, 1997;Brunel and Hakim, 1999;Brunel, 2000). In this mean fi eld 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 andSompolinsky, 1996, 1998). Further studies of fi nite networks found stable dynamics in weakly diluted networks of inhibi-tory coupled neurons (Zillmer et al., 2006), as well as in globally coupled networks with dominating inhibition (Jin, 2002). Recent analytical evidence confi rmed the existence of stable dynamics in inhibitory coupled networks of integrate-and-fi re 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 fi eld models. Thus, one cannot make strict statements about stability in the limit of infi nitely 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 coeffi cients 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 fi ring 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.

APPENDIX
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 specifi c neuron l 0 through the entire network.
All neurons spike at least once in a suffi ciently large but fi nite time interval T. Moreover, after τ max = max i,j (τ ij ) all spikes in transit have thus a neuron l i ∈ Post i (l 0 ) is connected to l 0 by a directed path of length i (cf. Figure 12). Further, we defi ne Post 0 (l 0 ) := {l 0 }. We estimate the bounds of the perturbation following one specifi c 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 j via neurons 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 Δ i ( ) . 0 The fi rst spiking of neuron l 0 after n = 0 is labeled by s 0 ≤ M. After a delay time τ l l 1 0 this spike is received by the postsynaptic neuron l 1 ∈ Post(l 0 ), we call the event r 1 ≤ 2M. After at most M further events, at s 1 ≤ 3M, the neuron l 1 emits a spike. In general, we recursively defi ne s i as the fi rst spiking event of neuron l i ∈ Post i (l 0 ) after r i and r i as the event when the spike generated by l i − 1 ∈ Post i − 1 (l 0 ) at s i − 1 is received (cf. Figure 12). Due to the defi nition of M, the relations s i ≤ (2i+1)M and r i ≤ 2iM hold.
First, we prove by induction that the perturbation of the neuron l i , before sending of a spike at s i , is bounded from above by Post ( ) contains all neurons of the network (cf. Figure 12). After the s i th event, neuron l i can still receive spikes. Before the s N − 1 th event, taken as reference, it receives in the worst case scenario (s N − 1 − s i ) inputs with maximal initial perturbation Δ max ( ) 0 and minimal averaging factor c min . Using Eq. 65 we repeatedly apply Eq. 14 (s N − 1 − s i ) times which leads to At 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 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 := 2NM ≥ s N − 1 + M events are bounded by