We compute bounds on network spike-response entropy that give novel answers. In particular, we show that the noise entropy of multi-cell spike responses is at least an order of magnitude lower than would be naively extrapolated from from single-cell measurements, under the assumption that spike variability is independent from cell to cell. The direction of the comparison between noise entropy of single cell and multi-cell spike responses agrees with intuition provided by the shape of the Lyapunov spectrum, which indicates time-dependent chaotic attractors of lower dimension than phase space. Thus, the phase space dynamics of each neuron are not independent. What we quantify explicitly is the order-of-magnitude size of the effect, as it is manifested in the binary *spiking* outputs of the system—a fact which might seem especially striking given that pairs of spike trains appear to be very weakly correlated on average.

If one considers the level of noise entropy as an indicator of potential information contained in spike patterns, we show that balanced networks may be able to encode inputs stimuli using spike timing if these inputs contain strong enough temporal structure. This mechanism takes root in the complex noiseinteractions that chaos induces between neurons. The extensive nature of this phenomenon suggests that this mechanism is scalable with network size. Moreover, the strong dependence of entropy on the input signal's mean and variance indicate that

a network can operate in different "regimes" modulating the repeatability of spike patterns. This is in addition to known advantages of balanced networks, such as efficiently tracking changes in common, mean inputs with firing rates (van Vreeswijk and Sompolinsky, 1998)—which may encode coarser statistics about inputs at the population level.

parametrization. Markers: square: η = −1, ε = 0.69, star: η = −0.5, ε = 0.5, circle: η = 0.07, ε = 0. **(B)** Lyapunov spectra along contour

To formalize these notions, future work could seek to compute the mutual information between an input ensemble and a system's response. In order to estimate this quantity, one needs to compute the *total entropy* (Rieke et al., 1996) of spike patterns in addition to the noise entropy computed in this paper– which captures how many distinct spike outputs can be produced by the network, for any input *I*. This quantity can be thought of as noise entropy marginalized over the set of possible inputs. Estimating the total entropy in large networks is a difficult problem since it depends on the evolution of ensembles of trajectories driven by ensembles of inputs. In other words, one needs to capture the entropy of trajectories when system (1) is treated as a stochastic differential equation rather than a RDS, a distinction that introduces a variety of challenges.

Our results complement prior work on the behavior of sparse, balanced networks in the large *N* limit. Seminal results use meanfield approaches (e.g., van Vreeswijk and Sompolinsky, 1998), deriving successful estimates of population activity statistics such as the mean and the variance of firing rates. In this approach, self consistent equations are derived for representative single cells based on the assumption that, when *N* is sufficiently large and *k*/*N* is sufficiently small, the inputs to each neuron in the network can be approximated by independent gaussian noise. In contrast, we derive estimates for the impact of correlations among these individual cells. Interestingly, in both the classical and the present work, noise entropy scales extensively with *N*; here, the predicted rate of scaling would be lower, as even weak correlations between cells combine to create statistical dependencies—especially when network activity is conditioned on an input.

and *HKS* extrapolated from a network with *N* = 500, as are quantities

Finally, we expect that the *HKS* bound can be adapted to other neuron models, provided a state space partition linking dynamics to spike patterns can be derived. This could prove to be a powerful tool to investigate stimulus encoding as a function of many network attributes, such as spike-generating dynamics, connectivity, learning rules and input correlations.

## **FUNDING**

from all other panels.

This work was supported in part by an NIH Training Grant from University of Washington's Center for Computational Neuroscience (5R90DA03346103, 5T90DA03243602), a Bernstein Fellowship from the Bernstein Center for Computational Neuroscience, a postdoctoral scholarship from the *Fonds de Recherche du Québec*, the Burroughs Wellcome Fund for Scientific Interfaces, the NSF under grant DMS CAREER-1056125 and NSERC Discovery and CIHR operating grants. Numerical simulations were performed on NSF's XSEDE supercomputing platform.

## **ACKNOWLEDGMENTS**

The authors thank Fred Wolf, Yu Hu, and Kevin K. Lin for helpful insights. We also thank two anonymous reviewers for comments and suggestions that improved the manuscript.

## **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fncom. 2014.00123/abstract

## **REFERENCES**

**Conflict of Interest Statement:** 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.

*Received: 16 May 2014; accepted: 11 September 2014; published online: 02 October 2014.*

*Citation: Lajoie G, Thivierge J-P and Shea-Brown E (2014) Structured chaos shapes spike-response noise entropy in balanced neural networks. Front. Comput. Neurosci. 8:123. doi: 10.3389/fncom.2014.00123*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2014 Lajoie, Thivierge and Shea-Brown. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Commentary on Structured chaos shapes spike-response noise entropy in balanced neural networks, by Lajoie, Thivierge, and Shea-Brown

## *Peter J. Thomas\**

*Department of Mathematics, Applied Mathematics and Statistics, Case Western Reserve University, Cleveland, OH, USA \*Correspondence<EMAIL_ADDRESS>
#### *Edited by:*

*David Hansel, University of Paris, France*

#### *Reviewed by:*

*Demian Battaglia, Max Planck Institute for Dynamics and Self-Organization, Germany Guillaume Lajoie, Max Planck Institute for Dynamics and Self-Organization, Germany*

**Keywords: spike-response noise entroy, balanced neural networks, information theory, chaos, mutual information**

#### **A commentary on**

**Structured chaos shapes spike-response noise entropy in balanced neural networks**

*by Lajoie, G., Thivierge, J.-P., and Shea-Brown, E. (2014). Front. Comput. Neurosci. 8:123. doi: 10.3389/fncom.2014.00123*

Animals with nervous systems generate complex adaptive behaviors in part through the computational capabilities arising from large networks of interconnected neurons in their brains (Churchland and Sejnowski, 1992). Although a full description of the nervous system would take into account the interactions of central circuits with sensory and motor systems (Chiel and Beer, 1997) it is more common to consider central circuitry in isolation. The individual nerve cells and synaptic junctions that comprise biological neural networks are spatially extended structures with fundamentally stochastic dynamics on a range of spatial and temporal scales (Andersen et al., 2006; Carnevale and Hines, 2006). Nevertheless, much progress has been made in understanding the repertoire of neural behavior through simplified deterministic one dimensional "phase" models such as the Ermentrout-Kopell canonical model (Ermentrout, 1996; Brown et al., 2004; Ermentrout, 2008)1 .

Even if we restrict attention to isolated networks of deterministic, instantaneously coupled phase models, we confront significant challenges. The behavior of such networks can be chaotic, as evidenced by the divergence of nearby trajectories (positive Lyapunov exponents). If we consider such a "chaotic network" driven by a collection of input signals, it is natural to ask how the intrinsic variability related to the chaotic dynamics impacts the networks' computational capabilities. It is equally natural to view the system as a communications channel. With the input signals drawn from some specified ensemble, and the output taken as the spike trains of (some or all of) the neurons, the mutual information between the input and output ensembles would be of great interest. However, this quantity is difficult to obtain, either analytically or numerically.

In Lajoie et al. (2014), the authors further the analysis of information processing in chaotic deterministic networks by formulating a computationally tractable upper bound on the spike-train noise entropy, building on Monteforte and Wolf (2010) and Lajoie et al. (2013). They study a network of deterministic canonical Ermentrout-Kopell "theta" neurons (Ermentrout and Kopell, 1986) with an *ad-hoc* interaction function. The network connectivity is fixed, sparse and random. Each neuron is driven by a quenched white noise injected current input of the form *Ii*(*t*) = η + *dWi*,*t*/*dt*. As the authors (and others) have shown previously, the spontaneous activity (i.e., with = 0) in this class of networks exhibits chaotic behavior. It has been observed that applying an input to such networks (i.e., setting  > 0) can reduce the apparent irregularity of the spike train ensemble. The *spike train entropy* quantifies this reduction in variability; the authors obtain an upper bound on this quantity through a state space partitioning construction that takes advantage of the Kolmogorov-Sinai entropy, which is given in turn by the Lyapunov spectrum, which the authors estimate numerically. They show convincingly that the KS entropy of the spike trains is roughly an order of magnitude smaller than what one would expect from a naive estimate based on the single-cell noise entropy. Their results help make rigorous the observation that the application of a driving stimulus reduces the variability of the resulting spike trains, although the networks remain chaotic.

While this result is a substantive contribution, it is still some steps removed from telling us the mutual information *I*(*X* : *Y*) = *H*(*Y*) − *H*(*Y*|*X*) between an ensemble of inputs, *X*, and the corresponding ensemble of outputs, *Y*. The authors' result gives a bound on *H*(*Y*|*x*) for a specific realization of the frozen noise inputs *x* ∈ *X*. Because the system is ergodic, this estimate applies as well to the mean entropy *H*(*Y*|*X*) [as discussed in Lajoie et al. (2013)]. However, as the authors point out, one cannot replace the entropy *H*(*Y*) with *H*(*Y*|0), the entropy when the input fluctuations are switched off, since (as they convincingly demonstrate) turning on the input ( > 0) significantly changes the entropy. The entropy that would be needed for calculating the mutual information would be the spike train entropy

1As pointed out in Wolf et al. (2014), the choice of 1D model can strongly influence the entropy production properties of the resulting network.

for the ensemble unconditioned on a specific input—but with an ensemble of different white noises all with fixed > 0. It would be very interesting if one could investigate how *I*(*X* : *Y*) varied as a function of ; for instance, whether the mutual information changes smoothly or whether there is evidence for some kind of information processing phase transition. The authors' contribution provides a valuable step along the way to a deeper understanding of the impact of chaotic dynamics on computations in deterministic neural networks.

## **REFERENCES**

Wolf, F., Engelken, R., Puelma-Touzel, M., Weidinger, J. D. F., and Neef, A. (2014). Dynamical models of cortical circuits. *Curr. Opin. Neurobiol.* 25, 228–236. doi: 10.1016/j.conb.2014. 01.017

**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 02 October 2014; accepted: 08 February 2015; published online: 10 March 2015.*

*Citation: Thomas PJ (2015) Commentary on Structured chaos shapes spike-response noise entropy in balanced neural networks, by Lajoie, Thivierge, and Shea-Brown. Front. Comput. Neurosci. 9:23. doi: 10.3389/fncom. 2015.00023*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2015 Thomas. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks

## *Jorge F. Mejias 1,2\* and André Longtin2,3*

*1 Center for Neural Science, New York University, New York, NY, USA*

*2 Department of Physics, University of Ottawa, Ottawa, ON, Canada*

*3 Department of Cellular and Molecular Medicine, University of Ottawa, Ottawa, ON, Canada*

#### *Edited by:*

*Mark D. McDonnell, University of South Australia, Australia*

#### *Reviewed by:*

*Alberto Mazzoni, Scuola Superiore Sant'Anna, Italy Fabiano Baroni, Monash University, Australia*

#### *\*Correspondence:*

*Jorge F. Mejias, Center for Neural Science, New York University, 4 Washington Place, New York, NY 10003, USA e-mail<EMAIL_ADDRESS>
Recent experimental and theoretical studies have highlighted the importance of cell-to-cell differences in the dynamics and functions of neural networks, such as in different types of neural coding or synchronization. It is still not known, however, how neural heterogeneity can affect cortical computations, or impact the dynamics of typical cortical circuits constituted of sparse excitatory and inhibitory networks. In this work, we analytically and numerically study the dynamics of a typical cortical circuit with a certain level of neural heterogeneity. Our circuit includes realistic features found in real cortical populations, such as network sparseness, excitatory, and inhibitory subpopulations of neurons, and different cell-to-cell heterogeneities for each type of population in the system. We find highly differentiated roles for heterogeneity, depending on the subpopulation in which it is found. In particular, while heterogeneity among excitatory neurons non-linearly increases the mean firing rate and linearizes the f-I curves, heterogeneity among inhibitory neurons may decrease the network activity level and induces divisive gain effects in the f-I curves of the excitatory cells, providing an effective gain control mechanism to influence information flow. In addition, we compute the conditions for stability of the network activity, finding that the synchronization onset is robust to inhibitory heterogeneity, but it shifts to lower input levels for higher excitatory heterogeneity. Finally, we provide an extension of recently reported heterogeneity-induced mechanisms for signal detection under rate coding, and we explore the validity of our findings when multiple sources of heterogeneity are present. These results allow for a detailed characterization of the role of neural heterogeneity in asynchronous cortical networks.

**Keywords: heterogeneity, asynchronous state, gain control, mean-field, cortical networks, signal detection**

## **1. INTRODUCTION**

Mathematical models of neurons and neural circuits have become, in the last couple of decades, a highly valuable tool to analyze and understand real neural systems, from single cell to behavior. Models are commonly used to test hypotheses or to support experimental observations, and their potential usefulness increases as their predictions are refined to account for the actual behavior of neurons (Gerstner and Naud, 2009). While it is not uncommon to see a high level of biophysical detail in singleneuron models, most of these details are usually neglected when modeling larger systems, such as neural circuits of thousands of neurons, for the sake of simplicity.

