Drifting States and Synchronization Induced Chaos in Autonomous Networks of Excitable Neurons

The study of balanced networks of excitatory and inhibitory neurons has led to several open questions. On the one hand it is yet unclear whether the asynchronous state observed in the brain is autonomously generated, or if it results from the interplay between external drivings and internal dynamics. It is also not known, which kind of network variabilities will lead to irregular spiking and which to synchronous firing states. Here we show how isolated networks of purely excitatory neurons generically show asynchronous firing whenever a minimal level of structural variability is present together with a refractory period. Our autonomous networks are composed of excitable units, in the form of leaky integrators spiking only in response to driving currents, remaining otherwise quiet. For a non-uniform network, composed exclusively of excitatory neurons, we find a rich repertoire of self-induced dynamical states. We show in particular that asynchronous drifting states may be stabilized in purely excitatory networks whenever a refractory period is present. Other states found are either fully synchronized or mixed, containing both drifting and synchronized components. The individual neurons considered are excitable and hence do not dispose of intrinsic natural firing frequencies. An effective network-wide distribution of natural frequencies is however generated autonomously through self-consistent feedback loops. The asynchronous drifting state is, additionally, amenable to an analytic solution. We find two types of asynchronous activity, with the individual neurons spiking regularly in the pure drifting state, albeit with a continuous distribution of firing frequencies. The activity of the drifting component, however, becomes irregular in the mixed state, due to the periodic driving of the synchronized component. We propose a new tool for the study of chaos in spiking neural networks, which consists of an analysis of the time series of pairs of consecutive interspike intervals. In this space, we show that a strange attractor with a fractal dimension of about 1.8 is formed in the mentioned mixed state.

The study of balanced networks of excitatory and inhibitory neurons has led to several open questions. On the one hand it is yet unclear whether the asynchronous state observed in the brain is autonomously generated, or if it results from the interplay between external drivings and internal dynamics. It is also not known, which kind of network variabilities will lead to irregular spiking and which to synchronous firing states. Here we show how isolated networks of purely excitatory neurons generically show asynchronous firing whenever a minimal level of structural variability is present together with a refractory period. Our autonomous networks are composed of excitable units, in the form of leaky integrators spiking only in response to driving currents, remaining otherwise quiet. For a non-uniform network, composed exclusively of excitatory neurons, we find a rich repertoire of self-induced dynamical states. We show in particular that asynchronous drifting states may be stabilized in purely excitatory networks whenever a refractory period is present. Other states found are either fully synchronized or mixed, containing both drifting and synchronized components. The individual neurons considered are excitable and hence do not dispose of intrinsic natural firing frequencies. An effective network-wide distribution of natural frequencies is however generated autonomously through self-consistent feedback loops. The asynchronous drifting state is, additionally, amenable to an analytic solution. We find two types of asynchronous activity, with the individual neurons spiking regularly in the pure drifting state, albeit with a continuous distribution of firing frequencies. The activity of the drifting component, however, becomes irregular in the mixed state, due to the periodic driving of the synchronized component. We propose a new tool for the study of chaos in spiking neural networks, which consists of an analysis of the time series of pairs of consecutive interspike intervals. In this space, we show that a strange attractor with a fractal dimension of about 1.8 is formed in the mentioned mixed state.

