A single neuron is described by the standard leaky integrate-and-fire neuron model (Gerstner et al., 2014). The membrane voltage v evolves according to:

When v(t) > v_th, a spike is emitted and, after a refractory period of τ_R = 2 ms, the voltage is reset to v(t) = v_r. The parameter τ_m in Equation (1) is the membrane time constant, which may be different depending on whether the neuron belongs to the excitatory or the inhibitory population. The input current (scaled by the membrane resistance R), is denoted by RI.

The statistics inspected in this work are based on spike-trains, which are defined as sums of Delta functions

where t_i is the time instant of the ith spike. The instantaneous firing rate ν is the (generally time-dependent) average of the spike train, ν = 〈x(t)〉, and can be determined for a specific neuron within the network by an average over different runs with randomized initial conditions. We are only interested in the asynchronous state, in which ν does not depend on time. In practice, we often average the rate over the population (if appropriate, i.e., if the neurons are statistically equivalent) which is indicated by 〈.〉 (ensemble average) and over time:

For the calculation of spectral measures, we define the Fourier transform by

where T is our time window and is set in our simulations to T = 2 s if not mentioned otherwise. In all simulations we neglect a transient period of 1 s before extracting the statistics over the next T = 2 s. The power spectrum of a spike train is then defined by

where x~* is the complex conjugate of x~. We note that the power spectrum saturates for infinite frequency at the firing rate, limf→∞Sxx(f)=ν.

Two important statistical measures can be extracted from the power spectrum. The first is the Fano factor FF which is defined as the variance of the spike count N=∫0Tdtx(t) over its mean, an expression that can be related to the power spectrum at zero frequency:

The second statistical measure is the correlation time τ_c. Following Neiman et al. (2007) and Wieland et al. (2015), we consider the spike train's correlation function c(τ) = 〈x(t)x(t + τ)〉−〈x(t)〉〈x(t + τ)〉 (note that here 〈.〉 indicates a time average) and its continuous part ĉ(τ) = c(τ)−νδ(τ) to define the correlation time as an integral over the squared and normalized ĉ(τ)

an integral which in turn can be related to an integral over the power spectrum via the Parseval theorem on the right side.

Different network compositions are studied, many of which are based on the work of Brunel (2000), specifically on his Model B, a heterogeneous random network with fixed in-degree. In contrast to Brunel (2000), we use a larger number of neurons, i.e., an excitatory population size NE=105 instead of NE=104. Independently of the number of populations, there is always a mixture of excitatory to inhibitory neurons with a ratio of 4:1, i.e., N_I = γN_E where γ = 0.25. Therefore, the total network size is N = N_E + N_I.

The ℓth neuron from the network has the dynamics

The external input current to each neuron and its membrane time constant depend on the population it belongs to which is here indicated by index α. The ℓth neuron receives a fixed number of CℓE (CℓI) excitatory (inhibitory) randomly selected neurons connections from population α = {E, I}. The local input is described by:

where g_α is the ratio between average inhibitory and average excitatory synaptic weights, which depends via α on the target neuron (α ∈ E, I) The number of presynaptic neurons CℓE,I will be constant in some cases (fixed in-degree, in sections 2.3 and 2.4 as well as from 3.1 to 3.3) and random in others (as a consequence of a Erdős-Rényi topology in sections 2.5 and 3.4). The excitatory (inhibitory) input neurons are picked randomly from the network and the set of the neuron indexes is denoted by m_ℓ,k and n_ℓ,i. The synaptic coupling strength (also called synaptic weight or synaptic efficacy) will be either constant, J_ℓj = J, (in sections 2.3–3.4) or exponentially distributed with mean value 〈J_ℓj〉 = J (in section 3.4). We fix the transmission delay at τ_D = 1.5 ms unless otherwise indicated, and K(t) is a an optional synaptic filter. In most cases, the filter is not used, which means K(t) = δ(t). Otherwise, it is a simple exponential filter:

where θ(t) is the Heaviside function and τ_s is the synaptic filter time. Note that in the limit τ_s → 0, the case without synaptic filter is recovered. If not explicitly stated otherwise, we use the parameter values that are given in Table 1.