A particularly interesting case is the natural intrinsic variability found in the biophysical properties of neurons, which is averaged out in most theoretical and computational modeling studies. Real neural systems display a significant level of cell-tocell diversity at the neuron level, even among same-class neurons, as well as other differences at the subcellular or synaptic level (Bannister and Larkman, 1995a,b; Reyes et al., 1998; Hausser and Mel, 2003; Jinno et al., 2007). Experimental observations indicate that this type of structural heterogeneity has non-trivial effects on several neural information processing mechanisms. For instance, neural heterogeneity has been shown to have an impact in burst coding *in vivo* (Avila-Akerberg et al., 2010) and in envelope coding and non-linear responsiveness of the electroreceptors of weakly electric fish (Savard et al., 2011). The presence of a certain level of heterogeneity at the cell-to-cell level has also been recently reported to have a beneficial role for population coding (Marsat and Maler, 2010; Tripathy et al., 2013), and it can also induce the decorrelation of neuronal firing and the optimization of information content (Padmanabhan and Urban, 2010; Angelo et al., 2012; Urban and Tripathy, 2012). These experimental observations can not be explained by neural circuit models where, for instance, any given pyramidal neuron is perfectly identical to all the other pyramidal neurons in the system. Models which take into account the intrinsic heterogeneity of neural systems are, therefore, necessary to understand neural coding.

In response to this increasing body of evidence, a significant number of theoretical and computational studies, especially in the last years, have contributed to explaining the properties and dynamics of networks of heterogeneous neurons. In particular, the role of heterogeneity on synchronization has been extensively studied (Golomb and Rinzel, 1993; White et al., 1998; Neltner et al., 2000; Golomb et al., 2001; Denker et al., 2004; Talathi et al., 2008, 2009; Luccioli and Politi, 2010; Olmi et al., 2010; Brette, 2012; Mejias and Longtin, 2012). More recently, the effect of neural heterogeneities on neuronal correlations (Chelaru and Dragoi, 2008; Yim et al., 2013), detection of weak signals (Tessone et al., 2006; Perez et al., 2010) and different types of neural coding (Chelaru and Dragoi, 2008; Savard et al., 2011; Mejias and Longtin, 2012; Hunsberger et al., 2014) have drawn special attention as well. Novel approaches and mean-field approximations to tackle the problem of heterogeneity have also been recently proposed (Nicola and Campbell, 2013; Yim et al., 2013). These studies, however, typically focus on one type of cell (such as pyramidal neurons) and consider the presence of heterogeneity on this specific population. The possible—and potentially relevant interplay between populations of different cell types, each one of them presenting its own heterogeneity level, has remained a goal for future work and has not been fully addressed yet.

In this work, we analytically and numerically study the properties of a typical cortical circuit with cell-specific heterogeneity levels. Our basic circuit is constituted by a population of excitatory (i.e., pyramidal) neurons and a population of inhibitory neurons (i.e., interneurons), connected in a sparse manner. Both the excitatory and the inhibitory populations have their own independent level of intrinsic heterogeneity. This allows us to quantitatively study the effect of heterogeneity of a given population (excitatory, or inhibitory, or both) and to characterize its impact on the whole network dynamics. Our results indicate highly differentiated roles for heterogeneity, depending on the population in which it is introduced. In particular, heterogeneity among excitatory neurons (which we call here *excitatory heterogeneity*) non-linearly increases the mean firing rate of the whole network and linearizes the input/output f-I curves. On the other hand, heterogeneity among inhibitory neurons (*inhibitory heterogeneity*) may decrease the network activity level and induce divisive gain effects in the f-I curves of the excitatory population, providing an effective gain control mechanism to influence the flow of information across the network. We also compute the conditions for stability of the network activity and provide an extension of the recently described heterogeneityinduced mechanism which optimizes information transmission under rate coding (Mejias and Longtin, 2012). Our novel meanfield approach extends our previous theoretical results for fully connected excitatory neurons (Mejias and Longtin, 2012) to cortical-like sparsely connected networks of heterogeneous excitatory and inhibitory cells, providing a strong analytical tool to characterize the role of neural heterogeneity in cortical networks.

#### **2. MATERIALS AND METHODS**

We consider a sparsely connected network of *N* integrate-andfire neurons (see **Figure 1**), where any two given neurons are

unidirectionally connected with a probability - (the average number of synapses onto a given neuron is then *K* ≡ -*N*). A subset of this population is constituted by γ *N* excitatory neurons, while the remaining (1 − γ )*N* neurons in the network are inhibitory. A given neuron *i* is governed by the dynamics

$$
\tau\_m \frac{dV\_i(t)}{dt} = -V\_i(t) + RI\_i^{\text{ext}}(t) + RI\_i^{\text{net}}(t), \tag{1}
$$

where τ*m* is the neuron membrane time constant, *Vi* is the membrane potential of the *i* − *th* neuron in the network, *R* is the membrane resistance, and *Iext i* , *Inet i* are the external and recurrent input to the *i* − *th* neuron, respectively. Each neuron *i* is assumed to fire an action potential (AP) every time *Vi* reaches a certain firing threshold, and after that the membrane potential is reset to *Vr* for a time period τ*ref* . The external and recurrent input to the *i* − *th* neuron are given, respectively, by

$$RI\_i^{\text{ext}}(t) = \mu\_i + \sigma \sqrt{\tau\_m} \xi\_i(t), \tag{2}$$

$$RI\_i^{net}(t) = \mathfrak{r}\_m \sum\_j \sum\_k f\_{ij} \,\delta(t - t\_j^k), \tag{3}$$

where μ*i* is a constant input, ξ*i*(*t*) is a Gaussian white noise of zero mean and unitary variance, σ is the noise intensity, *Jij* is the coupling strength of the synapse from neuron *j* to neuron *i* (considered zero if there is not such a synapse between both neurons), and the *k* − *th* spike from neuron *j* arrives at neuron *i* at *t k j* . The synaptic coupling strength between two neurons *i*, *j* takes the value *Jij* = *J*αβ, where α = {*E*, *I*} is a label denoting the population to which the postsynaptic neuron belongs, and β = {*E*, *I*} denotes the population to which the presynaptic neuron belongs. We define the external input to the network as μ*i* = μ (arriving at all excitatory neurons), and we also define a constant bias μ*i* = μ0 for all inhibitory neurons.

In this framework, we assume that each neuron of the network is characterized by a different distance-to-threshold value, which may be related with several biophysical properties (such as the membrane resistance, the firing threshold, or extra non-linear considerations). We assume here that such heterogeneity in the distance-to-threshold corresponds to heterogeneity in firing threshold values, although heterogeneities in distance-tothreshold can be translated into heterogeneities in other kind of parameters in more sophisticated neuron models. In particular, each excitatory neuron has a firing threshold θ*E*,*i* which is randomly distributed following a Gaussian profile *PE*(θ*E*) with mean θ and standard deviation *wE*. Equivalently, each inhibitory neuron has a firing threshold θ*I*,*i* randomly chosen from a Gaussian distribution *PI*(θ*I*) with mean θ and standard deviation *wI* (see **Figure 1**). Such heterogeneity serves to reflect some of the variability in the individual excitability properties of neurons found in actual neural systems, while treating separately the heterogeneity of excitatory and inhibitory cells will allow us to discern the effects caused by each population. For convenience, we will define a *low-threshold neuron* as a neuron whose firing threshold value (or, more precisely, its distance between the threshold for spiking and its resting state in absence of input) is below the average for its population (i.e., excitatory or inhibitory). A *high-threshold neuron* will therefore have a firing threshold value which is higher than the average for its population.

In the following, and unless specified otherwise, we choose *K* = 200 connections (in simulations, - = 0.2 for *N* = 1000 neurons), γ = 0.8, τ*m* = 20 ms, *Vr* = 10 mV, σ = 3 mV, θ = 20 mV, τ*ref* = 5 ms, *JEE* = *JIE* = 0.05 mV, and *JEI* = *JII* = −0.08 mV. These parameter values are within the physiological range for cortical neurons, and similar values have been used in previous modeling studies (Brunel and Hakim, 1999; Brunel, 2000). When computing the response of the system (for instance, the network mean firing rate for a given heterogeneity value), we average over 10 trials (or simulation runs on a random realization of the connectivity matrix) of 10 s each. The results presented in this work hold for these other parameter choices as well.

Together with the numerical simulations of the neural network described above, we have obtained an analytical meanfield solution of the model, which is described in detail in the Supplementary Material (Section Mean-field approach). In short, we have employed the diffusion approximation in the input to a single IF neuron to compute its mean firing rate in steady state conditions (Tuckwell, 1989; Brunel, 2000). Since the input to any given neuron will depend on the activity of the whole network (due to the recurrent nature of the system), we can average over the heterogeneity and obtain a mean-field description of the excitatory and inhibitory network mean firing rate, which will be given, respectively, by ν*E* = *E*(ν*E*, ν*I*, *wE*, *wI*) and ν*I* = *I*(ν*E*, ν*I*, *wE*, *wI*) (see Section Mean-field approach in Supplementary Material for an explicit form of these functions). An analytical estimation of the stability of this solution has been obtained as well (see Supplementary Material). The heterogeneity parameters *wE* and *wI* have an important effect on the mean firing rates, and allow us to use this mean-field solution, together with numerical simulations of the network, to explore the properties of the system.

## **3. RESULTS**

## **3.1. EFFECT OF HETEROGENEITY ON MEAN FIRING RATE**

Our first step is to understand the effect of an increase of the level of cell-to-cell heterogeneity on the stationary firing rate of the neurons in the network. Due to input noise and the sparseness of the network (which leads to a different number of incoming connections for each neuron), the neurons in the network are not characterized by a common unique mean firing rate (even in the absence of threshold heterogeneity), but rather each neuron has an individual mean firing rate, distributed around an average value following a Gaussian-like profile. This can be seen in **Figure 2A** (in light red, for excitatory neurons) and **Figure 2B** (in light blue, for inhibitory neurons). When we consider some degree of heterogeneity in the neuron firing thresholds, this original distribution of firing rates becomes wider and spans over a large range of firing rate values. For instance, a heterogeneity level of only *wE* = 2 mV for excitatory neurons lead the excitatory firing rate distribution from the previous narrow, Gaussian-like profile to a broad, long-tailed distribution which contains firing rates from zero to even tens of Hertz (**Figure 2A**, dark red). The same effect is observed for the inhibitory population: an increase in *wI* from 0 to 2 mV leads from a narrow, peaked distribution (**Figure 2B**, light blue) to a long-tailed one (dark blue). Excitatory (inhibitory) heterogeneity has also an effect on the shape of the inhibitory (excitatory) firing rate distribution, although it is not as strong as the effect of heterogeneity in a given population on the firing rate distribution of that same population (not shown).

Heterogeneity has a significant effect not only on the shape of the distributions, but on the mean firing rate of the populations in the network as well. The case of a single excitatory population was already considered in a previous work (Mejias and Longtin, 2012), where it was shown that an increase of neural heterogeneity triggered the appearance of a group of low-threshold neurons with higher firing rates (similar to the long tail in the dark-colored distributions in **Figures 2A,B**), and this group produced an extra recurrent input on the high-threshold neurons forcing them to increase their firing rate. The effect had a strong collective component, since a simple firing rate increase in the low-threshold neurons would have been at least partially compensated by a decrease in the high-threshold neurons, were they not connected to each other. The overall recurrent activity generated by the low-threshold neurons contributed to avoid a sudden drop in the firing rate of high-threshold neurons, yielding an overall quadratic-like increment in the network mean firing rate as the heterogeneity level increased.

The situation is more complex in the present case, where we have two different and interconnected populations of neurons (the excitatory and the inhibitory population), and also one heterogeneity parameter for each population. The results from the mean-field approach as well as from the numerical simulations can be seen in **Figures 2C,D**. **Figure 2C** shows the effect of increasing the excitatory heterogeneity in the activity level of the system. Our mean-field prediction, which agrees very well with numerical simulations, shows that increasing the heterogeneity level of the excitatory population leads to a rise in both excitatory and inhibitory activity. This can be easily understood by considering that the increase of *wE* produces the effect in the excitatory

firing rate described above, and this in turn increases the input entering from the excitatory to the inhibitory neurons, rising the inhibitory firing rate as well. The increment in the inhibitory rate also modulates back the excitatory population, which implies that the effect of *wE* on the excitatory population is not as pronounced as for the case of an isolated excitatory population.