INTRODUCTION
The study of collective synchronization has attracted the attention of researchers across fields for now over half a century (Winfree, 1967;Kuramoto, 1975;Peskin, 1975;Buck, 1988;Pikovsky and Rosenblum, 2015). Kuramoto's exactly solvable mean field model of coupled limit-cycles (Kuramoto, 1975), formulated originally by Winfree (1967), has helped in this context to establish the link between the distribution of natural frequencies and the degree of synchronization (Gros, 2010). Moreover, the functional simplicity of this model, and other extensions, has permitted to analytically study the collective response of the system to external perturbations in the form of phase resets (Levnajić and Pikovsky, 2010). Networks of phase coupled oscillators may show, in addition, attracting states corresponding to limit cycles, heteroclinic networks, and chaotic phases (Ashwin et al., 2007;Dörfler and Bullo, 2014), with full, partial, or clustered synchrony (Golomb et al., 1992), or asynchronous behavior (Abbott and van Vreeswijk, 1993) Different degrees of collective synchronization may occur also in networks of elements emitting signals not continuously, such as limit-cycle oscillators, but via short-lived pulses (Mirollo and Strogatz, 1990;Abbott and van Vreeswijk, 1993;Strogatz and Stewart, 1993). Networks of pacemaker cells in the heart (Peskin, 1975), for instance, synchronize with high precision, acting together as a robust macroscopic oscillator. Other well-known examples are the simultaneous flashing of extended populations of southeast Asian fireflies (Hanson, 1978;Buck, 1988) and the neuronal oscillations of cortical networks (Buzsáki and Draguhn, 2004). In particular, the study of synchronization in the brain is of particular relevance for the understanding of epileptic states, or seizures (Velazquez et al., 2007).
The individual elements are usually modeled in this context as integrate and fire units (Kuramoto, 1991;Izhikevich, 1999), where the evolution (in between pulses, flashes, or spikes) of a continuous internal state variable V is governed by an equation of the type: Here τ is the characteristic relaxation timescale of V, with f representing the intrinsic dynamics of the unit, and I the overall input (both from other units and from external stimuli). Whenever V reaches a threshold value V θ , a pulse is emitted (the only information carried to other units) and the internal variable is reset to V rest . These units are usually classified either as oscillators or as excitable units, depending on their intrinsic dynamics. The unit will fire periodically even in the absence of input when f (V) > 0 (∀ V ≤ V θ ). Units of this kind are denoted pulse-coupled oscillators. The unit is, on the other hand, an excitable unit, if an additional input is required to induce firing.
A natural frequency given by the inverse integration time of the autonomous dynamics exist in the case of pulse-coupled oscillators. There is hence a preexisting, albeit discontinuous limit cycle, which is then perturbed by external inputs. One can hence use phase coupling methods to study networks of pulse coupled oscillators (Mirollo and Strogatz, 1990;Kuramoto, 1991;Izhikevich, 1999), by establishing a map between the internal state variable V and a periodic phase φ given by the state of the unit within its limit cycle. From this point of view systems of pulse-coupled units share many properties with sets of coupled Kuramoto-like oscillators (Kuramoto, 1975), albeit with generally more complex coupling functions (Izhikevich, 1999). For reviews and examples of synchronization in populations of coupled oscillators see Strogatz (2000) and Dörfler and Bullo (2014).
These assumptions break down for networks of coupled excitable units as the ones here described. In this case the units will remain silent without inputs from other elements of the system and there are no preexisting limit cycles and consequently also no preexisting natural frequencies (unlike rotators (Sonnenschein et al., 2014), which are defined in terms of a periodic phase variable, and a count with a natural frequency). The firing rate depends hence exclusively on the amount of input received. The overall system activity will therefore forcefully either die out or be sustained collectively in a self-organized fashion (Gros, 2010). The respectively generated spiking frequencies for a given mean network activity could be considered in this context as self-generated natural frequencies.
The study of pulse coupled excitable units is of particular relevance within the neurosciences, where neurons are often modeled as spike emitting units that continuously integrate the input they receive from other cells (Burkitt, 2006). The proposal (Shadlen and Newsome, 1994;Amit and Brunel, 1997), and later the empirical observation that excitatory and inhibitory inputs to cortical neurons are closely matched in time (Sanchez-Vives and McCormick, 2000;Haider et al., 2006), has led researchers to focus on dynamical states (asynchronous states in particular) in networks characterized by a balance between excitation and inhibition (Abbott and van Vreeswijk, 1993;van Vreeswijk and Sompolinsky, 1996;Hansel and Mato, 2001;Vogels and Abbott, 2005;Kumar et al., 2008;Stefanescu and Jirsa, 2008). This balance (E/I balance) is generally presumped to be an essential condition for the stability of states showing irregular spiking, such as the one arising in balanced networks of integrate and fire neurons (Brunel, 2000). The type of connectivity usually employed in network studies however, is either global, or local consisting of either repeated patterns, or random connections drawn from identical distributions (Kuramoto and Battogtokh, 2002;Abrams and Strogatz, 2004;Ashwin et al., 2007;Alonso and Mindlin, 2011). Our results show, however, that only a minimal level of structural variability is necessary for excitatory networks to display wide varieties of dynamical states, including stable autonomous irregular spiking. We believe that these studies are not only interesting on their own because of the richness of dynamical states, but also provide valuable insight into the role of inhibition.
Alternatively, one could have built networks of excitatory neurons with high variability in the connection parameters, reproducing realistic connectivity distributions, such as those found in the brain. The large number of parameters involved would make it however difficult to fully characterize the system from a dynamical systems point of view, the approach taken here. An exhaustive phase-space study would also become intractable. We did hence restrict ourselves in the present work to a scenario of minimal variability, as given by a network of globally coupled excitatory neurons, where the coupling strength of each neuron to the mean field is non-uniform. Our key result is that stable irregular spiking states emerge even when only a minimal level of variability is present at a network level.
Another point we would like to stress here is that asynchronous firing states may be stabilized in the absence of external inputs. In the case here studied, there is an additional "difficulty" to the problem, in the sense that the pulse-coupled units considered are in excitable states, remaining quiet without sufficient drive from the other units in the network. The observed sustained asynchronous activity is hence self-organized.
We characterize how the features of the network dynamics depend on the coupling properties of the network and, in particular, we explore the possibility of chaos in the here studied case of excitable units, when partial synchrony is present, since this link has already been established in the case of coupled oscillators with a distribution of natural frequencies (Miritello et al., 2009), while other studies had also shown how stable chaos emerges in inhibitory networks of homogeneous connection statistics (Angulo-Garcia and Torcini, 2014).

