Mechanisms underlying the response of mouse cortical networks to optogenetic manipulation
Abstract
GABAergic interneurons can be subdivided into three subclasses: parvalbumin positive (PV), somatostatin positive (SOM) and serotonin positive neurons. With principal cells (PCs) they form complex networks. We examine PCs and PV responses in mouse anterior lateral motor cortex (ALM) and barrel cortex (S1) upon PV photostimulation in vivo. In ALM layer five and S1, the PV response is paradoxical: photoexcitation reduces their activity. This is not the case in ALM layer 2/3. We combine analytical calculations and numerical simulations to investigate how these results constrain the architecture. Twopopulation models cannot explain the results. Fourpopulation networks with V1like architecture account for the data in ALM layer 2/3 and layer 5. Our data in S1 can be explained if SOM neurons receive inputs only from PCs and PV neurons. In both fourpopulation models, the paradoxical effect implies not too strong recurrent excitation. It is not evidence for stabilization by inhibition.
Introduction
Local cortical circuits comprise several subclasses of GABAergic interneurons which together with the excitatory neurons form complex recurrent networks (Goldberg et al., 2004; Jiang et al., 2015; Karnani et al., 2016; Markram et al., 2004; Moore et al., 2010; Pfeffer et al., 2013; Tasic et al., 2018; Tremblay et al., 2016). The architecture of these networks depends on the cortical area and layer (Beierlein et al., 2003; Jiang et al., 2013; Rudy et al., 2011; Xu et al., 2013; Xu and Callaway, 2009).
Optogenetics is now classically used to reversibly inactivate a particular cortical area or neuronal population to get insights into their functions (Atallah et al., 2012; Guo et al., 2014b; Lee et al., 2012; Li et al., 2015; Svoboda and Li, 2018). Optogenetics has also been applied to isolate the different components (e.g. feedforward vs. recurrent) of the net input into cortical neurons (Lien and Scanziani, 2018; Lien and Scanziani, 2013). It can also be used to experimentally probe the architecture of local cortical circuits (Moore et al., 2018; Xu et al., 2013). However, because of the complexity of these networks and of their nonlinear dynamics, qualitative intuition and simple reasoning (e.g. ‘boxandarrow’ diagrams) are of limited use to interpret the results of these manipulations.
‘Paradoxical effect’ designates the phenomenon that stimulation of a GABAergic interneuron population not only decreases the average activity of the principal cells (PCs) but also decreases the activity of the stimulated population (Murphy and Miller, 2009; Ozeki et al., 2009; Tsodyks et al., 1997). Intuitively, paradoxical effect arises when the stimulation induces a strong activity suppression in the PCs (Kato et al., 2017; Moore et al., 2018), such that the overall (synaptic+stimulus) excitation to the stimulated population decreases. However, the precise conditions under which the paradoxical effect occurs are difficult to establish without mathematical modeling.
In simple models consisting of only two populations (one excitatory and one inhibitory) these conditions have been mathematically derived. The paradoxical effect occurs when the networks operates in the regime known as inhibition stabilized (inhibition stabilized networks, ISN) in which the total the total recurrent excitation is so strong that inhibition is necessary to prevent a blow up in the activity (Murphy and Miller, 2009; Ozeki et al., 2009; Tsodyks et al., 1997). Networks, with several inhibitory populations have been recently investigated (Garcia Del Molino et al., 2017; LitwinKumar et al., 2016; Sadeh et al., 2017). These studies considered network models with synaptic currents small compared to neuronal rheobase currents (Gerstner et al., 2014; Lapicque, 1909). However, interactions in cortex are stronger than what is assumed in these studies (Shadlen and Newsome, 1994).
Simple networks with strong interactions comprising one excitatory and one inhibitory population have been studied extensively. In a broad parameter range not requiring finetuning, such networks dynamically evolve into a state in which strong excitation is balanced by strong inhibition such that the net input into the neurons is comparable to their rheobases (van Vreeswijk and Sompolinsky, 1998; van Vreeswijk and Sompolinsky, 1996). The theory of balanced networks has been developed for a variety of single neuronal models including binary neurons (van Vreeswijk and Sompolinsky, 1998; van Vreeswijk and Sompolinsky, 1996), rate models (Harish and Hansel, 2015; Kadmon and Sompolinsky, 2015), leakyintegrateand fire neurons (Hansel and Mato, 2013; Mongillo et al., 2012; Rosenbaum and Doiron, 2014; Roxin et al., 2011; Van Vreeswijk and Sompolinsky, 2005) and conductancebased models (Hansel and van Vreeswijk, 2012; Pattadkal et al., 2018).
In the present study, we investigate experimentally the effects of the photostimulation of PV interneurons on the anterior lateral motor cortex (ALM) and barrel cortex (S1) of the mouse. We show that twopopulation network models do not suffice to account for these effects. To overcome this limitation, we develop a theory for the paradoxical effect in balanced networks that takes into account the multiplicity of GABAergic neuronal populations. Combining analytical calculations and numerical simulations, we study the responses of these networks at population and single neuron level. For twopopulation balanced networks it has been shown that the paradoxical effect only occurs when the network is inhibition stabilized (Pehlevan and Sompolinsky, 2014; Wolf et al., 2014). Here we show that in contrast, in fourpopulation networks, the paradoxical effect can occur even if the network is not inhibition stabilized. We conclude with prescriptions for experiments that according to the theory can be informative about network architectures in cortex.
Results
ALM layer 5 and S1 exhibit paradoxical effect but not ALM layer 2/3
We expressed a redshifted channelrhodopsin (ReaChR) in PV interneurons to optogenetically drive local inhibition in the barrel cortex (S1) and anterior lateral motor cortex (ALM) of awake mice (Hooks et al., 2015). We used orange light (594 nm) to illuminate a large area of ALM or S1 (2 mm diameter), photostimulating a large proportion of PV interneurons (Figure 1A). We measured the lightinduced effects on neural activity using silicon probe recordings. In both brain areas, putative PCs and putative PV neurons were identified based on spike width (Methods). Neurons with wide spikes were likely mostly PCs. Units with narrow spikes were fast spiking (FS) neurons and likely expressed parvalbumin (Cardin et al., 2009; Guo et al., 2014b; Olsen et al., 2012; Resulaj et al., 2018). We investigated the responses of these neurons as a function of the photostimulation intensity in ALM layer 2/3 and layer 5, and in S1.
We found that in all recorded layers and areas, the population average activity of the PCs decreased with the optogenetic drive (Figure 1B, Figure 2). In contrast in ALM, the PV population exhibited a behavior which depended on the recorded layer.
In ALM layer 2/3, the population average firing rate of PV neurons monotonically increased with the photostimulation intensity. However, individual neuron responses were heterogeneous. Most PV neurons increased their spike rates from baseline with increased photostimulation intensity. Some PV neurons initially decreased their spike rates below baseline for low light intensity.
In ALM layer 5, the response of the PV population was nonmonotonic. For low laser intensity, their activity paradoxically decreased with the optogenetic drive. The slope of the normalized firing rate v.s. laser intensity was significantly different from zero for both the PC and PV populations (Figure 1F). The ratio of their slopes was 0.62 ± 0.28. At high photostimulation intensity, the activity of the PV population increased. At intermediate photostimulation intensity (0.5 mW/mm^{2}), the response of the PV neurons was significantly different between layer 2/3 and layer 5 (Figure 1E, p<0.005, unpaired ttest, twotailed test).
Paradoxical decrease in PV neurons activity with the optogenetic drive was also observed in S1. Remarkably, the concomitant decrease of the PC and the PV population activities was proportional (Figure 1G, ratio of slopes PV/PC, mean ± SEM; S1, 1 ± 0.29).
In both ALM layer 5 and S1, there was also a large diversity of responses. Most PV neurons decreased their activity at low photostimulation intensity. At high laser intensity (5 mW/mm^{2}), a fraction of PV neurons (6/12 in ALM layer 5 and 6/10 in S1) had a larger response than baseline, while the rest remained suppressed. Figure 2 shows the spike rates of PCs and PV neurons at an intermediate light intensity (0.5 mW.mm^{2}).
Network models
To assess the network mechanisms which may account for the experimental data from ALM and S1, we first considered models consisting of one excitatory and one inhibitory population. Since it is well established that cortical circuits involve a variety of inhibitory subpopulations, we later extended the theory to network models of four populations of neurons representing PCs and three subtypes of GABAergic interneurons in cortex. In all our models, neurons are described as integrateandfire elements. The data we seek to account for, were obtained in optogenetic experiments in which the laser diameter was substantially larger than the spatial range of neuronal interactions and comparable to the size of the area in which activity was recorded. Therefore, in all our models, we assume for simplicity that the connectivity is unstructured. We modeled the ReachRoptogenetic stimulation of the PV population as an additional external input, ${I}_{opto}$ , into PV neurons. We assumed that it depends on the intensity of the laser, ${\Gamma}_{opto}$ , as ${I}_{opto}={I}_{0}log\left(1+\frac{{\Gamma}_{opto}}{{\Gamma}_{0}}\right)$ where ${I}_{0}$ and ${\Gamma}_{0}$ are parameters (Figure 3—figure supplement 1; Hooks et al., 2015).
Twopopulation model
The twopopulation network is depicted in Figure 3A. It is characterized by four recurrent interaction parameters, ${J}_{\alpha \beta}$ , and two feedforward interaction parameters, ${J}_{\alpha 0}$ , $\alpha ,\beta \in \{E,I\}$ (see Materials and methods).
Results from numerical simulations of the model are depicted in Figure 3B and C where, the dependence of the population activities normalized to baseline, are plotted against the intensity of the laser, ${\Gamma}_{opto}$ . Figure 3B shows the response of the network where the recurrent excitation, J_{EE} , is non zero. The activity of the PV population, r_{1} varies nonmonotonically with the laser intensity. For small intensities, r_{1} paradoxically decreases together with the activity of the PCs, r_{E} . This paradoxical effect stems from the fact that the decrease in the activity of the PCs yields a reduction in the excitation to PV neurons which is not compensated for by the optogenetic drive. As a result, the net excitation to PV neurons diminishes yielding a decrease in r_{I} . When r_{E} becomes very small, this mechanism does not operate anymore and consequently, r_{I} increases as ${\Gamma}_{opto}$ is increased further. In Figure 3C, J_{EE} is zero, r_{I} monotonically increases with the light intensity whereas r_{E} monotonically decreases. For small intensities, r_{I} is close to a constant. It starts to increase appreciably only when ${r}_{E}\simeq 0$ . Therefore, the PV response is not paradoxical.
Qualitatively this model seems to account for our experimental data from ALM layer 2/3, ALM layer 5 and S1. It would imply that in layer 5, J_{EE} is sufficiently large to generate the paradoxical effect, while in layer 2/3 this is not the case. On closer inspection however, there are major discrepancies between the simulation results and the experimental data. In our recordings in both ALM layer 5 and S1, the PV population activity reaches a minimum while the PCs are still significantly active: relative to baseline the activity is 40% in ALM and 25% in S1. In contrast, in the twopopulation model, the minimum of the PV activity is reached (Appendix 1B) when excitatory neurons are virtually completely silenced (Figure 3B, Figure 3—figure supplement 2A). In fact one can show that for sufficiently large K, when r_{I} is minimum, the activity of the excitatory population is exponentially small in K. As a result, to account for the data one needs to assume that $K\simeq 10$ .
In addition, in the experimental data the activities of the PC and PV populations in S1 decrease in equal proportions before the minimum of the PV activity (Figure 1B). This cannot be accounted for in a twopopulation model unless parameters are finetuned (Figure 3—figure supplement 3). Analytical calculations (Appendix 1B) supplemented with numerical simulations show that this proportional decrease only happens when the determinant of the interaction matrix, J _{αβ}, is close to zero. Moreover, the external input must also be finetuned so that the neurons have biologically realistic firing rates (Figure 3—figure supplement 3).
The experimental data from ALM layer 2/3 show that for already small light intensity the activity of PV neurons increases appreciably. This is in contrast with Figure 3C. In Figure 3—figure supplement 2B, we show that the twopopulation model can account for this feature only if the recurrent excitation is very weak in that layer and the connectivity is extremely sparse.
These discrepancies prompted us to investigate whether models with several populations of inhibitory neurons can account for our experimental data without finetuning. We focus on two fourpopulation network models. Both consist of three populations representing PCs, PV and SOM neurons and a fourth population representing other inhibitory neurons. The main difference between the two models lies in the inhibitory populations from which SOM neurons receive inputs.
A fourpopulation model with V1like architecture (Model 1)
We first investigated the dynamics of a fourpopulation network with an architecture that is similar to the one reported in layer 2/3 in V1 (Pfeffer et al., 2013) and S1 (Lee et al., 2013) (Figure 4A). The model consists of four populations representing PCs, PV, SOM and VIP neurons. SOM neurons do not interact with each other (Adesnik et al., 2012; Gibson et al., 1999; Hu et al., 2011). VIP neurons only project to the SOM population (Jiang et al., 2015; Pfeffer et al., 2013). All neurons except SOM receive inputs from sources external to the network (e.g. thalamus) (Beierlein et al., 2003; Beierlein et al., 2000; Cruikshank et al., 2010; Ma et al., 2006; Xu et al., 2013). The same architecture was considered in LitwinKumar et al. (2016).
Following Pfeffer et al. (2013), the PV population does not project to the SOM population. Other studies have reported such a connection (Jiang et al., 2015). However, adding such a connection to Model 1 does not qualitatively affect the PC and PV responses (see Appendix 1C).
We considered parameter sets such that: 1) At baseline, the network is operating in the balanced state with all populations active; 2) the activity of the PC population decreases with the laser intensity as observed in our experiments.
Theory in the large N, K limit
It is instructive to consider the limit in which the number of neurons in the network, N, and the average number of connections per neuron, K, go to infinity. In this limit, the analysis of the stationary state of the network simplifies (see Materials and methods). This stems from the fact that when interactions are numerous, excitatory and inhibitory inputs are strong and only populations for which excitation is balanced by inhibition have a finite and nonzero activity. The average activities of the four populations are then completely determined by four linear equations, the balance equations, which reflect this balance. Solving this system of equations yields the population activities, r_{α} , α = E, I, S, V, as a function of the external inputs to the network. In particular, when the laser intensity is sufficiently small, the four populations are active and their firing rates vary linearly with the current induced by the photostimulation (Appendix 1C).
Figure 4 plots the activities of the populations vs. the optogenetic input into PV neurons, I_{opto} , for two sets of interaction parameters. In Figure 4B, the activity of the PV population, r_{I} , increases with I_{opto} . In contrast, in Figure 4C, r_{I} decreases with I_{opto} : the response of the PV population is paradoxical.
To characterize for which interaction parameters the PV response is paradoxical, we consider the 4 × 4 susceptibility matrix $\left[{\chi}_{\alpha \beta}\right]$ . The element ${\chi}_{\alpha \beta}\left(\alpha ,\beta =E,I,S,V\right)$ is the derivative of the population activity, ${r}_{\alpha}$ , with respect to a small additional input, into population $\beta $ , ${I}_{\beta}$ . Evaluated for small ${I}_{\beta}$ , ${\chi}_{\alpha \beta}$ characterizes by how much r_{α} varies with an increasing but weak extra input into population β. Its sign indicates whether r_{α} increases or decreases with I _{β}. The elements of the susceptibility matrix can be decomposed in several terms corresponding to the contributions of different recurrent loops embedded in the network (Appendix 1C). Using this decomposition one can show whether the PV response is paradoxical or not depends on the interplay between two terms. One is the gain of the disinhibitory feedback loop PCVIPSOMPC and the other is the product of the recurrent excitation, J_{EE} , with the gain of the disinhibitory feedback loop VIPSOMVIP (Figure 4—figure supplement 1). Remarkably, PV neurons are not involved in these two terms. A straightforward calculation (Equation A37) then shows that the response of PV neurons increases with I_{opto} if the recurrent excitation is sufficiently strong, namely if
The denominator in ${J}_{EE}^{*}$ is the strength of the connection from the SOM population to the VIP population. The numerator is the gain of the pathway which connects these two populations via the PCs. When $J}_{EE}>{J}_{EE}^{\text{*}$ the negative contribution of the disinhibitory loop PCVIPSOMPC dominates in the expression of ${\chi}_{II}$ . It is the opposite when $J}_{EE}<{J}_{EE}^{\ast$ . The stability of the balanced state provides other necessary conditions that the interactions must satisfy (see Materials and methods). In particular, the determinant of the interaction matrix, J, must be positive.
The difference between the behaviors in Figure 4B and C can now be understood as follows: in Figure 4B, $J}_{EE}>{J}_{EE}^{\text{*}$ and ${\chi}_{II}=1.6>0$ , thus, ${r}_{I}$ increases with I_{opto} ; in Figure 4C, $J}_{EE}<{J}_{EE}^{\ast$ and ${\chi}_{II}=5.1<0$ and thus, r_{I} decreases. Remarkably, in both cases the activities of the PC and VIP populations normalized to baseline, are always equal (Figure 4B–C, right panel). This is a consequence of the balance of excitatory and inhibitory inputs into the SOM population which implies that r_{E} and r_{V} are proportional (see Materials and methods, Equation 19.3).
In Figure 4B, the activity of the SOM population decreases with the laser intensity. This also stems from the fact that $J}_{EE}>{J}_{EE}^{\text{*}$ (Appendix 1C, Equations A3134). This qualitative behavior is therefore independent of parameter sets, provided that inequality (1) is satisfied. In contrast, for parameters for which $J}_{EE}<{J}_{EE}^{\ast$ the activity of the SOM population either decreases or increases with I_{opto} depending on other parameters. Moreover, it is straightforward to prove that if $J}_{EE}>{J}_{EE}^{\text{*}$ , the product ${\chi}_{EI}{\chi}_{IE}$ is positive (Appendix 1C). Since we assumed that r_{E} decreases upon photostimulation of PV neurons, namely ${\chi}_{EI}<0$ , this implies that ${\chi}_{IE}$ is also negative. In other words, in Model 1, a nonparadoxical response of the PV population upon PV photostimulation implies that the PV activity decreases when PCs are photostimulated.
When I_{opto} is sufficiently large, the solution of the four balance equations will contain one or more populations for which r_{α} < 0. Obviously such a solution is inconsistent. Instead, other solutions should be considered where at least one population has a firing rate which is zero and the firing rates of the other populations is determined by a new system of linear equations with lower dimensions (see Materials and methods, Appendix 1C). Consistency requires that in these solutions the net input is hyperpolarizing for the populations with r_{α} = 0. As a consequence, the network population activities are in general piecewise linear in I_{opto} (Figure 4—figure supplement 2).
The large N, K analysis provides precious insights into the dynamics of networks with reasonable size and connectivity. In particular, we will show that the criterion for the paradoxical effect, Equation 1, remains valid up to small corrections. Although it is possible to treat analytically the dependence of r_{α} on I_{opto} for finite K, these calculations are very technical and beyond the scope of this paper. Instead here, we proceed with numerical simulations.
Numerical simulations for $J}_{EE}>{J}_{EE}^{\text{*}$
Figure 5 depicts the results of our numerical simulations of Model 1 for the same parameters as in Figure 4B (see Materials and methods, Tables 3–4). The response of PV neurons is nonparadoxical: the activity of the PV population increases monotonically with ${\Gamma}_{opto}$ in the whole range (Figure 5A). Concurrently, the population activities of PC, SOM and VIP neurons monotonically decrease with ${\Gamma}_{opto}$ (Figure 5AB). For sufficiently large ${\Gamma}_{opto}$ , PCs become very weakly active and the SOM and VIP populations dramatically reduce their firing rates. The variations with ${\Gamma}_{opto}$ of r_{E}, r_{I}, r_{S} and r_{V} and are robust to changes in the average connectivity, K (Figure 5—figure supplement 1) and in qualitative agreement with the predictions of the large N, K limit (Figure 4B Appendix 1C, Figure 4—figure supplement 2).
To test the robustness of our results with respect to changes in the interaction strengths, we generated 100 networks with J _{αβ} chosen at random within a range of ±10% of those of Figure 4B. All the networks exhibited a balanced state which was stable with respect to slow rates fluctuations in the large N, K limit. We simulated those networks with K = 500 and computed the population activity at baseline and for ${\Gamma}_{opto}=0.07mW.m{m}^{2}$ . For all these networks, the results were consistent with the one of the control set: for ${\Gamma}_{opto}=0.07mW.m{m}^{2}$ , ${r}_{I}$ was larger and r_{E}, r_{S}, r_{V} were smaller than baseline (Figure 5—figure supplement 2). However, a small percentage of these networks (10%) exhibited oscillations with at most an amplitude 20% of their mean in the firing rates. Apart from that, the results were robust to changes in J _{αβ}.
In contrast to what happens in the large N, K limit (Figure 4B, right panel), in the results depicted in Figure 5 the activity of the PC and VIP populations are not proportional. Moreover, in the large K limit, PC and VIP neurons are inactivated before the SOM population is. For K = 500, VIP is the first population to be silenced followed by the SOM and finally the PC population. Simulations with increasing values of K show that these differences are due to substantial finite K effects (Figure 5—figure supplement 1).
Figure 5 also depicts the changes in the firing rates (normalized to baseline) with ${\Gamma}_{opto}$ for several example neurons. These changes are highly heterogeneous across neurons within each population. Whereas the population average varies monotonically, individual cells activity can either increase or decrease and the response can even be nonmonotonic with ${\Gamma}_{opto}$ .
The heterogeneity in the single neuronal responses are also clear in Figure 6A–B that plots, for two different light intensities, the perturbed firing rate vs. baseline for PCs and PV neurons. Remarkably, in both populations a significant fraction of neuron exhibits a response which is incongruous with the population average. The pie charts in Figure 6 depict the fraction of PCs and PV neurons which increased, decreased, or did not change their firing rates. The fraction of neurons whose activity is almost completely suppressed, is also shown. Remarkably, even for ${\Gamma}_{opto}=1.0mW.m{m}^{2}$ , some of the PCs show an activity increase. Moreover, the fraction of PV neurons whose firing rate increases is less for ${\Gamma}_{opto}=1.0mW.m{m}^{2}$ than ${\Gamma}_{opto}=0.5mW.m{m}^{2}$ . It should be noted that in the model all PV neurons receive the same optogenetic input, therefore, the heterogeneity in the response is not due to whether or not the PV neurons were “infected”. This heterogeneity is solely due to the randomness in the connectivity.
Numerical simulations for $J}_{EE}<{J}_{EE}^{\ast$
Figure 7 depicts the results of our numerical simulations of Model 1 when $J}_{EE}<{J}_{EE}^{\ast$ . Parameters are the same as in Figure 4C (see Materials and methods, Tables 3–5). The population activities of PCs and VIP neurons, r_{E} and r_{V} , decrease monotonically with the laser intensity, ${\Gamma}_{opto}$ . Conversely, the variations of the activities of the PV and SOM populations, r_{I} and r_{S} , are nonmonotonic with ${\Gamma}_{opto}$ . For small light intensities, r_{I} decreases and then abruptly increases with larger ${\Gamma}_{opto}$ ; r_{S} exhibits the opposite behavior. Remarkably, when r_{I} is minimum, r_{S} is maximum for nearly the same value of ${\Gamma}_{opto}$ . We show in Figure 7—figure supplement 1 that this proportional decrease only happens in a small region of parameter space when the determinant of the interaction matrix, ${J}_{\alpha \beta}{\u03f5}_{\beta}$ , is close to zero.
This behavior is qualitatively similar to the one derived in the large N, K limit (Figure 4—figure supplement 2). As suggested by the large N, K analysis, the paradoxical response of the PV neurons in the simulations, is driven by the positive feedback loop PCVIPSOMPC (Figure 4—figure supplement 1). Remarkably, when the activity of the PV neurons is minimum, the PCs are still substantially active (40% of baseline level). This is due to finite K corrections to the large N, K predictions (Figure 7—figure supplement 2). These corrections are strong and scale as $\frac{1}{\sqrt{K}}$ (Appendix 1C). Indeed, even for K as large as 2000, r_{E} is still 25% of the baseline when ${r}_{I}$ is minimum.
We checked the robustness of these results with respect to changes in the interaction parameters as we did for $J}_{EE}>{J}_{EE}^{\text{*}$ . We found that for small light intensity all the 100 simulated networks were operating in the balanced state and exhibited the paradoxical effect (Figure 7—figure supplement 3).
Finally, the single neuron responses are highly heterogeneous. Figure 8 plots the perturbed activities of PCs and PV neurons vs. their baseline firing rates for two light intensities. In Figure 8A, the PV response is paradoxical. This is not the case in Figure 8B. Interestingly, the fraction of PV neurons incongruous with the population activity is larger for ${\Gamma}_{opto}=0.5mW.m{m}^{2}$ than for ${\Gamma}_{opto}=1.0mW.m{m}^{2}$ . For both light intensities the activity of almost all the PCs is decreased.
Fourpopulation network: Model 2
In S1, in the range of laser intensities in which the PV response is paradoxical, the decrease of the PC and PV activity is proportional. This feature of the data can be accounted for in Model 1 but only with a finetuning of the interaction parameters (Figure 7—figure supplement 1 and Figure 7—figure supplement 4). This prompted us to investigate whether a different architecture could account robustly for this remarkable property. Our hypothesis is that this property is a direct consequence of the balance of excitation and inhibition.
Theory in the large $N,K$ limit
We first considered the threepopulation model depicted in Figure 9A. It consists of the PC, PV and SOM populations. SOM neurons receive strong inputs from PCs and PV neurons, but do not interact with each other and do not receive feedforward external inputs. In the large N, K limit, the balance of excitation and inhibition of the SOM population reads (see Materials and methods, Equation 20.2).
Therefore, the activities of the PC and PV populations are always proportional. However, as we show in (Appendix 1D) a threepopulation network with such an architecture cannot exhibit the paradoxical effect.
We therefore considered a network model in which a third inhibitory population, referred to as ‘X’, is added without violating Equation (3) (Figure 9B). This requires that SOM neurons do not receive inputs from X neurons (Appendix 1D). This network exhibits the paradoxical effect if and only if $J}_{SE}{J}_{EX}{J}_{XS}>{J}_{XX}{J}_{ES}{J}_{SE$ , that is if the gain of the positive feedback loop, SOMXPCSOM, is sufficiently strong (Appendix 1D). Obviously, this condition simplifies and reads
Remarkably, this inequality does not depend on J_{EE} . This is in contrast to what happens in Model 1 where the paradoxical effect occurs only if J_{EE} is small enough (see Equation (2)).
As in Model 1, we further required that the activity of the PC population increases with its feedforward external input. This adds the constraint (Appendix 1D):
Equations (35) do not depend on J_{XI} . For simplicity, we take J_{XI } =0 and refer to the resulting architecture as Model 2.
In Figure 9C, the slope of the PV population activity changes from negative to positive while PCs are still active. This is because if SOM neurons are completely suppressed, the loop SOMXPCSOM which is responsible for the paradoxical effect, is not effective anymore. Interestingly, the analytical calculations also show that, when the SOM population activity vanishes, the activity of the X population is maximum. Since the SOM population is inactive before PCs, there is a range of laser intensities where the activity of the latter keeps decreasing while the activity of the PV population increases. Once PCs are inactive, the activity of the X population do not vary with I_{opto} . This is because then they only receive a constant feedforward excitation from outside the network which is balanced by their strong recurrent mutual coupling, J_{XX} .
Simulations for finite K
These features are also observed in our simulations depicted in Figure 10. For small laser intensities, the network exhibits a paradoxical effect where the activities of the PC and PV populations decrease with ${\Gamma}_{opto}$ and in a proportional manner (Figure 10A), until the SOM neurons become virtually inactive (Figure 10B). At that value, r_{I} is minimum and r_{X} is maximum. For larger ${\Gamma}_{opto}$ , r_{I} increases while r_{E} keeps decreasing and is still substantial. After r_{E} has vanished, r_{X} saturates but r_{I} continues to increase. All these results are robust to changes in the connectivity, K (Figure 10—figure supplement 1) as well as to changes in the interaction parameters (Figure 10—figure supplement 2). Single neuron responses are more heterogeneous than in the experimental data (Figure 11). It should be noted however that we did not tune parameters to match the experimental heterogeneity.
Discussion
We studied the response of cortex to optogenetic stimulation of parvalbumin positive (PV) neurons and provided a mechanistic account for it. We photostimulated the PV interneurons in layer 2/3 and layer 5 of the mouse anterior motor cortex (ALM). In layer 2/3 photostimulation increased PV activity and decreased the response of the PCs on average. In contrast, in layer five the response of the PV population was paradoxical: both PC and PV activity decreased on average. This is similar to what we found in the mouse somatosensory cortex (S1) (Li et al., 2019). To account for these results, we first investigated the dynamics of networks of one excitatory and one inhibitory population of spiking neurons. We showed that twopopulation network models of strongly interacting neurons do not fully account for the experimental data. This prompted us to investigate the dynamics of networks consisting of more than one inhibitory population.
We considered two network models both consisting of one excitatory and three inhibitory populations. Interneurons are known to be unevenly distributed throughout the cortex. For instance, SOM neurons have been reported to be most prominent in layer five whereas VIP neurons are mostly found in layer 2/3 (Tremblay et al., 2016). Instead of giving a complete description of these layers and all neuronal populations they include, we propose here models with the minimal number of inhibitory populations that can account for the data.
The three inhibitory populations in Model 1 represent PV, somatostatin positive (SOM) and vasoactive intestinal peptide (VIP) interneurons with a connectivity similar to the one reported in primary visual cortex (Pfeffer et al., 2013) and S1 layer 2/3 (Lee et al., 2013). In Model 2, the first two inhibitory populations likewise represent PV and SOM neurons and the third population, denoted as X, represents an unidentified inhibitory subtype. The main difference with Model one is that here, the third population does not project to SOM neurons.
Depending on network parameters, the response of PV neurons in Model one can be paradoxical or not. To have equal relative suppression of the PCs and PV activities, however, interaction parameters have to be finetuned. In Model 2, the relative changes in the PC and PV activity are the same independent of interaction parameters.
For a twopopulation network, the paradoxical effect only occurs when it is inhibition stabilized (Pehlevan and Sompolinsky, 2014; Wolf et al., 2014). This is because the mechanism requires strong recurrent excitation. In the fourpopulation networks we studied, however, the mechanism responsible for paradoxical effect is different. It involves a disinhibitory loop. In fact, strong recurrent excitation prevents the paradoxical effect in these networks. Therefore, the observation of the paradoxical effect upon PV photoexcitation is not a proof that the network operates in the ISN regime.
Strong vs. weak interactions
Cortical networks consist of a large number (N) of neurons each receiving a large number of inputs (K). Because N and K are large, one expects that a network behaves similar to a network where N and K are infinite. In this limit the analysis is simplified and the mechanisms underlying the dynamics are highlighted. When taking the large K limit one needs to decide how the interaction strengths scale with K. Two canonical scalings can be used: in one the interactions scale as 1/K (Hansel and Sompolinsky, 1992; Hennequin et al., 2018; Knight, 1972; Rubin et al., 2015), in the other as $1/\sqrt{K}$ (Darshan et al., 2017; Renart et al., 2010; Rosenbaum et al., 2017; van Vreeswijk and Sompolinsky, 1996). These differ in the strength of the interactions. For instance, for K = 900 interactions are weaker by a factor 30 in the first scaling than in the second. Importantly, these two scalings give rise to qualitatively different dynamical regimes.
When interactions are strong, the excitatory and inhibitory inputs are both very large (of the order of $K.\frac{1}{\sqrt{K}}=1$ ). They, however, dynamically balance so that the temporal average of the net input and its spatial and temporal fluctuations are comparable to the rheobase (Van Vreeswijk and Sompolinsky, 2005; van Vreeswijk and Sompolinsky, 1998), Appendix 1A). In this balanced regime, the average firing rates of the populations are determined by a set of linear equations: the “balance equations”. These do not depend on the neuronal transfer function. For large but finite $K$ , the network operates in an approximately balanced regime. In this regime, the population activities are well approximated by the balance equations, interspike intervals are highly irregular and firing rates are heterogeneous across neurons.
When the interactions are weak, excitatory and inhibitory inputs are both comparable to the rheobase even when K is large, but their spatial and temporal fluctuations vanish as K increases. The activity of the network is determined by a set of coupled nonlinear equations which depends on the neuronal transfer function. For large but finite K, the firing of the neurons is weakly irregular and heterogeneities mostly arise from differences in the intrinsic properties of the neurons.
In which of these regimes does cortex operate invivo? This may depend on the cortical area and on whether the neuronal activity is spontaneous or driven (e.g. sensory, associative, or motor related). There are, however, several facts indicating that the approximate balanced regime may be ubiquitous. Many cortical areas exhibit highly irregular spiking (Shinomoto et al., 2009) and heterogeneous firing rates (Hromádka et al., 2008; Roxin et al., 2011). Excitatory and inhibitory postsynaptic potentials (PSPs) are typically of the order of 0.2 to 2mV or larger (Levy and Reyes, 2012; Ma et al., 2012; Pala and Petersen, 2015; Seeman et al., 2018). Model networks with PSPs of these sizes and reasonable number of neurons and connections exhibit all the hallmarks of the balanced regime (Amit and Brunel, 1997; Hansel and Mato, 2013; Hansel and van Vreeswijk, 2012; Lerchner et al., 2006; Pehlevan and Sompolinsky, 2014; Argaman and Golomb, 2018; Rao et al., 2019; Roudi and Latham, 2007; Roxin et al., 2011 Van Vreeswijk and Sompolinsky, 2005). Moreover, there is experimental evidence of covariation of excitatory and inhibitory inputs into cortical neurons (Haider et al., 2006; Shu et al., 2003). Finally, in cortical cultures synaptic strengths have been shown to approximately scale as $1/\sqrt{K}$ (Barral and D Reyes, 2016). Therefore in this paper we focused on cortical network models in which interactions are strong, that is of the order of $1/\sqrt{K}$ .
Model 1 accounts for the responses in ALM layer 2/3 and layer 5
In Model 1, whether the network exhibits a paradoxical effect depends on the value of the ratio $\rho ={J}_{EE}/{J}_{EE}^{*}$ where ${J}_{EE}^{*}\equiv {J}_{VE}{J}_{ES}/{J}_{VS}$ . Here, ${J}_{\alpha \beta},\alpha ,\beta \in \{E,S,V\}$ , is the strength of the connection from population β to population $\alpha $ . When ρ > 1, the PV response is nonparadoxical and its activity increase can be substantial well before suppression of the PC activity. On the other hand when ρ > 1, the PV response is paradoxical and the PV activity reaches its minimum for light intensities at which the PCs are still substantially active.
In ALM layer 2/3, the activity of the PV population increases with the light intensity while the activity of the PC decreases on average. Remarkably, our experiments showed that the increase in the PV activity was already substantial for small light intensities, where the PCs were still significantly active. In ALM layer 5 the activity of the PV population initially decreased with the light intensity together with the activity of the PC population. As the light intensity is further increased, the PV activity reaches a minimum after which it increases. At this minimum, the PC activity is still substantial.
Thus, Model 1 accounts for our experimental findings in ALM layer 2/3 provided that J_{EE} is sufficiently large. It accounts for the paradoxical effect in layer 5 provided that J_{EE} is sufficiently small. Note that this does not mean that J_{EE} , is larger in the former layer as compared to the latter. The interactions J_{VE}, J_{ES} and J_{VS} are likely to be layer dependent (Jiang et al., 2015) and therefore so is the value of ${J}_{EE}^{*}$ .
Model 2 accounts for the paradoxical effect in S1 while model 1 would require finetuning
Similar to ALM layer 5, the PV response in S1 is paradoxical. Remarkably however, in S1 the relative suppression of the PC and PV activities is the same for low light intensity. Model 1 can account for this feature only when the interaction parameters are finetuned. In contrast, in Model 2 the comodulation of the PC and PV activities stems from the architecture and therefore occurs in a robust manner. Furthermore, it can equally well account for the fact that in S1 the PV activity reaches its minimum when the PC population is active.
Note that in ALM layer 5 the difference between the slopes of the PC and PV population activities is not significantly different (p>0.05). Therefore, we cannot exclude that Model 2 describes ALM layer 5.
The main difference between Models 1 and 2 is that in Model 1, the third inhibitory population (VIP) projects to SOM neurons while in Model 2, the third population (X) does not. This suggests that population X is not the VIP population. For example, X could be chandelier cells that do not express the PV marker (Jiang et al., 2015) Alternatively, population X could describe the effective interaction of several inhibitory populations with PC and PV neurons.
Models 1 and 2 account for the heterogeneity of single neuron responses
The responses of PCs and PV neurons in the experimental data are highly heterogeneous across cells. Indeed in ALM layer 5 and S1, PV neurons on average show a paradoxical response but at the single neuron level the effect of the laser stimulation is very diverse. Moreover, the firing rate of a neuron can vary monotonically or nonmonotonically with the laser intensity. For instance, when stimulated, the firing rates of many PV neurons increase, although, on average the activity is substantially smaller than baseline. Conversely, for some PV neurons the paradoxical effect is so strong that the laser completely suppresses their activity.
We observed an even larger diversity in single neuron responses in our simulations of Model 1 and 2. We should emphasize that in the simulated networks all the neurons were identical and the cells in the same population received the same feedforward constant external input. The only possible source of heterogeneity therefore comes from the randomness in the network connectivity. The effect of this randomness on the network recurrent dynamics is however nontrivial: one may think that the effect of the fluctuations in the number of connections from neuron to neuron should average out since in the models the number of recurrent inputs per neuron is large (K = 500 or more). This is not what happens because in our simulations populations which are active operate in the balanced excitation/inhibition regime (Roxin et al., 2011; van Vreeswijk and Sompolinsky, 1998; van Vreeswijk and Sompolinsky, 1996). In this state, relatively small homogeneity in the number of connections per neuron is amplified to a substantial inhomogeneity in the response. Thus, strong heterogeneity in the response of neurons is not a prima facie evidence for the heterogeneity of the level of Channelrhodopsin expression in the cells nor is it for the diversity of the single neuron intrinsic properties.
Limitations
We give here a qualitative account for the mechanisms underlying the responses of different cortical areas to optical stimulation. A quantitative analysis of the data, in particular of the heterogeneity is beyond our scope. Such an analysis would require a much larger number of PV neurons. Moreover, it would necessitate the use of more complicated neuronal models making the mathematical analysis intractable, limiting the investigation to simulations only and thus obscuring the mechanisms.
In our experiments, we expressed ReaChR in all PV neurons and in all layers in ALM. In particular, all PV neurons in layer 2/3 and layer five were simultaneously affected by the photostimulus. PCs in layer 2/3 project to layer 5 and receive feedback from the latter (Hooks et al., 2013; Naka and Adesnik, 2016). Interlaminar interactions are likely to also contribute to the effect of the photostimulation.
In our models, we did not take into account such interactions. Including strong connections from layer 2/3 PCs to neurons in layer 5 and/or feedback connections from layer 5 neurons to layer 2/3, could alter our interpretations. In the absence of data that reveal the nature of interlaminar interactions, extending our model to incorporate these is impractical given the large number of parameters to vary. Experiments in ALM and S1 where the optogenetic marker is expressed in only one layer at a time would constraint models which include interlaminar interactions and facilitate their analysis (Moore et al., 2018).
There is a large amount of experimental evidence indicating that different synapses can exhibit diverse dynamics depending on their pre and postsynaptic populations (Ma et al., 2012). For instance, recent studies have shown that PCs to PV synapses are depressing while the PCs to SOM synapses are highly facilitating (Karnani et al., 2016; Xu et al., 2013). Synaptic facilitation and depression mechanisms could give rise to dynamics which will make the network responses depend on the duration of the photostimulation. Here, we did not take into account short term plasticity. Mice neocortex mostly comprises PV, SOM and 5HT3aR expressing interneurons. There is a growing amount of experimental evidence indicating that these populations include different subtypes which may have distinct connectivity patterns (Naka and Adesnik, 2016; Nigro et al., 2018; Tremblay et al., 2016). In the present work, we only considered three populations of identical interneurons: PV, SOM and VIP or X. As the number of populations increases, the number of interaction parameters increases quadratically, making it a great challenge to uncover even simple mechanisms that could underlie the network responses.
Comparison with previous theoretical work
The paradoxical effect was first described in Tsodyks et al. (1997) and Ozeki et al. (2009) for weak interactions using coarse grained twopopulation rate models (Wilson and Cowan, 1972). These models were extended in Rubin et al. (2015) to a spatially structured network to explain centersurround interactions and other contextual effects in primary visual cortex. They found that these effects can be accounted for if the neuronal transfer function is supralinear and the network is operating in the inhibition stabilized regime (ISN). With supralinear transfer functions, whether or not the network exhibits a paradoxical effect depends on the background rate of the inhibitory neurons. These models were further extended by LitwinKumar et al. (2016) to networks consisting of PC, PV, SOM and VIP neurons with an architecture similar to Pfeffer et al. (2013). They studied the effect of photostimulation of the different inhibitory populations on the responses and orientation tuning properties of the neurons. In a recent study (Sadeh et al., 2017) have investigated the effects of partial activation of PV neurons upon photostimulation in an ISN. They argued that depending on the degree of viral expression, the average response of the infected neurons can decrease or increase with the light intensity: it decreases only if a large proportion of the population is infected. (Garcia Del Molino et al., 2017) showed that due to the nonlinearity in the neuronal transfer function, the response of the network to stimulation can be different for different background rates. In particular, they showed that it can reverse the response of SOM neurons to VIP stimulation.
All these works considered inhibition stabilized networks in which the total recurrent excitation is so strong that the activity would blow up in the absence of inhibitory feedback. With our notations, this means that ${G}_{E}{j}_{EE}>1/K$ , where G_{E} is the gain of the noise average transfer function (fI curve) of the excitatory neurons. In fact, in these models all the interactions j _{αβ} are of order 1/K so they are weak in our sense. Moreover, these studies considered networks that are so small that it is impossible to extrapolate their results to mouse cortex size networks. Here we studied large network models (N = 76800) with strong interactions, that is j _{αβ} are of order $1/\sqrt{K}$ operating in the balanced regime. Note that such networks are ISNs provided that ${j}_{EE}\ne 0$ . We showed that paradoxical effect can be present or not depending on the interaction parameters.
Since we used static synapses, changes in the background rates cannot reverse the paradoxical effect in our models. This is because with static synapses the balance equations are linear. One can recover this reversal if one introduces shortterm plasticity which will make the balance equations nonlinear. We did not consider partial expression of channelrhodopsin in the PV population because our goal was to account for experimental data where virtually all neurons were infected. These effects have been studied in Gutnisky et al. (2017); Sanzeni et al. (2019) in strongly coupled networks of two populations yielding to the same conclusions as (Sadeh et al., 2017).
Predictions
Our theory (Model 1) predicts that in ALM layer 2/3 the activity of the SOM and VIP populations will decrease upon PV photostimulation (Figure 4B). It also predicts that upon PC photoinhibition, the PV activity will increase whereas the activity of the SOM and VIP populations will decrease (Figure 12A). This is because in Model 1 when the PV response is nonparadoxical ( ${\chi}_{II}>0$ ) the product X_{EI }X_{IE} is also positive (see Appendix 1C). Furthermore, in ALM layer 2/3 the population activity of PCs decreases upon PV photostimulation, X_{EI} < 0. Hence, X_{IE} is negative. The balance of the PC and the VIP inputs into SOM neurons implies that VIP and PC activity covary. Finally, in Appendix 1C we show that if X_{EE } > 0 and X_{IE} < 0 then necessarily X_{SE} > 0. Thus, in ALM layer 2/3, the SOM population activity should decrease upon PC photoinhibition (Figure 12A).
In auditory and prefrontal cortex (Pi et al., 2013) as well as in S1 (Lee et al., 2013), photostimulation of VIP neurons, activates them (X_{VV} > 0) and disinhibits the PCs (X_{EV} > 0) through an inhibition of the SOM population (X_{SV} > 0). If this is also true in ALM layer 2/3, our model predicts that photostimulation of VIP neurons should increase the PV activity (X_{IV} > 0) (Appendix 1C, Figure 12B).
In S1 our theory (Model 2) predicts that the PC and PV activities will proportionally decrease upon PC photoinhibition (Equation (3), Appendix 1D, Figure 12C). Photostimulation of the SOM neurons modifies Equation (3) and consequently, the changes in PC and PV activity no longer covary (Figure 12D). Thus, our theory can be tested by photostimulating PV neurons as in our experiment, while also photostimulating SOM neurons with a second laser with constant power. In this case, the model predicts that S1 will still exhibit the paradoxical effect but that the responses of the PC and PV populations will no longer be proportional (Figure 12E).
Perspectives
We only considered response of the neurons for a large radius of the laser beam. In a recent study Li et al. (2019), have investigated the spatial profile of the response and its dependence on the light intensity. Our theory can be extended to incorporate spatial dependencies. Studying the interplay between the connectivity pattern and laser beam width in the response profile of the networks will provide further constraints on cortical architectures.
Due to the strong interactions in our models, the nonlinearity of the single neuron fI curves hardly affects the population average responses. However, it influences the response heterogeneity that naturally arises in our theory (Figures 6–8). An alternative model for the paradoxical effect is the supralinear stabilized network (SSN) (Rubin et al., 2015) which relies on an expansive nonlinearity of the inputoutput transfer function of the inhibitory populations. Whether this mechanism can account for our experimental data is an issue for further study. In particular, it would be interesting to know whether the SSN scenario can account for the strong heterogeneity in the responses and for the proportionality of the PC and PV population activities in S1. Answering these questions may provide a way to discriminate between the balance network and SSN theory.
Materials and methods
Animals and surgery
Request a detailed protocolThe experimental data are from 9 PVIresCre x R26CAGLSLReaChRmCitrine mice (age >P60, both male and female mice) (Hooks et al., 2015). three mice were used for photoinhibition in somatosensory cortex (S1). six mice were used for photoinhibition in anterior lateral motor cortex (ALM). All procedures were in accordance with protocols approved by the Janelia Research Campus and Baylor College of Medicine Institutional Animal Care and Use Committee.
Mice were prepared for photostimulation and electrophysiology with a clearskull cap and a headpost (Guo et al., 2014a; Guo et al., 2014b). The scalp and periosteum over the dorsal surface of the skull were removed. A layer of cyanoacrylate adhesive (Krazy glue, Elmer’s Products Inc) was directly applied to the intact skull. A custom made headbar was placed on the skull (approximately over visual cortex) and cemented in place with clear dental acrylic (Lang Dental Jet Repair Acrylic; Part# 1223clear). A thin layer of clear dental acrylic was applied over the cyanoacrylate adhesive covering the entire exposed skull, followed by a thin layer of clear nail polish (Electron Microscopy Sciences, Part# 72180).
Photostimulation
Request a detailed protocolLight from a 594 nm laser (Cobolt Inc, Colbolt Mambo 100) was controlled by an acoustooptical modulator (AOM; MTS110A3VIS, Quanta Tech; extinction ratio 1:2000; 1µs rise time) and a shutter (Vincent Associates), coupled to a 2D scanning galvo system (GVA002, Thorlabs), then focused onto the brain surface (Guo et al., 2014a). The laser at the brain surface had a diameter of 2 mm. We tested photoinhibition in barrel cortex (bregma posterior 0.5 mm, 3.5 mm lateral) and ALM (bregma anterior 2.5 mm, 1.5 mm lateral).
To prevent the mice from detecting the photostimulus, a ‘masking flash’ pulse train (40 1 ms pulses at 10 Hz) was delivered using a LED driver (Mightex, SLA1200–2) and 590 nm LEDs (Luxeon Star) positioned near the eyes of the mice. The masking flash began before the photostimulus started and continued through the end of the epoch in which photostimulation could occur.
The photostimulus had a near sinusoidal temporal profile (40 Hz) with a linear attenuation in intensity over the last 100–200 ms (duration: 1.3 s including the ramp). The photostimulation was delivered at ~7 s intervals. The power (0.5, 1.2, 2.2, 5, 12 mW for S1 photostimulation; 0.3, 0.5, 1, 1.5, 2, 3.3, 5, 8, 15 mW for ALM photostimulation) were chosen randomly. Because we used a timevarying photostimulus, the power values reported here reflect the timeaverage.
Electrophysiology
Request a detailed protocolAll recordings were carried out while the mice were awake but not engaged in any behavior. Extracellular spiking activity was recorded using silicon probes. We used 32channel NeuroNexus silicon probes (A4 × 8–5 mm100200177) or 64channel Cambridge NeuroTech silicon probes (H2 acute probe, 25 μm spacing, two shanks). The 32channel voltage signals were multiplexed, digitized by a PCI6133 board at 400 kHz (National Instruments) at 14 bit, demultiplexed (sampling at 25,000 Hz) and stored for offline analysis. The 64channel voltage signals were amplified and digitized on an Intan RHD2164 64Channel Amplifier Board (Intan Technology) at 16 bit, recorded on an Intan RHD2000Series Amplifier Evaluation System (sampling at 20,000 Hz) using OpenSource RHD2000 Interface Software from Intan Technology (version 1.5.2), and stored for offline analysis.
A 1 mm diameter craniotomy was made over the recording site. The position of the craniotomy was guided by stereotactic coordinates for recordings in ALM (bregma anterior 2.5 mm, 1.5 mm lateral) or barrel cortex (bregma posterior 0.5 mm, 3.5 mm lateral).
Prior to each recording session, the tips of the silicon probe were brushed with DiI in ethanol solution and allowed to dry. The surface of the craniotomy was kept moist with saline. The silicon probe was positioned on the surface of the cortex and advanced manually into the brain at ~3 µm/s, normal to the pial surface. The electrode depth was inferred from manipulator depth and verified with histology. For ALM recordings, putative layer 2/3 units were above 450 µm and putative layer 5 units were below 450 µm (Hooks et al., 2013). For S1, our recording did not distinguish layers.
Data analysis
Request a detailed protocolThe extracellular recording traces were bandpass filtered (300–6 kHz). Events that exceed an amplitude threshold (four standard deviations of the background) were subjected to manual spike sorting to extract single units (Guo et al., 2014a).
Our final data set comprised of 204 single units (S1, 95; ALM, 109). For each unit, its spike width was computed as the trough to peak interval in the mean spike waveform (Guo et al., 2014a). We defined units with spike width <0.35 ms as FS neurons (31/204) and units with spike width >0.45 ms as putative pyramidal neurons (170/204). Units with intermediate values (0.35–0.45 ms, 3/204) were excluded from our analyses.
To quantify photoinhibition strength, we computed ‘normalized spike rate’ during photostimulation. For each neuron, we computed its spike rate during the photostimulus (1 s time window) and its baseline spike rate (500 ms time window before photostimulus onset). The spike rates under photostimulation were divided by the baseline spike rate. The ‘normalized spike rate’ thus reports the total fraction of spiking output under photostimulation. For normalized spike rate of individual neurons, each neuron’s spike rate with photostimulation was normalized by dividing its baseline spike rate (Figure 1B–D, top). For normalized spike rate of the neuronal population (Figure 1B–D, bottom), the spike rates with photostimulation were first averaged across the population (without normalization) and then normalized by dividing the averaged baseline spike rate.
Bootstrap was performed over neurons to obtain standard errors of the mean. For each round of bootstrapping, repeated 1000–10000 times, we randomly sampled with replacement neurons in the dataset. We computed the means of the resampled datasets. The standard error of the mean was the standard deviation of the mean estimates from bootstrap.
Network models
Request a detailed protocolAll the models we consider consist of strongly interacting leaky integrateandfire neurons. We first study networks of one excitatory (E) and one inhibitory (I) population. We then investigate two models comprising three inhibitory populations, namely parvalbumin positive (PV or I), somatostatin positive (SOM or S) and a third population either corresponding to the vasoactive intestinal peptide positive (VIP or V) neurons (Model 1) or to an unidentified population denoted by X (Model 2).
In all models the total number of neurons is N = 76800. In the two population model, 75% are excitatory and 25% inhibitory. In the fourpopulation networks, 75% are excitatory and the number of cells is the same, N/12, for all GABAergic inhibitory population.
The data we seek to account for, were obtained in optogenetic experiments in which the laser diameter was substantially larger than the spatial range of neuronal interactions and comparable to the size of the cortical area were the recordings were performed. Therefore, in all models we assume for simplicity that the connectivity is unstructured: neuron (i, α), (α = E, I, S, V/X), is postsynaptically connected to neuron (j) (j, β) with probability
For simplicity, we take ${K}_{\alpha \beta}$ the same for all populations, ${K}_{\alpha \beta}=K$ .
Neuron dynamics: The dynamics between spikes of the membrane potential of the neuron (i, α) is given by
Here, ${I}_{rec}^{\alpha i}\left(t\right)$ is the net recurrent input into neuron $\left(i,\alpha \right)$ , ${\Lambda}_{ext}^{\alpha}$ represents inputs from outside the circuit (e.g. thalamic excitation) to population α, and ${\Lambda}_{opto}^{\alpha i}$ is the optogenetic input into neuron (i, α).
We assumed that the capacitance, C_{M} , is identical for all neurons and the leak conductance, ${g}_{leak}^{\alpha}$ , is identical for all the cells in the same population. We take ${C}_{M}=1\mu F.c{m}^{2}$ , ${g}_{leak}^{I}=0.1mS.c{m}^{2}$ and ${g}_{leak}^{E}$ = ${g}_{leak}^{S}$ = ${g}_{leak}^{V/X}=0.05mS.c{m}^{2}$ .
Equation (2) has to be supplemented by a reset condition: if at time $t$ the membrane potential of the neuron (i, α) crosses the threshold ${V}_{i}^{\alpha}({t}^{})\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}{V}_{th}\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}50mV$ , the neuron fires a spike and its voltage is reset to the resting potential ${V}_{i}^{\alpha}({t}^{+})\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}{V}_{R}\phantom{\rule{thinmathspace}{0ex}}=\phantom{\rule{thinmathspace}{0ex}}70mV$ .
Recurrent inputs: The net recurrent input into neuron (i, α) is
where C ^{αβ} is the connectivity matrix between (presynaptic) population β and (postsynaptic) population α, such that ${C}_{ij}^{\alpha \beta}=1$ if neuron (j, β) projects to neuron (i, α) and ${C}_{ij}^{\alpha \beta}=0$ otherwise. The parameter j _{αβ} is the strength of the interaction from neurons in population β to neurons population α. We assumed it to depend on the pre and postsynaptic populations only. The polarity (excitation or inhibition) of the interaction is denoted by ε_{β}. Therefore if β = E, ε_{β} = 1 and ε_{β} = 1 otherwise.
The function ${S}_{j}^{\alpha \beta}\left(t\right)$ is
where ${t}_{\beta j}^{k}$ is the time at which neuron (j, β) has emitted its k ^{th} spike, the sum is over all the spikes emitted by neuron (j, β) prior to time t and
where τ_{αβ} is the synaptic time constant of the interactions between neurons in population β and α.
External and optogenetic inputs: The feedforward input, ${\Lambda}_{ext}^{\alpha}$ , into the neurons in population $\alpha $ is described by inputs from 2K external neurons with constant firing rate r _{0} = 5 Hz and an interaction strength j _{α0}, therefore, ${\Lambda}_{ext}^{\alpha}={2Kj}_{\alpha 0}{r}_{0}$ .
We model the ReachR photostimulation as an additional external constant input to the stimulated population. For simplicity, we assume that this input, ${\Lambda}_{opto}^{\alpha i}={\Lambda}_{opto}^{\alpha}$ , is the same for all stimulated neurons. Unless specified otherwise, we only consider ${\Lambda}_{opto}^{I}={\Lambda}_{opto}$ and ${\Lambda}_{opto}^{\alpha}=0$ for $\alpha \ne I$ .
In qualitative agreement with Figure 3, and Figures 5, 7, Figure 7—figure supplement 1, Figure 10; (Hooks et al., 2015) we take
where ${\Gamma}_{opto}$ is the laser intensity and ${\Lambda}_{0}$ and ${\Gamma}_{0}$ are parameters.
Architectures of the fourpopulation models
Request a detailed protocolThe network of Model one is depicted in Figure 4A. In line with the results of Pfeffer et al. (2013), there are no connections from PV to SOM, VIP to PC and VIP to PV neurons. There is no mutual inhibition between SOM as well as between VIP neurons. All the populations except SOM receive feedforward external input.
The interaction matrix of the network is
The network of Model two is depicted in Figure 9B. SOM only receives projections from PCs and PV neurons. X neurons are recurrently connected and project to PCs and PV neurons. The PC and SOM populations project to the population X. All the populations except SOM receive feedforward external input.
The interaction matrix is
Numerical simulations: The dynamics of the models was integrated numerically using a secondorder RungeKutta scheme (Press et al., 1986) without spike time interpolation. Unless specified otherwise the time step was Δt = 0.01 ms and the temporally averaged firing rates were estimated over 100s.
The balance equations
Request a detailed protocolWe consider recurrent networks of strongly interacting neurons (van Vreeswijk and Sompolinsky, 1996) in which order $\sqrt{K}$ excitatory synaptic inputs are sufficient to bring the voltage above threshold. To understand the behavior of such networks, it is imperative to analyse how it behaves when K goes to infinity. To this end, we scale the interactions as
where J _{αβ} does not depend on K. Since a neuron receives on average K inputs from each of its presynaptic populations, the total interaction from population β to a neuron in population α is ${J}_{\alpha \beta}\sqrt{K}$ . To keep the relative strength of the optogenetic input, ${\Lambda}_{opto}^{\alpha}$ , as $K$ increases we take
where ${I}_{opto}^{\alpha}$ depends on the intensity of the laser:
We take: ${I}_{0}^{\alpha}={I}_{0}=8nA$ and ${\Gamma}_{0}^{\alpha}={\Gamma}_{0}=0.5mW.m{m}^{2}$ .
The net input into the neurons must remain finite in the infinite K limit. This implies that up to corrections which are of the order of $\frac{1}{\sqrt{K}}$ ,
In a npopulation network, these $n$ equations determine the $n$ firing rates, ${r}_{\alpha},\alpha \in \{1,...,n\}$ .
This set of linear equations express the fact that, for the population activities to be finite, excitatory and inhibitory inputs to the neurons must compensate. These 'balance' equations have a unique solution (unless the determinant of the matrix ${J}_{\alpha \beta}{\u03f5}_{\beta}$ is zero). To be meaningful the solution must be such that all population activities are positive. This constrains the feedforward and recurrent interaction parameters.
The stability of this balanced solution further constraints the interaction parameters and synaptic time constants. A necessary condition for the stability is that $det\left[{J}_{\alpha \beta}{\u03f5}_{\beta}\right]>0$ . This condition guarantees that the 'balanced state' is stable with respect to divergence of the firing rates. A complete study of these constraints for our LIF networks is beyond the scope of this paper.
In all the models, we study parameter ranges in which, at baseline ( ${I}_{opto}^{\alpha}=0$ ), the network operates in a stable balanced state where distributions of rates exhibit a quasilognormal shape and spikes are emitted irregularly as in a Poisson process (Figure 5—figure supplement 3; Figure 7—figure supplement 5; Figure 10—figure supplement 3). For ${I}_{opto}^{\alpha}$ sufficiently large, it may happen that one or more population activity reaches zero. In this case, the network evolves to a partially balanced state in which the rates of the populations that remain active satisfy a reduced set of balanced equations. For example, if we consider a solution were the rate of population $\gamma $ , ${r}_{\gamma}$ is zero and all other rates are positive, the reduced balance equations are
Consistency of this solution leads to the requirement that the input into population $\gamma $ is hyperpolarizing.
Note that they may be multiple selfconsistent solutions which are partially balanced.
Upon photostimulation of PV, in Model 1, the balanced equations are
In particular, Equation (19.3) implies that ${r}_{E}$ and ${r}_{V}$ are always proportional ( ${J}_{SE},{J}_{SV}>0$ ).
Similarly, in Model 2, the balanced equations are
Equation (20.3) implies that in this network ${r}_{E}$ and ${r}_{I}$ are always proportional $\left({J}_{SE},\text{}{J}_{SI}0\right)$ .
Appendix 1
Mean field theory
Let us consider a network consisting of n populations (e.g. n = 4) receiving feedforward input, ${\Lambda}_{ext}^{\alpha}$ , from an external population with constant firing rate, r _{0} ${r}_{0}$ , and an optogenetic input, ${\Lambda}_{opto}^{\alpha}$ (Materials and Methods). The total input into neuron (i, α) is
If the size of the network, N, and mean connectivity, K are large and the synaptic time constants are sufficiently small compared to the membrane time constants, one can take the diffusion approximation and neglect the temporal correlations and write
where ${\zeta}_{i}^{\alpha}$ is an i.i.d. Gaussian with zero mean and unit variance, and ${\eta}_{i}^{\alpha}\left(t\right)$ is a Gaussian white noise with zero mean and unit variance. The mean input, ${u}_{\alpha}$ , is
where the population average firing rate of population β is ${r}_{\beta}=\text{[}{r}_{j}^{\beta}\text{]}$ and ${r}_{j}^{\beta}$ is the firing rate of the neuron (j, β). Here <.> denotes temporal average (i.e. over ${\eta}_{i}^{\alpha}\left(t\right)$ ) and $\left[.\right]$ is the average over the quenched disorder ( ${\zeta}_{i}^{\alpha}$ ). The latter stems from heterogeneities in the indegree of the inputs into the neurons.
In Equation (A2), ${A}_{\alpha}$ is the variance of the quenched disorder which is given by
while ${B}_{\alpha}$ is the variance of the temporal fluctuations (Van Vreeswijk and Sompolinsky, 2005; Roxin et al., 2011)
In Equation (A4), ${q}_{\beta}\text{=[(}{r}_{j}^{\beta}{\text{)}}^{2}\text{]}$ .
Equations (A4A5) have to be supplemented with the expression of the inputoutput transfer function which relates the average firing rate, ${r}_{i}^{\alpha}$ , to the statistics of ${I}_{tot}^{\alpha i}\left(t\right)$ ,
where $D\zeta$ = $\frac{1}{\sqrt{2\pi}}{e}^{{\zeta}^{2}/2}$ , and $\mathrm{\Phi}}_{\alpha$ is given by Capocelli and Ricciardi (1971)
where $X}_{\alpha}^{}=\frac{x{g}_{leak}^{\alpha}{V}_{R}}{\sqrt{y}$ , $X}_{\alpha}^{+}=\frac{x{g}_{leak}^{\alpha}{V}_{Th}}{\sqrt{y}$ and $\tau}_{\alpha}=\frac{{C}_{M}}{{g}_{leak}^{\alpha}$ is the membrane time constant of the neurons in population $\alpha$ .
With $j}_{\alpha \beta}=\frac{{J}_{\alpha \beta}}{\sqrt{K}$ , $\mathrm{\Lambda}}_{ext}^{\alpha}=2\sqrt{K$ and $\mathrm{\Lambda}}_{opto}^{\alpha}={I}_{opto}^{\alpha}\sqrt{K$ (see Materials and methods), we obtain
For finite, but large K, the average activity of population $\alpha $ is
where ${\Psi}_{\alpha}$ is the right handside of Equation (A7).
In the limit where ${u}_{\alpha}\to \infty $ , it can be shown that
In the large K limit, the activities, r_{a} , have to satisfy a set of n linear balance equations (Equation (12), Materials and methods) and are given by
We define the susceptibility matrix, X _{αβ} , as the derivative of the activity, r_{a} , with respect to ${I}_{opto}^{\beta}$ ,
At baseline $\left({I}_{opto}^{\beta}=0\right)$ , the positivity of ${r}_{\alpha},\forall \alpha $ imposes conditions on the recurrent and feedforward interaction strengths, ${J}_{\alpha \beta}$ and ${J}_{\alpha 0}$ . The requirement that there are no 'partially' balanced solutions for which one or more of the $n$ populations is inactive or saturates and the stability of the balanced solution imposes further constraints.
Twopopulation model
Large $K$ limit
For a twopopulation (one excitatory E and one inhibitory I) network, solving Equation (13) gives for a perturbation, ${I}_{opto}$ , upon I,
where ${\Delta ={J}_{EI}{J}_{IE}J}_{EE}{J}_{II}$ .
The requirement that at baseline the network state is fully balanced and stable implies that
Therefore, $\mathrm{\Delta}>0$ .
The susceptibilities with respect to a perturbation of I are
which both are negative. Therefore, ${r}_{E}$ and ${r}_{I}$ decrease linearly with ${I}_{opto}$ , that is the response of the I population is paradoxical.
It is useful to consider the susceptibilities normalized to baseline rate
Equation (A19) implies that, $\left{\overline{\chi}}_{EI}\right$ is larger than $\left{\overline{\chi}}_{II}\right$ .
Moreover, whereas ${\overline{\chi}}_{EI}$ is independent of ${J}_{EE}$ , ${\overline{\chi}}_{II}$ depends on ${J}_{EE}$ . When ${J}_{EE}=0$ , ${\overline{\chi}}_{II}$ is zero: the PV activity is insensitive to ${I}_{opto}$ .
The identity of the two normalized susceptibilities can only be achieved with a finetuning of the interaction parameters such that $\Delta \simeq 0$ for
Concurrently, as ${J}_{EE}\to {J}_{EI}{J}_{IE}/{J}_{II}$ , the activity of the two populations diverge as $\frac{1}{\Delta}$ with a constant ratio equal to $\frac{{J}_{IE}}{{J}_{II}}$ . Thus, to keep the activities finite, $2\left({J}_{II}{J}_{E0}{J}_{EI}{J}_{I0}\right){r}_{0}$ and $2\left({J}_{IE}{J}_{E0}{J}_{EE}{J}_{I0}\right){r}_{0}$ must also tend to zero.
Finally, if ${I}_{opto}\text{=}{I}_{opto}^{\text{*}}\equiv 2\left({J}_{E0}{J}_{II}/{J}_{EI}{J}_{I0}\right){r}_{0}$ , ${r}_{E}$ vanishes (Figure 3—figure supplement 1). When $I}_{opto}\text{>}{I}_{opto}^{\text{*}$ , the balance between the total external excitatory (optogenetic+feedforward) and recurrent inhibitory inputs into I implies that ${r}_{I}$ linearly increases with ${I}_{opto}$ and the slope is $1/{J}_{II}$ .
Finite K corrections to ${r}_{E}$ and ${r}_{I}$ near ${I}_{opto}^{\text{*}}$
When $K$ is finite, ${r}_{I}$ starts to increase with ${I}_{opto}$ when ${r}_{E}$ is exponentially small in $K$ . To show that, we have to derive the leading order correction to the activities near ${I}_{opto}^{\text{*}}$ .
We make the ansatz that when ${I}_{opto}={I}_{opto}^{\text{*}}+\delta I\sqrt{\frac{log\left(K\right)}{K}}$ , ${r}_{E}={\nu}_{E}\frac{\sqrt{log\left(K\right)}}{K}$ and ${r}_{I}={r}_{I}^{\infty}+{\nu}_{I}\sqrt{\frac{log\left(K\right)}{K}}$ , where ${\nu}_{E}$ and ${\nu}_{I}$ are $O\left(1\right)$ and ${r}_{I}^{\infty}=2{J}_{E0}{r}_{0}/{J}_{EI}$ is the inhibitory firing rate at ${I}_{opto}\text{=}{I}_{opto}^{\text{*}}$ in the large $K$ limit.
To leading order:
where ${A}_{\alpha}^{\infty}$ and ${B}_{\alpha}^{\infty}$ , $\alpha \in \{E,I\}$ , are the variance of the temporal and quenched noise in the large $K$ limit (Equations (A11A12)).
Equation (A25.1) implies that
Together with Equation (A25.2) one obtains
where $\Delta ={J}_{EI}{J}_{IE}{J}_{EE}{J}_{II}$ .
For large $K$ ,
where $Q=\frac{1}{{\tau}_{m}^{E}\sqrt{\pi}}\frac{{B}_{E}^{\infty}}{{\left(2{A}_{E}^{\infty}+{B}_{E}^{\infty}\right)}^{3/2}}$ .
Since ${\nu}_{E}$ must be positive, $\left({J}_{EI}\delta I+{\nu}_{E}\Delta \right)$ must also be positive, Equation (A28) then implies that to leading order
Hence, ${\nu}_{I}$ is
Therefore, both ${\nu}_{E}$ and ${\nu}_{I}$ decrease with $\delta I$ . This holds for $\delta I\le \frac{{J}_{II}}{{J}_{EI}}\sqrt{{A}_{E}^{\mathrm{\infty}}+\frac{{B}_{E}^{\mathrm{\infty}}}{2}}$ . Beyond this range ${r}_{E}$ is exponentially small, ${\nu}_{I}=\frac{\delta I}{{J}_{II}}$ and ${r}_{I}$ increases with ${I}_{opto}$ .
In conclusion, when the response of the I population is minimum the firing rate of the excitatory population is exponentially small in $K$ .
Fourpopulation model: Model 1
Large $K$ limit
In Model 1, the population susceptibilities in response to a perturbation of the PV population are given by Equation (A16)
where $\Delta =det\left(\left[{J}_{AB}{\u03f5}_{B}\right]\right)$ .
Note, in this model we do not take into account any PV to SOM connections. Nevertheless even If one includes these, the expressions of the PC and PV susceptibility will only differ by a scaling factor from the ones in A31 and A32 (because of $\Delta $ ) and therefore their sign will depends on the same conditions than A31 and A32.
Interestingly, for stable solutions $\left(\mathrm{\Delta}>0\right)$ , then ${\chi}_{II}>0$ implies that $J}_{EE}\text{}{J}_{VS}{J}_{ES}\text{}{J}_{VE$ . while ${\chi}_{EI}<0$ implies that $J}_{ES}\text{}{J}_{VI}{J}_{EI}\text{}{J}_{VS$ . Therefore, $J}_{EE}\text{}{J}_{VS}\text{}{J}_{VI}{J}_{VE}\text{}{J}_{ES}\text{}{J}_{VI$ . and $J}_{ES\text{}}{J}_{VI}\text{}{J}_{VE}{J}_{EI}\text{}{J}_{VS}\text{}{J}_{VE$ . Combining the latter one has $J}_{EE}\text{}{J}_{VS}\text{}{J}_{VI}{J}_{EI}\text{}{J}_{VS}\text{}{J}_{VE$ . Therefore, $J}_{EE}\text{}{J}_{VI}{J}_{EI}\text{}{J}_{VE$ which is equivalent to ${\chi}_{SI}<0$ .
Similarly one can show that if ${\chi}_{EE}>0$ and ${\chi}_{IE}<0$ necessarily ${\chi}_{SE}>0$ .Let us consider a particular set of parameters for which a stable balanced solution exists when ${J}_{EE}=0$ $\left(\mathrm{\Delta}\left(0\right)>0\right)$ .
The susceptibility ${\chi}_{II}$ as a function of ${J}_{EE}$ is
where ${\hat{\chi}}_{EE}\equiv {\chi}_{EE}.\mathrm{\Delta}\left({J}_{EE}\right)={J}_{SV}\left({J}_{VI}{J}_{IS}{J}_{II}{J}_{VS}\right)$ , is the numerator in the susceptibility ${\chi}_{EE}$ .
In our models, we assumed ${\chi}_{EE}>0$ . When ${J}_{EE}=0$ , $\Delta \left(0\right)$ is positive thus, ${\chi}_{II}\left(0\right)<0$ . As ${J}_{EE}$ increases, the sign of ${\chi}_{II}\left({J}_{EE}\right)$ depends on the order relationship between two quantities. The first one, ${J}_{EE}^{\text{*}}$ , is the value of ${J}_{EE}$ for which the numerator in Equation (A35) changes sign
The second one, ${J}_{EE}^{c}$ , is defined by $\Delta \left({J}_{EE}^{c}\right)=0$
Therefore, for $J}_{EE}>{J}_{EE}^{c$ , the dynamics is unstable. Two cases can be distinguished:
If $J}_{EE}^{\text{*}}<{J}_{EE}^{c$ , then ${\chi}_{II}$ is an increasing function of ${J}_{EE}$ . It is negative if $J}_{EE}<{J}_{EE}^{\text{*}$ and becomes positive for $J}_{EE}>{J}_{EE}^{\text{*}$ .
If $J}_{EE}^{\text{*}}<{J}_{EE}^{c$ , ${\chi}_{II}$ is a decreasing function of ${J}_{EE}$ and is negative in all the region where the dynamics is stable.
The derivative of ${\chi}_{II}$ , (Equation (A35)), with respect to ${J}_{EE}$ , has the same sign as ${\chi}_{EI}{\chi}_{IE}$ . Therefore, ${\chi}_{EI}{\chi}_{IE}$ is positive in the first case and negative in the second.
Experimental data shows that the activity of the PC population decreases upon PV photostimulation, i.e., ${\chi}_{EI}<0$ . Therefore, if ${\chi}_{II}>0$ as in ALM layer 2/3, ${\chi}_{IE}$ must be negative, i.e., the activity of the PV population decreases upon PC photostimulation.
Finite $K$
When ${I}_{opto}$ is sufficiently strong, a fully balanced solution $\left({r}_{\alpha}>0,\mathrm{\forall}\alpha \right)$ no longer exists in our case ${r}_{E}={r}_{V}=0$ for $I}_{opto}>{I}_{opto}^{\text{*}$ where
To understand the network behavior after this point we need to consider finite $K$ corrections.
Since the PC and VIP population activities decrease with ${I}_{opto}$ , when ${I}_{opto}$ is sufficiently large and due to the balance of the SOM input, ${r}_{E}$ and ${r}_{V}$ will both be at most $O\left(\frac{1}{\sqrt{K}}\right)$ . Let us write: ${r}_{E}\equiv \frac{{\nu}_{E}}{\sqrt{K}}$ and ${r}_{V}\equiv \frac{{\nu}_{V}}{\sqrt{K}}$ where ${\nu}_{E}$ and ${\nu}_{V}$ are at most $O\left(1\right).$
One should consider four cases:
1) ${\nu}_{E}$ and ${\nu}_{V}$ are $O\left(1\right)$
In this case, the average net input into the SOM population, ${u}_{S}={J}_{SE}{\nu}_{E}{J}_{SV}{\nu}_{V}$ , is $O\left(1\right)$ and the temporal fluctuations, ${B}_{S}$ , and heterogeneities, ${A}_{S}$ , are negligible. If ${u}_{S}$ is larger than the rheobase, $\left({V}_{th}{V}_{R}\right)/{g}_{leak}^{S}$ , ${r}_{S}$ is also $O\left(1\right)$ . Otherwise, ${r}_{S}=0$ .
Because ${\nu}_{E}$ and ${\nu}_{V}$ are $O\left(1\right)$ , ${u}_{E}$ and ${u}_{V}$ are $o\left(\frac{1}{\sqrt{K}}\right)$ . Thus, to leading order,
Moreover, the balance of the PV population implies that
Thus, there are three linear equations (Equations (A39A41)) for two unknowns $\left({r}_{I}\text{and}{r}_{s}\right)$ . These cannot be satisfied and hence, in this case, there is no consistent solution.
2) ${\nu}_{E}=o\left(1\right)$ and ${\nu}_{V}=O\left(1\right)$
Here, to leading order, ${u}_{S}={J}_{SV}{\nu}_{V}<0$ , while ${A}_{S}={B}_{S}=0$ . As a result, to leading order, ${r}_{S}=0$ . The activity of the PV population is then
Because ${\nu}_{V}$ is $O\left(1\right)$ ,
Equations (A42, A43) cannot both be satisfied. This solution is also inconsistent.
3) ${\nu}_{E}=O\left(1\right)$ and ${\nu}_{V}=o\left(1\right)$
In this case ${u}_{S}={J}_{SE}{\nu}_{E}>0$ and therefore ${r}_{S}$ can be $O\left(1\right)$ . Equations (A39) and (A41) imply
which determine ${r}_{I}$ and ${r}_{S}$ as ${r}_{I}=\frac{\left({J}_{ES}{J}_{I0}{J}_{IS}{J}_{E0}\right){r}_{0}+{J}_{ES}{I}_{opto}}{{J}_{ES}{J}_{II}{J}_{EI}{J}_{IS}}$ and ${r}_{S}=\frac{\left({J}_{II}{J}_{E0}{J}_{EI}{J}_{E0}\right){r}_{0}{J}_{EI}{I}_{opto}}{{J}_{ES}{J}_{II}{J}_{EI}{J}_{IS}}$ .
Provided that the parameters are such that they are positive, ${\nu}_{E}$ is given by
Finally, since ${\nu}_{V}=o\left(1\right)$ consistency implies that
This solution is valid for a finite range of ${I}_{opto}$ . It exists as long as ${r}_{s}>0$ which implies that $J}_{E0}\frac{{J}_{II}}{{J}_{EI}}{J}_{I0}>{I}_{opto}>{I}_{opto}^{\text{*}$ .
4) ${\nu}_{E}=o\left(1\right)$ and ${\nu}_{V}=o\left(1\right)$
Here, ${u}_{S}={A}_{S}={B}_{S}=0$ and thus, ${r}_{S}=0$ . This solution exists only for sufficiently large ${I}_{opto}$ such that ${u}_{E}$ and ${u}_{V}$ are $O\left(\sqrt{K}\right)$ and negative. Therefore, PV is the only active population and ${r}_{I}$ is given by Equations (A40).
In conclusion, in this model at the minimum of ${r}_{I}$ , ${r}_{E}$ is of order $\frac{1}{\sqrt{K}}$ in contrast to the twopopulation case where ${r}_{E}$ is exponentially small in $K$ .
Fourpopulation model: Model 2
Large $K$ limit
To get insights on the network architecture that could explain the proportional paradoxical effect observed in layer 5 of ALM and S1, we first considered a threepopulation network consisting of the PC, PV and SOM populations (Figure 9A).
In this network, the population activities are
where $\mathrm{\Delta}=\left({J}_{II}{J}_{SE}{J}_{IE}{J}_{SI}\right){J}_{ES}+\left({J}_{EE}{J}_{SI}{J}_{EI}{J}_{SE}\right){J}_{IS}>0$ .
The full balance of the network activities implies
The inequality on the left side stems from the positivity of the rates. The inequality on the right side stems from the fact that the balanced state is the only solution of the dynamics, namely that no partially balanced solution (in particular, ${r}_{E}=0$ , ${r}_{I}=O\left(1\right)$ and ${r}_{S}=0$ and ${r}_{E}=0$ , ${r}_{I}=O\left(1\right)$ and ${r}_{S}=O\left(1\right)$ ) exists.
${r}_{E}$ and ${r}_{I}$ are proportional Equations (A49) and increase with ${I}_{opto}$ . As a consequence, the network never exhibits the paradoxical effect.
In this threepopulation network, the proportionality of ${r}_{E}$ and ${r}_{I}$ stems from the balance of inputs into the SOM population. To account for the proportional paradoxical effect, we consider a network model with an additional inhibitory population, denoted X (Figure 9B). Because in this network the SOM neurons only receive inputs from PCs and PV neurons, here, the balance of the SOM input also ensure the proportionality of ${r}_{E}$ and ${r}_{I}$ .
The susceptibilities upon PV stimulation are
where $\Delta =det\left(\left[{J}_{AB}{\u03f5}_{B}\right]\right)$ (see Material and methods).
Paradoxicality implies that
The susceptibilities upon PC stimulation are
Therefore, the PC population activity increases upon PC stimulation if
One can find a range of parameters (e.g. Figure 9C) such that:
The relative decrease in the SOM population is larger than that in the E and I populations. As a consequence, as ${I}_{opto}$ is increased, ${r}_{S}$ approaches zero when the PC and PV activities are still finite.
As ${I}_{opto}$ is increased further, the network settles into a partially balanced state where ${r}_{E}$ , ${r}_{I}$ and ${r}_{X}$ are finite and ${r}_{I}$ increases with ${I}_{opto}$ , while ${r}_{E}$ continues to decrease.
Thus, ${r}_{I}$ reaches its minimum value when ${r}_{E}$ is finite even in the large $K$ limit.
Data availability
All simulation, raw, and processed data software is available on GitHub (https://github.com/Amahrach/Paper4pop; copy archived at https://github.com/elifesciencespublications/Paper4pop).
References

Two dynamically distinct inhibitory networks in layer 4 of the neocortexJournal of Neurophysiology 90:2987–3000.https://doi.org/10.1152/jn.00283.2003

A canonical neural mechanism for behavioral variabilityNature Communications 8:15415.https://doi.org/10.1038/ncomms15415

BookNeuronal Dynamics: From Single Neurons to Networks and Models of CognitionCambridge University Press.https://doi.org/10.1017/CBO9781107447615

Mechanisms underlying a thalamocortical transformation during active tactile sensationPLOS Computational Biology 13:e1005576.https://doi.org/10.1371/journal.pcbi.1005576

ShortTerm plasticity explains irregular persistent activity in working memory tasksJournal of Neuroscience 33:133–149.https://doi.org/10.1523/JNEUROSCI.345512.2013

Synchronization and computation in a chaotic neural networkPhysical Review Letters 68:718–721.https://doi.org/10.1103/PhysRevLett.68.718

The mechanism of orientation selectivity in primary visual cortex without a functional mapJournal of Neuroscience 32:4049–4064.https://doi.org/10.1523/JNEUROSCI.628411.2012

Asynchronous rate Chaos in spiking neuronal circuitsPLOS Computational Biology 11:e1004266.https://doi.org/10.1371/journal.pcbi.1004266

Organization of cortical and thalamic input to pyramidal neurons in mouse motor cortexJournal of Neuroscience 33:748–760.https://doi.org/10.1523/JNEUROSCI.433812.2013

DualChannel circuit mapping reveals sensorimotor convergence in the primary motor cortexJournal of Neuroscience 35:4418–4426.https://doi.org/10.1523/JNEUROSCI.374114.2015

The organization of two new cortical interneuronal circuitsNature Neuroscience 16:210–218.https://doi.org/10.1038/nn.3305

Transition to Chaos in random neuronal networksPhysical Review X 5:041030.https://doi.org/10.1103/PhysRevX.5.041030

Dynamics of encoding in a population of neuronsThe Journal of General Physiology 59:734–766.https://doi.org/10.1085/jgp.59.6.734

A disinhibitory circuit mediates motor integration in the somatosensory cortexNature Neuroscience 16:1662–1670.https://doi.org/10.1038/nn.3544

Response variability in balanced cortical networksNeural Computation 18:634–659.https://doi.org/10.1162/neco.2006.18.3.634

Tuned thalamic excitation is amplified by visual cortical circuitsNature Neuroscience 16:1315–1323.https://doi.org/10.1038/nn.3488

Inhibitory stabilization and visual coding in cortical circuits with multiple interneuron subtypesJournal of Neurophysiology 115:1399–1409.https://doi.org/10.1152/jn.00732.2015

Distinct subtypes of somatostatincontaining neocortical interneurons revealed in transgenic miceJournal of Neuroscience 26:5069–5082.https://doi.org/10.1523/JNEUROSCI.066106.2006

Interneurons of the neocortical inhibitory systemNature Reviews Neuroscience 5:793–807.https://doi.org/10.1038/nrn1519

Inhibitory circuits in cortical layer 5Frontiers in Neural Circuits 10:35.https://doi.org/10.3389/fncir.2016.00035

Diversity and connectivity of layer 5 SomatostatinExpressing interneurons in the mouse barrel cortexThe Journal of Neuroscience 38:1622–1633.https://doi.org/10.1523/JNEUROSCI.241517.2017

The spatial structure of correlated neuronal variabilityNature Neuroscience 20:107–114.https://doi.org/10.1038/nn.4433

A balanced memory networkPLOS Computational Biology 3:e141.https://doi.org/10.1371/journal.pcbi.0030141

On the distribution of firing rates in networks of cortical neuronsJournal of Neuroscience 31:16217–16226.https://doi.org/10.1523/JNEUROSCI.167711.2011

Three groups of interneurons account for nearly 100% of neocortical GABAergic neuronsDevelopmental Neurobiology 71:45–61.https://doi.org/10.1002/dneu.20853

Noise, neural codes and cortical organizationCurrent Opinion in Neurobiology 4:569–579.https://doi.org/10.1016/09594388(94)900590

Relating neuronal firing patterns to functional differentiation of cerebral cortexPLOS Computational Biology 5:e1000433.https://doi.org/10.1371/journal.pcbi.1000433

Neural mechanisms of movement planning: motor cortex and beyondCurrent Opinion in Neurobiology 49:33–41.https://doi.org/10.1016/j.conb.2017.10.023

Paradoxical effects of external modulation of inhibitory interneuronsThe Journal of Neuroscience 17:4382–4388.https://doi.org/10.1523/JNEUROSCI.171104382.1997

Chaotic balanced state in a model of cortical circuitsNeural Computation 10:1321–1371.https://doi.org/10.1162/089976698300017214

BookIrregular Activity in Large Networks of neuronsIn: Chow C. C, Gutkin B, Hansel D, Meunier C, Dalibard J, editors. Methods and Models in Neurophysics, 80. Elsevier. pp. 341–406.https://doi.org/10.1016/S09248099(05)800150

Dynamical models of cortical circuitsCurrent Opinion in Neurobiology 25:228–236.https://doi.org/10.1016/j.conb.2014.01.017
Article and author information
Author details
Funding
Agence Nationale de la Recherche (14NEUC000101)
 Carl van Vreeswijk
Agence Nationale de la Recherche (13BSV4001402)
 David Hansel
Agence Nationale de la Recherche (09SYSC00201)
 David Hansel
Helen Hay Whitney Foundation
 Nuo Li
Robert and Janice McNair Foundation
 Nuo Li
Alfred P. Sloan Foundation
 Nuo Li
National Institutes of Health (NS104781)
 Nuo Li
Pew Charitable Trusts
 Nuo Li
Simons Foundation (543005)
 Nuo Li
Agence Nationale de la Recherche (ANR17NEUC0005)
 Alexandre Mahrach
 Carl van Vreeswijk
 David Hansel
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Karel Svoboda for illuminating discussions and comments on the manuscript. We are also thankful to Ran Darshan and Tohar Yarden for discussions. DH thanks Svoboda’s lab. and Janelia Research Campus for their warm hospitality. This work was supported by ANR grants ANR14NEUC0001–01 (CvV and DH), ANR13BSV4001402 (DH, CvV), the ANR09SYSC002–01 (DH, CvV), ANR17NEUC0005 (DH, CvV, AM), the Janelia Research Campus visiting program (DH), the Helen Hay Whitney Foundation fellowship (NL), the Robert and Janice McNair Foundation (NL), Whitehall Foundation (NL), Alfred P Sloan Foundation (NL), Searle Scholars Program (NL), NIH NS104781 (NL), the Pew Charitable Trusts (NL), and Simons Collaboration on the Global Brain (#543005, NL). Work performed in the framework of the FranceIsrael Center for Neural Computation (CNRS/Hebrew University of Jerusalem).
Ethics
Animal experimentation: All procedures were in accordance with protocols approved by the Janelia Research Campus and Baylor College of Medicine Institutional Animal Care and Use Committee.
Copyright
© 2020, Mahrach et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 3,492
 views

 609
 downloads

 53
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Immunology and Inflammation
 Neuroscience
Acute retinal ischemia and ischemiareperfusion injury are the primary causes of retinal neural cell death and vision loss in retinal artery occlusion (RAO). The absence of an accurate mouse model for simulating the retinal ischemic process has hindered progress in developing neuroprotective agents for RAO. We developed a unilateral pterygopalatine ophthalmic artery occlusion (UPOAO) mouse model using silicone wire embolization combined with carotid artery ligation. The survival of retinal ganglion cells and visual function were evaluated to determine the duration of ischemia. Immunofluorescence staining, optical coherence tomography, and haematoxylin and eosin staining were utilized to assess changes in major neural cell classes and retinal structure degeneration at two reperfusion durations. Transcriptomics was employed to investigate alterations in the pathological process of UPOAO following ischemia and reperfusion, highlighting transcriptomic differences between UPOAO and other retinal ischemiareperfusion models. The UPOAO model successfully replicated the acute interruption of retinal blood supply observed in RAO. 60 min of Ischemia led to significant loss of major retinal neural cells and visual function impairment. Notable thinning of the inner retinal layer, especially the ganglion cell layer, was evident postUPOAO. Temporal transcriptome analysis revealed various pathophysiological processes related to immune cell migration, oxidative stress, and immune inflammation during the nonreperfusion and reperfusion periods. A pronounced increase in microglia within the retina and peripheral leukocytes accessing the retina was observed during reperfusion periods. Comparison of differentially expressed genes (DEGs) between the UPOAO and high intraocular pressure models revealed specific enrichments in lipid and steroid metabolismrelated genes in the UPOAO model. The UPOAO model emerges as a novel tool for screening pathogenic genes and promoting further therapeutic research in RAO.

 Computational and Systems Biology
 Neuroscience
Diffusional kurtosis imaging (DKI) is a methodology for measuring the extent of nonGaussian diffusion in biological tissue, which has shown great promise in clinical diagnosis, treatment planning, and monitoring of many neurological diseases and disorders. However, robust, fast, and accurate estimation of kurtosis from clinically feasible data acquisitions remains a challenge. In this study, we first outline a new accurate approach of estimating mean kurtosis via the subdiffusion mathematical framework. Crucially, this extension of the conventional DKI overcomes the limitation on the maximum bvalue of the latter. Kurtosis and diffusivity can now be simply computed as functions of the subdiffusion model parameters. Second, we propose a new fast and robust fitting procedure to estimate the subdiffusion model parameters using two diffusion times without increasing acquisition time as for the conventional DKI. Third, our subdiffusionbased kurtosis mapping method is evaluated using both simulations and the Connectome 1.0 human brain data. Exquisite tissue contrast is achieved even when the diffusion encoded data is collected in only minutes. In summary, our findings suggest robust, fast, and accurate estimation of mean kurtosis can be realised within a clinically feasible diffusionweighted magnetic resonance imaging data acquisition time.