The effects of increasing the inhibitory heterogeneity are, however, qualitatively different from those produced by its excitatory counterpart. As **Figure 2D** shows, increasing the inhibitory heterogeneity produces a rise in the inhibitory activity but decreases the excitatory activity. The origin of this effect is that increasing *wI* leads to the appearance of low-threshold inhibitory neurons with high firing rates, which increases the firing rate of the inhibitory network. This in turn induces more inhibition in the excitatory population, which lowers its level of activity as a consequence. It is interesting to note that, due to the negative character of the feedback within the inhibitory population, the increment in the inhibitory firing rate with the inhibitory heterogeneity is only due to the appearance of low-threshold inhibitory neurons, which pull the average firing rate up. On the other hand, we have three different factors that pull this average down: (i) the appearance of high-threshold inhibitory neurons, (ii) the decay in the positive contribution of the excitatory firing rate, and (iii) the presence of negative feedback within the inhibitory population. Because of this, the increase in the inhibitory firing rate with *wI* is not as strong as the increase of the excitatory firing rate with *wE*, where the feedback is positive. In particular, the effect of *wE* on the excitatory firing rate is about twice as strong as that of *wI* in the inhibitory firing rate, as one can see from the differences in the scale of the vertical axis in **Figures 2C,D**.

## **3.2. HETEROGENEITY AS A GAIN CONTROL MECHANISM**

After observing the strong effect that heterogeneity has on the mean firing rate of a cortical network for a given external input, an immediate question follows: how does neural heterogeneity influence the general input–output properties of cortical circuits? A first approach to answering this question is to analyze the effect of heterogeneity on the input–output dependence, or f-I curve, of the neural network. The f-I curve of a given neural system gives the relationship between a slow (usually considered constant) input to the circuit and the *readout* or mean firing rate of that circuit. There are a number of biophysical mechanisms which are able to modify or control the shape of this curve (a strategy commonly referred to as gain control). One can typically distinguish between several types of gain control, the most common ones being subtractive, divisive, or non-monotonic gain control (see Mejias et al., 2014 for an example of a system able to display these three regimes).

In this section, we will use our mean-field approach, together with numerical simulations, to address the role of neural heterogeneity as a relevant factor for gain control. Our results are shown in **Figure 3**, where different f-I curves are plotted, for both the excitatory and the inhibitory populations and varying either *wE* or *wI*. The effect of the excitatory heterogeneity on the f-I curves of the system is shown in the top panels of the figure, where we can see that increasing *wE* linearizes the f-I curve for both the excitatory (**Figure 3A**) and the inhibitory (**Figure 3B**) populations. The effect of the excitatory heterogeneity is similar to the noiseinduced linearization effect in a recurrent spiking neural network (Sutherland et al., 2009). As in the case of temporal fluctuations in noisy input currents, the effect of cell-to-cell heterogeneity is particularly important around the onset of the f-I curve, when most of the neurons lie between the fluctuation-driven and the mean-driven regime.

The effect of increasing the heterogeneity level in the inhibitory population is, however, more complex. As we can see in **Figure 3C**, the increase of *wI* leads to a decrease in the excitatory firing rate for any given input value. Interestingly, the effect of *wI* on the excitatory f-I curve is of a divisive nature, meaning that the inhibitory heterogeneity can be used as a divisive gain control parameter to perform multiplicative and divisive operations in cortical computations. Such a divisive gain control effect holds

**the f-I curves of the system.** The input in the horizontal axis corresponds to the external excitatory input μ. The inhibitory bias remains at μ0 = 17 mV at all times. Left panels show the effect on the excitatory population, while right panels do the same for the inhibitory population. In **(C)**, one can observe a clear divisive gain control of the excitatory f-I curve when *wI* varies. This divisive effect is not present in correspond to the *wI* = 0 case (red line) but rescaled by a constant factor to fit the other cases, indicating that the effect of *wI* can be described as divisive gain control. The constant factor ζ used to rescale the *wI* = 0 f-I curve was obtained for each case by minimization of the squared distance (see main text) and the resulting values are, for increasing *wI* , the following: 0.772, 0.531, 0.355, and 0.223.

for small and moderate values of the f-I curve, although as we move to strongly mean-driven conditions (i.e., μ ≥ 20 mV) the divisive control is lost (not shown). The main effect of increasing *wI* on the inhibitory population is an increase in the mean firing rate for a wide range of input values, as **Figure 3D** shows; no clear subtractive, divisive or linearization effect is apparent for the inhibitory cells as their heterogeneity is increased.

Due to the relative complexity of the model (and, in particular, its recurrent connectivity nature), it is not easy to obtain, even approximately, a theoretical proof of the divisive nature of the gain control effect of *wI*. Indirect measurements, however, can be used to test this hypothesis. For instance, in **Figure 3E** we have computed, using our mean-field approach, the f-I curve for the excitatory population and different values of *wI*. We have set μ0 = 12 mV and *JEI* = −0.4 mV so that the divisive effect is stronger and more easily identifiable than in **Figure 3C** for the range of biases shown. By taking the f-I curve for *wI* = 0 (in red) and multiplying it by a given constant factor ζ , one can obtain an f-I curve for *wI* > 0. To do this systematically, we have defined the squared distance between the rescaled *wI* = 0 f-I curve [namely ζ *r*0(*I*)] and a given *wI* > 0 f-I curve [namely *rw*(*I*)] as

$$\Delta = \frac{1}{n} \sum\_{i=1}^{n} \left[ \zeta \left[ r\_0(i) - r\_{\text{w}}(i) \right] \right]^2,\tag{4}$$

where the subindex *i* runs over all the input values considered in the numerical evaluation of the curve, and *n* = 25 is the total number of these values. By systematically varying the rescaling constant factor ζ , we find the value of this factor that minimizes the squared distance between both curves. This fitting is possible for all values of *wI* considered, and the squared distance at the optimal rescaling factor is always small (<0.003). For instance, by multiplying the original (*wI* = 0) f-I curve by a factor of ζ = 0.772, we obtain an f-I curve that fits very well ( < 0.0022) the f-I curve for *wI* = 2 mV. The good overlap between the rescaled *wI* = 0 curve (circles in **Figure 3E**) using different multiplicative constants and the original mean-field solutions (in colored lines), as demonstrated by the small values of obtained and graphically displayed in **Figure 3E**, indicates that the observed effect is indeed divisive. Simulation results are not displayed for **Figure 3E** for an easier visual comparison with the rescaled curves, although simulations agree very well with the mean-field predictions as in the previous set of parameters (see **Figure 3C** for a reference). We have further assessed the goodness of fit by checking that the residuals for each fit are distributed around zero, with approximately two thirds of the data points falling within one standard deviation of the data distribution, as expected for zero-mean Gaussian statistics. Other quantities for measuring the goodness of fit, like a normalized version of the quantity used here (which prevents our fit to depend on the average firing rate), also give the same results.

#### **3.3. STABILITY AND PHASE DIAGRAM**

So far, we have described the behavior of our cortical network model by assuming fixed point conditions, which led us to asynchronous steady-state solutions of the dynamics. Spiking neural networks are known to display other non-linear dynamics for certain conditions, such as multistability (Compte et al., 2000; Wang, 2001), fast global oscillations (Brunel and Hakim, 1999; Brunel, 2000), or winner-take-all dynamics (Wang, 2002; Wong and Wang, 2006). Although our aim in this study is not to fully characterize this kind of behavior in heterogeneous cortical networks– which would require more advanced calculations– we can study the local stability of the dynamics of the network to map the regions in the parameter space where our conclusions hold. To accomplish this, we compute the Lyapunov exponents of our system (see Equation 15 in Supplementary Material) and estimate, for a given value of the heterogeneity parameters, the maximum external input μ for which the asynchronous steady-state is stable (i.e., all eigenvalues have a negative real part). This limit would give us a clear frontier between asynchronous (top panel in **Figure 4A**) and synchronous (bottom panel in **Figure 4A**) network mean firing rate.

**Figure 4B** shows this maximum external input as a function of the excitatory heterogeneity, for different values of inhibitory heterogeneity. We can observe that, for each *wI* value, the maximum external input decreases linearly with the level of excitatory heterogeneity, as in the previously studied case of a purely excitatory heterogeneous network (Mejias and Longtin, 2012). This indicates that networks with highly heterogeneous excitatory neurons are able to enter the synchronous regime with less external stimuli than for the homogeneous case. The observed early synchronization in heterogeneous networks arises due to the presence of low-threshold excitatory neurons. This subset of neurons has a higher firing rate, and therefore they generate a stronger recurrent input that makes them closer to the bifurcation point from asynchrony to synchrony. As a consequence, low-threshold neurons become synchronous with less external input and they in turn contribute to the early synchronization of the rest of the neurons in the network, as we observe in **Figure 4B**.

The effect of the inhibitory heterogeneity is much less significant (see **Figure 4C**), although one can distinguish a small increase in the maximum external input for large enough values of *wI*. This is to be expected, since a large *wI* value would increase the inhibitory firing rate, inducing a decrease in the excitatory firing rate that must be compensated with a higher external input. Therefore, the synchronization onset will be located at a higher external input value. For both panels (**Figures 4B,C**), numerical simulations (points) agree very well with our mean-field predictions (lines). For large heterogeneity values (*wE*, *wI* ≥ 2.5 mV), the quenched disorder together with the stochasticity of the system make it difficult to accurately detect the synchronization onset. To avoid this problem, we have restricted our analysis to situations in which both the excitatory and inhibitory heterogeneity levels were small (*wE*, *wI* < 2.5 mV).

We have also used our mean-field approximation to compute a *wE* − *wI* phase diagram of the behavior of the system, which is shown in **Figure 4D**. For both *wE* and *wI* small, the system is in the asynchronous regime. The asynchronous state continues being stable for increasing *wI*, since the subsequent increment in the inhibitory firing rate contributes to stabilize the network dynamics as explained above (see **Figure 4C**). Only when *wI* takes moderate values and *wE* is significantly increased, synchronous behavior appears in the network dynamics. As the inhibitory

**FIGURE 4 | Stability of networks of heterogeneous neurons. (A)** Examples of the excitatory network mean firing rate in the asynchronous (top) and synchronous (bottom) regimes. **(B)** Critical external input as a function of the excitatory heterogeneity *wE* , for μ0 = 12 mV and different values of *wI* . An external input higher than the critical one will induce spontaneous synchronization in the network. Mean-field predictions (lines) agree with numerical simulations (symbols). **(C)** Critical external input as a function of *wI* , for μ0 = 12 mV and different values of *wE* . **(D)** Phase diagram in the *wE* − *wI* space, obtained using the mean-field approach from **(B,C)**, for an external input of 17 mV and different levels of the inhibitory population bias μ0. For all panels, σ = 1 mV.

bias μ0 takes larger values, the stabilizing effect of increasing *wI* reduces the area of the regime of synchronous dynamics. This is due to the fact that *wI* allows for a stronger modulation of the inhibitory firing rate when μ0 is larger, as the activity of the inhibitory low-threshold neurons will be higher in this case.

It is interesting to highlight that the results presented here (together with other recent studies such as Mejias and Longtin, 2012) provide counter-intuitive situations where heterogeneity promotes synchronization rather than impede it (Borgers and Kopell, 2003; Denker et al., 2004). A comprehensive study of the contrast between our results and the dynamical mechanisms previously reported is, however, beyond the scope of this study and will be addressed in future work.

## **3.4. SIGNAL DETECTION**

Since both excitatory and inhibitory heterogeneity have a significant impact on the input–output characteristics for constant input, it is convenient to extend our analysis to consider the effect of heterogeneity in the transmission of more realistic, time-varying signals. In particular, previous work showed that the presence of a certain level of heterogeneity can optimize the transmission of slow signals under rate coding in excitatory populations (Mejias and Longtin, 2012). The phenomenon was also present in cortical-like networks with sparseness and inhibition, although in this more realistic case, no theoretical approximations were provided to support these claims.

The mean-field approximation presented in the Supplementary Material (see Section Mean-field approach) constitutes a useful tool to investigate these heterogeneityinduced resonances in cortical-like network models, and to evaluate the likelihood of this phenomenon to occur in real cortical circuits. We consider an external input constituted by

$$
\hat{\mu}(t) = \mu + \text{S}\_0 \sin \left(2\pi f\_\text{\textdd} \right),
\tag{5}
$$

with the first part of the r.h.s. being a constant input and the second part being a slow and weak modulation. Such a weak input modulation is able to drive the excitatory mean firing rate of the network under certain circumstances, a situation which is shown in the inset of **Figure 5B**. A convenient measurement to quantify this behavior is the zero-lag input–output covariance function, which is given by

$$C \equiv \langle \hat{\mu}(t) \nu\_E(t) \rangle - \langle \hat{\mu}(t) \rangle \langle \nu\_E(t) \rangle. \tag{6}$$