THE MODEL
In the current work we study the properties of the selfinduced stationary dynamical states in autonomous networks of excitable integrate-and-fire neurons. The neurons considered are characterized by a continuous state variable V (as in Equation 1), representing the membrane potential, and a discrete state variable y that indicates whether the neuron fires a spike (y = 1) or not (y = 0) at a particular point in time. More precisely, we will work here with a conductance based (COBA) integrate-and-fire (IF) model as employed in Vogels and Abbott (2005) (here however without inhibitory neurons), in which the evolution of each neuron i in the system is described by: where E ex = 0 mV represents the excitatory reversal potential and τ = 20 ms is the membrane time constant. Whenever the membrane potential reaches the threshold V θ = −50 mV, the discrete state of the neuron is set to y i = 1 for the duration of the spike. The voltage is reset, in addition, to its resting value of V rest = −60 mV, where it remains fixed for a refractory period of t ref = 5 ms. Equation (2) is not computed during the refractory period. Except for the times of spike occurrences, the discrete state of the neuron remains y i = 0 (no spike). The conductance g i in Equation (2) integrates the influence of the time series of presynaptic spikes, decaying on the other side in absence of inputs: where τ ex = 5 ms is the conductance time constant. Incoming spikes from the N − 1 other neurons produce an increase of the conductance g i −→ g i + g i , with: Here the synaptic weights w ij represent the intensity of the connection between the presynaptic neuron j and the postsynaptic neuron i. We will generally employ normalized synaptic matrices with j w ij /(N − 1) = 1. In this way we can scale the overall strength of the incoming connections via K i , retaining at the same time the structure of the connectivity matrix.

