Beyond excitation/inhibition imbalance in multidimensional models of neural circuit changes in brain disorders
Abstract
A leading theory holds that neurodevelopmental brain disorders arise from imbalances in excitatory and inhibitory (E/I) brain circuitry. However, it is unclear whether this onedimensional model is rich enough to capture the multiple neural circuit alterations underlying brain disorders. Here, we combined computational simulations with analysis of in vivo twophoton Ca^{2+} imaging data from somatosensory cortex of Fmr1 knockout (KO) mice, a model of FragileX Syndrome, to test the E/I imbalance theory. We found that: (1) The E/I imbalance model cannot account for joint alterations in the observed neural firing rates and correlations; (2) Neural circuit function is vastly more sensitive to changes in some cellular components over others; (3) The direction of circuit alterations in Fmr1 KO mice changes across development. These findings suggest that the basic E/I imbalance model should be updated to higher dimensional models that can better capture the multidimensional computational functions of neural circuits.
https://doi.org/10.7554/eLife.26724.001eLife digest
In many brain disorders, from autism to schizophrenia, the anatomy of the brain appears remarkably unchanged. This implies that the problem may reside in how neurons communicate with one another. Unfortunately, neuroscientists know little about how brain activity might differ from normal in these disorders, or how specific changes in activity give rise to symptoms. One leading theory, first proposed over a decade ago, is that these disorders reflect an imbalance in the activity of excitatory and inhibitory neurons. Excitatory neurons activate their targets, whereas inhibitory neurons suppress or silence them. While studies in mice have lent support to this theory, they have not yet culminated in new treatments for brain disorders.
One limitation of the excitationinhibition imbalance theory is that it is onedimensional. It assumes that there is an optimal balance of excitation and inhibition, and that brain disorders can be arranged in an imaginary line on either side of this optimum. Disorders to the right of the optimum, such as epilepsy and some forms of autism, feature too much excitation. Disorders to the left, such as the developmental disorder Rett syndrome, feature too much inhibition. But can diverse brain disorders really be classified on the basis of a single property, or do scientists need to consider other factors?
To find out, O’Donnell et al. analyzed recordings of brain activity from genetically modified mice with the mutation that causes fragile X syndrome, the most common form of inherited learning disability and autism. The mice showed changes in their overall brain activity compared to control animals. Their neurons also tended to fire in a more synchronized manner. A computer simulation revealed that an imbalance in excitation and inhibition alone could not explain these changes. Yet, a more complex simulation incorporating extra properties of neural circuits did a better job of explaining the altered neural activity seen in the mice.
O’Donnell et al. propose that this more advanced multidimensional model of changes in neural circuits could be used to screen candidate drugs before testing them in patients. In principle, the model could even help with designing drugs or other interventions by making it easier for researchers to target more precisely the changes in neural circuits that occur in brain disorders.
https://doi.org/10.7554/eLife.26724.002Introduction
The nervous system shows complex organization at many spatial scales: from genes and molecules, to cells and synapses, to neural circuits. Ultimately, the electrical and chemical signaling at all these levels must give rise to the behavioral and cognitive processes seen at the wholeorganism level. When trying to understand prevalent brain disorders such as autism and schizophrenia, a natural question to ask is: where is the most productive level of neuroscientific investigation? Traditionally, most major disorders are diagnosed entirely at the behavioral level, whereas pharmaceutical interventions are targeted at correcting alterations at the molecular level. However, even for the most successful drugs, we have little understanding of how pharmaceutical actions at the molecular level percolate up the organizational ladder to affect behavior and cognition. This classic bottomup approach may even be further confounded if phenotypic heterogeneity in disorders such as autism turn out not to reflect a unique cellular pathology, but rather ‘a perturbation of the network properties that emerge when neurons interact’ (Belmonte et al., 2004). These considerations imply that a more promising level of analysis might be at the level of neural circuits, since the explanatory gap between circuits and behavior is smaller than the gap between molecules and behavior. This circuitlevel viewpoint argues for a reverseengineering approach to tackling brain disorders: rather than start at the molecular level and working up, we should instead start by asking how cognitive and behavioral symptoms manifest as alterations at the circuit level, then interpret these changes at the levels of cells, synapses, and molecules as appropriate.
One prominent circuitlevel hypothesis for brain disorders has been the idea of an imbalance in excitatory and inhibitory signaling. First proposed as a model for autism (Rubenstein and Merzenich, 2003), the concept has since been applied to many other brain disorders, including Schizophrenia, Rett syndrome, fragileX syndrome, tuberous sclerosis, and Angelman Syndrome. However, a major drawback of this model is that it only considers overall activity, which is onedimensional. It implies that either too much excitation or too much inhibition is unhealthy (Figure 1A). Although several studies have found evidence that the E/I balance is indeed upset in multiple brain disorders (Bateup et al., 2011; Dani et al., 2005; Gibson et al., 2008; Kehrer et al., 2008; Wallace et al., 2012), a model’s usefulness should not be judged on whether it is nominally true or false, but on its explanatory and predictive powers as compared with competing alternative models. In this study, we argue that that even if the E/I imbalance model proves correct, its unidimensionality might ultimately limit its applicability, for three reasons.
First, by placing all disorders on the same single axis, the E/I imbalance model implicitly lumps together some vastly different disorders, such as epilepsy, schizophrenia and autism (Figure 1A) because they share an excess of excitation. By extension it implies that the symptoms of diverse disorders could be normalized solely by either enhancing or reducing the level of, say, GABAergic signaling as appropriate. Although clinical trials for such GABAergicbased interventions are ongoing (Braat and Kooy, 2015), no treatment for a neurodevelopmental disorder based on this principle has yet been approved.
A second issue with the unidimensionality of the E/I imbalance model is that it lumps together all excitatory and inhibitory neural circuit components. In Figure 1B, we show a schematic diagram of a generic neural circuit with excitatory components colored red and inhibitory components colored blue. The E/I imbalance model implies that varying any of the excitatory components, such as the strength of recurrent excitatory synapses or the input resistances of excitatory neurons, would have the same overall effect on circuit function. In contrast, theorists have found that these equivalences often do not hold even in very simple circuit models (Wilson and Cowan, 1972).
Third, because the standard E/I imbalance model is given in terms of circuit components, not circuit function, it does not specify which aspect of a neural circuit’s activity should be maintained for healthy performance. For example, it leaves unclear which of neuronal firing rates, synchrony, or reliability of responses might be altered if E/I balance is upset.
To motivate our study, we began by investigating which circuit activity properties are altered in a model brain disorder. We reanalyzed published in vivo twophoton Ca^{2+} imaging data we previously recorded from somatosensory cortex in Fmr1 knockout mice (Gonçalves et al., 2013), a wellstudied animal model for fragileX syndrome (The DutchBelgian Fragile X Consortium, 1994). We compared the data from wildtype (WT) mice with Fmr1 KO mice, across three different developmental time points: just before (P9–11) and after (P14–16) the critical period for heightened activitydependent synaptic plasticity in L2/3 barrel cortex, and a more mature timepoint (P30–40). Example ΔF/F raster plots from each group are shown in Figure 1C, top left. We binned the data into 1 s timebins (originally imaged at 4 Hz), then transformed each neuron’s timeseries of ΔF/F values into a probabilistic sequence of binary ON/OFF values by assuming a Poisson firing model (Materials and methods). We then summarized the neural population activity from each animal with three statistics: the mean ON probability across all recorded neurons, the standard deviation (s.d.) in ON probability across neurons, and the mean correlation between all pairs of neurons (Figure 1C, bar charts right and scatter plots lower left). Together these measures capture both the statistics of the bulk population activity and some indication of the heterogeneity across neurons.
For mean firing rates, the only change we detected was a decrease in firing probability in WT between P14–16 and P30–40 (p=0.027), which was coupled with an increased s.d. of firing rates (p=0.015). We also detected a higher firing rate s.d. in P9–11 KO animals than WT (p=0.031). Finally, as previously reported (Golshani et al., 2009; Gonçalves et al., 2013; Rochefort et al., 2009), we found a substantial decrease in pairwise correlations in both genotypes across development, with slightly higher correlations in KO animals than WT at P9–11 (p=0.029) and P14–16 (p=0.047).
These results show that multiple statistics of cortical circuit activity are altered in Fmr1 KO mice. However, two questions remain: (1) Which circuit components are responsible for these activity alterations? (2) What aspects of these activity alterations impact circuit computation? In the remainder of this study, we used computational simulations and further data analysis to ask whether the E/I imbalance model could help address these questions.
We first built a detailed spiking neural circuit model of mouse L2/3 somatosensory cortex to explore how its various excitatory and inhibitory components affect the circuit’s spiking output, and found that the E/I imbalance model was not flexible enough to capture key aspects of this relationship. We then derived an abstract twodimensional circuit model that captures more features of the circuit function than the 1dimensional E/I imbalance model did. Using this new model, we found that certain sets of circuit components have redundant effects on circuit function, and that circuit function is vastly more sensitive to changes in some components over others. To ask how this 2D model could help interpret brain circuit abnormalities in a particular test case, we fit a version of the model to the Ca^{2+} imaging data from FragileX mouse models presented above. We found that the model predicts opposite changes in Fmr1 KO circuit properties at different developmental ages. Finally, we applied a new largescale neural population analysis method (O'Donnell, 2017b) to the same Ca^{2+} imaging data, and found systematic shifts in the distribution of neural circuit activity patterns in FragileX that were not predictable from neural firing rates or correlations alone.
Results
Neural circuits consist of many components that typically interact nonlinearly to generate complex circuit activity dynamics. Although many properties of cortical circuit components have been found to be altered in animal models of brain disorders, it remains extremely difficult to predict the net effect of varying any one particular parameter on circuit activity. The E/I imbalance model seeks to simplify this problem by projecting all circuit alterations onto a onedimensional axis (Figure 1A) where the goal is to achieve a ‘healthy’ balance of excitation and inhibition. Under this model, either too much excitation or too much inhibition leads to improper circuit function.
To explicitly test whether the E/I imbalance model can account for the effects of cellular component changes on circuit function, we built a detailed computational model of layer (L) 2/3 mouse somatosensory cortex. This circuit has been studied in detail by neurophysiologists, and several of its properties have been found to be altered during development in mouse models of FragileX syndrome, including parvalbuminpositive interneuron number (Selby et al., 2007), GABA receptor reversal potential (He et al., 2014), dendritic spine dynamics (CruzMartín et al., 2010), and L4 excitatory input (Bureau et al., 2008), reviewed by Contractor et al., 2015. We used numerical computer simulations to perform hypothetical experiments where we perturbed various parameters of the circuit model and observed the resulting changes in circuitlevel activity. Although we focused on this particular brain circuit for tractability, our general conclusions and methodology should be readily applicable to other brain circuits (Frye and MacLean, 2016).
The L2/3 computational model we built (Figure 2A, see Materials and methods for details) consisted of four interconnected populations of leaky integrateandfire neurons: one group of 1700 excitatory (E) pyramidal neurons and three groups of inhibitory neurons: 115 5HT_{3A}Rexpressing neurons, 70 parvalbuminexpressing (PV) neurons, and 45 somatostatinexpressing (SOM) neurons. This L2/3 circuit was driven by a separate population of 1500 L4 excitatory neurons. Cellular numbers, intrinsic properties, synaptic strengths, and connectivity statistics were taken from published in vitro data from P17–22 wildtype mice (Avermann et al., 2012; Lefort et al., 2009; Tomm et al., 2014). We chose this level of detail for the model in order to relate experimentally measureable biophysical properties of neurons to their putative role in the circuit at large.
L2/3 neurons of the rodent somatosensory neocortex respond only sparsely to sensory stimulation in vivo. For example, twitching a whisker activates, on average, only ~20% of L2/3 neurons in its corresponding barrel, each of which typically emits only one action potential (Clancy et al., 2015; Kerr et al., 2007; Sato et al., 2007). Hence any individual neuron carries very little information about the stimulus on its own, implying that information must instead be encoded at the circuit level as the identities of the subset of neurons that respond.
To model whisker stimulation, we simulated a volley of spikes arriving from L4 as input to the population of L2/3 cells. We chose a random subset of L4 neurons as ON, then sent a single spike from each of these L4 cells to their target neurons in L2/3, and recorded the responses of all neurons in L2/3, some of which spiked and some of which did not. We repeated this identical stimulation multiple times, in order to get an average response probability for each L2/3 neuron, given the probabilistic vesicle release at synapses in the model. Then, we chose a different random subset of L4 neurons as ON, and repeated the entire procedure. Finally, we varied the fraction of L4 cells active and plotted the probability of response for each individual L2/3 neuron as a function of L4 activity level (Figure 2B–C).
We found that the mean response probability of each neuron in the simulation increased from zero to one monotonically with increasing L4 activity level. This sigmoidalshaped response profile of simulated L2/3 neurons mimics the spiking response of mouse L2/3 pyramidal cells to extracellular L4 stimulation in vitro (Elstrott et al., 2014), while the sparse, noisy and distributed network responses were reminiscent of in vivo activity following whisker stimulation (Clancy et al., 2015; Kerr et al., 2007). Neurons of all four cell types in the simulation responded to the L4 stimulus, including the SOM interneurons which did not receive direct L4 input, but were instead activated by disynaptic connections via L2/3 excitatory neurons. The detailed shape of the response curve varied systematically across cell types and was heterogeneous for different neurons of a given cell type. To quantify these differences, we used logistic regression to fit the response profile of each neuron with a sigmoid function (Figure 2D, inset), which has just two parameters: the slope (representing the steepness of the response curve) and threshold (representing the minimal fraction of L4 neurons needed to activate the cell). When we plotted the fitted slope and threshold values for each neuron against each other, we found that each cell type falls into a distinct cluster in this twodimensional space. For example, all PV inhibitory neurons had a low slope and low threshold, whereas SOM inhibitory neurons had a steep slope and moderate threshold. We then used these slopethreshold measurements to summarize the circuitlevel inputoutput function of this ‘default’ model of L2/3 somatosensory cortex. This 2D logistic model has two benefits over the 1D E/I imbalance model: first, its extra degree of freedom allows for richer and more flexible fits to data, and second, by describing an inputoutput mapping for the L2/3 circuit it can capture some aspects of the computation that the circuit performs for the animal. In contrast, the E/I imbalance model is specified purely in terms of circuit components, and so is agnostic to the circuit’s computational function.
The biophysical circuit model contained 100 parameters (Materials and methods). How sensitive is the circuit’s macroscopic inputoutput function to alterations in its lowlevel components? To test this, we varied 76 of the model parameters in turn by ±20%, and repeated the entire set of simulations for each case (Materials and methods). For each individual parameter alteration, we fit a new logistic response function for each model neuron. We summarize the results by plotting the logistic slope and threshold parameters and comparing their values to those found with the default model. The outcomes were hugely varied. We show three examples from the set of 76 in Figure 3A, chosen to illustrate three qualitatively different effects that neural parameter changes can have on circuit function. First, when we increased the amplitude of postsynaptic potentials (PSPs) of excitatory synapses from L4 to L2/3 excitatory neurons, we found that the logistic threshold parameter of all cell types shifted leftwards to lower values (Figure 3A left), implying that fewer L4 neurons were needed to activate the entire L2/3 circuit. In contrast, when we increased the PSP amplitude of a different excitatory synapse, the recurrent connections between L2/3 excitatory neurons, we found (Figure 3A center) that excitatory and SOM inhibitory neurons had increased slope parameters relative to default, with little change in their threshold parameters. At the same time, 5HT_{3A}R inhibitory neurons had decreased slopes and thresholds, while PV neurons had little change at all. As a third example, we increased the probability of inhibitory synaptic connections from L2/3 PV interneurons to L2/3 excitatory neurons (Figure 3A right). In this case, we found that excitatory neurons had a marginally lower slope and increased threshold, SOM inhibitory neurons had a lower slope, 5HT_{3A}R inhibitory neurons had both an increased slope and threshold, and PV neurons showed little change, even though their outgoing synapses were the parameter that was altered.
To synthesize the findings from all simulations, for each altered parameter we plotted the shift in mean slopethreshold values for L2/3 excitatory neurons from the mean values found with the default model (Figure 3B). We focused on excitatory neurons because they constitute 90% of the neurons in this layer (Lefort et al., 2009) and are the primary output to downstream circuits (Mao et al., 2011; Petreanu et al., 2007). Overall, we found a very heterogeneous picture. First, the magnitude of the shift in circuit response varied from parameter to parameter (Figure 3B–C). Varying some parameters, such as the first two examples given above, had large effects, whereas varying other parameters such as w_{I5htE} (the strength of synapses from 5HT_{3A}R inhibitory neurons to E neurons) or τm_{Isom} (the membrane time constant of SOM inhibitory neurons) had little effect. Second, the direction of shift in 2D slopethreshold space also depended on parameter (Figure 3C). Increasing some parameters changed either circuit slope or threshold in isolation, while other parameters changed both slope and threshold together. All four quadrants of the slopethreshold plane could be reached by various subsets of the model parameters. Third, basic knowledge of whether a component was ‘excitatory’ or ‘inhibitory’ was insufficient to predict the direction of slopethreshold change. For example, the two glutamatergic projections considered in Figure 3A had distinct effects on circuit function.
In summary, these simulations indicate that the L2/3 somatosensory cortex circuit has extremely varied sensitivities to changes in its cellular components and that the eventual circuitlevel consequences cannot be predicted from knowledge of the class of the perturbed neurotransmitter alone. Since the E/I imbalance model groups all excitatory and inhibitory components as respective equals, it cannot account for these results.
Firing rates and correlations from the logistic model
In the above analysis, we investigated how lowlevel circuit components affect a highlevel circuit inputoutput function, as parameterized by the slope and threshold of fitted logistic functions. But how is this logistic inputoutput function related to more common measures of neural population activity, such as firing rates and pairwise correlations between neurons? To investigate this, we considered the following reduced statistical model of cortical activity. We assumed for simplicity that the magnitude of the total input to the L2/3 circuit can be described by a Gaussian distributed random variable, with zero mean and unit standard deviation (Figure 4A lower left). Then we described each L2/3 neuron’s inputoutput as a logistic function as before (Figure 4A upper left), with threshold and slope defined relative to the Gaussian input’s mean and standard deviation, respectively. Given this model, we can numerically calculate the probability distribution over a neuron’s firing probability, which in general is skewed and nonGaussian (Figure 4A upper right). From this function, we compute (Materials and methods) both the neuron’s mean firing probability and the pairwise correlation of two identical neurons following this profile (Figure 4A lower right). Example samples from the model are illustrated in Figure 4B.
Neural firing rates and correlations had qualitatively different dependencies on the underlying logistic model’s slope and threshold. Neural firing rate was greatest when threshold was low and slope was high (top left of phase plot, Figure 4C left), whereas correlations were greatest when both threshold and slope were high (top right of phase plot, Figure 4C right). This implies that any change in the circuit’s inputoutput function slope or threshold will in general have distinct effects on firing rate versus correlations, and so could not be captured by a 1dimensional E/I balance model that sought to, for example, normalize firing rates alone. To illustrate this fact, we plot the calculated correlation values along a contour where firing probability is fixed at 0.1 (Figure 4D). In the region of parameter space where both the slope and threshold are low (Figure 4C bottom left), correlations are low,~0.01. However, as we move along the contour for firing rate = 0.1 toward the region of parameter space where slope and threshold are high (Figure 4C top right), the pairwise correlations increase to ~0.4. This shows that a onedimensional E/I balance rule that exclusively sought to normalize neural firing rates would leave neural correlations free to achieve arbitrary values.
Previous studies have found evidence for an E/I imbalance in ASD (Lee et al., 2017; Nelson and Valakh, 2015). FragileX syndrome is the leading inherited cause of ASD, and also carries alterations in excitability (Contractor et al., 2015). We aimed to interpret our Fmr1 KO Ca^{2+} imaging data (Figure 1C) via the 2D logistic model. Since our earlier analysis found substantial withinanimal heterogeneity in neural activity levels in both genotypes (Figure 1C), we extended the 2D logistic model for single neurons to a 5D neural population version that captured celltocell heterogeneity. The three extra parameters represented the standard deviations and correlation in slope and threshold parameters across the neural population (see Materials and methods for details). We fit the parameters of the logistic model to reproduce the same neural population Ca^{2+} imaging data presented in Figure 1C. Given the three summary statistics from each animal (firing rate mean and s.d., and mean pairwise correlation), we used a gradient descent algorithm to find the five parameters of the populationlevel version of the logistic model that best matched the activity statistics (see Materials and methods). The output statistics of the fitted models matched well those of the target data (Figure 5–figure supplement 1). Example neural population activity patterns drawn from the mean model fits for each group are shown in Figure 5A, along with the fitted slopethreshold functions (Figure 5A insets), to be compared with the Ca^{2+} imaging rasters in Figure 1C. We also plot the full 5D parameter fits for all animals in Figure 5–figure supplement 2. For the rest of the analysis, we focus on the mean slope and mean threshold parameters, which showed the most prominent changes. In Figure 5B, we plot the mean slope and mean threshold fits on top of the previously calculated (Figure 4C) 2D slopethreshold maps of firing rate and correlation. We found that in young animals, P9–11, most points were scattered at high values of both slope and threshold (Figure 5B left). With age, the parameter fits for both genotypes moved southwest toward the low slope and low threshold region of parameter space (Figure 5B center and right). The mean location of the cloud of points at each developmental age differed between WT and KO. We plot the direction of shift in group mean from WT to KO in Figure 5C. In young animals, P9–11 and P14–16, the KO group had both higher slope and higher threshold than WT, whereas in adult animals, P30–40, the KO group had a lower slope and lower threshold than WT. These results demonstrate an opposite direction of circuit parameter change in young FragileX mice compared to adults, which was not be uncovered by measures of neural firing rates and correlations (Figure 1C).
Earlier we asked how sensitive the logistic model slope and threshold parameters were to alterations in the many underlying neural circuit components (Figure 3). In a similar way, we can also ask how sensitive the neural firing rates and correlations are to alterations in the logistic slope and threshold parameters. This is important since inspection of the twodimensional maps in Figure 3C shows that these sensitivities will differ depending on starting location within the slopethreshold space. To quantify this effect, we calculated the sensitivity of both the firing rate and correlations to small changes in the slope and threshold (Figure 6A–B, see Materials and methods), quantified as the partial derivatives local to the fitted logistic parameter values for each animal (black and red circles in Figure 5B). In general, increasing the slope or decreasing the threshold always increased both firing rates and correlations, as can be predicted from Figure 5B. However, the magnitude of sensitivities varied across animals. We found only minor differences in sensitivities between genotypes (Figure 6C–D), and as a result we pooled the sensitivity measurements across genotypes to test for statistical differences in sensitivity with developmental age. In young animals, P9–11, changes in the logistic threshold (solid bars in Figure 6) had substantial effects on both firing rates and correlation. This sensitivity decreased with age (p≤0.013 for firing rates, p<0.01 for correlations from P9–11 to P14–16), so that in adult animals, P3040, changes in threshold had relatively little effect on neural activity statistics. A different picture emerged for the logistic slope parameter (striped bars in Figure 6). There, the firing rate sensitivity increased from P9–11 to P14–16 (p<1e6) (Figure 6C), while correlation sensitivity stayed approximately constant (p≥0.18) (Figure 6D). These results show that the quantitative relationships between neural activity statistics and the parameters of the logistic model, and perhaps also the underlying circuit components, are not fixed across development.
What are the functional implications of these alterations in firing rates and correlations in FragileX mice across development? To address this, we calculated the entropy of the neural population activity for the data from each animal. Entropy is a quantity from information theory, measured in bits, that puts a hard upper bound on the amount of information that can be represented by any coding system (Cover and Thomas, 2006). Intuitively, the entropy measures how uniform the neural population activity pattern distribution is: it is large if the circuit exhibits many different activity patterns over time, and small if only a few activity patterns dominate. Entropy is an appealing measure for the present problem because it is sensitive both to neural firing rates and to correlations at all orders. It is typically highest when firing rates are high and correlations are low. Although entropy is notoriously difficult to calculate for large neural populations because most estimation methods require impractically long data recordings (Quian Quiroga and Panzeri, 2009), we recently developed a new statistical method for this purpose, called the population tracking model, that scales well to large numbers of neurons, even for limited data (O'Donnell, 2017b). This model matches both the synchrony distribution for the number of neurons simultaneously active, and the variations in individual celltocell firing rates. We fit this population tracking model to the same Ca^{2+} imaging data as analyzed above (Figure 7). An intermediate step in estimating the neural entropy involves calculating a lowparameter approximation of the entire probability distribution over all 2^{N} neural population activity patterns, where N is the number of neurons. The cumulatives of these probability distributions calculated for 50neuron subsets of the recordings are shown in Figure 7A. In young animals P9–11, a small number of activity patterns accounts for a large fraction of the probability mass (Figure 7A left). For example, based on these curves, 50% of the time we would expect to see the same 1000–10,000 patterns out of a possible total 2^{50} ≈ 10^{15} patterns. In contrast, in older animals P14–16 and P30–40 the cumulative distributions shift rightwards so that more patterns are typically observed (Figure 7A center, right). In these cases, around 1,000,000 patterns are needed to account for 50% probability mass.
Instead of attempting to quantify these shifts by asking how many patterns are needed to cross an arbitrary threshold of probability mass, we instead calculated the entropy, which takes into account the shape of the entire probability distribution. The entropy depends on the number of neurons analyzed, so we normalized all estimates to calculate the entropy per neuron (Figure 7B–C). Since we are treating neurons as binary, the entropy/neuron was bounded between 0 and 1 bits. For all age groups, and for both WT and Fmr1 KO animals, entropy/neuron progressively decreased with the number of neurons analyzed (Figure 7B). Because each imaging session captured a different number of neurons (range 40–198, median 97), we fit the entropy/neuron versus number of neurons data with a double exponential function (see Materials and methods) and use the fit to provide a standardized estimate of the entropy/neuron for 100neuron populations (Figure 7C). In WT animals, entropy/neuron showed a nonmonotonic trajectory across development (O'Donnell, 2017b). At P9–11 it was low, 0.38 bits (95% c.i. [0.35:0.41]), before increasing at P14–16 (p<0.001) to 0.50 bits (95% c.i. [0.48:0.52]), before decreasing again at P30–40 (p=0.028) to 0.45 bits (95% c.i. [0.42:0.48]). We found a different entropy trajectory in Fmr1 KO animals. There, although entropy/neuron also began low at 0.34 bits (95% c.i. [0.30:0.39]), not different from WT (p=0.19), when it increased at P14–16 (p<0.001) to 0.465 bits (95% c.i. [0.45:0.48]) it remained lower than for WT (p=0.048). Finally, instead of decreasing as in the WT case, entropy continued to increase in P30–40 Fmr1 KO animals (p=0.033) to 0.51 bits (95% c.i. [0.47:0.55]), higher than WT (p=0.034). These entropy values estimated directly from Ca^{2+} imaging data agreed well with entropy estimates for synthetic data sampled from the previously fit logistic models (Figure 7–figure supplement 1). In summary, unlike WT animals, Fmr1 KO mice showed a monotonically increasing entropy/neuron from P9–11 to P30–40. Furthermore, the direction of change in entropy between P14–16 and P30–40 was opposite for WT and Fmr1 KO animals, decreasing in the former and increasing in the latter.
Discussion
The onedimensional E/I imbalance model has been widely used for interpreting neural circuit changes observed in animal models of diverse brain disorders (Bateup et al., 2011; Dani et al., 2005; Gibson et al., 2008; Kehrer et al., 2008; Wallace et al., 2012). In the case of FragileX syndrome, the hyperexcitability prediction of the E/I imbalance model is consistent with many of the symptoms of the disease (e.g. seizures, hyperarousal, hyperactivity, hypersensitivity to sensory stimuli) and the known pathogenic defects implicated in Fmr1 KO mice (diminished GABA signaling, exaggerated intrinsic excitability, increased neuronal firing rates; reviewed by Contractor et al., 2015). Here, we tested the hypothesis that the E/I imbalance model can account for alterations in other neural activity statistics beyond the mean firing rates; however, our results demonstrated that it was inadequate. The model was too inflexible to account for the joint alterations in both neural firing rates and correlations observed in FragileX model mice. This suggests that future studies of brain disorders may need to consider higherdimensional models of neural circuit dysfunction.
To test how cellular components affect their circuit function, we built computational models of mouse L2/3 somatosensory cortex at two levels of abstraction: a detailed, 100parameter biophysical model, and a twoparameter logistic response model. The purpose of the detailed model was to build a representation of the circuit where each parameter has a onetoone mapping with something that could be experimentally measured in a real animal  indeed many of these parameters have been shown to be altered in Fmr1 KO mice. The purpose of the logistic model was different: it simple enough to be both derivable from the complex model, and provide a direct link with measurable activity variables in our in vivo Ca2 +imaging data, firing rates and correlations. The disadvantages of detailed models are that they contain many parameters, and so are hard to constrain to data – in this case it was only possible because of the large dataset from Petersen et al., for P1722 WT mice (Avermann et al., 2012; Lefort et al., 2009; Tomm et al., 2014). The disadvantages of the simple 2D model is that its logistic inputoutput structure implies a very strong and specific assumption about the functional purpose of the circuit – to generate single spikes across a subset of neurons. Although this may be a physiologically relevant computation for this particular brain circuit (Clancy et al., 2015; Kerr et al., 2007; Sato et al., 2007), it is not immediately obvious how to extend this approach to include temporal correlations, for example, or to apply it to other brain circuits where we may have less insight into their natural computations. Nevertheless, our approach demonstrates a new way to tackle such problems.
After building the detailed computational model of L2/3 of mouse somatosensory cortex (Figure 2), we asked how sensitive the spiking responses of the overall circuit were to changes in its underlying neural components, many of which are known to be altered in Fmr1 KO mice (Bureau et al., 2008; Gibson et al., 2008; Gonçalves et al., 2013; Harlow et al., 2010; Hays et al., 2011; Paluszkiewicz et al., 2011; Patel et al., 2013; TestaSilva et al., 2012). We found that while changing some neural parameters did have a large effect, changing other parameters had little or no effect on circuit function (Figure 3B). This redundancy property has been reported as widely prevalent in computational models of biological systems (Gutenkunst et al., 2007; O'Leary et al., 2015). Its existence has two important implications for studies of brain disorders: first, many of the physiological component changes discovered in animal models may be entirely benign at the circuit level. Second, any treatment designed to correct circuit function is free to push the system by arbitrary amounts along insensitive directions in parameter space without consequence, as long as it makes the correct perturbations along the sensitive directions. The insensitive directions form a null space, which is a subspace of the parameter space.
An important caveat to our parameter sensitivity analysis is that it was linear and local to a particular point in the highdimensional model parameter space, corresponding to WT P1722 mice. Since the circuit dynamics are nonlinear, it is likely that the particular parameter sensitivities would be different in other parts of parameter space, especially near bifurcations where qualitatively different dynamics emerge (Hirsch et al., 2013). However, as long as the redundancy property is widely preserved, as suggested by studies on computational models of other biological systems (Fisher et al., 2013; Gutenkunst et al., 2007; Machta et al., 2013; Panas et al., 2015), then our conclusions for brain disorders remain valid.
In addition to the varying magnitudes of circuit components’ effect on circuit function, we also found that different components shifted the circuit inputoutput function in different directions, as defined by our 2D logistic response model (Figures 1 and 2). Even circuit parameters that are nominally of the same type, such as the strength of glutamatergic synapses between excitatory (E) neurons in L4 to E neurons in L2/3 or synaptic strength between E neurons within L2/3, had qualitatively different effects on the circuit response to stimulation (Figure 3). According to the standard E/I imbalance model (Rubenstein and Merzenich, 2003), both of these parameters should have similar effects on circuit function; but according to the logistic response model we studied, their differing effects on slope and threshold parameters must necessarily lead to different magnitudes of change in neural firing rates and correlations (Figure 4C). Indeed, no 1dimensional model of circuit function could ever capture the heterogeneity in parameter sensitivities that we observed (Figure 3B).
Next, we fit the parameters of the logistic response model to match the in vivo firing statistics of neural populations from WT and Fmr1 KO mice of varying age (Figure 5). Previous studies had found that neural correlations decrease during development (Golshani et al., 2009; Rochefort et al., 2009), and that early postnatal Fmr1 KO mice had higher correlations and firing rates than WT mice (Gonçalves et al., 2013; La Fata et al., 2014). Circuit hypersynchrony may be a general defect in autism disorders, as it is also found in mouse models of Rett syndrome (Lu et al., 2016). However, the relationship between these changes in firing statistics and the underlying neural circuit components were unclear. Our logistic model helps bridge this gap, leading to two findings: first, the direction of circuit parameter change from WT to KO was opposite in young (P9–11 and P14–16) versus mature (P30–40) animals (Figure 5C). Similar opposing switches in sensory cortex properties with age were also recently reported in Fmr1 KO and WT rats (Berzhanskaya et al., 2016). Second, we found that the sensitivity of neural firing rates and correlations to changes in underlying circuit components depends on developmental age (Figure 6). Taken together, these findings imply that qualitatively different interventions may be needed at different stages of development in FragileX, and perhaps other neurodevelopmental disorders, to shift cortical circuit function towards typical wildtype operation.
Spontaneous, intrinsic activity is ubiquitously present in mammalian cerebral cortex. It is highly structured at multiple spatiotemporal scales (Mitra et al., 2015; Ringach, 2009) and interacts strongly with the signals evoked by sensory stimulation (Ringach, 2009). Cellularresolution recordings in animals have shown that the patterns of spontaneous activity in neural populations are representative of the ensemble of activity patterns used by the brain to represent sensory stimuli (Berkes et al., 2011; Luczak et al., 2009; Miller et al., 2014). Here, we found that the entropy of spontaneous activity in WT mouse somatosensory cortex follows an invertedU shaped trajectory across development, and that this trajectory is dramatically altered in the Fmr1 KO mouse model of FragileX (Figure 7). Although we saw no reliable differences across genotypes in early postnatal animals (P9–11), Fmr1 KO animals showed lower entropy than WT after the second postnatal week (P14–16), while surprisingly switching to show higher entropy than WT in adult (P30–40). Notably, this switch in the direction of entropy change from WT to KO during development mirrors the reversing we saw in logistic model parameter changes in Figure 5C. Together, these findings suggest a perturbed trajectory of cortical development during the critical period in Fmr1 KO mice (Meredith et al., 2012). However, our results cannot distinguish whether the observed perturbation in L2/3 activity statistics reflects a developmental delay, or a permanently altered developmental trajectory. Further studies at later developmental time points are needed.
What is the functional significance of these shifts in population entropy? Previous work suggested that the entropy of neural circuit activity may be optimally tuned at intermediate levels as a tradeoff between maximizing representational capacity at high entropy, versus maintaining error correction and regularization at low entropy (Schneidman et al., 2006). These properties can also be thought of as trading off between discrimination and generalization, respectively (Qian and Lipkin, 2011). If we assume that WT mice are optimally tuned, our findings predict that young Fmr1 KO mice should show poorer somatosensory discrimination in behavioral tasks than wildtype animals, while in contrast adult Fmr1 KO mice should perform more poorly on tasks involving generalization across somatosensory stimuli.
If the unidimensional E/I imbalance model is not sufficiently rich to capture the circuit changes observed in neurodevelopmental disorders, what should replace it? How many dimensions or degrees of freedom should a working model for a brain disorder have? Theoretical neuroscientists have long studied E/I balance in generic models of recurrent neural circuits (Brunel and Brunel, 2000; Tsodyks and Sejnowski, 1995). These models have uncovered important distinctions between ‘loose’ balanced regimes, where E and I inputs to a neuron are equal only on average, and finetuned ‘tight’ balanced regimes where E and I inputs to a neuron track each other closely on fast timescales (Denève and Machens, 2016; Hennequin et al., 2017). In principle, these generic models could be used to investigate multidimensional E/I imbalances in brain disorders (Vogels and Abbott, 2007). However, it is currently difficult to directly fit these manyparameter network models to data (although see Arakaki et al., 2017; Fisher et al., 2013; Stringer et al., 2016), and they are agnostic to circuit function. Instead we suggest an alternative, complementary approach: start by assuming a computational function for the particular neural circuit under study, then work backwards to design a model that is both sophisticated enough to capture the key information processing features of the circuit, but simple enough to interpret and link to physiological data. In this study, we considered a twoparameter model of L2/3 somatosensory cortex’s inputoutput function, which could account for both neural firing rates and correlations. Other brain circuits may demand models with more degrees of freedom. Crucially, the most informative models need not be those that include the highest level of physiological detail. All models are ultimately wrong in the sense that they make abstractions about their underlying parts, and detailed models carry the additional burden of fitting many parameters, which may be difficult to adequately constrain (O'Leary et al., 2015). Nonetheless, some models are useful (Box, 1979).
One potential use of simple parametric circuit models such as the ones we employed here may be as a tool for rationally designing candidate intervention compounds and then screening their effects on neural population activity. For example, the current study could have been extended to fit the logistic model to neural activity data from another cohort of Fmr1 KO mice that had received a candidate treatment, then ask if the fitted model parameters were closer in value to those from WT animals or Fmr1 KO controls. Approaches like this could complement the traditional strategy of designing drugs based on reversing molecular deficits and then assessing the drug’s impact on animal model behavior. Indeed, our results suggest that given the multidimensionality of circuit properties, it may prove difficult or impossible to find a single compound that can correctly reverse deficits at any age. This scenario might require a combination of drugs chosen to push circuitlevel properties towards the ‘correct’ region of parameter space. The framework we have introduced in this study can facilitate this type of highdimensional intervention analysis for diverse neurodevelopmental disorders.
Materials and methods
Mouse in vivo calcium imaging recording
Request a detailed protocolAll Ca^{2+} imaging data were published previously (Gonçalves et al., 2013). Briefly, data were collected from male and female C57Bl/6 wildtype and Fmr1 KO mice at P9–40. For each group the animal numbers were: P911, n = 13 WT and n = 9 Fmr1 KO; P1416, n = 8 WT and n = 10 Fmr1 KO; P3040, n = 7 WT and n = 6 Fmr1 KO. There were variations in the number of cells recorded from each animal. The range of cell numbers for each group were: P911, 49–198 cells in WT and 84–144 cells in Fmr1 KO; P1416, 65–119 cells in WT and 40–149 cells in Fmr1 KO; P3040, 60–114 cells in WT and 69–105 cells in Fmr1 KO. Mice were anesthetized with isoflurane, and a cranial window was fitted over primary somatosensory cortex by stereotaxic coordinates. Mice were then transferred to a twophoton microscope and headfixed to the stage while still under isoflurane anesthesia. 2–4 injections of the Ca^{2+} sensitive OregonGreen BAPTA1 (OGB) dye and sulforhodamine101 (to visualize astrocytes) were injected 200 um below the dura. Calcium imaging was performed using a TiSapphire Chameleon Ultra II laser (Coherent) tuned to 800 nm. Imaging in unanesthetized mice began within 30–60 min of stopping the flow of isoflurane after the last OGB injection. Images were acquired using ScanImage software (Pologruto et al., 2004) written in MATLAB (MathWorks; RRID:SCR_001622). Wholefield images were collected using a 20 × 0.95 NA objective (Olympus) at an acquisition speed of 3.9 Hz (512 × 128 pixels). Several 3 min movies were concatenated and brief segments of motion artifacts were removed (always <10 s total). Data were corrected for x–y drift. Cell contours were automatically detected and the average ΔF/F signal of each cell body was calculated at each time point. The baseline F for each cell was calculated as the mean ROI fluorescence across the entire 3 min timeseries. Neuronal and neuropil signals were analyzed separately and astrocytic signals were excluded from analysis. Each ΔF/F trace was lowpass filtered using a Butterworth filter (coefficient of 0.16) and deconvolved with a 2 s singleexponential kernel (Yaksi and Friedrich, 2006). To remove baseline noise, the standard deviation of all points below zero in each deconvolved trace was calculated, multiplied by two, and set as the positive threshold level below which all points in the deconvolved trace were set to zero. Estimated firing rates of the neurons, r_{i}(t), were then obtained by multiplying the deconvolved trace by a factor of 78.4, previously derived empirically from cellattached recordings in vivo (Golshani et al., 2009).
Computational methods
Request a detailed protocolData analysis and logistic model calculations were done using MATLAB (Mathworks; RRID:SCR_001622). All simulation codes are available online at https://github.com/cianodonnell/ODonnelletal_2017_imbalances (copy archived at https://github.com/elifesciencespublications/ODonnelletal_2017_imbalances), and the population tracking model code (O'Donnell et al., 2017) is available at https://github.com/cianodonnell/PopulationTracking. (copy archived at https://github.com/elifesciencespublications/PopulationTracking)
Detailed layer 2/3 model simulations
Request a detailed protocolLayer 2/3 model simulations (Figures 1 and 2) were implemented with the Pythonbased simulator Brian 2 (http://briansimulator.org; RRID:SCR_002998) (Goodman and Brette, 2009), and results analyzed with MATLAB (Mathworks; RRID:SCR_001622). The model consisted of four populations of reciprocally connected leaky integrateandfire neurons representing a L2/3 somatosensory barrel circuit: 1700 excitatory neurons, 70 PV inhibitory neurons, 115 5HT_{3A}R inhibitory neurons, and 45 SOM inhibitory neurons, driven by a separate population of 1500 excitatory spike sources representing input from L4. Cell numbers were estimated by combining layerspecific excitatory and inhibitory cell count information from (Lefort et al., 2009) with the approximate percentages of the three inhibitory cell groups given by (Petersen and Crochet, 2013). The voltage V of each neuron evolved as
where $R}_{in$ is the input resistance, $E}_{rev,e}\phantom{\rule{thinmathspace}{0ex}}and\phantom{\rule{thinmathspace}{0ex}}{E}_{rev,i$ are the excitatory and inhibitory synaptic reversal potentials respectively, $\tau}_{m$ is the membrane time constant, and $g}_{e}\phantom{\rule{thinmathspace}{0ex}}\mathrm{a}\mathrm{n}\mathrm{d}\phantom{\rule{thinmathspace}{0ex}}{g}_{i$ are the summed excitatory and inhibitory synaptic input conductances, respectively. Between input events the total excitatory synaptic condunctance $g}_{e$ evolved in time according to the equation
where τ_{syn,e} is the excitatory synaptic time constant. Similar equations governed the inhibitory conductances. When a spike arrived at a synapse, a Bernoulli random number was drawn with release probability set according to the particular synaptic connection type. If this number was equal to one, then the total synaptic conductance for that neuron was instantaneously incremented by the specific amplitude of the chosen conductance for that individual synapse, indexed j: ${g}_{e}\to {g}_{e}+{\stackrel{}{g}}_{j}$.
All synaptic connections were formed probabilistically by drawing independent random Bernoulli variables with connection typespecific probabilities. Synaptic PSP amplitudes were drawn independently for each synapse from a lognormal distribution constrained by the experimentally reported mean and median values for each particular connection type. The maximum postsynaptic potential amplitude was set to 8 mV. Synapses in the model were conductancebased, but since synaptic strengths reported in the literature were typically in terms of EPSP/IPSP amplitude, in accordance with how the experiments were performed (Avermann et al., 2012), we set each maximal synaptic conductance as the value needed to generate a PSP of the desired amplitude when the target neuron started at resting potential in the case of EPSPs or −55 mV in the case of IPSPs, which we computed analytically. Refractory periods were calculated as the interspikeinterval corresponding to the maximal experimentally reported firing rate. Release probability and synaptic strength values for unconnected neurons are excluded from Table 1. Excitatory synaptic time constants were set at 2 ms, which is typical for the fast component of AMPA receptor responses, but could not be estimated from the PSP statistics in (Avermann et al., 2012) because of masking by the slower membrane time constant. The mathematical form of our model meant that inhibitory synaptic time constants needed to be equal for all incoming inhibitory synapses to a neuron. We set these to 40 ms for E, 5HT_{3A}R and SOM neurons and 16 ms for PV neurons, which were the typical values of the IPSP decay time constants in the (Avermann et al., 2012) dataset. Due to lack of direct data for this circuit, connection probabilities for synapses from L4 E neurons to E, PV and SOM L2/3 neurons was set to a reasonable cortical value of 0.15, while 5HT_{3A}R neurons did not receive any input from L4 (Gentet et al., 2012). Similarly due to a lack of direct data, we set synaptic release probabilities for connections from L4 to L2/3 neurons to a typical cortical value of 0.25, while mean and median L4 excitatory PSP amplitudes onto L2/3 PV and SOM were set to 0.8 and 0.48 mV, respectively, to match reported data for L4 EPSP amplitudes onto L2/3 E neurons (Lefort et al., 2009). The differential equations were solved using the forward Euler method with an integration timestep of 0.01 ms. Each simulation run was 50 ms long, during which we recorded whether or not each neuron responded. In the rare cases where a neuron spiked more than once, we disregarded the extra spikes. L4 neuron dynamics were not explicitly simulated, but instead modeled only as a set of output spike trains. After selecting the subset of active L4 neurons, spike times were drawn randomly from a Gaussian distribution with standard deviation of 2 ms. We repeated the simulations 10 times for this identical input pattern to average over the noise due to probabilistic vesicle release. We repeated this procedure further 10 times for different random allocations of the ‘ON’ inputs. Then, a neuron’s ON probability was defined as the fraction of these $10\times 10=100$ simulations for which it responded with one or more spikes. Finally, we repeated the entire procedure for varying levels of L4 input sparsity.
For the simulations presented in Figure 3 we varied only 76 model parameters, which is 24 less than the total number of 100 model parameters listed in Table 1. We excluded the four neuronal refractory periods (because in almost all simulations each neuron spiked a maximum of once, making the refractory period irrelevant), and the six connection probabilities that were fixed at zero. Finally, we grouped together the mean and median PSP amplitudes for each of the fourteen nonzero synaptic connections, so that both parameters were increased or decreased by the same fraction in tandem. Together these choices reduced the number of test parameters from 100 to 76.
For all parameters that naturally range from 0 upwards, such as the number of neurons or release probability, we increased or decreased their values during testing in the most intuitive way, by adding ±20% of the baseline value. However, this method was less useful for other parameters, such as cell resting voltage, for which we reasoned it made more sense to scale relative to another parameter, such as spike threshold. As a result, we varied (1) resting voltage relative to its difference from spike threshold; (2) spike threshold relative to its difference with resting voltage; (3) excitatory synaptic reversal potentials relative to resting voltage; (4) inhibitory synaptic reversal potentials relative to spike threshold.
Source [1] is (Lefort et al., 2009), [2] is (Avermann et al., 2012), [3] is (Fanselow et al., 2008), [4] is (Kinnischtzke et al., 2012), [5] is (Fino and Yuste, 2011). N is number of neurons, Vrest is resting potential, Vth is spike voltage threshold, Rin is input resistance, tref is refractory period, τm is the membrane time constant, τsyn is the synaptic time constant with the first subscript indicating the postsynaptic neuron type and the second subscript the neurotransmitter type of the presynaptic neuron (e or i), Erev is the synaptic reversal potential, pcon is the synaptic connection probability, prel is the synaptic release probability, w is the mean or median postsynaptic potential amplitude as indicated. For all neuronal parameters, the subscript indicates the neuron type: E is L2/3 excitatory neurons, Ipv is PV neurons, I5ht is 5HT_{3A}R neurons, Isom is SOM neurons, and EL4 is L4 excitatory neurons. For synaptic parameters, the first and second subscripts indicate the pre and postsynaptic neuron types, respectively.
An important caveat is that although this model may be considered detailed by some measures, it also simplifies many aspects of L2/3 circuit. For example, we assumed that all 5HT_{3A}R cells were homogeneous, even though they likely separate into different subclasses with typespecific connectivity (Gentet, 2012; Petersen and Crochet, 2013). Layer 2 and layer 3 may also consist of distinct cell populations (Petersen and Crochet, 2013). Not all likely connections were included in the model (Dalezios et al., 2002; Pfeffer et al., 2013), and connectivity was assumed to be random, even though it is likely nonrandom (Tomm et al., 2014). Although these choices will likely not affect the conclusions of the current study, they may be important to consider for future work that seeks to understand the biological function of the L2/3 somatosensory microcircuit.
Logistic model
Request a detailed protocolFrom the L2/3 circuit model simulations, we numerically estimated the probability q that each neuron in the model fires a spike as a function of the fraction of L4 inputs that were active, f. We then used the generalized linear model regression tool ‘glmfit’ in MATLAB to find the best fit of the two logistic model parameters for each neuron:
where the parameter β represents the slope, and the parameter $f}_{1/2$ represents the fraction of active L4 neurons at which the response probability $q=0.5$. For clarity of presentation, in the main text we converted this $f}_{1/2$ parameter to what we termed the ‘threshold’, $f}_{thresh$, which we defined as the fraction of L4 neurons needed to reach a specified spike probability, $q}_{thresh$. Throughout the study, we fixed ${q}_{thresh}=0.01$. The threshold is related to q_{thresh} via the inverse of the logistic function
We computed firing rates and pairwise correlations from the logistic model (Figures 4–5) in the following way. First, we assumed that the fraction of active L4 neurons is described by a normally distributed random variable with zero mean and unit variance:
We defined the $\beta $ and ${f}_{1/2}$ parameters relative to the mean and standard deviation of the input distribution. Since $q$ is a monotonically increasing function of $f$, the probability distribution for $q$ is
where $f\left(q\right)$ is the inverse of the logistic function $q\left(f\right)$ and
$\frac{df}{dq}=\frac{{\left(\mathrm{exp}\left(\beta \left(f{f}_{1/2}\right)\right)+1\right)}^{2}}{\beta \mathrm{exp}(\beta (f{f}_{1/2}\left)\right)}$. We calculate a neuron’s mean firing rate $\mu $ as the expectation of q,
We calculate the pairwise covariance of two homogeneous neurons driven by a common input f as
$\text{cov}=\mathbb{E}\left[{q}^{2}\right]{\left(\mathbb{E}\left[q\right]\right)}^{2}=\mathbb{E}\left[{q}^{2}\right]{\mu}^{2}={\int}_{0\text{}}^{1}\left[{q}^{2}\times p\left(f\left(q\right)\right)\left\frac{df}{dq}\right\right]dq{\mu}^{2}$, then find the pairwise correlation by normalizing the covariance by the neurons’ shared variance,$\text{var}=\mu (1\mu )$.
For fitting the logistic model to the recorded neural firing rates and correlations (Figure 5), we considered a population model where the joint probability distribution across threshold and slope was specified by a 2D Gaussian, which has five parameters: threshold mean and s.d., slope mean and s.d., and slopethreshold correlation. The three constraint statistics we considered from the neural population data were the mean neural ON probability, the s.d. of neural ON probabilities, and the mean pairwise correlations. We found the bestfit model parameters for each dataset using stochastic gradient descent (code available at https://github.com/cianodonnell/ODonnelletal_2017_imbalances). Briefly, the fitting procedure followed: (1) initialize the parameters at a starting guess points, (2) compute the predicted three output firing statistics via numerical integration over the model’s probability distributions, (3) compute the fitting error as the summed squared difference between the model output predictions and the target values, (4) generate a new set of parameter values by adding a small perturbation of a zeromean Gaussian random number to each parameter, (5) compute the new output statistics, (5) recompute the fitting error, (6) if the new error is smaller than the old error, accept the updated parameter values, otherwise reject them and revert to the old parameters, (7) return to step 4 unless the error is lower than the desired tolerance. We checked for fit convergence by sampling a large number of logistic model parameters from the fitted 2D Gaussian, drawing binary samples from these logistic model neurons and computing the ON probability mean and s.d., and mean pairwise correlation from the synthetic binary samples, and comparing the computed statistical values to the original data statistics (Figure 5 – figure supplement 1). For the sensitivity analysis presented in Figure 6, we numerically computed the partial derivative in mean firing rate and pairwise correlation with respect to the mean slope and mean threshold parameters in the population logistical model, using standard finite difference methods.
Statistical tests
Request a detailed protocolTo avoid parametric assumptions, all statistical tests were done using standard bootstrapping methods with customwritten MATLAB scripts. For example, when assessing the observed difference between two group means $\mathrm{\Delta}{\mu}_{obs}$ we performed the following procedure to calculate a pvalue. First, we pool the data points from the two groups to create a null set $\mathrm{S}}_{\mathrm{n}\mathrm{u}\mathrm{l}\mathrm{l}$. We then construct two hypothetical groups of samples S_{1} and S_{2} from this by randomly drawing n_{1} and n_{2} samples with replacement from $\mathrm{S}}_{\mathrm{n}\mathrm{u}\mathrm{l}\mathrm{l}$, where n_{1} and n_{2} are the number of data points in the original groups 1 and 2 respectively. We take the mean of both hypothetical sets $\mathrm{\Delta}{\mu}_{1}$ and $\mathrm{\Delta}{\mu}_{2}$ and calculate their difference $\mathrm{\Delta}{\mu}_{null}=\mathrm{\Delta}{\mu}_{1}\mathrm{\Delta}{\mu}_{2}$. We then repeat the entire procedure 10^{7} times to build up a histogram of $\mathrm{\Delta}{\mu}_{null}$. This distribution is always centered at zero. After normalizing, this can be interpreted as the probability distribution $f(\mathrm{\Delta}{\mu}_{null})$ for observing a group mean difference of $\mathrm{\Delta}{\mu}_{null}$ purely by chance if the data were actually sampled from the same null distribution. Then the final pvalue for the probability of finding a group difference of at least $\mathrm{\Delta}{\mu}_{obs}$ in either direction is given by $p={\int}_{\infty}^{\Delta {\mu}_{obs}}f\left(\Delta {\mu}_{null}\right)d\Delta {\mu}_{null}+{\int}_{\Delta {\mu}_{obs}}^{\infty}f\left(\Delta {\mu}_{null}\right)d\Delta {\mu}_{null}$.
For Figure 5C, we estimated twodimensional 95% confidence ellipses for the shift in mean slopethreshold parameters between Fmr1 KO and WT by computing the sample error variances and covariance through bootstrapping. Then, the 95% confidence ellipse can be computed using the Chisquared distribution. We plotted the confidence interval ellipse using the MATLAB function error_ellipse.m, downloaded from https://www.mathworks.com/matlabcentral/fileexchange/4705errorellipse.
Conversion from firing rate to ON/OFF probabilities for Ca^{2+} imaging data
Request a detailed protocolFor the Ca^{2+} imaging data, we began with estimated firing rate time series r_{i}(t) for each neuron i recorded as part of a population of N neurons. For later parts of the analysis we needed to convert these firing rates to binary ON/OFF values. This conversion involves a choice. One option would be to simply threshold the data, but this would throw away information about the magnitude of the firing rate. We instead take a probabilistic approach where rather than deciding definitively whether a given neuron was ON or OFF in a given time bin, we calculate the probability that the neuron was ON or OFF by assuming that neurons fire action potentials according to an inhomogeneous Poisson process with rate r_{i}(t). The mean number of spikes λ_{i}(t) expected in a time bin of width Δt is λ_{i}(t)=r_{i}(t)Δt. We choose Δt = 1 s. Under the Poisson model the actual number of spikes m in a particular time bin is a random variable that follows the Poisson distribution P(m = k) = λ^{k} exp(λ)/k!. We considered a neuron active (ON) if it is firing one or more spikes in a given time bin. Hence, the probability that a neuron is ON is p_{on}(t) = 1 P(m=0) = 1exp(λ). This approach has two advantages over thresholding: (1) it preserves some information about the magnitude of firing rates and (2) it acts to regularize the probability distribution for the number of neurons active by essentially smoothing nearby values together.
Entropy estimation for large numbers of neurons
Request a detailed protocolEntropy was estimated by fitting a statistical model we recently developed, called the population tracking model (O'Donnell, 2017a), to the binarized Ca^{2+} imaging data. Briefly, the population tracking model fits two aspects of the data: the probability distribution for the number of neurons synchronously active in the population, and also the conditional firing probability that each individual neuron is active given the population count. Hence, the model captures both some aggregate statistics of the population activity, and some aspects of the heterogeneity across neurons. See O'Donnell, 2017b) for complete details and validation of the method. Code for fitting the model to data is available at https://github.com/cianodonnell/PopulationTracking.
The entropy/neuron generally decreased with the number of neurons considered as result of the population correlations (Figure 7B), so we needed to control for neural population size when comparing data from different experimental groups. On the one hand, we would like to study as large a number of neurons as possible, because we expect the effects of collective network dynamics to be stronger for large population sizes and this may be the regime where differences between the groups emerge. On the other hand, our recording methods allowed us to sample only typically around ~100 neurons at a time, and as few as 40 neurons in some animals. Hence, we proceeded by first estimating the entropy/neuron in each animal by calculating the entropy of random subsets of neurons of varying size from 10 to 100 (if possible) in steps of 10. For each population size we sampled a large number of independent subsets, calculated the entropy of each. Finally, for each dataset we fit a double exponential function to the estimated entropy/neuron as a function of the number of neurons: H/N = A*exp(b*N)+C*exp(d*N)+e, and used this fit to estimate H/N for 100 neurons.
Data availability
References

Microcircuits of excitatory and inhibitory neurons in layer 2/3 of mouse barrel cortexJournal of Neurophysiology 107:3116–3134.https://doi.org/10.1152/jn.00917.2011

Loss of Tsc1 in vivo impairs hippocampal mGluRLTD and increases excitatory synaptic functionJournal of Neuroscience 31:8862–8869.https://doi.org/10.1523/JNEUROSCI.161711.2011

Autism and abnormal development of brain connectivityJournal of Neuroscience 24:9228–9231.https://doi.org/10.1523/JNEUROSCI.334004.2004

BookRobustness in the strategy of scientific model buildingIn: Launer R. L, Wilkinson G. N, editors. Robustness in Statistics. Academic Press. pp. 201–236.https://doi.org/10.1016/B9780124381506.500182

Dynamics of sparsely connected networks of excitatory and inhibitory spiking neuronsJournal of Computational Neuroscience 8:183–208.https://doi.org/10.1023/A:1008925309027

Circuit and plasticity defects in the developing somatosensory cortex of FMR1 knockout miceJournal of Neuroscience 28:5178–5188.https://doi.org/10.1523/JNEUROSCI.107608.2008

Structure of a single whisker representation in layer 2 of mouse somatosensory cortexJournal of Neuroscience 35:3946–3958.https://doi.org/10.1523/JNEUROSCI.388714.2015

Delayed stabilization of dendritic spines in fragile X miceJournal of Neuroscience 30:7793–7803.https://doi.org/10.1523/JNEUROSCI.057710.2010

Efficient codes and balanced networksNature Neuroscience 19:375–382.https://doi.org/10.1038/nn.4243

Cellular mechanisms for response heterogeneity among L2/3 pyramidal cells in whisker somatosensory cortexJournal of Neurophysiology 112:233–248.https://doi.org/10.1152/jn.00848.2013

Selective, statedependent activation of somatostatinexpressing inhibitory interneurons in mouse neocortexJournal of Neurophysiology 100:2640–2652.https://doi.org/10.1152/jn.90691.2008

Spontaneous activations follow a common developmental course across primary sensory areas in mouse neocortexJournal of Neurophysiology 116:431–437.https://doi.org/10.1152/jn.00172.2016

Functional diversity of supragranular GABAergic neurons in the barrel cortexFrontiers in Neural Circuits 6:52.https://doi.org/10.3389/fncir.2012.00052

Internally mediated developmental desynchronization of neocortical network activityJournal of Neuroscience 29:10890–10899.https://doi.org/10.1523/JNEUROSCI.201209.2009

Circuit level defects in the developing neocortex of Fragile X miceNature Neuroscience 16:903–909.https://doi.org/10.1038/nn.3415

Universally sloppy parameter sensitivities in systems biology modelsPLoS Computational Biology 3:e189–1878.https://doi.org/10.1371/journal.pcbi.0030189

The developmental switch in GABA polarity is delayed in fragile X miceJournal of Neuroscience 34:446–450.https://doi.org/10.1523/JNEUROSCI.444713.2014

Inhibitory plasticity: balance, control, and codependenceAnnual Review of Neuroscience 40:557–579.https://doi.org/10.1146/annurevneuro072116031005

BookDifferential Equations, Dynamical Systems, and an Introduction to Chaos (3rd ed)Academic Press.

Altered excitatoryinhibitory balance in the NMDAhypofunction model of schizophreniaFrontiers in Molecular Neuroscience 1:6.https://doi.org/10.3389/neuro.02.006.2008

Spatial organization of neuronal population responses in layer 2/3 of rat barrel cortexJournal of Neuroscience 27:13316–13328.https://doi.org/10.1523/JNEUROSCI.221007.2007

Excitation/inhibition imbalance in animal models of autism spectrum disordersBiological Psychiatry 81:838–847.https://doi.org/10.1016/j.biopsych.2016.05.011

Sensitive timewindows for susceptibility in neurodevelopmental disordersTrends in Neurosciences 35:335–344.https://doi.org/10.1016/j.tins.2012.03.005

Computational models in the age of large datasetsCurrent Opinion in Neurobiology 32:87–94.https://doi.org/10.1016/j.conb.2015.01.006

Impaired inhibitory control of cortical synchronization in fragile X syndromeJournal of Neurophysiology 106:2264–2272.https://doi.org/10.1152/jn.00421.2011

Sloppiness in spontaneously active neuronal networksJournal of Neuroscience 35:8480–8492.https://doi.org/10.1523/JNEUROSCI.442114.2015

Channelrhodopsin2assisted circuit mapping of longrange callosal projectionsNature Neuroscience 10:663–668.https://doi.org/10.1038/nn1891

Monitoring neural activity and [Ca2+] with genetically encoded Ca2+ indicatorsJournal of Neuroscience 24:9572–9579.https://doi.org/10.1523/JNEUROSCI.285404.2004

A learningstyle theory for understanding autistic behaviorsFrontiers in Human Neuroscience 5:77.https://doi.org/10.3389/fnhum.2011.00077

Extracting information from neuronal populations: information theory and decoding approachesNature Reviews Neuroscience 10:173–185.https://doi.org/10.1038/nrn2578

Spontaneous and driven cortical activity: implications for computationCurrent Opinion in Neurobiology 19:439–444.https://doi.org/10.1016/j.conb.2009.07.005

Model of autism: increased ratio of excitation/inhibition in key neural systemsGenes, Brain and Behavior 2:255–267.https://doi.org/10.1034/j.1601183X.2003.00037.x

Fmr1 knockout mice: a model to study fragile X mental retardation. The DutchBelgian Fragile X ConsortiumCell 78:23–33.

Connectiontypespecific biases make uniform random network models consistent with cortical recordingsJournal of Neurophysiology 112:1801–1814.https://doi.org/10.1152/jn.00629.2013

Rapid state switching in balanced cortical network modelsNetwork: Computation in Neural Systems 6:111–124.https://doi.org/10.1088/0954898X_6_2_001

Gating deficits in model networks: a path to schizophrenia?Pharmacopsychiatry 40:S73–S77.https://doi.org/10.1055/s2007992130
Decision letter

Frances K SkinnerReviewing Editor; University Health Network, Canada
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Multidimensional imbalances in cortical circuit activity in FragileX Syndrome mice" for consideration by eLife. Your article has been reviewed by four peer reviewers, one of whom, Frances Skinner, is a member of our Board of Reviewing Editors and the evaluation has been overseen by Huda Zoghbi as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Bill Lytton (Reviewer #2); Mark D Humphries (Reviewer #4).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This study combines computational modeling and analysis of twophoton Ca^{2+} imaging data from in vivo FragileX model mice to test the excitatoryinhibitory (EI) imbalance theory of circuit dysfunction. The microcircuit model highlights the complex nonlinear effects of perturbations to cellular components on circuit function, and combined with their logistic response model, indicates that EI imbalance theory is not sufficient to explain dissociated effects of firing rate and correlations. Analysis of data from Fmr1 knockout mice demonstrated differences in inputoutput circuit parameters across development, which could not be captured by examining firing rate and pairwise correlations between neurons alone.
Overall, the model and analysis are able to show relationships between firing statistics and underlying circuit parameters.
The demonstration of capturing networkwide changes in dynamics by a neuron's or population's P(spike) function is innovative, and opens up a range of possible applications in both modelling and data analysis – both of which are demonstrated here. A nice bonus is the characterization of the developmental trajectory of the WT and Fmr1 mice. Collectively, the work attains its goal: to evaluate the relevance of E/I balance, and suggest an alternative.
All of the reviewers found the work to be interesting and felt positive in different ways, but all of the reviewers also felt that extensive rewriting and reorganization was required as several points were confusing and/or unclear along with clarifying various issues. There are three essential aspects that the authors need to address as well as the 5 further comments.
Essential revisions:
A) A reorganization and rewrite is required so that the motivation and rationale is clear at various junctures. Some reviewers thought it would be better to present the experiment first (as described at the end of the Introduction, "In this study we compare in vivo….") rather than the model first as it is at present. That is, starting with the problem and working forward rather than starting with some solutions and seeming to look for a problem.
For example, presenting the data first would allow it to be used as a motivator for the model, and make it more accessible to a general neuro audience. However, this then would perhaps seem like P(spike) comes out of the blue. If the order is left as is with the model first, then the motivation needs to be more clearly described and the work should perhaps be presented with less emphasis on FragileX and instead a focus on the how models can highlight the simplicity of EI imbalance theories etc. With this latter option, a title change should be considered.
We leave it to the author to decide how they prefer to reorganize and rewrite to present their work in a clear, motivating fashion through the various steps, so that the reader can appreciate and grasp what are the main exciting results. That is, there is a need to bring out the logic of certain choices, details of the analyses, and the motivation for the analyses.
B) Caveats and limitations, as well as advantages and disadvantages of their model should be brought forth in the main text for the reader to appreciate up front, and not only given in the Materials and methods section.
For example, explicitly state why and how a simpler model could/should not be considered, and why they do not expect existing model choices/limitations to not necessarily affect the outcome (e.g., "varying the strength of synapses from 5HT_{3A}R inhibitory neurons to E neurons.… has little effect" might possibly be troublesome given the model setup).
Even within the Materials and methods section, more rationale for model details should be provided. Model equations need to be provided.
For the above points, additional subsections in the Results section could be helpful.
C) Models are only reliably replicable if downloaded. The model (microcircuit, logistic response model) should be made available via ModelDB, GitHub or other repository.
The following comments from the four reviewers also need to be addressed.
1) Explicit equations for population model and circuit model should be provided. Integration scheme and time constant for easy reproduction needed.
Also, several connectivities are fixed or not examined at all (e.g., 'typical cortical value of 0.15?' – what is that based on if there is no data? etc.).
While robustness is tackled via a +/20% adjustment and found to be 'sloppy', this is essentially because a dynamical systems understanding does not exist (e.g., there would be more sensitivity near a bifurcation) as the authors would know. While it is clearly not possible to do a dynamical systems analyses (and certainly true in general for highdimensional systems), I wonder if the authors could provide some understanding of the essence of their network I/O in some way to understand the sloppiness and parameter fixing or not including certain connectivities. While this could be challenging, I think this is important to bring forth in some way since the authors go on to analyse their model as if it is 'good enough' to represent the biological situation, so it seems glossed over.
Further, this forms the basis of the linkage they are making here, and it was not clear to me. That is, the interpretation of model (and its parameters) analyses relative to the calcium imaging analyses that was being compared was obscured and/or unclear.
2) I am somewhat perplexed by the decision to lead with the model and then introduce fragile X data late in the paper. In general, I would be more interested if the paper could present some problems with understanding fragile X pathophysiology, likely by interpreting it initially by using various levels of strawman simplistic singlefactor thinking with a simple or 2stage (I/E) logisticfunction explicit model. From there, the paper could demonstrate how the more sophisticated model can at least permit assessment of features which are inaccessible in a logistic model. Granted that this final L2/3 model also does not directly predict treatments fragile X, it at least points in the direction from which future solutions will come.
3) Issues regarding the construction of the microcircuit model. The authors constrained their PV and 5HT_{3A}R model parameters from the experimental work of Avermann et al., 2012. However, this paper identified two GABAergic populations: fastspiking (shown to be PVexpressing), and nonfastspiking (NFS; i.e., everything else). Although a large percentage of these NFS cell are likely the heterogeneous 5HT_{3A}R population, they would also include SOMexpressing cells. This assumption that NFS are strictly 5HT_{3A}R has led to some perhaps misleading choices in model construction and choices for intrinsic and synaptic parameters. For example, ~40% 5HT_{3A}R cells are VIP+ in S1 (Lee et al., 2010), and a primary target of VIP+ interneurons are SOM+ (Pfeffer et al., 2013; Dalezios et al., 2002), but the model has no 5HT_{3A}R→ SOM connections. Also, it has 5HT_{3A}R→5HT_{3A}R connections, but not SOM→PV (Pfeffer et al., 2013)?
Overall, these model choices will likely not affect the main results – an overall varied sensitivity to cellular changes. However, it makes the specifics of how distinct populations are affected in Figures 1 and 2 (and the associated text) difficult to interpret from a biological standpoint. Along those lines, Table 1 describes the circuit model variables and could serve as a resource for modelers and experimentalists alike. However, for clarity, readers should be aware of the assumptions made so that they can use with caution.
4) The authors found that the direction of circuit parameters change in young to mature mice is opposite in KO vs WT, which could not be captured by examining neural activity statistics. Given that there is so much redundancy and variability at the cellular level, how does knowing about these circuit properties help us in terms of mechanisms of change?
The logistic response model also fits the threshold s.d., the slope s.d., and the slopethreshold correlation. How were these parameters affected in development/in KO?
5) It is unclear what we learn from the populationlevel P(spike) model that we could not learn from the single neuron model. Two things it would be good for the authors to address here:
i) The populationlevel model lacks motivation – I would expect the fourth paragraph of subsection “Firing rates and correlations from the logistic model” to start with something like "In order to characterise the population, we fitted… because…". Why fit a 5parameter model to the whole population: why not just fit each neuron with a 2parameter model, and average over their fits? [Just as was done for the model]. Then one would also obtain the population variation in the slope and threshold parameters, and that variation need not conform to a Gaussian model as is enforced by the 5parameter model. I suspect there is an underlying issue with having enough data perneuron. Please elaborate.
ii) Having fitted a population model, not much is done with it. It seems plausible that the variation in slope and threshold (quantified by the SD) could change with development, as the underlying parameter ranges contract or expand. So what is the variation over development, and does it differ between WT and Fmr1 mice?
https://doi.org/10.7554/eLife.26724.017Author response
Thank you for submitting your article "Multidimensional imbalances in cortical circuit activity in FragileX Syndrome mice" for consideration by eLife. Your article has been reviewed by four peer reviewers, one of whom, Frances Skinner, is a member of our Board of Reviewing Editors and the evaluation has been overseen by Huda Zoghbi as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Bill Lytton (Reviewer #2); Mark D Humphries (Reviewer #4).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This study combines computational modeling and analysis of twophoton Ca^{2+} imaging data from in vivo FragileX model mice to test the excitatoryinhibitory (EI) imbalance theory of circuit dysfunction. The microcircuit model highlights the complex nonlinear effects of perturbations to cellular components on circuit function, and combined with their logistic response model, indicates that EI imbalance theory is not sufficient to explain dissociated effects of firing rate and correlations. Analysis of data from Fmr1 knockout mice demonstrated differences in inputoutput circuit parameters across development, which could not be captured by examining firing rate and pairwise correlations between neurons alone.
Overall, the model and analysis are able to show relationships between firing statistics and underlying circuit parameters.
The demonstration of capturing networkwide changes in dynamics by a neuron's or population's P(spike) function is innovative, and opens up a range of possible applications in both modelling and data analysis – both of which are demonstrated here. A nice bonus is the characterization of the developmental trajectory of the WT and Fmr1 mice. Collectively, the work attains its goal: to evaluate the relevance of E/I balance, and suggest an alternative.
All of the reviewers found the work to be interesting and felt positive in different ways, but all of the reviewers also felt that extensive rewriting and reorganization was required as several points were confusing and/or unclear along with clarifying various issues. There are three essential aspects that the authors need to address as well as the 5 further comments.
Essential revisions:
A) A reorganization and rewrite is required so that the motivation and rationale is clear at various junctures. Some reviewers thought it would be better to present the experiment first (as described at the end of the Introduction, "In this study we compare in vivo….") rather than the model first as it is at present. That is, starting with the problem and working forward rather than starting with some solutions and seeming to look for a problem.
For example, presenting the data first would allow it to be used as a motivator for the model, and make it more accessible to a general neuro audience. However, this then would perhaps seem like P(spike) comes out of the blue. If the order is left as is with the model first, then the motivation needs to be more clearly described and the work should perhaps be presented with less emphasis on FragileX and instead a focus on the how models can highlight the simplicity of EI imbalance theories etc. With this latter option, a title change should be considered.
We leave it to the author to decide how they prefer to reorganize and rewrite to present their work in a clear, motivating fashion through the various steps, so that the reader can appreciate and grasp what are the main exciting results. That is, there is a need to bring out the logic of certain choices, details of the analyses, and the motivation for the analyses.
On reflection we agree with the reviewers – much of the original manuscript’s logic may have proven difficult for a reader to follow. To address this we have now totally reworked the Introduction, inserted a new Figure 1 that includes the motivating FragileX data, and added new text to the Results in several places to help clarify the paper’s logic. We believe that the new version is much improved.
B) Caveats and limitations, as well as advantages and disadvantages of their model should be brought forth in the main text for the reader to appreciate up front, and not only given in the Materials and methods section.
For example, explicitly state why and how a simpler model could/should not be considered, and why they do not expect existing model choices/limitations to not necessarily affect the outcome (e.g., "varying the strength of synapses from 5HT_{3A}R inhibitory neurons to E neurons.… has little effect" might possibly be troublesome given the model setup).
We have added new sentences to the Results to motivate both the biophysical model (“We chose this level of detail for the model in order to relate experimentally measureable biophysical properties of neurons to their putative role in the circuit at large.”) and logistic model (“This sigmoidalshaped response profile of simulated L2/3 neurons mimics the spiking response of mouse L2/3 pyramidal cells to extracellular L4 stimulation in vitro (Elstrott et al., 2014), while the sparse, noisy and distributed network responses were reminiscent of in vivo activity following whisker stimulation (Clancy et al., 2015; Kerr et al., 2007).”), and a section comparing the logistic model with the E/I imbalance model: “This 2D logistic model has two benefits over the 1D E/I imbalance model: first, its extra degree of freedom allows for richer and more flexible fits to data, and second, by describing an inputoutput mapping for the L2/3 circuit it can capture some aspects of the computation that the circuit performs for the animal. In contrast, the E/I imbalance model is specified purely in terms of circuit components, and so is agnostic to the circuit’s computational function.” And finally, we added a new paragraph to the Discussion on the same issue.
Even within the Materials and methods section, more rationale for model details should be provided. Model equations need to be provided.
We have now added more rationale and all the model equations to the Materials and methods section.
For the above points, additional subsections in the Results section could be helpful.
C) Models are only reliably replicable if downloaded. The model (microcircuit, logistic response model) should be made available via ModelDB, GitHub or other repository.
We agree that this is an important step. We have now put all the code online at https://github.com/cianodonnell/ODonnelletal_2017_imbalances and have indicated so at the beginning of the Methods section on computational modelling methods.
The following comments from the four reviewers also need to be addressed.
1) Explicit equations for population model and circuit model should be provided. Integration scheme and time constant for easy reproduction needed.
For the detailed circuit model we used a timestep of 0.01 ms and the forward Euler integration scheme. These details and the model equations are now mentioned in the Methods.
Also, several connectivities are fixed or not examined at all (e.g., 'typical cortical value of 0.15?' – what is that based on if there is no data? etc.).
We are unsure what the reviewer’s point is here – we varied 76 out of the 100 total circuit model parameters, including those for which there was no direct data, such as the synaptic connections with probability of 0.15 that the reviewer highlights. The logic for not varying the remaining 24 parameters was given in the original Materials and methods section: “We excluded the four neuronal refractory periods (because in almost all simulations each neuron spiked a maximum of once, making the refractory period irrelevant), and the six connection probabilities that were fixed at zero. Finally, we grouped together the mean and median PSP amplitudes for each of the fourteen nonzero synaptic connections, so that both parameters were increased or decreased by the same fraction in tandem. Together these choices reduced the number of test parameters from 100 to 76.”
While robustness is tackled via a +/20% adjustment and found to be 'sloppy', this is essentially because a dynamical systems understanding does not exist (e.g., there would be more sensitivity near a bifurcation) as the authors would know. While it is clearly not possible to do a dynamical systems analyses (and certainly true in general for highdimensional systems), I wonder if the authors could provide some understanding of the essence of their network I/O in some way to understand the sloppiness and parameter fixing or not including certain connectivities. While this could be challenging, I think this is important to bring forth in some way since the authors go on to analyse their model as if it is 'good enough' to represent the biological situation, so it seems glossed over.
This is a very good point. We have added a new paragraph to the Discussion section to address this issue:
“An important caveat to our parameter sensitivity analysis is that it was linear and local to a particular point in the highdimensional model parameter space, corresponding to WT P1722 mice. Since the circuit dynamics are nonlinear, it is likely that the particular parameter sensitivities would be different in other parts of parameter space, especially near bifurcations where qualitatively different dynamics emerge (Hirsch et al., 2013). However, as long as the redundancy property is widely preserved, as suggested by studies on computational models of other biological systems (Fisher et al., 2013; Gutenkunst et al., 2007; Machta et al., 2013; Panas et al., 2015), then our conclusions for brain disorders remain valid.”
Further, this forms the basis of the linkage they are making here, and it was not clear to me. That is, the interpretation of model (and its parameters) analyses relative to the calcium imaging analyses that was being compared was obscured and/or unclear.
This is also a good point. We consider this to be the purpose of the logistic model, which acts as a bridge between the biophysical model and the data. This bridge is necessary because it is not tractable to fit the biophysical model to our calcium imaging data. We have tried hard to clarify this relationship in the new version of the manuscript.
2) I am somewhat perplexed by the decision to lead with the model and then introduce fragile X data late in the paper. In general, I would be more interested if the paper could present some problems with understanding fragile X pathophysiology, likely by interpreting it initially by using various levels of strawman simplistic singlefactor thinking with a simple or 2stage (I/E) logisticfunction explicit model. From there, the paper could demonstrate how the more sophisticated model can at least permit assessment of features which are inaccessible in a logistic model. Granted that this final L2/3 model also does not directly predict treatments fragile X, it at least points in the direction from which future solutions will come.
This reviewer seems to have misinterpreted the point of the paper: our goal was not to “understand fragilex pathophysiology” directly, but to test whether the E/I imbalance model is flexible enough to account for the neural activity alterations in any brain disorder, taking Fmr1 KO mice as an example case. We do appreciate that this logic was not clear in the original manuscript and have taken care to emphasize our exact line of thinking in the new version.
3) Issues regarding the construction of the microcircuit model. The authors constrained their PV and 5HT_{3A}R model parameters from the experimental work of Avermann et al., 2012. However, this paper identified two GABAergic populations: fastspiking (shown to be PVexpressing), and nonfastspiking (NFS; i.e., everything else). Although a large percentage of these NFS cell are likely the heterogeneous 5HT_{3A}R population, they would also include SOMexpressing cells.
When we were designing the model we were unclear about these same issues, so contacted Carl Petersen directly by email to ask if in the Avermann et al., 2012 study, the cells labelled NFS were both SOM and 5HT_{3A}R, or just 5HT_{3A}R. He replied to say, “We think that our NFS cells are almost exclusively 5HT_{3A}Rexpressing GABAergic neurons (although this is not proven)”. We then based our model on this assumption.
This assumption that NFS are strictly 5HT_{3A}R has led to some perhaps misleading choices in model construction and choices for intrinsic and synaptic parameters. For example, ~40% 5HT_{3A}R cells are VIP+ in S1 (Lee et al., 2010), and a primary target of VIP+ interneurons are SOM+ (Pfeffer et al., 2013; Dalezios et al., 2002), but the model has no 5HT_{3A}R→ SOM connections. Also, it has 5HT_{3A}R→5HT_{3A}R connections, but not SOM→PV (Pfeffer et al., 2013)?
We agree that the model could be made more detailed, and of course makes approximations by lumping together all 5HT_{3A}R cells into one population. However, splitting the population would involve multiplying the number of model parameters. Designing such models necessarily involves some tradeoff between adding more detail versus availability of constraining data, so we settled at some point on this tradeoff. A second difficulty with designing such models is the fact that different studies from different labs, or conducted in different brain areas/animal ages, etc often result in mutually inconsistent findings.
Overall, these model choices will likely not affect the main results – an overall varied sensitivity to cellular changes. However, it makes the specifics of how distinct populations are affected in Figures 1 and 2 (and the associated text) difficult to interpret from a biological standpoint. Along those lines, Table 1 describes the circuit model variables and could serve as a resource for modelers and experimentalists alike. However, for clarity, readers should be aware of the assumptions made so that they can use with caution.
This is a fair point and we have now added a paragraph to the Methods making clear these simplifications: “An important caveat is that although this model may be considered detailed by some measures, it also simplifies many aspects of L2/3 circuit. For example, we assumed that all 5HT_{3A}R cells were homogeneous, even though they likely separate into different subclasses with typespecific connectivity (Gentet, 2012; Petersen and Crochet, 2013). Layer 2 and layer 3 may also consist of distinct cell populations (Petersen and Crochet, 2013). Not all likely connections were included in the model (Dalezios et al., 2002; Pfeffer et al., 2013), and connectivity was assumed to be random, even though it is likely nonrandom (Tomm et al., 2014). Although these choices will likely not affect the conclusions of the current study, they may be important to consider for future work that seeks to understand the biological function of the L2/3 somatosensory microcircuit.
4) The authors found that the direction of circuit parameters change in young to mature mice is opposite in KO vs WT, which could not be captured by examining neural activity statistics. Given that there is so much redundancy and variability at the cellular level, how does knowing about these circuit properties help us in terms of mechanisms of change?
We would argue that the redundancy at the cellular level in fact strengthens the case for looking at the circuit level, instead of worrying so much about the neural component changes. We hope that this point is made clearer in the new version of the manuscript.
The logistic response model also fits the threshold s.d., the slope s.d., and the slopethreshold correlation. How were these parameters affected in development/in KO?
We have now added Figure 5—figure supplement 2 that displays all five fitted parameters across genotypes and development. The most obvious changes across development were in the mean slope and threshold. This was why we focused on these in the original manuscript.
Nevertheless, we also found evidence for subtle changes in slope and threshold s.d. between groups, and interestingly also a possible increase in slopethreshold correlation in Fmr1 KO P3040 animals. We have added a clarifying note to the results text that points to this supplemental material: “For each animal, we plot the full 5D parameter fits for all animals in Figure 5—figure supplement 2. We will focus on the most prominent parameter changes which were in mean slope and mean threshold, but we also note an increase in the slope s.d. between P1416 and P3040 WT animals that was not observed in Fmr1 KO, mirroring the increased heterogeneity in firing rates in the same animals (Figure 1C).
5) It is unclear what we learn from the populationlevel P(spike) model that we could not learn from the single neuron model. Two things it would be good for the authors to address here:
i) The populationlevel model lacks motivation I would expect the fourth paragraph of subsection “Firing rates and correlations from the logistic model” to start with something like "In order to characterise the population, we fitted… because…". Why fit a 5parameter model to the whole population: why not just fit each neuron with a 2parameter model, and average over their fits? [Just as was done for the model]. Then one would also obtain the population variation in the slope and threshold parameters, and that variation need not conform to a Gaussian model as is enforced by the 5parameter model. I suspect there is an underlying issue with having enough data perneuron. Please elaborate.
The reviewer’s suggestion to fit the logistic model for each neuron independently is a good idea, but unfortunately it is a difficult statistical problem. The logistic model has two parameters per neuron (slope and threshold) so fitting it separately for each of N neurons would imply fitting 2N parameters total. However, if the goal was to match the firing rates and pairwise correlations of all the neurons, then we would be faced with N + N*(N1)/2 constraints, which is much greater than 2N, for all N> 3. Hence the problem would be vastly overconstrained and therefore likely could not be solved exactly. It might be possible to attempt to find the 2N parameters that come “closest” to matching the constraints, but that is in general a difficult problem that we are interested in pursuing but is beyond the scope of this paper. Instead by using the 5parameter model we flipped the problem to an underconstrained version where we use only 3 target constraints.
ii) Having fitted a population model, not much is done with it. It seems plausible that the variation in slope and threshold (quantified by the SD) could change with development, as the underlying parameter ranges contract or expand. So what is the variation over development, and does it differ between WT and Fmr1 mice?
As discussed above in response to a previous reviewer comment, the changes in s.d. across development/genotypes were nonzero but not as widespread as the changes in mean slope and threshold. To show these to the reader we have now added a new Figure 5—figure supplement 2.
https://doi.org/10.7554/eLife.26724.018Article and author information
Author details
Funding
FRAXA Research Foundation (Postdoctoral fellowship)
 Cian O'Donnell