**Figure 5A** shows this input–output covariance as a function of the excitatory heterogeneity. For networks of homogeneous neurons, the modulation part of the signal is typically too weak to be noticed. For higher values of the excitatory heterogeneity, however, the sensitivity of the network to small inputs increases, due to the existence of a larger number of low-threshold excitatory neurons. As a consequence, the small input modulation is now able to strongly drive the output (i.e., the excitatory network mean firing rate). For even higher values of *wE* the overall activity of the network increases drastically and the variations due to input modulations become slightly weaker. This makes the signaldriven firing rate modulations too small compared to the baseline firing rate, and as a consequence the quality of signal transmission to a given linear readout system decreases. The overall effect is a bell-shape dependence of the input–output covariance with *wE* shown in **Figure 5A**, indicating that a certain non-zero level of

by *S*0 = 0.5 mV and *fs* = 2 Hz. **(B)** Same as in **(A)**, but as a function of excitatory heterogeneity optimizes signal transmission, as in the

findings (symbols). Other parameters are *JEE* = *JIE* = 0.043 mV, *JEI* = *JII* = −0.06 mV, and μ = μ0 = 15 mV, and the signal is characterized

simpler case of a purely excitatory population studied in Mejias and Longtin (2012). Furthermore, our novel mean-field approach for excitatory and inhibitory sparse populations closely follows the numerical results.

The effect of inhibitory heterogeneity on signal transmission is notably different from the situation explained above. As we can see from **Figure 5B**, the input–output covariance tends to be weaker in networks whose inhibitory neurons are highly heterogeneous. This is due to the fact that inhibitory heterogeneity causes an increase in the inhibitory firing rate, which reduces the sensitivity of the excitatory population to weak stimuli and therefore hinders its capacity for signal detection. This is consistent with recent experimental and theoretical findings which show that correlations between two neurons decrease as their firing rate decrease (de la Rocha et al., 2007), and it suggests that heterogeneity in inhibitory neurons may have an important role in decorrelation between input signals and neural activity.

We can also see that the detection of the signal is frequencydependent. In **Figure 5C** the input–output covariance as a function of *wE* is computed for different values of the signal frequency *fs*, with a small shift of the peak toward more heterogeneous networks as the input frequency is increased. This behavior was also observed for the case of one isolated excitatory population (Mejias and Longtin, 2012), and suggests that the ability of neural networks to efficiently detect and transmit signals of a given frequency range depends on the heterogeneity level of the network. The decrease of the signal detection as a function of the inhibitory heterogeneity also depends on the input frequency considered (**Figure 5D**), although for large enough *wI* the signal

since the mean-field predictions does not hold for high input frequencies. **(D)** Same as **(B)**, but for different input frequencies and only considering

the numerical results. **(C,D)** share the same color code.

## **3.5. MULTIPLE HETEROGENEITY SOURCES**

is not detected regardless of the frequency.

As a final remark, it should be noted that, in all of the simulations and analyses presented so far, either the excitatory or the inhibitory heterogeneity was varied, while the heterogeneity of the other population was kept fixed at a very low level (0.1 mV, which would correspond to a almost homogeneous population). In the situations in which the system under study behaves in a linear fashion and the effects caused by parameter variations are independent, this approach is convenient to systematically characterize the behavior of the system. The neural network under study, however, is known to display multiple kinds of non-linear behavior (such as, for instance, the non-linear dependence of the mean firing rate with the heterogeneity level, as shown in **Figures 2C,D**). It is, therefore, unclear whether one can infer the response of the system for arbitrary combinations of heterogeneity parameters from the curves and results presented in previous sections, in which mainly only one type of heterogeneity was analyzed at a time.

In order to test the validity of our results in more complex scenarios, we have jointly increased both heterogeneity levels (*wE* and *wI*) at the same time, and we numerically computed the excitatory mean firing rate of the network as a function of this *combined* heterogeneity level *W* (with *W* = *wE* = *wI*). To test the linearity of the system to the presence of multiple sources of heterogeneity, we also compute, using the mean-field solution, the changes in the excitatory firing rate due to only *wE* or only *wI* (while keeping the other heterogeneity level at zero value), and we add these two contributions together. The comparison is shown in **Figure 6**, where we can see that the simulation results closely follows the linear prediction up to values of the combined heterogeneity of ∼ 4 mV. This finding implies that the effects of multiple types of heterogeneity can add up linearly in some parameter regimes, such as the one in this preliminary investigation. Therefore, the results of the present work are also valid for more realistic situations in which different types of heterogeneity (i.e., excitatory and inhibitory) are simultaneously present in the system.

## **4. DISCUSSION**

The importance and roles of intrinsic neuronal heterogeneities on the dynamics of neural networks is starting to being uncovered in recent years (Marsat and Maler, 2010; Padmanabhan and Urban, 2010; Savard et al., 2011; Angelo et al., 2012; Tripathy et al., 2013). Although commonly disregarded in most modeling studies, an increasing level of attention has been drawn to the subject by theoretical and computational models as well (Golomb et al., 2001; Denker et al., 2004; Luccioli and Politi, 2010; Mejias and Longtin, 2012; Nicola and Campbell, 2013). In particular, novel theoretical frameworks for addressing the heterogeneity of neural systems

have been proposed recently. Yim and colleagues, for example, propose a theoretical approach especially useful for addressing the relationship between neural heterogeneity and neural correlations (Yim et al., 2013), and they sketch a possible explanation for recent evidence of a positive role of heterogeneity on population coding (Padmanabhan and Urban, 2010). In another recent work, Nicola and Campbell provide a set of mean-field approaches used to shed light onto a heterogeneity-induced change on the nature of the Hopf bifurcation responsible for burst generation (Nicola and Campbell, 2013).

The theoretical understanding of the effects of heterogeneity on neural systems is still a young problem, though, and only simple situations have been considered up to now. In this work, we have analytically and computationally studied the interplay between population-specific levels of cellular heterogeneity, an important problem that has not been properly addressed to date. Interestingly, the effects that excitatory heterogeneity produces on neural networks are quite different from the ones produced by inhibitory heterogeneity. Excitatory heterogeneity, as we have shown, non-linearly increases the network mean firing rate with respect to that of a homogeneous network, and the f-I curves of the system are linearized as well. In this sense, excitatory heterogeneity may be viewed as a classical *quenched disorder* in excitable systems, with similar effects on the f-I curve than that of pure noise (Doiron et al., 2001). On the other hand, the introduction of inhibitory heterogeneity induces an increase (with respect to homogeneous networks) in the inhibitory firing rate and a decrease in the excitatory firing rate, and a divisive modulation of the f-I curve as a result. Divisive gain control mechanisms is often assumed as a key operation for neural computations (Carandini and Heeger, 1994; Chance and Abbott, 2000), but biophysical mechanisms for such a modulation have been hard to identify, regardless of being network-based mechanisms (Holt and Koch, 1997; Doiron et al., 2001; Chance et al., 2002; Mejias et al., 2014) or cell-based mechanisms (Prescott and De Koninck, 2003; Mehaffey et al., 2005). The identification of neuronal heterogeneity in inhibitory populations as a biophysically realistic mechanism for multiplicative and divisive gain control constitutes one of the key achievements of the present study.

The analysis of the stability of the fixed point solutions of heterogeneous networks also provides useful information about the effects of heterogeneity on neural networks. Again, the effects of neural heterogeneity heavily depend on the population in which it is found. Excitatory heterogeneity leads to an easier spontaneous synchronization of the neural network, while inhibitory heterogeneity has a weaker effect and tends to slightly increase the robustness of the asynchronous state. This produces a rich repertoire of stability behaviors in neural networks, with the stability conditions of a given particular network depending on its balance between excitatory and inhibitory heterogeneity.

In recent works, an optimal information transmission has been shown to occur for heterogeneous populations of neurons (Marsat and Maler, 2010; Padmanabhan and Urban, 2010), and the presence of short-term synaptic plasticity has been suggested to increase the efficiency of coincidence detection in the presence of heterogeneity via the appearance of *optimal frequencies* (Mejias and Torres, 2008). These findings indicate that heterogeneity may have an important role in information transmission in neural systems. In the present study, we have demonstrated here that the optimization of signal detection by networks of heterogeneous neurons under rate coding, first described in Mejias and Longtin (2012), holds for the more realistic cortical-like network used here, by means of both numerical simulations and mean-field approaches. The improvement of signal detection in heterogeneous neural and excitable systems has been a recent focus of interest. For instance, Tessone et al. (2006) found that global synchronized events in response to weak, slowly modulated external signals can be optimized in heterogeneous networks, a result that has been also obtained in neural networks with electrical and chemical synapses (Perez et al., 2010). Global synchronized events in heterogeneous networks can also work at very short time scales, being triggered by fast input and allowing for an efficient temporal coding (Mejias and Longtin, 2012). Recent experimental work has also highlighted an optimization of population coding in networks of heterogeneous neurons (Marsat and Maler, 2010; Padmanabhan and Urban, 2010; Savard et al., 2011; Angelo et al., 2012), establishing a solid ground for neural heterogeneity as a key ingredient of neural coding.

Finally, it is worth noting that, although we have studied exclusively the case of heterogeneity in the distance-to-threshold of LIF neurons, our mean-field approach can be used to study heterogeneity in other parameters as well. Indeed, there are many potential biophysical sources of heterogeneity in neural systems, both at the network level (i.e., heterogeneity in the network connectivity, as in Olmi et al., 2010) and at the neuron level. In this second group, possible heterogeneity sources can be defined in terms of anatomical and morphological properties, or also at a functional level, including neuronal excitability (Tessone et al., 2006; Perez et al., 2010), different degrees of spike frequency adaptation (Hemond et al., 2008; Nicola and Campbell, 2013), or other biophysical properties (Padmanabhan and Urban, 2010; Tripathy et al., 2013), to name a few. Understanding their individual or joint role in neural dynamics will require future modeling work at different scales and levels of detail, for which mean-field approaches could be of great help.

#### **ACKNOWLEDGMENTS**

This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC) and in particular its Discovery Accelerator program, and the Canadian Institutes for Health Research (CIHR). Jorge F. Mejias would like to thank Lidia Roman-Gonzalez for help with some figures.

#### **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fncom.2014. 00107/abstract

## **REFERENCES**

**Conflict of Interest Statement:** 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.

*Received: 23 April 2014; accepted: 21 August 2014; published online: 12 September 2014.*

*Citation: Mejias JF and Longtin A (2014) Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks. Front. Comput. Neurosci. 8:107. doi: 10.3389/fncom.2014.00107*

*This article was submitted to the journal Frontiers in Computational Neuroscience. Copyright © 2014 Mejias and Longtin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Heterogeneity of heterogeneities in neuronal networks

## *Fabiano Baroni 1,2,3\* and Alberto Mazzoni 4*

*1 School of Psychological Sciences, Monash University, Clayton, VIC, Australia*

*2 NeuroEngineering Laboratory, Department of Electrical and Electronic Engineering, University of Melbourne, Parkville, VIC, Australia*

*3 Centre for Neural Engineering, University of Melbourne, Carlton, VIC, Australia*

*4 The BioRobotics Institute, Scuola Superiore Sant'Anna, Pontedera, Italy*

*\*Correspondence<EMAIL_ADDRESS>
#### *Edited by:*

*Mark D. McDonnell, University of South Australia, Australia*

#### *Reviewed by:*

*Roberto Fernández Galán, Case Western Reserve University, USA Benjamin Lindner, Bernstein Center for Computational Neuroscience, Germany*

**Keywords: integrate-and-fire neurons, excitation/inhibition balance, neuronal network dynamics, biophysical heterogeneity, dynamical heterogeneity, signal modulation, high-frequency oscillations**

**A commentary on**

**Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks**

*by Jorge F. Mejias, André Longtin (2014). Front. Comput. Neurosci. 8:107. doi: 10.3389/fncom.2014.00107*

Neurons in the brain exhibit a broad spectrum of heterogeneities even within a given morphological or physiological class. In a recent modeling study, Mejias and Longtin investigated the effects of heterogeneity in the voltage threshold for spike generation on the dynamics of random networks of excitatory and inhibitory neurons (Mejias and Longtin, 2014), hence extending their previous results on purely excitatory networks (Mejias and Longtin, 2012).

The authors focused on the different effects of heterogeneity when incorporated in either the excitatory or the inhibitory population. A greater heterogeneity in the excitatory population increases the average firing rate of both neuron types, because the subset of most excitable excitatory neurons provides a positive feedback to the whole network. Conversely, when heterogeneity is included in the inhibitory population, only the average firing rate of inhibitory neurons increases, while the average firing rate of excitatory neurons decreases. This result can be explained by the presence of lowthreshold, highly excitable inhibitory neurons which tend to shift the average firing rate of the inhibitory population to higher levels, while increasing inhibitory currents in the network. While the silencing effect of an increase in inhibitory currents in the heterogeneous inhibitory population is overcome, on average, by the presence of highly excitable neurons, it dominates the dynamics of the homogeneous excitatory population, hence reducing the average firing rate of excitatory neurons.