The iterative self-consistent scheme developed by Lerchner et al. (2006) and Dummer et al. (2014) is able to reproduce the single-spike-train power spectrum for homogeneous populations close to the balanced regime. In this procedure, in one generation a single neuron is stimulated with a colored noise over many trials, the power spectrum of its spike train is estimated, and using this spectrum and the output firing rate, a new surrogate colored Gaussian noise is generated which is used as the stimulus in the next generation. This procedure is repeated iteratively until the mean value and the spectrum of the driving noise matches in a self-consistent manner approximately the firing rate and the power spectrum of the resulting spike train. We present in detail the heterogeneous iterative self-consistent scheme in section 2.4, for further details of the homogeneous scheme we refer to Dummer et al. (2014).

The version of the scheme by Lerchner et al. (2006) and Dummer et al. (2014) is unable to reproduce self-consistently the statistics of single neurons in a recurrent network with strong relative inhibition g. More specifically, in cases where the inhibition is high, the scheme loses stability and the measured firing rate ν oscillates as a function of the generations (Dummer et al., 2014) (a numerical instability in the balanced case is reported in Lerchner et al., 2006, which is unrelated to the instability at strong recurrent inhibition). As a result, network regimes of low firing rate (such as those seen in cortex) cannot be captured.

Here, we propose a method to ensure convergence even with strong recurrent inhibition. We observed that the firing rate oscillations are around the target firing rate, therefore we can use the average firing rate over all past n generations as input to the next generation:

This procedure stabilizes the scheme, see Figure 3 for a numerical example. Note that averaging over a higher number of past generations can yield a faster convergence (cf. Figures 3B,C). The effect can be visualized using a similar approach as in Dummer et al. (2014): a map from the input rate to the output rate. We calculate the output rate from the input rate with the approximation for synaptically filtered white noise (Brunel and Sergi, 1998). The effect of the averaging over resulting output rate and input can be captured by the functions ν_out,0 = ν_out(ν_in), νout,1=νout(νin+νout(νin)2),…, which are shown in Figure 3D. These functions display an increasingly flatter shape in the dependence on the initial firing rate illustrating the stabilizing effect of the averaging.

The procedure of averaging the rate over the past generations will be used only in cases of unstable convergence. Typically, if excitation and inhibition are (nearly) balanced, the scheme is stable and we do not need to apply the averaging procedure.

The self-consistent scheme for a homogeneous population can be generalized in different ways. First of all, real networks consist of several types of neurons, that all differ with respect to their physiological parameters. A first important step is to distinguish between excitatory and inhibitory neurons not solely with respect to their postsynaptic effect but to endow inhibitory neurons also with other cellular parameters (membrane time constant, leak potential, mean input current) than excitatory cells. Generally, we distinguish between P^E excitatory and P^I inhibitory populations. In the self-consistent scheme each population is represented by one neuron.

Here, we generalize the iterative scheme to the case where the number of input neurons and the synaptic weights are not the same for the same type of neuron. We deal with a fixed connection probability (Erdős-Rényi topology) and an exponential distribution of synaptic weights, resulting in heterogeneity of firing rates and spike-train spectra. For simplicity, we use the same membrane time constant for inhibitory and excitatory neurons τ_E = τ_I = τ = 20 ms and do not use a synaptic filter. We represent the neurons of the network in the iterative scheme by a number of M neurons in each generation. The connection parameters of each neuron in each generation are randomly drawn from the assumed distributions. The input to each neuron in the subsequent generation is determined independently from the other neurons taking into account all power spectra of the previous generation.

In case of an Erdős-Rényi topology the numbers of presynaptic neurons are binomially distributed. Likewise, in the scheme we draw binomially distributed random numbers CkE and CkI as numbers of presynaptic neurons connected to neuron α with power spectrum S_k(f) from the previous generation. The mean of these random numbers is 〈CkE〉=CE/M and 〈CkI〉=CI/M respectively. In case of exponentially distributed synaptic weights we draw CkE excitatory synaptic weights Jk,lE and CkI inhibitory synaptic weights Jk,l′I for each neuron in the previous generation from the exponential distribution. The input current to the αth neuron is constructed with the mean input

and with the power spectrum of the effective input that is now taken to be

We assume here that E and I cells have the same spectrum if their connection parameters are equal. With given spectrum and mean input we calculate the voltage dynamics of neuron α as in section 2.4. To adapt the stabilization method described in section 2.3 to the case of a single heterogeneous population described by M single neurons, we sum a constant to all firing rates in the n-th generation, in order to set its mean firing rate to the average of all generations

