Imbalanced amplification: A mechanism of amplification and suppression from local imbalance of excitation and inhibition in cortical circuits

Understanding the relationship between external stimuli and the spiking activity of cortical populations is a central problem in neuroscience. Dense recurrent connectivity in local cortical circuits can lead to counterintuitive response properties, raising the question of whether there are simple arithmetical rules for relating circuits’ connectivity structure to their response properties. One such arithmetic is provided by the mean field theory of balanced networks, which is derived in a limit where excitatory and inhibitory synaptic currents precisely balance on average. However, balanced network theory is not applicable to some biologically relevant connectivity structures. We show that cortical circuits with such structure are susceptible to an amplification mechanism arising when excitatory-inhibitory balance is broken at the level of local subpopulations, but maintained at a global level. This amplification, which can be quantified by a linear correction to the classical mean field theory of balanced networks, explains several response properties observed in cortical recordings and provides fundamental insights into the relationship between connectivity structure and neural responses in cortical circuits.


Introduction
Information about a sensory stimulus is passed along a hierarchy of neural populations, from subcortical areas to the cerebral cortex where it propagates through multiple cortical areas and layers. Within each layer, lateral synaptic connectivity shapes the response to synaptic input from upstream layers and populations. In a similar manner, lateral connectivity shapes the response of cortical populations to artificial, optogenetic stimuli. The densely recurrent structure of local cortical circuits can lead to counter-intuitive response properties [1][2][3][4][5], making it difficult to predict or interpret a population's response to natural or artificial stimuli. This raises the question of whether there are underlying arithmetical principles through which one can understand the relationship between a local circuit's connectivity structure and its response properties.
In principle this relationship could be deduced from detailed computer simulations of the neurons and synapses that comprise the circuit. In practice, such detailed simulations can be computationally expensive, depend on a large number of unconstrained physiological parameters, and their complexity can make it difficult to pinpoint mechanisms underlying observed phenomena. In many cases, however, one is not interested in how the response of each neuron is related to the detailed connectivity between every pair of neurons. Relevant questions are often more macroscopic in nature, e.g. "How does increased excitation to population A affect the average firing rate of neurons in population B?" For such questions, it is sufficient to establish a relationship between macroscopic connectivity structure and macroscopic response properties.
One such approach is provided by the mean-field theory of balanced networks [6][7][8][9][10], which is derived in the limit of a large number of neurons and a resulting precise balance of strong excitation with strong inhibition. This notion of precise balance implies a simple relationship between the macroscopic structure of connectivity and firing rates, and balanced network models naturally produce the asynchronous, irregular spiking activity that is characteristic of cortical recordings [6,7,11,12]. However, classical balanced network theory has some critical limitations. While cortical circuits do appear to balance excitation with inhibition, this balance is not always as precise and spike trains are not as asynchronous as the theory predicts [13][14][15][16][17][18][19]. Moreover, precise balance is mathematically impossible under some biologically relevant connectivity structures [8][9][10], implying that the classical theory of balanced networks is limited in its ability to model the complexity of real cortical circuits.
We show that cortical circuits with structure that is incompatible with balance are susceptible to an amplification mechanism arising when excitatory-inhibitory balance is broken at the level of local subpopulations, but maintained at a global level. This mechanism of "imbalanced amplification" can be quantified by a linear, finite-size correction to the classical mean field theory of balanced networks that accounts for imperfect balance and local imbalance. Through several examples, we show that imbalanced amplification explains several experimentally observed cortical responses to natural and artificial stimuli.