While firing rates increase with heterogeneity of excitatory cells, potentially leading to run-away excitation in the absence of saturation or adaptation mechanisms, other features of network dynamics exhibit a non-monotonic dependence on heterogeneity. For example, the encoding of an oscillatory input signal is optimal when the network exhibits an intermediate level of heterogeneity, consistently with a recent experimental study (Tripathy et al., 2013). This behavior is reminiscent of stochastic resonance, a general phenomenon observed in excitable systems, whereby intermediate levels of noise enable optimal information encoding. In fact, heterogeneity can be considered as a form of spatial noise.

When discussing about heterogeneity in the nervous system, it is important to distinguish between biophysical heterogeneity, which relates to neuronal parameters (in simulations) or quantities that are static in the time scales of interest (in experiments), and dynamical heterogeneity, which refers to measures of ongoing neuronal activity such as firing rates and correlations. The relationships between the two can be usefully explored in both directions: while Mejias and Longtin explored the dynamical consequences of different levels of biophysical heterogeneity (bottom-up), others started from experimental observations of dynamical heterogeneity, and investigated neuronal models that are consistent with the observed dynamics (topdown, Koulakov et al., 2009; Roxin et al., 2011).

Crucially, different biophysical sources of heterogeneity can yield similar effects at the level of network dynamics. For example, the strongly skewed, lognormallike distribution of firing rates typically observed in large-scale neuronal recordings (recently reviewed in Buzsáki and Mizuseki, 2014) can be explained by models that include nonrandom connectivity among linear neurons (Koulakov et al., 2009), as well as by homogeneous networks with random connectivity among more realistic nonlinear neurons, due to the expansive nonlinearity of the *f-I* curve (that is, the superlinear increase in output firing rate *f* with increasing input current *I*) in the presence of noise (Roxin et al., 2011). Similarly, delay and synaptic time scale diversity yield equivalent effects (Biggio et al., 2013).

Further, neuronal heterogeneity can arise from different biophysical substrates, and how the effects of different sources of heterogeneity interact is unclear. In fact, heterogeneity has been observed in virtually every aspect of neuronal physiology where it has been investigated. These include intrinsic neuronal properties, which are mostly determined by ion channels' dynamics, such as neuronal excitability, propensity to bursting, amplitude and time course of spikefrequency adaptation, post-inhibitory rebound, and many more (e.g., Marder, 2011; Angelo et al., 2012); as well as properties related to the connectivity among neurons, such as dynamics and efficacy of synaptic transmission, density and size of dendritic spines, thickness and myelination of axons (e.g., Dobrunz and Stevens, 1997; Song et al., 2005; Wang et al., 2008).

Mejias and Longtin reported that, in their simulations, the effects of heterogeneity in the distance-to-threshold of excitatory and inhibitory neurons summed linearly when combined. However, we should expect that, in general, the effects of heterogeneity in different biophysical parameter might interact in a complex manner. For example, Roxin et al. reported that heterogeneity in synaptic weights can decrease the variability in firing rates caused by the expansive nonlinearity of the *f-I* curve, contrary to expectations of a linear interaction between these two sources of variability (Roxin et al., 2011). The effects of neuronal heterogeneity are expected to be more complex in more physiologically relevant settings, where different sources of heterogeneity coexist and are distributed in a highly non-random fashion (Yassin et al., 2010).

While biophysical heterogeneity is widespread at all levels of description, we believe that important insights into the relationships between biophysical and dynamical heterogeneities can be obtained using reductionist approaches, where the degree of biophysical diversity can be described by few parameters. For example, neuronal network studies that investigated the role of connectivity heterogeneity yielded important insights by focusing on random and scale-free connectivity (described by a single parameter), both of which display important dynamical differences with respect to homogeneous all-to-all connectivity. We propose that a similar approach could be fruitfully applied to other forms of biophysical heterogeneity, and ultimately result in useful taxonomies of the different sources of biophysical heterogeneity, describing the dynamical heterogeneities they result in and the interactions between their effects. This level of understanding would facilitate the conceptual integration of different results and eventually lead to basic functional principles of neuronal processing beyond area- or species- specific details.

## **REFERENCES**

**Conflict of Interest Statement:** 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.

*Received: 16 October 2014; accepted: 21 November 2014; published online: 10 December 2014.*

*Citation: Baroni F and Mazzoni A (2014) Heterogeneity of heterogeneities in neuronal networks. Front. Comput. Neurosci. 8:161. doi: 10.3389/fncom.2014.00161*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2014 Baroni and Mazzoni. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## A Markov model for the temporal dynamics of balanced random networks of finite size

## *Fereshteh Lagzi\* and Stefan Rotter*

*Bernstein Center Freiburg and Faculty of Biology, University of Freiburg, Freiburg, Germany*

#### *Edited by:*

*Mark D. McDonnell, University of South Australia, Australia*

#### *Reviewed by:*

*Maurizio Mattia, Istituto Superiore di Sanità, Italy Jonathan D. Touboul, Collège de France, France Richard Hardstone, Vrije Universiteit Amsterdam, Netherlands*

#### *\*Correspondence:*

*Fereshteh Lagzi, Bernstein Center Freiburg and Faculty of Biology, University of Freiburg, Hansastraße 9a, 79104 Freiburg, Germany e-mail: fereshteh.lagzi@ bcf.uni-freiburg.de*

The balanced state of recurrent networks of excitatory and inhibitory spiking neurons is characterized by fluctuations of population activity about an attractive fixed point. Numerical simulations show that these dynamics are essentially nonlinear, and the intrinsic noise (self-generated fluctuations) in networks of finite size is state-dependent. Therefore, stochastic differential equations with additive noise of fixed amplitude cannot provide an adequate description of the stochastic dynamics. The noise model should, rather, result from a self-consistent description of the network dynamics. Here, we consider a two-state Markovian neuron model, where spikes correspond to transitions from the active state to the refractory state. Excitatory and inhibitory input to this neuron affects the transition rates between the two states. The corresponding nonlinear dependencies can be identified directly from numerical simulations of networks of leaky integrate-and-fire neurons, discretized at a time resolution in the sub-millisecond range. Deterministic mean-field equations, and a noise component that depends on the dynamic state of the network, are obtained from this model. The resulting stochastic model reflects the behavior observed in numerical simulations quite well, irrespective of the size of the network. In particular, a strong temporal correlation between the two populations, a hallmark of the balanced state in random recurrent networks, are well represented by our model. Numerical simulations of such networks show that a log-normal distribution of short-term spike counts is a property of balanced random networks with fixed in-degree that has not been considered before, and our model shares this statistical property. Furthermore, the reconstruction of the flow from simulated time series suggests that the mean-field dynamics of finite-size networks are essentially of Wilson-Cowan type. We expect that this novel nonlinear stochastic model of the interaction between neuronal populations also opens new doors to analyze the joint dynamics of multiple interacting networks.

**Keywords: Markov process, self-consistent noise model, balanced random network, nonlinear dynamics, networks of finite size, Wilson-Cowan type model**

## **1. INTRODUCTION**

Cortical neurons of behaving animals show highly irregular patterns of activity. One hypothesis for the source of such irregularity is the balance of excitation and inhibition in the steady state activity of the network (Softky and Koch, 1993; Bell et al., 1994; Shadlen and Newsome, 1994, 1998; Tsodyks and Sejnowski, 1995; van Vreeswijk and Sompolinsky, 1996). Experimental evidence in favor of this hypothesis suggest that excitation-inhibition balance is the principle of brain dynamics (Sanchez-Vives and McCormick, 2000; Shu et al., 2003; Haider et al., 2006; Okun and Lampl, 2008). The balanced state is an emergent self-consistent and stable solution of the temporal dynamics of the network (van Vreeswijk and Sompolinsky, 1996, 1998; Amit and Brunel, 1997b; Brunel, 2000). In other words, in a recurrent balanced network, both excitatory and inhibitory activity are shaped such that in cooperation with each other, they generate a stationary, selfconsistent input-output behavior on the level of the mean and the fluctuations. The collective activity of the involved neuronal populations include weakly correlated and irregular spike trains. Due to its stochastic appearance, this feature is referred to as "self-generated noise" or simply "noise" in this paper. In fact, these fluctuations are generated mainly by the complex recurrent interactions in the network, even in absence of any external source of noise (van Vreeswijk and Sompolinsky, 1996; Kriener et al., 2008). As in our model there is no external source of noise, the fluctuations are very likely due to deterministic chaos in a high-dimensional system (for details see Jahnke et al., 2009).

Temporal fluctuations of neuronal activity reflect brain processes. Transient activity of neuronal networks, for instance, correspond to different neural computations at different stages of a cognitive task (see for example Churchland et al., 2011) or the representation of information in the brain (Destexhe and Contreras, 2006). Fluctuations also influence sensory perception in the case of ambiguous input. This phenomenon has been modeled by a multi-stable noise-driven dynamical system in which activity fluctuations cause transitions between meta-stable fixed points (Moreno-Bote et al., 2007; Deco and Romo, 2008). It was shown that dynamic noise in a neuronal network also gives rise to different dynamical states of the network. Both theoretical and simulation studies have addressed the role of noise for dynamic stability, or for the emergence of oscillations in network dynamics (Brunel and Hansel, 2006; Ghosh et al., 2008; Touboul et al., 2011; Cai et al., 2012) which may have implications for brain functions. Therefore, to understand the functional properties of neuronal networks, it is essential to understand the dynamics of the time dependent variability in such systems. Theoretical studies of balanced random networks indicate that the fluctuations are essentially determined by two factors: neuronal correlations and the finite size of the network (Ginzburg and Sompolinsky, 1994; Brunel and Hakim, 1999; Brunel, 2000).

To address the fluctuations of the activity of a complex high-dimensional system such as a spiking neural network, a reduced low-dimensional description of network activity is needed. However, to compensate for the loss of degrees of freedom, an analytical treatment is needed such that the essential properties of the stochastic dynamics of the system fluctuations are preserved to an acceptable degree. Finding a suitable stochastic model to replace the spiking dynamics of a network is a challenge. One reason is that such networks are hybrid systems, as the membrane potential of each neuron is a continuous variable, and the spiking activities are discrete quantities. Secondly, an appropriate time scale has to be defined because the amplitude and the dynamics of the fluctuations depend on the temporal resolution. Thirdly, the stationary activity of the reduced model should maintain the statistical properties of the population activities, as their statistics are crucial for the switching dynamics of a network with two or more interacting populations (Bressloff and Newby, 2013). Finally, each neuron is a highly nonlinear element in a random network, due to its threshold and refractoriness. It has been hypothesized that in the thermodynamic limit, when the number of neurons becomes very large, the global dynamics is linearized due to the negative feedback from the inhibitory population (van Vreeswijk and Sompolinsky, 1996, 1998; Tetzlaff et al., 2012; Helias et al., 2013). However, in networks of finite size, as we show in this paper, the dynamics are not fully described by a linear framework.

An influential study on population dynamics of excitatory and inhibitory neurons has been performed more than 40 years ago (Wilson and Cowan, 1972). The authors considered an infinitely large number of neurons in each population such that a fraction of all neurons in the population are in the refractory state. They derived an *ad hoc* response function for the non-refractory neurons and using coarse-graining of activities in time, they derived a set of coupled ordinary differential equations. However, correlations and finite size fluctuations were not captured by this model. In essence, without fluctuations in the input, there are no output fluctuations in this model. In fact, most theoretical studies of population interactions using a mean-field approach (see e.g., Gerstner, 1995; Amit and Brunel, 1997a,b; Brunel, 2000; Aviel and Gerstner, 2006; Kriener et al., 2008; Toyoizumi et al., 2009; Cardanobile and Rotter, 2011; Ledoux and Brunel, 2011; Ostojic and Brunel, 2011) have not addressed the consequences of the finite size of the network and/or pairwise correlations for the temporal dynamics of the population activities in a self-consistent way.

There are a few explicit or implicit suggested approaches to study finite size population dynamics. The first type of studies are based on deterministic equations using mean-field approaches that are derived by an external source of noise (Kriener et al., 2008; Toyoizumi et al., 2009; Tetzlaff et al., 2012). Seminal studies on balanced random networks show that a stochastic input is not needed for network fluctuations and the noise in the system is self-generated as a result of recurrent activity and stochastic spiking of neurons (van Vreeswijk and Sompolinsky, 1996, 1998). Therefore, this approach cannot describe the temporal dynamics in a self-consistent way. Moreover, a deterministic set of equations with additive noise (Kriener et al., 2008) cannot reproduce the state dependent fluctuations of the network activity, as we will show in our study.

A second class of studies considered independent neurons with Poisson statistics (Brunel and Hakim, 1999; El Boustani and Destexhe, 2009). This approach could lead to unrealistic number of spike counts in a short time bin. Moreover, with the assumption of uncorrelated neurons, neural Poisson statistics result in network Poisson statistics. As we show in this article, the statistics of the population activities in balanced random networks are not Poissonian. Exploiting a general escape noise model in a point process framework, Spiridon and Gerstner (1999) derived an integral equation for the population activity of a fully coupled network. According to their approach, the finite size effect would show up as a multiplicative noise term in the original equation. As we show later, our analysis supports these results to some extent.