where ν_k,n is the firing rate of the k-th neuron in generation n and ν^k,n is the transformed firing rate that should be used to calculate the input of the next generation. Note that the mean of the first two terms is zero. These ν^k,n are used in place of ν_k in Equation (16) to stabilize the scheme in case of strong inhibition.

For the simulation of these networks, we use the Brian spiking network simulator (Goodman and Brette, 2009).

We would like to start with results for the inhibition-dominated network (g > 4), in which firing rates are low. In this regime, the iterative scheme as proposed by Dummer et al. (2014) is highly unstable and we only obtain convergence with the averaging procedure described in section 2.3. To demonstrate that the averaging procedure works in such a situation, we consider in Figure 4 the network studied by Ostojic (2014) who found two contrasting asynchronous states when varying the synaptic strength J.

In Figure 4A spike-train spectra for strong recurrent inhibition (g = 5) for different values of J and different network sizes are shown. The power spectra of these network simulations are close to those of the iterative scheme in most cases. For weak coupling, the agreement between spectra is always good; discrepancies for large J become smaller with increasing network size because cross-correlations become less important in this limit. An additional reason for discrepancy is that the Gaussian approximation becomes less accurate for strong synaptic strength. The change in spike-train power spectra upon increase of the synaptic coupling does not hinge on the specific nature of the subthreshold function in the IF model. If we replace the leaky IF model by an exponential IF model (Fourcaud-Trocmé et al., 2003) in the network and in the recurrent scheme, we observe a similarly drastic change in low-frequency power if the synaptic strength is doubled (inset of Figure 4A). Also for this single-neuron model the agreement between spectra from the network and from the self-consistent scheme is fairly good.

When the coupling strength J increases, the firing rate first decreases and then increases (Figure 4D). More interestingly, with increasing coupling we see a transition from Poisson-like irregular firing (Figure 4B) to bursty firing of single neurons (Figure 4C), i.e., periods of strong firing are separated by pauses. In the latter state one can observe a broad distribution of spike counts, and that is why this state has been referred to as heterogeneous asynchronous state (Ostojic, 2014). In terms of the power spectrum this transition becomes manifest as an amplification at low frequencies (Wieland et al., 2015); correspondingly the Fano factor increases (Figure 4F). Together with the minimum in the correlation time (Figure 4E) (attained at a coupling where the Fano factor is about unity), our results confirm that the transition described by Ostojic (2014) in the inhibition-dominated regime is essentially the same as the one observed and explained by Wieland et al. (2015) for the balanced case g = 4.

In summary, the results in Figure 4 indicate that the emergence of a new heterogeneous asynchronous state for strong synaptic coupling can be explained only using the properties of a single neuron and the self-consistency condition, here demonstrated by our iterative single-neuron scheme.

We now investigate the effect of a finite synaptic filter, Equation (10). Not surprisingly, a pronounced synaptic filter (large τ_s) leads to a long time scale in the network dynamics, as revealed by the increased power at low frequencies (Figure 5). The synaptic filter, Equation (10) is scaled such that the total charge per input spike remains constant. Therefore, an increased time constant for the exponential decay renders the postsynaptic response smaller in amplitude and longer in duration. This longer duration of postsynaptic responses extends the range of temporal correlations in the input to the neuron, which in turn causes the slow fluctuations in the neuron's activity. The resulting power spectrum (Figure 5A), especially for long synaptic time constant, looks similar to that of a colored-noise driven perfect IF model (see Figure 9 in Middleton et al., 2003). We emphasize that the emergence of the slow time scale is here imposed by the long-lasting synaptic filter which is in marked contrast to the network amplification of slow fluctuations for strong synaptic coupling discussed before in Figure 4.

We also verified that the synaptic filter does not change qualitatively the emergence of slow fluctuations for strong coupling (i.e., the heterogeneous asynchronous state discussed above). Using τ_s = 10 ms we still see the characteristic strong increase in Fano factor (Figure 5B) and a minimum in the correlation time (Figure 5C).

In the following, we return to the limit of instantaneous synapses τ_s → 0, i.e., K(t) = δ(t), and introduce different parameter values for excitatory neurons (E-cells) and inhibitory neurons (I-cells). First of all, in order to test whether the applicability of the scheme hinges on the exact value of crucial parameters, we choose a small change of the membrane time constant between E and I-cells: τ_I = 19 ms and τ_E = 20 ms. Secondly, we make the relative strength of recurrent inhibition, g_E and g_I, different for the two populations in order to see whether the generalized iterative scheme with two neurons can cope with this heterogeneous situation.