The arithmetic of imprecise balance in cortical circuits
We begin by reviewing and demonstrating the classical mean-field theory of balanced networks and a linear correction to the large network limit that the theory depends on. A typical cortical neuron receives synaptic projections from thousands of neurons in other cortical layers, cortical areas or thalamus. These long range projections are largely excitatory and provide enough excitation for the postsynaptic neuron to spike at a much higher rate than the sparse spiking typically observed in cortex. The notion that excitation to cortical populations can be excessively strong has been posed in numerous studies and is typically resolved by accounting for local, lateral synaptic input that is net-inhibitory and partially cancels the strong, net-excitatory external synaptic input [2,6,[20][21][22][23]. Balanced network theory takes this cancellation to its extreme by considering the limit of large external, feedforward synaptic input that is canceled by similarly large local, recurrent synaptic input. In this limit, a linear mean-field analysis determines population-averaged firing rates in terms of the macroscopic connectivity structure of the network [6,7].
To demonstrate these notions, we first simulated a recurrent network of N E = 4000 excitatory (population E) and N I = 1000 inhibitory spiking neurons (population I) receiving synaptic connections from an "external" population (X) of N X = 4000 excitatory neurons modeled as Poisson processes. Cortical circuits are often probed using optogenetic methods to stimulate or suppress targeted neuronal sub-populations [24,25]. As a simple model of optogenetic stimulation of cortical pyramidal neurons, we added an extra inward current to all neurons in population E halfway through the simulation (Fig 1a). Neurons in the local population (E and I) were modeled using the adaptive exponential integrate-and-fire (AdEx) model, which A population of recurrently connected excitatory (red) and inhibitory (blue) spiking neuron models receive synaptic input from an external population (X; green) of Poisson-spiking neurons. Optogenetic stimulation of excitatory neurons was modeled by an extra inward current to the excitatory population at 5s. b) Spike rasters from 50 randomly selected excitatory (red) and inhibitory (blue) neurons from recurrent network. c) Average firing rate of excitatory (red) and inhibitory (blue) neurons in the recurrent network from simulations (light solid), from the balanced network approximation (Eq (3); solid dark) and from the corrected approximation (Eq (4); dashed). d) Mean synaptic currents to 200 randomly selected excitatory neurons in the recurrent network from external inputs (X; green), from the local population (E + I; purple) and the total synaptic current (black). Currents are measured in units of the neurons' rheobase (rheobase/C m = 10.5 V/s). e) Mean firing rates plotted against mean input currents to all neurons in populations E and I (gray dots) and a rectified linear fit to their relationship (black line). f) Mean firing rates from identical simulations without stimulation except the total number of neurons, N, in the recurrent network was modulated while scaling synaptic weights and connection probabilities so that $ 1= ffiffiffiffi N p (see Methods). Solid light curves are from simulations, solid dark from Eq (3), and dashed from Eq (4). accurately captures the responses of real cortical neurons [26][27][28]. Connectivity was random with each neuron receiving 800 synaptic inputs on average and postsynaptic potential amplitudes between 0.19 and 1.0 mV in amplitude. The recurrent network produced asynchronous, irregular spiking (Fig 1b), similar to that observed in cortical recordings [11,21,29,30]. Firing rates in populations E and I were similar in magnitude to those in population X and were increased by optogenetic stimulation (Fig 1c). As predicted by balanced network theory, local synaptic input (from E and I combined) was net-inhibitory and approximately canceled the external input from population X and artificial stimulation combined (Fig 1d).
A review of the mean field theory of balanced networks. To capture the notion that the net external synaptic input to neurons is strong, we define the small number, where K EX = p EX N X is the average number of external synaptic projections received by each neuron in E from all neurons in X, p EX is connection probability, and J EX is the synaptic strength of each connection. Specifically, J EX is the total postsynaptic current induced in a postsynaptic neuron in E by a single spike in a presynaptic neuron in X. Hence, 1/ quantifies the synaptic current that would be induced in each neuron in E (on average) if every neuron in X spiked once simultaneously. For the simulations in Using this convention, the mean synaptic input to each neuron in populations E and I from all sources can be written in vector form as A population of recurrently connected excitatory (red) and inhibitory (blue) spiking neuron models receive synaptic input from an external population (X; green) of Poisson-spiking neurons. Optogenetic stimulation of excitatory neurons was modeled by an extra inward current to 20% of the excitatory population at 5s. b) Spike rasters from 10 randomly selected ChR2-expressing and 40 non-expressing excitatory (red) neurons and 50 inhibitory (blue) neurons from the recurrent network. c) Average firing rate of ChR2-expressing excitatory neurons from simulations (light solid) and from the corrected approximation (Eq (4); dashed). d) Mean synaptic currents to 200 randomly selected ChR2-expressing excitatory neurons from external inputs (X; green), from the local population (E + I; purple) and the total synaptic current (black). e) Same as c, but for non-expressing excitatory neurons (red) and inhibitory neurons (blue). f) Same as d, but for non-expressing excitatory postsynaptic neurons. g) Same as c and d, but averaged over all excitatory neurons (expressing and non-expressing). Currents are measured in units of the neurons' rheobase (rheobase/C m = 10.5 V/s). Firing rates predicted by Eq (3)  where I = [I E I I ] T (superscript T denotes the transpose) is the vector of mean synaptic input to neurons in each population and similarly for their mean rates, r = [r E r I ] T . The rescaled external synaptic input, X = [X E X I ] T , is given by where r X is the average rate of neurons in population X and s/ is the strength of the inward current induced by optogenetic stimulation (s = 0 when stimulation is off). The recurrent and feedforward mean-field connectivity matrices are given by respectively where w ab = K ab J ab /(K EX J EX ) quantifies the relative number, K ab = p ab N b , and strength, J ab , of synaptic connections from population b to a. To achieve moderate firing rates when is small, local input, Wr, must be net-inhibitory and partially cancel the strong external excitation, X, in Eq (1). Balanced network theory [6,7] takes this cancellation to its extreme by considering the limit of large number of neurons, N = N E + N I , while scaling connection strengths and probabilities in such a way that $ Oð1= ffiffiffiffi N p Þ ! 0. Under this scaling, Eq (1) would seem to imply that mean synaptic currents diverge in the limit, but this divergence is avoided in balanced networks by a precise cancellation between external and recurrent synaptic input. To achieve this cancellation, firing rates must satisfy the mean-field balance equation, Wr þ X ¼ 0 in the large N limit, so that [6][7][8][9][10] Hence, balanced network theory provides a closed form, linear expression for firing rates in the large network limit. Generally speaking, the firing rate of a neuron depends nonlinearly on the mean and variance of its input current [22,31,32]. Notably, however, the derivation of the fixed point in Eq (3) did not require us to specify the exact form of this dependence. Instead, Eq (3) represents the unique fixed point firing rates for which synaptic currents remain bounded as N ! 1. More specifically, if Eq (3) is not satisfied as N ! 1 then kIk ! 1 (where k Á k is the Euclidean norm). The existence of this fixed point does not guarantee that it is stable. Precise, general conditions on the accuracy of Eq (3) for spiking network models are not known and the investigation of such conditions is outside the scope of this study. However, the approximation tends to be accurate in the N ! 1 limit whenever all eigenvalues of W have negative real part, the solution in Eq (3) is strictly positive, and inhibitory synaptic kinetics are sufficiently fast [6-10, 32, 33]. Indeed, Eq (3) provides a reasonable, but imperfect approximation to firing rates in our spiking network simulation (Fig 1c, compare light and dark solid). Balanced network theory has some critical limitations. Local cortical circuits are, of course, finite in size so the N ! 1 (equivalently ! 0) limit may not be justified. Moreover, excitation and inhibition in cortex may not be as perfectly balanced and spike trains not as asynchronous as predicted by balanced network theory [13][14][15][16][17][18][19]34]. More importantly, under many biologically relevant connectivity structures, precise cancellation cannot be realized so Eq (3) cannot even be applied [8][9][10]. We next review a simple, linear correction to Eq (4) that partially resolves these issues.
A linear correction to precise balance. A correction to Eq (3) can be obtained by considering non-zero, but this requires making assumptions on the relationship between neurons' input statistics and firing rates. A simple approximation is obtained by assuming that population-averaged firing rates, r, depend only on population-averaged mean inputs, I, yielding the fixed points problem r = f(I) = f([Wr + X]/) where f is the population-level f-I curve. When f is an increasing function over relevant ranges of I, this fixed point equation can be re-written as Wr þ X ¼ f À 1 ðrÞ: Hence, in strongly coupled networks ( small), the shape of f-I curves has a small effect on steady-state firing rates under such an approximation. Indeed, in the ! 0 limit, the f-I curve has no effect and firing rates are determined by Eq (3). This conclusion easily generalizes to the case where f also depends on the average temporal variance of neurons' inputs.
A simple case of this approximation is obtained by using a rectified linear approximation, r = g[I] + where [Á] + denotes the positive part. We fit such a function to the relationship between neurons' mean synaptic inputs and firing rates from our spiking network simulation (Fig 1e). Assuming that the average firing rates of all populations are positive, this rectified linear approximation produces a linear rate model [35] with mean firing rates given by solving Wr + X = /g r to obtain The AdEx neuron model used in our simulations has a nonlinear f-I curve (Fig 1e; gray dots) and its firing rate depends on all statistics of its input, not just the mean [26,36]. Nevertheless, the linear approximation in Eq (4) was accurate in predicting firing rates in our simulations (Fig 1c, dashed), outperforming the balanced network approximation from Eq (3). This can be explained by the fact that the balanced approximation in Eq (3) is already somewhat accurate and the linear approximation in Eq (4) corrects for some of the error introduced by imperfect balance, even though the true dependence of r on I is nonlinear.
To further investigate the relative accuracy of Eqs (3) and (4), we repeated the spiking network simulations from Fig 1a-1d while proportionally scaling the number of neurons (N E , N I , and N X ) in each population and scaling connection weights and probabilities in such a way that $ 1= ffiffiffiffi N p (see Methods). As predicted by balanced network theory, excitatory and inhibitory firing rates increased toward the limit in Eq (3) (Fig 1f, compare light and dark solid lines). The linear correction in Eq (4) tracks this increase in firing rates and is more accurate than the approximation in Eq (3), particularly for smaller N (Fig 1f,  dashed). It is worth noting that, in applying Eq (4) to obtain the dashed curve in Fig 1f, we fixed the value of g to the one obtained from the simulation in Fig 1a-1e. Hence, a single estimate of the gain yields an accurate approximation even under different parameter values.
The predictive power of Eq (4) in these examples is, of course, limited by the fact that it was only applied after estimating the gain of the neurons using firing rates obtained in simulations.
Moreover, highly nonlinear f-I curves could introduce additional error. However, the purpose of Eq (4) in this work is to provide a first-order approximation to and general understanding of firing rates in networks under which Eq (3) cannot be applied. For these purposes, Eq (4) is sufficient.