A third approach tries to describe either the dynamics of each neuron (Ohira and Cowan, 1993, 1995; Soula and Chow, 2007; El Boustani and Destexhe, 2009) or the network (Touboul and Ermentrout, 2011; Buice and Chow, 2013a) by a Markov process. The latter has the problem that possible jumps of the Markov process are limited to the immediate neighbors of each state; meaning that the number of active neurons in each population can either increase or decrease by one. The former approach seems to be able to better capture the dynamics and statistics of the network. Soula and Chow (2007) assumed a Markov model in discrete time with a time scale in the range of membrane time constant in the case of instantaneous synapses. The transition rate of a typical neuron in the network is calculated from the stationary firing rate and includes the net amount of excitation in the system. However, balanced random networks operate on a much faster time scale compared to that of a single neuron. Also, the interaction between excitation and inhibition, and the consequences of the negative feedback on temporal dynamics are not analyzed in this study. In another study, following a similar approach in continuous time, a master equation for the activity of a balanced network with current-based and conductance-based synapses was derived (El Boustani and Destexhe, 2009). This method needs the static transfer function of a single neuron that maps input rates to output rates.

Finally, population density approaches based on a conservation law imposed on the probability flux (the number of neurons is constant) (for details see Knight, 1972; Abbott and van Vreeswijk, 1993; Treves, 1993; Knight et al., 1996; Omurtag et al., 2000; Sirovich et al., 2000; Haskell et al., 2001; Nykamp and Tranchina, 2001; Mattia and Del Giudice, 2002) are yet another way of deriving stochastic dynamics from deterministic equations. In a study on the temporal dynamics of the interaction between excitation and inhibition in networks of finite size, an eigenfunction expansion of the essentially nonlinear Fokker-Planck equation resulted in a set of coupled ordinary differential equations (Mattia and Del Giudice, 2002, 2004). However, the noise term included in the model to account for the finite size of the network, was a white noise of fixed amplitude. Recently, Buice and Chow (2013a,b), using a population density method and moment hierarchies of the equation, derived path integrals of a moment generating functional. To get a time dependent correlation function of the system, by introducing a small perturbation, they applied a linear expansion of the equation of moments (Buice and Chow, 2013a,b).

In this paper, we aim at describing statistics and dynamics of finite-size fluctuations in a balanced random network of excitatory and inhibitory neurons self-consistently, such that the temporal dynamics of the network is driven by the interactions in the network. Moreover, a high correlation between the excitatory and inhibitory population, in their stochastic representation in the model, has to be preserved. To this end, a typical neuron in the network is modeled by a two-state Markov system that its transition probabilities needs to be derived dependent on the network activity states. Our model is based on the state space analysis of numerical simulations of interactions between the excitatory and the inhibitory population in a large balanced network, in the regime dominated by inhibition.

Systematic analysis of the two-dimensional population spike counts shows that no dynamic model with additive Gaussian noise can fully describe the temporal dynamics of the network activity. Specifically, the more excitation and the less inhibition is recruited in the network at any given point in time, the higher is the variance of the self-generated noise. It will be demonstrated in this paper that a stationary external input results in nonlinear interactions between excitation and inhibition. Moreover, we will show that the self-generated and state dependent noise emerges naturally as a result of the finite number of neurons in the network. Furthermore, the suggested two-state Markov model is capable of producing a heavy-tailed (positively skewed) distribution of excitatory and inhibitory spike counts, a property of balanced random networks that is considered here for the first time. We show that this heavy-tailed distribution can be well approximated by a log-normal distribution.

#### **2. MATERIALS AND METHODS**

In this section, we first describe our assumptions about the structure and the parameters of the network. Then, we show how to reconstruct the dynamic flow based on the results of spiking network simulations. In Section "Estimation of the dynamic flow underlying the mean-field dynamics" we first suggest a Markov model to represent the collective activity of neuronal populations, with transition probabilities inferred from numerically simulated networks of spiking neurons. The objective of this study is to identify a suitable bin size to represent the temporal dynamics of the network fluctuations. Then, in Section "Markov model for the mean-field and stochastic dynamics" we introduce a two-state Markov model for each neuron, termed the "Active-Refractory Markov" (ARM) model. A method to find the mean-field equations and a self-consistent noise model based on the Markovian single-neuron dynamics is introduced at the end of this section.

#### **2.1. NETWORK STRUCTURE AND PARAMETERS**

The network under study is composed of 10,000 excitatory and 2500 inhibitory neurons, similar to the network studied by Brunel (2000). Each neuron receives local inputs from randomly chosen fraction of the excitatory and the inhibitory population (10% each). An external Poisson process of rate 25 spikes/ms mimics the input from other brain areas. The neurons are modeled by a leaky integrate-and-fire (LIF) dynamics with pulse-like post-synaptic currents (PSC), and exponential post-synaptic potentials (PSP). Therefore, the dynamics of neuron *i* in the network, regardless of whether it is excitatory or inhibitory, satisfies

$$\tau \dot{\nu}\_i(t) = -\nu\_i(t) + \tau \sum\_{j=1}^{N} J\_{\vec{l}\vec{j}} \mathbf{S}\_{\vec{j}}(t - t\_d) \tag{1}$$

where *Sj*(*t*) = *k* δ(*t* − *t j k*) is the spike train of neuron *j*, which is seen by the postsynaptic neuron with a delay of *td*, and integrated by it with a membrane time constant τ = 20 ms. *Jij* is the amplitude of the post synaptic potential (PSP), *J* = 0.1 mV for the available excitatory synapses to each neuron and −*gJ* = −0.6 mV for inhibitory synapses. *N* is the total number of neurons in the network. The membrane potential of each neuron, once reached the threshold at θ = 20 mV, is reset to *v*reset = 10 mV and a spike will be generated. The membrane potential remains at *v*reset for a refractory period of 2 ms. These parameters are identical to those studied in Brunel (2000) and represent a simplified model for one cubic millimeter of neocortical tissue. The results of our study are valid for a wide range of biologically realistic parameters, as far as the network is in the inhibition dominated regime. Numerical simulations of the network were all conducted in NEST (Gewaltig and Diesmann, 2007) with a time resolution of 0.01 ms and a minimal synaptic delay *td* equal to the time resolution. The network simulation time was 100 s. As a result, a total number of 107 data points for each population in a histogram with a bin size of 0.01 ms was available for further data analysis.

## **2.2. ESTIMATION OF THE DYNAMIC FLOW UNDERLYING THE MEAN-FIELD DYNAMICS**

We reconstructed the flow corresponding to the mean-field dynamics from a simulated time series of the spiking network activities. To this end, for each possible combination (*i*, *j*) of excitatory spike count *i*, and inhibitory spike count *j*, the "state" of the system, we first collected the corresponding derivatives by computing the increments in successive time bins divided by the time bin, for both excitatory and the inhibitory spike counts. Taking the mean of the encountered derivatives for each state gave the average direction into which the system moved forward in time. This way, for each state visited in the simulation, a velocity vector is obtained, yielding a vector field that approximates the flow in a two-dimensional state space. Calculating the variance of the state dependent derivatives results in the state dependent variance of the self generated noise (fluctuations) in the system (see Results for more details).

## **2.3. MARKOV MODEL FOR THE MEAN FIELD AND STOCHASTIC DYNAMICS**

As mentioned in the beginning of this section, we considered two different Markov models for different purposes. First, for the selected time bin *dt*, chosen such that a Markov model is capable of reproducing the power spectrum of the network activity with the highest fidelity, a state transition matrix was inferred entirely based on the data obtained from a spiking network simulation. It is important to stress that the Markov model for the selected bin size *dt* provides a compressed description of the network dynamics, assuming that only the most recent bin, and no longer account of the history of the network activity, determines the dynamics of the network at any given point in time. The analysis described here explores the limits of such a description of the the large-scale dynamics, and in this paper we refer to it as "Markov chain analysis." However, deriving the transition matrix from a self-consistent analysis is difficult. Therefore, we looked for a more analytically tractable approach in a second step. This approach, which is the main focus of this study, is based on a two-state Markov description of individual neurons. As shown in this article, this model is able to describe both the dynamics and the statistics of the collective activity of the network in a self-consistent way. The transition rates of this model are also estimated form the spiking network simulations.

#### *2.3.1. Markov chain analysis*

In this section we address the question whether we can identify a time scale such that a Markov model can describe the time dependent dynamics of the network sufficiently well. To that end, based on the joint time series of the two population spike counts, a matrix of transition probabilities for any possible jump from state (*i*, *j*) at time *t* to state (*i* - , *j* - ) at time *t* + *dt* was inferred. Starting from an arbitrary initial condition, a stochastic signal was then generated which was able to match the dynamics and statistics of the two populations PSTH sufficiently well (**Figures 1B–D**). However, as the power spectra of the stochastic implementation of the Markov model were not exactly identical to those of the spiking network simulations, even for the optimal step size *dt* (see next subsection), we conclude that the Markov property is only approximately satisfied. This approximation to the complicated dynamics of the system provides a simple model for the stochastic behavior of the network. It is important to check which properties of the system can be explained by a Markov model of the type discussed here. Some limits of the model due to this assumption are described in the results and will be further analyzed in the Discussion Section.

#### *2.3.2. Selection of the bin size*

To extract the activity of the two neuronal populations, we consider a variant of the Peri-Stimulus Time Histogram (PSTH) of the spike trains of the excitatory and inhibitory neuronal populations. To find a suitable time resolution such that the dependency of the dynamics on the past history is essentially reduced to the most recent time bin, we looked at the autocorrelation functions of the PSTH of each population. As commonly done, we considered the first zero-crossing of the autocorrelation function as a first estimate for the time scale of the dynamics. As **Figure 1A** shows, however, the two populations have slightly deviating time scales. In the range of these two time scales, we explored a set of bin sizes that were integer multiples of the temporal resolution for the spiking network simulations (0.01 ms). For each value, we reconstructed the transition probabilities from the data, as explained in Section "Markov chain analysis." Using these transition probabilities, a stochastic signal for both populations was generated and the power spectra of the signals were compared with those of the spiking network simulations extracted at the same bin size. We found that a bin size of 0.1 ms provides power spectral densities best matching to those obtained from spiking network simulations for both populations (**Figure 1D**). This yields a reasonable time basis for the dynamics, based on numerical experiments, in line with the Markov assumption. Due to the discrete nature of spike counts in each bin, the dynamics in this paper is analyzed in discrete time, with a discrete noise model corresponding to the finite size of both neuronal populations.

## *2.3.3. Two state Active-Refractory Markov model*

The goal is to derive a Markov model the parameters of which can be interpreted in terms of neural dynamics, and that captures the finite size effects. The model should also capture the properties of the spiking network simulation on the level of mean-field as well as the transient fluctuations. We came up with a twostate Markov model, where each neuron is assumed to generate its spikes independently of the other neurons, given the input from the rest of the network. Moreover, it is assumed that each neuron's membrane potential falls into one of two classes: close to threshold (active state) or far away from threshold (refractory state). Hence, we call this process the "Active-Refractory Markov" (ARM) model. A schematic of the neuron model is depicted in **Figure 2**. Transitions from the active to the refractory state can be of two different types. Either the membrane potential of a neuron decays due to the membrane leak [Equation (1)], or the neuron receives inhibitory input spikes, or the neuron fires a spike itself and the membrane potential is reset. The former is described by the β branch; the latter is due to the γ branch in the model (**Figure 2**). More specifically:

In general, α, β, and γ are nonlinear functions of the state (spike counts). In our model, it turns out that α is an exponential function of a linear combination of spike counts; β and γ are constants estimated from the data. Note that the definition of a state in the ARM model deviates from the definition in the Markov chain model described in the previous section. In the latter, states are just spike counts, and transition probabilities are

curves) corresponding to a spiking network simulation with *v*reset = 0 mV. The low frequency dynamics are not very well captured by the Markov model; however, the overall shape of the spectrum, as well as the frequency beyond which the power drops, are well preserved. There are several peaks in the low frequency regime, corresponding to slow network oscillations. **(D)** Power spectral density for a simulation with *v*reset = 10 mV. In contrast to the case of stronger reset, **(C)**, there are fewer peaks in the low frequency range.

also defined in terms of spike counts. In the ARM model, however, the actual states are the occupation numbers of the active and refractory pools. These can neither be directly observed nor inferred from numerical simulations of the spiking neuronal network (SNN), but for our analysis we found a way work around this problem. All neurons in the network under study have the same in-degree and, consequently, share the same input statistics. Therefore, the transition rates are assumed to be identical for excitatory and inhibitory neurons. Neurons are assumed to make their state transitions, in particular spike firing, independently of each other, given their inputs.