In Figure 6 we show power spectra obtained from simulations of the recurrent network and of the iterative scheme for different combinations (g_E,g_I). In Figure 6A the two populations are statistically rather different with an E-cell firing rate of ν_E = 3.2 Hz whereas I-cells fire at ν_I = 9.7 Hz. Both spectra are well reproduced by the iterative scheme and show a “green” shape (in the colored noise lingo, this is white minus red noise). That means, the spectra exhibit a dip at low frequencies, but this is much more pronounced for the I-cells. Even when we increase the difference in recurrent inhibition and the two types of neurons fire at lower frequencies of ν_E = 0.1 Hz and ν_I = 7.4 Hz, the agreement of the spectra from the iterative scheme and from the network simulations is excellent (Figure 6C). If we choose the relative recurrent inhibition to be the same, the neural dynamics differ only by the small difference of the membrane time constants, which does neither cause significant differences in the firing rates (ν_E = 128.9 Hz and ν_I = 129.9 Hz) nor in the shape of the power spectra (cf. Figure 6B).

In order to explore the quality of the approximation systematically, we evaluated the discrepancy using the relative integrated error

where f_cut = 2ν_I (we use the inhibitory firing rate because it is usually higher).

In our scheme the assumption of weak cross-correlations among neurons in the network is crucial - indeed we assume an infinitely sparse system that is in a perfectly asynchronous state. This is, of course, a somewhat artificial limit and thus it is interesting how, for a fixed number of connections (about 10³), the squared deviation as well as important statistics such as the Fano factor depend on the system size. In Figure 7 this dependence is illustrated for the case where g_E = 4.2 and g_I = 4.0, the same parameters as in Figure 6A. For the chosen connectivity, a minimal number of N_E = 20, 000 E-cells seems to be required to reach a good approximation (relative error below 1% for both E and I cells) with the self-consistent scheme. This plot illustrates that although sparsity is an important assumption for the self-consistent determination of spike-train power spectra, it does not lead to the necessity to consider exorbitantly large networks.

In principle, the proposed iterative scheme is applicable to any number of populations. As long as the resulting activity is sufficiently asynchronous (implying weak cross-correlations) and the synaptic strength is not excessively large (needed for the Gaussian approximation), the iterative scheme should converge to a self-consistent result. Here we demonstrate that the extended scheme also works for networks with more than two populations and study two cases: a network with three distinct populations and a modular network.

An example of three populations is given by a combination of one excitatory and two inhibitory populations (Figure 8A), biologically inspired by a cortical network with excitatory regular spiking neurons (RS), inhibitory fast-spiking (FS), and low-threshold spiking (LTS) neurons (see Izhikevich, 2003; Tomov et al., 2014 and references therein). Their firing rates are ordered such that ν_FS > ν_LTS > ν_RS. This heterogeneous situation is achieved by changing both membrane time constants, which are chosen to be τ_FS = 21 ms, τ_LTS = 20 ms, and τ_RS = 19 ms, and making one of the synaptic weights in the network (connecting RS neurons to FS neurons) 1.4 times stronger (indicated by the thick arrow in Figure 8A, left). This setting illustrates how heterogeneity of connectivity and membrane time constants shape the power spectra statistics.

The resulting spectra are well-captured by the iterative scheme; they all display the effect of neural refractoriness by the dip at low frequencies (Bair et al., 1994; Franklin and Bair, 1995) but to a different degree. The dip is most pronounced for the fast spiking neurons; the regular spiking neurons fire with a statistics that is closest to a Poisson process with a flat power spectrum.

According to a common view, the cortex possesses a modular structure (Boucsein et al., 2011; Tomov et al., 2014, 2016), a feature that we take into account in the next setup. We consider two different modules as shown in Figure 8B. The two modules are equal to each other with respect to the population size and each consists of an E-I network with N_E = 100, 000 and membrane time constants τ = 20 ms, requiring the simulation of two neurons in total in the self-consistent scheme. For module 1 and 2 we choose J = 0.1 mV and J = 0.27 mV, respectively. With this choice, the two modules operate in different regimes: module 1 in a fast-fluctuation mode with low Fano factor and peaked power spectrum, module 2 in a regime of dominating slow fluctuations (cf. section 3.1). In module 2, we rewired 50% of the connections so that they come from module 1, i.e., it receives 0.5C_E excitatory inputs and 0.5C_I inhibitory inputs from module 1. This results in a highly heterogeneous situation which is reflected in the power spectrum of module 2: in contrast to the behavior observed in Figure 4, the power spectrum of module 2 contains an additional hump around 80 Hz. The power spectra of all different neurons in this setup are well represented by the iterative scheme, cf. Figure 8B.