Global Couplings
Different connectivity structures are usually employed in the study of coupled oscillators, ranging from purely local rules to global couplings (Kuramoto and Battogtokh, 2002;Abrams and Strogatz, 2004;Ashwin et al., 2007;Alonso and Mindlin, 2011). We start here by employing a global coupling structure, where each neuron is coupled to the overall firing activity of the system: which corresponds to a uniform connectivity matrix without self coupling. All couplings are excitatory. The update rule (Equation 4) for the conductance upon presynaptic spiking then take the form: wherer =r(t) represents the time-dependent mean field of the network, viz the average over all firing activities.r is hence equivalent to the mean field present in the Kuramoto model (Kuramoto, 1975), resulting in a global coupling function as an aggregation of local couplings. With our choice (Equation 5) for the coupling matrix the individual excitable units may be viewed, whenever the mean fieldr is strong enough, as oscillators emitting periodic spikes with an "effective" natural frequency determined by the afferent coupling strength K i . The resulting neural activities determine in turn the mean fieldr(t).

Coupling Strength Distribution
We are interested in studying networks with non-uniform K i , We mostly consider here the case of equidistant K i , defined by: for the N neurons, whereK represents the mean scaling parameter, and K, the maximal distance to the mean. It is possible, alternatively, to use a flat distribution with the K i drawn from an interval [K − K,K + K] around the meanK. For large systems there is no discernible difference, as we have tested, between using equidistant afferent coupling strengths K i and drawing them randomly from a flat distribution.

RESULTS
Several aspects of our model, in particular the asynchronous drifting state, can be investigated analytically as a consequence of the global coupling structure (Equation 5), as shown in Section 3.1. All further results are obtained from numerical simulations, for which, if not otherwise stated, a timestep of 0.01 ms has been used. We have also set the spike duration to one time-step, although these two parameters can be modified separately if desired, with our results not depending on the choice of the time-step, while the spike width does introduce minor quantitative changes to the results, as later discussed.

Stationary Mean-Field Solution for the Drifting State
As a first approach we compute the response of a neuron with coupling constant K i to a stationary mean fieldr, as defined by Equation (6), representing the average firing rate of spikes (per second) of the network. This is actually the situation present in the asynchronous drifting state, for which the firing rates of the individual units are incommensurate. Withr being constant we can combine the update rules (3) and (4) for the conductances g i to where we have denoted with g * i the steady-state conductance. With the individual conductance becoming a constant we may also integrate the evolution Equation (2) for the membrane potential, obtaining the time t * i it takes for the membrane potential V i to reach the threshold V θ , when starting from the resting potential V rest : with: We note, that the threshold potential V θ is only reached, if For the t * i to be finite we hence have (from Equation 9) The spiking frequency is r i = T −1 i , with the intervals T i between consecutive spikes given by T i = t * i + t ref , when (Equation 12) is satisfied. Otherwise the neuron does not fire. The mean fieldr is defined as the average firing frequencȳ of the neurons. Equations (10) and (13) describe the asynchronous drifting state in the thermodynamic limit N → ∞. We denote this self-consistency condition forr the stationary mean-field (SMF) solution.