Howard Hughes Medical Institute
 Cian O'Donnell
 Terrence J Sejnowski
SloanSwartz
 Cian O'Donnell
 Terrence J Sejnowski
Dana Foundation
 J Tiago Gonçalves
 Carlos PorteraCailliau
John Merck Fund (20160969)
 Carlos PorteraCailliau
Simons Foundation (295438)
 Carlos PorteraCailliau
National Institute of Neurological Disorders and Stroke (RC1NS068093)
 J Tiago Gonçalves
 Carlos PorteraCailliau
Eunice Kennedy Shriver National Institute of Child Health and Human Development (R01HD054453)
 J Tiago Gonçalves
 Carlos PorteraCailliau
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Timothy O'Leary, Hannes Saal, and Alex Williams for comments on earlier versions of the manuscript. This study was supported by funding from FRAXA Research Foundation, Howard Hughes Medical Institute, SloanSwartz Foundation, the Dana Foundation, Developmental Disabilities Translational Research Program grant #20160969 (The John Merck Fund), SFARI grant 295438 (Simons Foundation) and the NIH (NICHD R01HD054453 and NINDS RC1NS068093).
Ethics
Animal experimentation: All experiments were conducted according the US National Institutes of Health guidelines for animal research, under an animal protocol (ARC#2007035) approved by the Chancellor's Animal Research Committee and the Office for the Protection of Research Subjects at the University of California, Los Angeles.
Reviewing Editor
 Frances K Skinner, University Health Network, Canada
Version history
 Received: March 11, 2017
 Accepted: October 4, 2017
 Accepted Manuscript published: October 11, 2017 (version 1)
 Accepted Manuscript updated: October 13, 2017 (version 2)
 Version of Record published: October 31, 2017 (version 3)
Copyright
© 2017, O'Donnell 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

 7,022
 Page views

 1,071
 Downloads

 46
 Citations
Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.
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

 Chromosomes and Gene Expression
 Neuroscience
Mathys et al. conducted the first singlenucleus RNAseq (snRNAseq) study of Alzheimer’s disease (AD) (Mathys et al., 2019). With bulk RNAseq, changes in gene expression across cell types can be lost, potentially masking the differentially expressed genes (DEGs) across different cell types. Through the use of singlecell techniques, the authors benefitted from increased resolution with the potential to uncover cell typespecific DEGs in AD for the first time. However, there were limitations in both their data processing and quality control and their differential expression analysis. Here, we correct these issues and use bestpractice approaches to snRNAseq differential expression, resulting in 549 times fewer DEGs at a false discovery rate of 0.05. Thus, this study highlights the impact of quality control and differential analysis methods on the discovery of diseaseassociated genes and aims to refocus the AD research field away from spuriously identified genes.

 Neuroscience
The strength of a fear memory significantly influences whether it drives adaptive or maladaptive behavior in the future. Yet, how mild and strong fear memories differ in underlying biology is not well understood. We hypothesized that this distinction may not be exclusively the result of changes within specific brain regions, but rather the outcome of collective changes in connectivity across multiple regions within the neural network. To test this, rats were fear conditioned in protocols of varying intensities to generate mild or strong memories. Neuronal activation driven by recall was measured using cfos immunohistochemistry in 12 brain regions implicated in fear learning and memory. The interregional coordinated brain activity was computed and graphbased functional networks were generated to compare how mild and strong fear memories differ at the systems level. Our results show that mild fear recall is supported by a wellconnected brain network with smallworld properties in which the amygdala is wellpositioned to be modulated by other regions. In contrast, this connectivity is disrupted in strong fear memories and the amygdala is isolated from other regions. These findings indicate that the neural systems underlying mild and strong fear memories differ, with implications for understanding and treating disorders of fear dysregulation.