The result in Figure 8B demonstrates that the iterative scheme can capture complex situations involving the interaction among different modules. The simulated network contained in total 250,000 neurons and the iterative scheme reproduced the single neuron correlation statistics with high accuracy using only two neurons.

Here we use the iterative scheme to represent network dynamics for distributed connection parameters in a single population. As a first example, we consider a network in which connections are made with a fixed probability (Erdos and Rényi, 1960) instead of a fixed in-degree (as in all previous examples). This Erdős-Rényi topology yields a binomially distributed in-degree which is certainly a more realistic scenario than the fixed in-degree. Furthermore, a distributed in-degree leads to a distribution of firing rates, as seen in cortical networks (Griffith and Horn, 1966; Koch and Fuster, 1989; Shafi et al., 2007; Hromádka et al., 2008; O'Connor et al., 2010; Roxin et al., 2011). In addition to a random number of connections, in cortical networks also the synaptic weights are not fixed but follow a long-tailed distribution (Song et al., 2005; Gilson and Fukai, 2011). This feature can be approximated in the recurrent network model by drawing the weights from an exponential distribution.

In this version of the iterative scheme, we have M neurons that represent independent samples of the distributions of excitatory and inhibitory in-degrees and weights. Consequently, the firing rates and power spectra of the M neurons will differ and reflect the heterogeneity of these measures in the recurrent network (for further details, see section 2.5).

The first striking feature for a randomized connectivity is a broad distribution of firing rates. The stronger source for the rate variability seems to origin in the random number of connections (the case shown in Figure 9A); an exponential distribution of the synaptic strength makes the histograms only somewhat broader (case with additional weight variability, shown in Figure 9B). The histograms obtained from our iterative scheme with a modest number of M = 50 representative neurons (data are collected over several generations to improve the sampling) agree well with those from the network simulations. Note that analytical expressions for self-consistent distributions of the firing rate have been found previously using the diffusion approximation (Amit and Brunel, 1997; Roxin et al., 2011).

In Figure 10 we show the mean Figure 10A and the standard deviation Figure 10B of the firing rate histograms of the iterative scheme as a function of the generation. Because the number M of neurons per generation (i.e., the number of samples) is finite we consider the standard error σs=σ/M-1 represented by the colored areas in Figure 10A, where σ is the standard deviation of the firing rates measured from the network. The mean firing rate and standard deviation of the iterative scheme fluctuate around the network values within the expected standard error (colored areas), which indicates that the scheme works even for rather limited numbers of neurons M ≪ N. For the considered parameter sets, we observe fast convergence within five iterations for all used numbers of neurons except in the case of strong recurrent inhibition with g = 6 for which we observe convergence only after the 10th generation.

Power spectra are different for all neurons in the network and depend most strongly on the mean input that the respective cell receives. In the network we group neurons with similar firing rates (within a 1 Hz interval) and average their spectra within each group. We compare these spectra to those resulting from the self-consistent scheme having firing rates in the same interval (Figure 10C); we find a good agreement for a modest number of neurons M = 50 used in the scheme. Note that averaged spectra are very different from each other and from the average over all neurons in the network (Figure 10D). Even for a low number M = 10, the iterative scheme reproduces this average spectrum well.

If we use in addition to the random number of connections also randomly distributed synaptic weights, the overall picture does not change qualitatively (Figure 11). Compared to the case of equal weights, the variance in the firing rates goes up by about 50% (Figure 11B) whereas the mean firing rate increases only slightly (Figure 11A). Power spectra look similar to the previous case and the agreement between scheme and network simulation is again good, except for the case of strong inhibition g = 6, for which we find a discrepancy between the averaged power spectra for the scheme and the network (Figure 11D) even for a large M (similarly for the single neuron spectra in Figure 11C). Due to the very low firing rates and some high synaptic weights, the Gaussian approximation becomes inaccurate. We calculated the relative error defined in Equation (19) for the averaged power spectra for M = 10 up to M = 100, using f_cut = 2〈ν〉. All observed relative errors were smaller than five percent.

The results of this section demonstrate how a relatively small number of neurons in the self-consistent scheme can be used to capture the first- and second-order statistics of neurons in a network with a variable number of inputs and synaptic weights.

The iterative scheme assumes neuronal input that obeys a Gaussian distribution. This assumption is only accurate for a large number of pre-synaptic spikes with small amplitude. A low firing rate or an increase of (some) synaptic weights (as observed for stronger recurrent inhibition, i.e., larger values of g) shifts the statistics from approximately Gaussian to a distribution with pronounced non-Gaussian features (skewed and with fat tails). In fact, even the most basic statistics of the output, the firing rate, can systematically deviate from the Gaussian setting if synaptic amplitudes are sufficiently large (Richardson and Swarbrick, 2010). We can thus expect that at large values of g the iterative scheme (that assumes Gaussianity) may provide power spectra that deviate from those measured in the network; indeed, this is what we observed above. Here we want to illustrate the relation between the deviations of power spectra in the network and for the isolated neuron on the one side and non-Gaussian features of the network input statistics on the other side.

We simulate the heterogeneous network from section 3.4 with exponentially distributed weights and either g = 4 (as in Figure 11, middle column) or g = 6 (as in Figure 11, right column). We pick out 1,000 representative neurons randomly. For each of these neurons we measure the input current spectrum and the output spike-train spectrum. Input currents are integrated over a time bin of 2 ms (10% of the membrane time constant), in which the membrane voltage does not change too drastically but we can expect to collect many input spikes and thus Gaussian statistics. Clearly, for a much smaller time bin of the order of [(CE+CI)ν̄]-1≈ 0.01 ms (where ν̄ is the average firing rate of the presynaptic neurons), one single spike at most would typically fall into one time bin. Thus, on a very fine temporal scale, the input has always a highly non-Gaussian statistics. Nevertheless, the Gaussian approximation works well if the temporally integrated input noise has approximately Gaussian statistics for an effective time step, in which the voltage does not change much. As the deterministic drift causes already a change by several mV over a time of 2 ms, our choice of the effective time step is conservative in the sense that non-Gaussian features of the input in this time step will certainly become noticeable in the output statistics of the driven cell.

For each of the 1,000 selected neurons, we first measure the skewness and excess kurtosis of its input current to quantify non-Gaussianity. Then, we generate a Gaussian noise with the same power spectrum as the input current to this specific neuron and use it to drive an isolated LIF neuron. Finally, we calculate the relative error, Equation (19), between the output power spectrum of this isolated neuron and the corresponding neuron from the recurrent network. Note that for this comparison, we do not iterate the self-consistent scheme but rather use a Gaussian version of the (neuron-specific) network input noise for the isolated cell. In this way, we isolate non-Gaussianity of the input noise as one potential source of deviation between network spectra and spectra from the self-consistent scheme.

Plots of the relative error vs. skewness and excess kurtosis are presented in Figures 12A,B, respectively. For g = 4 (balanced case, moderate inhibitory synaptic amplitudes), the relative error is generally small and skewness and excess kurtosis deviate only little from zero. In contrast, for stronger inhibitory synaptic weigths g = 6 (inhibition-dominated recurrent feedback), the input current has a pronounced negative skew and fat tails because of large inhibitory amplitudes, resulting in stronger deviations between the spectra of isolated neurons and the corresponding network neurons. Consequently, we find a correlation between non-Gaussianity and relative error (cf. the Pearson correlation coefficients given atop of the figure). To gauge the relative error, we need a lower limit. To this end, we consider the deviation between two measurements of the same power spectrum for a given neuron of the network. The average of this lower bound over the 1,000 selected neurons is shown in Figure 12 for comparison. It is close to the mean error for g = 4 but is significantly lower than the mean error for g = 6, and it barely depends on the skewness and kurtosis, again confirming the relation between non-Gaussianiaty of the input noise and the relative error between isolated neuron and network neuron.

For both values of g we obtain a high relative error for exceptionally large firing rates (>100Hz) (because the time step of 2 ms becomes too large in comparison to the mean ISI) or exceptionally low firing rate (because spectra are very noisy in this case). If we consider only a range of moderate firing rates (Figure 13) then the correlation between the deviation from Gaussianity and the discrepancy between network and isolated becomes even clearer (cf. Pearson correlations coefficients stated atop of the figure).

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.