Numerical Simulations
We studied our model, as defined by Equations (2) and (3), numerically for networks with typically N = 100 neurons, a uniform coupling matrix (see Equation 5) and coupling parametersK and K given by Equation (7). We did not find qualitative changes when scaling the size of the network up to N = 400 for testing purposes (and neither with downscaling), see Figure 7. Random initial conditions where used. The network-wide distribution of firing rates is computed after the system settles to a dynamical equilibrium. Three examples, forK = 2.0 and K/K = 0.9, 0.6, and 0.2, of firing-rate distributions are presented in Figure 1 in comparison with the analytic results obtained from the stationary mean field approach (SMF), as given by Equation (13). The presence or absence of synchrony is directly visible. In all of the three parameter settings presented in Figure 1 there is a drifting component, characterized by a set of neurons with a continuum of frequencies. These neurons fire asynchronously, generating a constant contribution to the collective mean field.
The plateau present in the case K/K = 0.2, corresponds, on the other hand, to a set of neurons firing with identical frequencies and hence synchronously. Neurons firing synchronously will do so however with finite pairwise phase lags, with the reason being the modulation of the common mean field r through the distinct afferent coupling strengths K i . We note that the stationary mean-field theory (Equation 13) holds, as Frontiers in Computational Neuroscience | www.frontiersin.org expected, for drifting states, but not when synchronized clusters of neurons are present.
In Figure 2 we systematically explore the phase space as a function ofK and K. For labeling the distinct phases we use the notation I : inactive, I+D : partially inactive and drifting, D : fully drifting (asynchronous), D+S : mixed, containing both drifting and synchronized components, and S : fully synchronized. Examples of the rate distributions present in the individual phases are presented in Figure 2C.
The phase diagram is presented in Figure 2A. The activity dies out for a low mean connectivity strengthK, but not for largerK. Partial synchronization is present when bothK and the variance K are small, taking over completely for larger values ofK and small K. The phase space is otherwise dominated by a fully drifting state. The network averager of the neural firing rates, given in Figure 2B, drops only close to the transition to the inactive state I, showing otherwise no discernible features at the phase boundaries.
The dashed lines in Figures 2A,B represent the transitions between the inactive state I and active state I+D, and between states I+D and D, as predicted by the stationary mean field approximation (Equation 13), which becomes exact in the thermodynamic limit. The shaded region in these plots indicates the co-existence of attracting states S and S+D. As a note, we found that the location of this shaded region depends on the spike width, shifting to higherK values for narrower spikes. While real spikes in neurons have a finite width, we note from a dynamical systems point of view, that this region would most likely vanish in the limit of delta spikes.
For a stable (non-trivial) attractor to arise in a network composed only of excitatory neurons, some limitation mechanism needs to be at play. Otherwise one observes a bifurcation phenomenon, similar to that of branching problems, in which only a critical network in the thermodynamic limit could be stable (Gros, 2010). In this case, the limiting factor is the refractory period. Refractoriness prevents neurons from firing continuously, and prevents the system activity from exploding. Interestingly, this does not mean that the neurons will fire at the maximal rate of 1/t ref which would correspond in this case to 200 Hz. The existence of this refractory period allows for self-organized states with frequencies even well bellow this limit, as seen in Figure 2B. We have tested these claims numerically by setting t ref = 0, observing that the neural activity either dies out or the neurons fire continuously.
In order to study the phase transitions between states D and D+S and between D+S and S, we will resort in the  following section to adiabatic paths in phase space crossing these lines.

Adiabatic Parameter Evolution
Here we study the nature of the phase transitions between different dynamical states in Figure 2. To do so, we resort to adiabatic trajectories in phase space, crossing these lines. Beginning in a given phase we modify the coupling parameters K and K (Equation 7), on a timescale much slower than that of the network dynamics. Along these trajectories, we then freeze the system in a number of windows in which we compute the rate distribution as a function of the K rank (see Figure 1). During these observation windows the parameters do not change. In this way, we can follow how the rate distribution varies across the observed phase transitions. The results are presented in Figure 3. We observe that the emergence of synchronized clusters, the transition D→(D+S), is completely reversible. We believe this transition to be of second order and that the small discontinuity in the respective firing rate distributions observed in Figure 3A are due to finite-size effects. The time to reach the stationary state diverges, additionally, close to the transition, making it difficult to resolve the locus with high accuracy.
The disappearance of a subset of drifting neurons, the transition S→(D+S) is, on the other hand, not reversible. In this case, whenK is reduced, the system tends to get stuck in metastable attractors in the S phase, producing irregular jumps in the rate distributions. Furthermore, when we increaseK, we observe that the system jumps back and forth between states D+S and S in the vicinity of the phase transition, indicating that both states may coexist as metastable attractors close to the transition. We note that a similar metastability has been observed in partially synchronized phase of the Kuramoto model (Miritello et al., 2009).