The network is comprised of an excitatory and an inhibitory population. A population of size *N* will have *A* neurons in the active state and *N* − *A* neurons in the refractory state. Considering identical and independent transition rates for all neurons, our approach was to take neurons out of a state based on binomial distributions with the number of neurons in each state corresponding to the occupation numbers. The transition probability is given by multiplying the transition rate with *dt*, assuming that this is a small number. This results in the following stochastic description of the system

$$
\Delta A(t) = A(t + dt) - A(t) = \Delta A^+(t) - \Delta A^-(t) \tag{2}
$$

In equation (2), *A*+(*t*) and *A*−(*t*) are the increment and decrement from the active pool *A* that indicate the number of incoming neurons from the refractory pool, and the number of outgoing neurons to the refractory pool, respectively. *A*−(*t*) has a binomial distribution with parameters *A*(*t*) and (β + γ )*dt* as the number of available neurons and the probability of selection, respectively. Therefore,

$$\varrho\left(\Delta A^{-}(t) = \mathbf{x}\right) = \binom{A(t)}{\mathbf{x}} \left( (\boldsymbol{\beta} + \boldsymbol{\nu})dt \right)^{\mathbf{x}} \left( 1 - (\boldsymbol{\beta} + \boldsymbol{\nu})dt \right)^{A(t) - \mathbf{x}} \tag{3}$$

where *A*(*t*) *x* = *Cx A*(*t*) is the binomial coefficient. In this model, the spike count *S* is determined by the number of neurons in the active state and the transition rate γ in two steps: First, the total number of outgoing neurons from the active pool is calculated from a binomial distribution with probability (β + γ )*dt*. This quantity is exactly equal to *A*−(*t*). Second, from *A*−(*t*), the neurons that are actually firing a spike are drawn with another binomial distribution with probability γ γ +β . Therefore,

$$p\left(\mathbf{S}(t) = \mathbf{z}\right) = C\_{\Delta A^{-}(t)}^{z} \left(\frac{\mathcal{V}}{\mathcal{V} + \beta}\right)^{z} \left(\frac{\beta}{\mathcal{V} + \beta}\right)^{\Delta A^{-}(t) - z} \tag{4}$$

A similar expression as the one given in equation (3), a binomial distribution with rate α, describes the number of incoming neurons to the active pool at time *t*.

$$p\left(\Delta A^{+}(t) = \mathbf{y}\right) = C\_{N-A(t)}^{\mathcal{Y}} (\alpha dt)^{\mathcal{Y}} (1 - \alpha dt)^{N-A(t)-\mathcal{Y}} \tag{5}$$

With *A* neurons in the active state and *N* − *A* neurons in the refractory state at time *t*, the dynamic equation describing the expected value of the dynamics of *A* is

$$\frac{d}{dt}\mathbb{E}[A(t)] = \alpha(N - \mathbb{E}[A(t)]) - (\beta + \gamma)\mathbb{E}[A(t)]\tag{6}$$

It is important to point out that the number of neurons in the active state, *A*(*t*), is not observable, as we record only spikes. To describe the dynamics of the network it is, therefore, easier to describe the temporal dynamics in terms of observables, like the number of spiking neurons *S*(*t*). As mentioned before, *S* is the integer number of spikes generated in the γ branch, where E[*S*] = E[*A*]γ *dt*. It is straightforward to rewrite equation (6) as a function of *S* and get an equation describing the temporal dynamics of the spike counts

$$\frac{d}{dt}\mathbb{E}[\mathcal{S}(t)] = \alpha(\mathcal{N}\boldsymbol{\chi}\,dt - \mathbb{E}[\mathcal{S}(t)]) - (\boldsymbol{\beta} + \boldsymbol{\gamma})\mathbb{E}[\mathcal{S}(t)] \tag{7}$$

For better readability, we will drop the E[. ] operator that indicates the expected values of *A* and *S*. As all neurons are assumed to have the same membrane potential dynamics, determined by equation (1), and receive the same number of inputs regardless of their identity, the two populations have identical transition rates and therefore equation (7) holds for both populations. The only difference between the two populations, however, is the total number of neurons included in each of them. We denote the number of excitatory neurons by *Ne* and the number of inhibitory neurons by *Ni*. Excitatory and inhibitory spike counts generated at time *t* are given by *Se*(*t*) and *Si*(*t*), respectively. A nonlinear regression analysis of equation (7) applied to data from a spiking network simulation shows that α is an exponential function of recent spike counts, and that β and γ are well approximated by constant rates (see the Result Section for more details). Therefore, the two-dimensional mean-field has the following dynamics

$$\begin{cases} \dot{S}\_{\epsilon}(t) &= \exp\left(c\_{0} + c\_{1}S\_{\epsilon}(t) + c\_{2}S\_{i}(t)\right) \left(N\_{\epsilon}\gamma \, dt - S\_{\epsilon}(t)\right) \\ &- (\beta + \gamma \,)S\_{\epsilon}(t) \\ \dot{S}\_{i}(t) &= \exp\left(c\_{0} + c\_{1}S\_{\epsilon}(t) + c\_{2}S\_{i}(t)\right) \left(N\_{i}\gamma \, dt - S\_{i}(t)\right) \\ &- (\beta + \gamma \,)S\_{i}(t) \end{cases} \tag{8}$$

All the unknown parameters of the ODE system (8) can be estimated from the vector field extracted from simulated data. For this purpose, the regression analysis was performed with Python, using general purpose least-square optimization available in the SciPy library (Jones et al., 2001).

The exponential shape of α can be qualitatively justified. As pointed out in previous studies (Brunel and Hakim, 1999; Ricciardi et al., 1999; Brunel, 2000), the time dependent distribution of the membrane potential of a typical integrate and fire neuron, under the assumption of stochastic input and small PSPs, follows the Fokker-Planck equation with a drift and a diffusion term. In general, these two terms are functions of the recurrent activity of the network and therefore result in a nonlinear partial differential equation. However, under the assumption of Gaussian white noise input that is independent of the activity of the network, the steady state solution of equation (9), with appropriate boundary conditions, characterizes the stationary distribution of the membrane potential of a typical leaky integrate-and-fire (LIF) neuron. The dynamics of the distribution of the membrane potential is therefore

$$\tau \frac{\partial}{\partial t} p(\nu, t) = -\frac{\partial}{\partial \nu} \left[ \left( -\nu + \tau \sum\_{k} r\_{k}(t) f\_{k} \right) p(\nu, t) \right] \qquad (9)$$

$$+ \frac{1}{2} \left[ \tau \sum\_{k} r\_{k}(t) f\_{k}^{2} \right] \frac{\partial^{2}}{\partial \nu^{2}} p(\nu, t)$$

$$= -\frac{\partial}{\partial \nu} \left[ (-\nu + \mu) p(\nu, t) \right] + \frac{\sigma^{2}}{2} \frac{\partial^{2}}{\partial \nu^{2}} p(\nu, t)$$

where *rk*(*t*) is the firing rate of the pre-synaptic source *k* and *Jk* is the corresponding PSP to source *k*. μ and σ are the first and second moments, respectively, of the external input to the neuron. Equation (9) is based on the conservation of the probability flux of the membrane potential and could be rewritten in the following form:

$$
\pi \frac{\partial}{\partial t} \mathfrak{p}(\nu, t) = -\frac{\partial}{\partial \nu} \Phi(\nu, t), \tag{10}
$$

where (*v*,*t*) represents the probability flux and shows the probability mass crossing any arbitrary *v* per unit of time *t*. The stationary time independent solution of *p*(*v*) without considering refractoriness and with appropriate boundary conditions at threshold θ, *v*reset and *v* = −∞ is (Gerstner and Kistler, 2002)

$$p(\nu) = \begin{cases} \frac{c\_1}{\sigma} \exp\left(-\frac{(\nu-\mu)^2}{\sigma^2}\right) & \nu < \nu\_{\text{reset}}\\ \frac{c\_2}{\sigma} \exp\left(-\frac{(\nu-\mu)^2}{\sigma^2}\right) \int\_{\nu}^{\theta} \exp\frac{(\mathbf{x}-\mu)^2}{\sigma^2} \,d\mathbf{x} & \nu\_{\text{reset}} < \nu < \theta \end{cases} (11)$$

Intuitively, active and refractory states of a neuron are linked to high and low membrane potentials, respectively. The probability to encounter a potential exceeding a certain value, therefore, should be related to the rate α that describes a transition from refractory to active in the ARM model. The transition rate α, therefore, is proportional to the time dependent flux of the membrane potential crossing the border between the active and the refractory state. For simplicity we can approximate the time (state) dependent flux with the steady state flux plus some fluctuations. These fluctuations are determined by the spike counts of the activity: the larger the excitatory spike counts and the smaller the inhibitory spike counts are, the larger the probability flux crossing the border between the active and the refractory state is. As in equation (10), is related to the voltage integral of *p*(*v*, *t*), the time dependent flux will be proportional to the local behavior of the Cumulative Density Function (CDF) of the membrane potential distribution. In our Result section, it is argued that for a wide range of the membrane potential between *v*reset and θ, the stationary Cumulative Density Function (CDF) is locally well approximated by an exponential function, matching the functions we saw in numerical simulations of SNNs based on the LIF neuron model. A formal derivation of the link between the LIF and the ARM model, however, is mathematically involved and beyond the scope of the present paper.

#### **3. RESULTS**

In this section, numerical results from the simulation of large but finite balanced random networks are illustrated, which all imply either the nonlinear dynamics of interactions, or the non-Gaussian and state dependent nature of the self-generated noise. A comparison between the spiking network simulation results and simulation of the Active-Refractory Markov (ARM) model is made to show the limits and strengths of the model.

#### **3.1. MARKOV CHAIN INFERRED FROM TIME SERIES**

A numerical study of the network shows that with a time bin of 0.1 ms the most essential features of the system under study are recovered by a Markov chain model where the probabilities for state transitions are all extracted from the data. The vector field obtained from the spike counts (**Figure 1B**) and the power spectral density of both populations modeled with a Markov chain (**Figure 1D**) are similar to those obtained from the simulation of the spiking network.

The autocorrelation function (ACF) of the excitatory and inhibitory spike counts (**Figure 1A**) reflect the memory of the system. It has been reported previously (Tetzlaff et al., 2012), and is confirmed again in our study, that the time constant of the decay is smaller for the inhibitory population, although the input to excitatory and inhibitory neurons are statistically the same. This behavior was hypothesized to be related to the negative feedback contribution of the inhibitory population in the large-scale dynamics of the network (Tetzlaff et al., 2012). To examine the Markovian nature of the dynamics, a good choice for the time bin of the histograms is the point where the autocorrelation has its first zero-crossing. For a simulation with a temporal resolution of 0.01 ms, this point is roughly 0.07 ms for the inhibitory population and 0.20 ms for the excitatory population, according to **Figure 1A**. Since a unique time scale for the model and simulations is needed, we chose 0.1 ms to construct the PSTH of both populations.

The low frequency power spectra of the populations are not very well captured by the Markov model. This may be caused by the fact that there is a dependence on the past spiking activity of the two populations, which cannot be reflected by the Markov process employed here. We wanted to find out which features of the system are nevertheless recoverable by a Markov process. The low frequency behavior of the system depends on the distance between *v*reset and the spiking threshold θ (**Figures 1C,D**). The bigger this distance is, the more peaks in the low frequency part of the system appear. Refractoriness can also change the shape of the power spectrum in these frequency ranges (Franklin and Bair, 1995; Mar et al., 1999; Spiridon and Gerstner, 1999). It is obvious that a Markov process due to its lack of memory cannot capture this phenomenon.

## **3.2. TRANSITION PROBABILITIES OF THE ACTIVE-REFRACTORY MARKOV (ARM) MODEL**

All neurons in the network under study are statistically the same. Therefore, identical transition probabilities α, β, and γ were imposed for excitatory and inhibitory neurons. A nonlinear regression to estimate the parameters of equation (8) given the excitatory and inhibitory spike counts and their temporal derivatives from the time series results in an exponential link function for α. The function β exhibits a slight negative dependency on the excitatory spike counts (**Figure 3C**-dots). Ignoring this does not visibly affect the simulation results of the model (data not shown). Therefore, it was assumed to be a constant parameter. An exponential function for α and constant parameters for β (solid lines in **Figure 3**) and γ result in *c*0 = −0.046,

*c*1 = 0.032, *c*2 = −0.152, β = 7.78 and γ = 0.325 (all quantities in spikes/ms) as the estimated parameters of the model in equation (8). Interestingly, *c*2, the coefficient of the inhibitory spike counts is roughly 5 times bigger in amplitude compared with *c*1. This was expected since IPSP = −*g* EPSP. In the stochastic implementation of the ARM model, the transition rate α is a function of random variables *Se* and *Si*. Therefore, the ARM model can be interpreted as a type of "doubly stochastic" point process.

