Inhibition stabilization is a widespread property of cortical networks

Many cortical network models use recurrent coupling strong enough to require inhibition for stabilization. Yet it has been experimentally unclear whether inhibition-stabilized network (ISN) models describe cortical function well across areas and states. Here, we test several ISN predictions, including the counterintuitive (paradoxical) suppression of inhibitory firing in response to optogenetic inhibitory stimulation. We find clear evidence for ISN operation in mouse visual, somatosensory, and motor cortex. Simple two-population ISN models describe the data well and let us quantify coupling strength. Although some models predict a non-ISN to ISN transition with increasingly strong sensory stimuli, we find ISN effects without sensory stimulation and even during light anesthesia. Additionally, average paradoxical effects result only with transgenic, not viral, opsin expression in parvalbumin (PV)-positive neurons; theory and expression data show this is consistent with ISN operation. Taken together, these results show strong coupling and inhibition stabilization are common features of the cortex.


Introduction
Extensive recurrent connectivity between nearby neurons is a ubiquitous feature of the cerebral cortex (Braitenberg and Schüz, 2013;Lefort et al., 2009;Binzegger et al., 2004;Thomson and Lamy, 2007), and theoretical work has shown that the strength of recurrent coupling has a major impact on several computational properties of networks of excitatory (E) and inhibitory (I) neurons. Strong excitatory recurrent coupling can increase the speed of network response to external stimuli (van Vreeswijk and Sompolinsky, 1996;van Vreeswijk and Sompolinsky, 1998), allow a network to sustain persistent activity (Amit and Brunel, 1997), and increase the capacity and robustness of memory storage (Rubin et al., 2017). Further, strong recurrent coupling can allow networks to amplify specific input patterns, as well as generate complex spatiotemporal activity patterns in response to briefer inputs (Murphy and Miller, 2009;Hennequin et al., 2018).
Strong excitatory recurrent coupling, however, can lead to unstable dynamics unless stabilized by inhibition. When recurrent connectivity is weak, excitatory cells can show stable firing rates independent of the activity of inhibitory cells. But in networks with strong recurrent connections, excitatoryto-excitatory (E-E) connections amplify responses so that the excitatory network is unstable if the firing rates of inhibitory neurons are kept fixed. Stable excitatory network operation across a range of firing rates can be restored if inhibitory recurrent connections are sufficiently strong, allowing inhibition to track and balance excitation (Amit and Brunel, 1997;van Vreeswijk and Sompolinsky, 1996;van Vreeswijk and Sompolinsky, 1998;Renart et al., 2010). Such network models, with strong recurrent connections rendering the excitatory cells self-amplifying and thus unstable, and requiring inhibition for stability, are called inhibition-stabilized networks (ISNs) (Ozeki et al., 2009).
Whether cortical networks function in the ISN regime, in which conditions they do so, and which cortical areas may operate as ISNs remain the subject of debate. One key open question has been whether cortical networks function as ISNs only during high levels of network activity (as when strong sensory stimuli are used to drive sensory cortical regions), or also operate as ISNs for weak input or even during spontaneous activity states. Merely because a network is stabilized by inhibition during one network state does not imply it is inhibition-stabilized for all network states (Ahmadian et al., 2013). For example, the fact that the loss of cortical inhibition can lead to epileptic seizures (Steriade and Contreras, 1998) may seem at first sight to imply that cortical networks are inhibition-stabilized. However, this observation leaves open the possibility that cortical networks are inhibition-stabilized only in some network states but not all, and that those inhibition-stabilized states are the ones that generate seizure activity.
One area, cat primary visual cortex (V1), shows behavior clearly consistent with the ISN regime (Ozeki et al., 2009). This was established in the presence of sensory stimuli, and thus could not determine whether cat V1 operates as an ISN in the absence of visual stimuli when network activity is lower (i.e. at rest, during spontaneous activity). Based on these data and others, Miller and coauthors later developed a model called the stabilized supralinear network (SSN) (Ahmadian et al., 2013;Rubin et al., 2015), which explains cat and ferret V1 responses to visual sensory stimuli of different sizes and intensities. The SSN shows ISN-regime operation for strong visual stimuli, but predicts that as sensory stimuli decrease in strength, the network should transition into a non-ISN state.
To examine whether cortical networks operate as ISNs or non-ISNs in different activity levels, optogenetic stimulation can be used to test predictions of ISN models. A strong prediction of inhibition stabilization is paradoxical inhibitory suppression -a counterintuitive decrease of inhibitory firing rate as inhibitory cells receive direct input ( Figure 1A-B, Figure 1-figure supplement 1; . The paradoxical effect is a result of the strong synaptic coupling between excitatory and inhibitory cells. In an ISN, due to the strong recurrent coupling, more recurrent excitation is withdrawn after inhibitory stimulation than the stimulation adds, leading to a net decrease in input to stimulated inhibitory cells (Figure 1-figure supplement 1). This decrease in recurrent excitation is the cause of the paradoxical effect. Thus, examining whether inhibitory cells show paradoxical suppression when stimulated should be able to reveal whether cortical areas are operating in the ISN regime.
However, optogenetic ISN studies, largely performed in the mouse, have produced conflicting data about ISN operation during spontaneous activity states. While some paradoxical changes in  inhibitory currents have been observed in mouse primary auditory cortex (A1) (Kato et al., 2017;Moore et al., 2018), several studies have found, in contrast, non-paradoxical effects in inhibitory firing rate or intracellular currents in mouse V1, A1, and primary somatosensory cortex (S1) (Atallah et al., 2012;Moore et al., 2018;Gutnisky et al., 2017). One study that observed paradoxical changes in firing in mouse A1 via stimulation with ArchT (Moore et al., 2018) found that those paradoxical effects were not produced via the ISN mechanism, but instead mediated by a circuit where L4 inhibitory cells suppressed L2/3 inhibitory cells (i.e. using feedforward inhibition). Possible explanations for these varying observations include differences between anesthetized and awake states, differences between measurements of firing and intracellular currents, or, as we show experimentally below, differences in which, or how many, inhibitory cells are stimulated (Litwin-Kumar et al., 2016;Sadeh et al., 2017;Gutnisky et al., 2017). Additionally, in mouse V1, one study (Adesnik, 2017) found that excitatory and inhibitory currents scale as predicted in the SSN with increasing visual stimulation. But this study of SSN phenomena also leaves open the question whether mouse V1 operates as an ISN without sensory stimuli (see also Litwin-Kumar et al., 2016, their Discussion).
Here, we examine whether cortical networks operate as ISNs at rest, by combining optogenetic stimulation of inhibitory cells, extracellular recordings in awake animals from several cortical areas, and theoretical analyses. We find clear evidence for ISN operation even without sensory stimulation in multiple mouse cortical areas and both in superficial and deep layers. By providing multiple experimental and theoretical lines of evidence for ISN operation, our results argue against non-ISN explanations (e.g. against paradoxical suppression arising from one inhibitory population inhibiting another).
We also address whether stimulating a single genetic subclass of inhibitory neurons in V1, the parvalbumin-positive (PV) cells, can produce paradoxical inhibitory effects. We find differences in paradoxical responses between viral and transgenic expression strategies. These differences are explained by an ISN model where the fraction of stimulated PV cells is varied (Gutnisky et al., 2017;Sadeh et al., 2017), and are supported by histological measurements, supporting the idea that cortical responses to inhibitory input can be highly dependent on the number of stimulated cells (Sadeh et al., 2017). Additionally, while synaptic plasticity (Varela et al., 1999) such as depression at excitatory synapses can impact network stability, we show via theoretical work (Appendix 1) that the paradoxical effects we observe imply stabilization by inhibition, even in the presence of synaptic plasticity.
Together, these results support the idea that the cortex operates in the ISN regime across a range of areas and brain states.

Mouse primary visual cortex is inhibition stabilized
We first describe experiments where all inhibitory neurons were optogenetically stimulated. We expressed an excitatory opsin (ChannelRhodopsin-2; ChR2) in all inhibitory neurons using a transgenic mouse line (VGAT-ChR2). We delivered blue light to the surface of the cortex, while recording activity extracellularly with multi-site silicon recording arrays ( Figure 1, Figure 1-figure supplement 2). Mice were given drops of reward approximately once per minute (Materials and methods) to keep them awake and alert. Because neurons in the superficial cortical layers showed the largest responses to the stimulation light (Figure 2-figure supplement 1), we restricted our analyses to units within » 400 mm of the cortical surface (Materials and methods), which primarily includes neurons in layer 2/3 and upper layer 4 (Harris et al., 2018). We sorted units into single-(separated from noise and other units) and multi-units (lower SNR or containing several apparent units) by waveform (Materials and methods; using a more restrictive threshold for unit inclusion does not affect our conclusions, Figure 2-figure supplement 2), and we excluded the multi-units from further analyses.
The majority of recorded units were suppressed by light (146/167 single units, 87%; Figure 2B-D; 4 animals, 7 recording days). Given that approximately 80% of cortical neurons are excitatory (DeFelipe et al., 2013;Tremblay et al., 2016), and that extracellular recordings can show several biases in sampling cortical neurons (Margrie et al., 2002;Olshausen and Field, 2005), this measured percentage could be consistent with either the presence or absence of paradoxical  We first recorded neurons' responses to optogenetic stimulation, and then applied blockers of excitatory and inhibitory synapses (CNQX, AP5, bicuculline; which block AMPA, kainate, NMDA, and GABA-A synapses; Methods). When recurrent synapses are substantially suppressed, inhibitory cells expressing ChR2 should increase their firing rate to optogenetic input, while excitatory cells' firing should not increase. Indeed, by measuring responses to inhibitory stimulation in the presence of E and I blockers, we were able to classify units into two groups ( Figure 2E). Units in the first class (putative excitatory) were silent for strong stimulation and units in the second class (putative inhibitory) increased their firing rate in response to stimulation. To quantitatively classify cells, we measured the change in activity induced by light when the blockers had been applied ( Figure 2F-G); cells were labeled inhibitory (33%; 56/167 units) if the change compared with baseline produced at maximum laser intensity was positive (according to Welch's t-test, p<0:01), and excitatory otherwise. Choosing different classification thresholds did not affect our results (Figure 2-figure supplement 4), examining a subset of cells with the most stable waveforms over time did not affect our results ( Figure 2-figure supplement 5), and, finally, examining a subset of inhibitory units where waveforms had high signal-to-noise ratio also did not qualitatively change the results (Figure 2-figure supplement 2). Therefore, we used these response differences in the presence of blockers to identify recorded inhibitory units.
Pharmacologically identified inhibitory and excitatory units showed differences in waveform width and spontaneous firing rate ( Figure 2H), although these two factors were insufficient alone to completely classify inhibitory cells (i.e. insufficient to predict whether units would increase firing rate in the presence of blockers). Identified excitatory cells had a wider waveform (0.58 ± 0.01 ms) and lower firing rate (6.0 ± 0.8 spk/s) than inhibitory cells (width 0.46 ± 0.02 ms, rate 10.7 ± 1.7 spk/s). However, inhibitory neurons' widths were broadly distributed, supporting the idea that that some inhibitory classes can have broader waveforms than others (McCormick et al., 1985;Neske et al., 2015).
Once recorded units were classified as excitatory or inhibitory, we could ask whether the same units before adding blockers (i.e. with the animal in the awake state, with recurrent connections intact) showed paradoxical inhibitory effects. Indeed, inhibitory neurons showed strong suppression when stimulated ( Figure 2I-J, Figure 2-figure supplement 7). The mean response over all recorded inhibitory units was suppressed ( Figure 2L), and also, the large majority of inhibitory units were suppressed (distributions of response slopes shown in Figure 2J-K; responses for all inhibitory cells as power is varied shown in Figure 2-figure supplement 7). This paradoxical inhibitory suppression is one piece of evidence that an ISN model is a good description of the upper layers of primary visual cortex. (Below we show that response dynamics and responses at high laser intensity provide additional evidence for an ISN.)       These experiments were done without sensory stimulation, as animals viewed an unchanging neutral gray screen. Because visual cortical neurons' firing rates increase with increasing contrast across a range of overall luminance levels, we expected that we would see little difference between data collected with a neutral gray screen and data collected in the dark. Confirming that expectation, in one experimental session we measured inhibitory responses to stimulation while animals were in the dark (with the blue optogenetic light shielded with light-blocking materials), and found paradoxical suppression (N = 6/6 inhibitory units with initial negative slope, median initial slope À6.9 ± 3.5 spk/s; p = 0.046, median < 0, Mann-Whitney U test).
In an ISN, paradoxical inhibitory suppression arises because exciting inhibitory cells in turn suppresses excitatory cells, which withdraw excitation from the stimulated inhibitory cells. Then, because more excitatory input from recurrent sources is withdrawn from inhibitory cells than excitatory input from optogenetic stimulation is added, the net steady-state response of inhibitory cells is suppression. For this mechanism to hold, inhibitory cells must transiently increase their firing, even if only a small number of inhibitory spikes are generated, before the paradoxical suppression Ozeki et al., 2009). Therefore, we looked for transient increases in inhibitory firing with short latencies after stimulation, and found a small increase in inhibitory firing rate immediately after stimulation, before the larger paradoxical suppression (Figure 3). Another ISN prediction is that inhibitory cells should increase their firing rates at high laser intensities. In a recurrent network operating as an ISN, stimulating inhibitory cells first decreases inhibitory firing rates, but stronger stimulation eventually suppresses excitatory cells until they reach a firing rate where the excitatory network is stable without reactive inhibition. Then, increasing stimulation of inhibitory cells beyond this point produces increases in inhibitory firing rates, as the network moves into a non-ISN regime ( Figure 1C). Indeed, we observed this in our excitatory and inhibitory identified populations ( Figure 2H).
The link between paradoxical responses and inhibition stabilization was initially derived  assuming recurrent connections do not change their strength. Short-term synaptic plasticity (Tsodyks and Markram, 1997;Markram et al., 1998) could modify this link, as synaptic plasticity can influence network stability, for example by reducing the self-amplification of the excitatory network as excitatory rates rise. We thus set out to determine analytically how synaptic plasticity might affect both responses and dynamics in a network model of excitatory and inhibitory neurons (Appendix 1). We find that, even when short-term synaptic plasticity is included in the model, paradoxical responses imply inhibitory stabilization. Therefore, even when our data are interpreted in a modeling framework that includes short-term plasticity, our experimental observations provide evidence for inhibition stabilization.
While our extracellular recording approach has advantages, such as allowing precise measurements of spike rate dynamics, and recording a large number of E and I neurons in the same experiment, one concern might be that misclassification of excitatory waveforms into inhibitory units could mask inhibitory cell responses. However, confirming that our recording procedures accurately measured inhibitory waveforms, (1) paradoxical suppression is maintained when we make our classification thresholds more stringent (restricting analysis to units with highest waveform signal-to-noise ratio, Figure  Also, as shown below, in parvalbumin-positive neuron stimulation with viral transfection, using the same waveform-sorting procedures, we observed no average paradoxical effects. Together, these observervations argue that paradoxical effects in VGAT-ChR2 animals do not arise due to waveform sorting, but instead due to the properties of the cortical inhibitory network.
Taken together, these results show that mouse V1 responses are consistent with those of a strongly-coupled network whose activity is stabilized by inhibition.

Model-based inference of network coupling strength and stability
We next used these data to infer network connection parameters in a standard two-population model that describes the dynamics of populationaveraged firing rates of excitatory and inhibitory neurons (see Equation 6; Materials and methods). This parameter inference is possible because we obtained E and I response measurements in three pharmacological conditions (no blockers; with E blockers; with E and I blockers; Figure 4A), and with data from these three conditions, there are more experimental observations than model parameters (Materials and methods: Model degrees of freedom). The model uses a rectified-linear single-neuron transfer function; using nonlinear transfer functions give similar results (Figure 4-figure supplement 1 and Appendix 2).
First, we find that the model, despite being overconstrained by the data, is a good description of V1 responses ( Figure 4A). The model also makes a number of predictions that are verified in the data. It predicts that the inhibitory firing rate slope for high laser intensities should be the same in the no-blocker and E-blocker conditions, as seen in the data ( Figure 4A, blue: left vs. middle columns; see also Appendix 2-figure 1). The model makes three other predictions: that the reversal point of inhibitory cells (L 0 , the point at which initial negative inhibitory slopes become positive), in all three experimental conditions, should match the point where excitatory cells change their slope to become nearly silent. These effects are also seen in the data (e.g. clearly visible in the no-blocker case, Figure 4A, left column). Note that, although these predictions have been derived in a threshold-linear model (Materials and methods), they are general features of network models which show a transition from non-ISN to ISN as the activity level of the excitatory population increases.
One prediction of balanced-state models with strong recurrent coupling (van Vreeswijk and Sompolinsky, 1996) is that the total excitatory and inhibitory currents should each be large, although they cancel to produce a small net input. To check this, we computed from the model the currents flowing into excitatory cells ( Figure 4B). The excitatory and inhibitory currents are each approximately ten times larger than threshold (medians are 10.9 and À8.63, respectively), and their sum is small compared to the magnitudes of E and I currents (median net current/threshold = 2.2). The model also allows  Figure 3. Inhibitory neurons show a small increase in firing before paradoxical suppression. (A) Average timecourse of neural responses from units classified as excitatory (N = 111) and inhibitory (N = 56). Light pulse is constant-intensity and lasts for 800 ms (high-intensity pulse, 2.6 ÁL 0 ; stimulation strength larger than L 0 is predicted by ISN models to produce larger transients than at L 0 , though steady-state I rates are larger than the minimum; see also Figure 3-figure supplement 1). Data from VGAT-ChR2 animals. Heavy lines: population mean rate, smoothed with a LOWESS filter. Shaded regions: SEM (Materials and methods). The onset of the light pulse produces an artifact in the 0 ms time bin; we drop spikes at that time point (broken lines, time = 0) to remove the artifact. Inhibitory neurons show a brief increase in firing before suppression, and the increase above baseline is statistically significant (inhibitory firing rates over 12 ms window after pulse are greater than a matchedduration interval before pulse, p<0:01, KS test). (B) Same data as (A), enlarged to show the initial transient Figure 3 continued on next page us to infer the value of the key parameter for inhibition stabilization, W EE (the excitatory selfamplification). The excitatory network is unstable, requiring inhibition for stabilization, and yielding an ISN, when W EE >1. As expected when paradoxical effects are present, the fitted value is greater than 1 ( Figure 4C: median 4.7, mode 2.5). This value is not much larger than unity, suggesting recurrent coupling is not as strong as some balanced-state models might predict (van Vreeswijk and Sompolinsky, 1996;van Vreeswijk and Sompolinsky, 1998;Amit and Brunel, 1997) but is more consistent with a moderately-coupled network (Ahmadian et al., 2013;Rubin et al., 2015;Hennequin et al., 2018). In Appendix 2, we extend this model analysis to examine networks of spiking neurons which feature a fixed level of input noise, and find qualitatively similar results. The model shows that measured activity is consistent with the mean input to cells being below threshold. This implies the recorded cells operate in the fluctuation-driven regime, matching expectations from cortical response variability (Shadlen and Newsome, 1994;Shadlen and Newsome, 1998;van Vreeswijk and Sompolinsky, 1996;van Vreeswijk and Sompolinsky, 1998;Amit and Brunel, 1997). Moreover, when input fluctuations are taken into account in the model, mean excitatory and inhibitory currents are found to be of order of the distance between rest and threshold, in agreement with recent estimates suggesting the cortex operates in the loosely balanced regime (Ahmadian and Miller, 2019).
Finally, we computed the parameter ranges over which the network is stable, as a function of the time constants of E and I firing rate dynamics, t E and t I . We find, by computing the eigenvalues of the connection matrix (Materials and methods) that the network is stable for t I =t E in the interval [0, 5.3], and the observed dynamics are consistent with t I =t E ¼ 4:4 and t E ¼ 7:8 ms ( Figure 4E); that is, with excitatory and inhibitory synaptic time constants of the same order. The values of t ðE;IÞ ( Figure 4D) inferred from the response dynamics are within the stable region, providing further support that the ISN model is a good description of the effects we measure. (Note that the time constants of the opsin may increase these estimates, so the values of t E and t I , while small, are only upper bounds on the true dynamic response of the network to an instantaneous conductance change.) The ratio between E and I time constants is close to the point at which the network becomes unstable through an oscillatory instability ( Figure 4E). The fact that the inferred network model describes the data and is stable with W EE >1 supports the experimental observation that the underlying cortical network is inhibition-stabilized.

Inhibition stabilization is present in other cortical areas
The data above show that the superficial layers of mouse V1 are inhibition stabilized. Another open question is how general the ISN regime is across cortical areas. We performed experiments to look for signatures of inhibitory stabilization in two other cortical areas ( Figure 5): somatosensory (bodyrelated primary somatosensory, largely medial to barrel cortex) and motor/premotor cortex. (See Figure 1-figure supplement 2 for locations of all recording sites.) As in V1, we recorded units extracellularly, restricted analysis to the superficial layers of the cortex, and examined responses to optogenetic stimulation of all inhibitory cells (as above, using the VGAT-ChR2 line).
In these experiments, we classified neurons as inhibitory based on their activity at high laser intensity. In our V1 experiments, paradoxical effects are clear whether neurons are classified in this way, or with the pharmacological method ( Figure 2-figure supplement 3). This high-laser-power classification is likely more stringent (i.e. may reject some inhibitory cells if they are suppressed by other (time range here is indicated in (A) by gray shaded region). Inhibitory initial positive transient has peak amplitude 3.6 spk/s above baseline, latency to peak 7 ms, full width at half maximum 7.1 ms, and inhibitory rate crosses baseline into suppression at 13.1 ms. (C) Size of transient vs. size of steady-state suppression for all I units. The majority of units show steady-state paradoxical suppression relative to both baseline rate and to transient rate. y-axis: steady-state rate change (normalized to baseline; rateÀbaseline baseline ). x-axis: normalized transient; measured from 1 to 12 ms after light onset. Units below the horizontal line at zero show paradoxical suppression; units in the blue shaded region are those whose normalized steady-state rate is lower than their normalized transient rate.   Figure 3. Best-fit values are t I ¼ 34:3 ms and t E ¼ 7:8 ms; note that network response dynamics shows faster time constants due to recurrent network effects. The full model provides a good approximation to the dynamics even though it is constrained to simultaneously fit the time constants and the responses at different intensities. (E) Stability analysis. Real (black) and imaginary (purple) parts of the eigenvalues of the Jacobian matrix as a function of the ratio t I =t E . Imaginary part (purple) greater than zero signifies the network can show damped oscillations when being driven to a new stationary point. Note that these damped oscillations are not seen in (D) because of rectification. When the real part (black) is greater than zero, the network is unstable (shaded red area). The online version of this article includes the following figure supplement(s) for figure 4: Figure supplement 1. Comparison of linear and nonlinear rate models, and global optimization method. Figure 4 continued on next page inhibitory cells, due e.g. to heterogeneity of recurrent connections) than pharmacological classification. Supporting that this classification may be more stringent, some pharmacologically-classified inhibitory cells show paradoxical effects but little increase at high power (Figure 2-figure supplement 8).
In both the recorded areas ( Figure 5), inhibitory neurons showed paradoxical effects, supporting that these cortical areas do also operate as an ISN. Similar to the V1 data, and as expected in an ISN, both areas showed a transition to a non-paradoxical response at large laser intensity. Compared to the somatosensory data, the motor/premotor recordings show a smaller average suppression ( Figure 5B) and more inhibitory units with non-paradoxical effects (rate increases with stimulation; Figure 5D). However, in both areas, means and medians are significantly negative (see legend for statistical tests). Thus, the superficial layers of both sensory cortical areas we examined (visual and somatosensory), and one non-sensory area (motor/premotor) showed responses consistent with strong coupling and ISN operation.

Differences in paradoxical suppression with viral or transgenic opsin expression can be explained by different numbers of stimulated cells
Up to this point, we have studied paradoxical suppression by stimulating an opsin expressed in all inhibitory cells via the VGAT-ChR2 mouse line. A remaining question is whether stimulation of a subclass of inhibitory neurons also yields paradoxical suppression. Even in an ISN, stimulation of any single subclass of inhibitory cells need not produce paradoxical suppression (see Materials and methods and Rubin et al., 2015;Litwin-Kumar et al., 2016;Sadeh et al., 2017;Gutnisky et al., 2017).
To study this, we stimulated parvalbumin-positive (PV) neurons, which provide strong inhibitory input to other cells ( Figure 6). PV basket cells are the most numerous class of cortical inhibitory cells    Figure 2), weaker paradoxical suppression is seen (blue line; mean rate is not significantly suppressed: p>0:05, paired t-test between rate at 0, rate at 1 Á L 0 ) as initial response slope is near zero. (E) A model with a subset of inhibitory cells stimulated (60%) can recapitulate the data in (B). Shaded region: ±1 s.d. via data bootstrap as in Figure 4. Other than splitting inhibitory population into two subsets, network parameters (Materials and methods) are as inferred in Figure 4. (F) Population responses for transgenic expression of ReaChR in PV neurons. Unlike the viralexpression data (D), the identified inhibitory cells in these experiments show paradoxical suppression; mean firing rate is significantly suppressed (p<10 À6 , paired t-test on initial slopes). (G) Cell counts in histological sections. Solid/dotted black lines: means across two independent human counters; upper and lower gray boundaries give results from each counter. X-axis shows variation as counters were asked to use a high or low threshold for accepting an opsin positive cell. Cells were counted across 400 mm depth and are expressed (y-axis) as a percent of DAPI-positive nuclei. Transgenic Figure 6 continued on next page and make strong synapses near the somata of excitatory cells (for review, see Tremblay et al., 2016), and PV stimulation effectively suppresses network firing rates (e.g. Glickfeld et al., 2013;Sparta et al., 2014).
We first used viral methods to express an opsin in PV neurons, injecting a Cre-dependent adenoassociated virus (AAV) encoding an excitatory opsin (Chronos; Klapoetke et al., 2014), into a transgenic mouse line (PV-Cre; Hippenmeyer et al., 2005). In these experiments, the network neurons can be divided into three populations ( Figure 6A): (1) excitatory, (2) Chronos-expressing PV inhibitory (PV-Chronos), and (3) remaining inhibitory neurons: non-PV (e.g. somatostatin-positive, etc.) or non-Chronos-expressing PV neurons. We identified PV-Chronos cells by measuring whether cells' firing is increased at high laser intensity (statistically significant increase at maximum laser power, as in Figure 5; see Materials and methods). In V1 VGAT-ChR2 experiments, this classification by response at high laser power produces qualitatively similar measures of paradoxical suppression as pharmacology-based classification, Figure 2-figure supplement 3. Non PV-Chronos cells (without increase in firing at high laser intensity), are likely either excitatory cells, non-PV, or non-expressing cells. Supporting the idea that our classification approach identifies PV-Chronos cells, the majority of classified inhibitory cells have narrow waveforms ( Figure 6I). While not all PV-positive cells are basket cells (Taniguchi et al., 2013), the large proportion of narrow-waveform units we found suggest our inhibitory classification detects a number of fast-spiking cells.
Stimulating the PV cells produced no significant paradoxical effect on average ( Figure 6D, blue). From a theoretical point of view, this average inhibitory response is an important measure, as with the standard two-population ISN model, e.g. , it is the average response (averaged over inhibitory cells) that is paradoxically suppressed when a network is inhibition-stabilized and strongly coupled. A second important measure, e.g. for experimentalists wishing to identify inhibitory neurons, is how many individual inhibitory cells show paradoxical effects. Examining individual units showed that PV-Chronos units often showed no paradoxical effect. Initial slope was often positive and thus non-paradoxical (e.g. Figure 6H; median not sig. dif. from zero; summarized in Figure 6K-L; see legend for statistical tests).
One reason why stimulation with viral expression might produce no paradoxical effect on average is that viral expression could target only a subset of the PV cells (Mathematical methods; [Sadeh et al., 2017]). We first examined the effect of inhibitory cell subset stimulation in a model, the rate model of Figure 4, using network parameters inferred there. Stimulating a subset of inhibitory neurons in the model (60%) reduced the magnitude of the average paradoxical effect ( Figure 6E, analytical derivation in Materials and methods, and Figure 6-figure supplement 1), and described the data ( Figure 6D) well. We also considered a model with multiple subclasses of inhibitory neurons, using the connectivity structure measured in Pfeffer et al., 2013, and found there also that paradoxical effects were not clear when approximately half of the PV cells were stimulated ( Figure 6-figure supplement 1).
Why, intuitively, can stimulation of a subset of inhibitory cells eliminate paradoxical effects? In networks of strongly-coupled excitatory and inhibitory neurons, increased inhibitory activity produced by stimulation suppresses excitatory activity and results in a withdrawal of recurrent excitation; this, in turn, drives suppression of stimulated cells (the paradoxical effect). When only a fraction of inhibitory cells is stimulated, the withdrawal of excitation coexists with other effects produced by Figure 6 continued expression gives about twice as many opsin expressing cells as viral expression, in addition to the differences in expression heterogeneity seen in (B-C). (H) Example viral (PV-Chronos) units show diversity of responses to stimulation. Some narrow spiking units (blue, rightmost two panels) show nonparadoxical initial increases, and some show paradoxical initial suppression. (I,J) Distribution of spike widths (I) and initial slopes (J) for experiments using transgenic PV-ReaChR (bottom). Same conventions as in Figure 5, except here colors are as shown in panel (A). (K,L) Same as (I,J) for viral expression. Although inhibitory units are classified by response at high laser power, differences in spike width are visible in both datasets. Viral mean and median slopes are zero or positive (t-test for negative mean p>0:05, Mann-Whitney U for negative median p>0:05; 22/42 (53%) negative I slopes); transgenic mean and medians are negative (negative mean p<10 À7 , negative median p<10 À7 ; 24/27 (89%) negative I slopes). (M) Summary of mean and median inhibitory cell initial slopes for these data and data from V1, somatosensory, and motor cortex. The online version of this article includes the following figure supplement(s) for figure 6:  recurrent interactions. In particular, increased activity of stimulated inhibitory cells tends to suppress non-stimulated inhibitory cells and to produce a withdrawal of inhibition to excitatory cells and to stimulated inhibitory cells. This withdrawal of inhibition increases with the fraction of non-stimulated inhibitory cells and, if large enough, can overcome the withdrawal of excitation and prevent paradoxical suppression of stimulated cells.
To test these models, and experimentally determine whether differences in number of cells expressing opsin could change the paradoxical effect, we used a transgenic approach to express a different excitatory opsin in PV cells (PV-Cre;ReaChR transgenic mice). In these mice, the opsin is expressed in most or all PV cells (Lin et al., 2013). (Since we study steady-state firing rates (Materials and methods), we do not expect differences in the onset or offset kinetics of the opsin to affect these measurements, and we verified that considering different time windows of the steadystate response does not change the results, Figure 6-figure supplement 2). Because we measured responses for a range of light intensities, differences in viral and transgenic mean levels of opsin expression would not affect paradoxical effects. On the other hand, differences in the number of expressing cells, or variability in opsin levels across cells, are predicted to change the paradoxical effect. Indeed, we found that, with transgenic expression, the average paradoxical effect was present ( Figure 6F). Further, the fraction of inhibitory cells that showed paradoxical responses was significantly larger in PV-ReaChR transgenic animals than in PV-Chronos animals ( Figure 6J,L), and was similar to the VGAT-ChR2 data in V1, somatosensory, and motor/premotor cortex ( Figure 6M).
Changes in the number of stimulated cells could be affected by the virus infecting only a subset of neurons (Watakabe et al., 2015), or by the virus yielding different levels of opsin expression in different neurons, so that only a subset are recruited strongly at any given light intensity. To examine differences in expression pattern, we counted neurons in histological sections from both the viral (PV-AAV-Chronos) and transgenic (PV-ReaChR) animals. For these comparisons, the Cre line, expressing Cre in PV neurons, was held constant. We observed both effects. First, viral expression created some cells with strong, and some cells with weak, opsin expression ( Figure 6B-C). Second, viral expression yielded opsin in about half as many cells as transgenic expression ( Figure 6G); positive-opsin cell percentage significantly different between viral and transgenic cases, both observers and both thresholds, p<0.05, c 2 test). The differences in the number of opsin-expressing cells roughly matched our model predictions (60% of neurons stimulated, Figure 6E).
Taken together, these data and our analysis support the idea that V1 and other areas of the mouse cortex are strongly coupled and operate as an ISN. Yet whether an average paradoxical effect is seen depends on the number of inhibitory cells stimulated. The heterogeneous responses we saw in PV-Chronos experiments (about half of inhibitory cells paradoxical, half non-paradoxical, Figure 6D,H,L, also may in part explain why paradoxical inhibitory suppression,and thus ISN operation, has not been more widely reported in optogenetic PV-stimulation experiments (also see Discussion).

Paradoxical effects are also seen in deep-layer recordings
To this point all the neurophysiological data we have reported is from the upper layers (L2/3 and 4) of the cortex, recorded within 400 mm of the cortical surface. We focused on the superficial layers because the blue light we used to activate ChR2 and Chronos does not penetrate more than a few hundred microns into the tissue (Yona et al., 2016), so the inhibitory neurons that receive direct optogenetic input are those in the superficial layers. However, blue light delivered to the cortical surface to stimulate inhibitory neurons has been seen to suppress activity across cortical layers, presumably due to polysynaptic effects (Li et al., 2019). Because we recorded data using silicon probes that span most of the depth of the cortex, we could also assess whether light delivered to the top of the cortex also produced paradoxical suppression in deeper inhibitory neurons.
We examined units (Figure 7) recorded ! 500 mm from the cortical surface (Materials and methods). In our multi-area recordings from superficial layers (e.g. Figure 5), we classified inhibitory cells in different brain areas based on increases in inhibitory response at strong light intensities, when excitatory cells are silent or largely suppressed. But in deep recordings, presumably due to attenuation of stimulation light from scattering and absorption, at our highest stimulation intensities we did not always see clear mean increases in inhibitory firing rates (e.g. Figure 7B). Therefore, for these deep recordings we classified neurons as inhibitory based on waveform width. As in superficial layers, we found bimodal distributions of waveform width in our deep-layer data ( Figure 7A-C, insets). In our upper-layer V1 data, we compared classification of V1 units via pharmacology, waveform width, and responses at high power (Figure 2-figure supplement 3). Based on this comparison, in these deep-layer data, we do not expect excitatory cells to be classified as inhibitory (i.e. few or no excitatory units have narrow waveforms), even though some inhibitory cells with wider waveforms will be missed with this approach (i.e. some inhibitory units have wide waveforms).
Examining cells recorded from deeper layers, we found that V1, motor, and somatosensory areas all showed initial paradoxical responses (Figure 7). However, the maximal suppression of inhibitory and excitatory cells was not as strong as in superficial layers (Figure 2, Figure 5). (Supporting the idea that this weaker effect of stimulation is due to attenuation of blue light, paradoxical effects in deep layers were smaller in these data than when using red light to stimulate ReaChR in V1 PV  Figure 2K. As in that panel, slopes here are normalized by baseline rate so that minimum slope is À1. Means and medians of individual inhibitory neurons' slopes are all negative (t-test, V1 p<0:05, others p<0:01) except for V1 median (blue error bar shows 95% confidence interval around median via bootstrap: upper CI, V1 0.03, motor À0.22, somato À0.33). Also, the population firing rate decrease is significant for all three areas (p<0:001, except V1 p<0:05, Mann-Whitney U on summed population counts, baseline vs. rate at L 0 ). The online version of this article includes the following figure supplement(s) for figure 7: neurons, Figure 7-figure supplement 1.) In all these cases ( Figure 7A-C), the mean population inhibitory responses showed initial suppression.
Thus, the deep-layer data shows no evidence for a different pattern of responses than seen in the upper layers -deep-layer inhibitory cells show paradoxical suppression to excitatory stimulation. However, because light applied to the cortical surface affects the superficial layer inhibitory neurons most directly, we do not make the claim that these data are final proof that deep layers operate in the ISN regime. On the other hand, the deep-layer data do provide evidence against alternative, non-ISN, models for L2/3 responses, such as the potential alternative where deeper inhibitory cells increase their firing and inhibit L2/3 cells. Instead, we find that most inhibitory cells we record, across layers, are paradoxically suppressed by stimulation.
No transition out of inhibition stabilization is seen during anesthesia, at lower network firing rates Theoretical studies have pointed out that a network can switch from non-ISN to ISN when external inputs are increased (Ahmadian et al., 2013). If the cortical network did transition out of the ISN state at lower activity levels, this would be computationally important, as the way inputs sum can change depending on whether the network operates as an ISN or not (Ahmadian et al., 2013;Rubin et al., 2015). Our results above (Figure 2, Figure 3) show that superficial layers of mouse V1 operate as an ISN at rest (i.e. without sensory stimulation), and also show the transition from ISN to non-ISN is below the level of spontaneous activity, as network activity is decreased by optogenetic stimulation of inhibitory neurons. (The transition point is where the slope of I rates vs stimulation intensity switches from being negative to positive.) However, we wished to determine whether a brain state with lower activity levels, as that produced by light anesthesia, might result in a transition to a non-ISN regime.
Under light isoflurane anesthesia, we found no evidence of a transition, and instead found that paradoxical inhibitory suppression was maintained ( Figure 8A). At 0.25% isoflurane (a low level, as surgical levels are often 1.0% and above), spontaneous firing rates of excitatory neurons are reduced ( Figure 8C; mean 6.4 spk/s reduced to 3.5 spk/s, p<0:02, Wilcoxon signed-rank test; inhibitory neurons' rates show a negative trend but are not statistically different, 14.5 spk/s awake to 11.6 spk/s anesth., p<0:10, Wilcoxon signed-rank tesk). Thus, the network changes induced by anesthesia do not cause the network to transition out of the ISN state. At this low level of anesthesia, we did not observe prominent up and down state slow oscillations. In one experiment, we used a higher level of anesthesia (0.5% isoflurane, Figure 8B), yielding even lower firing rates but still preserving the paradoxical effect. Further confirming the robustness of coupling to changes with anesthesia, the distribution of response slopes is roughly unchanged ( Figure 8C), suggesting the network is far from a transition into a non-ISN state. Anesthesia thus preserves the paradoxical inhibitory response, leaving the network still an ISN.

Discussion
These data show a signature of inhibitory stabilization, the paradoxical suppression of inhibitory cells to optogenetic excitation, in several different areas of mouse cortex, V1, S1, and motor cortex, during spontaneous activity in the absence of sensory stimuli. An important aspect of the work is that we stimulated all inhibitory cells together, using a mouse line in which all inhibitory cells express opsin. When single inhibitory subclasses are stimulated, paradoxical effects do not prove that the network's excitatory cells are unstable without inhibition (Mathematical Methods; also see Litwin- Kumar et al., 2016;Mahrach et al., 2020). However, our observation of paradoxical suppression when stimulating all inhibitory cells together, supported by a model that describes the data well (Figure 4), is evidence that these cortical networks do indeed have strong excitatory recurrent coupling, and do indeed operate in the ISN regime.
We found, as predicted by ISN models, that strong inhibitory stimulation causes the network to transition into a non-ISN state. In this non-ISN state, inhibitory responses are non-paradoxical and increase their firing when stimulated, and excitatory neurons are suppressed. This transition from ISN to non-ISN behavior is a prediction of ISN models. We saw paradoxical suppression of the mean inhibitory firing response in both superficial and deep layers of the three areas, and we found that a clear majority of inhibitory units (except in one of these six measurements, deep layers of V1) show paradoxical suppression when all inhibitory cells are stimulated. We used several ISN rate models to fit data from V1, aided by pharmacological manipulations that yield additional informative data about the network. We find that despite having fewer parameters than data degrees of freedom, the ISN model describes the data well. During light anesthesia, mean V1 network firing rates decrease, but we still observed ISN paradoxical suppression. Finally, we explored the effect of stimulating a subclass of inhibitory neurons with strong network effects: the PV neurons. We find that stimulating PV cells also produces paradoxical effects. We compared viral and transgenic expression strategies and found that viral expression did not produce average paradoxical effects. This viral/transgenic difference is consistent with theoretical observations that changing the fraction of stimulated inhibitory cells changes the strength of the paradoxical effect, and is supported by histological data.

Multiple lines of evidence rule out other circuits and argue for the strongly-coupled ISN explanation
Our data and theoretical results provide the strongest evidence currently available to rule out non-ISN interpretations. Three prominent features of our data supporting a strongly-coupled recurrent origin for the results are: (1) that inhibitory response dynamics match ISN predictions (Figure 3), (2) that an ISN model describes the data well and makes predictions verified in the data (Figure 4), and (3) that many inhibitory cells show paradoxical effects at low stimulation intensity but increase their firing rate at higher stimulation intensity (e.g. Figure 2, Figure 5).
These observations allow us to rule out several potential alternatives to the ISN explanation for the paradoxical effects we observed. One potential alternative might be a model with multiple inhibitory populations and weak recurrent coupling, perhaps driven to fire spontaneously by strong external input such as thalamic or corticocortical input. In principle, paradoxical inhibitory responses could be produced in such a non-ISN model via feedforward, inhibitory-onto-inhibitory, interactions. Specifically, in this alternative, stimulation of one inhibitory subclass (e.g. PV cells) would inhibit a second set of inhibitory cells whose firing rate would decrease (e.g. somatostatin or SOM cells; the same explanation could hold with the PV and SOM classes swapped). However, in a feedforward inhibitory model, there would be little reason to predict inhibitory neurons to change from decreasing to increasing rate as stimulation intensity increases. ISN models, on the other hand, are predicted to transition into a non-ISN state for strong inhibitory stimulation where inhibitory cells begin to increase their firing, as our data show (e.g. Figure 2L).
The large fraction of inhibitory cells that we found to be suppressed (paradoxical) at low stimulation intensity ( Figure 2J-K)   Population average response with (gray) and without (blue;red) anesthesia. Spontaneous firing rates are reduced for both the E (top) and I (bottom) populations below the awake firing rates, yet inhibitory paradoxical suppression is preserved with anesthesia (lower panel: gray line shows negative initial slope). In these experiments, inhibitory responses are weaker at high powers, but this does not affect paradoxical suppression. (C) Distribution of initial slope and spontaneous rate before and after anesthesia. Initial slope is largely preserved while excitatory and inhibitory spontaneous rates are suppressed. Colored horizontal bars within each panel show mean ± SEM for each distribution. Initial slopes are plotted here non-normalized (units spk/s/mW) to show the slopes are quantitatively similar across firing rate changes; normalized slopes are shown in created by multiple differentially-responding inhibitory cell populations. While excitatory opsins expressed in inhibitory neurons can sometimes cause synaptic release at terminals (Zhao et al., 2011;Babl et al., 2019), we found that many inhibitory cells fired in response to optogenetic stimulation in the presence of E and I blockers, showing that our stimulation did produce substantial neural spiking. Further, ISN-generated paradoxical effects are network effects that rely on inhibitory cells affecting other local neurons, so even if stimulation did cause terminal release, this should not affect our conclusions. Finally, ISN models predict, and we observed, that the transition point from inhibitory suppression to excitation occurs when excitatory cells are largely suppressed. It might be possible for that to occur in a non-ISN case, such as in the feedforward scenario, but it would require fine-tuning of parameters, and our ISN model's accurate description of the pharmacology data makes such fine tuning even more unlikely.
Other non-ISN scenarios would also rely on feedforward inhibitory effects made unlikely by our data. For example, alternative models could have feedforward inhibition arising from rare, powerful inhibitory cells, perhaps so rare that we could not record them, or input from deeper-layer inhibitory cells. These feedforward inhibition scenarios are also argued against by the observations outlined above (e.g. the transition from suppression to excitation showed by our recorded inhibitory cells, etc.). And our recordings across depth showed that a majority of inhibitory cells' firing rates were paradoxically suppressed by stimulation, independent of the depth of recording, making it unlikely that our effects were produced by feedforward inhibition arising from a deeper cortical layer. The ISN model remains a consistent explanation of all the data we observed.

Potential explanations why prior evidence for paradoxical suppression has been mixed
In the last decade, multiple experiments have performed optogenetic stimulation of cortical interneurons, but these studies have yielded mixed results on inhibition stabilization. In layer 2/3 of mouse V1, optogenetic excitation of PV cells with viral expression was shown to generate non-paradoxical modulation of inhibition (Atallah et al., 2012). In mouse auditory cortex (A1), in contrast, Kato et al., 2017 found that suppression of inhibitory cells with viral ArchT produced paradoxical effects in intracellular currents. ISN models make similar predictions for somatic intracellular currents as for the firing rates we measured: if a network is strongly coupled and thus inhibition-stabilized, stimulating all inhibitory neurons will yield paradoxical changes in inhibitory current, as well as a transient increase in inhibitory current before suppression (Ozeki et al., 2009;Litwin-Kumar et al., 2016). Another A1 study observed mixed paradoxical effects but linked them to a non-ISN mechanism. Moore et al., 2018 found that stimulation of PV cells using viral transfection produced paradoxical effects in some PV cells and not others. and they reported evidence for feedforward, non-ISN, mechanisms in L2/3 via L4 inhibitory input. In somatosensory and motor cortex, a recent study observed paradoxical suppression for PV subclass stimulation in some layers and not others (Mahrach et al., 2020). Finally, many experiments have used excitation of inhibitory cells (e.g. Glickfeld et al., 2013) to suppress excitatory activity, without reporting paradoxical inhibitory effects, although a recent survey of such methods does report paradoxical effects (Li et al., 2019) in somatosensory and motor cortex.
Our results, combining theoretical analyses with a range of experimental conditions, suggest some explanations why paradoxical inhibitory suppression has not been previously widely reported.
First, it appears a significant fraction of inhibitory cells must be stimulated to produce large paradoxical changes in firing rate (Sadeh et al., 2017). Stimulating PV cells with viral methods (which can yield subsets of cells with weaker or no expression; see also Sadeh et al., 2017;Gutnisky et al., 2017 for a similar explanation for non-paradoxical effects with viral expression), produced in our hands mixed effects on single cells and weaker mean paradoxical effects ( Figure 6). Further, stimulating subtypes of inhibitory cells even with transgenic expression can produce different inhibitory responses compared to stimulating all inhibitory neurons (Mahrach et al., 2020). Another factor which can, in principle, generate non-paradoxical effects even when the network is an ISN is reducing the area of tissue illuminated by stimulation light. This reduces the number of inhibitory cells stimulated and therefore produces a situation similar to stimulating a subset of network inhibitory cells. However, this effect of reducing optogenetic light area seems to apply mainly for fairly small optogenetic light spots. Sadeh et al., 2017 found via modeling that paradoxical effects were diminished when illumination spot sizes fell to around 100 mm, well below the 500+ mm diameter spots we used here to reveal ISN-related paradoxical suppression. This effect of illumination area may, however, impact experiments using optical fibers, which can have small diameters. Beyond the fraction of neurons stimulated, a second explanation for differences, especially in experiments using extracellular recording, is that it can be difficult to definitively identify inhibitory cells. Some inhibitory cells' waveforms are broad (Figure 2, Figure 2-figure supplement 8), and in an ISN, excitatory and inhibitory cells' firing rate changes are often similar (inhibitory neurons may only show increased activity for a few milliseconds before paradoxical suppression begins). Third, for strong enough drive to inhibitory neurons, ISNs transition to a non-ISN state. In this non-ISN state with low excitatory activity, tracking of excitatory fluctuations is not needed to stabilize the network, and some past experiments may have used stimulation large enough for parts of the network to enter into this second phase of activity, hiding the signature of inhibition stabilization.

PV cells may be the main source of stabilizing input
While paradoxical effects in a single inhibitory subclass do not imply the overall network is inhibition stabilized ( Figure 6E; Mahrach et al., 2020), PV neurons may well be the principal providers of stabilizing input. In an ISN, strong bidirectional E-I synaptic coupling allows inhibitory cells to stabilize an unstable excitatory network, by tracking and responding to changes in local excitation. Thus, one major parameter that would predict which inhibitory cells would show paradoxical effects is the strength of their synaptic influence on other network neurons. PV neurons receive input from many diverse local excitatory cells (Bock et al., 2011) show responses that reflect an average of local excitatory responses (Sohya et al., 2007;Kerlin et al., 2010), cause network instability when strongly suppressed (an effect which does not result from SOM suppression) (Veit et al., 2017), and target virtually all nearby excitatory cells with strong peri-somatic synapses (Fino and Yuste, 2011;Packer and Yuste, 2011). Other cell types of the diverse inhibitory classes in the cortex (e.g. SOM or VIP+ interneurons; Tremblay et al., 2016) also could contribute to stabilization (Sadeh et al., 2017), though if their effect on the whole local network is less strong than PV cells, they may be unable to produce paradoxical suppression when stimulated alone.
Our data support 'loosely balanced' cortical models, and show evidence for ISN operation even at rest Several influential models of cortical function include inhibition-stabilized parameter regimes, and our data provide new constraints on such models. The 'balanced network' model (Amit and Brunel, 1997;van Vreeswijk and Sompolinsky, 1996) predicts that both excitatory and inhibitory inputs to single neurons should be large but approximately cancel each other, leading to a balance between excitation and inhibition. This scenario accounts in a parsimonious way for multiple ubiquitous properties of activity in cortex, such as the highly irregular nature of neuronal firing, and the broad distributions of firing rates across neurons. Balanced network models, however, can vary in their recurrent coupling strength, and thus the size of the total excitatory and total inhibitory currents. The analytical work of van Vreeswijk and Sompolinsky (van Vreeswijk and Sompolinsky, 1996;van Vreeswijk and Sompolinsky, 1998) was performed in the limit of very large numbers of synaptic inputs per cell (K), leading to very large E and I currents (whose leading order in K cancel each other -that is, the sum of the E and I currents is small while the sum's variability can be substantial). Other balanced network models (Amit and Brunel, 1997;Brunel, 2000) used finite, moderate-to-large coupling strength. The network parameters that best fit our data produce moderately large (i.e. within an order of magnitude of threshold) excitatory and inhibitory inputs (Figure 4), which are on the same order of magnitude as the ones used in Brunel, 2000. This moderate recurrent coupling, or 'loose balance', is also consistent with the SSN model (for review, see Ahmadian and Miller, 2019).
The SSN also makes a specific prediction about the ISN regime that our data constrains. The SSN predicts that cortical networks are inhibition-stabilized for strong inputs, but not an ISN for sufficiently weak inputs. Thus, the SSN predicts a transition between two operating regimes as network activity increases. Finding this transition point is important for understanding computation, as the two regimes show different (supralinear or sublinear) modes of input summation. Prior ISN studies (Ozeki et al., 2009) have left open whether this transition occurs only with sensory stimulation. A study of responses to two combined sensory stimuli (Britten and Heuer, 1999) found a linear-tosublinear transition for increasing strength (contrast) of visual stimuli, but did not resolve whether the non-ISN to ISN transition was above or below the level of spontaneous activity. Our data do resolve this, showing several cortical networks are in the ISN state even without sensory stimulation. But a second question our data address is whether during spontaneous activity the network is on the edge of this transition point, or well into the ISN regime. Our work gives two pieces of evidence that without sensory activity the network is far above a transition into a non-ISN state: first, excitatory rates must be substantially suppressed before inhibitory responses switch from paradoxical suppression to firing rate increases (e.g. Figure 2E), and second, under light anesthesia, we find network activity is lowered but ISN behavior is preserved. Therefore, our data provides support for the idea that cortical areas generally operate above this transition point, and that the weakly-coupled, non-ISN regime predicted by the SSN is not a regime mouse cortex commonly enters during normal processing in awake behavior.

Conclusion
Together, the results reported here suggest that inhibition stabilization is a ubiquitous property of cortical networks. This data is consistent with a network that is in the strong coupling regime, but is not extremely strongly coupled as some balanced network models would predict (Amit and Brunel, 1997;van Vreeswijk and Sompolinsky, 1996). It is tempting to speculate that, while strong coupling allows the network to perform non-trivial computations on its inputs, increasing coupling further might not be optimal for several reasons. First, in the strong coupling limit, only linear network responses are available to a balanced network (unless synaptic non-linearities are included, Mongillo et al., 2012), while networks with moderate coupling can combine inputs in a non-linear fashion (Ahmadian et al., 2013). Second, maintaining strong connections might be metabolically expensive, so that a moderately-coupled network might represent an ideal compromise for cortical computations.

Experimental methods Animals
All procedures were conducted in accordance with the guidelines of the National Institutes of Health. Seven VGAT-ChR2 (Zhao et al., 2011; ChR2 targeted at the Slc32a1 locus), three PV-Cre (Hippenmeyer et al., 2005; Cre targeted at the Pvalb locus), and two PV-Cre;ReaChR (Lin et al., 2013; ReaChR targeted at the Gt(ROSA)26Sor locus) mice were used (JAX stock n 014548, 008069, and 024846; 2 females and 10 males; singly housed on a reverse light/dark cycle).

Cranial window implants
Mice were implanted with a titanium headpost and a transparent window (optical glass, 0.8 mm thickness, 3 or 5 mm diameter) in the left cerebral hemisphere. The windows provided access to the primary visual (V1) and somatosensory cortex, or motor cortices, for imaging and silicon neural probe recordings with optogenetics. Mice were given dexamethasone (3.2 mg/kg, i.p.), 2 hr before surgery. Animals were anesthetized during surgery with isoflurane (1.0%-4% in 100% O 2 ). Using aseptic technique, a headpost was affixed using Metabond (Parkell), and a 3 or 5 mm diameter craniotomy made.

Hemodynamic intrinsic imaging
To determine the location of V1, we delivered small visual stimuli to animals at different retinotopic positions and measured changes in the absorption of 530 nm light resulting from cortical hemodynamic responses (Ma et al., 2016). We evenly illuminated the brain with a 530 nm fiber coupled LED (M530F2, Thorlabs) passed through a 532 nm-center bandpass filter (Thorlabs). Images were collected on a stereo microscope (Discovery V8, Zeiss) through a green long-pass emission filter using a 1x objective (PlanApoS 1.0x, Zeiss) onto a Retiga R3 camera (Q Imaging, captured at 2 Hz with 4 Â 4 binning). For retinotopic mapping we presented upward-drifting square wave gratings (2 Hz, 0.1 cycles/degree) masked with a circular window (10˚diameter) for 5 s with 10 s of mean luminance preceding each trial. Stimuli were presented in random order at four positions in the right monocular field of view. The response to a stimulus was calculated as the fractional change in intensity between the average of the 10 frames immediately preceding the stimulus (as baseline) and a 6-10 frame window 1 s after stimulus onset (as response) to match the timecourse of the hemodynamic response (Chen-Bee et al., 2007;Heimel et al., 2007).

Viral injections
We used viral injections to express Chronos (Klapoetke et al., 2014) in PV-Cre animals. Mice were anesthetized (isoflurane 1-1.5%) and the cranial window implant removed. We used a stereotaxic injection system (QSI, Stoelting Inc) to deliver 500 nL of a 10:1 mixture of AAV1-hSyn-FLEX-Chronos-GFP (UNC Vector Core Stock) and 100 mM sulforhodamine (SR101, Invitrogen, for visualization) at a depth of 200-400 mm. We made 3 to 4 injections spaced 0.5-1.0 mm apart to cover the visual cortex. After the injections a new cranial window was affixed. Viral expression was monitored over the course of days by imaging GFP fluorescence, and allowed to develop for greater than 4 weeks before electrophysiological recordings.

Electrophysiological recording with optogenetic activation of inhibitory neurons
For recording experiments, we first affixed a 3D printed ring to the cranial window to retain fluid. With the ring in place, we removed the cranial window and flushed the craniotomy site with sterile normal saline to remove debris. Kwik-Sil silicone adhesive (World Precision Instruments) was used to seal the craniotomy between recording days. Using a stereomicroscope on an articulating arm we positioned the end of an optical fiber (600 mm diameter, Doric Lenses) fitted with a light-tight coupler to an optical cannula (400 mm diameter, Thorlabs CFMLC14L02) over our target cortical area, with a slight 10-30˚angle from vertical to provide space for the electrodes. We used fiber coupled LED light sources (M470F3for Chronos and ChR2, or M625F2 for ReaChR, Thorlabs) to deliver illumination with peak at 470 nm or 625 nm to the brain. We calibrated total intensity at the entrance of the cannula using a power meter and photodiode (meter model 1918R, photodiode model 918D-SL-OD3R, Newport Corp). The cannula distance from the dura was adjusted to provide a light spot with a full width at half maximum intensity of 0.8-1.2 mm (as measured with a small digital macro documentation microscope, Opti-TekScope); thus, our reported spot irradiance in mW/mm 2 is within a factor of two of the reported power (mW). We targeted a multisite silicon probe electrode (Neuro-Nexus; 32-site model 4 Â 8-100-200-177; four shank, eight sites/shank; sites were electrochemically coated with PEDOT:PSS [poly(3,4-ethylenedioxythiophene): poly(styrenesulfonate)], Xiao et al., 2006) to the center of the LED spot using a micromanipulator (MPC-200, Sutter Instruments). With the fiber optic cannula and electrode in place, we removed the saline buffer using a sterile absorbent triangle (Electron Microscopy Sciences, Inc) and allowed the dura to dry for 5 min. After insertion, we waited 30-60 min without moving the probes to reduce slow drift and provide more stable recordings. We isolated single and multiunit threshold crossings (three times RMS noise) by amplifying the site signals filtered between 750 Hz and 7.5 Khz (Cerebus, Blackrock microsystems). During recordings, animals were awake and passively viewing a gray screen. To keep animals awake and alert, animals were water-scheduled (Histed and Maunsell, 2014), and a 1 ml water reward was randomly provided on 5% of the stimulus trials; we verified animals were licking in response to rewards during the experiments. Superficial units were those recorded from a site within 400 mm of the cortical surface, identified by monitoring spike and local field potential activity (LFP) on sites at different depths; we typically observed desynchronized cortical LFP activity on any site below the surface and found the first substantial unit activity 100 mm below the surface. Deep units were those recorded 500-800 mm below the surface. Optogenetic stimuli were square light pulses with 2 ms linear ramps at start and end to reduce recording artifacts. Pulses were on for 600 ms and off for 1000 ms at a range of power levels (0.3-10 mW), presented in random order, with 100 repetitions per power level. For population plots, we normalized the light intensity in each experiment to control for fluctuations across experiments due to e.g. changes in dural thickness or tissue light absorption. We found a minimum value of the inhibitory responses (L 0 , with one value for each experiment, used for each unit recorded in that experiment) by fitting a 2-segment piecewise-linear function to the average inhibitory response (Figure 2-figure supplement 6). To avoid biases due to light attenuation at deep sites, L 0 was computed for each experiment on the superficial sites only, and this value was used for the deep site data (Figure 7). Population firing rate std. dev. and SEM (Figure 3) were calculated on the sum of all unit counts in a given time bin.

Pharmacological blocking of excitatory and inhibitory synapses
To classify cells using pharmacological weakening of cortical synapses, we divided the experiment into three phases, each 30-45 min in duration: first no blockers, then excitatory blockers, then excitatory plus inhibitory blockers. For each section, we delivered the same optogenetic stimulation protocol (above). To apply the pharmacological agents, we made a hole near the recording electrodes in the agarose on top of the brain, removed the normal saline covering the agarose and hole, and replaced it with a solution containing the agent in saline. After 15 min, the initial application of blocking solution was removed via aspiration and refreshed. We waited a total of 20 min for the pharmacological agent(s) to take effect before recording. Excitatory synaptic blockers (2 mM CNQX, 6 mM APV) affected AMPA, kainate, and NMDA synapses. Inhibitory blockers (1 mM bicuculline) affected GABA-A synapses. Source: Sigma-Aldrich (#C239, #A8054, and #14343).

Histology and electrophysiology probe tracking
To determine the location and depth of recording electrodes, we coated the shanks with 1,1'-Dioctadecyl-3,3,3',3'-tetramethylindocarbocyanine perchlorate (DiI) (Thermo Fisher #D282, 50 mg/mL solution in ethanol), prior to insertion into the brain. We dipped the electrode tips 5 times into the DiI solution and allowed the coating to dry for 30 s between immersions. Fluorescent electrode tracks were then visualized in coronal sections of fixed brain tissues. Mice were anesthetized with isoflurane and injected intraperitoneally with pentobarbital sodium (150 mg/kg). They were then perfused transcardially with cold (4˚C) PBS followed by cold 4% paraformaldehyde. Brains were extracted and fixed in 4% paraformaldehyde for 6-12 hr and then cryoprotected in a 30% (w/v) sucrose solution in PBS until they sank. Brains were sectioned at 50 mm on a freezing microtome (Leica), mounted on glass slides, and coverslipped with mounting media containing DAPI (Fluoromount-G with DAPI, Electron Microscopy Sciences). Slides were imaged using an Olympus slidescanner (Olympus BX61VS, Japan).

Counting opsin-expressing cells
To quantify opsin expression differences between viral (PV-Cre,AAV-Chronos-YFP) and transgenic (PV-Cre::ReaChR-mCitrine) animals, we used a confocal microscope (Zeiss LSM780) to image fluorescent neurons in L2/3 using 50 mm thick coronal brain tissue sections. (We used the same Cre line for both cases, so that any Cre expression patterns would affect neurons in both groups, and >90% of cortical PV neurons are reported to express Cre in this line [Hippenmeyer et al., 2005]). For each animal, we chose four 100 mm x 100 mm areas in layer 2/3, and constructed 3D stacks by aligning images across adjacent sections. 3D volumes from 3 viral and three transgenic animals were anonymized and manually counted by two independent observers. Observers were instructed to count cells once with a high threshold for accepting a cell as positive, and once with a low threshold for accepting a cell. Percentages of opsin-expressing cells were calculated as the average number of fluorescently labeled cells divided by the number of DAPI-labeled nuclei.

Anesthesia during electrophysiology
To test the effects of lowering overall activity on the stability of the cortical network we fitted animals with an isoflurane inhalation mask system (model V-1, VetEquip Inc) and provided anesthesia (0.25-0.5% isoflurane in 100% O 2 ) to put the animal into a lightly anesthetized but awake state during recordings. At the lower 0.25% isoflurane concentration, recordings in V1 displayed no synchronized oscillatory Up/Down state activities and the animals eyes were open throughout. The anesthesia experiments were divided into four sections: awake animal, anesthesia with no blockers, anesthesia with excitatory blockers, anesthesia with inhibitory blockers. Anesthesia was delivered for 20 min before starting recording.

Spike sorting
Spike waveforms were sorted after the experiment using OfflineSorter (Plexon, Inc). Single units were identified as waveform clusters that showed clear and stable separation from noise and other clusters, unimodal width distributions, and inter-spike interval histograms consistent with cortical neuron absolute and relative refractory periods. Multiunits were clusters that were distinct from noise but did not meet one or more of those criteria, and thus these multiunits likely group together a small number of single neurons. Signal-to-noise ratio (SNR) (Kelly et al., 2007;Histed, 2018) of single unit waveforms (median ± std): in visual cortex, 3.6 ± 1.7 (N = 167), see Figure 2-figure supplement 2; in other datasets: motor, 3.7 ± 1.2 (N = 103); somatosensory, 3.7 ± 1.7; in V1 PV data ( Figure 6) 2.91 ± 1.3. Only single-units were analyzed; multiunits were discarded. We also repeated our single-unit analyses using both more-and less-stringent criteria to qualify a unit as a single unit, and found no qualitative differences in the results (Figure 2-figure supplement 2). Supporting the idea that single units did not group together multiple units, we found no significant correlation between SNR and baseline firing rate in any of the four datasets (all p>0:2 by linear regression).

Model equations
To quantitatively analyze the network response, we used a two population rate model (Wilson and Cowan, 1972), in which the average firing rates of excitatory (E) and inhibitory (I) populations (r E and r I , respectively) evolve according to where f A , I AX and t X are the static transfer function (f-I curve), external input and time constant of population A (=E,I), respectively, while W AB is the strength of connections from population B to A. The optogenetic stimulation is described by the parameters l and L which represent the efficacy and the intensity of the stimulation light. Results shown in the main text ( Figure 4) have been obtained with a rectified-linear transfer function which is zero for input x smaller than the threshold x 0A , and increases linearly, with a gain a A , otherwise. This transfer function is the simplest one that can describe the data, which shows an approximately piece-wise linear dependence of firing rates on stimulation intensity (see Figure 4). We also fit the data using a different transfer function that smoothes the threshold non-linearity of the rectified linear function, using an additional parameter b A that controls the width of the exponential region around threshold. This function reduces to the rectified-linear transfer function when b A ! 0. We find that the nonlinear transfer function provides a minor improvement in describing the data and does not significantly affect the values of the inferred parameters (Figure 4-figure supplement 1).
In the experiments, recordings of the response are done in three separate phases: (1) A 'normal' phase in which recurrent interactions are intact; (2) A phase with E synaptic blockers; and (3) A phase with both E and I synaptic blockers. The addition of blockers in the model is described with two parameters, E;I 2 0; 1 ½ , representing the decrease in strength of excitatory and inhibitory synapses. After the addition of excitatory blockers, connectivity is modified as W EE ; I EX ; W IE ; I IX ! E W EE ; E I EX ; E W IE ; E I IX ; (4) and with inhibitory blockers, connectivity is modified as Combining Equations 1, 4 and 5 with the transfer function of Equation 2, the model uses 15 parameters (2 time constants, 2 thresholds, 2 gains, 4 connectivity strengths, 2 external inputs, 1 stimulation efficacy, 2 blocker efficacies) to describe simultaneously the three phases of the experiment.
To compare the model response with the data, we find the equilibrium solution of Equation 1, shown in Equation 6. Since t E and t I do not affect the equilibrium solution, the number of relevant independent variables in the model reduces to 13. Moreover, the values of a E and a I can be reabsorbed in the definitions of W, I and l. It follows that all the parameters can be inferred up to a proportionality constant and, without loss of generality, we can fix a E ¼ a I ¼ 1. This reduces the number of independent parameters to 11.
Note that this number of parameters is smaller than the number of parameters needed for piecewise-linear fits of the data (average firing rates vs light intensity for both E and I neurons). Such piecewise-linear fits of the data need 15 parameters: 3 phases times five parameters per phase (three for inhibitory neurons -two per linear region, minus one for the continuity constraint; and two for excitatory neurons, for the single linear region at low intensities). Thus, the model does not have enough parameters to fit successfully the data, unless the model structure accurately describes the cortical responses. Indeed, Figure 4 shows that the model gives a good description of the data. The additional 4 degrees of freedom present in the data lead to four parameter free model predictions described in the main text.
The fixed point solutions of Equation and For every value of the laser intensity, the excitatory and inhibitory nullclines are defined as the functions r I ¼ r I ðr E Þ which solve the first and the second line of Equation 6, respectively. Every intersection between the two nullclines is a solution of Equation 6 and gives a stationary state of the network dynamics.
Equation 6 has two dimensions that can be inferred only up to a multiplicative constant: ; g E W EI ; g E I EX ; g E x 0;E ; W IE ; W II þ 1; I IX ; x 0;I ; l ! g I W IE ; g I ðW II þ 1Þ; g I I IX ; g I x 0;I ; g I l : Because of these invariances, using data from a single experimental phase, only ratios between parameters can be inferred. This is no longer true when the three phases are considered; in this case all the 11 relevant model parameters can be found.

Parameter inference
Given the dataset of excitatory and inhibitory responses, the best set of parameters describing the data is found via global optimization as follows: 1. Select random initial parameters in the interval [0,10]. 2. Find the optimal set of parameters through a least-squares optimization (python function 'curve_fit', variables are constrained to be positive and, for E;I , in the interval [0,1]) of the difference between observed rates and model predictions given by Equation 6. 3. Repeate the procedure 10 4 times, and select the solution with minimal error as the optimal parameter set describing the data.
The optimization procedure is applied simultaneously on data from all phases and recording sessions in V1. For the first phase (awake), data from all recording sessions from VGAT-ChR2 animals in V1 (9 days, four animals) were pooled together for fitting. For the second and third phases (E blockers, E+I blockers), in order to describe the variability from day to day of the efficacy of the blockers, we use a different set of E;I for each recording session (6 days, two animals). For these two phases, three out of the nine recording sessions are excluded from the analysis because blockers were added after anesthetizing the animals. The optimal parameters found by this approach are: These parameters are used in Equation 6 to generate the best model description of the data in the first phase of the experiment ( Figure 4A, first column, black line). For the second and third phases, predictions also depend on the efficacy of synaptic blockers, which are different in the different recording sessions (see Figure 4-figure supplement 2). Data and model predictions for all sessions are shown in Figure 4-figure supplement 2; in the main text, we showed averages of the two computed across sessions for each phase ( Figure 4A, second and third columns, black line).
We use a bootstrap approach, combined with the global optimization described above, to estimate the precision with which model parameters can be inferred from the data. The optimal (minimum-error) estimates obtained from 10 4 random resamplings (random with replacement) are shown in Figure 4-figure supplement 2. Despite the non-uniqueness of the solutions, the bootstrap shows (dashed lines in Figure 4: 1 s.d. via boostrap) that model rates are clustered around the data in all three phases of the experiment.

Stability of the solution
To analyze the stability of the system, we compute the eigenvalues of the dynamical matrix defined by Equation 1. For the linear transfer function we find with r ¼ t E aI WIE t I aE WEI This ratio cannot be inferred from the data since t E;I do not appear in the static solution. We thus determined t E and t I (shown in Figure 4E) from the response dynamics ( Figure 4D).

With multiple inhibitory subclasses, ISN operation need not imply paradoxical suppression
To understand how the network response to PV stimulations depends on its structure, here we analyze a three-population model ( Figure 6A) with one excitatory and two inhibitory populations, referred hereafter as E (pyramidal cells), P (ChRonos-expressing PV cells), and I (non-expressing PV , SOM, VIP. . .). The network dynamics is described by t P drP dt ¼ Àr P þ f P W PE r E À W PI r I À W PP r P þ I PX þ lL ð Þ: The model predictions in Figure 6 have been obtained from Equation 10 where frac 2 ð0; 1Þ represents the fraction of ChRonos-expressing PV cells, W VGAT AB and I VGAT A represent the parameters of Equation 8, obtained using the V1 data from the VGAT-ChR2 mouse line.
Using the rectified linear transfer function of Equation 2, and assuming a A ¼ 1 to simplify expressions, the equilibrium response to stimulation of the P population is given by shows that, both for W EE <1 and W EE >1, depending on the network connectivity, the ratio on the r.h.s. can be either positive or negative. It follows that, in a network with two inhibitory populations, an unstable excitatory subnetwork does not imply paradoxical suppression when one of the populations is stimulated. This can be seen more easily in a simplified model in which connectivity depends only on the identity of the presynaptic unit. In particular, using In the large W limit, in which the excitatory population is unstable without inhibition, the network shows the paradoxical effect only if k I <1 and k P >1 À k I . (Note that, for the network to be stable in this condition, the matrix of the coefficients of Equation 10, linearized around the fixed point, must have eigenvalues with a negative real part, a requirement which also involves the magnitude of t E , t I , and t P .)

Mathematical description of STP
To investigate the effects of STP on the relationship between paradoxical effect and inhibitory stabilization, we use the phenomenological description of STP developed in Markram et al., 1998;Tsodyks et al., 1998. This framework is accurate enough to capture quantitatively both short term depression (STD) and short term facilitation (STF), and is simple enough to allow analytical investigation. The state of each synapse is described by two variables: the fraction x (0 x 1) of vesicles available for release, and the fraction u (0 u 1) of available vesicles that release neurotransmitter after a presynaptic spike (release probability). In the model, a presynaptic spike arriving at time t opens calcium channels in the presynaptic terminal, which generates an increases the value of uðt À Þ by an amount U½1 À uðt À Þ (0<U<1 is a fixed parameter) and produces the release of a fraction uðt þ Þxðt À Þ of vesicles. After the spike, channels close with a time constant t F and vesicles recover with a time constant t D . As shown in Tsodyks et al., 1998, when presynaptic spikes are Poissondistributed with rate r, the dynamics of x and u is given by In this model, synaptic efficacy and postsynaptic current are proportional to Wx and Wxr, respectively. The modulation of synaptic efficacy by the presynaptic rate can be understood as follows. When u is fixed, x and the synaptic efficacy decrease with r, leading to synaptic depression. When x is fixed, the increase in u with r augments synaptic efficacy and leads to synaptic facilitation. When both u and x are allowed to vary, facilitation dominates at low r and depression dominates at high r Tsodyks et al., 1998). A stationary input rate rðtÞ ¼ r 0 gives and, since U<1, produces a synaptic current which increases monotonically with r 0 ; this monotonic increase will be important in the upcoming sections.
In what follows, we include the description of STP discussed above in our network model (Equation 1 of the main text) by adding variables x AB and u AB for all projections W AB , so that the synaptic current W AB r B becomes x AB u AB W AB r B . Therefore, the network dynamics is given by Here A; B 2 ½E; I and, to simplify the analysis, we assumed threshold-linear single neuron transfer function.

Conditions for paradoxical response in networks with STP
In this section, we investigate conditions for paradoxical inhibitory response in networks with STP; the role of inhibition in stabilizing activity is discussed in the next section.
As discussed in Ozeki et al., 2009 and in Figure 1-figure supplement 1, paradoxical inhibitory response can be studied using the r E and r I nullclines of the model in the r E =r I plane. The difference with those previous studies is that we have additional dynamical variables for synaptic strengths, and these should be set to their equilibrium values. Using dx AB =dt ¼ 0 and du AB =dt ¼ 0, we express W AB and x AB as a function of r B (as done in Equation A2) and define f E ðr I Þ ¼ x EI ðr I ÞW EI ðr I ÞW EI r I ; f I ðr I Þ ¼ r I þ x II ðr I ÞW II ðr I ÞW II r I : Since each term x AB ðr B Þu AB ðr B ÞW AB r B is a monotonically increasing function of the presynaptic rate r B , f E;I can be inverted and, in the region of rates r E;I >0, we can write the excitatory and inhibitory nullclines as Using the fact that f À1 I ½: is an increasing function of its argument, we obtain that the inhibitory nullcline r inhibitory I ðr E ; LÞ increases with r E and L. It follows that the response of the inhibitory population is paradoxical, that is the stationary value of r I decreases with L, only if the excitatory nullcline r excitatory I ðr E Þ has positive slope, that is if Equation A6 must be satisfied by any network of excitatory and inhibitory neurons with STP which shows paradoxical effect. In the next section, we show that this condition is violated in networks which are stable with fixed inhibition.
To conclude, we note that the condition of Equation A6 is only necessary but not sufficient, and in order for the paradoxical response to emerge, one also needs df À1 E ðr E Þ=dr E <df À1 I ðr E Þ=dr E . Without STP, this condition is met any time the network is dynamically stable (Ozeki et al., 2009). With STP, proving an analogous result is more complicated, as it involves the stability of Equation A3, that is of a ten-dimensional dynamical system.

Inhibition stabilization in networks with STP
In this section, we investigate the relation between paradoxical effect and inhibitory stabilization in networks with STP. As discussed above, plasticity can be included in all the synapses of the model (i.e. at all projections W AB in Equation 1) but only STP at recurrent excitatory synapses (W EE ) can influence the role of inhibition in the stabilization of dynamics. In fact, a network is said to be inhibition stabilized if instabilities emerge once inhibitory rates are forced to be constant in time Ozeki et al., 2009). This condition is not affected by plasticity involving inhibitory neurons (both pre-and post-synaptically), therefore we can limit our investigation to STP in W EE . In what follows, first we study analytically networks in which only one between STD or STF is present, then we investigate numerically networks in which both STD and STF are included simultaneously. In all these cases, we find that paradoxical inhibitory response implies that the network is unstable with fixed inhibition. Therefore, if a stable network shows paradoxical response, the excitatory network must be unstable on its own and has to be stabilized by inhibition. It follows that the conclusions derived in the main text about cortical networks remain valid also when STP is taken into account in the analysis.

Effects of STD
With STD in W EE synapses, the network dynamics is given by The effect of STD on the stationary response can be seen by looking at the r E and r I nullclines of the system (Appendix 1- figure 1A) in the r E =r I plane, setting x to its steadystate value. Expressing x as a function of r E , we find The above expression shows that STD decreases the self-excitation of excitatory cells as the excitatory rate increases. The stationary solutions of the Equation A7 (here referred to as x Ã , r Ã E , r Ã I ) are the points at which the two nullclines of Equation A8 cross. As discussed in general terms in the previous section, increasing L in Equation A8 moves the inhibitory nullcline upward and Equation A6 provides a necessary condition for the paradoxical inhibitory response to emerge; this condition is met any time that is not satisfied, the network response is not paradoxical. The stability of the stationary solutions is found with a standard perturbation analysis around the stationary state. With frozen inhibition, we can write perturbation around the stationary solutions as x ¼ x Ã þ d x and r E ¼ r Ã E þ d E . In the case of r Ã E >0, the dynamics of these perturbations is given by Equation A10 shows that, with fixed inhibition and fixed STD (d x ¼ 0), the excitatory network is unstable if which, since 0<x<1, also implies that W EE >1. With fixed inhibition and dynamic STD, the dynamics is stable if If Equation A12 is verified, the network is stable without dynamic inhibition. Since this condition cannot be satisfied simultaneously with Equation A9, stabilization in a network with STD in W EE synapses which shows paradoxical effect has to be realized by inhibition.
Appendix 1-figure 1. Network dynamics with short term depression. (A) Excitatory (red) and inhibitory (blue) nullclines of Equation A7. Different columns correspond to different values of W IE (from left to right, W IE ¼ 11:10, 3.42, 0.60); symbols correspond to the stationary solutions of the Equation A7. STD bends the excitatory nullcline downward, with a decrease that starts when r E~1 =Ut D . Since increasing L moves the inhibitory nullcline to the left, paradoxical inhibitory response emerges only if the nullclines meet in the region of positive slope of the excitatory nullcline. (B) As in (A) but for Excitatory (red) and x (black) nullclines of Equation A7. Dashed lines correspond to threshold for paradoxical response (purple, Equation A9) and stability of the excitatory population with fixed x and r I (green, Equation A12). Note that paradoxical response can appear only if the network with fixed r I is unstable. Therefore, if a network is stabilized by STD, it cannot show paradoxical response. (C) Dynamics of excitatory neurons in response to perturbations applied at time t ¼ 0. Different colors correspond to different constraints on the dynamics. Different columns are as in panels A and B. (D) Values of x Ã as a function of W EE for which the network: has paradoxical response (Equation A9, purple), is stable with fixed r I (Equation A12, green), is stable with fixed x and r I (Equation A11, red). Paradoxical inhibitory response emerges only if the excitatory subnetwork is unstable on its own (W EE >1) and with STD and fixed inhibition (Equation A12). Symbols correspond to the stationary solutions of panels A and B. Simulation parameters are: W EE ¼ 3:33, W EI ¼ 1:77, W II ¼ 7:11, I EX ¼ 8:51, I IX ¼ 32:45, l ¼ 6:3, x 0E ¼ 1:19, x 0I ¼ 8:65, t E ¼ 7:8 ms, t I ¼ 34:3 ms, t D = 500 ms, U = 0.6.

Effects of STF
In the presence of STF in W EE synapses, the network dynamics is given by Proceeding as in the case of STD, we find that inhibitory response is paradoxical, that is Equation A6 is verified, if On the other hand, dynamics with fixed inhibition is stable if Therefore, as in the case of STD, stability with frozen inhibition is incompatible with paradoxical effect and paradoxical inhibitory response implies inhibitory stabilization.

Analysis of networks with STF and STP in E!E synapses
To conclude our analysis, we investigate the relation between paradoxical inhibitory response and inhibitory stabilization in networks with both STF and STD in W EE synapses. In these networks, the dynamics is given by t E drE dt ¼ Àr E þ x W W EE r E À W EI r I þ I EX À x E0 ½ þ t I drI dt ¼ Àr I þ W IE r E À W II r I þ I IX À x I0 þ lr I ½ þ du dt The analytical approach we used to study networks with either STD or STF cannot be applied here, as finding the stationary solutions of Equation A16 requires solving a system of equations with polynomials of degree higher than two in r E . Therefore, to investigate the relation between paradoxical inhibitory response and inhibitory stabilization, we used a numerical approach.
We considered networks with randomly generated parameters and, for each realization, we analyzed numerically the response to external excitation of inhibitory neurons and the role of inhibition in stabilizing the dynamics. We considered networks with t E ¼ t I ¼ 10 ms, l ¼ 1, x 0E ¼ x 0I ¼ 0 and randomly generated W AB ; I AX 2 f0; 10g, t F;D 2 f0:01s; 0:99sg and U 2 f0:01; 0:99g, using uniform distributions in the specified intervals. For each network, that is for each realization of the randomly-generated parameters, we simulated the dynamics induced by Equation A16 for a time time T ¼ 10 s (Euler method, integration time step dt ¼ 0:1 ms) and measured the stationary values of r E , r I , u and x. We assessed numerically that each solution was stationary, and that rates were not oscillating or running away, by measuring the standard deviation of r E for t>5 s; we considered the dynamics stationary if, for t>5 s, std ðr E Þ= mean ðr E Þ<10 À3 . For each network with stationary dynamics, we ran another simulation (T ¼ 10 s) and measured the response of inhibitory neurons to external stimulation. Specifically, we increased the value of L from to 0.1 for t 2 f5; 10g s and considered the response to be paradoxical if ðr I ðL ¼ 0Þ À r I ðL ¼ 0:1ÞÞ=ðr I ðL ¼ 0Þ þ r I ðL ¼ 0:1ÞÞ>10 À3 , where r I ðL ¼ zÞ represents the rate measured at the end of the simulation period (t ¼ T ¼ 10 s) for L ¼ z. To investigate the role of inhibition in stabilizing the dynamics, we ran yet another simulation (T ¼ 10 s) in which, for time t>5 s, we clamped the value of inhibitory rates and applied a perturbation to excitatory neurons (increase in 0.05 spk/s of r E for t>5:00 s, t<5:01 s). We considered the network to be inhibition-stabilized if jr E ðP ¼ 0:05 spk=sÞ À r E ðP ¼ 0 spk=sÞj=ðr E ðP ¼ 0:05spk=sÞ þ r E ðP ¼ 0spk=sÞÞ>10 À3 , where r E ðP ¼ zÞ represents the rate measured at the end of the simulation period (t ¼ T ¼ 10 s) with a perturbation of z spk/s. Applying this approach to N ¼ 998845 networks, we found that 98% (979227/998845) of them had stationary dynamics. Of these, 65.5% (639245/979227) were inhibition stabilized and 9% (88559/979227) had paradoxical inhibitory response. Of the inhibition stabilized networks, only 14% (88559/639245) showed paradoxical inhibitory response. This is consistent with our analytical results in networks with either STD or STP, where we found that inhibitory stabilization does not imply paradoxical inhibitory response. Of the networks with paradoxical response, consistent with what we found in networks with either STD or STP, 100% (88559/88559) were inhibition stabilized.
Our numerical results suggest that, also in networks in which STD and STF coexist in W EE synapses, paradoxical inhibitory response implies inhibitory stabilization.
where the integral is evaluated at the inferred baseline activity. The quantity W EE is the generalization of the parameter W EE analyzed in the main text, which takes into account the amplification of coupling strength due to transfer function nonlinearities. Values of W EE inferred with data bootstrap are shown in Appendix 2- figure 1D. Results are consistent with W EE larger than one, indicating that indeed the excitatory subnetwork is unstable without inhibition, and in the range 5-15, consistent with the result obtained in the main text using a linear model.
To summarize, the analysis derived in this section shows that results obtained in the main text are preserved when a more biologically accurate, spiking, model is used to describe the data. Moreover, this approach also shows that inferred parameters are consistent with known biophysical properties of cortical networks.