Balanced neural architecture and the idling brain

A signature feature of cortical spike trains is their trial-to-trial variability. This variability is large in the spontaneous state and is reduced when cortex is driven by a stimulus or task. Models of recurrent cortical networks with unstructured, yet balanced, excitation and inhibition generate variability consistent with evoked conditions. However, these models produce spike trains which lack the long timescale fluctuations and large variability exhibited during spontaneous cortical dynamics. We propose that global network architectures which support a large number of stable states (attractor networks) allow balanced networks to capture key features of neural variability in both spontaneous and evoked conditions. We illustrate this using balanced spiking networks with clustered assembly, feedforward chain, and ring structures. By assuming that global network structure is related to stimulus preference, we show that signal correlations are related to the magnitude of correlations in the spontaneous state. Finally, we contrast the impact of stimulation on the trial-to-trial variability in attractor networks with that of strongly coupled spiking networks with chaotic firing rate instabilities, recently investigated by Ostojic (2014). We find that only attractor networks replicate an experimentally observed stimulus-induced quenching of trial-to-trial variability. In total, the comparison of the trial-variable dynamics of single neurons or neuron pairs during spontaneous and evoked activity can be a window into the global structure of balanced cortical networks.


INTRODUCTION
The rich structure of neural firing patterns provides ample challenges for the systems neuroscience community. The complexity of in vivo neural responses depends, in part, on detailed single neuron biophysics, synaptic dynamics, and network interactions. To make sense of such complexity, it is often necessary to treat brain activity as probabilistic, similar to how statistical physics treats large ensembles of particles. In this spirit, collecting the responses of a neuron (or a population of neurons) over many trials of an experiment permits a statistical characterization of neural activity. A widespread observation is that the trial averaged response of a neuron's spike output is stimulus or action tuned (Decharms and Zador, 2000), and this has grounded many theories of neural coding (Dayan and Abbott, 2001). However, only considering the trial averaged response of a neuron glosses over many complex response dynamics that may give further insight into neural function.
Large trial-to-trial variability of the neural response is a general characteristic of cortical dynamics (Britten et al., 1992;Shadlen and Newsome, 1994;Averbeck et al., 2006;Cohen and Kohn, 2011;Ponce-Alvarez et al., 2013). In stimulus or task evoked states, spiking responses are well modeled as an inhomogeneous Poisson process, whose firing rate dynamics can be obtained from trial averaged responses. In this case, the response variability is equal to the trial averaged evoked response, and little is learned about cortical processing by considering neural variability. Neural activity during spontaneous dynamics (dynamics in the absence of a driving stimulus or action) is more complex (Ringach, 2009). In particular, the trial-to-trial spike variability during spontaneous states is larger than that predicted by a Poisson process (Churchland et al., 2006(Churchland et al., , 2010, and exhibits fluctuations over a range of timescales (Smith and Kohn, 2008). The state of cortex during spontaneous conditions is both an influence (Arieli et al., 1996;Fukushima et al., 2012) and a reflection of stimulus evoked activity (Tsodyks et al., 1999;Luczak et al., 2009;Luczak and MacLean, 2012). Thus, understanding the mechanics behind spontaneous activity may give new insight into the structure and function of cortical circuits.
Recent theoretical studies have shown that cortical models with balanced excitation and inhibition and clustered excitatory connectivity capture the trial-to-trial variability of both spontaneous and evoked conditions (Deco and Hugues, 2012;Litwin-Kumar and Doiron, 2012). In this article, we present a unified view of spontaneous dynamics in a variety of neural architectures. The core mechanic is that networks with a large number of stable states allow spontaneous fluctuations to stochastically "sample" the various states (Goldberg et al., 2004), leading to high spiking variability in spontaneous conditions. For stimuli to quench spiking variability, such models require that the recurrent architecture that supports network metastability be coherent with feedforward afferent projections. Under this assumption, we predict that pairwise spiking correlations as measured in the spontaneous state will be related to signal correlations in evoked states. Our article extends the framework of balanced excitation and inhibition Sompolinsky, 1996, 1998;Renart et al., 2010) to more structured network architectures, leading to a theory of variability that is consistent with spontaneously active cortex and its relation to the evoked cortical response.

NEURAL VARIABILITY IN BALANCED EXCITATORY AND INHIBITORY NETWORKS
Cortical circuits have recurrent excitatory and inhibitory connections that are sparse and random ( Figure 1A). As a first model, we consider no underlying global cortical architecture and rather assume wiring that obeys a pairwise independent model (meaning that a single connection between two neurons occurs with a probability that is independent of the rest of the network). Despite sparse wiring, the high number of cortical neurons ensures that a postsynaptic cell receives a large number (∼10 3 ) of synaptic inputs from other cortical neurons (Binzegger et al., 2004). Furthermore, the membrane potential deflection from an individual excitatory or inhibitory input can be on the order of 1 mV (Lefort et al., 2009). The convergence of a high number of large amplitude excitatory or inhibitory inputs could overwhelm a postsynaptic cell. However, in practice, the firing rates of cortical neurons are neither negligible nor driven to saturation (Hromádka et al., 2008). To explain this seeming contradiction between anatomy and physiology, Newsome (1994, 1998) proposed that large excitation and inhibition effectively balance one another (on average), so that the net mean input is small and output rates can be moderate ( Figure 1B). This conjecture has experimental support from in vivo whole cell recordings from cortical neurons, in which the magnitudes of isolated excitatory and inhibitory inputs are large yet balanced when neurons are near their resting state (Haider et al., 2006).
A key feature of balanced cortical networks is that spiking activity is dominated by the fluctuations in synaptic input (Shadlen and Newsome, 1994). These fluctuations can be sizable, since synaptic strengths are relatively large (only a small fraction of afferent inputs are necessary to drive the voltage to its threshold value). Synaptic input fluctuations provide a mechanism for the large dynamic and trial-to-trial variability of spiking activity that is characteristic of cortex ( Figure 1B, top). Using techniques from statistical physics, van Vreeswijk and Sompolinsky (1996Sompolinsky ( , 1998 provided a formal link between balance in recurrent networks and spiking variability. Specifically, they considered the statistics of the synaptic input into a representative cell in the network that received on average K inputs from recurrent excitatory and inhibitory wiring. With minimal assumptions, the mean M and variance V of the input obey: Here J α and r α are the synaptic strength and firing rates of the α ∈ {E, I} population and I ext is an external drive to the network. For a reasonable theory we require that M and V to be independent of system size for large K. For the input variance V to neither vanish nor explode as K → ∞, we must scale J α = j α / √ K where the parameter j α ∼ O(1). If we let the external input scale as I ext = i ext √ K, then we have the following expression for M: For the mean input M to neither vanish nor explode as K → ∞, the interior of the bracketed term must, as a whole, scale as O(1/ √ K). This only occurs when the recurrent and external excitation are balanced by recurrent inhibition, and with i ext > 0 this can occur robustly over a wide range of network parameters.
Using this mean field approach van Vreeswijk and Sompolinsky showed that balanced networks were dynamically stable over a wide range of parameters, with asynchronous spiking dynamics (Figure 1C), and single neuron statistics consistent with Poisson-like variability (Figure 1C, inset). An important component of this variability is that it is self-generated through network interactions, making it akin to high dimensional chaos. Such chaotic network-based variability has in vivo experimental support (London et al., 2010), and been the focus of recent theoretical investigations (Banerjee et al., 2008;Wolf, 2010, 2012). Furthermore, the fact that the variability is self-generated shows that the high degree of variability of cortical dynamics may be an intrinsic property of a cortical circuit, rather than being inherited from external sources to the network.
In these analyses it is often assumed that connectivity is sparse. However, later work by Renart et al. (2010) relaxed this assumption and showed that asynchronous activity could be achieved in densely coupled balanced networks. More recent work has discussed the impact of balance in working memory (Boerlin and Denève, 2011;Lim and Goldman, 2013) 1 , and the formation of orientation selectivity in visual cortices that lack functional topography (Hansel and van Vreeswijk, 2012). In total, balanced networks have been a successful framework to probe the mechanisms through which cortical dynamics are irregular, both in time and over repeated presentations of a fixed input.

SPONTANEOUS AND EVOKED DYNAMICS IN METASTABLE POPULATION DYNAMICS
The spike patterns from model neurons in balanced networks capture the irregular and asynchronous spike dynamics of evoked cortical response. However, balanced networks fail to capture the dynamics of cortical networks in spontaneous conditions (i.e., when the cortex is "idling"). In both sensory and motor cortices, the variability of spiking activity during spontaneous conditions is larger than that of a homogeneous Poisson process of the same firing rate (Churchland et al., 2006(Churchland et al., , 2010. A useful simplification is to treat a spike train in spontaneous conditions as a "doubly stochastic" process (Byron et al., 2009;Churchland and Abbott, 2012;Ponce-Alvarez et al., 2013), with one process modeling fast Poisson-like spike discharge and the other process capturing a slow fluctuation in the firing rate. Balanced networks of model neurons with independent wiring do not show the slow firing rate variability, and rather only capture the fast spiking variability ( Figure 1C). In this section we show how higher order architectural structure in excitatory-excitatory connectivity allow balanced networks to replicate the variability reported in both spontaneous and evoked conditions.

Excitatory networks with clustered architectures
Over the past decade, several experimental groups have performed careful analysis of the microstructure in cortical networks. One clear finding is an architectural structure that is beyond that of a simple independent wiring model. In particular, small interconnected clusters of excitatory neurons are overrepresented in the cortex (Song et al., 2005;Perin et al., 2011). Cluster membership is often stimulus Ko et al., 2011), circuit input , or activity (Yassin et al., 2010) dependent. Motivated by these anatomical findings, we subdivided excitatory cells into small clusters, and let the wiring probability be higher for two cells within the same cluster, compared to cells that are in distinct clusters (Figure 2A top, Materials and Methods). The chaotic dynamics of the network (Figure 2A, middle) is very distinct from that of an unclustered balanced network ( Figure 1C). Specifically, cells within a cluster show periods of coordinated low and high firing activity, and the transitions between activity states occur over long timescales. These firing rate transitions are in addition to the spike time variability caused by the balanced state input fluctuations, so that spiking dynamics are "doubly stochastic." Indeed, the spike count Fano factor over trials (normalized measure of variability, see Materials and Methods) is significantly above unity (Figure 2A, bottom), consistent with experimental measurements across cortex (Churchland et al., 2010). When the neurons in a cluster are stimulated together they have a higher likelihood of being in a high firing rate state when compared to spontaneous conditions. Furthermore, recurrent inhibition permits only a fraction of the clusters to have high firing rates at any time, meaning that neurons in unstimulated clusters have a higher likelihood of being in a low firing rate state compared to spontaneous conditions. Thus, by stimulating a subset of the clusters, the network is biased to be in a specific firing rate configuration. This removes the trial-to-trial rate variability characteristic of spontaneous dynamics (Figure 2A, middle), leaving only the Poisson-like spike dynamics due to balanced excitation and inhibition. As a result the Fano factor is reduced to near one in evoked conditions (Figure 2A, bottom), again matching experiments (Churchland et al., 2010). In our analysis we have corrected for reductions in Fano factor that would be caused by an increase in firing rate. This was done through a meanmatching algorithm, originally presented in Churchland et al. (2010), that discards data in evoked and spontaneous conditions so that the firing rate distributions overlap (see Methods). Thus, the stimulus-induced reduction in Fano factor cannot be trivially due to an increase in the mean spike count that exceeds an increase in spike count variance, and is rather because of an active reduction in spike count variance.
The mean field analysis that was applied to unclustered balanced networks can be applied to clustered networks. Consider the mean synaptic input M given to a representative neuron. Split the recurrent excitation into the number of connections coming from within the neuron's cluster K in , and the remaining connections from outside the cluster K out = K − K in . This gives the mean input as: The Fano factor for a homogenous (i.e., unstructured) balanced network (see Figure 1) network is also shown (dashed lines). For the homogeneous networks a random subset of neurons were excited during the stimulation period. In all networks, the inhibitory connections were not selective for any excitatory topology, and rather connected to all targets with a uniform probability.

Frontiers in Computational
For excitatory cells that are within and outside the cluster we distinguish their connection strengths, j in E and j out E , and firing rates, r in E and r out E . The first term on the right hand side is O(1) when the large recurrent excitation from outside the cluster and external drive are balanced by large recurrent inhibition. This is identical to the balance condition of the unclustered network. The second term is a perturbation from the balanced state that is induced by the clustered architecture. For this term to neither vanish nor overwhelm the network we require that the number of inputs that drive a cell from within its cluster K in scale as √ K (van Vreeswijk and Sompolinsky, 2005). For our network we set K ≈ 800 so that K in ≈ 30, meaning that the excitation from within the cluster is a small fraction of the total recurrent excitation. Equivalently, for each neuron to be in a single cluster, we require the network to be divided into Np in EE /K in ≈ 50 clusters. Overall, this means that the network will have a large number of stable states that provide a rich repertoire of dynamics in spontaneous conditions. The interested reader can see a full exposition of the trial-totrial variability of balanced networks with clustered connections in Litwin- Kumar and Doiron (2012).
The above cluster-based mechanism for a stimulus-induced reduction in spiking variability suggests a general framework. Highly variable spontaneous dynamics should occur in balanced networks with a large number of symmetric stable firing rate configurations. When population driven fluctuations cause stochastic transitions between distinct rate configurations, individual neurons will inherit these transitions as a new source of long timescale firing variability. Stimuli can break the network symmetry and bias specific firing rate patterns, reducing trial-to-trial firing rate variability. As a proof of principle, we show stimulus-induced reductions in variability in balanced networks with two global architectures that are distinct from clustered networks, yet satisfy the core requirements of our framework.

Excitatory networks with ring and feedforward architectures
Networks with local excitation and global (or lateral) inhibition arranged on a "ring" (Figure 2B, top) have been well studied as models of cortical orientation tuning (Ben-Yishai et al., 1995) as well as working memory (Compte et al., 2000). Ring networks support neutrally stable solutions where elevated rates are localized in space on the ring, often labeled a "bump" of activity (Ermentrout, 1998). In the absence of a stimulus, the network has a continuum of stable bump solutions, and population fluctuations allow a slow timescale wandering of the bump (Figure 2B, middle). Over realizations of the network, the random initial position and stochastic wandering of the bump drive slow firing rate variability so that single neuron Fano factors are large (Figure 2B, bottom). When a specific region on the ring is stimulated it becomes globally stable and other regions lose their stability, so that a bump of activity appears in that region and is pinned in space. Consequently, the firing rate variability is removed and single unit Fano factors drop to levels below spontaneous conditions (Figure 2B, bottom).
Another well studied architecture in theoretical neuroscience is networks with dominant feedforward coupling ( Figure 2C top) 2 . Feedforward networks support the creation and propagation of synchronous waves of activity along the network (Kumar et al., 2010), as well as provide a mechanism for selective amplification of activity that produces long timescale dynamics (Goldman, 2009;Murphy and Miller, 2009). The symmetry in the network supports wave initiation at any position in the chain, so that population fluctuations spuriously cause and destroy wave dynamics over time ( Figure 2C, middle). Over trials of the network, these wave events drive single unit firing rate variability so that Fano factors are elevated (Figure 2C, bottom). Stimuli at a fixed point along the chain break the translational symmetry along the network and effectively pin activity in that region since global network inhibition prevents wave propagation. Thus, over repeated trials the wave induced dynamics are absent and the evoked Fano factor are below what is observed in spontaneous conditions (Figure 2C, bottom).
In total, networks with balanced excitation and inhibition with ring or feedforward architectures also capture high variability in spontaneous conditions and the reduction of variability in evoked states.

COHERENCE BETWEEN STIMULUS DRIVE AND RECURRENT ARCHITECTURES
In the balanced network with clustered excitation, we chose feedforward stimulus drive to be aligned with cluster membership (Figure 2). In other words, neuron pairs that belonged to the same cluster received coordinated stimulation, while neuron pairs in different clusters did not. In this section we compare the recruitment of clustered network activity when stimulation respects cluster membership (cluster matched, Figure 3A, left) to when a random assortment of neurons are stimulated, with no relation to the underlying cluster membership (cluster interleaved, Figure 3A, right). The stimulus-response gain of the trial-averaged single neuron firing rate is higher for the cluster matched than for the cluster interleaved protocol ( Figure 3B). Furthermore, the cluster interleaved stimulation did not cause a reduction in the spike count Fano factor compared to the spontaneous state, in contrast to the cluster matched case ( Figure 3C). Thus, for the stimulus intensities studied, only cluster matched stimulation cause an evoked response that differs from the spiking statistics of the spontaneous state. This result is expected, as we describe below.
By design, the cluster of neurons driven by cluster matched inputs share more recurrence with one another than the neurons driven by cluster interleaved inputs. Because of this the positive feedback recruited by the stimulus is larger for the cluster matched protocol than for the interleaved protocol. Positive feedback increases response gain in spiking networks (Sutherland et al., 2009), allowing the response to weak inputs to be amplified. Furthermore, when the stimulation is interleaved, each cluster receives inputs to only a fraction of its neurons, distributed roughly evenly over the clusters. Thus, stimulus interleaved stimulation does not bias any one firing rate configuration more than any other, and the mechanism for firing rate variability is left unaffected. Hence so that the Fano factor remains high throughout stimulation. In total, these results show that a coherence between how the stimulus is distributed over the network and the global recurrent architecture is needed for evoked spiking variability to differ from those in the spontaneous state.

RELATING SPIKING CORRELATIONS IN THE SPONTANEOUS AND EVOKED STATES
Balanced networks with unstructured connectivity have a stable asynchronous state (Renart et al., 2010). As such, in these networks the spike count correlation coefficient from pairs of neurons is on average close to zero, and for large networks (N > 1000) the likelihood of correlation coefficients above 0.2 is negligible. The clustered network is also balanced and, similar to the unclustered network, has a density of pairwise correlation that is approximately Gaussian with a small spread about a near zero mean ( Figure 4A, black curve). While asynchronous dynamics are consistent with some cortical recordings (Ecker et al., 2010;Renart et al., 2010), there are many studies that show an average positive spike count correlation, with some pairs having significant (>0.2) correlations (Cohen and Kohn, 2011). However, in these experimental studies the correlation is measured from evoked, and not spontaneous, states.
The trial-to-trial co-variability of pairwise response in evoked states is termed "noise" correlations (Averbeck et al., 2006), to be distinguished from correlations in the trial-averaged activity due to the signal preferences of the neuron pair, i.e., "signal" correlations. In the clustered network we have assumed that stimulus preference and cluster membership are related (Figure 3). When pairwise correlation is conditioned on each neuron of the pair being a member of the same cluster, the mean correlation is far from zero, and the density has a heavy tail at positive correlations ( Figure 4A, green curve). This is due to neurons that are members of the same cluster being subject to the same firing rate fluctuations in spontaneous conditions, manifesting as a positive correlation coefficient. Thus, clustered architectures in balanced networks provide a mechanism for positive spike count correlations between select pairs. Further, because of the large number of clusters in the network, the number of pairs that are in the same cluster is far less than the number of total pairs in the network. This fact allows the clustered architecture to support a rough asynchronous state across the entire network, a requirement of balanced solutions (Renart et al., 2010).
The combination of these results suggests that a pair of neurons in the same cluster will have both a net positive correlation and an assumed common feedforward stimulus drive. This produces a clear and testable prediction: the signal correlation between pairs of neurons should be positively related to their correlation in the spontaneous state. This is certainly the case in our network simulations ( Figure 4B); however, this has yet to be tested from in vivo data (to our knowledge).

SPONTANEOUS ACTIVITY "SAMPLES" THE RANGE OF EVOKED STATES
Our analysis of spontaneous and evoked network activity has focused on its reflection in single neuron (Figures 2, 3) or pairwise (Figure 4) spike count statistics. However, the high correlation for pairs of neurons within the same cluster suggests that an analysis of the population dynamics of pairs of clusters will provide a network perspective on the dynamics of spontaneous and evoked activity.
In the spontaneous state, cluster activity randomly transitions between low and high firing rate regimes. The recurrent inhibition regulates the overall network activity such that if neurons in one cluster are firing at high rates then neurons in another cluster are likely firing at low rates. This is seen in the joint activity of a pair of clusters in the network, since their joint firing rate trajectory shows an anti-correlation ( Figure 5, gray curve). However, when a stimulus is given that drives either one of the clusters, the joint firing rate trajectory is effectively clamped so that the driven cluster's firing rate is high and the other cluster's is low. In the space of joint firing rates this collapse of network dynamics is such that the evoked state is a subset of the spontaneous state (Figure 5, red and blue curves compared to the gray curves).

STIMULUS QUENCHED VARIABILITY IN CHAOTIC FIRING RATE NETWORKS
We have presented a model for spontaneous cortical dynamics with high variability that is based on either a discrete (Figure 2A) or continuous (Figures 2B,C) attractor structure. Attractor dynamics require that recurrent cortical architecture provides a large set (or continuum) of stable network states. Such a constraint confers a fragility to the model dynamics, with unstructured heterogeneity in network architecture causing a collapse of network dynamics to only a few stable states. This criticism of attractor dynamics is often countered by assuming appropriate homeostatic mechanisms which provide a robustness to the attractor structure (Renart et al., 2003;Vogels et al., 2011). However, more research is needed to properly understand how attractor dynamics depends on the combination of homeostasis and rich spontaneous dynamics. An alternative model of cortex is a recurrent excitatoryinhibitory network built from phenomenological firing rate models (Ermentrout and Terman, 2010). When recurrent coupling is strong, these models show chaotic solutions (Sompolinsky et al., 1988;Sussillo and Abbott, 2009;Rajan et al., 2010), where firing dynamics are trial-to-trial variable due to the extreme sensitivity of the state of the network at the time of stimulation. In contrast to attractor networks, these networks are not dependent on any fine structured architecture, so that the chaotic solutions are robust to parameter variations. However, by design these models consider only output firing rates, and we must assume a suitable mechanism for how spike trains are generated from a rate dynamic. It is natural to consider the rate solutions as generating an inhomogeneous Poisson process. In this case, the trial variable chaotic dynamics of the rate model would produce spike trains that are "doubly stochastic" point processes, with an expected Fano factor larger than unity (Churchland and Abbott, 2012). Rajan et al. (2010) observed that global, dynamic stimulation quenched internally generated chaotic dynamics in firing rate models. Thus, network firing rates in stimulated conditions lacked the dynamic variability characteristic of the spontaneous state. If a network produced spikes according to an inhomogenous Poisson process with firing rates consistent with this model, then we expect a stimulus-induced reduction in the Fano factor, similar to that reported in our attractor model (Figure 2). However, it is still unclear if a network composed of actual spiking units, rather than interacting firing rates, is capable of reproducing the results of Rajan et al. (2010). This reproduction is not trivial because in most recurrent networks firing rates are stable ( Figure 6A, top) (van Vreeswijk and Sompolinsky, 1996Sompolinsky, , 1998Brunel, 2000;Renart et al., 2010), as opposed to chaotic.
Fortunately, Ostojic (2014) has recently shown how strong recurrent coupling in excitatory-inhibitory networks of spiking neuron models can destabilize firing rates, causing chaotic dynamics reminiscent of those predicted in firing rate networks (Figure 6A, bottom). This destabilization produces longer timescale spike train fluctuations reminiscent of a burst-like discharge that is absent in standard balanced networks (compare inset auto-covariance functions in Figure 6A). The network nevertheless shows roughly asynchronous spiking behavior and is a reasonable candidate model for rich spontaneous cortical dynamics. We implemented this spiking network in an effort to test the alternative theory for stimulus induced reduction in spike count variability proposed by Rajan et al. (2010).
The rate variability in the strongly coupled network produces high trial-to-trial variability in spiking dynamics in the spontaneous state, as characterized by a sizable Fano factor (∼4; Figure 6B, bottom panel before stimulation). When a static depolarizing input stimulus was applied to a half of the excitatory neurons, the firing rates of the stimulated neurons naturally increased (Figure 6B, top and middle panels). However, the mean matched Fano factor also increased, because the elevated firing rates recruited more burst behavior (Figure 6B, bottom panels). This is in contrast to the drop in variability seen in the attractor networks presented earlier (Figure 2), as well recorded spike trains from sensory and motor cortex (Churchland et al., 2010).
We also tested if periodic stimulation could quench the spiking variability in strongly coupled networks, as suggested by chaotic rate networks (Rajan et al., 2010) as well as data from rodent visual cortex (White et al., 2012). Periodic input was reflected in a periodic modulation of both firing rate and spike count Fano factor ( Figure 6C); however, the average Fano factor was not reduced by stimulation. This weak modulation of Fano factor, despite large changes in evoked rates, occurred over a range of stimulation frequencies ( Figure 6D). Thus, further work is needed to determine whether certain stimulus or network structures can reproduce the quenching of chaos seen in continuous firing rate networks.

DISCUSSION
The relation between spontaneous and evoked cortical dynamics is a popular topic of study (Ringach, 2009). Voltage sensitive dye imaging of visual cortex in cats has shown the spontaneous dynamics of large regions of cortex can mirrors evoked dynamics (Arieli et al., 1996;Tsodyks et al., 1999;Luczak and MacLean, 2012). Population recordings in the auditory system of rodents (Luczak et al., 2009) and primates (Fukushima et al., 2012) show that spontaneous dynamics contain a large repertoire of activity, with an evoked state being merely a subset of that repertoire. Our study contributes to this literature by identifying cortical network architecture which supports rich spontaneous dynamics, that nevertheless has a clear relation to evoked dynamics. We expand on our previous study (Litwin- Kumar and Doiron, 2012) and argue that the excess variability in spontaneous conditions (Churchland et al., 2010) is a characteristic of balanced networks with a large number of stable firing rate states.

BALANCED ATTRACTOR NETWORKS AND SCALING
The networks discussed in this paper are a merging of ideas from the balanced network literature Sompolinsky, 1996, 1998;Renart et al., 2010), which proposes models of trialto-trial variability, and the Hopfield network literature (Hopfield, 1982), which proposes models of neuronal assemblies. This connection has been made in several past studies (Amit and Brunel, 1997;van Vreeswijk and Sompolinsky, 2005;Renart et al., 2007;Roudi and Latham, 2007). Our contribution is specifically to study the implications of such attractor structures for trial-to-trial variability. Below, we discuss more thoroughly the plausibility of merging these ideas in the context of scaling arguments (see Renart et al., 2007 for an additional discussion of balance and scaling).
One requirement needed to merge these two lines of research was a specific scaling of the size of the neuronal assemblies. In balanced networks, if neurons receive K connections, then the connection strength J must be proportional to 1/ √ K. We considered the case of dense connectivity, so that K ∝ N. Further, the perturbation away from a homogeneous balanced network due to the attractor structure was set to be O(1), to prevent large attractors from dominating the network activity. Let the size of a neuronal assembly be K in neurons and define the coding fraction as K in /N. Then we must have: where J in E − J out E is how much stronger connections within the assembly are and p in − p out is how much more probable they are. There are three scenarios: 1. K in ∝ K, that is, each assembly is composed of a macroscopic number of neurons and the coding fraction is constant. This is the case of a classic Hopfield network (Hopfield, 1982). However, to satisfy Equation (1), we must have J in E p in − J out E p out ∝ 1/K. This presents a fine-tuning problem, as the J's themselves are proportional to 1/ √ K, much larger than the desired difference between J in E and J out E . Classic Hopfield networks avoided this problem as J was proportional to 1/K. 2. K in ∝ √ K, so the coding fraction is small. This was the case considered in our study and does not require fine-tuning as setting J in E = aJ out E for some constant a ensures Equation (1) is satisfied. In such networks, the number of neurons in a local assembly is smaller by a factor of √ K than the total number of inputs received. If stimuli are assumed to excite a few assemblies at once, this corresponds to a network exhibiting sparse coding, which is frequently encountered in cortex (Hromádka et al., 2008). 3. K in ∝ 1. This corresponds to a case of extremely sparse coding, as assembly size does not grow at all with K. Further, in this case we must have J in corresponding to a few extremely strong local connections. For these reasons, we do not consider this case.
Typically, spiking implementations of attractor networks fall into the second category (Amit and Brunel, 1997;van Vreeswijk and Sompolinsky, 2005;Renart et al., 2007), as was true for our study. Roudi and Latham (2007) took an alternative approach and studied the first category, differentiating between "background weights" that were O( √ K) and "foreground weights" (connections between neurons in the same assembly) that were O(1/K). In order to obtain stable assembly activation, these foreground weights needed to be tuned to within 6% in a network of 10,000 neurons, and it is unclear that such fine-tuning is plausible.
It is therefore an open question how to robustly merge balanced and attractor networks with a non-vanishing coding fraction. Unlike primary sensory regions, higher-order associative regions such as prefrontal cortex have a large degree of "mixed selectivity," with high coding fractions and assembly overlap (Rigotti et al., 2013). The fine-scale architecture of such regions has yet to be determined.

SYMMETRIC AND FINITE SIZE CLUSTERS
Despite their name, balanced networks are robust to changes in recurrent synaptic strengths J and admit a stable asynchronous solution over a wide range of network firing rates Sompolinsky, 1996, 1998;Renart et al., 2010). However, the mechanism for spontaneous dynamics presented in this study requires a strong assumption of symmetry in the network architecture. More to the point, rich spontaneous dynamics requires that each firing rate configuration be roughly equally stable. If this is not the case then some stable configurations will be overrepresented in the spontaneous state, and the population firing rates will tend to reside in those states. In the limit where one specific state is very stable, then the spiking variability in spontaneous and evoked dynamics would be similar. In our model, symmetric metastability is achieved by assuming that all clusters are of the same size, introducing a certain fragility into the framework. This is not an exact symmetry since quenched connection variability is present in any realization of the network. Nevertheless, our framework cannot support clusters of widely different sizes, as might be expected in real cortical networks. Future work should leverage powerful homeostatic mechanisms (Turrigiano and Nelson, 2000;Vogels et al., 2011) to prevent "winner-take-all" dynamics when the size and recurrence are heterogeneous across attractors.
The fluctuations in our model are an emergent feature of recurrent wiring, as opposed to being imposed by an external source. Balanced networks have a well defined solution as the number of inputs K → ∞ (van Vreeswijk and Sompolinsky, 1996Sompolinsky, , 1998Renart et al., 2010). However, in the clustered network the internal variability that promotes stochastic transitions between stable firing rate states requires K < ∞. Otherwise the time between firing rate transitions diverges, and the spontaneous state is no longer dynamically rich. Estimates for the number of neurons that are afferent to a cortical cell are imperfect, since anatomically and functionally defined architectures are often distinct, making a concrete value for K is difficult to obtain. Further, estimates for the size of a cluster are also difficult; however, (Perin et al., 2011) give an indirect number of ∼100 neurons per cluster, consistent with our network model. For models of cortical networks to improve, techniques for functional microcircuit measurements of large numbers of cells must be developed.

CORTICAL ARCHITECTURE DETERMINES VARIABILITY IN SPONTANEOUS AND EVOKED STATES
In our model, feedforward afferent projections are assumed to be coherent with recurrent architecture (Figure 3A). Without this feature stimulus drive does not quench trial-to-trial variability (Figure 3C), and our model would then be at odds with experimental data across sensory and motor cortex (Churchland et al., 2010). Also, if stimulus input recruited distinct populations than those recruited during spontaneous activity, then the patterning of spontaneous activity would not resemble that in evoked states, again in disagreement with cortical data (Tsodyks et al., 1999;Luczak et al., 2009). This circuit assumption is well justified from recent circuit recordings in mouse visual cortex Ko et al., 2011), where the stimulus correlation and connection probability of pairs of neurons are positively correlated with one another. Thus, a specific stimulus should drive neurons that are members of a specific cluster to a larger degree than other neurons in the network. Finally, such a wiring rule is consistent with standard Hebbian plasticity rules, where repeated co-activation of neuron pairs would engage plasticity rules that strengthen recurrent coupling between them. Understanding the trial-to-trial variability in evoked conditions is important for a complete description of a population code (Dayan and Abbott, 2001;Averbeck et al., 2006;Josić et al., 2009;Ponce-Alvarez et al., 2013). For this reason, a large number of studies focus on the trial-to-trial "noise" correlations of neuron pairs in evoked states (Cohen and Kohn, 2011). In both visual (Smith and Kohn, 2008;Ponce-Alvarez et al., 2013) and auditory (Rothschild et al., 2010) cortex, signal and noise correlations are themselves correlated, so that neurons that have similar stimulus preference, i.e., positive signal correlation, also have positive noise correlation. This suggests that the underlying circuits that establish stimulus preference also contribute to common fluctuations. Our model extends this idea to the spontaneous state, and predicts a positive relation between signal correlations and the correlations measured during spontaneous conditions. However, balanced networks (including the ones used in this study) fail to capture the positive average pairwise correlation that is consistent with some (Cohen and Kohn, 2011), yet not all (Ecker et al., 2010), data sets ( Figure 4A, black curve). A spiking network framework that internally generates non-zero mean correlation consistent with data is, to our knowledge, an open problem.
In sum, our model shows that a coherence between feedforward and recurrent architectures produces cortical dynamics that is consistent with a wide array of spontaneous and evoked data, and makes a clear and testable prediction for pairwise correlations in the spontaneous state.

SPIKING NETWORK SIMULATIONS
Neurons were modeled as leaky integrate-and-fire units whose voltages obeyed:V When neurons reached a threshold V th = 1, a spike was emitted and they were reset to V re = 0 for an absolute refractory period of 5 ms. The membrane time constant τ was 15 and 10 ms for excitatory and inhibitory neurons, respectively. The bias μ was chosen according to a uniform random distribution between 1.1 and 1.2 for excitatory neurons and between 1 and 1.05 for inhibitory neurons. While these values are superthreshold, balanced dynamics ensured that the mean membrane potentials were subthreshold (van Vreeswijk and Sompolinsky, 1998). In Figure 1, the non-dimensionalized voltages were transformed so that V th = −50 mV and V re = −65 mV. Synapses between neurons were modeled as differences of exponentials, and the total synaptic input to neuron i in population x was: where * denotes convolution. Here x, y ∈ {E, I} denote excitatory or inhibitory populations of N E = 4000 and N I = 1000 neurons each, J xy ij is the strength of synaptic connections from neuron j in population y to neuron i in population x, F y (t) is the synaptic filter for projections from neurons in population y, and s y j (t) is the spike train of neuron j in population y, a series of delta-functions at the time points where the neuron emitted a spike. F y (t) = 1 τ 2 − τ 1 e −t/τ 1 − e −t/τ 2 , with τ 2 = 3 ms for excitatory synapses and 2 ms for inhibitory synapses while τ 1 = 1 ms.
Connection probabilities p xy from neurons in population y to x were p EI = p IE = p II = 0.5, while p EE was on average 0.2. In clustered networks, excitatory neurons were partitioned into 50 assemblies of 80 neurons each and the connection probability was set to p EE in for neurons in the same assembly and p EE out for neurons in different assemblies. These were chosen so that p EE remained on average 0.2, but p EE in /p EE out = R EE = 2.5. For ring and feedforward networks, periodic boundary conditions were applied by identifying i = 4000 with i = 0. For ring networks, neurons i and j were said to belong to the same assembly if |i − j| < 40. For feedforward networks, connections from neuron i to neuron j were drawn with probability p EE in if i − j ∈ [−35, 45] and p EE out otherwise, thus biasing connections in one direction along the ring. If a connection from neuron j in population y to neuron i in population x existed, J xy ij = J xy (unless x, y = E and neurons were in the same assembly, then connection strength was multiplied by 1.9), otherwise J xy ij = 0. Synaptic strengths were J EE = 0.024, J EI = −0.045, J IE = 0.014, and J II = −0.057. These parameters, multiplied by 15 mV, would give the deflection of the membrane potential of the post-synaptic target, neglecting leak, in our dimensionalized units. When clusters of excitatory neurons were stimulated, stimulation was accomplished by increasing μ for neurons in those clusters by 0.07. Simulations were performed using Euler integration with a timestep of 0.1 ms.

SPIKE TRAIN STATISTICS
Spike train statistics were computed for excitatory neurons. We denote the spike times of neuron i as {t i1 , t i2 , t i3 , . . .}. We can then define neuron i's spike train: y i (t) = k δ(t ik ). The number of spikes emitted by the neuron between times t and t + t is The firing rate of a neuron over an interval (t, t + t) was defined as For the networks studied, firing rates and other statistics for the spontaneous state were calculated with t = 1.5 s to prevent effects due to initial conditions and t = 1.5 s. We also computed the Fano factor F i (t, t + t) for neuron i by evaluating F i (t, t + t) = Var (N i (t, t + t)) N i (t, t + t) where the expectations are over repeated trials of the same network with random initial conditions. When computing the Fano factor as a function of time relative to stimulus onset, we computed the mean-matched Fano factor described in Churchland et al. (2010) to control for changes in firing rate. Fano factors were computed over 100 ms windows.
We computed correlation coefficients for the spike counts of neuron pairs. The correlation between neurons i and j was given by where the covariances and variances were computed over overlapping windows within each trial and then averaged across trials. For spontaneous activity trial averaging was replaced with temporal averaging.