As mentioned in the Methods Section, for a wide range of membrane potentials *v* the CDF of the membrane potential can be approximated by an exponential function (**Figure 4B**). For values of the membrane potential close to threshold, however, the CDF does not match the fitted exponential function very well. This sub-exponential behavior is also observed in the data (dots in **Figures 3A,B**). This similarity suggests that a better approximation for α might come from the integral of the analytical solution of equation (9), but a formal mathematical analysis of this idea is beyond the scope of this paper. Moreover, β can be assumed as the leak term in the dynamics of leaky integrate and fire model due to its role in taking the value of the membrane potential away from the membrane threshold when there is no other inputs to the neuron. Inhibitory inputs in general can influence this parameter as well, however, for the sake of simplicity, we assumed it to be a constant. γ is a constant rate determining the number of spiking neurons at any given time interval *dt*. This term is proportional to the outgoing probability flux from the threshold θ in equation (11).

In comparison to the well-known Wilson-Cowan model (Wilson and Cowan, 1972), the ARM model assumes an exponential, instead of a sigmoidal, function as a transfer function for low input level. The reason is that in the balanced network, the firing rates of the neurons are low and therefore the activity of the populations are far from saturation. The ARM model can be considered as a special case of Wilson-Cowan model that was suggested for the dynamics of fluctuations around the mean firing rates of the populations for time-dependent inputs. The advantage of our model is that it can generate the statistics of the noise from the mean-field dynamics of the system in a self-consistent way.

#### **3.3. STATE SPACE ANALYSIS OF SELF-GENERATED NOISE**

Each pair of excitatory and inhibitory spike counts that are observed in the same time bin define a state in the state space. Spiking network simulations show that the future evolution of states during a particular trajectory of the system is highly dependent on the current state of the network. Therefore, the increment or decrement of the spike counts, as well as the derivative estimated from this (which are basically the difference between the spike counts at successive points in time divided the time interval 0.1 ms) are state dependent (**Figures 5A,B**). For each state, the distribution is heavy-tailed (positively skewed) and has a higher variance in excitation (compare the color-bars of **Figures 6A,B**). For a particular value of inhibitory spike counts, the variance increases as excitatory spike counts increase. Also, for a particular value of excitatory spike counts, increasing inhibitory spike counts results in less variability of the derivative of the spike counts for both excitatory and inhibitory populations. In other words, the more net excitation is in the system, the higher is the variance of the derivatives of the population spike counts. This statistical property of the network dynamics indicates that the noise model cannot be that of additive Gaussian white noise, as otherwise the variability in derivatives would be identical for the entire state space. The reason is that for any stochastic system governed by an equation of type *x*˙ = *f*(*x*) + ξ , if the variance of *x*˙ is the same for the entire state space, ξ could be considered as an additive white noise where the variance of the noise ξ is not state dependent. Otherwise, the noise term has to be state dependent, or non-additive. Multiplicatively interacting point processes are an explicit model for the interaction between neurons (Cardanobile and Rotter, 2010), with non-additive and non-Gaussian noise.

**hint for state probability transition** *α***. (A)** Stationary PDF of the membrane potential of a leaky integrate and fire neuron. The value of the membrane potential that separates the active from the refractory state is somewhere

In the ARM model the transition probabilities are generally state dependent. Therefore, each state has a different distribution of the spike rate increments (derivatives). The model shows the same pattern of state dependent variance, however, the variance is not as high as that of spiking network simulations (**Figures 6C,D**). Particularly, for the excitatory spike counts the variance generated in the model is smaller by a factor of 10 (**Figures 6A,C**). However, this fact does not affect the distribution of the spike counts in the populations drastically (**Figure 9**) and the normalized correlation functions (**Figure 8**) are still recovered with a high accuracy.

**Figure 5** shows that for some states in the state space, the distribution of the derivatives of the spike counts is reproduced by the ARM model with very high fidelity. Specifically for excitation, however, the model does not seem to generate enough variability in the excitatory spike counts. In **Figures 5C,D** the overlap between the distributions of the derivatives generated from SNN and ARM in the entire state space are represented by a number between 0, indicating no overlap, and 1, corresponding to a complete overlap, respectively. To calculate the overlap, first the state dependent distribution of the derivatives both for the model and the simulation data were normalized. Then the two distributions were compared, and for each possible value of the derivative, the minimum between the two distributions was determined. The integral over the minima is a number between 0 and 1 corresponding to the overlap. It is clear that the performance of the model in terms of variability of the inhibitory activity is very good (**Figure 5D**), however, the model is not as good in generating large enough variability in the derivative of the excitatory population. This might be related to the larger number of excitatory

is observed in the case of inhibition. **(C)** Variance of spike count derivatives for the excitatory population in the stochastic implementation of the ARM model. The same pattern as **(A)** is observed but with a 10 times reduced variance. **(D)** Variance of spike count derivatives for the inhibitory population in the ARM model. The scale of the state dependent variance of noise is the same as in **(B)**.

neurons compared to inhibitory ones, and therefore the bigger influence of the pairwise correlations among excitatory neurons that is ignored in the ARM model. The ARM model is based on the assumption that neurons, given the network is in a certain state, perform spike transition independently. This leads us to use binomial distributions to describe the transitions between the two neuronal states. This results in a linear scaling of the variance with the population size. As shown in **Figures 6C,D** (note the different scales of the color bars), the variance of the derivatives in both populations differ by a factor of 4, which is exactly the ratio of the excitatory and inhibitory population sizes. However, in the spiking network simulation, the variance of the derivatives for the excitatory population is much bigger than the inhibitory population (**Figures 6A,B**). We conclude that due to the linear relationship between the variance and the population size, the ARM model systematically underestimates the variance of the activity for the larger subpopulation in the network. This is a drawback of the model and in the discussion section we will suggest ways to overcome this problem. In the next section it is shown that the portion of the state space visited in a stochastic simulation of the model is less spread, and we attribute this fact to the reduced variability in the increments of the excitatory spike counts in the model.

### **3.4. NONLINEAR ISOCLINES**

The state space of the system reconstructed from the simulated data (vector field shown in **Figure 7**) is a good representation of the mean-field dynamics that represents smooth transition of the average fluctuations toward the fixed point. Starting from an arbitrary initial condition in the state space, the mean-field dynamics leads the trajectory to a stable fixed point. However, due to recurrent activity and the finite size of the network, the trajectory driven by the mean-field is continuously perturbed and the result is a quasi-stochastic signal that only on average follows the

**FIGURE 7 | Nonlinear isoclines and dynamic flow (vector-field) in the state space of excitatory and inhibitory spike counts. (A)** Left: Mean of the excitatory spike count derivatives; three levels of contour lines, −200, 0, 200 are shown in black. The isoclines are clearly nonlinear. Right: Mean of the inhibitory spike count derivatives; the contour lines are shown for −50, 0, 50. For inhibition, the nonlinearity is not dominant. **(B)** Mean of the derivatives of the excitatory and inhibitory spike counts in the ARM model. The isoclines

are similar to those in a spiking network simulation. **(C)** Vector field extracted from a numerical simulation of spiking neuronal network (SNN; left) and simulation of the ARM Markov model (ARM); (middle: analytical model, right: stochastic simulation of the model). Parameters of the model were chosen such that the vector fields on the left and in the middle are identical. There is a good match between the vector fields and the nullclines in the simulation and in the model.

mean-field dynamics. In other words, the vector field describes how a trajectory evolves on average.

The vector field has one component for each population representing the mean increment (or derivative) of the spike counts, given the two-dimensional state of the joint population. The excitatory and inhibitory components of the average state dependent derivatives are shown in **Figure 7A**. Isoclines in the state space represent contour lines of a particular value of the derivative. In **Figure 7A**, the isoclines corresponding to values −200, 0, 200 of the excitatory spike count derivative, and the −50, 0, 50 isoclines for the inhibitory spike count derivative are depicted in black. The inhibitory isoclines seem to be a linear function of excitatory and inhibitory spike counts. The isoclines for the derivatives of the excitatory spike counts are non-linear with a negative dependency of inhibition on excitation when there are only few excitatory spikes and relatively more inhibitory spikes are available. The slope of the dependency becomes positive when the number of excitatory spikes increases. The nullclines (0 isocline) of the system are the solutions of the two dimensional mean-field equations. An immediate conclusion from the shape of the nullclines is that the dynamics of a finite size network of excitatory and inhibitory neurons is nonlinear. More theoretical evidence for the nonlinearity of the dynamics, based on the system characteristic equation analysis is provided in the Supplementary Material.

The ARM model with parameters estimated form the vector field of the spiking network simulation reproduces the flow with very high fidelity (**Figure 7C**, middle). The stochastic implementation of the ARM model generates a comparable vector field with a similar shape of the nullclines (**Figure 7C**, right). However, the spread of the excitatory spike counts is reduced and therefore, compared to the spiking network simulation results, the variance of the excitatory spike counts is smaller (**Figure 9B**). The low variance of the excitatory spike counts may be explained by the state dependent distribution of the derivatives of the excitatory spike counts which the model generates (**Figure 5A**, red curve compared to the blue). In comparison with the spiking network simulation, the variance of this distribution is reduced. We conjecture that if the excitatory spikes are not only generated from a binomial distribution, but are also correlated with the recently generated inhibitory spike counts, the variance of this distribution will increase. Another way of generating excitatory spike counts from the occupation number would be to consider the two dimensional log-normal distribution of the spike counts in the state space.

#### **3.5. CORRELATION FUNCTIONS OF THE TWO POPULATIONS**

The balanced random network with stationary input is a highly recurrent system. A high correlation between excitation and inhibition nevertheless results in a decorrelated input to single neurons, because the recurrent inhibitory input cancels the effect of recurrent excitation (Renart et al., 2010). In an inhibition-dominated network, this provides an input to each neuron such that the mean membrane potential is below threshold and only fluctuations of the input cause threshold crossing and therefore result in a very low rate and irregular spiking activity of a neuron in the network. The high correlation between excitation and inhibition manifests itself in the cross-correlation between the spiking activity of the excitatory and inhibitory population (**Figure 8C**). The ARM model is capable of producing a high correlation function between excitation and inhibition, and the time scale of the correlation function is the same as for the network simulations.

The autocorrelation function of the excitatory and inhibitory population activities in the simulation of the spiking neuronal network has a very sharp decay in the range of the time scale selected for the Markov model (0.1 ms). For the excitatory population, the autocorrelation is close to zero after 0.2 ms and the ARM model reproduces the same correlation function (**Figure 8A**). We conclude that the model is able to capture the temporal dynamics of the excitatory population on such a short time scale. The autocorrelation of the inhibitory spiking activity decays faster and exhibits a small undershoot after 0.1 ms. The ARM model shows the same time-scale of correlation decay, however, the undershoot is more pronounced (**Figure 8B**). The shorter time scale of the inhibitory auto-correlation function compared to the excitatory one is also represented in the power spectral density of the excitatory and inhibitory populations, illustrated in **Figures 1C,D**. The excitatory population concentrates most of its power in a relatively narrow low frequency band compared to the inhibitory population.

### **3.6. HEAVY TAILED DISTRIBUTION OF SPIKE COUNTS**

Some previous studies of the temporal dynamics of interacting populations assume a Gaussian distribution of the activity around the fixed point solution of a low-dimensional system (Kriener et al., 2008; Tetzlaff et al., 2012; Helias et al., 2013), however, our numerical study shows that a log-normal distribution provides a good fit to the spike counts observed in the the spiking neuronal network (SNN) simulation (**Figure 9A**). Note that the population spike counts are bounded from above (due to the finite size of the system), a property that is not reflected by the

log-normal distribution. Nevertheless, it was found to provide a good approximation for the range of data actually observed in our network simulations. Beyond the upper bound, the probability of spike counts is zero and the log-normal distribution is not valid. Therefore, strictly speaking, the distribution of such a bounded random variable is not heavy-tailed but we use this term for the positively skewed distribution of spike counts. The parameters μ and σ of the fitted log-normal distribution are 3.214, 0.664 and 1.938, 0.582 for the excitatory and inhibitory spike counts, respectively. The ARM model with the parameters extracted from the mean-field flow which was reconstructed from simulated data can capture this statistical property of the system (**Figure 9B**). The corresponding parameters of the ARM model are 3.405, 0.405 and 2.036, 0.539 for the excitatory and inhibitory spike counts, respectively. A log-normal distribution fits to the distribution of the spike counts generated in spiking network simulation and in the stochastic implementation of the ARM model very well, much better than possible alternative right-skewed distributions (gamma distribution, negative binomial distribution; data not shown). Insets of **Figure 9** show a fitted normal distribution to the envelope of the distribution of the logarithm of the spike counts. The fact that spike counts are integer numbers needs to be taken into account in the fitting process. We conclude from these results that the spike counts of the population activity embedded in a network do not follow Poisson statistics, as the mean and the variance of the spike counts are not identical.