Time Structure
In networks of spiking neurons, it is essential to characterize not only the rate distribution of the system, but also the neurons' interspike-time statistics (Perkel et al., 1967a,b;Chacron et al., 2004;Lindner, 2004;Farkhooi et al., 2009). In this case, we have computed the distribution p i (s) of the interspike intervals s (ISI) of the individual neurons respectively for full and partial drifting and synchronized states. The distribution of inter-spike intervals in Figure 4 shows the network average of the p i (s), normalized individually with respect to the average T i = s p i (s)ds spiking intervals.
• D : The input received by a given neuron i tends to a constant, as discussed in Section 3.1, in the thermodynamic limit N → ∞. The small but finite width of the ISI for the fully drifting state D evident in Figure 4 is hence a finite-size effect. • D+S : The input received for drifting neuron i in a state where other neurons form a synchronized subcluster is intrinsically Of interest here are the finite CVs of the drifting units in the D+S state, which are considerably larger than the CVs of the drifting neurons when no synchronized component is present. This phenomenon is a consequence of the interplay between the periodic driving of the drifting neurons by the synchronized subcluster in the D+S state, where the driving frequency will in general be in mismatch with the effective, self-organized natural frequency of the drifting neurons. The firing of a drifting neuron is hence irregular in the mixed D+S state, becoming however regular in the absence of synchronized drivings.

