A Balanced Memory Network

A fundamental problem in neuroscience is understanding how working memory—the ability to store information at intermediate timescales, like tens of seconds—is implemented in realistic neuronal networks. The most likely candidate mechanism is the attractor network, and a great deal of effort has gone toward investigating it theoretically. Yet, despite almost a quarter century of intense work, attractor networks are not fully understood. In particular, there are still two unanswered questions. First, how is it that attractor networks exhibit irregular firing, as is observed experimentally during working memory tasks? And second, how many memories can be stored under biologically realistic conditions? Here we answer both questions by studying an attractor neural network in which inhibition and excitation balance each other. Using mean-field analysis, we derive a three-variable description of attractor networks. From this description it follows that irregular firing can exist only if the number of neurons involved in a memory is large. The same mean-field analysis also shows that the number of memories that can be stored in a network scales with the number of excitatory connections, a result that has been suggested for simple models but never shown for realistic ones. Both of these predictions are verified using simulations with large networks of spiking neurons.


Introduction
A critical component of any cognitive system is working memory-a mechanism for storing information about past events, and for accessing that information at later times. Without such a mechanism, even simple tasks, such as deciding whether to wear a heavy jacket or a light sweater after hearing the weather report, would be impossible. Although it is not known exactly how storage and retrieval of information is implemented in neural systems, a very natural way is through attractor networks. In such networks, transient events in the world trigger stable patterns of activity in the brain, so by looking at the pattern of activity at the current time, other areas in the brain can know something about what happened in the past.
There is now considerable experimental evidence for attractor networks in areas such as inferior temporal cortex [1][2][3], prefrontal cortex [4][5][6][7][8][9], and hippocampus [10,11]. And from a theoretical standpoint, it is well understood how attractor networks could be implemented in neuronal networks, at least in principle. Essentially, all that is needed is an increase in the connection strength among subpopulations of neurons. If the increase is sufficiently large, then each subpopulation can be active without input, and thus ''remember'' events that happened in the past.
While the basic theory of attractor networks has been known for some time [12][13][14], moving past the ''in principle'' qualifier, and understanding how attractors could be implemented in realistic, spiking networks, has been difficult. This is because the original Hopfield model violated several important principles: neurons did not obey Dale's law; when a memory was activated, neurons fired near saturation, much higher than is observed experimentally in working memory tasks [1,15]; and there was no null background state-no state in which all neurons fired at low rates.
In spite of these advancements, there are still two fundamental open questions. One is: how can we understand the highly irregular firing that is observed experimentally in working memory tasks [24]? Answering this question is important because irregular firing is thought to play a critical role both in how fast computations are carried out [25] and in the ability of networks to perform statistical inference [26]. Answering it is hard, though, because, as pointed out in [27], with naive scaling the net synaptic drive to the foreground neurons (the neurons that fire at elevated rate during memory) is proportional to the number of connections per neuron. Consequently, because of the high connectivity observed in cortex, the mean synaptic drive is much larger than the fluctuations, which implies that the foreground neurons should fire regularly. Moreover, as pointed out by Renart et al. [28], even for models that move beyond the naive scaling and produce irregularly firing neurons, the foreground neurons still tend to fire more regularly than the background neurons, something that is inconsistent with experiments [24].
Several studies have attempted to get around this problem, either directly or indirectly [22,[27][28][29]. Most of them, however, did not investigate the scaling of the network parameters with its size (i.e., with the number of neurons and connections). So, although parameters were found which led to irregular activity, it was not clear how those parameters should scale as the size of the network increased to realistic values. In the two that did investigate scaling [27,28], irregular firing was possible only if a small fraction of neurons was involved in each memory; i.e., only if the coding level was very small. Although there have been no direct measurements of the coding level during persistent activity, at least to our knowledge, experiments in superior temporal sulcus [30] suggest that it is much larger than the one used in these models. We should point out, though, that the model of Renart et al. [28] is the only one in which the foreground neurons are at least as regular as the background neurons.
The second open question is: what is the storage capacity of realistic attractor networks? That is, how many different memories can be stored in a single network? Answering this is critical for understanding the highly flexible and seemingly unbounded memory capacity observed in animals. For simple, albeit unrealistic, models the answer is known: as shown in the seminal work of Amit, Gutfreund, and Sompolinsky [31], the number of memories that can be stored in a classical Hopfield network [12] is about 0.14 times the number of neurons. For slightly more realistic networks the answer is also known [16,19,21,27,[32][33][34][35][36][37][38]. However, even these more realistic studies lacked biological plausibility in at least one way: connectivity was all-all rather than sparse [19,21,33,38], the neurons were binary (either on or off, with nothing in between) [16,19,21,32,33,37], there was no null background [16,32,33,35,37,38], the firing rate in the foreground state was higher than is observed experimentally [16,27,32,33,36,37], or the coding level was very small [27,36].
Here we answer both questions: we show, for realistic networks of spiking neurons, how irregular firing can be achieved, and we compute the storage capacity. Our analysis uses relatively standard mean-field techniques, and requires only one assumption: neurons in the network fire asynchronously. Given this assumption, we first show that neurons fire irregularly only if the coding level is above some threshold, although a feature of our model is that the foreground neurons are slightly more regular than the background neurons. We then show that the maximum number of memories in our network-the capacity-is proportional to the number of connections per neuron, a result that is consistent with the simplified models discussed above. These predictions are verified with simulations of biologically plausible networks of spiking neurons.

Results Model
To address analytically the issues of irregularity and storage capacity in attractor networks, we consider a model in which neurons are described by their firing rates. Although firing rate models typically provide a fairly accurate description of network behaviour when the neurons are firing asynchronously [39,40], they do not capture all features of realistic networks. Therefore, we verify all of our predictions with large-scale simulations of spiking neurons.
Our network consists of two populations, one excitatory and one inhibitory, with N E neurons in the former and N I in the latter. (In general we use E for excitation and I for inhibition.) We represent the firing rate of the ith neuron in pool Q(¼E,I) by v Qi . As we show in the section ''Fast fluctuations,'' and discuss below, the time evolution equations for the firing rates are given by where s E and s I are the excitatory and inhibitory time constants, h Qi is the synaptic input to the ith neuron in pool Q, and F Q (h) is a function that tells us the steady state firing rate of a neuron receiving synaptic input h. This function, which has a relatively stereotyped quasi-sigmoidal shape, can be determined analytically (or semi-analytically) for specific noise models [41][42][43], and numerically for more realistic models [40]. The synaptic drive, h Qi , is related to the activity of the presynaptic neurons via whereJ QR ij is the synaptic weight from the jth neuron in pool R to the ith neuron in pool Q, andh Qex is the external, purely excitatory, input to neurons in pool Q. Finally, the steadystate firing rate of each neuron is determined by setting dv Ei / dt and dv Ii /dt to zero, yielding the equation The bulk of our analysis focuses on solving Equation 3; we use the dynamics, Equation 1, only when investigating stability. Our goal is to determine the conditions that support retrieval states-states such that subpopulations of neurons have elevated firing rates.
Since the gain functions, F Q (h), that we use in Equation 1 play such a central role in our analysis, we briefly justify them here; for additional details, see the section ''Fast fluctuations.'' These gain functions come from an average over the fast temporal fluctuations of the synaptic input-basically, filtered spikes. Calculating the temporal fluctuations selfconsistently is a hard problem [44], but, fortunately, it's not a

Author Summary
A critical component of cognition is memory-the ability to store information, and to readily retrieve it on cue. Existing models postulate that recalled items are represented by self-sustained activity; that is, they are represented by activity that can exist in the absence of input. These models, however, are incomplete, in the sense that they do not explain two salient experimentally observed features of persistent activity: low firing rates and high neuronal variability. Here we propose a model that can explain both. The model makes two predictions: changes in synaptic weights during learning should be much smaller than the background weights, and the fraction of neurons selective for a memory should be above some threshold. Experimental confirmation of these predictions would provide strong support for the model, and constitute an important step toward a complete theory of memory storage and retrieval.
problem we have to solve. As we show in the section ''Fast fluctuations,'' in the limit that each neuron receives a large number of connections, the temporal fluctuations experienced by all the excitatory neurons have the same statistics, as do the temporal fluctuations experienced by all the inhibitory neurons. Thus, we can use a single function, F E (h), for the excitatory neurons, and another function, F I (h), for the inhibitory ones. Of course, we won't be able to calculate the shape of F Q without knowing the structure of the temporal fluctuations. However, as we show below, the precise shapes of the gain functions don't play a strong role in our analysis.
Connectivity. The main determinant of network behaviour, at least in this model, is the set of connection strengths, thẽ J QR ij . To choose connection strengths that will lead to attractors, we build on the model proposed by Hopfield more than two decades ago [12]. In that model, random patterns are stored via a Hebbian learning rule, so connection strengths among neurons have the form where A ij is the strength of the connection from neuron j to neuron i, n l i ¼ 1 if neuron i participates in pattern l and n l i ¼ 0 otherwise,b is a constant that determines the memory strength, and p is the number of patterns. For each neuron, the probability of participating in a given pattern, l, is equal to the coding level, which we denote a. Thus, With this definition, the term ðn l j À aÞ in Equation 4 ensures that, on average, P j A ij is zero. Thus, the learning rule does not change the total synaptic weight onto a neuron, a form of postsynaptic normalisation that has been observed experimentally in cultured networks [45,46].
While Equation 4 produces a network that exhibits attractors, it is inconsistent with biology in at least two important ways. First, the neurons can exhibit both excitatory and inhibitory connections (for fixed presynaptic neuron j, A ij can be positive for some postsynaptic targets i and negative for others), which violates Dale's law. Second, connectivity is all to all, which is inconsistent with the sparse connectivity seen in cortex [47]. Both can be fixed by introducing sparse, random background connectivity among excitatory and inhibitory neurons, and adding a threshold so that neurons are either excitatory or inhibitory, but not both. This yields a set of connection strengths of the form where theJ QR set the background connection strengths (with, of course,J EE andJ IE positive andJ EI andJ II negative), and c QR ij tells us whether neuron j of type R is connected to neuron i of type Q. We assume that the connection probability is independent of type, so With this connectivity matrix, every neuron in the network projects to, on average, K E excitatory and K I inhibitory neurons, and every neuron receives, on average, K E excitatory and K I inhibitory connections, where The probability of connection, c, is assumed to be much smaller than 1, leading to a sparsely connected network [47], and it is independent of the size of the network unless otherwise stated. While we could have made the connectivity scheme more general by letting the connection probability between neurons depend on their type and/or by letting the nonzero c QR ij in Equation 7 have some variability, this would merely add complexity without changing any of our conclusions.
Although we are including the threshold-linear operator in Equation 6 (and also in the simulations), we neglect it in the forthcoming theoretical analysis. This is because A ij tends to be small: its mean is zero and, as we discuss in the sections ''Storage capacity'' and ''Mean-field equations,'' its variance is O(p/K E ). Thus, as long as p is sufficiently small compared with K, the threshold-linear operator can be neglected. For our model, we find that p/K is at most about 0.01, which means that the threshold-linear operator is unlikely to have much effect. Importantly, even if p/K were large, the scaling relation that we derive for storage capacity, i.e., p max } K, would still be correct; the only effect would be a slight modification to the precise value of p max /K [16].

Network Equilibria
As discussed above, much of our focus in this paper is on solving Equation 3. For even moderate size networks, this corresponds to solving thousands of coupled, highly nonlinear equations, and for large networks that can number into the millions. We do not, therefore, try to find a particular solution to this equation, but instead look for a statistical description-a description in terms of probability distributions over excitatory and inhibitory firing rates. The main tool we use is self-consistent signal-to-noise analysis [48,49]. The idea behind this analysis is to treat the synaptic input (h Ei and h Ii in Equation 3) as Gaussian random variables. Solving Equation 3 then reduces to finding, self-consistently, their means and variances.
Because h Ei and h Ii consist of 2K (very weakly) correlated terms, where naive central limit arguments tell us that the standard deviations of these quantities should be smaller than their means by a factor of K 1/2 . It would seem, then, that in the kinds of high connectivity networks found in the brain, where K is on the order of 5,000-10,000, neuron-to-neuron fluctuations in firing rate would be small, on the order of K À1/2 . By the same reasoning, temporal fluctuations in the firing rates would also be small, again on the order of K À1/2 . Neither of these, however, are observed in biological networks: there are large fluctuations in firing rate both across neurons and over time [24,[50][51][52][53].
To resolve this apparent contradiction, one need only notice that h Ei and h Ii consist of both positive and negative terms (the first and third terms in Equation 2 are positive; the second is negative). If these terms approximately cancel-to within O(K À1/2 )-then both the mean and standard deviation of the synaptic drive will be on the same order, and network irregularity will be restored. As showed by van Vreeswijk and Sompolinsky in a groundbreaking set of papers [25,54], under fairly mild conditions this cancellation occurs automatically, thus placing networks very naturally in what they called the balanced regime. In this regime, fluctuations across both neurons and time are large. Whether networks in the brain really operate in the balanced regime is not completely clear, although recent experimental evidence has come down strongly in favour of this hypothesis [55,56].
While the work of van Vreeswijk and Sompolinsky was extremely important in shaping our understanding of realistic recurrent networks, their focus was primarily on random connectivity. The situation, however, is more complicated in attractor networks. That's because these networks consist of three classes of neurons rather than two: background excitatory neurons and background inhibitory neurons, as found in randomly connected networks, but also foreground excitatory neurons. Our goal in the next several sections is to understand how all three classes can be balanced, and thus fire irregularly.
Strong synapses and the balanced condition. A reasonable constraint to place on our theoretical framework is that, in the large K limit, our results should be independent of K. This suggests that the synaptic strength, theJ QR in Equation 6, should scale as K À1/2 . With this scaling, the mean value of the positive and negative terms in h Ei and h Ii become O(K 1/2 ); with cancellation these terms are O(1), and the variance is also O(1). Thus, if the gain functions, the F Q (h) in Equation 3, are also O(1), our results will be independent of the number of connections. To make the K À1/2 scaling explicit, we define a new set of synaptic strengths and external input, which we denote J QR and h Qex , respectively, where J QR and h Qex are both O(1) and, recall, tells us how to scale the background connectivity, but it does not directly apply to the part of the connection matrix associated with memories, A ij . To determine how A ij should scale, we need only note that the mean contribution from the memories should be O(1)-sufficiently large to have an effect, but not so large as to overwhelm the background. Consequently, A ij should scale as 1/K (see the section ''Mean-field equations'' for details), which we can guarantee by defining a new variable, b, via the relationb where b is O(1) and the factor a(1Àa) is for convenience only.
Mean-field equations for the retrieval states. Now that we have the ''correct'' scaling-scaling that makes our results independent of network size and ensures that the mean and variance of the synaptic input are both O(1)-we can apply self-consistent signal-noise analysis to Equation 3. The first step is to divide the excitatory and inhibitory synaptic currents (h Ei and h Ii ) into two pieces: one that is nearly independent of index, i (the ''mean''), and one that is a random variable with respect to i (the fluctuating piece). To do that, we rewrite the synaptic current in terms of our new variables, J QR and b, rather thanJ QR ij andb. Combining Equations 4, 6, 9, and 10 with Equation 2, we have Note that Equation 11 is identical to Equation 2; it is just expressed in different variables. For the terms in brackets, the mean and fluctuating pieces are easy to compute: the mean comes from replacing c QR ij by its average, c, and the fluctuating piece comes from replacing c QR ij by the residual, c QR ij À c. For the second term in Equation 11a, separating the mean from the fluctuating piece is harder, as there is a nontrivial dependence on i associated with the p memories. Ultimately, however, we are interested in the case in which only one memory is retrieved, so when computing the mean we can consider only one term in this sum on l; the other p À 1 terms contribute only to the fluctuations. Assuming, without loss of generality, that the first memory is retrieved, averaging over the randomness associated with the sparse connectivity allows us to replace c EE ij with c, and we find that the mean of the last term in Equation 11a is proportional to n 1 i . Putting all this together, we arrive at the eminently reasonable result that the mean excitatory and inhibitory synaptic currents are linear in the mean excitatory and inhibitory firing rates, and the mean excitatory current has an extra, memory-induced, dependence proportional to n 1 i . Dropping the superscript ''1'' (a step taken only to simplify the equations), we find that the synaptic current may be written where h E and h I are the averages of the terms in brackets on the right-hand side of Equation 11, n i bm is the mean contribution from the first memory, and dĥ Ei and dh Ii contain everything else. More specifically, the terms in Equation 12 are as follows. First, h E and h I are given by To solve these equations, we use the fact that there are a large number of neurons; this allows us to turn the sum over i into an integral over the probability distributions of dĥ E and dh I , denoted pðdĥ E Þ and pðdh I Þ, respectively. Replacing the sum by an integral in Equation 16, and also averaging over n i , the mean-field equations become where the subscript on the angle brackets indicates an average with the statistics given in Equation 5. Because both dĥ E and dh I are Gaussian random variables (see the section ''Fast fluctuations''), these integrals are reasonably straightforward; what makes them at all difficult is that the variance of dĥ E and dh I must be found selfconsistently. This results in two more equations, for a total of five (see Equation 47). This is still far fewer than our original set of thousands or more. And the situation gets even better: it turns out that we really only need to consider three, at least if all we want to do is gain qualitative insight into how attractor networks function. That's because the integrals are simply Gaussian convolutions, so all they do is smooth the gain functions. Using a bar to denote the smoothed functions, and performing the average over n (which is straightforward because it has simple 0/1 statistics; see Equation 5), we have These equations-which are identical in form to the ones derived in [23]-are oversimplified versions of the full meanfield equations. Basically, the bar over F hides a dependence on two additional order parameters-the second moments of the excitatory and inhibitory firing rates-which in turn depend on our main order parameters, m E , m I , and m. While these dependencies are important for making detailed predictions, for an intuitive picture of what the mean-field equations mean they can be ignored. Consequently, in the next several sections, we focus on Equations 17a-17c, which we refer to as the reduced mean-field equations. At the end of the next section, we argue that, under very general conditions, all the conclusions we draw based on the reduced mean-field equations apply to the full set (which are given in Equation 47).
Reduced mean-field equations in the infinite K limit. When solving the reduced mean-field equations, we have a choice: we can think of them as functions of m E , m I , and m, or as functions of h E , h I , and m. Since m E and m I are related to h E and h I via an invertible transformation-Equation 13-the two prescriptions are identical. The latter, however, turns out to be more convenient, especially in the infinite K limit. To see why, we need only solve Equation 13 for the mean firing rates, which yields are the mean firing rates in the infinite K limit and is the determinant of the background connectivity matrix; as shown in [54] and the section ''Stability analysis,'' D must be positive for the background to be stable. Since we are in the balanced regime, h E and h I are O(1). Consequently, in the infinite K limit, the mean excitatory and inhibitory firing rates are simply given by m E0 and m I0 , respectively, independent of h E and h I . Using this fact, the reduced mean-field equations, Equation 17, become, in the K ! ' limit, An important feature of these equations is that h I decouples from h E and m. This greatly simplifies the analysis, since it means we can find the equilibrium value of h I simply by inverting F I : Our approach to finding the equilibrium values of h E and m is a graphical one: we plot, in h E À m space, the two curves that correspond to the solutions to Equations 21a and 21b-the h E and m nullclines, respectively-and look for their intersections. The goal is to determine the conditions under which there are multiple intersections, with at least one of them corresponding to an equilibrium with m . 0, and thus to a retrieval state.
To be as general as possible, we make only two assumptions: F E ðhÞ is monotonicaly increasing, and it is quasisigmoidal, where we use ''quasi-sigmoidal'' to mean convex  Figure 1A for two values of the coding level, a. Note that the nullcline curves upward in this plot, a consequence of the fact that we use Àh E rather than h E on the y-axis.
To find the m-nullcline-the set of points in h E À m space that satisfy Equation 21b-we proceed in two stages. First, we plot the right-hand side of Equation 21b versus m and look for intersections with the 458 line; these intersections correspond to points on the m-nullcline. Second, we vary h E and sweep out a curve in h E À m space; this curve is the full mnullcline. A typical plot versus m with h E fixed is shown in Figure 1B. There are three intersections with the 458 line, which means that the m-nullcline consists of three points at this particular value of h E : one with m ¼ 0 and two with m . 0. To find out how these three points move as we vary h E , we compute dm(h E )/dh E where the derivative is taken along the mnullcline; using Equation 21b, this is given by We are primarily interested in the sign of dm/dh E , which can be found by examining the signs of the numerator and denominator separately. For the denominator, note that the derivative of the term in square brackets is the slope of the curve in Figure 1B. Consequently, the denominator is negative for the intermediate intersection (where the slope is greater than 1) and positive for the upper intersection (where the slope is less than 1). The sign of the numerator depends primarily on the size of h E . If h E is small, so that both F E ðh E þ bmÞ and F E ðh E Þ lie on the convex part of the sigmoid, then the numerator is positive. If, on the other hand, h E is large, so that F E ðh E þ bmÞ and F E ðh E Þ lie on the concave part, then it is negative (see insets in Figure 1C). This gives us the following picture: when h E is small, so that the numerator in Equation 22 is positive, decreasing h E causes the two intersections in Figure 1B to move closer, and eventually to annihilate. When h E is large, on the other hand, so that the numerator is negative, increasing, rather than decreasing, h E causes the intersections to move closer, and eventually annihilate, this time for sufficiently large h E . Filling in the points away from the extrema, we see that the mnullcline is topologically equivalent to a circle ( Figure 1C). Finally we note that the line m ¼ 0 is also part of the nullcline, as can easily be seen from Equation 21b; this line is also included in Figure 1C.
In Figure 1D, we combine the h E -nullclines from Figure 1A and the m-nullcline form Figure 1C. Clearly there is always an equilibrium at m ¼ 0, corresponding to no active memories; i.e., corresponding to a null background. There are also two equilibria at m . 0, corresponding to active memories. In the section ''Stability analysis,'' we show that the one at larger m is stable. Importantly, this equilibrium can occur at small m, and thus low firing rate, something we will see more quantitatively in the next section, where we consider a specific example. Although not shown in Figure 1, the m-nullcline can shift far (D) The m-and h E -nullclines on the same plot. The intersections correspond to network equilibria. There are three equilibria: one at m ¼ 0, corresponding to the background state, and two at m . 0, corresponding to potential retrieval states. The one at m ¼ 0 and the one at large m are stable; the intermediate one is not. Consequently, only the large m equilibrium is observed during retrieval. Note that when the coding level, a, is small (dashed blue line), the retrieval state occurs at large m, and thus has a high firing rate. Only when a is large (solid blue line) is it possible to have a low firing rate during retrieval. doi:10.1371/journal.pcbi.0030141.g001 enough up so that m can be negative at equilibrium. When this happens, m ¼ 0 becomes unstable, which in turn implies that the background becomes unstable. We see this in the simulations: when b becomes too large, memories are spontaneously activated.
We can now see the critical role played by the coding level, a. In the limit a ! 0, the right-hand side of Equation 21a becomes almost independent of m. This makes the h Enullcline almost horizontal (dashed line in Figure 1D), so the only stable retrieval state occurs at large m, and thus high firing rate (the peak of the m-nullcline typically occurs near the maximum firing rate of the neurons, about 100 Hz; see next section). If, on the other hand, a is reasonably large, then the h E -nullcline can curve up and intersect the m-nullcine to the left of its highest point (solid blue line in Figure 1D). As just discussed, this intersection corresponds to the intermediate intersection in Figure 1B, which means it corresponds to low firing rate, and thus a biologically realistic retrieval state.
We end this section by discussing the conditions under which the nullclines in Figure 1D, which were derived from Equation 17, are the same as the nullclines for the full meanfield equations, Equation 47. The primary effect of the full set of equations is to couple h I to h E and m. One could, however, solve for h I in terms of h E and m, insert that solution into the equations for h E and m, and derive a new coupled set of equations that again involve only h E and m. This would, effectively, replace F E in Equation 17 with a more complicated function of h E and m. Examining Equations 47 and 48, we see that these manipulations would result in the following replacements, Retracing the steps that led us to Figure 1D, we see that if F E ðh E ;r E ðh E ; mÞÞ and F E ðh E þ bm;r E ðh E ; mÞÞ are quasi-sigmoidal functions of h E and m, we recover the nullclines in Figure 1D. Both of these conditions are likely to hold for real neurons, since increasing h E and m correspond to increasing excitatory drive. Thus, for neurons with reasonable gain functions, we expect Figure 1D to fully capture the shape of the nullclines.
An example: Nullclines for a simple gain function. As an illustrative example, we consider a specific form for the gain functions (the F Q ), and compute the resulting nullclines numerically. The form we choose is a rather standard one, where m max is the maximum firing rate of both excitatory and inhibitory neurons, which without loss of generality we take to be 100 Hz, H(x) is given by and r Q is an approximate standard deviation based on Equation 44, Before computing the nullclines for these gain functions, we introduce a transformation that changes the nullclines without changing the equilibria. Combining Equations 21a and 21b, we see that Equation 21a can be written Note that the right-hand side of Equation 25 is an increasing function of both h E and m, so the h E -nullcline based on Equation 25 has the same qualitative shape as the h Enullcline based on Equation 21a. This form is more useful than the one in Equation 21a, however, because we can immediately write down an expression for h E (m), Computing the nullclines is now a straightforward numerical task, and in Figure  The first thing we notice is that when b is sufficiently small (Figure 2A), the m-nullcline consists only of a line at m ¼ 0, which means that the only possible equilibria are at m ¼ 0, and so retrieval states are not possible. When b is slightly larger ( Figure 2B), the m-nullcline gains a second piece away from the line m ¼ 0. However, this second piece lies below both h Enullclines, so the only intersections are again at m ¼ 0, and retrieval is again not possible. The fact that there is no memory retrieval when b is small makes sense: b controls the connection strength among the neurons within each memory, so if it is too small there will not be enough recurrent connectivity to produce elevated firing.
For still larger b, there is an intersection with one of the h Enullclines-the one corresponding to low coding level ( Figure  2C). The stable equilibrium, which is the equilibrium with larger m, corresponds to memory retrieval (see the section ''Stability analysis''). Finally, at sufficiently large b, the system acquires an intersection with the h E -nullcline corresponding to high coding level ( Figure 2D). Again, the stable equilibrium is the one with larger m.
An important point is that the value of m at the retrieval state, and thus the firing rate of the foreground neurons, depends strongly on the coding level, a. For small a (dashed blue line), retrieval occurs near saturation, and thus at an unrealistically high firing rate. For larger a (solid blue line), the retrieval occurs at low firing rate, consistent with experiments (when a ¼ 0.05 and b ¼ 1.2, the equilibrium value of m is 20 Hz). This is exactly the behaviour we saw in the previous section.
As can be expected from these figures, increasing b even further would shift the intermediate intersection to negative values of m. In this regime the background becomes unstable. Again this makes sense: if the coupling among the neurons within a memory is too strong, they become spontaneously active. Examining Figure 1B, we see that this occurs when the slope of F E ðh E þ bmÞ À F E ðh E Þ with respect to m is 1 at m ¼ 0 (and, of course, h E is at its equilibrium value). The value of b at which this happens, denoted b max , is given by (see Equations 21b and 26). For the sigmoidal gain function used in this example (Equation 24), b max is given by The phase diagram for this model-a plot showing stability and, in the stable region, the firing rate of the foreground neurons-is shown in Figure 3. Storage capacity. In the above analysis, there was no way to determine how many memories could be embedded in a network, and thus no way to determine storage capacity. That's because we hid all effects of the quenched noise-the noise associated with the random elements of the connectivity matrix-in F E and F I (see Equation 17). However, the quenched noise can have a nontrivial effect, in two ways. First, within the context of the self-consistent signal-to-noise analysis, it changes both F E and F I , and thus modifies the nullclines. Second, and potentially more important, as we add memories we increase the number of preferred modes that can be activated in the network, and thus we increase the quenched noise. Either effect could cause memories to be active when they should not be, and inactive when they should be.
To quantify these effects, we note that both scale with the fluctuations associated with the memories that are not recalled on any particular trial. The size of these fluctuations can be found by computing the contribution of the memories to dĥ E , the fluctuating piece in Equation 12a. Examining the memory portion of the connectivity matrix, A ij , which is given in Equation 4, and noting thatb is proportional to K À1 E (Equation 10), we show in the section ''Mean-field equations'' that the variance of the quenched fluctuations associated with this term scale as p/K E (Equation 45). Intuitively, that is because when we sum the right-hand side of Equation 4 on j and l, there are (p À 1)K E terms: K E that come from the sum on j, and p À 1 that come from the non-activated memories in the sum on l. Each of these terms has variance that is Oð1=K 2 E Þ. Central limit type arguments then tell us that the variance of such a sum is on the order of ðp À 1ÞK E =K 2 E ' p=K E , where the approximation is valid if p is large. Consequently, there is a critical value of p/K E above which none of the stored patterns could be retrieved. Thus, the maximum number of memories in a network should scale linearly with K E . This is what we found in our simulations (see the section ''Computer Simulations'').
Unfortunately, the scale factor we found in our simulations ). The parameters were The h E -nullcline with a ¼ 0.001 is essentially a straight line, so memory retrieval occurs at a firing rate that is too high to be biologically realistic. The h E -nullcline with a ¼ 0.05, on the other hand, has strong upward curvature, so memory retrieval occurs at a much lower, and thus biologically plausible, firing rate. doi:10.1371/journal.pcbi.0030141.g002 was small, in that the maximum number of memories scaled as 0.01 K E . A natural question to ask, then, is: can the scale factor be improved by, for example, using different parameters in our network? In the rest of this section, we focus on the effect of the coding level, a, on the storage capacity. We choose the coding level because, at least in simple models, the storage capacity is inversely proportional to a [33,34,37]. We have already shown that as the coding level decreases, the foreground firing rate becomes large, so we cannot make a arbitrarily small. However, the minimum allowable value of a depends on the model. What we show below, though, is that even for models which exhibit realistic foreground firing rate at relatively low coding levels, the 1/a scaling of the storage capacity does not hold. This suggests that decreasing the coding level cannot be used to increase the storage capacity in realistic networks. Examining Equations 12a and 46, we see that the background neurons receive an input drawn from a Gaussian distribution with mean h E and standard deviationr E , while the foreground neurons receive input with larger mean, h E þ bm, and the same standard deviation,r E . When the standard deviation of these distributions,r E , is smaller than the separation between the means, the two populations are well separated ( Figure 4A) and memory recall is possible. The standard deviation, however, is an increasing function of p; see Equation 47d and note that p enters this equation only through the storage load, a, which is defined to be When a, and thus p, becomes large enough, the standard deviation is on the same order as the separation. At this point, the two distributions have a significant overlap with each other ( Figure 4B), and memory recall fails. Using this intuitive picture and Equation 47d, we can find the value of a for whichr E is on the order of the separation between the means; this should give us an estimate of the storage capacity, a max . Using Equation 47d and the fact that the means are separated by bm (see Figure 4), we see that this happens when where ; the lower boundary is determined numerically by finding the minimum value of b (for a given a) such that the m-nullcline and h E -nullcline intersect. The parameters are the same as in Figure 2: The standard deviation is much smaller than the distance between the means of the two distributions. In this regime, the two populations are well separated, there is no interference between them, and memory retrieval is supported.
Solving Equation 29 for a max then leads to If the background synaptic weights, J EE and J EI , were zero and there was zero background firing so that c 2 vanished, we would recover the 1/a scaling (in the small a limit) found in simpler models [33,3437]. With nonzero background synaptic weights, however, the capacity no longer scales as 1/a. Consequently, we expect that the maximum capacity cannot be improved much by using sparser codes.

Computer Simulations
Our mean-field analysis gave us two predictions. The first is that if the background synaptic weights, theJ, scale as K À1/2 , the foreground weights, A, scale as K À1 , and the coding level, a, is sufficiently high, then both the background and foreground neurons should operate in the balanced regime and the neurons should fire irregularly. The second prediction is that the number of memories that can be stored is proportional to the number of excitatory connections per neuron, K E .
To test these predictions, we perform simulations with large networks of spiking neurons. We start by finding, for a particular network size, parameters such that both foreground and background neurons exhibit irregular activity. We then increase the size of the network while scaling the synaptic weights according to the above prescriptions. If the larger networks continue to exhibit irregular activity, then our predicted scalings are correct. To test the relation between storage capacity and number of connections per neuron, we calculate the storage capacity for networks with different sizes. A linear relation would indicate a scaling consistent with our predictions.
Network model. Each neuron is modeled as a conductancebased quadratic integrate and fire (QIF) neuron. Dendritic trees and axonal arborizations are not considered. The spikes generated in any neuron immediately affect all the postsynaptic neurons connected to it. The membrane potential of neuron i of type Q, denoted V Qi , evolves according to Here s is the membrane time constant, s s is the synaptic time constant, V r and V t are the nominal resting and threshold voltages, V 0i determines the actual resting and threshold voltages, (V 0i is constant for each i, but as a function of i it's a Gaussian random variable with mean V 0i and standard deviation DV 0 ),J QR ij is the connection strength from cell j in population R to cell i in population Q, e E and e I are the excitatory and inhibitory reversal potentials, respectively, the notation j 2 R means sum over only those cells of type R, d(Á) is the Dirac d-function, t k j is the kth spike emitted by neuron j, andh Qex;i is the external input to neuron i of type Q. The external input is modeled as where the t k Qex are the times of the external spikes. These are taken to be Poisson at constant rate v Qex .
There are two features of these equations that are worth commenting on. First, the connection strengths,J have the same form as in Equation 6, except that we introduce an extra scaling factor so that connection strengths can be directly related to postsynaptic potential (PSP) size. Specifically, we use the fact that if neuron j spikes and neuron i is at rest, then the PSP generated at neuron i will have peak amplitudeJ QR ij V R where V R ¼ e R À V r ðs=s s Þexp½lnðs=s s Þ=ðs=s s À 1Þ ; see [58] for a derivation of this expression. This suggests that we should scale our connection strengths by V R , so we writẽ where c QR ij is the same binary random variable defined in Equation 7, d Q,R is the Kronecker delta, andJ QR andb in Equation 32 correspond to, but typically have different numerical values than, the ones in Equations 6 and 10. If V R is in mV, thenJ QR is the peak PSP, in mV, that occurs in a neuron in pool Q when a neuron in pool R fires (assuming the two are connected, the postsynaptic neuron is at rest, and b ¼ 0).
Our analytical results have been derived by assuming current-based neurons. However, it is possible to extend such analysis to a more realistic network of conductance-based neurons by noting that the effective connection strength in a conductance-based model is proportional to the PSP size [40,44,57]. Thus, for the network to operate in the balanced regime, we should have the following scalings, p}K; Note that the mean external excitatory input must be proportional to K 1/2 . Therefore, given Equation 31a and the scaling ofJ Qex in Equation 33b, the firing rate of the neurons that provide external input, v Qex , should scale as K. We performed simulations using three different networks, called Networks 1, 2, and 3, that differ in the number of neurons (they contain a total of 10,000, 20,000, and 30,000, respectively). In all three networks, c ¼ 0.15, so K is proportional to the total number of neurons in the network. Because of the scaling in Equation 33, the values ofJ QR ,b, v Qex , and p also differ. The parameters for the three networks are given in Table 1. Our goal in these simulations is to determine whether, as predicted by our mean-field analysis, the above scaling leads to behaviour that is independent of K and the firing of both foreground and background neurons is irregular.
Building a balanced network. Our first step in assessing our mean-field predictions is to build a network that operates in the balanced regime and supports retrieval states. To test whether a network is operating in the balanced regime, we rely on two indicators. One is that it exhibits irregular firing, quantified by the coefficient of variation (CV)-the ratio of the standard deviation to the mean interspike interval-and that the CV is independent of K. The second is that the mean excitatory and inhibitory firing rates scale linearly with the external input, as predicted by Equation 19. To test whether a network supports retrieval states, we simply activate a memory by bombarding all the neurons within a memory with excitatory input, and ask whether the memory stays active for several seconds. Very little fine-tuning was required to find a network that exhibited both balance and retrieval states: we simply chose reasonable peak PSPs, set the coding level, a, to 0.1, and increasedb until at least one memory was stored.
In Figure 5A, we show an example of the retrieval of a stored pattern for Network 1. The first 2 s in this figure consists of background firing; at t ¼ 2 s, neurons selective for one of the patterns receive an excitatory external input lasting for 100 ms; and at t ¼ 27.3 s, the same neurons receive an inhibitory external input, which again lasts for 100 ms. The blue line is the mean firing rate of the foreground neurons, the black line is the mean firing rate of the excitatory neurons (both foreground and background), and the red line is the mean firing rates of the inhibitory neurons.
Two points are worth mentioning. One is that the background firing rate in our simulations is lower than the background firing rate observed in studies of delay activity, which range from 1.5 to 8 Hz [15], although we should point out that the firing rates determined from extracellular recordings may be overestimated due to selection bias [59]. We could, however, achieve a higher background rate by increasing the excitatory external input; an example is shown in Figure 6, for which the network parameters are the same as Network 1 ( Figure 5A) except that the external input to excitatory and inhibitory neurons is five times higher, b is a factor of about two higher, and there is just one stored pattern instead of five. With the higher input, the background and foreground rates are in the range reported from neurons in, for example, anterior ventral temporal cortex [1,3] and entorhinal cortex [15].
The second point is that during retrieval, the mean firing rates of the excitatory and inhibitory neurons differ from the background rates; i.e., from the rates when no memories are activated. This appears to be inconsistent with the balance condition, which predicts that the mean firing rate during the activation of a memory is the same as that when the network is in the background state (see Equation 19). However, this prediction holds only in the limit of infinite connectivity. For finite connectivity, there are corrections, and they are particularly important when the firing rate is low [54]. For example, in Figure 5A the average excitatory activity increased from 0.28 Hz in the background to 1.07 Hz during retrieval (an increase of about 400%), whereas in Figure 6, where the background is higher, it increased from 1.06 Hz to 1.73 Hz (an increase of 60%). Thus, the increase in the mean excitatory firing rate during retrieval is reduced when the firing rate is higher. However, this is accompanied, at least in the parameter range we looked at, by a decrease in the storage capacity. Since we would like to study the scaling of storage capacity, we operate in the lower firing rate regime. A detailed search of parameter space is required to determine whether both high storage capacity and high background firing can be achieved.
In Figure 7A we show the CV versus firing rate, again for Network 1. Here and in what follows, the CV is calculated only for those neurons that emit at least five spikes during the 25 s period that the pattern is active. The data in Figure 7A fall into two clusters, one (blue dots) corresponds to  background neurons and the other (red crosses) to foreground neurons. The distributions of CVs and firing rates are shown in Figure 7B and 7C. The CV of both background and foreground neurons are on the order of 0.8, which indicates irregular firing. This suggests that the network is operating in the balanced regime. To further test for balance, in Figure 8A we plot the average excitatory and inhibitory firing rates versus the external input. As predicted by Equations 18 and 19, the relation is approximately linear.
Scaling of the parameters. To test our predicted scaling with the number of connections, we considered networks with two and three times the number of neurons and connections as in Network 1; these are Networks 2 and 3. At the same time, we scaledJ QR by K À1=2 , m Qex by K, and p by K (see Equations 33a-33c). The value ofb was set, as in Network 1, to the minimum value that results in retrieval of a single stored pattern. The 1/K scaling of b (Equation 10) gives us the values reported asb predicted in Table 1. The values found from simulations (b in Table 1) do not exactly follow the expected 1/K scaling:b is 20% too large in Network 2 and 40% too large in Network 3. As discussed in the section ''Retrieval states in the finite connectivity regime,'' this is because of finite K effects, and the trends we see here follow the trends predicted in that section.
Examples of stored memories in Networks 2 and 3 are shown in Figures 5B and 5C, the CV versus firing rate is shown in Figures 7D and 7G, and the distribution of background and foreground CV and firing rates during the 25 s period that the memory is active are shown in Figure 7E and 7F for Network 2 and Figure 7H and 7I for Network 3. These plots show that when the connection strengths are scaled properly, both the background and foreground neurons exhibit irregular firing, just as in Network 1. Finally, Figure 8B and 8C show the relationship between the external input and the firing rate of the inhibitory and excitatory populations. As we saw for Network 1, the firing rate of excitatory and inhibitory neurons are linearly related to external input, further evidence for the balanced regime. In theory, the lines should lie on top of each other; however, due to finite size effects, this does not happen. The fact that finite size effects are responsible for this deviation from the theory can be seen by noting that the lines corresponding to Network 2 and Network 3 are much closer to each other than Network 1 and Network 2.
Scaling of the maximum number of memories. Our last prediction is that the maximum number of memories should be linear in the number of excitatory connections, K E . To test this, for each of our three networks we increased the number of patterns, p, until the network failed to exhibit retrieval states. Specifically, we performed simulations as described in Figure 5, except that the memory was active for 4 s rather than 25 s. For each value of p, we activated, one at a time, either all p memories (if p was smaller than 20) or 20 memories (if p was larger). If the mean activity of the foreground neurons during the activation period was at least three times larger than the activity averaged over all the excitatory neurons, then that memory was said to be successfully retrieved.
The results of these simulations are shown in Figure 9A, where we plot the fraction of successful retrievals versus p/ KE for the three networks. Consistent with our predictions, the transition to a regime where none of the patterns could be retrieved occurs at approximately the same value of p/ KE for all three networks. Moreover, as one would expect, the transition for the largest network is sharper than for the others.
Although Figure 9A shows that p max scales linearly with K E , in these simulations N E also scales with K E , so this does not rule out the possibility that p max is proportional to N E rather than K E . To test for this, in Figure 9B we plot the fraction of successful retrievals versus p/ KE , but this time with K E fixed and N E varied. This figure shows that p max is proportional to K E , not N E , ruling out the N E scaling.

Discussion
In this paper we addressed two questions. First, can all the neurons in an attractor network-both background and foreground-exhibit irregular firing? And second, what is the storage capacity in networks of realistic spiking neurons? To answer these questions, we applied self-consistent signal-tonoise analysis to large networks of excitatory and inhibitory neurons, and we performed simulations with spiking neurons to test the predictions of that analysis.
Our primary finding is that two conditions must be met to guarantee irregular firing of both foreground and background neurons. The first is proper scaling with the number of connections per neuron, K: the strength of the background weight matrix must scale as K À1/2 and the strength of the structured part of the weight matrix (the part responsible for the memories) as K À1 . What this scaling does is guarantee ''balance,'' meaning the network dynamically adjusts its firing rates so that the mean input to a neuron is on the same order as the fluctuations, independent of K. This in turn guarantees that the degree of irregular firing is independent of K.
While balance is a necessary condition for irregular firing, it is not sufficient. That's because balance ensures only that the mean and fluctuations are independent of K, but does not rule out the possibility that the mean is much larger than the fluctuations, which would result in regular firing. To ensure that this does not happen, a second condition must be satisfied: the coding level, a, must be above some (Kindependent) threshold. This condition is needed to ensure that the coupling between background and foreground neurons is sufficiently strong to stabilize a low firing rate foreground state on the unstable branch of the m-nullcline (see Figure 1).
The analysis that led to predictions of irregular firing also quite naturally provided us with information about the capacity of attractor networks-the maximum number of patterns that could be stored and successfully retrieved. What we found, under very general conditions, was that this maximum, denoted p max , is linear in the number of excitatory connections per neuron, K E . This scaling relation has been observed in studies of simplified attractor networks [16,32,34], but, as discussed in the Introduction, those models did not include all the features that are necessary for a realistic recurrent networks. Thus, the analysis performed here is the first to show that the number of memories is linear in K E in biophysically plausible networks.

Scaling in Other Models, and the Importance of O(1) Input to the Foreground Neurons
Note that there are other types of scaling, different from what we proposed, which can result in irregular firing of both foreground and background neurons. What is critical is that the net input a foreground neuron receives from the other foreground neurons should be O(1). We achieved this by letting the structured part of the connection matrix (the second term in Equation 11a) be O(1/K), and using a coding level, a, that was O(1). However, this is not the only possible combination of connection strengths and coding levels, and in the two other studies that address both scaling and irregularity in memory networks [27,28], different combinations were used. In the model proposed by van Vreesjwik and Sompolinsky [27], the structured part of their connection matrix was a factor of K 1/2 larger than ours; to balance that, the coding level was a factor of K 1/2 smaller. In the model proposed by Renart et al. [28], the structured part of the synaptic weights was K times larger than ours, so their coding level had to scale as O(1/K). Whether such low coding levels are consistent with reality needs further investigation; however, data from studies conducted on selectivity of neurons to visual stimuli suggests that it is too low [30]. In addition to the very low coding level that these two models require, they also exhibit non-biologically high foreground firing rate. Nevertheless, the model of Renart et al. [28] does have one advantage over others: the foreground neurons are as irregular as, or even more irregular than, the background neurons, something our model does not achieve (see next section).

Not as Irregular as It Could Be
Although our simulations showed irregular activity, we found that the mean CV was only about 0.8. This is smaller than the values measured in vivo, which are normally close to, or slightly above, one [24,[50][51][52][53]. In addition, in our simulations the CV showed a small, but consistent, decrease with firing rate (see the left column in Figure 7). This is due to the fact that with the scaling that we chose, the fluctuations in the input current to foreground and background neurons are the same but the mean current to the foreground neurons is higher (see the section ''Fast fluctuations''). This decrease in the CV disagrees slightly with a study by Compte et al. [24], who found that the CV in prefrontal cortex does not depend on the mean firing rate, at least in a spatial memory task. While there are many possible reasons for this discrepancy, a likely one arises from the fact that the neurons in our network contained only two time scales, the membrane and synaptic time constants, and both were short: 10 ms for the former and 3 ms for the latter. Real neurons, however, have a host of long time scales that could contribute to irregularity [60]. In addition, in vivo optical imaging [61][62][63] and multielectrode [64] studies indicate that the background activity varies coherently and over long time scales, on the order of seconds, something we did not model. Both of these would increase the CV, although how much remains to be seen.
Although multiple time scales could certainly increase irregularity, it is not the only possible way to do this. As discussed in the Introduction and in the previous section, the model proposed by Renart et al. [28] also increases irregularity, and is consistent with the experimental results of Compte et al. [24]. However, it requires a very small coding level (a } 1/K), and fine-tuning of the parameters.

Subthreshold versus Suprathreshold Persistent Activity
In conventional models of persistent activity [14,22,29], the foreground activity necessarily lies on the concave part of the excitatory gain function, F E (h E ), whereas the background activity lies on the convex part. Since the inflection point of realistic gain functions is typically near the firing threshold [42,43], this type of bistability is called suprathreshold bistability [22,28]. Because the concave part of the gain function is typically at high firing rate, with suprathreshold bistability it is hard to have either low foreground firing rate or high CV. Consequently, there has been interest in understanding whether it is possible to have subthreshold bistability; that is, whether it is possible for both foreground and background solutions to lie on the subthreshold part of the gain function [28].
The model presented here can in fact show subthreshold bistability: as discussed in the section ''Reduced mean-field equations in the infinite K limit,'' increasing the coding level, a, brings the foreground firing rate very close to the background rate. Therefore, for sufficiently large a, the foreground state would be on the convex part of the transfer function. Our model, and the recently proposed model by Renart et al. [28], are the only ones that can show subthreshold bistability.

Bimodal Distribution of Firing Rates
One rather striking feature of our networks is that they all produce a highly bimodal distribution of firing rates: as can be seen in the first and third columns of Figure 7, the background neurons fire at a much lower rate than the foreground neurons-so much lower, in fact, that they form a distinct, and easily recognizable, population. This occurs because the patterns we store-the n l i -are binary, which makes the average input current to every neuron in the foreground exactly the same. This feature is potentially problematic, as the distinction between foreground and background rates observed in experiments is not nearly as striking as the one in Figure 7 [65]. However, this feature is not essential to our analysis, for two reasons. First, as discussed in the section ''Building a balanced network'' (see especially Figure 6), we deliberately made the background firing rate low to increase the capacity. Second, it is easy to extend our analysis to real valued patterns in which the elements of the n l i are drawn from a continuous distribution [34]. Under this, more realistic, scenario, it should be possible to match the statistics of the response seen in the cortex. This will be the subject of future work.

Fine-Tuning of the Weights
In our model, every time a new pattern is learned, the weights change by an amount proportional to K À1 . This is a factor of K À1/2 smaller than the background weights. Since weight changes are unlikely to be under such fine control, it is natural to ask whether errors during learning will lead to a major reduction in storage capacity. The answer, of course, depends on the size of the errors. In the section ''Fine-tuning in the learning rule,'' we show that errors can be larger than the weight changes by a factor of (K/p) 1/2 , with only a small change in storage capacity. More specifically, every time a pattern is learned, noise of O((Kp) À1/2 ) can be added to the synaptic strength, and the network will retain its ability to store and recall patterns.
Although this result tells us that the noise in the weight changes can be large compared with the structured part, the fine-tuning problem is not entirely eliminated: the noise must still be a factor of p 1/2 smaller than the background weights. Because of the low storage capacity found in these networks (at most 2.5% [23]), even when K is as large as 10,000, 1/p 1/2 is on the order of 6%. It seems plausible that biological machinery has evolved to achieve this kind of precision. However, for networks with larger capacity, the requirements on the precision of the weight would be more stringent.
It is also possible to have a probabilistic learning rule for which the changes in the weights are on the same order as the background weight, but this decreases the capacity significantly, by a factor of ffiffiffiffi K p (see the section ''Fine-tuning in the learning rule,'' Equation 78; we thank Carl van Vreeswijk for pointing this out). Although this probabilistic learning rule guarantees a balanced state with irregular background and foreground firing, it has the drawback that the storage capacity scales as ffiffiffiffi K p rather than K.

Low Storage Capacity
Although we showed that p max } K E , we did not compute analytically the constant of proportionality. In our simulations, this constant was small: from Figure 9, p max is about 0.01 K E , which means that for K E ¼ 10,000 we can store only about 100 patterns. It is important, though, to note that we made no attempt to optimize our network with respect to other parameters, so the constant of proportionality 0.01 is unlikely to be a fundamental limit. In fact, Latham and Nirenberg [23] were able to store about 50 patterns in a network with 2,000 excitatory connections, 2.5 times larger than our capacity. Interestingly, the only substantial difference between their network and ours was that in theirs the  There is a critical value of a, above which the fraction of successful runs is zero; this is the storage capacity a max . The transition at a max is sharp for K E ¼ 3,600 but smoother for K E ¼ 2,400 and K E ¼ 1,200, due to finite size effects. The fact that a max is almost the same for all three values of K E implies that the maximum number of patterns that could be stored and retrieved, p max , is linear in K E . (B) The fraction of successful runs versus the storage load, a ¼ p/K E , for three networks with all parameters, except for the total number of neurons in the network, is equal to those of Network 1. This figure shows that increasing the size of the network does not change p max . doi:10.1371/journal.pcbi.0030141.g009 background activity was generated by endogenously active neurons rather than by external input.
Can we further increase the scaling factor? One potential mechanism is to decrease the coding level, a, since, at least in simple models [33,34,37], the maximum number of patterns that could be stored and retrieved is inversely proportional to the coding level. But, as we showed in the section ''Storage capacity,'' realistic networks do not exhibit this 1/a scaling. Consequently, sparse coding cannot be used as a way to improve the storage capacity in our network. Simplified models also suggest that one can increase the storage capacity by a factor of 3-4 by using other schemes, such as non-binary patterns [34], or spatially correlated patterns [66]. Whether these techniques can be extended to the kind of network we have studied here is not clear, and requires further investigation. However, an increase beyond a factor of 3-4, to a capacity above about 0.1, seems unlikely within this class of networks.
In any case, there is a limit to the number of memories that can be stored in a single attractor network with a fixed number of connections per neuron, no matter how many neurons in the network. This suggests that, in order to make the best use of the existing connections, realistic working memory systems must be composed of interconnected modules. In this paradigm, each module would consist of an attractor network [67][68][69]. Such modular structure naively suggests a combinatorial increase in storage capacity; however, understanding how to achieve such an increase has proved difficult. For simple models whose storage capacity could be calculated analytically, either no increase in the storage capacity [67] or a modest increase [69] was found. It is yet to be determined how modular networks could be implemented in realistic networks of spiking neurons, and what their storage capacity would be.