Imbalanced amplification under partial optogenetic stimulation
We next show that a more realistic model of optogenetic stimulation breaks the classical balanced state, providing a demonstrative and experimentally relevant example of imbalanced amplification and suppression that explains phenomena observed in recordings from mouse somatosensory cortex.
Firing rates are increased by stimulating fewer neurons. The model of optogenetic stimulation considered in Fig 1 is somewhat inaccurate since optogenetic stimulation of excitatory neurons is often incomplete. For example, only a fraction of cortical pyramidal neurons express the channelrhodopsin 2 (ChR2) protein targeted in many optogenetic experiments [3,24,37,38]. To more accurately model optogenetic stimulation, we modified the example above so the extra inward current was provided to a randomly chosen 20% of the excitatory neurons (Fig 2a), modeling ChR2-expressing pyramidal cells. This change produced surprising results. The ChR2-expressing neurons increased their firing rates by a larger amount than they did when all excitatory neurons received the current (Fig 2b and 2c; compare to Fig 1b and 1c). Hence, counterintuitively, stimulating fewer neurons actually amplifies the effects of stimulation on the targeted cells. In contrast, non-expressing excitatory neurons were suppressed during stimulation and inhibitory neurons increased their rates, but by a smaller amount than they did under complete stimulation (Fig 2e; compare to Fig 1c).
Similar effects were observed in experiments by Adesnik and Scanziani [3]. In that study, pyramidal neurons in layers (L) 2/3 of mouse somatosensory cortex (S1) were stimulated optogenetically, but only about 23% of the pyramidal neurons expressed ChR2. During stimulation, non-expressing L2/3 pyramidal neurons were suppressed and inhibitory synaptic currents increased, implying an increase in inhibitory neuron firing rates.
To understand these effects, we first extended the mean-field theory above to account for multiple subpopulations by defining to be the vector of average firing rates for the ChR2-expressing (exp), non-expressing (nexp) excitatory neurons and inhibitory (I) neurons. The vector of average input to the network is again given by Eq (1) where ; and q = 0.2 represents the proportion of neurons that express ChR2. Importantly, W is singular (i.e., not invertible), so classical balanced network theory fails for this example since Eq (3) cannot be evaluated. More specifically, it is impossible for I in Eq (1) to remain finite as ! 0 since there is no vector, r, such that Wr = −X. Intuitively, this can be understood by noting that expressing and non-expressing excitatory neurons receive the same local input on average (Fig 2d and 2f, purple), since local connectivity is not specific to ChR2 expression, but they receive different external input during stimulation (Fig 2d and  2f, green). Therefore, local synaptic input cannot simultaneously cancel the external input to both sub-populations, so the precise cancellation required by classical balanced network theory cannot be achieved (Fig 2d and 2f, black). A similar mechanism has been used to explain a lack of cancellation between positive and negative correlations in balanced networks [12,39].
Amplification in the nullspace: A general analysis. We now give a general analysis of network responses when W is singular. The example of partial optogenetic stimulation is then considered as a special case. If W is a singular matrix then only vectors, X, that are in the column space of W admit solutions to Wr + X = 0. The column space of W is defined as the linear space of all vectors, u, such that u = Wr for some r. The column space of a matrix, W, is the orthogonal complement of the nullspace of W T . We can therefore decompose where X 0 ¼ proj NðW T Þ X is the projection of X onto the nullspace of W T and X 1 = proj C(W) X is the projection onto the column space of W. Moreover, note that proj NðW T Þ Wr ¼ 0 since Wr is in the column space of W. Therefore, the projection of the total input onto the nullspace of Hence, the projection of the total synaptic input onto the nullspace of W T is Oð1=Þ whenever X has an Oð1Þ component in the nullspace of W T . Note that, despite the 1/ term in Eq (1), the total synaptic input, I, is Oð1Þ when balance is realized due to cancellation (as in Fig 1d). Hence, the singularity of W introduces large, Oð1=Þ synaptic currents where they would not occur if W was non-singular. In other words, external input in the nullspace of W T produces strong synaptic currents in the network. Importantly, this conclusion does not rely on any assumptions about neurons' f-I curves or other properties. This result is a fundamental property of balanced networks or, more generally, networks receiving strong feedforward input.
To understand the implications of this result on firing rates in the network, however, we must specify an f-I curve. We again consider the linear rate approximation quantified by Eq (4). Importantly, unlike Eq (3) from classical balanced network theory, the approximation in Eq (4) is applicable to this example because it accounts for imperfect cancellation between local and external inputs. Specifically, the regularized matrix, D − W, is invertible so Eq (4) can be evaluated even though Eq (3) cannot. The resulting firing rate solution from Eq (4) agrees well with spiking network simulations (Fig 2c and 2e). Hence, Eq (4) provides an accurate approximation to firing rates in networks to which classical balanced network theory is not applicable at all.
Eq (4) also provides a concise mathematical quantification of firing rates when W is singular. Namely, if X 0 , X 1 $ Oð1Þ then firing rates can be expanded as where r 0 is in the nullspace of W and r 0 , r 1 $ Oð1Þ. To derive this result, first note that Eq (4) can be rewritten as If X has components in the nullspace of W T then we can project both sides of this equation onto this nullspace to obtain where we again used the fact that proj NðW T Þ Wr ¼ 0 since Wr is in the column space of W. Since X 0 ¼ proj NðW T Þ X and D are assumed Oð1Þ, this equation is only consistent when r $ Oð1=Þ. We can therefore decompose r = (1/)r 0 + r 1 where r 0 , r 1 $ Oð1Þ. We next show that r 0 is in the nullspace of W. From Eq (7), we have Isolating the Oð1=Þ terms gives Wr 0 = 0 and therefore r 0 is in the nullspace of W. In summary, components of external input in the nullspace of W T partially break balance to evoke amplified firing rates in the nullspace of W.
In the special case that W has a one-dimensional nullspace, a more precise characterization of r 0 is possible. Let v 0 be in the nullspace of W with kv 0 k = 1. Note that W T also has a onedimensional nullspace (since W is a square matrix). Let v 2 be in the nullspace of W T with kv 2 k = 1. Since r 0 is in the nullspace of W, we can write r 0 = av 0 for some scalar, a. Now, dot product both sides of Eq (7) by v 2 to obtain where we have used that v 2 Á Wr = 0 since v 2 is in the nullspace of W T , which is orthogonal to Wr in the column space of W. Keeping only Oð1Þ terms and making the substitution r 0 = av 0 , we get yielding a concise expression for the amplified component of firing rates when W has a onedimensional nullspace. Note that v 0 is in the nullspace of W, so this result is consistent with the more general conclusions above.
Amplification in the nullspace under partial optogenetic stimulation. For the specific example of partial optogenetic stimulation considered in Fig 2, (8) gives Hence, ChR2-expressing neurons are amplified and non-expressing neurons are suppressed by optogenetic stimulation, as observed in simulations. A more precise result is given by expanding the full approximation from Eq (4) to obtain þ OðÞ: Here, OðsÞ is a constant proportional to s, c = |w IE /w II | and r off is the vector of firing rates in the balanced, ! 0, limit when stimulation is off (s = 0). Specifically, r off is the unique vector that satisfies Wr off + W XrX = 0, which is solvable even though W is singular because W X r X is in the column space of W, so balance can be maintained when s = 0.
The gs/ term in Eq (9) quantifies the amplification and suppression observed in simulations: Non-expressing neurons are suppressed by stimulation since −q < 0 and the response of ChR2-expressing neurons is amplified since 1 − q > 0 and s/ is large. Stimulating fewer excitatory neurons (q smaller), increases the rate of the stimulated neurons, as indicated by the 1 − q factor. The OðsÞ term shows why inhibitory neurons increase their rates by a smaller amount. In summary, the optogenetically induced suppression observed experimentally by Adesnik and Scanziani [3] is a generic feature of balanced or strongly coupled networks under partial stimulation.
Local imbalance with global balance explains intralaminar suppression and interlaminar facilitation. Interestingly, despite the break of balance at the level of ChR2-expressing and non-expressing subpopulations, global balance is maintained in this example. This can be understood by repeating the mean-field analysis above without partitioning neurons into ChR2-expressing and non-expressing sub-populations, thereby quantifying the global average of firing rate of all excitatory neurons. In particular, the average synaptic input, I = [I E I I ] T , to excitatory and inhibitory neurons is given by Eq (1) where W and W X are as in Eq (2), and to account for the fact that only a proportion q of the excitatory neurons receive the inward current from optogenetic stimulation. In this case, W is non-singular so the balanced solution in Eq (3) is applicable. Indeed, the average firing rates of all excitatory neurons in our spiking network simulation is close to the prediction from Eq (3) and even closer to the prediction from Eq (4) (Fig 2g; compare to Fig 1c). The average feedforward input to all excitatory neurons is canceled by net-inhibitory local input (Fig 2h; compare to Fig 1d). Hence, balance is maintained globally even though the network is imbalanced at the level of ChR2-expressing and non-expressing populations.
In the same study by Adesnik and Scanziani considered above [3], recordings were made in L5, which was not directly stimulated optogenetically, but receives synaptic input from L2/3. Interestingly, despite the fact that most excitatory neurons in L2/3 were suppressed during stimulation, firing rates in L5 increased.
To model these experiments, we interpreted the recurrent network from Fig 2 as a local neural population in L2/3, which sends synaptic projections to L5 (Fig 3a). We modeled a neural population in L5 identically to the L2/3 population, except its feedforward input came from excitatory neurons in the L2/3 network, instead of from Poisson-spiking neurons. As in experiments [3], L5 neurons increased their firing rates during stimulation (Fig 3b) and approximate balance was maintained (Fig 3c). This can be understood by noting that, in our model, L5 receives synaptic input sampled from all excitatory neurons in L2/3. Hence, the feedforward excitatory current to L5 neurons increases proportionally to the average excitatory firing rates in L2/3 during stimulation. As we showed above, this average rate increases (Fig 2e), despite the fact that most excitatory neurons in L2/3 are suppressed by stimulation. Hence, the combination of intralaminar suppression and interlaminar facilitation observed during optogenetic stimulation in experiments [3] results from the fact that the stimulated layer is locally imbalanced, but globally balanced during partial stimulation. These conclusions rely on an assumption that synaptic projections to L5 come from both ChR2-expressing and non-expressing excitatory neurons in L2/3. Imbalanced amplification of weak stimuli. Sufficiently small or large s would introduce negative rates in Eq (9), representing a regime in which non-expressing neurons cease spiking and the firing rate of ChR2-expressing neurons saturate at a high value. In this sense, firing rates do not truly have a Oð1=Þ component for very small. However, smaller values of allow weak stimuli (small s) to be strongly amplified. Strictly speaking, if one takes s $ OðÞ, then under the linear approximation in Eq (4), partial optogenetic stimulation would have an Oð1Þ effect on the average firing rate of stimulated and unstimulated subpopulations, but an OðÞ effect on globally averaged firing rates. In practical terms, this means that, in strongly coupled networks ( small), partial optogenetic stimuli can have a moderate effect on the firing rates of stimulated neurons while having a negligible effect on the average firing rates of all excitatory neurons.
To demonstrate this idea, we repeated the simulations from Fig 2 in a network with four times as many neurons (N = 2 × 10 4 ) where synaptic weights and probabilities were scaled so that $ 1= ffiffiffiffi N p (as in Fig 1f) and we reduced stimulus strength, s, as well. In this simulation, ChR2-expressing neurons' firing rates nearly doubled (Fig 4a) and non-expressing neurons were noticeably suppressed (Fig 4b). However, the change in the average firing rate of all excitatory neurons was nearly imperceptible (Fig 4c) and similarly for the firing rates of inhibitory neurons (Fig 4b and 4c). As a result, the firing rates in a downstream layer were unnoticeably modulated during stimulation (Fig 4d; compare to Fig 3). This effect could mask the effects of optogenetic stimulation in recordings.
Imbalanced amplification with nearly singular connectivity matrices. An apparent limitation of the results above is that they rely on the singularity of the connectivity matrix, W. Singularity is a fragile property of matrices that arises from structural symmetries. In the example above, singularity arises from our implicit assumption that local synaptic connectivity is independent of whether neurons express ChR2. Even a slight difference in connectivity to or from ChR2-expressing neurons would make W non-singular so that its nullspace would be , which receives external input from Poisson-spiking neurons (green X) and from partial optogenetic stimulation, sends excitatory synaptic input to an identical network (downstream layer). b) Average firing rate of ChR2-expressing excitatory neurons in the stimulated layer from simulations (light solid) and from the corrected approximation (Eq (4); dashed). c) Same as b, but for non-expressing excitatory (red) and inhibitory (blue) neurons in the stimulated layer. d) Same as b, but averaged over all excitatory neurons in the stimulated layer. e) Same as d, but for the downstream layer. Mean firing rates from simulations of the stimulated layer changed from 5.8 Hz before stimulation to 10.0 Hz during stimulation for ChR2-expressing neurons, from 5.9 to 5.1 Hz for non-expressing excitatory neurons, from 5.9 to 6.1 Hz averaged over all excitatory neurons, and from 7.8 to 8.0 Hz for inhibitory neurons. Mean firing rates from simulations of the downstream layer changed from 7.2 Hz to 8.1 Hz for excitatory neurons and from 8.5 Hz to 9.5 Hz for inhibitory neurons.
A matrix, W, is singular if it has λ = 0 as an eigenvalue. A matrix can therefore be considered approximately singular if it has an eigenvalue with small magnitude. Specifically, let λ be an eigenvalue of W with |λ| ( 1. Note that λ is also an eigenvalue of W T . Now let v be the associated eigenvector so that W T v = λv and assume that kvk = 1 without loss of generality. Take the projection of each term in Eq (7) onto the subspace spanned by v to get If proj v X $ Oð1Þ and proj v [Dr] * proj v r then this implies ðjlj þ Þ proj v r $ proj v X: where δ = |λ| + . This generalizes Eq (6) to the case where W is only approximately singular. In summary, the mechanism of imbalanced amplification is a general property of strongly coupled networks with singular or nearly singular connection matrices. We next show that networks with connection probabilities that depend on continuous quantities like distance or tuning preference necessarily have singular or nearly singular connectivity kernels and are therefore naturally susceptible to the amplification and suppression mechanisms described above.

Imbalanced amplification and suppression in continuously indexed networks
So far we considered networks with discrete subpopulations. Connectivity in many cortical circuits depends on continuous quantities like distance in physical or tuning space. To understand how the amplification and suppression mechanisms discussed above extend to such connectivity structures, we next considered a model of a visual cortical circuit. We arranged 2 × 10 5 AdEx model neurons (80% excitatory and 20% inhibitory) on a square domain, modeling a patch of L2/3 in mouse primary visual cortex (V1). Neurons received external input from a similarly arranged layer of 1.6 × 10 5 Poisson-spiking neurons, modeling a parallel patch of L4 (Fig 5a). We additionally assigned a random orientation preference to each neuron, modeling the "salt-and-pepper" distribution of orientation preferences in mouse V1. Connectivity was probabilistic and, as in cortex [40][41][42], inter-and intralaminar connections were more numerous between nearby and similarly tuned neurons. Specifically, connection probability decayed like a Gaussian as a function of distance in physical and orientation space (Fig 5b), where distance in both spaces was measured using periodic boundaries, i.e. wrapped Gaussians were used in place of regular Gaussians.
Amplification and suppression in simulations with spatially narrow stimuli. An oriented stimulus localized in the animal's visual field (Fig 5c) was modeled by imposing firing rate profiles in L4 that were peaked at the associated location in physical and tuning space, again with a Gaussian profile (Fig 5d and 5e). This produced external input to L2/3 that was similarly peaked, but was nearly perfectly canceled by net-inhibitory lateral input (Fig 5f and  5g). Excitatory and inhibitory firing rate profiles in L2/3 were also peaked at the associated location in physical and tuning space (Fig 5h and 5i), demonstrating that neurons in L2/3 were appropriately tuned to the stimulus.
A smaller visual stimulus was modeled by shrinking the spatial profile of firing rates in L4 while leaving the orientation-dependence of L4 rates unchanged (Fig 5j and 5k). As above, synaptic inputs and firing rate profiles were appropriately peaked in physical and orientation tuning space (Fig 5l-5o). However, the smaller stimulus produced a surprising change to firing rates in L2/3. Despite the fact that L2/3 neurons at all locations received less excitation from L4 (Fig 5l), peak firing rates in L2/3 increased and a surround suppression dynamic emerged (Fig  5n). Hence, a more localized external input produced an amplification and suppression dynamic similar to the one observed in our model of optogenetic stimulation (compare to Fig  2). On the other hand, responses in orientation tuning space were mostly unchanged by the smaller stimulus size (Fig 5m and 5o).
Mean-field theory of balance in two-dimensional spatial networks with orientation-tuning-specific connectivity. The mean-field theory of balanced networks was previously extended to continuously indexed networks in one and two dimensions [8,12,43]. We now review a straightforward extension to two spatial dimensions and one orientation dimension.  β a . b) Neurons are assigned random orientations and connection probability also depends on the difference, dθ, between neurons' preferred orientation. c) An oriented stimulus in the animal's visual field. d,e) The location of the stimulus is modeled by firing rates in L4 that are peaked at the location of the stimulus in physical and orientation space. f,g) Synaptic current to neurons in population E from the external network (green), the local network (purple) and total (black) as a function distance from the receptive field center and as a function of neurons' preferred orientation. h,i) Firing rate profiles of excitatory (red) and inhibitory (blue) neurons in the local network from simulations (light curves), classical balanced network theory (solid, dark curves; from Eq (13)) and under the linear correction (dashed; from Eq (17)) in physical and orientation space. j-o) Same as (d-i) except for a smaller visual stimulus, modeled by a narrower spatial firing rate profile in L4. Firing rates from Eq (13) are not shown in panels n and o because balance cannot be realized and Eq (13) cannot be applied when external input is narrower than recurrent connectivity (see main text, S1 Text, and [8] Eq (1) generalizes naturally to where r(x, θ) = [r E (x, θ) r I (x, θ)] T is the vector of mean firing rates of excitatory and inhibitory L2/3 neurons near spatial coordinates x = (x, y) with preferred orientation near θ, and similarly for the neurons' external input, X(x, θ), and total input, I(x, θ). The external input is given by θ) is the profile of firing rates in L4 The connectivity kernels, W and W X , are convolution integral operators defined by Here, w ab (x, θ) is the mean-field connection strength between neurons separated by x in physical space and θ in orientation space (see Methods), and [w ab Ã r b ](x, θ) denotes circular convolution with respect to x and θ, i.e., convolution with periodic boundaries. These convolution operators implement low-pass filters in orientation and physical space, capturing the effects of synaptic divergence and tuning-specific connection probabilities. Similar filters describe feedforward connectivity in artificial convolutional neural networks used for image recognition [44].
Taking ! 0 in Eq (10) shows that that firing rates must satisfy This is an analogue to Eq (7) for spatial networks. From here, one may be tempted to invert the integral operator W to obtain a spatial analogue of Eq (3). However, integral operators are never invertible [45]. Specifically, since Eq (11) is an integral equation of the first kind, there necessarily exist external input profiles, X(x, θ), for which Eq (11) does not admit a solution so that the classical balanced state cannot be realized [8]. This implies that there always exist inputs that prevent a continuously indexed network from maintaining excitatory-inhibitory balance. To better understand why this is the case, we follow previous work [8,12,43,46,47] in transitioning to the spatial Fourier domain to rewrite Eq (11) as Here, e rðn; kÞ ¼ ½e r E ðn; kÞ e r I ðn; kÞ T is a Fourier coefficient of r(x, θ) and similarly for e Xðn; kÞ ¼ f W X ðn; kÞe r X ðn; kÞ where n = (n 1 , n 2 ) is the two-dimensional spatial Fourier mode and k is the Fourier mode in tuning space. Importantly, the convolution operators above become ordinary matrices in the Fourier domain. Specifically, where e w ab ðn; kÞ is a Fourier coefficient of w a (x, θ). Note that going from Eq (11) to Eq (12) requires that W is a convolution operator and that the boundaries of the network are treated periodically, i.e., the convolutions are circular.
Solving Eq (12) gives an analogue to Eq (3) for spatial networks in the Fourier domain, This equation gives all Fourier coefficients, e rðn; kÞ. However, this solution is only viable when the inverse transform exists, i.e., when the Fourier series of e rðn; kÞ converges, which in turn requires that k e Xðn; kÞ k converges to zero faster than k f W ðn; kÞ k as n ! 0 and k ! 0. More specifically, e rðn; kÞ in Eq (13) must be square-summable. Hence, balance can only be realized when recurrent connectivity, quantified by f W ðn; kÞ, has more power at high spatial frequencies than external input, e Xðn; kÞ. In other words, for balance to be realized, external input, X (x, θ), cannot have "sharper" spatial features than the recurrent connectivity kernels, w ab (x, θ) for a, b = E, I.
The use of Fourier series to solve Eq (11) relies on the translation invariance of the integral operator, W, but the existence of external input profiles for which Eq (11) does not admit a solution is a more general feature that is true for any integral operators with Hilbert-Schmidt kernels. Specifically, the convolutions, w ab Ã r b , that define W and W X can be replaced by integrals of the form R w ab (u, v)r b (v)dv and the result that there exist X(u) for which Eq (11) does not admit a solution is still true whenever the kernels, w ab (u, v) are finite (no delta functions) and square-integrable [45].
Balance and imbalance in networks with Gaussian-shaped connectivity kernels. A more intuitive understanding of when and why balance is broken is provided by considering the Gaussian-shaped connectivity profiles used in our spiking network simulations. This explanation applies equally to the spatial profile of firing rates and connectivity in physical and orientation space, so we do not distinguish between the two in this discussion. Similar calculations were performed previously for spatial networks [8], so we only review the results here and discuss some of their implications here. Details of the calculations are provided in S1 Text.
When L4 firing rates and all connectivity kernels are Gaussian-shaped, all firing rates and synaptic input profiles are also Gaussian-shaped in the balanced state (see [8] and S1 Text), so the analysis can be performed solely in terms of the widths of each of these profiles.
Let σ a be the width of the firing rate profile in population a, α a the width of outgoing synaptic connections from the presynaptic neurons in population a, and β a the width of the spatial profile of synaptic input from population a = X, E, I (Fig 5a, 5d, 5f and 5h). Since L4 is the feedforward population in this context, we use X subscripts for L4 and E, I for L2/3 neurons.
Synaptic divergence broadens the profile of synaptic currents so that for a = E, I, X. This is due to the Gaussian shape of firing rates and connectivity, as well as the convolution that describes the resulting mean-field synaptic inputs (the convolution of a Gaussian with a Gaussian is Gaussian). For balance to be maintained, feedforward synaptic input from L4 must be precisely canceled by lateral synaptic input in L2/3. This, in turn, Combined with Eq (14), this implies that balance requires the widths of firing rate profiles in L2/3 to satisfy [8] This approximation accurately predicted firing rate profiles in our first spiking network simulation (Fig 5h and 5i, solid, dark curves have widths given by Eq (15)). Hence, by Eq (15), the requirement of cancellation in balanced networks implies that recurrent connectivity sharpens neurons' tuning, both in physical and orientation space. Interestingly, Eq (15) implies that the amount by which excitatory and inhibitory firing rate profiles are sharpened in balanced networks is determined by the width of their outgoing synaptic projections. Pyramidal neurons in L2/3 of mouse V1 preferentially target similarly tuned neurons in L2/3, but the tuning of these lateral connection probabilities is much broader than the tuning of pyramidal neurons' firing rates [41] (α E > σ E in orientation space). This observation is consistent with Eq (15): Excitatory neuron tuning curves are sharpened precisely because their outgoing connections are broadly tuned. Hence, sharpening of excitatory neuron tuning curves in L2/3 is naturally achieved in balanced networks with lateral excitation, without requiring lateral inhibition. Following the same line of reasoning, the broader orientation tuning of inhibitory neurons [40] (σ I larger) suggests that they project more locally in orientation tuning space than pyramidal neurons (α I < α E in orientation space).
Eq (15) also clarify when and why balanced network theory fails for continuously indexed networks. If external inputs are sharper than lateral connectivity (β X < α E or β X < α I ) in physical or orientation space, then Eq (15) do not yield solutions for σ E or σ I . In other words, balance requires that a E < b X and a I < b X because Eq (11) does not admit a solution when these inequalities are broken (see [8] and S1 Text). Intuitively, when the inequalities above are broken, recurrent connectivity paints with too broad a brush to cancel the sharper feedforward input, so balance cannot be realized. As a result, balanced network theory cannot be applied to the example in Fig 5j-5o with a smaller visual stimulus.
A linear correction to balance quantifies amplification and suppression in continuously indexed networks. We next derive a linear correction to Eq (13) that accounts for imperfect cancellation and, in doing so, gives firing rate approximations where classical balanced network theory fails. Specifically, we generalize the derivation of Eq (4) to continuously indexed networks. Under this linear approximation, firing rate profiles are given by solving This is an integral equation of the second kind, which generically admits firing rate solutions, r, even when Eq (11) does not [45]. We again transition to the Fourier domain so Eq (16) becomes From Eq (17), firing rates, r(x, θ), can be computed numerically through an inverse transform (the Fourier series over n and k), yielding an accurate approximation to firing rates from spiking network simulations even where classical balanced network theory fails (Fig 5n and 5o). Interestingly, the linear correction is slightly less accurate at predicting firing rates than the classical theory in Fig 5i. However, the main purpose of the corrected theory is not increased accuracy, but that it is applicable in situations where the classical theory is not. The amplification and suppression caused by the smaller visual stimulus can be roughly explained by the balanced amplification mechanism discussed previously. Since W is a lowpass filter, it approximately cancels high frequency components of firing rate profiles. Hence, high frequency components are in the approximate nullspace of the local connectivity operator, W, and are therefore amplified by the network through the same mechanism discussed for discrete networks previously.
A more precise explanation is given by first averaging firing rates over orientation preference by setting k = 0 in Eq (17) to give e rðnÞ that depends only on spatial frequency, and similarly for e XðnÞ and f W ðnÞ. The convolution operator, W, implements a low-pass filter, so f W ðnÞ is Oð1Þ in magnitude at low spatial frequencies and converges to zero at higher frequencies (large knk). The regularized inverse, ½D À f W ðnÞ À 1 , is therefore Oð1Þ in magnitude at low frequencies and Oð1=Þ at higher frequencies (Fig 6a, purple). When external input, X(x), has sharp features, e XðnÞ has power at higher spatial frequencies (Fig 6a, green), which are amplified by the Oð1=Þ component of ½D À f W ðnÞ À 1 while low frequencies remain Oð1Þ. The result is that the magnitude of e rðnÞ has a Oð1=Þ peak at a non-zero spatial frequency (Fig 6a,  black), introducing a high-amplitude, non-monotonic rate profile (as in Fig 5n; see [47] for a similar analysis). When X(x) has spatially broad features, e XðnÞ has little power at high spatial frequencies so that this amplification dynamic is weak or absent (as in Fig 5h). An identical argument applies in orientation space. In summary, high-frequency components of external input profiles are transmitted more strongly than low-frequency components in strongly coupled networks, and the cutoff frequency is determined by the width (α E or α I ) of lateral synaptic projections.
It is worth noting that the average firing rates (over all orientations and spatial positions) are given by the zero Fourier coefficient e rð0; 0Þ. When balance is broken by sharp external input features, the zero Fourier mode is not affected as long as mean firing rates, r(x, θ), input ( e XðnÞ, green) and the resulting firing rate profile (e rðnÞ, black) as a function of the spatial frequency, k n k¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 2 þ n 2 p , from the simulation in Fig 5j-5o. Magnitude is measured by the Frobenius norm for ½D À f W ðnÞ À 1 . Curves normalized by their peaks. b) Firing rates of excitatory (red) and inhibitory (blue) neurons with receptive fields at the center of a grating stimulus plotted as the width of the stimulus increases (represented by increasing σ X ) using parameters from Fig 5j-5o. c) Same as b, but the excitatory rate is plotted for different widths of the excitatory synaptic projection width, α E . d) Same as c, but firing rates in L4 are shaped like a disc with radius σ X instead of a Guassian with width parameter σ X .
https://doi.org/10.1371/journal.pcbi.1006048.g006 remain non-zero at all locations and orientations. Hence, sharp input features can break balance locally without breaking global, network-averaged, balance. This is analogous to the global balance obtained in the optogenetic example when local balance was broken at the level of subpopulations (Fig 2).
Implications of imbalanced amplification on receptive field tuning. We next considered a study by Adesnik et al. [48] in which drifting grating stimuli of varying sizes were presented to mice while recording from neurons in L2/3 of V1. In that study, pyramidal neurons' firing rates first increased then decreased as the stimulus size was increased. On the other hand, somatostatin-expressing (SOM) neuron's firing rates increased monotonically with stimulus size. Intracellular recordings combined with optogenetic stimulation in that study showed that SOM neurons project locally and pyramidal neurons form longer range projections.
To test our model against these findings, we applied Eq (17) to a network with local inhibition and longer-range excitation (α E > α I ) with increasing size of a visual stimulus (increasing σ X ). Our results are consistent with recordings in Adesnik et al., 2012 [48]: Excitatory neuron firing rates changed non-monotonically with stimulus size, while inhibitory neuron firing rates monotonically increased (Fig 6b). The non-monotonic dependence of excitatory firing rates on stimulus size in Fig 6b is explained by the mechanism of imbalanced amplification. When σ X is sufficiently small, balance is broken so imbalanced amplification introduces a large peak firing rate surrounded by suppression (as in Fig 5n). However, the total amount of external excitation introduced by the stimulus is proportional to the size of the stimulus, so a very small σ X introduces very little excitation and peak firing rates are small. As σ X increases, more excitation is recruited and the network is still imbalanced, which leads to increasingly large peak firing rates (as in Fig 5n). Once σ X becomes large, balance begins to be restored and the peak excitatory firing rate decreases to moderate values (as in Fig 5h).
The degree to which excitatory neurons suppress depends on the spatial width, α E , of lateral excitatory projections (Fig 6c) and suppression of inhibitory neurons similarly depends on the spatial width, α I , of lateral inhibition (not pictured). Specifically, suppression occurs when lateral connectivity is broader than feedforward input (α E > β X or α I > β X ) because this is when the balanced solution in Eq (15) disappears. When a sub-population's lateral connectivity is more localized than feedforward connectivity from L4 (α E < α X as in the lightest gray curve in Fig 6c; or α I < α X ), that sub-population cannot exhibit suppression since feedforward input width (b 2 X ¼ a 2 X þ s 2 X ) is always larger than lateral connectivity, regardless of the stimulus size (σ X ).
A similar line of reasoning explains why peak inhibitory neuron firing rates increase monotonically with stimulus size in Fig 6b. Inhibitory neurons in that example project locally (α I = α X ), so the inequality α I < β X is always satisfied because b X ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi a 2 X þ s 2 X p > a I whenever α I α X . Whenever α I < β X , inhibitory firing rates reflect their balanced state values which increase monotonically with the increase in total excitation induced by a larger stimulus.
Unlike SOM neurons, parvalbumin-expressing (PV) neurons were found to exhibit suppression by Adesnik et al., 2012 [48]. Hence, our theory predicts that PV neurons project more broadly in space than SOM neurons. Indeed, PV interneurons in L2/3 are primarily basket cells whose axons project to larger lateral distances than other inhibitory neuron subtypes such as Martinotti cells that comprise most SOM neurons [49].
We observed a unimodal dependence of firing rate on stimulus size (Fig 6c, all curves have a single peak). However, Rubin, et al. [47] observed a multi-modal, oscillatory dependence of firing rate on stimulus size in recordings and in a computational model. In that study, the drifting grating stimuli were disc-shaped with a sharp cutoff of contrast at the edges of the disc. Above, we considered a Gaussian-shaped contrast profile with soft edges (Fig 6c, inset). Repeating our calculations with a sharp-edged, disc-shaped stimulus (Fig 6d, inset) produced an oscillatory dependence of firing rate on stimulus size (Fig 6d), as observed by Rubin et al.. This oscillation only arose when lateral synaptic projections were narrower than the stimulus size (α E small). The oscillation results from a Gibbs phenomenon: The sharp edge in the stimulus produces high-frequency power in e X, which passes through the high-pass filter ½D À f W À 1 when α E is small.
We next explored the functional consequences of these results on receptive field tuning. We first considered a disc-shaped grating stimulus (Fig 7a), producing a disc-shaped firing rate profile in L4 (Fig 7b). Synaptic divergence causes the profile of synaptic input from L4 to L2/3 to be "blurred" at the edges (Fig 7c), as quantified by the low-pass filter, W X . This illustrates a fundamental problem in receptive field tuning: Synaptic divergence from one layer to another implements a low-pass filter that blurs sharp features. This problem is resolved by our observation above that lateral, recurrent connectivity implements a high-pass filter. If the width of lateral, excitatory connections in L2/3 is similar to that of feedforward connections from L4, the high-pass filter implemented by the recurrent network cancels the low-pass filter implemented by feedforward connectivity, effectively implementing a deconvolution that can recover the sharpness of firing rate profiles in L4 (Fig 7d-7f). Hence, counterintuitively, broader lateral excitation actually sharpens receptive field tuning. Broadening lateral connections further increases the sharpness of the firing rate profiles, but introduce oscillatory, Gibbs phenomena near sharp features (Fig 7f). These points are illustrated more clearly in an example with an asymmetrically shaped stimulus (Fig 7g-7l). Hence, the high-pass filter described above corrects the blurring caused by synaptic divergence between layers in V1.
In summary, imbalanced amplification and linear rate models provide a concise and parsimonious theoretical basis for understanding how suppression, amplification and tuning depends on the profile of neuron's incoming and outgoing synaptic projections in physical and orientation tuning space. Imbalanced amplification and suppression reverse the blurring introduced by interlaminar synaptic divergence. a) A disc-shaped grating stimulus gives rise to b) a disc-shaped firing rate profile, r X (x), in L4 with slightly blurred edges (achieved by convolving contrast from a with a Gaussian kernel). c) Input, X E (x), from L4 to excitatory neurons in L2/3 is blurred by synaptic divergence, which effectively applies a low-pass filter, W X , to the L4 rates. d) Excitatory firing rates in L2/3 are sharper than external input when lateral excitation is similar, but smaller, in width than interlaminar excitation (α E = 0.85α X ). e) Same as c, but lateral excitation is exactly as broad as interlaminar excitation (α E = α X ), which sharpens the edges further, making firing rates in L2/3 similar to those in L4. f) Same as c, but lateral excitation is broader than interlaminar excitation (α E = 1.1α X ), which sharpens the edges even further, but also introduced suppressed regions due to Gibbs phenomena. g-l) Same as a-f, but contrast was determined by the brightness of a photograph. Horizontal and vertical axes are neurons' receptive fields. https://doi.org/10.1371/journal.pcbi.1006048.g007

