Beyond excitation/inhibition imbalance in multidimensional models of neural circuit changes in brain disorders

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 one-dimensional model is rich enough to capture the multiple neural circuit alterations underlying brain disorders. Here, we combined computational simulations with analysis of in vivo two-photon Ca2+ imaging data from somatosensory cortex of Fmr1 knock-out (KO) mice, a model of Fragile-X 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.

Introduction 26 The nervous system shows complex organization at many spatial scales: from genes and 27 molecules, to cells and synapses, to neural circuits. Ultimately, the electrical and chemical 28 signaling at all of these levels must give rise to the behavioral and cognitive processes seen at 29 the whole-organism level. When trying to understand prevalent brain disorders such as autism 30 and schizophrenia, a natural question to ask is: where is the most productive level of 31 neuroscientific investigation? Traditionally, most major disorders are diagnosed entirely at the 32 behavioral level, whereas pharmaceutical interventions are targeted at correcting alterations at 33 the molecular level. However even for the most successful drugs, we have little understanding 34 of how pharmaceutical actions at the molecular level percolate up the organizational ladder to 35 affect behavior and cognition. This classic bottom-up approach may even be further 36 confounded if phenotypic heterogeneity in disorders such as autism turn out not to reflect a 37 unique cellular pathology, but rather "a perturbation of the network properties that emerge 38 when neurons interact" (Belmonte et al., 2004). These considerations imply that a more 39 promising level of analysis might be at the level of neural circuits, since the explanatory gap 40 between circuits and behavior is smaller than the gap between molecules and behavior. This 41 circuit-level viewpoint argues for a reverse-engineering approach to tackling brain disorders: 42 rather than start at the molecular level and working up, we should instead start by asking how 43 cognitive and behavioral symptoms manifest as alterations at the circuit level, then interpret 44 these changes at the levels of cells, synapses, and molecules as appropriate. 45 One prominent circuit-level hypothesis for brain disorders has been the idea of an imbalance in 46 excitatory and inhibitory signaling. First proposed as a model for autism (Rubenstein and 47 Merzenich, 2003), the concept has since been applied to many other brain disorders, including 48 Schizophrenia, Rett syndrome, fragile-X syndrome, tuberous sclerosis, and Angelman 49 Syndrome. However, a major drawback of this model is that it only considers overall activity, 50 which is 1-dimensional. It implies that either too much excitation or too much inhibition is 51 unhealthy ( Figure 1A). Although several studies have found evidence that the E/I balance is 52 indeed upset in multiple brain disorders (Bateup et  whether it is nominally true or false, but on its explanatory and predictive powers as compared 55 with competing alternative models. In this study we argue that that even if the E/I imbalance 56 model proves correct, its unidimensionality might ultimately limit its applicability, for three 57

reasons. 58
First, by placing all disorders on the same single axis, the E/I imbalance model implicitly lumps 59 together some vastly different disorders, such as epilepsy, schizophrenia and autism ( Figure  60 1A) because they share an excess of excitation. By extension it implies that the symptoms of 61 diverse disorders could be normalized solely by either enhancing or reducing the level of, say, 62 GABAergic signaling as appropriate. Although clinical trials for such GABAergic-based 63 interventions are ongoing (Braat and Kooy, 2015), no treatment for a neurodevelopmental 64 disorder based on this principle has yet been approved. 65 A second issue with the unidimensionality of the E/I imbalance model is that it lumps together 66 all excitatory and inhibitory neural circuit components. In Figure 1B we show a schematic 67 diagram of a generic neural circuit with excitatory components colored red and inhibitory 68 components colored blue. The E/I imbalance model implies that varying any of the excitatory 69 components, such as the strength of recurrent excitatory synapses or the input resistances of 70 excitatory neurons, would have the same overall effect on circuit function. In contrast, theorists 71 have found that these equivalences often do not hold even in very simple circuit models 72 (Wilson and Cowan, 1972). 73 Third, because the standard E/I imbalance model is given in terms of circuit components, not 74 circuit function, it does not specify which aspect of a neural circuit's activity should be 75 maintained for healthy performance. For example, it leaves unclear which of neuronal firing 76 rates, synchrony, or reliability of responses might be altered if E/I balance is upset. 77 To motivate our study, we began by investigating which circuit activity properties are altered in 78 a model brain disorder. We re-analyzed published in vivo 2-photon Ca 2+ imaging data we 79 previously recorded from somatosensory cortex in Fmr1 knockout mice (Gonçalves et al.,80 2013), a well-studied animal model for fragile-X syndrome (The Dutch-Belgian Fragile X 81 Consortium, 1994). We compared the data from wild-type (WT) mice with Fmr1 KO mice, across three different developmental time points: just before (P9-11) and after (P14-16) the 83 critical period for heightened activity-dependent synaptic plasticity in L2/3 barrel cortex, and a 84 more mature timepoint (P30-40). Example ΔF/F raster plots from each group are shown in 85 Figure 1C, top left. We binned the data into 1 s timebins (originally imaged at 4 Hz), then 86

120
Neural circuits consist of many components that typically interact non-linearly to generate 121 complex circuit activity dynamics. Although many properties of cortical circuit components 122 have been found to be altered in animal models of brain disorders, it remains extremely difficult 123 to predict the net effect of varying any one particular parameter on circuit activity. The E/I 124 imbalance model seeks to simplify this problem by projecting all circuit alterations onto a one-125 dimensional axis ( Figure 1A) where the goal is to achieve a 'healthy' balance of excitation and 126 inhibition. Under this model, either too much excitation or too much inhibition leads to improper 127 circuit function. circuit-level activity. Although we focused on this particular brain circuit for tractability, our 138 general conclusions and methodology should be readily applicable to other brain circuits (Frye 139 and Maclean, 2016). very little information about the stimulus on its own, implying that information must instead be 154 encoded at the circuit level as the identities of the subset of neurons that respond. 155 C: Probability of spiking as a function of the fraction of L4 neurons activated. Each curve represents the response probability of a single neuron, averaged over multiple trials and multiple permutations of active L4 cells.

D:
Each circle plots the fitted logistic slope and threshold values for a single neuron in the simulation. Circle color indicates cell type: red is excitatory, green is PV inhibitory, purple is 5HT 3A R inhibitory, blue is SOM inhibitory. Large black circles indicate mean for each cell-type.
Inset shows an example fitted logistic response function (orange) to the noisy simulation results from a single excitatory neuron (black).
To model whisker stimulation, we simulated a volley of spikes arriving from L4 as input to the 156 population of L2/3 cells. We chose a random subset of L4 neurons as ON, then sent a single 157 spike from each of these L4 cells to their target neurons in L2/3, and recorded the responses 158 of all neurons in L2/3, some of which spiked and some of which did not. We repeated this 159 identical stimulation multiple times, in order to get an average response probability for each 160 L2/3 neuron, given the probabilistic vesicle release at synapses in the model. Then we chose a 161 different random subset of L4 neurons as ON, and repeated the entire procedure. Finally, we 162 varied the fraction of L4 cells active and plotted the probability of response for each individual 163 L2/3 neuron as a function of L4 activity level ( Figure 2B-C). 164 We found that the mean response probability of each neuron in the simulation increased from 165 zero to one monotonically with increasing L4 activity level. This sigmoidal-shaped response 166 profile of simulated L2/3 neurons mimics the spiking response of mouse L2/3 pyramidal cells to 167 extracellular L4 stimulation in vitro (Elstrott et al., 2014), while the sparse, noisy and distributed 168 network responses were reminiscent of in vivo activity following whisker stimulation (Clancy et 169 al., 2015;Kerr et al., 2007). Neurons of all four cell types in the simulation responded to the L4 170 stimulus, including the SOM interneurons which did not receive direct L4 input, but were 171 instead activated by disynaptic connections via L2/3 excitatory neurons. The detailed shape of 172 the response curve varied systematically across cell types, and was heterogeneous for 173 different neurons of a given cell type. To quantify these differences, we used logistic 174 regression to fit the response profile of each neuron with a sigmoid function ( Figure 2D, inset), 175 which has just two parameters: the slope (representing the steepness of the response curve) 176 and threshold (representing the minimal fraction of L4 neurons needed to activate the cell). 177 When we plotted the fitted slope and threshold values for each neuron against each other, we 178 found that each cell type falls into a distinct cluster in this 2-dimensional space. For example, all PV inhibitory neurons had a low slope and low threshold, whereas SOM inhibitory neurons 180 had a steep slope and moderate threshold. We then used these slope-threshold 181 measurements to summarize the circuit-level input-output function of this 'default' model of 182 L2/3 somatosensory cortex. This 2-D logistic model has two benefits over the 1-D E/I 183 imbalance model: first, its extra degree of freedom allows for richer and more flexible fits to 184 data, and second, by describing an input-output mapping for the L2/3 circuit it can capture 185 some aspects of the computation that the circuit performs for the animal. In contrast, the E/I 186 imbalance model is specified purely in terms of circuit components, and so is agnostic to the 187 circuit's computational function.  In contrast, when we increased the PSP amplitude of a different excitatory synapse, the 201 recurrent connections between L2/3 excitatory neurons, we found ( Figure 3A center) that 202 excitatory and SOM inhibitory neurons had increased slope parameters relative to default, with 203 little change in their threshold parameters. 5HT 3A R inhibitory neurons had decreased slopes 204 and thresholds, while PV neurons had little change at all. As a third example we increased the 205 probability of inhibitory synaptic connections from L2/3 PV interneurons to L2/3 excitatory 206 neurons ( Figure 3A right). In this case we found that excitatory neurons had a lower slope and 207 increased threshold, SOM inhibitory neurons had a lower slope, 5HT 3A R inhibitory neurons had 208 both an increased slope and threshold, and PV neurons showed little change, even though 209 their outgoing synapses were the parameter that was altered. 210 'inhibitory' was insufficient to predict the direction of slope-threshold change. For example, the 225 two glutamatergic projections considered in Figure 3A had distinct effects on circuit function. 226 In summary, these simulations indicate that the L2/3 somatosensory cortex circuit has 227 extremely varied sensitivities to changes in its cellular components, and that the eventual 228 circuit-level consequences cannot be predicted from knowledge of the class of the perturbed 229 neurotransmitter alone. Since the E/I imbalance model groups all excitatory and inhibitory 230 components as respective equals, it cannot account for these results. 231

Firing rates and correlations from the logistic model
In the above analysis, we investigated how low-level circuit components affect a high-level 233 circuit input-output function, as parameterized by the slope and threshold of fitted logistic 234 functions. But how is this logistic input-output function related to more common measures of 235 neural population activity, such as firing rates and pairwise correlations between neurons? To 236 investigate this, we considered the following reduced statistical model of cortical activity. We 237 assumed for simplicity that the magnitude of the total input to the L2/3 circuit can be described 238 by a Gaussian distributed random variable, with zero mean and unit standard deviation ( Figure  239 4A lower left). Then we described each L2/3 neuron's input-output as a logistic function as 240 before ( Figure 4A upper left), with threshold and slope defined relative to the Gaussian input's 241 mean and standard deviation, respectively. Given this model, we can numerically calculate the 242 probability distribution over a neuron's firing probability, which in general is skewed and non-243   However we also note an increase in the slope s.d. between P14-16 and P30-40 WT animals 280 that was not observed in Fmr1 KO, mirroring the increased heterogeneity in firing rates in the 281 same animals ( Figure 1C). In Figure 5B we plot the mean slope and mean threshold fits on top 282 of the previously calculated ( Figure 4C) 2D slope-threshold maps of firing rate and correlation. 283 We found that in young animals, P9-11, most points were scattered at high values of both 284 slope and threshold (Figure 5B left). With age, the parameter fits for both genotypes moved 285 south-west towards the low slope and low threshold region of parameter space ( Figure 5B  286 center and right). The mean location of the cloud of points at each developmental age differed 287 between WT and KO. We plot the direction of shift in group mean from WT to KO in Figure 5C. 288 In young animals, P9-11 and P14-16, the KO group had both higher slope and higher 289 threshold than WT, whereas in adult animals, P30-40, the KO group had a lower slope and 290 lower threshold than WT. These results demonstrate an opposite direction of circuit parameter 291 change in young Fragile-X mice compared to adults, which was not be uncovered by 292 measures of neural firing rates and correlations ( Figure 1C). 293 Earlier we asked how sensitive the logistic model slope and threshold parameters were to 294 alterations in the many underlying neural circuit components (Figure 3). In a similar way, we 295 can also ask how sensitive the neural firing rates and correlations are to alterations in the 296 logistic slope and threshold parameters. This is important since inspection of the 2-dimensional 297 maps in Figure 3C shows that these sensitivities will differ depending on starting location within 298 the slope-threshold space. To quantify this effect, we calculated the sensitivity of both the firing 299 rate and correlations to small changes in the slope and threshold ( Figure 6A-B, see Methods), 300 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 302 threshold always increased both firing rates and correlations, as can be predicted from Figure  303 5B. However, the magnitude of sensitivities varied across animals. We found only minor 304 differences in sensitivities between genotypes ( Figure 5 -figure supplement 2), and as a result 305 we pooled the sensitivity measurements across genotypes to test for statistical differences in 306 sensitivity with developmental age. In young animals, P9-11, changes in the logistic threshold 307 (solid bars in Figure 6) had substantial effects on both firing rates and correlation. This 308 sensitivity decreased with age (p ≤ 0.013 for firing rates, p < 0.01 for correlations from P9-11 309 to P14-16), so that in adult animals, P30-40, changes in threshold had relatively little effect 310 on neural activity statistics. A different picture emerged for the logistic slope parameter (striped 311 bars in Figure 6). There, the firing rate sensitivity increased with from P9-11 to P14-16 (p < 312   What are the functional implications of these alterations in firing rates and correlations in 316 Fragile-X mice across development? To address this, we calculated the entropy of the neural 317 population activity for the data from each animal. Entropy is a quantity from information theory, 318 measured in bits, that puts a hard upper bound on the amount of information that can be 319 represented by any coding system (Cover and Thomas, 2006). Intuitively, the entropy 320 based on these curves, 50% of the time we would expect to see the same 1000-10,000 338 patterns out of a possible total 2 50 ≈ 10 15 patterns. In contrast, in older animals P14-16 and 339 P30-40 the cumulative distributions shift rightwards so that more patterns are typically 340 observed ( Figure 7A center, right). In these cases, around 1,000,000 patterns are needed to 341 account for 50% probability mass. 342   In addition to the varying magnitudes of circuit components' effect on circuit function, we also 425 found that different components shifted the circuit input-output function in different directions, 426 as defined by our 2D logistic response model (Figures 1 and 2). Even circuit parameters that 427 are nominally of the same type, such as the strength of glutamatergic synapses between 428 excitatory (E) neurons in L4 to E neurons in L2/3 or synaptic strength between E neurons 429 within L2/3, had qualitatively different effects on the circuit response to stimulation (Figure 3). for both neural firing rates and correlations. Other brain circuits may demand models with more 498 degrees of freedom. Crucially, the most informative models need not be those that include the 499 highest level of physiological detail. All models are ultimately wrong in the sense that they 500 make abstractions about their underlying parts, and detailed models carry the additional 501 burden of fitting many parameters, which may be difficult to adequately constrain (O'Leary et 502 al., 2015). Nonetheless, some models are useful (Box, 1979). 503 One potential use of simple parametric circuit models such as the ones we employed here may 504 be as a tool for rationally designing candidate intervention compounds and then screening their 505 effects on neural population activity. For example, the current study could have been extended 506 to fit the logistic model to neural activity data from another cohort of Fmr1 KO mice that had 507 received a candidate treatment, then ask if the fitted model parameters were closer in value to 508 those from WT animals or Fmr1 KO controls. Approaches like this could complement the 509 traditional strategy of designing drugs based on reversing molecular deficits and then 510 assessing the drug's impact on animal model behavior. Indeed, our results suggest that given the multi-dimensionality of circuit properties, it may prove difficult or impossible to find a single 512 compound that can correctly reverse deficits at any age. This scenario might require a 513 combination of drugs chosen to push circuit-level properties towards the 'correct' region of 514 parameter space. The framework we have introduced in this study can facilitate this type of 515 high-dimensional intervention analysis for diverse neurodevelopmental disorders. 516

Mouse in vivo calcium imaging recording 518
All Ca 2+ imaging data were published previously (Gonçalves et al., 2013). Briefly, data were 519 where R in is the input resistance, E rev,e and E rev,i are the excitatory and inhibitory synaptic 564 reversal potentials respectively, τ m is the membrane time constant, and g e and g i are the 565 summed excitatory and inhibitory synaptic input conductances respectively. Between input 566 events the total excitatory synaptic condunctance g e evolved in time according to the equation 567 where τ syn,e is the excitatory synaptic time constant. Similar equations governed the inhibitory 568 conductances. When a spike arrived at a synapse, a Bernoulli random number was drawn with 569 release probability set according to the particular synaptic connection type. If this number was 570 equal to one, then the total synaptic conductance for that neuron was instantaneously 571 incremented by the specific amplitude of the chosen conductance for that individual synapse, 572 indexed j: + → + + : . 573 All synaptic connections were formed probabilistically by drawing independent random 575 Bernoulli variables with connection type-specific probabilities. Synaptic PSP amplitudes were 576 drawn independently for each synapse from a log-normal distribution constrained by the 577 experimentally reported mean and median values for each particular connection type. The 578 maximum post-synaptic potential amplitude was set to 8 mV. Synapses in the model were 579 conductance-based, but since synaptic strengths reported in the literature were typically in 580 terms of EPSP/IPSP amplitude, in accordance with how the experiments were performed 581 lack of direct data for this circuit, connection probabilities for synapses from L4 E neurons to E, 594 PV and SOM L2/3 neurons was set to a reasonable cortical value of 0.15, while 5HT 3A R 595 neurons did not receive any input from L4 . Similarly due to a lack of direct 596 data, we set synaptic release probabilities for connections from L4 to L2/3 neurons to a typical 597 cortical value of 0.25, while mean and median L4 excitatory PSP amplitudes onto L2/3 PV and 598 SOM were set to 0.8 and 0.48 mV, respectively, to match reported data for L4 EPSP 599 amplitudes onto L2/3 E neurons (Lefort et al., 2009). The differential equations were solved 600 using the forward Euler method with an integration timestep of 0.01 ms. Each simulation run 601 was 50 ms long, during which we recorded whether or not each neuron responded. In the rare 602 cases where a neuron spiked more than once, we disregarded the extra spikes. L4 neuron 603 dynamics were not explicitly simulated, but instead modeled only as a set of output spike trains. 604 After selecting the subset of active L4 neurons, spike times were drawn randomly from a 605 Gaussian distribution with standard deviation of 2 ms. We repeated the simulations 10 times 606 for this identical input pattern to average over the noise due to probabilistic vesicle release. We 607 repeated this procedure further 10 times for different random allocations of the 'ON' inputs. 608 Then, a neuron's ON probability was defined as the fraction of these 10×10 = 100 simulations 609 for which it responded with one or more spikes. Finally, we repeated the entire procedure for 610 varying levels of L4 input sparsity. 611 For the simulations presented in Figure 3 we varied only 76 model parameters, which is 24 612 less than the total number of 100 model parameters listed in Table 1. We excluded the four 613 neuronal refractory periods (because in almost all simulations each neuron spiked a maximum 614 of once, making the refractory period irrelevant), and the six connection probabilities that were 615 fixed at zero. Finally, we grouped together the mean and median PSP amplitudes for each of 616 the fourteen non-zero synaptic connections, so that both parameters were increased or decreased by the same fraction in tandem. Together these choices reduced the number of test 618 parameters from 100 to 76. 619 For all parameters that naturally range from 0 upwards, such as the number of neurons or 620 release probability, we increased or decreased their values during testing in the most intuitive 621 way, by adding +/-20% of the baseline value. However, this method was less useful for other 622 parameters, such as cell resting voltage, for which we reasoned it made more sense to scale 623 relative to another parameter, such as spike threshold. As a result, we varied 1) resting voltage 624 relative to its difference from spike threshold; 2) spike threshold relative to its difference with 625 resting voltage; 3) excitatory synaptic reversal potentials relative to resting voltage; 4) inhibitory 626 synaptic reversal potentials relative to spike threshold. 627  (Fino and Yuste, 2011). N is number of neurons, Vrest is 629 resting potential, Vth is spike voltage threshold, Rin is input resistance, tref is refractory period, 630 τm is the membrane time constant, τsyn is the synaptic time constant with the first subscript 631 indicating the postsynaptic neuron type and the second subscript the neurotransmitter type of 632 the presynaptic neuron (e or i), Erev is the synaptic reversal potential, pcon is the synaptic 633 connection probability, prel is the synaptic release probability, w is the mean or median post-634 synaptic potential amplitude as indicated. For all neuronal parameters, the subscript indicates 635 the neuron type: E is L2/3 excitatory neurons, Ipv is PV neurons, I5ht is 5HT 3A R neurons, Isom 636 is SOM neurons, and EL4 is L4 excitatory neurons. For synaptic parameters, the first and 637 second subscripts indicate the pre-and post-synaptic neuron types, respectively. 638 An important caveat is that although this model may be considered detailed by some 639 measures, it also simplifies many aspects of L2/3 circuit. For example, we assumed that all 640 5HT 3A R cells were homogeneous, even though they likely separate into different subclasses 641 with type-specific connectivity (Gentet, 2012;Petersen and Crochet, 2013). Layer 2 and layer 642 3 may also consist of distinct cell populations (Petersen and Crochet, 2013). Not all likely 643 connections were included in the model (Dalezios et al., 2002;Pfeffer et al., 2013), and 644 connectivity was assumed to be random, even though it is likely non-random (Tomm et al., 645 2014). Although these choices will likely not affect the conclusions of the current study, they 646 may be important to consider for future work that seeks to understand the biological function of 647 the L2/3 somatosensory microcircuit. 648

Logistic model 649
From the L2/3 circuit model simulations, we numerically estimated the probability q that each 650 neuron in the model fires a spike as a function of the fraction of L4 inputs that were active, f. 651 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: 653

@ABCD (EF(GEG H/I ))
, 654 where the parameter β represents the slope, and the parameter f 1/2 represents the fraction of 655 active L4 neurons at which the response probability = 0.5. For clarity of presentation, in the 656 main text we converted this f 1/2 parameter to what we termed the 'threshold', f thresh , which we 657 defined as the fraction of L4 neurons needed to reach a specified spike probability, q thresh . 658 Throughout the study we fixed 4L-+3L = 0.01. The threshold is related to f 1/2 via the inverse of 659 the logistic function 660 We computed firing rates and pairwise correlations from the logistic model (Figures 3-4) in 662 the following way. First, we assumed that the fraction of active L4 neurons is described by a 663 normally distributed random variable with zero mean and unit variance: 664 . 665 We defined the and @/M parameters relative to the mean and standard deviation of the input 666 distribution. Since q is a monotonically increasing function of f, the probability distribution for q 667 where ( ) is the inverse of the logistic function ( ) and 669 . We calculate a neuron's mean firing rate as the expectation of q, 670 We calculate the pairwise covariance of two homogeneous neurons driven by a common input 672 For fitting the logistic model to the recorded neural firing rates and correlations ( Figure 5), we 677 considered a population model where the joint probability distribution across threshold and 678 slope was specified by a 2D Gaussian, which has five parameters: threshold mean and s.d., 679 slope mean and s.d., and slope-threshold correlation. The three constraint statistics we 680 considered from the neural population data were the mean neural ON probability, the s.d. of 681 neural ON probabilities, and the mean pairwise correlations. We found the best-fit model 682 parameters for each dataset using stochastic gradient descent (code available at 683 https://github.com/cianodonnell/ODonnelletal_2017_imbalances). Briefly, the fitting procedure 684 followed: 1) initialize the parameters at a starting guess points, 2) compute the predicted three

Statistical tests 701
To avoid parametric assumptions, all statistical tests were done using standard bootstrapping 702 methods with custom-written MATLAB scripts. For example when assessing the observed 703 difference between two group means Δμ obs we performed the following procedure to calculate 704 a p-value. First we pool the data points from the two groups to create a null set S null . We then construct two hypothetical groups of samples S 1 and S 2 from this by randomly drawing n 1 and 706 n 2 samples with replacement from S null , where n 1 and n 2 are the number of data points in the 707 original groups 1 and 2 respectively. We take the mean of both hypothetical sets Δμ 1 and Δμ 2 708 and calculate their difference Δμ null = Δμ 1 -Δμ 2 . We then repeat the entire procedure 10 7 709 times to build up a histogram of Δμ null . This distribution is always centered at zero. After 710 normalizing, this can be interpreted as the probability distribution f(Δμ null ) for observing a group 711 mean difference of Δμ null purely by chance if the data were actually sampled from the same 712 null distribution. Then the final p-value for the probability of finding a group difference of at 713 least Δμ obs in either direction is given by = . 714 For Figure 5C we estimated 2-dimensional 95% confidence ellipses for the shift in mean slope-715 threshold parameters between Fmr1 KO and WT by computing the sample error variances and 716 covariance through bootstrapping. Then the 95% confidence ellipse can be computed using 717 the Chi-squared distribution. We plotted the confidence interval ellipse using the MATLAB 718 function error_ellipse.m, downloaded from 719 https://www.mathworks.com/matlabcentral/fileexchange/4705-error-ellipse. 720

Conversion from firing rate to ON/OFF probabilities for Ca 2+ imaging data 721
For the Ca 2+ imaging data, we began with estimated firing rate time series r i (t) for each neuron 722 i recorded as part of a population of N neurons. For later parts of the analysis we needed to 723 convert these firing rates to binary ON/OFF values. This conversion involves a choice. One 724 option would be to simply threshold the data, but this would throw away information about the 725 magnitude of the firing rate. We instead take a probabilistic approach where rather than 726 deciding definitively whether a given neuron was ON or OFF in a given time bin, we calculate 727 the probability that the neuron was ON or OFF by assuming that neurons fire action potentials 728 according to an inhomogeneous Poisson process with rate r i (t). The mean number of spikes 729 λ i (t) expected in a time bin of width Δt is λ i (t)= r i (t)Δt. We choose Δt = 1 second. Under the 730 Poisson model the actual number of spikes m in a particular time bin is a random variable that 731 follows the Poisson distribution P(m=k) = λ k exp(-λ) / k!. We considered a neuron active (ON) if 732 it is firing one or more spikes in a given time bin. Hence the probability that a neuron is ON is preserves some information about the magnitude of firing rates, and 2) it acts to regularize the 735 probability distribution for the number of neurons active by essentially smoothing nearby 736 values together. 737  Black circles are WT, red are Fmr1 KO. Each data point corresponds to a recording from a single animal. Each plot shows model prediction versus raw data target value. Plots in each row correspond to data from a different age group (P9-11, P14-16, and P30-40), and each column corresponds to one of the three target activity statistics (mean firing rate, s.d. in firing rate, and mean pairwise correlation). Blue line is identity.