Materials and Methods
Fast fluctuations. The starting point for essentially all of our analysis is Equation 1, which, when combined with Equation 2, tells us that the time evolution of the firing rate of each neuron is purely a function of the firing rates of the other neurons. At a microscopic level, though, each neuron sees as input a set of spikes, not rates. However, for our model, rate-based equations do apply, as we show now.
In a spiking, current-based network, the input, h Qi (t), to the ith neuron in population Q has the form where t k j is the time of the kth spike on the jth neuron, f R (t), which mimics the PSP, is a non-negative function that integrates to 1 and vanishes for t , 0 and t large (greater than a few tens of ms). In a slight departure from our usual convention, R can refer to external input (R ¼ ex) as well as excitatory and inhibitory input (R ¼ E,I).
Our first step is to divide the input, h Qi , into a mean and a temporally fluctuating piece. The mean, which is found by timeaveraging the right-hand side of Equation 34a and using the fact that f R (t) integrates to 1, is simply where hÁÁÁi t represents a temporal average. The temporally fluctuating piece of the input can then be written The fluctuations, dh Qi , have zero mean by construction, and their correlation function, C Qi (s), is defined to be Assuming that h Qi is Gaussian (which is reasonable if there are a large number of neurons and they are not too correlated), then the firing rate depends only on the mean, hh Qi (t)i t , and the correlation function, C Qi (s). If the correlation function is independent of i, then the only i-dependence in the firing rate is through the mean input, and we recover Equation 1. What we now show is that, for our model, C Qi does not depend on i.
To understand the behaviour of C Qi , we express it in terms of dS R j ðtÞ; using Equation 36a, we have Under the assumption that the neurons are very weakly correlated, only the terms with j ¼ j9 survive, and this expression simplifies to Let us focus on the sum on j on the right-hand side of this expression. For Q 6 ¼ E or R 6 ¼ E, this sum is given by (see Equations For sparsely connected networks, c QR ij is independent of dS R j ðtÞ. Consequently, we can replace c QR ij on the right-hand side of Equation 38 by its average, c, and the right hand side becomes independent of i.
For Q ¼ R ¼ E, the situation is more complicated, asJ EE ij has an additional dependence on A ij , the structured part of the connectivity. Specifically using Equation 6a and again replacing c EE ij by its average, As discussed in the section ''Mean-field equations,'' A ij receives contributions from two sources: the p À 1 patterns that are not activated, and the one pattern that is. The non-activated patterns are not correlated with dS j , so they can be averaged separately in Equation 39, and thus do not produce any i-dependence. The activated pattern, on the other hand is correlated with dS j . However, the connection strength for the one activated pattern is smaller thanJ EE by a factor of K À1/2 (see the section ''Strong synapses and the balanced condition''). Consequently, in the high connectivity limit, we can ignore this contribution, and the right-hand side of Equation 39 is independent of i. This in turn implies that C Qi depends only on Q.
The upshot of this analysis is that the only i-dependence in the firing rate comes from hh Qi (t)i t . Moreover, comparing Equations 2 and 35, we see that hh Qi (t)i t is exactly equal to h Qi , the input current to the firing rate function, F Q , that appears in Equation 1. Thus, for the model used here, the rate-based formulation is indeed correct. What we do not do is compute F Q , as that would require that we compute the correlation function, C Q (s), self-consistently, which is nontrivial [44]. However, our results depend very weakly on the precise form of F Q , so it is not necessary to have an explicit expression for it.
Mean-field equations. In this section, we derive the mean-field equations for the model described in the section ''Model.'' As discussed in the main text, the derivation of these equations revolves around finding the distributions of dĥ Ei and dh Ii , the fluctuations around the mean excitatory and inhibitory synaptic input (both quantities are defined implicitly in . The main assumption we make is that dĥ Ei and dh Ii are zero mean Gaussian random variables, so all we need to do is find their variances self-consistently. In addition, primarily for simplicity (and because it is reasonable in large networks in the brain), we assume that the number of connections is small compared with the number of neurons, so c ( 1.
Our first step is to simplify the expressions for our main order parameters, m E , m, and m I . In the context of the self-consistent signalto-noise analysis, ''simplify'' means ''replace sums by Gaussian integrals.'' To see how to do this, note that, for any function g, where Var[Á] indicates variance, exact equality holds in the N E ! ' limit (but approximate equality typically holds when N E is only a few hundred), and where the average over n is with respect to the probability distribution given in Equation 5.
To complete Equation 40, we need the variance of dĥ Ei and dh Ii . It is convenient to break the former into two pieces, dĥ Ei ¼ dh Ei þ dh mi , where the first, dh Ei , is associated with the background neurons, and the second, dh mi , is associated with the foreground neurons (both will be defined shortly). Then, examining Equations 11-15, and performing a small amount of algebra, we find that and Here d l,v is the Kronecker delta; it is 1 if l ¼ m and zero otherwise. In addition, for notational convenience, we have returned the superscript ''1'' to n i . For the rest of the section, we will use n 1 i and n i interchangeably.
Let us focus first on the contribution from the background, Equation 41. Since c QR ij is equal to c on average, the mean of both terms on the right hand side of Equation 41 is zero. Moreover, these terms are uncorrelated, so their variances add. The variance of the QRth term is then where the angle brackets represent an average over the distribution of c QR ij . Because c QR ij and c QR ij 9 are independent when j 6 ¼ j9, only terms with j 6 ¼ j9 produce a nonzero average. Thus, all we need is the variance of c QR ij À c, which is given by (the last approximation is valid because, as mentioned above, we are assuming c ( 1). Performing the sums over j and j9 and collecting terms, we have The term on the right-hand side, hm 2 R i, is the second moment of the firing rate of the neurons in pool R. Inserting Equation 43 into 41, we find that The last quantity we need is the variance of dh m . A naive approach to computing it proceeds along lines similar to those described above: assume all the terms in the sum over j and l in Equation 42 are independent, so that the variance of dh m is just pN E (the number of terms in the sum) times the variance of each term. This yields, with rather loose notation for averages and ignoring the OðK À1 E Þ correction associated with l ¼ 1, All the averages in this expression are straightforward: hðc EE ij Þ 2 i ¼ c, hn 2 i ¼ a, h(n À a) 2 i ¼ a(1 À a), and hm 2 E i was defined in Equation 43. Putting all this together and defining q 2 to be the variance of dh m , we have While Equation 45 turns out to be correct, our derivation left out a potentially important effect: correlations between the patterns, ðn l j À aÞ, and the firing rates, m Ej in Equation 42. These correlations, which arise from the recurrent feedback, turn out to scale as c, and so can be neglected [32,35,70,71]. Rather than show this here, we delay it until the end of the section (see the section ''Loop corrections vanish in the small c limit'').
To write our mean-field equations in a compact form, it is convenient to define the total excitatory variance, Then, combining Equations 3, 40, 44, and 45, the mean-field equations become The functions F E and F I that we used in Equation 17 are equivalent to the ones defined in Equation 48, although we had suppressed the dependence on the standard deviation and dropped the superscript. Equation 47 constitutes our full set of mean-field equations. A key component of these equations is that the number of memories, p, enters only through the variable a, which is p/K E . Thus, the number of memories that can be embedded in a network of this type is linear in the number of connections.
Loop corrections vanish in the small c limit. To correctly treat the loop corrections in our derivation of the variance of dh m , we need to be explicit about the correlations between the patterns, ðn l j À aÞ, and the firing rates, m Ej , in Equation 42. We start by defining the idependent overlap, m l i , as Inserting this into Equation 42 leads to Each of the terms m l i is a Gaussian random variable whose variance must be determined self-consistently. This can be done by inserting Equation 3 into Equation 50 to derive a set of nonlinear equations for the m l i . There are two types of terms to consider: the activated memory, for which l ¼ 1, and the non-activated memories, for which l 6 ¼ 1. However, in the large p limit we can safely ignore the one term corresponding to l ¼ 1. Thus, considering the contributions from memories with l 6 ¼ 1, we have Taylor expanding around m l j ¼ 0 and defining where a prime denotes a derivative, we have We can write Equation 52 in matrix form as where I is the identity matrix, the ith component of m l is equal to m l i , and the matrices K l and N l are given by To solve Equation 53, we need to invert I À K, in general a hard problem. However, what we show now is that K has only one O(1) eigenvalue, with the rest OðK À1=2 E Þ. This allows us to write the inverse in terms of a single eigenvector and adjoint eigenvector, a simplification that allows us to perform the inversion explicitly.
The spectrum of the random matrix, K l , is determined primarily by the mean and variance of its components [72]. In the large N E limit, these are given by where hÁÁÁi ij indicates an average over i and j, and we used the fact that n l j and F l9 Ej are independent. Given that K l is an N E 3 N E matrix, the fact that the mean and variance of its elements are OðN À1 E Þ and O((K E N E ) À1 ), respectively, implies that it has one eigenvalue that is O(1) and N E À 1 eigenvalues that are OðK À1=2 E Þ [72]. Letting v k and v y k be the eigenvector and adjoint eigenvector of K l whose eigenvalue is k k , we can solve Equation 53 for m l , where ''Á'' represents dot product. Letting k ¼ 0 correspond to the O(1) eigenvalue and explicitly separating out this component, the expression for m l becomes Since v 0 and v y 0 are vectors whose components are all the same, without loss of generality we can choose v 0 ¼ (1,1,. . .1)/N E and v y 0 ¼ ð1; 1; :::; 1Þ. Combining this choice with Equation 55 and using Equation 54b for N l , we have We are now in a position to return to Equation 51 and compute the variance of dh m (which, recall, is denoted q 2 ). Treating, as usual, all the terms in Equation 51 as independent, we have To compute hm l2 i i n;z we use Equation 57 and the fact that the offdiagonal elements average to zero, and we find that To derive this expression, we again used h(n À a) 2 i ¼ a(1 À a). Our final step is to insert Equation 59 into Equation 58. Ignoring the two terms in brackets in Equation 59 that are a factor of c smaller than the first, and using the fact that hF 2 E ðh E þ bnm þr E zÞi n;z ¼ hm 2 E i n;z , this leads to the expression for q 2 given in Equation 45. Consequently, loop corrections vanish, and we can use our naive estimate for the variance of dh m .
Ignoring the two terms in brackets in Equation 59 is strictly correct for infinitely diluted networks; i.e., networks with c ! 0. When c is nonzero but small, the terms in the brackets can be ignored safely unless k 0 ! 1. However, as we now show, k 0 ! 1 is precisely the point where the background becomes unstable. Thus, it is not a regime in which we can operate.
The significance of the limit k 0 ! 1 can be seen by replacing When the largest eigenvalue of K l exceeds 1, the unactivated memories become unstable, and retrieval of just one memory is impossible. As discussed above, the largest eigenvalue of K l is k 0 . Consequently, loop corrections are necessarily important (no matter how dilute the network is) at precisely the point where the unactivated memories, and thus the background, become unstable.
Stability analysis. To determine stability, we need to write down time-evolution equations for the order parameters, and then linearize those around their fixed points. For m E , m I , and m, which are linear combinations of the firing rates, this is straightforward-we simply insert their definitions, Equation 16, into the time-evolution equations for the individual firing rates, Equation 1. For the variances, r 2 E and r 2 I , the situation is much more difficult, as these quantities do not admit simple time-evolution equations [73]. Fortunately, we expect the effects of the variances to be small-as discussed in the main text, their primary effect is to smooth slightly the gain functions, something that typically (although presumably not always) stabilizes the dynamics. Alternatively, if we assume that the variances are functions of m E , m I , and m, (meaning we give them instantaneous dynamics), we can rigorously neglect them. This is because derivatives of the gain functions with respect to m E and m I are large, on the order of K 1/2 , while derivatives with respect to the variances are O(1). Thus, as a first approximation, we will ignore these variables, and consider only the dynamics of m E , m I , and m. Because of this approximation, we expect our stability boundaries to be off by a small amount.
Combining Equation 1 and Equation 47, the time-evolution equations for m E , m I , and m may be written where the notation / a,b indicates a derivative of / a with respect to the argument specified by b (for example, / E,I ¼ @/ E /@m I and / I,M ¼ @/ I / @m). Since / I is independent of m (which means / I,m ¼ 0), the equation for the eigenvalues, denoted k, becomes 0 ¼ ðð/ E;E À 1Þ À s E kÞðð/ I;I À 1Þ À s I kÞ À / E;I / I;E Â Ã ðð/ m;m À 1Þ À s E kÞ þ / I;E / m;I À ðð/ I;I À 1Þ À s I kÞ/ m; Equation 63 is a cubic equation in k, and thus not straightforward to solve. However, in the large K limit it simplifies considerably. That's because derivatives with respect to m E and m I are O(K 1/2 ), which follows because the /'s depend on m E and m I through h E and h I , and the latter are proportional to K 1/2 (see Equation 13). Defining the O(1) quantities Examining Equation 64, it follows that if the eigenvalue, k, is O(K 1/2 ), then the term / m,m À 1 and the last term in brackets can be neglected. There are two such eigenvalues, and they are given by À 4½ðs E s I Þ À1 ð/ E;E / I;I À / E;I / I;E g 1=2 2 : Both eigenvalues are negative if Equation 66 reduces to where prime denotes a derivative.
Comparing Equations 62 and 49, we see that / E,m ¼ a/ m,m , which leads to This expression strongly emphasizes the role of the coding level, a: if it were zero, the only stable equilibria would be those with / m,m , 1, which would imply high firing rates for foreground neurons (see Figure 1B). Although Equation 67 tells us the stability of an equilibrium, it is not in an especially convenient form, as it does not allow us to look at a set of nullclines and determine instantly which equilibria are stable and which are not. However, it turns out that it is rather easy to determine the sign of k m for a given set of nullclines simply by looking at them. To see how, we make use of the expressions for / E and / m (Equations 62a and 62c) to reduce the right-hand side of Equation 67 to an expression with a single derivative. Our starting point is the definition where h E (m) is given by Equation 26; the solutions of the equation W(m) ¼ m correspond to network equilibria. The advantage of this one-dimensional formulation is that, as we show below, the condition k m , 0 is equivalent to dW/dm , 1. Thus, by plotting the function W(m) versus m and looking at its intersections with the 458 line, we can find the equilibrium values of m, and, more importantly, we can easily determine which of them is stable and which is unstable.
To show that dW/dm , 1 is equivalent to the condition k m , 0, we note first of all that where, recall, a prime denotes a derivative. By combining these expressions with Equation 67, and performing a small amount of algebra, the condition k m , 0 can be written To see how this compares to dW/dm, we use Equation 68 to write Then, using Equation 25, which tells us that Comparing Equations 69 and 70, we see that the condition dW/dm , 1 is equivalent to k m , 0. Thus, it is only when W(m) intersects the 458 line from above that, the equilibrium is stable. Since W(m) is bounded, if there are three equilibria, the smallest one must be stable, the middle one unstable, and the largest one again stable. Thus, we can look at the nullcline plots and immediately determine stability (see below and Figure 10).
As an example, we revisit Figure 2. In terms of our specific form for the gain functions, Equation 23, and with h E (m) given by Equation 26, the equation for m becomes This equation is solved graphically in Figure 10A where we plot W(m) versus m for the same values of b used in Figure 2 and with a ¼ 0.005. Intersections with the 458 line correspond to solutions of Equation 71, and thus to network equilibria. As we saw in the sections ''Reduced mean-field equations in the infinite K limit'' and ''An example: Nullclines for a simple gain function,'' the main factor that determines the number and location of the intersections, and thus the ability of the network to exhibit retrieval states, is b. For b ¼ 0.1 and 0.25, there is just one intersection at m ¼ 0, while for intermediate values of b, b ¼ 0.5 and 1.2, two additional intersections appear. Increasing b even further moves one of the solutions to negative m and destabilizes the background, but this is not shown. We can now easily see that the curves in Figure 10A with b ¼ 0.1 and 0.25 have a single stable intersection at m ¼ 0 (meaning that the solutions with m ¼ 0 in Figures 2A and 2B are stable); the curves with b ¼ 0.5 and b ¼ 1.2 have two stable intersections, one at m ¼ 0 and one at large m (and thus the solutions at m ¼ 0 in Figure 2C are stable, those at intermediate m are unstable, and those with large m are again stable).
Although we see bistability, the firing rate for the retrieval state is unrealistically high-on the order of 100 Hz, near saturation. As discussed in the main text, we can reduce the firing rate by increasing a. This is done in Figure 10B, where we plot W(m) versus m but this time for a ¼ 0.05 and b ¼ 1.2. Again there are three intersections (corresponding to the three intersections between the m-nullcline with b ¼ 1.2 and the h E -nullcline with a ¼ 0.05 in Figure 2C). With this higher value of a, the upper intersection is now in a biologically realistic range.
Retrieval states in the finite connectivity regime. When we performed network simulations, we found that the memory strength,b, did not exhibit exactly the predicted 1/K scaling. Here we ask whether the departure from predictions that we observed can be explained by finite K corrections. These corrections, as we will see shortly, are on the order of K À1/2 . Since in our simulations K is as small as 1,500, these corrections are potentially large.
Our starting point is the exact set of reduced mean-field equations, which is found by combining Equations 18 and 19, When K is large we can solve these equations by perturbing around the K ! ' solutions, which we denote h E0 , m 0 , and h I0 (these are the solutions to Equation 21). The zeroth step in this perturbation analysis is to replace h E and h I by h E0 and h I0 where they appear in brackets (and thus multiply K À1/2 ). This gives us a new set of equations, For the inhibitory firing rate, it is easy to see the effect of finite K: h I is shifted relative to h I0 by an amount proportional to dm I . Only slightly more difficult are h E and m, for which we have to consider how dm E affects the nullclines. Fortunately, only the h E -nullcline is affected, and we see that it shifts in a direction given by the sign of dm E . In particular, (We consider Àh E since, by convention, we plot our nullclines in a space with Àh E on the y-axis.) Thus, if dm E is positive then the h Enullcline shifts down relative to h E0 , while if it is negative the nullcline shifts up. In our simulations we set b to b min , the minimum value of b that allows retrieval of one memory. To determine how K affects b min , then, we need to know how to adjust b so that we keep the grazing intersection as K changes. Fortunately, the h E -nullcline depends on K but not on b, and the m-nullcline depends on b but not on K. Thus, all we need to know is how the m-nullcline changes with b. Using Equation 72b, it is easy to show that at fixed m, The numerator in this expression is clearly positive, and, for equilibria to the left of the peak of the m-nullcline, the denominator is also positive (see the section ''Stability analysis''). Thus, increasing b causes the m-nullcline to move up. Combining Equations 74 and 75, we have the following picture, dm E decreases ! h E À nullclline moves up ! b min increases where ''up'' corresponds to movement in the m À(Àh E ) plane. To complete the picture, we need to know how dm E depends on K. From Equation 73, we see that dm E } K À1/2 [J II h E0 À J EI h I0 ] ¼ K À1/2 [ÀjJ II j h E0 þ jJ EI jh I0 ]. Thus, whether dm E is an increasing or decreasing function of K depends on whether jJ II jh E0 is larger or smaller than jJ EI jh I0 . However, as we have seen, typically h E is negative. Thus, we expect dm E to be proportional to K À1/2 with a positive constant of proportionality, which means that dm E is a decreasing function of K. Combining that with the above picture, we conclude that when K increases, b min also increases. This is shown explicitly in Figure 11. Moreover, it was exactly what we saw in our simulations: b min (b in Table 1) was larger than predicted when we increased K (compareb withb predicted in Table  1). Fine-tuning in the learning rule. In the model described here, the structured part of the synaptic weights scale as K À1 , whereas the background scales as K À1/2 . This appears to require fine-tuning, since adjustments to the weights during learning of the attractors have to be a factor of K 1/2 times smaller than the background weights; a factor that can be as high as 100.
The first question to ask, then, is: exactly how big is the fine-tuning problem? In other words, how much noise can we add to the learning rule without having a huge effect on the storage capacity? This can be answered by considering a learning rule in which the weight changes during learning a pattern are not quite perfect. Specifically, let us consider the following modification of Equation 4, where the g l ij are zero-mean, uncorrelated random variables with variance r 2 g . The additional noise in this learning rule increases the variance of the quenched noise by an amount K E pr 2 g . As a result, if r 2 g ; OððK E pÞ À1 Þ; ð77Þ the effect on storage capacity is an O(1) increase in the quenched noise, and thus the storage capacity still scales as K E . With the scaling in Equation 77, weight changes during learning of each pattern is a factor of p 1/2 smaller than the background weights, and therefore the amount of fine-tuning depends on how many patterns are stored. Because of the low storage capacity found in these networks (at most 2.5% [23]), even when K is as large as 10,000, p À1/2 is on the order of 6%.
We should also point out that it is possible for the weight changes associated with the structured part of the connectivity to be on the same order as the background, although at the expense of storage capacity. Let us consider a third learning rule in which each synapse has a probability q of changing its value during learning, where the q l ij are Bernoulli variables; q l ij ¼ 1 with probability q and 0 with probability 1 À q. Let us define the coupling strength slightly differently than in Equation 10, where, as usual, b ; O(1). With this definition, the mean memory strength, hb9q l ij i, is again b/K E a(1 À a), as in Equation 10. But by setting q}K À1=2 E , the synaptic weight change-if there is one-is OðK À1=2 E Þ, just as it is for the background weights. However, there is a major drawback: as is easy to show, the variance associated with the structured part of the connectivity increases by a factor of K E , so the maximum number of patterns scales as p max }K 1=2 E rather than K E . We thus use Equation 4 for A ij in all of our analysis.