Discussion
We described a theory of amplification in cortical circuits arising from a local imbalance that occurs when recurrent connectivity structure cannot cancel feedforward input. We showed that this imbalanced amplification is evoked by optogenetic stimuli in somatosensory cortex and sensory stimuli in visual cortex, since these stimuli cannot be canceled by the connectivity structure in those areas. Our theoretical analysis of imbalanced amplification explains several observations from cortical recordings in those areas.
Even though firing rates in balanced networks in the large N limit do not depend on neurons' f-I curves (see Eq (3)), quantifying firing rates under imbalanced amplification relies on a finite size correction that requires an assumption on how firing rates depend on neurons' input. For simplicity, we used an approximation that assumes populations' mean firing rates depend linearly on their average input currents, giving rise to Eqs (4) and (17). In reality, neurons' firing rates depend nonlinearly on their mean input currents, and also depend on higher moments of their input currents. However, the salient effects of imbalanced amplification are not sensitive to our assumption of linearity. For instance, Eq (5), which quantifies the strong synaptic currents evoked under imbalanced amplification, does not depend on any assumption about neurons' f-I curves. However, the precise value of the firing rates elicited by this strong input does depend on neurons' f-I curves. We found that the linear approximation to f-I curves in Eqs (4) and (17) performed well at approximating firing rates in our spiking network simulations and also explained several observations from cortical recordings. This may be partly explained due to the fact that our spiking network simulations used neuron models that exhibit spike frequency adaptation, which is known to linearize f-I curves [50,51] and help networks maintain balance [10]. However, the linear approximation we used cannot explain some phenomena that rely on thresholding and other nonlinear transfer properties [47,52]. The notion of imbalanced amplification extends naturally to models with nonlinear transfer functions and future work will consider the implications of nonlinearities.
Balanced networks are related to, but distinct from, inhibitory stabilized networks (ISNs) [2,47,53] and stabilized supralinear networks that can transition between ISN and non-ISN regimes [47]. The primary distinction is that ISNs are defined by moderately strong recurrent excitation (strong E ! E) whereas balanced networks are defined by very strong external, feedforward excitation (strong X ! E) canceled by similarly strong net-inhibitory recurrent connectivity. Classical balanced networks are necessarily inhibitory stabilized at sufficiently large N (small ) unless w EE = 0. However, strongly coupled (approximately balanced) networks can be non-ISN at moderately large N (small ) if w EE is small. Cat V1 is believed to be inhibitory stabilized, which can be used to explain its surround suppression dynamic [2]. However, evidence from optogenetic and electrophysiological studies, suggests that mouse L2/3 V1 might not be inhibitory stabilized: Lateral connection probability is small between pyramidal neurons (small w EE ) [49], stimulation of PV neurons does not produce the paradoxical effects that characterize ISNs [54], and modulating pyramidal neuron firing rates only weakly modulates excitatory synaptic currents in local pyramidal neurons [48,54]. Nonetheless, pyramidal neurons and PV neurons in mouse V1 exhibit surround suppression [48], which we showed is explained by imbalanced amplification.
Despite the similarity in their names, the mechanism of imbalanced amplification studied here is fundamentally different from the mechanism of balanced amplification [55]. First, imbalanced amplification is related to steady-state firing rates, while balanced amplification is a dynamical phenomenon. Moreover, balanced amplification is intrinsic to the local, recurrent circuit: It produces large firing rate transients when local, recurrent inhibition is inefficient at canceling local, recurrent excitation. Imbalanced amplification, on the other hand, produces large steady state firing rates when local, recurrent input is unable to effectively cancel feedforward, external excitation.
The analysis of our spatially extended network model relied on an assumption of periodic boundaries in space, which are not biologically realistic, but approximate networks with more realistic boundary conditions [8]. Without periodic boundary conditions, the integral eqs (10), (11) and (16) are equally valid, but the integrals are defined by regular convolutions in space instead of circular convolutions. As a result, the spatial Fourier modes do not de-couple, so Eqs (12), (13) and (17) are no longer valid, though they should still offer a good approximation when connectivity is much narrower than the the spatial domain [8]. In addition, anisotropic connectivity statistics, arising for example from tuning dependent connectivity in visual cortical circuits with coherent orientation maps [56], would prevent the integral operator in Eqs (10), (11) and (16) from being a convolution operator, and therefore preclude the use of Fourier series for the solution. Future work will consider the effects of non-periodic boundaries and non-convolutional connectivity kernels on spatially extended balanced networks.
We only considered neuron models with current-based synaptic input. In S2 Text, we show that our numerical results and analysis extend naturally to models with more realistic conductance-based synapses. The analysis makes use of an approximation that relates a conductancebased model to a current-based model with similar membrane potential statistics [57][58][59].
We focused on firing rates, but sensory coding also depends on variability and correlations in neurons' spike trains. Our previous work derived the structure of correlated variability in heterogeneous and spatially extended balanced networks when connectivity structure prevents positive and negative correlations from cancelling, effectively providing an analogous theory of imbalanced amplification of correlated variability [12]. Combining those findings with the theory of steady-state firing rates presented here could yield a more complete theory of neural coding in cortical circuits and the effects of imbalanced amplification on coding.

