We study ring networks of N leaky integrate-and-fire (LIF) neurons with current-based synapses. N_E = βN, β ∈ [0, 1], of all neurons are excitatory, the residual N_I = N − N_E neurons are inhibitory. The membrane potential V_i(t) of neuron i is governed by the differential equation (1)τmdVi(t)dt=−Vi(t)+RIs,i(t−d)+RIx,i(t), where the membrane potential V_i(t) is reset to V_res whenever it reaches the threshold potential V_thr and a spike is emitted. The neuron is then refractory for a time τ_ref. τ_m denotes the membrane time constant, RIs,i(t−d)=τm∑j = 1, i ≠ jNWijsj(t−d), where I_s is the input current received from the local network, and s_j(t) = ∑_k δ(t − t_{j, k}) is the spike train produced by some neuron j. d is the transmission delay, and W is the coupling matrix, with entries that are either zero (no synapse present), J_E = J (excitatory synapse) or J_I = −gJ (inhibitory synapse). To keep spiking activity from entering a high rate state, we assume that the recurrent input I_s is net inhibitory by choosing g > N_E/N_I.

I_x denotes the external input current that is modeled as stationary Poisson noise injected via current-based synapses of strength J_x (for details, see below), and R is the membrane resistance. Such leaky integrate-and-fire neurons were extensively studied in various input (see e.g., Amit and Tsodyks, 1991; Amit and Brunel, 1997; Lindner, 2004; Fourcaud-Trocmé and Brunel, 2005; De la Rocha et al., 2007; Vilela and Lindner, 2009; Helias et al., 2010) and network settings, especially in the context of balanced random networks (Brunel, 2000; Tetzlaff et al., 2012). Here, we consider ring networks of LIF neurons, such that each fifth neuron is inhibitory and the coupling matrix W is given by (cf. Figure 3A) (2)Wij={Jif j excitatory, |i−j|modN≤κ/2 and i≠j−gJif j inhibitory, |i−j|modN≤κ/2 and i≠j.0otherwise

Thus, each neuron is connected to its κ nearest neighbors, where we assume κ = 0.1N. We note that the network layout chosen here is motivated by an actual embedding of individual neurons into some space, in contrast to earlier studies of pattern formation in spiking neuron networks where inhibitory and excitatory neurons are often distributed in equal numbers along two identical rings, and where the main difference lies in the spatial interaction ranges. These systems are usually motivated by population-density descriptions and can often be reduced to a two-dimensional model (see e.g., Ben-Yishai et al., 1995; Compte et al., 2000; Shriki et al., 2003; Roxin et al., 2005). As we will demonstrate in section 3.2, the network studied here leads to an effectively five-dimensional model. In section 3.4 we will discuss a reduction to a two-dimensional model that however already deviates noticeably in its details from the full five-dimensional network.

In the case of the mean-driven (μ[RI] ≥ V_thr − V_res) and intermediate regimes, which will be discussed later, the external current is purely excitatory, such that RI_x(t) = τ_mJ_xs_x(t) with J_x > 0 and rate ν_X = E[s_x(t)], where E[.] denotes the expectation value. To quantify the effective strength of the external input, we express ν_X = ηθ/J_xτ_m, i.e., the parameter η gives the strength of the external drive in terms of the excitatory rate necessary to reach threshold. Assuming stationary spiking activity with rate ν_o for both the excitatory and inhibitory population, i.e., ν_E = E[s_E(t)] = ν_o and ν_I = E[s_I(t)] = ν_o respectively, the mean and variance of the total input RI = R(I_s + I_x) are thus given by (3)μ[R(Is+Ix)]=∑j∈{E,I,X}τmνjJj=κτmνoJ(β−g(1−β))+τmνXJxσ2[R(Is+Ix)]=∑j∈{E,I,X}τmνjJj2=κτmνoJ2(β+g2(1−β))+τmνXJx2.

For stronger local coupling strength J the negative mean recurrent contribution μ[RI_s] to the total input is stronger, while the variance σ²[RI_s] increases. For the networks we consider here, for small J the network activity is thus typically driven by the mean μ[RI] of the total input, while for larger J the variance σ²[RI] has a notable impact. Moreover, the network firing rate ν_o depends on J.

The strongly fluctuation-driven regime is characterized by subthreshold mean total input, μ[RI] < V_thr − V_res, but high variance σ²[RI], such that spikes are initiated by transient input fluctuations. To guarantee that we are in the high variance regime for all J, we inject excitatory and inhibitory external currents with rates ν_Ex, ν_Ix, such that the total input current mean and variance, i.e., μ[RI] and σ²[RI], stay the same when varying the local coupling strength J. This is achieved by adapting μ[RI_x] = μ[RI] − μ[RI_s] and σ²[RI_x] = σ²[RI] − σ²[RI_s] as a function of ν_Ex, ν_Ix. The external current injected into each local neuron by current-based synapses of strengths J_Ex = J_x and J_Ix = −gJ_x is thus given by (4a)RIx(t)=τmJx(sEx(t)−gsIx(t)) (4b)with νEx:=E[sEx(t)]=1τmJx(1+g)(σ2[RIx]Jx+gμ[RIx]) (4c)and νIx:=E[sIx(t)]=1τmJxg(1+g)(σ2[RIx]Jx−μ[RIx]).

Because the variance of the recurrent input σ²[RI_s] increases quickly with increasing J, the total variance σ²[RI]− in order to keep it fixed for all J− usually needs to be chosen quite large in order to avoid non-sensical rates ν_Ix < 0. This moreover implies that this approach of choosing the external input is not useful for the mean-driven scenario.

For all neurons we set θ = V_thr − V_res = 20 mV, τ_m = 20 ms, and R = 80 MΩ. The excitatory synaptic coupling strengths J, J_x, the relative strength of inhibition g, as well as τ_ref and d will be specified individually. The model and parameters are listed in Table A1 in Appendix A1. All simulations were carried out with NEST (Gewaltig and Diesmann, 2007).

The analysis of the stability of spatially homogeneous activity dynamics to spatial perturbations follows a simple concept: first, the stationary, homogeneous state is determined in a mean-field approximation. In particular, in the diffusion limit, i.e., under the assumption that the coupling strength J is small compared to the distance between reset and threshold θ and that all neurons receive statistically identical input, the stationary firing rate ν_o of the neurons can be derived self-consistently as a function of mean and variance of the input current (a solution of the mean first-passage-time problem first derived by Siegert (Siegert, 1951), but also see e.g., Amit and Tsodyks, 1991; Brunel, 2000), such that (5)νo−1=τref+τmπ∫Vres−μoσoVthr−μoσoexp[x2](1+erf[x])dx, where μ_o = ∑_{i ∈ {E,I,X}} J_iν_{i, o}τ_m and σ²_o = ∑_{i ∈ {E,I,X}} J²_iν_{i, o}τ_m are the self-consistent input mean and variance. Because all neurons receive statistically identical input, in absence of symmetry-breaking the firing rate distribution is spatially homogeneous, i.e., ν_{o, i} ≡ ν_o for all i. To assess simple stability of the homogeneous state to spatial perturbations, we need to perform linear perturbation analysis. The critical eigenvector of the resulting linear stability operator then determines the spatial pattern that develops.

We will show that the exact nature of the appropriate linearization depends critically on the input current regime, and derive the respective linearization in the limits of subthreshold mean input with fluctuation-driven spiking in section 3.1.1, and of dominant suprathreshold mean total input current μ[R(I_x + I_s)] > θ in section 3.1.2. In the fluctuation-driven regime a quantitative approximation is obtained from the modulation of both mean and fluctuations of the synaptic input to the neurons, while in the mean-driven regime we recover the intuitive result that the system is governed by the noise-free input-output relation given by the f-I-curve of the individual neuron. In general, the full spiking dynamics will lie somewhere in between these two limiting cases, the case we call “intermediate.”

In this section we will analytically derive the eigenvalue spectrum of the ring networks under consideration that yields the critical eigenvalue and thus the critical coupling strength J_c with respect to the two linearizations Equation (9) and Equation (13). As described in section 2 we consider ring networks of size N with regular connectivity, in that each neuron is connected to its κ nearest neighbors, κ/2 on each side of the neuron but not to itself. Inhibitory neurons shall be distributed periodically across the network as illustrated in Figure 3A, where the ratio between excitation and inhibition is 4 : 1 (β = 0.8). This is in line with the observed frequencies of excitatory and inhibitory cells in cortex (Schüz and Braitenberg, 1994). The ring can thus be divided into structurally identical elementary cells of four excitatory neurons and one inhibitory neuron. We choose κ as an integer multiple of the size ℓ of an elementary cell (here ℓ = 5, κ mod ℓ = 0), such that it moreover divides the total number of neurons N (N mod κ = 0). The resulting coupling matrix is sketched in Figure 3B, where black squares depict inhibitory weights, white depict excitatory weights and gray denotes the lack of a connection. In this way all neurons receive the exact same amount of inhibition and excitation, and the translation invariant mode (1, 1, ···, 1)ᵀ is an eigenvector of the coupling matrix with eigenvalue μ = κJ(β − g(1 − β)). To keep network activity balanced, inhibitory synapses are assumed to be g-times stronger, with g > 4, cf. Equation (2). Hence, the local network is inhibition dominated.

To compute the eigensystem of the N-dimensional system we make use of its symmetry properties. The coupling matrix commutes with the unitary operator that shifts all neuron indices by multiples of ℓ = 5. Hence, we can diagonalize both in the same eigenbasis and moreover reduce the problem to an effectively five-dimensional one as outlined in Appendix A3. Figure 3C shows the exemplary eigenvalue spectrum of in this case the rescaled coupling matrix W/θ of a network of size N = 2500 with absolute coupling strength J = 1 mV and relative strength of inhibition g = 6. Thus, the eigenvalue λ_c of W/θ with largest real part will exceed unity, if the amplitude of the absolute coupling strength J exceeds J^md_c = 1/Re[λ_c], cf. Equation (16), here J^md_c ≈ 0.5.

The corresponding critical weight J^fd_c for the fluctuation-driven case is obtained by the identical eigenvalue decomposition of W˜ as defined by Equation (7). Since W˜ depends non-trivially on the self-consistent working point ν_o of the system, the critical coupling strength J^fd_c is the solution of (17)Jcfd:=J such thatmaxi[Re[λi(W˜(νo,J))]]≡1, where λ_i, i ∈ {1, …, N} are the N eigenvalues of W˜. If J > J^md/fd_c (depending on the regime), the spatially homogeneous rate distribution becomes unstable to perturbations in favor of the fastest growing eigenmode v_c, which only depends on the symmetry properties of the matrix as long as g and g˜ are larger than four, and it is hence the same irrespective of the input regime, cf. Figures 1C,F,I. The spatial frequency of the emerging spatial pattern (see Figure 3), i.e., the number of maxima of the rate distribution along the ring, is thus given by the wavenumber Λ (see Appendix A3) that corresponds to the respective λ_c := max_i[Re[λ_i]]. We note, that the eigenvalue spectra in the fluctuation-driven case look slightly different from those in the mean-driven case (cf. Figure 3C) because of the in general different scaling g˜ ≠ g between positive and negative entries in W˜. The maximum of the dispersion relation determining Λ_c and the eigenvectors are, however, usually the same.

Figure 4 shows the critical eigenvalue λ_c as a function of J for the linearizations Equation (9) (gray lines) and Equation (13) (red line) in the constant external drive regime, parameterized by η. The critical eigenvalue J_c is given by the intersection with the horizontal line at one. As was pointed out in section 3.1.2, in the noiseless approximation Equation (13) the eigenvalues of the linear operator do not depend on the working point, and thus the critical weight J^md_c is unique for all η such that ν_o > 0 (here, critical J^md_c = 0.506 indicated by red circle).

For the linearization Equation (9), and taking into account both μ and σ in Equation (7b) (solid lines), the critical weight J^fd_c is strongly dependent on the value of η = {3.5, 5, 10, 15, 20, 30, 40} (increasing from dark to light gray), and never comes close to the prediction of the noiseless linearization J^md_c. In particular, for larger η J^fd_c increases again.

If we neglect the σ-dependence in Equation (7b), and just take into account dνo/dμo = dνmdo/dRI = 1/τmθ as given by Equation (11), Equation (7) is linear in Wij and we recover the linear noiseless approximation Equation (13) for Equation (9). We might thus expect that Equation (9) in general will give a Jfdc-prediction close to Jmdc if we neglect the σ-dependent quadratic term in Equation (7b). Although the prediction of this Jfd|dνo/dμoc becomes smaller (cf. black triangle for Jfdc = 1.54 mV vs. the gray asterisk for Jfd|dνo/dμoc = 0.89 mV for η = 3.5) it never becomes as small as Jmdc. This is because every change in μo will also influence σo via its impact on νo, and thus dνo(μo, σo)/dμo ≠ dνmdo(RI)/dRI.

We note that the observation Jfdc > Jfd,dνo/dμoc is explained by the fact that W˜(Wij)|dνo/dμo is linear in Wij, and thus the ratio g˜ = |W˜(−gJ)/W˜(J)| = g, while for the full Equation (7), W˜(Wij) is quadratic in Wij, yielding a g˜ < g. Evaluating Equation (A20) in Appendix A3 shows that the critical eigenvalue λc in both linearizations is a linearly increasing function of g, g˜, respectively, and Jc is thus decreasing as a function of g, g˜.

For comparison, Figure 4 also shows the J-dependence of the real part of the critical eigenvalue λ_c in the strongly fluctuation-driven regime (μ[RI] = 5 mV, σ[RI] = 60 mV, blue line), with the critical J_c indicated by the crossing of the unity-line (blue triangle).

We will now compare the general performance of the linearizations Equations (9) and (13) and in particular the resulting predictions of the critical coupling strength J_c with the actual onset of pattern formation as observed in simulations. Figure 1 demonstrates how the identical ring network structure undergoes pattern formation at very different values of network coupling strength J depending on the input current regime. In particular, pattern formation sets in later if the system is in a fluctuation-driven regime. This can be understood in the two limit-cases that were analyzed in sections 3.1.1 and 3.1.2.

But what determines pattern formation onset in intermediate cases that comprise the most interesting and relevant cases? As we saw in section 3.1.2, the network will even undergo a crossover from the mean-driven to a more fluctuation-driven regime, if the external input is kept constant and not adapted to counterbalance the changes in recurrent input structure for changing J, see Figure 2C. Also, due to the spiking nature of the activity, even in a strongly mean-driven scenario there will always be a considerable amount of variance which might change the onset of pattern formation with respect to J^md_c.

To give more insight into the actual onset of pattern formation, useful reduced measures are the variance and kurtosis of the distribution of firing rates over neurons, cf. Figure 5. Figures 5A–C show the histogram of rates, and the variance and kurtosis of these histograms for various J-values in the mean-driven regime averaged over 10 trials each. Figures 5D–F show the same for the fluctuation-driven case and Figures 5G–I for the intermediate case, all for the same parameters as in Figure 1.

When the system is still well in the linear regime J < Jc the distribution is approximately Gaussian with small variance and a kurtosis close to zero (cf. Figures 5A,D,G, light gray curves), while once the critical weight Jc is surpassed the variance increases strongly (cf. Figures 5B,E,H) and the kurtosis becomes negative (platycurtic), indicating the broadening of the distribution due to the inhomogeneous rate per neuron (cf. Figures 5C,F,G and 5A,D,G, dark gray curves). The respective critical weights are indicated by vertical dashed lines in Figures 5B,E,H and C,F,I for visual guidance (red: Jmdc, black: Jfdc, gray: Jfd|dνo/dμoc).

The onset of pattern formation in the strongly fluctuation-driven regime at J^fd_c is more shallow than in the mean-driven regime at J^md_c. This might be due to deviations of the input spike statistics from the asynchronous-irregular activity assumption underlying the linear response derivation of J^fd_c: the standard-deviation of the input is extremely large at σ = 60 mV. This implies that the membrane-potential will be hyperpolarized for long times, while at other times it is highly depolarized, and several spikes are emitted in short succession. This typically leads to higher coefficients of variation, i.e., more irregular firing than expected for Poisson spike statistics [see also Brunel (2000) and supplementary material section S1]. Also, even in the highly fluctuation-driven regime there may be residual correlations between spike trains. Lastly, the increased network coupling strength might lead to deviations from the Gaussian white noise approximation of input currents underlying Equation (5) and the linear response theory yielding Equation (7) (see supplementary material section S1 for an analysis). How such deviations of the input statistics can indeed change the spike response behavior of LIF neurons was studied, e.g., in Moreno et al. (2002); Renart et al. (2007); Moreno-Bote et al. (2008); Helias et al. (2010).

However, another contributing factor, and more likely the reason for the shallow transition at J^fd_c, is the degeneracy of the maximal eigenvalue, see Figures 3D,E. In the fluctuation-driven regime the system is subject to strong perpetual perturbations that lead to a switching between the two dominant eigenvectors depicted in Figure 3D in the transition regime J^fd_c ± ΔJ. Indeed, as can be seen from Figure 1E, the periodic pattern is already clearly distinguishable at J^fd_c over several hundred milliseconds, but the average activity depicted in the histogram is still quite flat, yielding smaller variance and kurtosis of the respective rate distribution.

Finally, Figures 5D–F demonstrate clearly how pattern formation occurs for intermediate J-values J^md_c ≲ J^inter_c < J^fd_c if the system is neither clearly in the mean- nor strongly fluctuation-driven regime, and even changes from one to the other regime with increasing coupling strength, cf. Figure 2C. As we demonstrate in the supplementary material section S1, most input and network settings have pattern formation onsets that usually agree much better with the J^md_c prediction than with J^fd_c. This is a very consistent finding, even in cases where neurons are well in the fluctuation-driven regime in terms of subthreshold mean and pronounced input variance.

As we discuss in the supplementary material section S1 the reason for this finding lies in an asymmetry between the excitatory and inhibitory compound input current statistics. The excitatory subpopulation tends to be synchronized, even in the presence of strong balanced external noise, while inhibition actively decorrelates itself. The explanation for this effect lies in the local connectivity of ring networks, leading to a population-specific pattern of correlations already below the critical coupling: excess synchrony of the excitatory population will reinforce further spiking, while spikes from the inhibitory population decreases the instantaneous firing probability [see also Helias et al. (2013) for a related discussion of this effect in the framework of balanced random networks]. Near synchrony is thereby effectively enhanced in the excitatory population and suppressed in the inhibitory one. Volleys of excitatory input spikes act like compound pulses with large amplitude (on the order of up to θ) in an otherwise balanced asynchronous background activity. At these amplitudes the effective gain Equation (7) derived from linear response theory becomes linear in J with a slope proportional to 1/θ, i.e., the slope of the linear model Equation (11), see supplementary material section S1. This explains the fact that Equation (13) has better predictive power also in the intermediate or—comparably weakly—fluctuation-driven scenarios. Moreover, the fluctuation-driven prediction only becomes valid if the additional external noise sufficiently dilutes residual synchrony, such as it is the case for μ[RI] = 5 mV and σ[RI] = 60 mV.

In section 3.2 we saw that the full system can effectively be reduced to a five-dimensional system. However, the computation of eigensystems (cf. Appendix A3) is still quite involved. In this section we study how well a further simplified system that is formally identical to the well-known Ermentrout-Cowan networks [Ermentrout and Cowan (1979a,b, 1980) and supplementary material section S2] predicts the dynamics of the full network.

We coarse-grain the system and combine groups of ℓ = N/N_I neurons, such that they contain four excitatory and one inhibitory neuron (ℓ needs to be odd in the following), as indicated in Figure 6. The neurons i_c ∈ {0, …, ℓ − 1} within each cell c ∈ {0, …, C − 1}, C = N/ℓ, are connected to every other neuron within their cell c but not to themselves. Moreover, all neurons within one local cell c are connected to all neurons within the K < C/2 neighboring cells to the left and to the right, i.e., {(c − K) mod C, …, (c − 1) mod C, (c + 1) mod C, …, (c + K) mod C}. The computation of the eigensystem is outlined in Appendix A4. Due to the two mirror symmetries corresponding to exchange of neurons within the unit cell (with corresponding eigenvalues r, s ∈ {−1, 1}) and the translational symmetry with the eigenvalue determined by the wavenumber α, the system can be reduced to effectively two dimensions.

To conclude, we study how robust the findings of the previous sections are to effects of randomness in either structure or weight distribution. Though we again only show the eigenvalue spectra for the synaptic coupling matrix W/θ, that is the relevant linear operator in the mean-driven case, all results map straight-forwardly to the fluctuation-driven regime with appropriate rescaling of J and g.

The very nature of spiking neuron networks leads to considerable input current fluctuations, and hence in general both mean μ[RI] and standard deviation σ[RI] of the input will shape the output rate of the neuron. For leaky integrate-and-fire neurons the output rate as a function of μ[RI] and σ[RI], assuming stationary Poisson-like uncorrelated input spike trains and weak synaptic coupling J, is given by the Siegert equation (5) (Amit and Tsodyks, 1991; Amit and Brunel, 1997). Thus, in order to analyze the linear stability of this self-consistent solution, in general the derivatives with respect to both μ[RI] and σ[RI] need to be taken into account, cf. Equation (7).

We analyzed two different ways to drive the networks: in the first all neurons in the network were driven with the same excitatory Poisson-noise ν_X, independent on the recurrent coupling strength J. In this case the total input RI the neurons receive, i.e., the external input RI_x together with the recurrent input RI_s from the network, changes with increasing J, such that the mean μ[RI] decreases, while the standard deviation σ[RI] increases.

In the second input scenario, we corrected for the change in recurrent input contributions with changing J by administering compensating inhibitory and excitatory external Poisson input ν_Ix, ν_Ex. Thus, μ[RI] and σ[RI], and hence also ν_o are approximately constant with increasing J.

We find that in the first input scenario the critical J_c as observed in simulations deviates strongly from the one predicted by Equations (7),(9), which are derived from the self-consistent solution Equation (5) by linear response theory. In particular, for the mean-driven regime, where μ[RI] > V_thr, J^md_c becomes basically independent of the rate ν_o as well as of σ[RI], irrespective of its yet considerable magnitude. Instead, a very simple linearization derived from the noiseless f-I-curve Equation (10) applies. The J^md_c derived from this linearization is always considerably smaller than the prediction by Equations (7), (9). For intermediate cases where neurons are neither mean- nor strongly fluctuation-driven, the critical coupling strength lies in between the two predictions, and often closer to J^md_c.

Possible reasons for the failure of the prediction of Equations (7), (9) in the constant external drive scenario (parameterized by η, cf. section 2) are the break-down of the weak-coupling assumption, as well as the uncorrelated Poisson assumption for the input spike trains in the recurrent contribution. In line with earlier findings (Tetzlaff et al., 2012), we could not identify deviations with respect to coupling strength, however, we observed clear deviations from the Poisson assumption.

If spiking is Poissonian we expect the interspike intervals to be distributed exponentially, and the coefficient of variation (CV) of the interspike intervals (ISI) to be unity. A smaller CV indicates more regular, a larger one more irregular spiking. As can be seen from Figures 1A–C and 1G–I, activity is spatio-temporally structured, with locally clustered spiking that can induce significant pairwise correlations, and low coefficients of variation of the interspike intervals for small J (cf. also supplementary material section S1, Figures S1B,C). In effect, these deviations from uncorrelated Poisson spiking lead to deviations from the Gaussian white noise approximation underlying the derivation of Equations (5) and (9), especially for increasing J.

We observe that in particular the excitatory subpopulation tends to synchronize because of the highly recurrent local structure and positive feedback, while the inhibitory population actively desynchronizes itself due to negative feedback (see supplementary material section S1, Figures S3, S4, and Helias et al. (2013) for a related discussion). Thus, high amplitude volleys of excitatory input in a balanced background comprised of asynchronous inhibition and external Poisson-drive move the network in a regime that is better captured by the mean-driven low-noise approximation Equation (11).

Activity in the strongly fluctuation-driven regime, on the other hand, is much more asynchronous and irregular by construction, cf. Figures 1D–F with CV[ISI] of unity and larger (see supplementary material, section S1, Figures S1B,C). When in this scenario σ[RI] is very large, while μ[RI] is sufficiently subthreshold, the prediction for the critical coupling strength J^fd_c from Equations (7), (9) for linearity to break down indeed appears to be correct, or even seems to underestimate the actual onset of pattern formation.

This shallow transition around J^fd_c can be explained by slow noise-induced transitions between the two degenerate eigenvectors which grow quickest once pattern formation takes place. However, the fact that the CV is larger than unity will also lead to deviations from the prediction. In a series of papers Moreno et al. (2002); Renart et al. (2007) and Moreno-Bote et al. (2008) studied the impact of deviations of the Poisson-assumption of input spike trains on the firing rates of neurons. They studied the effect of positive and negative spike train auto- and cross-correlations parametrized by the Fano-factor F = σ²[counts]/μ[counts] of spike counts on the output firing rate of current-based LIF neurons in both the fluctuation- and the mean-driven regime. For renewal processes the CV of interspike intervals is directly related to the Fano-factor of counts by F = CV² (Cox, 1962). In Moreno et al. (2002); Renart et al. (2007); Moreno-Bote et al. (2008) it was shown that in the fluctuation-driven regime output rates are quite sensitive to input noise correlations (parameterized by the Fano factor), while in the mean-driven regime sensitivity is less pronounced.

In simulations of ring networks in the strongly fluctuation-driven case we observe that rates decrease with increasing correlations, while in otherwise equivalent random networks, rates increase (not shown). So the different recurrent input structure plays a crucial role for the effective input current and thus network dynamics properties. These differences are often, as here for ring and grid networks, a direct consequence of complex network topology, and it is thus important to understand how connectivity structure translates to single neuron activity, as well as collective network states. A more thorough analysis of such effects in the context of firing rate stability in random and spatially structured networks will be subject of subsequent research.

We assumed throughout the manuscript that the system has a negligible refractory time constant, small delays and instantaneous post-synaptic currents (δ-synapses). If these assumptions are dropped we observe differences especially for the mean-driven case: for larger delays d in the order of several milli-seconds, pronounced oscillations can occur that tend to synchronize the system locally on the order of the footprint (Kriener et al., 2009), or for very strong external drive also on a population-wide level (Brunel, 2000). A systematic study of the complex effects of delays on pattern formation in neural-field-type ring networks of spiking neurons was presented in Roxin et al. (2005). In the network we analyzed here, noisy delay oscillations mainly perturb developing patterns and can thus shift the onset of clear pattern formation to larger J_c.

Similar effects arise in the case of finite refractory time τ_ref. A larger τ_ref, especially in the presence of oscillations, i.e., highly coincidental input spiking, decreases the effective network input, in particular if τ_ref > d, and thus also shifts the input regime. This can be compensated for by explicitly taking into account τ_ref in the derivation of the noiseless linearization.

Finally, if synaptic time constants are finite, while total input currents stay the same, delay oscillations can be smoothened out, and the effect of absolute refractoriness is decreased. Increasing synaptic time constants can thus counteract the effect of finite τ_ref and d.

Pattern formation, such as activity bump formation or periodic patterns, is a well-studied phenomenon in ring and toroidal networks (see e.g., Ermentrout and Cowan, 1979a,b; Ben-Yishai et al., 1995; Usher et al., 1995; Bressloff and Coombes, 1998, 2000; Shriki et al., 2003; Marti and Rinzel, 2013). Earlier studies of periodic pattern formation in spiking neuron networks already showed that ring- and grid-networks of mean-driven neurons can undergo a Turing-bifurcation (Usher et al., 1995; Bressloff and Coombes, 1998, 2000). In particular, Usher et al. (1995) observed a dependence on the input level. For weaker, yet suprathreshold mean input current, activity patches diffuse chaotically across the spatial extend of the system, while for stronger inputs a stable spatial pattern develops that relates to the critical eigenmode. We observe similar dynamics in the networks considered here as well, if the system is in the intermediate regime. This can be understood as follows: if the noise component is relatively strong, it shifts the actual critical coupling strength J^md_c to higher values J^inter_c inbetween J^md_c and J^fd_c. If the system is thus effectively slightly sub-critical in terms of J^inter_c, input-noise can transiently excite eigenmodes close to effective criticality, due to Hebbian (Dayan and Abbott, 2001) or non-normal amplification (Murphy and Miller, 2009).

To compare our results to the classical work on neural field models on ring topologies (Ermentrout and Cowan, 1979a,b, 1980) we moreover studied a coarse-grained model that reduces the N-dimensional system to an effectively two-dimensional one that is structurally equivalent to the linearized two-population model presented, e.g., in (Ermentrout and Cowan, 1980) (see supplementary material section S2) which, however, allows for a direct quantitative mapping to the spiking network dynamics.

Because of the coarse-graining some of the connectivity details get lost and thus the critical eigenmode is not identical with that of the original ring network. Moreover, we observe the importance of excluding self-coupling in the coarse-grained model: if that is not excluded the resulting eigensystem looks dramatically different with nearly completely real-valued eigenvalues.

For the translation-invariant distribution of inhibitory neurons across a ring topology the eigensystem, and hence the respective linear stability, can be computed analytically. If inhibitory neurons are however randomly distributed, some parts of the ring will have a higher, others a lower density of inhibitory neurons, such that these networks form an ensemble of possible realizations, most of which are not invariant to translations. The resulting inhomogeneous rate per neuron in the subcritical regime can still be predicted well from the linear rate model. However, even when knowing the dominant eigenvector, it will usually not be periodic, and due to the inhomogeneous input per neuron different neurons will be at quite different working points. There is interference between the rate distribution and the dominant mode. Thus, the activity pattern that evolves in the supracritical regime is not straightforward to predict.

If the weights are distributed fully randomly across the whole ring topology irrespective of the identities of the neurons, i.e., in violation of Dale's principle, the resulting eigensystems look very similar to those of corresponding networks with random topologies. Only a few surviving singular eigenvalues serve to distinguish the underlying ring from the random topology. This directly translates to the spiking dynamics: the activity is much more asynchronous than that of Dale-conform networks and most of all pattern formation can not take place due to the effective lack of lateral inhibition.

If κ/N is kept constant when varying the network size N the critical coupling strength in both the mean- and the fluctuation-driven regime decreases. For example, for N = 10000 and κ = 1000 the critical coupling strength becomes J^md_c ≈ 0.2 mV for the mean-driven, and J^fd_c ≈ 0.32 mV for the fluctuation-driven regime. If the number of synapses per neuron κ is fixed when N increases, both J^md_c and J^fd_c hardly change because the maximal eigenvalue is approximately unaffected by network dilution.

We emphasize that the analysis presented here can be extended to two-dimensional grid networks in a straightforward manner, see supplementary material section S3. We conclude by remarking that the rate model Equation (13) is also a valuable tool to study the rate dynamics of spatially embedded spiking neuron networks with e.g., inhomogeneous neuron distribution and distance-dependent connectivity in the mean-driven regime, allowing for a treatment of more general networks than commonly studied random networks or neural field models. The presented model is in conclusion a further step in the efforts to find and understand appropriate descriptions of the high-dimensional spiking activity in structured networks.

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.