Balanced networks under spike-time dependent plasticity

The dynamics of local cortical networks are irregular, but correlated. Dynamic excitatory–inhibitory balance is a plausible mechanism that generates such irregular activity, but it remains unclear how balance is achieved and maintained in plastic neural networks. In particular, it is not fully understood how plasticity induced changes in the network affect balance, and in turn, how correlated, balanced activity impacts learning. How do the dynamics of balanced networks change under different plasticity rules? How does correlated spiking activity in recurrent networks change the evolution of weights, their eventual magnitude, and structure across the network? To address these questions, we develop a theory of spike–timing dependent plasticity in balanced networks. We show that balance can be attained and maintained under plasticity–induced weight changes. We find that correlations in the input mildly affect the evolution of synaptic weights. Under certain plasticity rules, we find an emergence of correlations between firing rates and synaptic weights. Under these rules, synaptic weights converge to a stable manifold in weight space with their final configuration dependent on the initial state of the network. Lastly, we show that our framework can also describe the dynamics of plastic balanced networks when subsets of neurons receive targeted optogenetic input.

We agree with this comment, and have rewritten the introduction to better relate our work to the existing literature on balanced networks and synaptic plasticity.
Correlated activity can also produce changes in synaptic weights [29,30]. For instance spike-time dependent plasticity (STDP), is driven by patterns in the timing of pre-and post-synaptic spikes [31,32]. However, we still lack a general theory that relates STDP to changes in neural activity, and the resulting neural computations. Hence, much of the theoretical analysis of the effects of STDP relies on simulations [29,[33][34][35]. Analytical treatments have been proposed for a number of cases, initially describing the impact of STDP on individual synaptic weights [36]. These results have been extended to small networks [34], and networks of Poisson neurons [37][38][39][40]. Other work provided analytical treatments of specific plasticity rules, such as homeostatic inhibitory plasticity [41,42]. Using linear response and motif resuming techniques [43], Ocker et al. developed a theory describing the evolution of mean weights in recurrent neural networks of noisy integrate-and-fire neurons under STDP [44]. This approach relies on the assumption that the input to individual cells is dominated by white noise and that local synaptic input is comparatively weak. Related results were obtained by treating neural firing as a Poisson process [37][38][39]45].
Here, we develop a complementary theory describing the evolution of synaptic weights and associated mean rates in tightly balanced networks in both correlated and asynchronous states. We combine the mean-field theory of firing rates and correlations in balanced networks [13,14,23,24,[46][47][48] with an averaging approach assuming a separation of timescales between changes in spiking activity, and the evolution of synaptic weights [30]. We show how the weights and the network dynamics co-evolve under different classical rules, such as Hebbian plasticity, Kohonen's rule, and a form of inhibitory plasticity [31,32,41,49,50]. In general, the predictions of our theory agree well with empirical simulations. We also explain when the mean-field theory fails, leading to disagreements with simulations, and we develop a semi-analytic extension of the theory that explains these disagreements. Our theory thus allows us to understand mean synaptic weight dynamics, and predicts how synaptic plasticity impacts network dynamics.
We find that spike train correlations, in general, have a mild effect on the synaptic weights and firing rates, in agreement with previous work [44,51]. We also show that for some STDP rules, synaptic competition can introduce correlations between synaptic weights and firing rates, resulting in the formation of a stable manifold of fixed points in weight space, and hence asymptotic weight distributions that depend on the initial state. Finally, we apply this theory to show how inhibitory STDP [41] can lead to a reestablishment of an asynchronous, balanced state that is broken by optogenetic stimulation of a neuronal subpopulation [52]. We thus extend the classical theory of balanced networks to understand how synaptic plasticity shapes their dynamics.

Referee A
Akil et al. develop a mean field theory (with some important semi analytical corrections) for the co-evolution of synaptic weights and first and second order firing rate statistics in spiking neuronal networks. This is potentially an important contribution to both theoretical and experimental neuroscience, where the understanding of the effects various plasticity rules have at the network level are lacking. This paper puts our understanding of these effects on more firm analytical grounds. However, I think that the results could be strengthened significantly if the authors leverage the theory they developed to address more fundamental problems than the applications included in the current version. Beyond that, the current version of the manuscript is missing some key references to and discussion of both theoretical and experimental literature that would help readers understand the importance of the results and their limitations. A fundamental problem in theoretical neuroscience is the stability of neuronal network activity subject to activity dependent plasticity. It is clear-and the authors mention it in their paper (l. 217)-that Hebbian excitatory plasticity leads to instabilities. The theory that the authors develop is exactly the analytical tool needed to understand what are minimal and/or realistic modifications of this rule which do lead to stable dynamics. More specifically-can the reduced, mean-field, description of the network+plasticity rule be leveraged to understand under what circumstances a Hebbian EE STDP rule could be stabilized? Alternatively, are there families of plasticity rules which are guaranteed to not stabilize the dynamics?
Thank you for the comments. We have now expanded the Discussion and the Appendix to bring out these points better. The following text now appears in the Discussion: One of the most important issues in understanding neural dynamics is establishing conditions under which the network remains active, yet stable as synaptic weights change. The theory we developed can help us address this questions, but it does have limitations. Since we used a mean-field approach, we can only capture first moments. While mean weight stability may not imply stable network dynamics (consider the case when weight variance diverges in Classical Hebbian STDP in S1 Appendix), instability in the mean weights does imply that the network is also unstable.
We have also expanded further on these results in S1 Appendix.
Classical Hebbian STDP leads to unstable dynamics. If we assume that the EE weights change according to: Mean-field theory predicts stable dynamics when the system has a stable fixed point determined by dJ ee /dt = 0. However, simulations show that network dynamics are unstable due to the divergence of individual weights to ±∞ (Supporting Fig. 10 A-D), so that the mean stays finite while the variance across the population diverges. Therefore, it is possible that the mean-field theory predicts stable dynamics, despite instabilities in higher moments that can cause weights to blow up. On the other hand, network dynamics will be unstable if the mean-field theory predicts that the first moment diverges.
In simulations Hebbian STDP is frequently stabilized by imposing hard lower and upper bounds on synaptic weights (Supporting Fig. 10 E,F). With appropriately chosen bounds, the network is stabilized and balance is achieved (Supporting Fig. 10 G,H). We can use mean-field theory to find alternative ways to modify learning rules in order to stabilize network dynamics: One such modification is to multiplying the LTP term in Eq. (1) by a constant (J max ∼ O(J ee jk )) and the LTD term by the current value of the synaptic weight, J ee jk , in what we call the "Weight-dependent Hebbian STDP rule" (See Table 1 in main text). This can be interpreted as a constraint on synaptic weights due to limitations on the size and strength of a synapse. Excitatory synaptic weights then converge to a fixed point and the network remains in the balanced regime (Supporting Fig. 10 I-L). This approach can be generalized to other plasticity rules while maintaining the underlying dependence on spike timing.
Does EE STDP combined with a specific form of inhibitory plasticity to do this job? (similar to the numerical work in Litwin-Kumar Doiron 2014 which is mentioned in the discussion)?
This is a good point, and we have now addressed it in the following paragraph that we added to the Discussion: As we mentioned, our theory can be used to show that small modifications to weight updates can stabilize different STDP rules. The question remains whether Hebbian EE plasticity can be stabilized through an interaction with STDP rules at different synapses? For instance, Litwin-Kumar and Doiron used a triplet voltage STDP rule that was stabilized by hard constraints and weight normalization to produce network assemblies [35]. This rule by itself lead to stable but pathological behavior, and they introduced iSTDP to restore a balanced, asynchronous network state. While such voltage-based triplet rules are outside the scope of the present study, we could use extensions of the mean-field theory to describe the impact of second and higher order moments on the evolution of weights, and network dynamics [53]. Our theory suggests that the classical pairwise Hebbian STDP cannot be stabilized by other STDP rules such as iSTDP.
Kohonen's rule that the authors investigate is stable-but I am not aware of this rule having an electrophysiological basis. The authors do not cite literature to that effect, and so I do not agree that the current results support their statement (l. 440): "Therefore, this framework provides mathematical tractability and biological relevance to a model of excitatory STDP [25], without the need for unrealistically fast inhibitory STDP."  Table 1). The mean synaptic weight is shown in black. Individual weights diverge to ±∞, destabilizing network dynamics. B: Distribution of synaptic weights right before the simulation is terminated. C: The mean synaptic input diverges as individual weights diverge. D: Positive feedback loops due to unconstrained Hebbian STDP lead to uncontrolled activity in the network. The simulation is terminated when input and rates grow with no bounds. E: Same as A-D, but for Classical Hebbian plasticity with hard constraints at 0 and 50. Mean synaptic weight (in black) remains unchanged. Network remains balanced and stable throughout the simulation. F: Same as A-D, but for weight dependent Hebbian STDP (See Table  1). The network maintains a stable, balanced state.
We agree with this comment. Our point here was to highlight the mathematical tractability of a nontrivial excitatory STDP rule such as Kohonen's. We are also not aware of electrophysiological evidence of such a rule, although it has been argued that it is biolog-ically plausible [54]. However, we still highlight generally desirable features of this rule such as stability and competition. We have now added the following statement to the revised manuscript.
108 "only pairwise interactions are meaningful". Please clarify what "meaningful" means here, since a few lines afterwards the authors say they will consider more general examples in the future (115 "the BCM rule... depend on interactions beyond second order ...").
Another class of plasticity rules that cannot immediately be written in terms of the current theory are calcium based rules (Shouval et al 2002, Graupner Brunel 2012 We agree that this statement was unclear, and have rephrased it. What we meant is that in our approach only pairwise interactions are accounted for. As the reviewer points out, higher order interactions can be essential in determining changes in synaptic weights. We have now rephrased this statement in the new version: Higher order interactions are at the heart of triplet rules [59][60][61], and other types of interactions may also be important, e.g. for calcium-based update rules [62,63]. For simplicity we here focus on pairwise interactions between spikes and eligibility traces, and leave extensions to more complex rules for future work.
124 Timescale assumption. The separation of timescales assumption should be introduced more clearly, on a number of levels: 1. Some readers may appreciate being given approximate values for the different timescales in the theory, perhaps in relation to the specific examples in Table 1.
We have now extended our analysis to better illustrate the importance of the separation of timescales in different STDP rules. This is now shown in the section entitled "Separation of timescales" in the Appendix. We also agree that it will be useful for readers to see the actual values of different timescales and we have now included this in the text.
We used 1/η ab = 10 5 ms, and τ STDP = 200ms, with correlation time window width T win = 250ms. Our derivations require τ STDP ∆T 1/η ab , however an exact numerical value for ∆T is neither used nor needed. The precise agreement can depend on the example. We used values for which the analytical expression agreed with numerical simulations qualitatively and quantitatively. As the separation of timescales is relaxed, both quantitative and qualitative agreements break down as the network becomes unstable due to increased fluctuations in the weights (See new Supplementary Figure 2). We now address this point in the Discussion.
The present theory relies on a separation of timescales between spiking dynamics and weight changes. Such timescale separation is supported by a number of experiments [30-32, 58, 64]. We show in the Appendix (see "S1 Appendix: What happens when timescales are not separated?") that reducing this timescale separation, and increasing weight updates leads to a breakdown of the theory, and can result in network instability.
In mammalian brains, timescales of weight changes may not always be separated from rate and correlation timescales. The size and timescale of weight updates is likely to depend on many factors that can modulated STDP, such as spiking patterns, synapses type, brain area, network state, neuromodulation, and others. Separation of timescales may not be pronounced in certain non-cortical areas, such as the hippocampus, which can be rapidly modified [58]. For example, Petersen et al., 1998 and Froemke et al., 2006 found significant changes in putative synaptic weights over short timescales in hippocampal CA1/CA3 slices [65] and in visual cortical slices subject to multispike pre-and post-synaptic bursts [66], respectively. However, it is possible that the rate of change of synaptic weights may be overestimated in vitro [58].
3. some computational studies indicate that the separation of timescales assumption must be violated for network stability, at least in relation to homeostatic plasticity (Zenke Gerstner Ganguli 2017). Does the theory provide a path for getting around these issues?
This is a good point to raise that needed clarification. We have added a paragraph addressing this in the Discussion.
How is our separation of timescales assumption affected when rapid compensatory processes are needed for homeostasis, given that experiments show that homeostasis is a process that is even slower than the timescale of STDP? Experimental evidence suggests that homeostatic processes can take hours or days [57,58,64]. On the other hand, theoretical models show that synaptic plasticity can be unstable in the presence of such slow homeostasis, and needs to be coupled with rapid compensatory processes such as inhibitory STDP [57,58]. The separation of timescales in our theory still puts synaptic dynamics on the "fast" side of the spectrum, as it separates network dynamics that occur over milliseconds from weight dynamics that take place over seconds or minutes. Hence, the assumption of timescale separation is still valid in our implementation of homeostatic inhibitory plasticity. A: Evolution of mean excitatory synaptic weight subject to weight-dependent Hebbian STDP in a balanced network. Fast weight timescales induce large fluctuations and deviations from the fixed point. Inset: When η = 1 weight updates are so large that the network becomes unstable. B: Same as A:, but for excitatory weights evolving according to Kohonen's rule. Convergence to the fixed point is very slow when η < 0.01, and the approximation breaks down for η = 1. C: Same as A:, but only inhibitory weights are plastic and follow inhibitory STDP updates. Deviations from the theory appear when timescales are not well separated (η ≥ 1/1000). In all panels, dashed lines represent theory (Eqs. (3,5,12)), and solid lines represent numerical results.
154 Intuitively I expect Neumann's approximation to break down for sparse matrices because the inverse W −1 0 is not well defined when the rank is not full. Is the 10% connectivity sufficiently large to avoid this issue? If yes, how small can the connection probability be before the approximation (and consequently the theoretical analysis) breaks down? Guzman, Schlogl, et al 2016 showed 1% connectivity in hippocampal region CA3-a region where recurrent connectivity and plasticity are thought to play important roles. The authors should comment on whether they expect their theory to be applicable in such a scenario.
This is a good point. It is crucial to keep in mind the difference between J (full connectivity matrix of size N × N , sparse) and W (mean-field connectivity matrix of size 2 × 2, full). We only assumed that W is invertible, which is generically true because W is full.
Previous work has considered mean-field theory of balance with multiple populations. In such cases W larger than 2 × 2, but still not sparse. In such cases it was shown that a singular connectivity matrix, W , can break balance. In the cases we consider, our assumption that W is non-singular is generic. We clarified these points with the following paragraph: The 2 × 2 mean-field connectivity matrix, W 0 , must be non-singular for the balanced state to exist and for Neumann's approximation to hold [13]. While the non-singularity of W 0 is a non-restrictive condition for two neural populations, W 0 can become singular in some models with several neural subpopulations [47,52].
158 ∆W can be defined such that the perturbation is done separately on EE, EI, IE, II connections. Does doing so within the current theory lead to the same results as Mongillo, Rumpel, Loewenstein 2018?
Mongillo et al., 2018 describe the remarkable stability of firing rate patterns to perturbations in excitatory, but not inhibitory synaptic weights. To understand the reason behind this observation they developed a numerical and theoretical approach to study the impact of randomly perturbing EE, EI, IE, and II weights. The fact that these perturbations are randomly generated, and not activity driven, is an important difference between their study and ours: While Mongillo et al. were agnostic about the particular rules that change the weights, we focused on modeling particular rules that drive changes in synaptic weights.
This difference in the approaches also makes it difficult to compare our results directly. For instance, the effect discussed by Mongillo et al. would likely be even stronger under the iSTDP rule that we discuss in the manuscript: The sensitivity of activity patterns to a change (rewiring in their case) of I-to-I and I-to-E, but not E-to-E connections is a result of the wider distribution of the I firing rates. Since the effect of iSTDP in our case is to reduce the variance in E firing rates further, we would expect such difference in sensitivity to be amplified in our network. On the other hand, applying the iSTDP rule to II synapses would reduce the distribution of I firing rates, as we discuss in the manuscript. We expect that this in turn would increase the sensitivity of firing rate patterns to perturbations in II synapses. It is possible that other plasticity rules, and their combination increase the variance of firing rates in either or both population resulting in different changes in sensitivity to perturbations. The answer to this question therefore depends on the learning rules, and their effect on the patterns of activity in the network.
We now mention this point in the Discussion where we write: We have shown that different plasticity rules can result in distinct firing rate distributions in different subpopulations. As shown by Mongillo et al. this can result in an increase or decrease in sensitivity of activity patterns and memories to perturbation of different synapse classes [67].
170 The writing in the first results paragraph should be improved to clarify this is an outline of the rest of the paper. Currently the different sections and applications of the theory come as a surprise. Otherwise I like the structure of the paper where introduction is mixed with development of the theory.
This is a good point. We have now rewritten this paragraph to make it clear that this is just an overview of the remainder of the paper.
In the following sections, we use the theory described in the Methods to show how synaptic weights and correlated activity co-evolve in tightly balanced networks. We first apply our theory to a balanced network with plastic excitatory synapses that evolve according to Kohonen's STDP rule [50]. The theory predicts that the network is in a stable, balanced state, and provides asymptotic values for rates, spike count covariances, and synaptic weights. We also show that correlated activity can have a moderate impact on the dynamics of synaptic weights. Next we show that the theory does not accurately predict weight dynamics for certain STDP rules for which heterogeneity in rates can lead to competition between synaptic weights. The resulting correlations between rates and weights are inconsistent with classical balanced network theory. We show how our framework can be extended to describe such networks. Finally, we demonstrate that our theory can describe the robust response of iSTDP to targeted optogenetic input to subsets of neurons.

How does the stability depend on beta?
This is an important question that we have now addressed in the Appendix. Under Kohonen's rule, the dynamics undergo a saddle-node bifurcation as β increases. We have now included these results in "S1 Appendix: Saddle-node bifurcation of excitatory weights in Kohonen's STDP rule", where we show the presence of a saddle-node bifurcation theoretically and confirmed it with simulations.
Saddle-node bifurcation of excitatory weights in Kohonen's STDP rule. Assume that EE weights evolve according to Kohonen's STDP rule [50,56]: Averaging gives the mean-field dynamics of the excitatory weights: In the asynchronous state, the fixed point of this equation is J * ee = βτ STDP r e .
How does the stability of the synaptic weights depend on β? Our meanfield theory shows that the stable network state disappears in a saddle-node bifurcation (Supporting Fig. 5 A). For the chosen set of parameters (See Supporting Table 1), the theory predicts a bifurcation at β ≈ 8 (Supporting Fig.  5 A). We confirmed the presence of a saddle-node bifurcation is confirmed in numerical simulations (Supporting Fig. 5 B). However, the bifurcation occurs at a lower value of β due to divergence of individual synapses (Supporting Fig.  5 C-F). This is not unexpected, as fluctuations can destabilize a system close to a saddle-node bifurcation by pushing it across the unstable fixed point.  Yes, they were only E cells. We have now included rasters of the full population color coding I cells as red and E cells as blue. The authors should indicate and discuss whether their simulations capture the corresponding experimentally observed properties of plastic neuronal networks beyond the averages.
We have now included more statistics and shown the dynamics of different realizations of networks under plasticity rules such as Kohonen's, iSTDP, and weight-dependent Hebbian in Supplementary Figs. 6-8. For the case of Kohonen's rule, we have shown this in "S1 Appendix: Statistics and dynamics of balanced networks undergoing pairwise STDP rules". The EE weights in Kohonen's rule follow a unimodal distribution. Distributions of rates are unimodal and asymmetric (Supporting Fig. 7 A-B) as seen in experiments. The distribution of covariances and that of correlations are unimodal centered close to zero in the asynchronous state (Supporting Fig. 7 E,F). We have added the following section to "S1 Appendix: Statistics and dynamics of balanced networks under pairwise STDP rules", and have included distributions of rates, covariances, correlations, and weights for networks undergoing weight-dependent Hebbian STDP, Kohonen's STDP, and iSTDP: Statistics and dynamics of balanced networks under pairwise STDP rules. Our theoretical results describe the evolution of averages (See Methods in main text). However, it is also important to understand changes in higher order statistics of individual, and collective activity under plasticity. Although we do not go beyond a mean-field theory, we provide the results of corresponding numerical simulations in three example balanced networks where: (1) excitatory weights were subject to weight-dependent Hebbian STDP; (2) EE weights changed according to Kohonen's rule; and (3) inhibitory weights changed according to an iSTDP rule (See Table 1 in main text for equations describing each rule).
Under weight-dependent Hebbian STDP individual synaptic weights converge to a unimodal distribution. This is in agreement with experimental results [68] (Supporting Fig. 6 A,B). Balance is preserved throughout the simulation (Supporting Fig. 6 C). Distributions of rates are unimodal and skewed to the right (Supporting Fig. 6 D and inset), also agreeing with experimental findings [68]. Spike count covariances and correlations are distributed around zero in the asynchronous state (Supporting Fig. 6 E,F).
EE weights under Kohonen's rule also converge to a unimodal distribution (Supporting Fig. 7 A,B). Again, the network is stable and balance is preserved (Supporting Fig. 7 C). The distributions of rates is again unimodal and skewed to the right (Supporting Fig. 7 D), and covariances and correlations are close to zero (Supporting Fig. 7 E,F). Lastly, EI weights subject to an inhibitory plasticity rule also follow a unimodal distribution (Supporting Fig. 8 A,B), distributions of I rates are unimodal and skewed to the right, and covariances and correlations are centered at zero in the balanced, asynchronous state (Supporting Fig. 8 C-F). However, the E rates are concentrated at the target rate (Supporting Fig. 8 D -inset), in disagreement with experimental findings [68]. In short, all three STDP rules show realistic distributions of weights, rates, and covariances, except for the iSTDP's distribution of E rates. B: EE synaptic weights follow a unimodal distribution at steady state. C: Evolution of synaptic currents showing that balance is achieved and maintained. Color scheme as in Fig. 1 in text: blue denotes mean excitatory input, and red the mean inhibitory input. The total mean input is shown in black. D: Steady state distribution of excitatory and inhibitory (inset) firing rates. In equilibrium most neurons fire at low rates, and a few neurons fire at high rates. E: Distribution of spike count covariances. The mean is close to zero indicating an asynchronous state. F: Same as E:, but for spike count correlations, again showing a mean that is close to 0. Balance is achieved and maintained. D: Equilibrium distribution of excitatory and inhibitory (inset) firing rates. Most neurons fire at low rates, and a few neurons fire at high rates. E: Distribution of spike count covariances shows that the mean is close to zero (asynchronous state). F: Same as E:, but for spike count correlations. The mean is close to zero, suggesting that the network is in an asynchronous state. Color scheme as in Supporting Fig. 6.
We believe there might be some confusion here between τ ST DP and 1/η ab . On one hand, τ ST DP is the time constant of eligibility traces, which determines the time window over which the eligibility traces, x a j (t), are affected by past spikes. It does not substantially affect the slope of dJ ab /dt. On the other hand, η ab is the learning rate of the synaptic weights, and hence 1/η ab is the time constant of STDP. The learning rules have the form dJ ab /dt = η ab · · · . It then follows that the slope of dJ ab /dt is proportional to 1/η ab .
There is a very strong slowing down as a function of input correlations (panel D). Does this mean that if rates are small, we should expect synaptic weights to always be away from steady state for correlated inputs? If not, can an experimental prediction be extracted from this or will the transient timescale always be swamped by firing rates? the rate contribution typically dominates the correlation term (even at the largest value of input correlations c x = 1), whether it be in determining the fixed point or modulating speed of convergence. Even when rates are small (low values of j ee in Fig. 2 D), we see that the rate term is still orders of magnitude larger than the correlation term. Thus our theory predicts that, generally, changes in synaptic weights in plasticity-inducing protocols (e.g., optogenetic stimulation) will be mainly driven by modulation in firing rates, rather than in correlations.
We have modified text to more clearly present these findings.
How large is the impact of correlations in plastic balanced networks more generally? To address this question, we assumed that only pairwise interactions affect EE synapses, as first order interactions depend only on rates after averaging. We thus set B α,β ≡ 1, and all other coefficients to zero in Eq. (10). While the network does not reach a stable state under this arbitrary plasticity rule, it allows us to estimate the contribution of rates and covariances to the evolution of synaptic weights. Here B α,β can have any nonzero value, since it scales both the rate and covariance terms. Under these conditions, our theory predicts that the rate term is at least an order of magnitude larger than the correlation term (even when rates themselves are small, i.e., when j ee is small), and so correlations only have a low impact on the dynamics of synaptic weights (Fig. 2 D). Therefore, our theory predicts that, in general, changes in synaptic weights will largely be driven by changes in firing rate patterns, rather than changes in pairwise correlations.

Fig. 2 C: is red E or I?
We have clarified this now in the legend of Figure 1. Unless otherwise stated, the color scheme remains the same throughout the paper.
Unless otherwise stated, colors carry the same meaning in all figures. Thank you for the suggestion. Yes, we have now included this in Fig. 3. See also our previous answer to the comment about convergence rates.  257: The authors should cite Vegue Roxin 2019 here. They did take some steps in extending the MFT to cases with heterogeneous connectivity which leads to changes in the firing rate distribution (which are predictable given the weights).
We have now added the reference and the following text: The failure is due to correlations between weights and pre-synaptic rates which are typically ignored [13,14,23,[46][47][48], but can cause the mean-field description of network dynamics to become inaccurate. This is similar to the breakdown of balanced state theory in the presence of correlations between in-and out-degrees discussed by Vegué and Roxin, 2019 [69].
304: From a theoretical point of view, it is nice that introduction of II plasticity (of the same form as EI) eliminates the heterogeneity of inhibitory firing rates. However, experimental data (e.g. Xue et al 2014, Chiu ,... Higley 2018) suggests that different inhibitory subtypes are subject to different inhibitory plasticity rules. I think it's fine to include this "trick" in a theoretical paper but I would at least mention the fact that inhibitory plasticity rules themselves are heterogeneous so in a real network we should not expect inhibitory firing rates to become similar due to homogeneity of the inhibitory plasticity rule.
This is an excellent point. Note that while one should not expect inhibitory firing rates to become uniform across subtypes, our mean-field approach can be extended to account for multiple subtypes, in which case we would only need (approximately) uniform firing rates within each subtype, or at least uncorrelated weights and rates within each subtype. But even that is probably not biologically very realistic. We have added the following text to the manuscript: This model of inhibitory plasticity is likely a large oversimplification. Synapses of different interneuron subtypes are likely subject to different plasticity rules operating on different timescales [17,70], and would therefore not lead to uniform inhibitory firing rates. The mean-field theory we presented here can be extended to account for multiple inhibitory subtypes with different plasticity rules.

331: missing citation
Thank your pointing that out, we have now added the missing citation.
Although classical mean-field theory produced singular solutions, Ebsch et al.
showed that the theory can be extended, and that a non-classical balanced state is realized: Balance at the level of population averages (E and I ) is maintained, while balance at the level of the three subpopulations is broken [52].
416: not sure that incapacity is a word. Consider using inability.
We agree and have modified the text now.
Partial stimulation of a population of E neurons has been shown to break balance due to the inability of the network in cancelling inputs when weights are static [52].

416: typo and missing ref. "I Ebsch ..."
Thank you for catching this. We have now added the citation and fixed the typo.
Ebsch et al. showed how classical balanced network theory can be modified to account for effects of input perturbations that break the classical balanced state [52]. Vogels et al. [41] (in addition to subsequent studies [35,42,[71][72][73][74]) showed empirically using simulations that inhibitory iSTDP can restore balance. We here provide a theoretical framework that describes the evolution of rates and weights before, during, and after a perturbation that breaks balance.

426:
The paragraph citing Ocker [57] should also cite Ravid-Tannenbaum Burak 2017 Thank you. We have now included the missing reference.
Ocker et al. also studied the dynamics of mean synaptic weights in recurrent networks and derived equations for mean synaptic weights in terms of rates and covariances [44,75].
Appendix 1, Generating correlated spike trains. I was not aware of this method. Its explanation is clear, but I lack intuition whether/why it is favorable over an alternative which seems simpler: 1. Generating a Gaussian process with the exact desired correlations/covariances 2. Computing the Gaussian CDF of the random number for each neuron at each time. The method described by the reviewer is very powerful, allowing one to specify pairwise covariances for the generated spike trains. While this method is computationally feasible for hundreds of spike trains (as demonstrated by Macke, et al. [76]), it becomes computationally intractable at the numbers required for our study. For example, the simulation in Figure 3 required the generation of 10 3 spike trains over 3 × 10 7 time bins. The covariance matrix for the Gaussian process would contain 10 6 elements, so evaluation of the Gaussian CDF would be computationally expensive.
The method we used is much simpler, but more constrained. It produces spike trains with homogeneous pairwise correlations (all pairs have the same covariance), but it is much more computationally efficient for generating a large number of spike trains. For example, the generation of spike trains for Figure 3 only took 24 seconds. Appendix 3, I think the authors should expand Appendix 3 to address my major concern 1: Are there minimal modifications of the classical EE STDP rule which imply a non-trivial solution to the fixed point equation (Eq. 2 in Appendix 3)?
We have addressed this in the previous comments. Please see our answers to major concern 1, as well as extensions of S1 Appendix, including the added figures.

Referee B
The authors developed a self-consistent theory of rates, covariance and synaptic weights in balanced recurrent network that receives correlated external input. The theory can help to understand coexistence of correlated activity and dynamic excitatory-inhibitory balance. Furthermore, they found a plasticity rule for inhibitory-to-excitatory connections which have interesting property of dependence on initial condition. Finally, they showed that their theory could be applied to "out of balance" scenario of optogenetic input. Although, these are interesting and important claims we have some reservation about the amount of proof they provide for such a general claims and consequently about biological relevance of the study. First claim is that authors developed a general theory of plasticity in balanced network which can address network dynamics with different plasticity rules and describe weight evolution with correlated activity. We have some concerns about generality, novelty and predictive power of the theory. The theory presented is far from being a general theory of plasticity in balance networks. The claim hold and have be demonstrated by authors fairly well only in the case of tightly balanced large neural network (with asymptotic 1/sqrt(N) scaling) for the class of pairwise STDP rules (ignoring other plasticity types i.e. bcm, metaplasticity, intrinsic plasticity and structural plasticity). This STDP class authors chose is one with multiplicative dependence on synaptic weights which naturally lead to "nice" fixed points and for which mean-field description works very well.
We agree with this point. We have rewritten the manuscript to clarify the limitations of the theory we develop. We believe that a similar theory can be developed for other types of balanced networks as long as a mean-field description holds, and timescales are separated. For learning rules that involve higher order interactions this will require extensions of the mean field theory.
Higher order interactions are at the heart of triplet rules [59][60][61]77], and other types of interactions may also be important, e.g. for calcium-based update rules [62,63]. For simplicity we here focus on pairwise interactions between spikes and eligibility traces, and leave extensions to more complex rules for future work.
We have also added the following text to the Discussion: Our theory assumes a large network, N 1, and weights that scale as 1/ √ N resulting in "tight balance". The theory is less accurate in smaller networks. However, our simulations show that the theory provides good approximations for biologically realistic network sizes and weights.
The main theoretical contribution of the paper is procedure to put together and solve self-consistently equations for rate, covariance and weights. This theory seems very computational efficient and simple, but it is not the only one which can threat this problem. Equations 1 and 2 for rate and covariance with tight balance condition and large networks are already derived earlier (partly by the same authors Baker et al 2019 "Correlated states in balanced neuronal networks"), and there is even a statistical field theory developments in this area e.g. Ocker et al 2016 paper "Linking structure and activity in nonlinear spiking networks". The equation for individual weight evolution with input correlation were derived already by Gutig et al 2003 in paper "Learning Input Correlations through Nonlinear Temporally Asymmetric Hebbian Plasticity". The self-consistent weight evolution was done by Ocker et al 2015 in the paper "Self-Organization of Microcircuits in Networks of Spiking Neurons with Plastic Synapses", for the case of pairwise STDP rule which allows for weight heterogeneity and structure in connectivity (they chose a rule which is very sensitive to correlation in inputs). There is also work with higher-order SPDP rule, i.e. Monangie et al 2020 paper "Autonomous emergence of connectivity assemblies via spike triplet interactions".
Thank you for the comment. We have rewritten our introduction to better place our work in the context of balanced networks and spike timing dependent plasticity. Please see the new introduction text at the beginning of the document. We have also included the missing citation to Gutig et al., 2003 and Montangie et al., 2020.
To my best of my knowledge this is the first theoretical study of effect of external input correlation on tightly-balance balanced plastic networks. Similar theoretical attempts have been done before for the cases of external and self-generated correlations, e.g. Gilson et al 2009 "Emergence of network structure due to spike-timing-dependent plasticity in recurrent neuronal networks. I. Input selectivity-strengthening correlated input pathways" and Gilson 2010 "STDP in recurrent neuronal networks". Above studies did not have condition for tight balance and this work is based on this assumption. It is not clear if prescribing a balance a priory could bring clarity to the issue, because temporal dynamics of correlation can play a big role in unstable dynamical states, e.g. trough large Fano factor. This is a good point. Our theory relies on the assumption that the system is in quasiequilibrium. It is likely that in some situations weight changes are driven by events that place the network out of equilibrium. On the other hand, we do point to evidence that, for instance, inhibitory plasticity may act over longer timescales to keep the network stable. There is also considerable evidence that many cortical networks are in a balanced state. As we stated above, we believe that our theory provides a basis for describing plastic balanced networks when a separation of timescale holds. These are the same conditions that underlie much of previous theoretical work. We believe that this is the first time this approach has been extended to balanced networks, as the reviewer notes.
We now point to the limitations of the theory in the following text added to the Discussion: It is likely that some changes in the weights are driven by transient events that place the network out of equilibrium. If plasticity is mainly driven by such fast changes, the timescale separation central to our theory does not hold. As we have shown in S1 Appendix, the predictions of the theory can break down in such cases.
One of the most important issues in understanding neural dynamics is establishing conditions under which the network remains active, yet stable as synaptic weights change. The theory we developed can help us address these questions, but it does have limitations. Since we used a mean-field approach, we can only capture first moments. While mean weight stability may not imply stable network dynamics (consider the case when weight variance diverges in Classical Hebbian STDP in S1 Appendix), instability in the mean weights does imply that the network is also unstable.
Second claim is that that correlations in the input mildly, but significantly affect the evolution of synaptic weights. This claim is well demonstrated, but it is a hardly surprising, given the type of the multiplicative pairwise rule authors use.
We agree that in retrospect this may not be surprising, and it is a general consequence of these learning rules. We now note that Ocker et al., 2015, andGraupner et al., 2016 have made similar observations. This conclusion is somewhat restricted to the types of rules that we used. For instance, as also shown by Ocket, et al., if we reduce the weight dependent terms by choosing a nearly anti-symmetric STDP function, correlations can drive plasticity [44]. Our theory can be applied to this case. We can also extend the main ideas to other types of plasticity rules as long as separation of timescales holds, although this may require describing the higher order statistics of the network. We now make a note of both of these points.
Higher order interactions are at the heart of triplet rules [59][60][61]77], and other types of interactions may also be important, e.g. for calcium-based update rules [62,63]. For simplicity we here focus on pairwise interactions between spikes and eligibility traces, and leave extensions to more complex rules for future work.
We have also commented about this in the Discussion: Ocker et al. also studied the dynamics of mean synaptic weights in recurrent networks and derived equations for mean synaptic weights in terms of rates and covariances [44,75]. They also found that the effect of correlations in synaptic weight dynamics is mild, however, their approach relied on linear response theory to obtain expressions of spike train covariances. Such a description is valid when neurons are driven by strong white noise in addition to weaker synaptic input from the network. Using networks of two neurons with varying firing patterns, Graupner et al. also found that changes in synaptic weights were dominated by firing rates, with correlations playing a secondary role [51].
Third claim is that their inhibitory-to-excitatory plasticity induce correlations between activity and synaptic weights, such that final weight depends on initial synaptic weight. This is certainly interesting result, and especially because it could explain one role of inhibitory-to-inhibitory plasticity. That being said the rule used is very peculiar. The RHS of the ODE is scaled with initial position and authors did not fully clarify what is the consequence of this choice? What is biological justification? It is very confusing for reader that initial position is a parameter of an ODE. We would be interested to see a comparison of this rule to the original Vogels at al 2011 "Inhibitory plasticity balances excitation and inhibition in sensory pathways and memory networks" to properly understand the scaling choice. More importantly temporal evolution of the weights is not shown in this case.
We understand why this can be confusing, especially since we show that weights in iSTDP depend on the initial condition later on. The only purpose of this normalization constant was to match the dimensions on both sides of the equation. The same could have been accomplished by absorbing the constant in the learning rate, η ab .
When we implemented the rule we used a single constant in all simulations, rather than an initial condition. Thus the equation we provided was in error and was not consistent with the simulations we presented, and the actual simulations are closer to what the reviewers suggest above. This holdover from a previous analysis has now been corrected in the text.
The normalization does not depend critically on the specific value of the constant, but it is required to be of the same order (and sign) as synaptic weights J ei jk , that is O(1/ √ N ) and negative. To clarify this point have now renamed the normalizing constant J norm , and addressed this in "Results: Balanced Network with Inhibitory Plasticity".
To illustrate this, we consider a balanced network subject to homeostatic plasticity [41]. This type of plasticity has been shown to stabilize the asynchronous balanced state and conjectured to play a role in the maintenance of memories [35,41,72]. Following [41] we assume that EI weights evolve according to: where α e is a constant that determines the target firing rates of E cells and J norm ∼ O(1/ √ N ) is a normalization constant. Note that J norm < 0 so the fraction in Eq. (4) is positive assuming J ei jk < 0. In a departure from the rule originally proposed by Vogels et al. [41], we chose to multiply the time derivative by the current weight value. This modification creates an unstable fixed point at zero, prevents EI weights from changing signs, and keeps the analysis mathematically tractable (See S1 Appendix for details). An alternative way to prevent weights from changing sign would be to place a hard bound at zero, but this would create a discontinuity in the vector field of J ei , complicating the analysis.
Final claim is that this framework could be used to treat case of optogenetic input to subpopulation. We think that this claim is shown well, but it is the question if the study brings more clarity to the field than original Vogels 2011 study.
We agree this may have not been very clear. Our main result here is not that iSTDP can help maintain balance. As the reviewer points out, this was already shown by Vogels et al., 2011. Rather we provide a theory that describes the evolution of weights and rates before, during, and after stimulation. While Vogels et al. demonstrated and explained how iSTDP can restore tight balance, as far as we know, ours is the first analytical description of the evolution of the actual values of mean firing rates and synaptic weights during and after a balance-breaking perturbation. In the manuscript, this serves as a demonstrative example of how our theoretical framework can be applied, using an example that arises naturally in previous work.
To explain this point more clearly we have added the following to the Discussion: Ebsch et al. showed how classical balanced network theory can be modified to account for effects of input perturbations that break the classical balanced state [52]. Vogels et al. [41] (in addition to subsequent studies [35,42,[71][72][73][74]) showed empirically using simulations that iSTDP can restore balance. We here provide a theoretical framework that describes the evolution of rates and weights before, during, and after a perturbation that breaks balance.

Referee C
The manuscript addresses the dynamics of plasticity under STDP in randomly connected neural networks in the balanced regime. The manuscript provides some analytical insights into this topic, backed up by numerical simulations. The topic is important, and results seem substantial and interesting. Nevertheless, there are some significant weaknesses in the analysis and presentation, which I recommend to address before the manuscript is considered again for publication. Since the theory is derived at the level of mean population activities, there is a hidden assumption that synaptic weights evolve uniformly. Is this true empirically in the simulations? How variable are the weights following plasticity?
Even though mean-field approximations are accurate under uniformity, they can also be accurate when weights and rates are not uniform. Our theory only requires that J ab jk r b k j,k ≈ J ab jk j,k r b k k where · j,k is the average over indices j and k. We provide an example of when and how this assumption breaks, and we describe a semi-analytic approach to restore the mean-field theory (Fig. 3).
We have now included supplementary figures that show the uniform evolution of randomly sampled individual weights in Kohonen's, weight-dependent Hebbian, and inhibitory STDP rules in (Supporting Figs. 6 A,7 A,8 A). A longer discussion of this new material is provided in the response to comments by Reviewer A.
Related to this point, there is also an implicit assumption that weights evolve independently of other quantities. One example is the correlation between presynaptic firing rate and the synaptic weight, which is discussed in the section on inhibitory plasticity. But there are other possible correlations that could develop between synaptic strength and, for example: the in-degree and out-degree of involved neurons; under some rules competition between synapses belonging to the same neuron could also break some of the assumptions of the mean field theory.
We agree with this comment. In all of our simulations, the in-degree and out-degree are almost equal, because p ab ≡ 0.1 for all a, b = e, i, x and N was large. Different extensions of the theory would have to account for other types correlations between different weights, degrees, and/or activity. We have now cited the work by Vegué and Roxin 2019, and commented on it in the Discussion: In plastic networks, correlations between weights and other features such as in-degrees, or out-degrees can emerge [69]. We have shown how the theory can capture the case in which synaptic weights and pre-synaptic rates are correlated. While we were not able to find analytical expressions for these correlations, we showed that a second-order correction is sufficient to explain the observed dynamics. Eventually, the mean-field theory would need to be extended to account for higher order network motifs and their potential correlations with synaptic weights and firing rates. This might be possible by extending our approach, but we leave these extensions for future work.
A full theory of the plasticity dynamics should address these questions, but I accept that there is value in deriving the predictions of a mean field theory that ignores such correlations, while comparing the predictions with numerical results. The presentation should acknowledge these limitations more explicitly. The appearance of a discussion of correlations between firing rates and synapses only in the section on inhibitory plasticity is unsatisfactory.
We agree with this point and have added it to the Discussion.
In plastic networks, correlations between weights and other features such as in-degrees, or out-degrees can emerge [69]. We have shown that the theory applies to the case in which synaptic weights and pre-synaptic rates are correlated. While we were not able to find analytical expressions for these correlations, we showed that a second-order correction is sufficient to explain the observed dynamics. Eventually, the mean-field theory would need to be extended to account for higher order network motifs and their potential correlations with synaptic weights and firing rates.
We have also added the following text to the Methods: Higher order interactions are at the heart of triplet rules [59][60][61]77], and other types of interactions may also be important, e.g. for calcium-based update rules [62,63]. For simplicity we here focus on pairwise interactions between spikes and eligibility traces, and leave extensions to more complex rules for future work.
Is there a way to understand why correlations did not play an important role in earlier sections?
This is an important question, and we have added a paragraph to the Discussion to address this more clearly.
A natural question that arises is why do correlations between weights and pre-synaptic rates only seem to play a role in iSTDP? In the examples of excitatory STDP we analyzed (Kohonen's rule and weight-dependent Hebbian rule), weights at equilibrium are determined by other parameters (Weightdependent Hebbian rule) or rates (Kohonen's rule). Therefore weights are updated until those steady state values are achieved, yielding values independent of initial conditions. On the other hand, in the case of the inhibitory plasticity rule, inhibitory weights at equilibrium are determined by the firing rates alone. Since the firing rate vectors are lower-dimensional than the weight matrices, the equilibrium solution does not fully determine the weight matrices. This is shown in Fig 3 C inset, where different distributions of weights can result in the same equilibrium firing rate when weights and pre-synaptic rates become correlated.
While the theory is derived for dynamics of the synaptic weights, the numerical analysis does not demonstrate that the dynamics indeed follow equation 6. Instead, simulations are used only to demonstrate agreement with the predictions on the steady state. Do the dynamics follow the prediction?
We agree. The system of equations can be solved iteratively over time as well, provided the timescales are separated. We have now included this in S1 Appendix: "Transient dynamics of synaptic weights".
Transient dynamics of synaptic weights. In much of this work we focused on steady state behavior. However, the system of equations (Eqs. (14,21,71)) relating rates, covariances, and weights can also be solved iteratively over time, I would also like to see some numerical support for statements that the network remains in an asynchronous state following plasticity.
The figures are included earlier in the letter, and we repeat the new text here for convenience: Statistics and dynamics of balanced networks under pairwise STDP rules. Our theoretical results describe the evolution of averages (See Methods in main text). However, it is also important to understand changes in higher order statistics of individual, and collective activity under plasticity. Although we do not go beyond a mean-field theory, we provide the results of corresponding numerical simulations in three example balanced networks where: (1) excitatory weights were subject to weight-dependent Hebbian STDP; (2) EE weights changed according to Kohonen's rule; and (3) inhibitory weights changed according to an iSTDP rule (See Table 1 in main text for equations describing each rule).
Under weight-dependent Hebbian STDP individual synaptic weights converge to a unimodal distribution. This is in agreement with experimental results [68] (Supporting Fig. 6 A,B). Balance is preserved throughout the simulation (Supporting Fig. 6 C). Distributions of rates are unimodal and skewed to the right (Supporting Fig. 6 D and inset), also agreeing with experimental findings [68]. Spike count covariances and correlations are distributed around zero in the asynchronous state (Supporting Fig. 6 E,F).
EE weights under Kohonen's rule also converge to a unimodal distribution (Supporting Fig. 7 A,B). Again, the network is stable and balance is preserved (Supporting Fig. 7 C). The distributions of rates is again unimodal and skewed to the right (Supporting Fig. 7 D), and covariances and correlations are close to zero (Supporting Fig. 7 E,F). Lastly, EI weights subject to an inhibitory plasticity rule also follow a unimodal distribution (Supporting Fig. 8 A,B), distributions of I rates are unimodal and skewed to the right, and covariances and correlations are centered at zero in the balanced, asynchronous state (Supporting Fig. 8 C-F). However, the E rates are concentrated at the target rate (Supporting Fig. 8 D -inset), in disagreement with experimental findings [68]. In short, all three STDP rules show realistic distributions of weights, rates, and covariances, except for the iSTDP's distribution of E rates.
Lines 359-365: What happens under the perturbation without plasticity? How does this depend on parameters? How can we see that the network is in a balanced state following plasticity (but possibly not before plasticity)?
After external input onset, the network initially exits the balanced state when some E neurons are silenced due to excess average inhibitory input (Fig. 4 D). However, while targeted cells still receive optogenetic input, inhibitory weights adjust to restore detailed balance in the network. The network will remain at or return to the balanced state so long as parameters are chosen such that W is non-singular. Thus a prolonged input will drive the network to a new balanced state, while a transient input may drive the network out of balance for a short period.
We now clarify this point in the following: We simulated a network of N = 10 4 neurons in an asynchronous state with c x 1 = c x 2 = 0. A subpopulation of 4000 E cells receives transient, correlated input. Solving Eqs. (27,28) predicts that inhibitory plasticity will alter EI synaptic weights so that the firing rates of both the E expr and the E non−expr approach the target firing rate before, during, and after stimulation. Once the network reaches steady state the mean inputs to each subpopulation cancel. Thus changes in EI weights restore balance at the level of individual subpopulations or "detailed balance," consistent with previous studies [41,72]. Simulations confirm these predictions (Fig. 4 B-D).
When the input is removed, the inhibitory weights onto cells in the E expr subpopulation converge to their pre-stimulus values, returning E expr rates to the target value, and reestablishing balance (Fig. 4 B-D). Correlations remain low (O ∼ 10 −4 ) before, during, and after stimulation (Fig. 4 E,F), suggesting that at equilibrium the network is in the asynchronous state. Our theory thus describes how homeostatic inhibitory STDP increases the stability and robustness of balanced networks to perturbations by balancing inputs at a level of individual cells, maintaining balance in regimes in which non-plastic networks cannot maintain balance. We presented an example in which only one subpopulation is stimulated. However, the theory can be extended to any number of subpopulations in asynchronous or correlated balanced networks receiving a variety of transient stimulus.
In comparisons between simulations and the theory, for how long were the simulations run? What was the criterion used to identify that the weights reached steady state?
We have added these details in the captions throughout the text. We applied the Kolmogorov-Smirnov 2-sample test to determine if synaptic weights had reached equilibrium. We added text clarifying this in "Appendix: Determining when synaptic weights reach equilibrium in simulations".
Determining when synaptic weights reach equilibrium in simulations. To determine if synaptic weights had reached steady state we sampled 1000 synaptic weights from the whole population (of plastic weights) at times t = 4000sec and t = 5000sec (end of simulation). We then compared the two distributions using the Kolmogorov-Smirnov 2-sample test. We found that in the three networks simulated in Supporting Figs. 6-8, the distributions were not distinguishable, and hence had reached equilibrium by the end of the simulation (Supporting Fig. 4 A-C)  Figure 4: General framework for STDP in balanced networks describes the dynamics of networks receiving optogenetic input. A: A recurrent network of excitatory, E, and inhibitory, I, neurons is driven by an external feedforward layer X 1 of uncorrelated Poisson neurons. Neurons that express ChR2 are driven by optogenetic input, which is modeled as an extra layer of Poisson neurons denoted by X 2 . B: Evolution of mean synaptic weights over the course of the experiment. C: Evolution of mean firing rates. Inhibitory STDP maintains E rates near the target, αe 2τ STDP . D: Evolution of mean excitatory, external, inhibitory, and total currents. Balance is transiently disrupted at stimulus onset and offset, but it is quickly restored by iSTDP. E: Mean spike count correlations before, during, and after stimulation remain very weak for all pairs. F: The distribution of spike count correlations also remains nearly unchanged with weak mean correlations before, during, and after stimulation. Solid lines represent simulations, dashed lines are values obtained from theory (Eqs. (3,5,27,28)).
Line 6: it will be appropriate to cite also Shaham and Burak (2017) and Darshan et al (2017) in the context of spatially correlated activity in the balanced state.
The introduction has been rewritten to more appropriately place the work in the balanced network and synaptic plasticity literature and these references have now been added.
In the discussion, it will be appropriate to cite Ravid-Tannenbaum and Burak (2016) in the context of the relationship between covariances and STDP dynamics.
Thank you. We have now included the missing reference.
Ocker et al. also studied the dynamics of mean synaptic weights in recurrent networks and derived equations for mean synaptic weights in terms of rates and covariances [44,75].
Lines 58-59: why is the assumption not essential? How does this agree with Ref. 34, which shows that heterogeneity in the connection probabilities can destroy the balanced state?
We meant that this assumption is not essential as long as the matrix W is non-singular. Balance can be maintained or destroyed depending on several factors, and heterogenous connectivity may not necessarily disrupt it (Pyle PRE, 2016). We have modified our statement to clarify this.
For simplicity we assume that both the probabilities, p ab , and weights, j ab ∼ O(1) are constant across pairs of subpopulations. This assumption is not essential as long as the balance condition remains valid. This is equivalent to having a non-singular mean-field connectivity matrix W , which is defined next.
PRESENTATION It is inconvenient that only some of the equations are numbered. This makes it difficult to refer to the unnumbered equations when discussion the manuscript.
We have now numbered all the equations for the readers' convenience.

Equation 2
: please define Gamma precisely, not only in limits, and either derive this result or provide a clear path from the results of previous papers to this expression, possibly in an appendix.
We now point to detailed derivations of Γ already presented in Baker et al., 2019. A closed form for finite N is not tractable for spiking neural networks.
The matrix Γ has the same structure as C and represents the covariance between external inputs (See Baker et al., 2019 Appendix A for a detailed derivation of this term [23]).
The role of Twin in equation 2 is not sufficiently clear. If I understood this correctly, Twin should be introduced already when defining the spike count covariance.
We agree with this comment and now define T win earlier.
We define the mean spike count covariance matrix as: where C ab is the mean spike count covariance between neurons in populations a = e, i and b = e, i respectively, counted over time windows of size T win . Throughout all simulations and theoretical predictions, we set T win = 250ms, however the theory is flexible to other time window sizes.
Lines 73-74: the statement that the second term in Equation 2 accounts for intrinsically generated variability is a bit confusing since the elements in the diagonal are negative. It will be helpful to clarify this.
We realize how this was confusing. We now remind the reader that in the spike count covariance matrix the diagonal terms are not variances, but the mean covariances between spike trains in the same population. We have clarified this in the following added text: In Eq. (5), F a is the Fano factor of the spike counts averaged over neurons in populations a = e, i over time windows of size T win . The second term in Eq. (5) is O(1/N ) and accounts for intrinsically generated covariability [23] within excitatory or inhibitory populations (note that this term does not refer to variances, but instead to mean covariances between spike trains in the same subpopulations).
The logic of how x and S combine to generate the STDP rule is mathematically simple, yet for many readers it will be difficult to understand the structure of the rule from the equations. It will be helpful to provide a graphical representation of the STDP rules in the different cases listed in Table 1.
We agree with this comment. We have now added Supporting Fig. 1 to show the STDP window for each rule presented in Table 1.  Figure 1: STDP windows of different plasticity rules. A: Change in synaptic weights as a function of the relative timing of pre-and post-synaptic spikes in Classical Hebbian STDP (same as weight-dependent Hebbian). B: Same as A, but for inhibitory STDP. C: Same as A, but for Kohonen's rule when weights are below parameter β. D: Same as C, but for the case when weights are above β. E: Same as A, but for Oja's rule when weights are below parameter β. F: Same as E, but for the case when weights are above β.
Lines 140-148: I didn't understand why root finding is involved in the numerical solution of the differential equations. Doesn't equation 6 simply provide the required update of J, based on the values of r anc C that are obtained from equations 1-2?
We realize this was not clear before and have addressed this now.
This process can be repeated iteratively at each time step to obtain the evolution of firing rates, spike count covariances, and synaptic weights provided our assumption about the separation of timescales holds (See S1 Appendix: "Transient dynamics of synaptic weights" for more details). Figure 3: "and the system predicts a stable fixed point". Should "system" be replaced by "theory"?
Yes, thank you for pointing that out. The rate of change of EI weights as a function of the weights themselves. The contributions of the covariance (blue) is considerably smaller than the contribution of the rate (red), and the theory predicts a stable fixed point.
Thank you for catching this.
We show that inhibitory STDP, as described by Eq. (22), can restore balance in the inputs to the stimulated and unstimulated subpopulations.
The last equation in page 1 of SI1 includes parameters whose meaning is not defined.
Thank you for noticing that. We have now pointed to the specific table containing these parameters.
See Supporting Table 1 for a description and values of the parameters of the voltage dynamics. In our simulations, this equation was integrated using Forward Euler's method with step size dt = 0.1ms.
In line 26 of SI1, units for dt are missing.
Thank you for pointing that out. We now specify units: dt = 0.1ms In line 53 of SI1, in the equation rm = r x /cx, should cx be c x ? (where stands for subscript).
Yes, thank you. We have now fixed this typo: To generate N x processes, all having firing rate r x and pairwise correlation c x , we first generate a "mother" process with rate r m = r x /c x .
Lines 53-56: It is not easy to understand from the description how the daughter process spikes are obtained. I suggest to clarify that each daughter process is obtained by first taking all the spikes of the mother process, followed by a deletion of spikes with probability 1 − c x , and to mention explicitly that deletion is performed independently in each one of the daughter processes.
We agree with this comment, and have made this point more clear in the text.
To generate N x processes, all having firing rate r x and pairwise correlation c x , we first generate a "mother" process with rate r m = r x /c x . Then, to generate each of the N x "daughter" processes, we first take all the spikes in the mother process and then delete each spike independently with probability 1 − c x . In other words, each spike time from the mother process is included in each daughter process independently with probability c x . The rates of the daughter processes are r m c x = r x , as desired. Also, two daughter processes share a proportion c x of spike times (since the probability that a spike time in one process also appears in the other process is c x ).
It will be more convenient to group all the supporting text in one file with several sections.
Thank you for the suggestion. We have now combined all Supporting Information into one file.