Materials and methods
We modeled recurrently connected networks with N neurons, composed of N E = 0.8N excitatory and N I = 0.2N inhibitory neurons. The recurrent network receives external input from a network of N X neurons that drive the recurrent network. The membrane potential of neuron j from the excitatory (a = E) or inhibitory (a = I) population has Adaptive Exponential integrate-and-fire dynamics, C m dV a j dt ¼ À g L ðV À E L Þ þ g L D T exp½ðV À V T Þ=D T þ I a j ðtÞ À w t w dw dt ¼ À w: Whenever V a j ðtÞ > V th , a spike is recorded, the membrane potential is held for a refractory period τ ref then reset to a fixed value V re , and w is incremented by B. Neuron model parameters for all simulations were τ m = C m /g L = 15ms, E L = −72mV, V T = −60mV, V th = −15mV, Δ T = 1.5mV, V re = −72mV, τ ref = 1ms, τ w = 150ms and B = C m 0.267 mV/ms. We write all currents in terms of C m , so C m can be any constant. Membrane potentials were also bounded below by V lb = −100mV. Synaptic input currents were defined by I a j ðtÞ ¼ ½X a j ðtÞ þ R a j ðtÞC m where X a j ðtÞ is the feedforward input and R a j ðtÞ the recurrent input to neuron j in population a = E, I. The recurrent input was defined by R a j ðtÞ ¼ where t b;k n is the nth spike time of neuron k in population b = E, I. The external input to the recurrent network is defined similarly by X a j ðtÞ ¼ where t X;k n is the nth spike time of neuron k = 1, . . ., N X in population X. Each coefficient, J ab jk , represents the synaptic weight from presynaptic neuron k in population b to postsynaptic neuron j in population a. For all simulations, we modeled synaptic kinetics using η b (t) = exp(−t/τ b )/τ b for t > 0 where τ E = 8ms, τ I = 4ms, and τ X = 10ms. Note that the integral of η b (t) over time is equal to 1 for all three kernels, so the choice of time constant, τ b , does not effect time-averaged synaptic currents. We used τ I < τ E < τ X to prevent excessive synchronous events that break the balanced state. While inhibition may be faster than excitation in many cortical circuits, excitatory neurons are more likely to contact distal dendrites and inhibitory neurons are more likely to contact the soma [60,61], which could make inhibition functionally faster than excitation. In any case, using fast inhibition is common practice in spiking network simulations with strong or dense connectivity [11,12,47,62,63] and a complete resolution of this issue is outside the scope of this study.
In Figs 1, 2 and 3 an extra term, S = 2 mV/ms, was added to X E j ðtÞ for stimulated neurons during the second half of the simulation to model optogenetic stimulation. We used N E = 4000, N I = 1000 and N X = 4000 (so N = 5000) except for Fig 1f where all N b values were scaled. Connections were drawn randomly with connection probabilities p EE = p IE = p IX = 0.1, p EI = p II = p EX = 0.2. Specifically, for each neuron in presynaptic population b = E, I, X, we sampled p ab N a postsynaptic targets from population a = E, I randomly and uniformly with replacement. Since outgoing connections were sampled with replacement, some neurons connected multiple times to other neurons. Synaptic weights were then defined by  For Figs 1f and 4, the values of J ab and the values of p ab were each multiplied by (5000/N) 1/4 so that they were unchanged at N = 5000 and so that $ 1= ffiffiffiffi N p . This is slightly different from the more common practice of fixing small connection probabilities and scaling J ab like 1= ffiffiffiffi N p . We instead fixed a relatively dense connectivity at N = 5000 and the network became increasingly sparse and weakly connected at increased N. Both approaches have the same mean-field (since the mean-field only depends on the product of p ab and J ab ), but our approach prevents excessively small synaptic weights at large N and prevents dense connectivity at large N, which is computationally expensive and susceptible to oscillatory and synchronous spiking.
Spike times in the external population were modeled as independent Poisson processes with r X = 5 Hz. In Fig 3,  , and connections probabilities were orientation space. The Fourier series of the convolution kernels defined above turns convolution into multiplication in the Fourier domain, from which Eq (10) gives e I ¼ ð1=Þ½ f W e r þ e X where e X, f W , and f W X are defined in Results with e w ab ðn; kÞ ¼ w ab exp ½À 2p 2 ðjnj 2 a 2 b þ k 2 a 2 b;y Þ, w ab ¼ e w ab ð0; 0Þ ¼ J ab p ab N b =ðJ EX p EX N X Þ, and k n k 2 ¼ n 2 1 þ n 2 2 . Using the linear approximation, r = gI then gives Eq (17). Firing rates for dashed curves in Fig 5 and all firing rates in Figs 6 and 7 were obtained by first computing Eq (17), then inverting the Fourier transform numerically using an inverse fast Fourier transform. Solid curves in Fig 5 were computed similarly, except using Eq (13) in place of Eq (17).
All simulations and numerical computations were performed on a MacBook Pro running OS X 10.9.5 with a 2.3 GHz Intel Core i7 processor. All simulations were written in a combination of C and Matlab (Matlab R 2015b, MathWorks). The differential equations defining the neuron model were solved using a forward Euler method with time step 0.1 ms.
Supporting information S1 Text. Derivation of firing rates in balanced networks with Gaussian-shaped connectivity kernels. This supplementary text contains details of the analysis of balanced networks with connection probabilities that decay like a Gaussian with distance in physical and orientation space. (PDF) S2 Text. Balance and imbalanced amplification in a model with conductance-based synaptic and optogenetic inputs. This supplementary text contains simulations and analysis of balance and imbalanced amplification in a network model with conductance-based synapses. (PDF)