Self Induced Chaos
The high variability of the spiking intervals observed in the mixed state, as presented in Figure 4, indicates that the firing may have a chaotic component in the mixed state and hence positive Lyapunov exponents (Gros, 2010).
Alternatively to a numerical evaluation of the Lyapunov exponents (a demanding task for large networks of spiking neurons), a somewhat more direct understanding of the chaotic state can be obtained by studying the relation between consecutive spike intervals. In Figure 5 we plot for this We note that the spiking would be regular, viz s i (n) constant, for all neurons either in the fully drifting state (D) or in the fully synchronized state (S). The plot of consecutive spike intervals presented in Figure 5 can be viewed as a poor man's approximation to Takens' embedding theorem (Takens, 1981), which states that a chaotic attractor arising in a d-dimensional phase space (in our case d = 2N) can be reconstructed by the series of d-tuples of time events of a single variable. With a blue line we follow in Figure 5 a representative segment of the trajectory, which jumps irregularly. A first indication that the attractor in Figure 5 may be chaotic comes from the observation that the trajectory does not seem to settle (within the observation window) within a limit cycle. The time series of consecutive spike-interval pairs will nevertheless approach any given previous pair s i (n), s i (n + 1) arbitrarily close, a consequence of the generic ergodicity of attracting sets (Gros, 2010). One of these close re-encounters occurs in Figure 5 near the center of the dashed circle, with the trajectory diverging again after the close re-encounter. This divergence indicates the presence of positive Lyapunov exponents.
We have determined the fractal dimension of the attracting set of pairs of spike intervals in the mixed phase by overlaying the attractor with a grid of 2 r × 2 r squares. For this calculation we employed a longer simulation with N spikes = 128, 000. The resulting box count, presented in Figure 6, yields a Minkowski or box-counting dimension D F ≈ 1.8, embeded in the 2D space of the plot, confirming such that the drifting neurons in the D+S phase spikes indeed chaotically. As a comparison, a limit cycle in where n denotes the spike number. In blue, we follow a representative segment of the trajectory. The system does not appear to follow a limit cycle, and these preliminary results would suggest the presence of chaos in the D+S state, consistent with studies of chaos in periodically driven oscillators (d 'Humieres et al., 1982). Indeed, if one looks at the two close points within the dashed circle, we observe how an initially small distance between them, rapidly grows in a few iterations steps, indicating a positive eigenvalue. For this simulation we have used a time-step dt = 0.001 ms, to improve the resolution of the points in the plot. We have evaluated the time the neuron needs to circle the attractor, finding it to be of the order of ∼ 5.3 spikes. Other drifting neurons take slightly longer or shorter. In Figure 6, we compute the fractal dimension of the here shown attractor.
this space, has a D F of 1. While we present here the result for one particular neuron, the same holds for every drifting neuron in this state, albeit with slightly different fractal dimension values. We note that the such determined fractal dimension is not the one of the full attractor in d = 2N phase space, for which tuples of 2N consecutive inter-spike intervals would need to be considered (Takens, 1981;Ding et al., 1993). Our point here is that a noninteger result for the single neuron D F strongly indicates that the full attractor (in the d-dimensional phase space) is chaotic.
We believe that the chaotic state arising in the mixed D+S state may be understood in analogy to the occurrence of chaos in the periodically driven pendulum (d 'Humieres et al., 1982). A drifting neuron with a coupling constant K in the D+S does indeed receive two types of inputs to its conductance, compare Equation (4), with the first input being constant (resulting from the firing of the other drifting neurons) and with the second input being periodic. The frequency r syn of the periodic driving will then be strictly smaller than the natural frequency r K of the drifting neuron as resulting from the constant input (compare Figure 1). It is known from the theory of driven oscillators (d 'Humieres et al., 1982) that the oscillator may not be able to synchronize with the external frequency, here r syn , when the frequency ratio r syn /r K is small enough and the relative strength of the driving not too strong.  Figure 5. N box denotes the number of squares occupied with at least one point of the trajectory of consecutive pairs of spike intervals, when a grid of 2 r × 2 r squares is laid upon the attracting set shown in Figure 5. The fractal dimension D F = log(N box )/log(r) is then ∼1.8. A time series with N spikes = 128000 spikes for the same drifting neuron in the same D+S state has been used. log 2 (N box ) saturates at log 2 (N spikes ) ≈ 16.97, observing that the linear range can be expanded further by increasing the number of spikes, albeit with an high cost in simulation time. Finally the resolution of the method is limited by the spike width.

Robustness
In order to evaluate the robustness and the generality of the results here presented, we have evaluated the effects occurring when changing the size of the network and when allowing for variability in the connectivity matrix w ij . We have also considered an adiabatically increasing external input, as well as Gaussian noise.
In Figure 7 (top half), the effect on the rate distribution of the network size is evaluated. Sizes of N = 50, 100, 200, and 400 have been employed. We observe that the plots overlap within the precision of the simulations. This result is on the one hand a consequence of the scaling K i ∼ 1/N for the overall strength of the afferent links and, on the other hand, of the regularity in firing discussed in Section 3.2.2. The neural activity is driven by the mean fieldr(t) which is constant, in the thermodynamic limit N → ∞, in the fully drifting state and non-constant but smooth (apart from an initial jump) in the synchronized states. Fluctuations due to finite network sizes are already small for N ≈ 100, as employed for our simulations, justifying this choice.
In the previous sections, we considered the uniform connectivity matrix described by Equation (5). This allowed us to formulate the problem in terms of a mean-field coupling. We now analyze the robustness of the states found when a certain degree of variability is present in the weight matrix, viz when an extra variability term η is present: Here we consider η to be drawn from a flat distribution with zero mean and a width W. Tests with W = 10, 20, and 50% were performed. In Figure 7 (lower half), the results for W = 50% are presented. We observe that the fully drifting state is the least affected by the variability in the weight matrix. On the other hand, the influence of variable weight matrices becomes more evident when the state is partially or fully synchronized, with the separation between the locked and the drifting neurons becoming less pronounced in the case of partial synchronization (lower panel of Figure 7B). The larger standard deviation evident for larger values of K rank in the lower panel of Figure 7C indicates the presence of drifting states in some of the ensemble realizations of weight matrices. Finally, we test the robustness of the drifting state when perturbed with an external stimulus. To determine the stability of the state, we adiabatically increase the external stimulus I ext and compute the firing rate as a function of the rank for several values of I ext . We do two excursions, one for positive values of I ext and another one for negative values. These plots are presented in Figure 8. We observe that the firing rates evolve in a continuous fashion, indicating that the drifting state is indeed stable. While positive inputs push the system to saturation, negative inputs reduce the average rate. We find, as is to be expected, that a large enough negative input makes the system silent. As a final test (not shown here), we have perturbed the system with random Gaussian uncorrelated noise, observing that the found attractors are all robust with respect to this type of noise as well.

DISCUSSION
In the present work we have studied a network of excitable units, consisting exclusively of excitatory neurons. In absence FIGURE 8 | Effect of an external input (either positive or negative), on the neural firing rates. The input is either increased (blue arrow) or decreased (red arrow) from zero (green curve) adiabatically. The drifting state remains stable for a wide range of inputs, with the activity disappearing only for I ext < −7.0 (gray curve, coinciding with the x-axis).
of external stimulus, the network is only able to remain active through its own activity, in a self-organized fashion. Below a certain average coupling strength of the network the activity dies out, whereas, if the average coupling is strong enough, the excitable units will collectively behave as pulse-coupled oscillators.
We have shown how the variability of coupling strengths determines the synchronization characteristics of the network, ranging from fully asynchronous, to fully synchronous activity. Interestingly, this variability, together with the neurons' refractoriness, is enough to keep the neural activity from exploding.
While we have initially assumed a purely mean field coupling (by setting all the synaptic weights w ij = 1), only regulating the intensity with which a neuron integrates the mean field with the introduction of a scaling constant K i , we have later shown how the here found states also survive when we allow the w ij s to individually vary up or down by up to a 50% value. We have also shown how the variability in coupling strength makes the asynchronous or drifting state extremely robust with respect to strong homogeneous external inputs.
Finally, we have studied the time structure of spikes in the different dynamical states observed. It is in the time domain that we find the main difference with natural neural networks. Spiking in real neurons is usually irregular, and it is often modeled as Poissonian, whereas in our network we found a very high degree of regularity, even in the asynchronous state. Only in the partially synchronous state we found a higher degree of variability, as a result from chaotic behavior. We have determined the fractal dimension of the respective strange attractor in the space of pairs of consecutive interspike intervals, finding fractional values of roughly 1.8 for the different neurons in the state.
While it has been often stated that inhibition is a necessary condition for bounded and uncorrelated activity, we have show here that uncorrelated aperiodic (and even chaotic) activity can be obtained with a network of excitatory-only connections, in a stable fashion and without external input. We are of course aware that the firing rates obtained in our simulations are high compared to in vivo activity levels and that the degree of variability in the time domain of spikes is far from Poissonian. We have however incorporated in this work only variability of the inter-neural connectivity, keeping other neural properties (such as the neural intrinsic parameters) constant and homogeneous.
In this sense, it would be interesting to study in future work, how intrinsic and synaptic plasticity (Triesch, 2007), modify these statistics, incorporating plasticity in terms of interspike-times (Clopath et al., 2010;Echeveste and Gros, 2015b), and in terms of neural rates (Bienenstock et al., 1982;Hyvärinen and Oja, 1998;Echeveste and Gros, 2015a). Here, instead of trying to reproduce the detailed connection statistics in the brain, which would in any case never be realistic without inhibitory neurons, we have shown how a minimal variability model in terms of non uniform link matrices is able to give rise to asynchronous spiking states, even without inhibition. Our results indicate therefore that further studies are needed for an improved understanding of which features of the asynchronous spiking state depend essentially on inhibition, and which do not.
We have shown here that autonomous activity (sustained even in the absence of external inputs) may arise in networks of coupled excitable units, viz for units which are not intrinsically oscillating. We have also proposed a new tool to study the appearance of chaos in spiking neural networks by applying a box counting method to consecutive pairs of inter-spike intervals from a single unit. This tool is readily applicable both to experimental data and to the results of theory simulations in general.

AUTHOR CONTRIBUTIONS
RE carried out the simulations and produced the figures. CG guided the project and provided the theoretical framework. Both authors contributed to the writing of the article.