Determining how synaptic coupling within and between regions is modulated during sensory processing is an important topic in neuroscience. Electrophysiological recordings provide detailed information about neural spiking but have traditionally been confined to a particular region or layer of cortex. Here we develop new theoretical methods to study interactions between and within two brain regions, based on experimental measurements of spiking activity simultaneously recorded from the two regions. By systematically comparing experimentally-obtained spiking statistics to (efficiently computed) model spike rate statistics, we identify regions in model parameter space that are consistent with the experimental data. We apply our new technique to dual micro-electrode array in vivo recordings from two distinct regions: olfactory bulb (OB) and anterior piriform cortex (PC). Our analysis predicts that: i) inhibition within the afferent region (OB) has to be weaker than the inhibition within PC, ii) excitation from PC to OB is generally stronger than excitation from OB to PC, iii) excitation from PC to OB and inhibition within PC have to both be relatively strong compared to presynaptic inputs from OB. These predictions are validated in a spiking neural network model of the OB–PC pathway that satisfies the many constraints from our experimental data. We find when the derived relationships are violated, the spiking statistics no longer satisfy the constraints from the data. In principle this modeling framework can be adapted to other systems and be used to investigate relationships between other neural attributes besides network connection strengths. Thus, this work can serve as a guide to further investigations into the relationships of various neural attributes within and across different regions during sensory processing.
As experimental tools advance, measuring whole-brain dynamics with single-neuron resolution becomes closer to reality [1–4]. However, a task that remains technically elusive is to measure the interactions within and across brain regions that govern such system-wide dynamics. Here we develop a theoretical approach to elucidate such interactions based on easily-recorded properties such as mean and (co-) variance of firing rates, when they can be measured in multiple regions and in multiple activity states. Although previous theoretical studies have addressed how spiking statistics depend on various mechanisms [5–8], these studies have typically been limited to a single region, leaving open the challenge of how inter-regional interactions impact the system dynamics, and ultimately the coding of sensory signals [9–13]. As a test case for our new theoretical tools, we studied interactions in the olfactory system. We used two micro-electrode arrays to simultaneously record from olfactory bulb (OB) and anterior piriform cortex (PC). Constrained by these experimental data, we developed computational models and theory to investigate interactions within and between OB and PC. The modeling framework includes two distinct regions: a network that receives direct sensory stimuli (here, the OB), and a second neural network (PC) that is reciprocally coupled to the afferent region. Each region contains multiple individual populations, each of which is modeled with a firing rate model [14]; thus even this minimal model involves several coupled stochastic differential equations (here, six) and has a large-dimensional parameter space. Analysis of this system would be unwieldy in general; we address this by developing a novel method to compute firing statistics that is computationally efficient, captures the results of Monte Carlo simulations, and can provide analytic insight. Thorough analysis of experimental data in both the spontaneous and stimulus-evoked states leads to a number of constraints on first- and second-order spiking statistics— many of which could not be observed using data from just one micro-electrode array. In particular, we find twelve (12) constraints that are consistent across different odorant stimuli. We use our theory and modeling to study an important subset of neural attributes (synaptic strengths) and investigate what relationships, if any, must be satisfied in order to robustly capture the many constraints observed in the data. We find that: i) inhibition within OB has to be weaker than the inhibition in PC, ii) excitation from PC to OB is generally stronger than excitation from OB to PC, iii) excitation from PC to OB and inhibition within PC have to both be relatively strong compared to inputs originating in OB (inhibition within OB and excitation from OB to PC). We validate these guiding principles in a large spiking neural network (leaky integrate-and-fire, or LIF) model, by showing that the many constraints from the experimental data are all satisfied. Finally, we demonstrate that violating these relationships in the LIF model results in spiking statistics that do not satisfy all of the data constraints. Our predictions provide insights into interactions in the olfactory system that are difficult to directly measure experimentally. Importantly, these predictions were inferred from spike rates and variability, which are relatively easy to measure. We believe that the general approach we have developed—using easy-to-measure quantities to predict hard-to-measure interactions—will be valuable in diverse future investigations of how whole-brain function emerges from interactions among its constituent components.
We first present our data from simultaneous dual micro-electrode array recordings in anesthetized rats. During each 30-second trial an odor was presented for roughly one second; recordings continued for a total of 30 seconds. This sequence was repeated for 10 trials with 2-3 minutes in between trials; the protocol was repeated for another odor. Recordings were processed to extract single-unit activity; the number of units identified was: 23 in OB and 38 in PC (first recording, two odors), 18 in OB and 35 in PC (second recording, another two odors). In total, there were four different odors presented. In this paper, we focus on the spike count statistics rather than the detailed temporal structure of the neural activity (Fig 1A and 1B). We divided each 30 s trial into two segments, representing the odor-evoked state (first 2 seconds) and the spontaneous state (remaining 28 seconds). We computed first- and second-order statistics for identified units; i. e. , firing rate νk, spike count variance, and spike count covariance (we also computed two derived statistics, Fano Factor and Pearson’s correlation coefficient, for each cell or cell pair). Spike count variances, covariances and correlations were computed using time windows Twin ranging between 5 ms and 2 s. In computing population statistics we distinguished between different odors (four total), different regions (OB vs. PC), and different activity states (spontaneous vs. evoked); otherwise, we assumed statistics were stationary over time. We then sought to identify relationships among these standard measures of spiking activity. For example, we found that mean firing rate of OB cells in the evoked state was higher than the mean firing rate in the spontaneous state, or ν O B E v > ν O B S p (although there is significant variability across the population, we focus on population-averaged statistics here). We found twelve (12) robust relationships that held across all odors. Table 1 summarizes the consistent relationships we found in our data, and Figs 1C and 1D, 2 and 3 show the data exhibiting these relationships when combining all odorant stimuli (see S1 Text for statistics plotted by distinct odors). Throughout the paper, when comparing activity states the spontaneous state is in black and the evoked state in red; when comparing regions the OB cells are in blue and PC cells in green. A common observation across different animals and sensory systems, is that firing rates increase in the evoked state (see, for example, Figure 3 in [15]). Indeed, we observed that average firing rates in both the OB and PC were higher in the evoked state than in the spontaneous state (Fig 1C and 1D). Furthermore, the firing rate in the OB was larger than the firing rate in the PC, in both spontaneous and evoked states (see mean values in Fig 1C and 1D). Stimulus-induced decorrelation appears to be a widespread phenomena in many sensory systems and in many animals [7]; stimulus-induced decorrelation was previously reported in PC cells under different experimental conditions [16]. Here, we found that in the PC, the average spike count correlation is lower in the evoked state (red) than in the spontaneous state (black), at least for time windows of 0. 5 s and above (Fig 2A). Although we show a range of time windows for completeness, we focus on the larger time windows because in our experiments the odors are held for 1 s; furthermore, our theoretical methods only address long time-averaged spiking statistics. Note that stimulus-induced decorrelation in the OB cells was not consistently observed across odors. Another common observation in cortex, is for variability to decrease at the onset of stimulus [15]: in Fig 2B we see that the Fano Factor of spike counts in PC cells decreases in the evoked state (red) compared to the spontaneous state (black); note that other experimental labs have also observed this decrease in the Fano factor of PC cells (see supplemental figure S6D in [16]). Fig 2C and 2D shows a comparison of PC and OB spike count correlation in the spontaneous state and evoked state, respectively. Spike count correlation in PC (green) is larger than correlation in OB (blue) in the spontaneous state, but in the evoked state the relationship switches, at least for time windows larger than 0. 5 sec. Fig 3 shows the four remaining constraints that are consistent for all odors and for all time windows. The Fano Factor in PC (green) is larger than in OB (blue), in the spontaneous state (Fig 3A); spike count variance in PC (green) is smaller than in OB (blue) in the evoked state (Fig 3B); spike count covariance in PC (green) is smaller than in OB (blue) in the evoked state (Fig 3C); and in OB the spike count variance in the evoked state (red) is larger than spontaneous (black, Fig 3D). Throughout the paper, we scale the spike count variance and covariance by time window for aesthetic reasons; this does not affect the relative relationships. We model two distinct regions (OB and PC) with a system of six (6) stochastic differential equations, each representing the averaged activity of a neural population [14] or representative cell (see Fig 4 for a schematic of the network). For simplicity, in this section we use the word “cell” to refer to one of these populations. Each region has two excitatory (E) and one inhibitory (I) cell to account for a variety of spiking correlations. We chose to include two E cells for two reasons: first, excitatory cells are the dominant source of projections between regions; we need at least two E cells to compute an E-to-E correlation. Moreover, in our experimental data, we are most likely recording from excitatory mitral and tufted cells (we do not distinguish between mitral vs tufted here, and therefore refer to them as M/T cells); therefore, the experimental measurements of correlations are likely to have many E-to-E correlations. The arrays likely record from I cell spiking activity as well, and the inclusion of the I cell is also important for capturing the stimulus-induced decreases in correlation and Fano factor [7,15] (also see [17] who similarly used these same cell types to analyze spiking correlations in larger spiking network models). We use j ∈ {1,2, 3} to denote three OB “cells” and j ∈ {4,5, 6} for three PC cells, with j = 1 as the inhibitory OB granule cell and j = 4 as the inhibitory PC cell. The equations are: τ d x j d t = - x j + μ j + σ j η j + ∑ k g j k F (x k) (1) where F (xk) is a transfer function mapping activity to firing rate. Thus, the firing rate is: ν j = F (x j). (2) We set the transfer function to F (X) = 1 2 (1 + tanh ( (X - 0. 5) / 0. 1) ), a commonly used sigmoidal function [14] for all cells; experimental recordings of this function demonstrate it can be sigmoidal [18–20]. All cells receive noise ηj, the increment of a Weiner process, uncorrelated in time but correlated within a region: i. e. 〈ηj (t) 〉 = 0, 〈ηj (t) ηj (t + s) 〉 = δ (s), and 〈ηj (t) ηk (t + s) 〉 = cjk δ (s). We set cjk to: c j k = { 0, if j ∈ { 1,2, 3 }; k ∈ { 4,5, 6 } 1, if j = k c O B if j ≠ k; j, k ∈ { 1,2, 3 } c P C if j ≠ k; j, k ∈ { 4,5, 6 } (3) The parameters μj and σj are constants that give the input mean and input standard deviation, respectively. Within a particular region (OB or PC), all three cells receive correlated background noisy input, but there is no correlated background input provided to both PC and OB cells. This is justified by the experimental data (see Fig S9 in S2 Text); average pairwise OB-to-PC correlations are all relatively small, and in particular, less than pairwise correlations within the OB and PC. Furthermore, anatomically there are no known common inputs to both regions that are active at the same time. We also set the background correlations to be higher in PC than in OB: i. e. , c P C > c O B. This is justified in part by our array recordings, where correlated local field potential fluctuations are larger in PC than in OB. Furthermore, one source of background correlation is global synchronous activity; Murakami et al. [21] has demonstrated that state changes (i. e. , slow or fast waves as measured by EEG) strongly affect odorant responses in piriform cortex but only minimally effect olfactory bulb cells. Finally, PC has more recurrent activity than the olfactory bulb; this could lead to more recurrent common input, if not cancelled by inhibition [22]. We constructed our model to have two distinct activity states, spontaneous and evoked. We modeled the evoked state by increasing the three parameters μ1, μ2, μ3, representing mean input to the olfactory bulb (values given in Table 2). All other parameters were the same for both states. While increasing the input to the I cells in OB in the evoked state (μ1) is not anatomically accurate because granule cells do not receive direct sensory input [23], overall this captures the net effect of stimulus input to granule cells (see section Generality of Firing Rate Model Predictions for how we apply this method to a specific olfactory system). The model we have described is less realistic than a large network of spiking models (such as Hodgkin-Huxley or leaky integrate-and-fire neurons). However, its simplicity permits fast and efficient evaluation of firing rate statistics, a necessity in exploring a large parameter space. Specifically, we calculate the statistics of the coupled network by solving a system of transcendental equations Eqs 28–45, rather than using Monte Carlo simulations. These equations were derived using an approximation based on asymptotic expansions (see Materials and methods: Approximation of Firing Statistics in the Firing Rate Model for details). This fast method allowed us to evaluate many parameter combinations, and therefore constrain the unknown coupling parameters, gjk, which would otherwise be an intractable problem. Comparisons of the firing statistics computed from our method and Monte Carlo simulations show that the mean activity and firing rates are very accurate; variance and covariance (and thus correlation) are not as accurate, for larger coupling strengths (see Fig S10 in S2 Text comparing 100 random parameter sets). Nonetheless, we will find that these reduced model results are replicated by more realistic and larger spiking network models. In principle, there can be up to 36 coupling strengths, which is intractable to explore in detail. We make the following assumptions: Within OB, there is also excitatory (M/T) input to the inhibitory (granule) cells: g12 = g13 = 0. 1—these values are small because feedforward inhibition is known to be a significant component in this circuit [28]. Within PC, we also include similar connections from E to I cells: g45 = g46 = 0. 1. Recurrent E to E connections in PC are omitted; such connections can cause problems for our reduction method, resulting in oscillatory firing rates that cannot be efficiently captured. We also make the following simplifying assumptions to limit the dimension of the parameter space of interest: The resulting network model is illustrated in Fig 4. Here we use non-standard notation for the 4 main connections of interest; instead of subscripts, we use two indicative capital letters (e. g. , gIP) so that readers can easily distinguish the connections we explore, vs. unexplored connections. Thus, we were left with four undetermined coupling strengths: gIO, gIP, gEO and gEP. We comprehensively surveyed a four-dimensional parameter space in which each coupling strength |gIO|, |gIP|, gEO, gEP was chosen between 0. 1 and 2, with a interval of 0. 1, giving us 204 = 1. 6 × 105 total models. Given each choice of 4-tuple {gIO, gIP, gEO, gEP}, we computed first- and second-order statistics of both activity xk and firing rates F (xk) using the formulas given in Eqs 28–45, and checked whether the results satisfied the constraints listed in Table 1—comparing the mean statistic across all 3 cells or all 3 possible pairs in various states and regions. We found that approximately 1. 1% of all 4-tuples satisfied the constraints; we display them in Fig 4, by projecting all constraint-satisfying 4-tuples onto a two-dimensional plane where the axes are two of the four coupling parameters. We show four out of six possible pairs (the other two show qualitatively similar patterns, see Fig S11 in S2 Text): |gIO| vs. |gIP| (Fig 4A), gEO vs. gEP (Fig 4B), |gIO| vs. gEP (Fig 4C), and |gIP| vs. gEO (Fig 4D). The results from the minimal firing rate model are: The statistics computed in Eqs 28–45 rely on the assumption that the activity distributions xk are only weakly perturbed from a normal distribution; this may be violated for larger coupling strengths. Thus, we used Monte Carlo simulations of Eq 1 to check the accuracy of this approximation; specifically we performed Monte Carlo simulations only on each 4-tuple of parameters for which the analytic approximation met our constraints. The resulting parameter sets that satisfied all 12 constraints are included as red dots in Fig 4A–4D (therefore a red dot indicates that all 12 constraints were satisfied both for the analytic approximation and for the Monte Carlo simulations). The result was a smaller set of parameters, but it is evident that the qualitative results derived from the fast analytic solver hold for the Monte Carlo simulations. Moreover, these results were robust to the choice of transfer function: in Fig. S12 of S2 Text, we show that the same constraints are obtained when using a “square root” transfer function, rather than a sigmoid. How do each of the 12 data constraints (Table 1) restrict the set of possible model parameters? Fig 5 addresses this question in two ways. In Fig 5A, we show, for each constraint, the fraction (as a percent) of all 204 parameter sets for which that constraint is satisfied, when statistics are computed via the reduction method (see Materials and methods, Approximation of Firing Statistics in the Firing Rate Model). Constraints have varying levels of restrictions, but the second order firing statistics in the evoked state appear more restrictive than the others. Together, only 1. 1% of the values in parameter space satisfy all 12 constraints. In Fig 5B, we show, for each constraint, the fraction of all 204 parameter sets for which that constraint is satisfied, in both the reduction method and in Monte Carlo simulations (recall that we took the relatively conservative approach of only testing the Monte Carlo simulations on the admissible set from the reduction method (1. 1%); this yielded only 0. 13% of parameter space. The constraint that ρ O B S p < ρ P C S p has the smallest percent by far in Fig 5B. We attribute this “mismatch” to inaccuracies in our method with stronger coupling (note that gIP and gEP are both relatively strong in the admissible set); the smaller percentages in Fig 5B compared to Fig 5A are likely due to errors in the Cov and Var calculations (see Fig S10 in S2 Text), as well as possible amplification of these errors when dividing by Var in the ρ calculation. Another way to succinctly examine the structure of the four neural attributes: gIO, gEO, gIP, gEP is to consider a matrix: A (j, :) = [ g I O (j), g E O (j), g I P (j), g E P (j) ] (4) where the jth row of A corresponds to a parameter set where all 12 constraints are satisfied. We first subtract the mean, finding that [gIO, gEO, gIP, gEP]T = [−0. 62,1. 11, −1. 38,1. 29]T, which is consistent with the results described in Fig 4. A standard singular value composition (SVD) of the mean-corrected matrix, A = U Σ V T, shows that two dimensions in the parameter space accounts for 82% of the remaining variance (as quantified by the singular values) and thus provide an approximation to the structure of the valid gIO, gEO, gIP, gEP values. The eigenvectors corresponding to the largest singular values are: [gIO, gEO, gIP, gEP]T = [−0. 05,0. 60, −0. 07,0. 79]T and [gIO, gEO, gIP, gEP]T = [0. 56,0. 05,0. 82,0. 08]T; that is, they reflect high positive correlations between the two inhibitory strengths gIP and gIO, and between the two excitatory strengths gEP and gEO. Therefore, with the minimal firing rate model we predict the connectivity strengths generally satisfy: | g I O | < g E O < g E P < | g I P |. We next asked whether the full set of data constraints were necessary; would we have seen a similar relationship between connectivity strengths, while using only a subset of the constraints outlined in Table 1? Because the admissible set is defined as the intersection over all constraints, removing any constraint would likely result in a different and (if different) larger parameter space. We considered i) keeping 8 of the 12 constraints in Table 1, neglecting the constraints on the Co-variability row, and ii) keeping only 4 of the 12 constraints in Table 1, neglecting both the Variability and Co-variability rows (i. e. , only with the firing rate). Briefly, the result is that i) 21. 5% of parameters in the analytic method satisfy the constraints; ii) 33. 4% of parameters in the analytic method satisfy the constraints; compare this to 1. 1% (and 0. 13% Monte Carlo) with all 12 constraints. The relationships of the connection strengths are different than when all 12 constraints are included: for example, it is no longer true that gEP > gEO, once the covariance constraints are omitted. In general, we should expect that if we change the wiring diagram of our simple firing rate model (Fig 4), then the same experimental constraints might result in different predictions. This could be a concern since our simple firing rate model is lacking many connections and cell types that exist in the real olfactory system [23]. However, we tested one alternative wiring diagram with different neurons receiving stimulus input, no E-to-I connections within OB, and no E-to-I connections within PC. Our predictions were robust to these changes. Second and most importantly, we tested whether our predictions held in a larger network of leaky integrate-and-fire neurons. This spiking network model also had more realistic network connectivity, more closely mimicking known anatomy of real olfactory systems. The following highlight the differences between the spiking model and the firing rate model: The parameter gEO will now refer to the strength of E-to-E connections, rather than E-to-I connections, from OB to PC. The next two sections demonstrate that our predictions hold for this LIF network model (also see S3 Text). Here we show that a general leaky integrate-and-fire (LIF) spiking neuron model of the coupled OB-PC system can satisfy all 12 data constraints. Rather than try to model the exact underlying physiological details of the olfactory bulb or anterior piriform cortex, our goal is to demonstrate that the results from the minimal firing rate model can be used as a guiding principle in a more realistic coupled spiking model with conductance-based synaptic input. The LIF model does not contain all of the attributes and cell types of the olfactory system, but is a plausible model that contains: i) more granule than M/T cells in OB (a 4-to-1 ratio, comparable to the 3-to-1 ratio used in [29]); ii) E-to-E connections from OB to PC that drive the entire network within PC; iii) E-to-I (granule cell) feedback from PC to OB; iv) lack of sensory input to granule I cells in OB. We also show that the minimal firing rate model results can be applied to a generic cortical-cortical coupled population (see S3 Text). We set the four conductance strength values to: g I O = 7 g E O = 10 g I P = 20 g E P = 15; (5) See Fig 6 or Eqs 65–67 for exact definitions of gXY; these conductance strength values are dimensionless scale factors. These values were selected to satisfy the relationships derived from the analysis of the rate model (see Fig 4). In contrast to the minimal firing rate model, here the conductance values are all necessarily positive; an inhibitory reversal potential is used to capture the hyperpolarization that occurs upon receiving synaptic input. With the conductance strengths in Eq 5, and other standard parameter values in a typical LIF model, we were able to easily satisfy all 12 constraints: see Table 3 and Fig 6 (Table 4 for LIF parameter values). While the firing rates in the LIF network (Table 3) do not quantitatively match with the firing rates from the experimental data, a few qualitative trends are apparent: (i) the ratio of mean spontaneous to evoked firing rates are similar to that observed in experimental data, for both OB and PC, (ii) the same is true of the standard deviation, (iii) the ratio of the mean OB firing rate to PC firing rate is similar to what is observed in the experimental data, in both spontaneous and evoked states. Therefore, the LIF network captures the mean firing rates reasonably well. One difference between the LIF spiking network and the minimal firing rate model is that in the evoked state, mean background input to both the OB and PC cells is increased, compared to the spontaneous state (recall that in the minimal firing rate model, only the mean input to the OB cells increased in the evoked state; this ensured that stimulus-induced changes in PC were due to network activity). When the mean input to the PC cells is the same in the spontaneous and evoked states, 10 of the 12 constraints were satisfied—the exception was the correlation of PC in the evoked state, which decreased but is still larger than the spontaneous correlation (see Fig. S13 in S2 Text). The reason is that as firing rates increase, the OB spiking is more variable and the synaptic input from OB to PC is noisier, so the input to PC activity is diffused. To capture the final two constraints, we allowed mean input drive to PC to increase in the evoked state. This has also been used in previous theoretical studies to achieve stimulus-induced decreases in spiking variability and co-variability [30]. Churchland et al. [15] used an extra source of variability in the spike generating mechanism, a doubly stochastic model, which was simply removed with stimulus onset. Thus, the mechanism we employ (increased mean input with lower input variability) is consistent with other studies that analyzed stimulus-induced changes in variability [15,30].
In computing statistics for the minimal firing rate model, we only considered equilibrium firing statistics, in which a set of stationary statistics can be solved self-consistently. More sophisticated methods might be used to address oscillatory firing statistics (see [60] where the adaptive quadratic integrate-and-fire model was successfully analyzed with a reduced method); capturing the firing statistics in these other regimes is a potentially interesting direction of research. The limitation to steady-state statistics is not unique, but is shared by other approximation methods. Some methods are known to have issues when the system bifurcates [61,62] because truncation methods can fail [63]. Several authors have proposed procedures to derive population-averaged first- and second-order spiking statistics from the dynamics of single neurons. The microscopic dynamics in question may be given by a master equation [61,62,64–66], a generalized linear model [67,68], or the theta model [69,70]. (Other authors have derived rate equations at the single-neuron level, by starting with a spike response model [74] or by taking the limit of slow synapses [73].) While we would ideally use a similar procedure to derive our rate equations, none of the approaches we note here is yet adapted to deal with our setting, a heterogeneous network of leaky integrate-and-fire neurons. Instead, we focused here on perturbing from a background state in which several populations (each population modeled by a single equation) receive correlated background input but are otherwise uncoupled. This allows us to narrow our focus to how spike count co-variability from common input is modulated by recurrent connections. We also note that other recent works have used firing rate models to explain observed patterns of correlated spiking activity in response to stimuli. Rosenbaum et al. [43] have studied the spatial structure of correlation in primate visual cortex with balanced networks [71]; Keane & Gong [72] studied wave propagation in balanced network models. Designing a spiking neural network model of two different regions that satisfies the many experimental data constraints we have outlined is a difficult problem that would often be addressed via undirected simulations. We have shown that systematic analysis of a minimal firing rate model can yield valuable insights into the relative strength of unmeasured network connections. Furthermore, these insights are transferable to a more complex, physiologically realistic spiking model of the OB–PC pathway. Indeed, incorporating the relative relationships of the four conductance strengths resulted in spiking network models that satisfied the constraints. Strongly violating the relative relationships of these conductance strengths led to multiple violations of the data constraints. Because our approach can be extended to other network features, we are hopeful that the general approach we have developed—using easy-to-measure quantities to predict hard-to-measure interactions—will be valuable in future investigations into how whole-brain function emerges from interactions among its constituent components.
All procedures were carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health and approved by University of Arkansas Institutional Animal Care and Use Committee (protocol #14049). Isoflurane and urethane anesthesia were used and urethane overdose was used for euthanasia. After the array recordings were spike sorted to identify activity from distinct cells, we further processed the data as follows: We divided each 30 s trial into two segments, representing the odor-evoked state (first 2 seconds) and the spontaneous state (remaining 28 seconds). In each state, we are interested in the random spike counts of the population in a particular window of size Twin. For a particular time window, the jth neuron has a spike count instance Nj in the time interval [t, t + Twin): N j = ∑ k ∫ t t + T w i n δ (t - t k) d t (6) The spike count correlation between cells j and k is given by: ρ T = Cov (N j, N k) Var (N j) Var (N k), (7) where the covariance of spike counts is: Cov (N j, N k) = 1 n - 1 ∑ (N j - μ (N j) ) (N k - μ (N k) ). (8) Here n is the total number of observations of Nj/k, and μ (N j): = 1 n ∑ N j is the mean spike count across Twin-windows and trials. The correlation ρT is a normalized measure of the the trial-to-trial variability (i. e. , noise correlation), satisfying ρT ∈ [−1,1]; it is also referred to as the Pearson’s correlation coefficient. For each cell pair, the covariance Cov (Nj, Nk) and variance Var (Nj) are empirically calculated by averaging across different time windows within a trial and different trials. A standard measure of variability is the Fano Factor of spike counts, which is the variance scaled by the mean: F F k = Var (N k) μ (N k). (9) In principle, any of the statistics defined here might depend on the time t as well as time window size Twin; here, we assume that Var, Cov, FF, and ρT are stationary in time, and thus separate time windows based only on whether they occur in the evoked (first 2 seconds) or spontaneous (last 28 seconds) state. Each trial of experimental data has many time windows (an exception is when Twin = 2s; in the evoked state, there is only 1 window per trial); the exact number depends on the state, the value of Twin, and whether disjoint or overlapping windows are used. In this paper we use overlapping windows by half the length of Twin to calculate the spiking statistics (e. g. , if the trial length is 2 s and Twin = 1s, then there are 3 total windows per trial: [0s, 1s], [0. 5s, 1. 5s], and [1s, 2s]). The results are qualitatively similar for disjoint windows and importantly, the relationships/constraints are the same with disjoint windows. We limit the size of Twin ≤ 2s because this is the maximum duration of the evoked state, within each trial. The average spike count μ (Nj) of the jth neuron with a particular time window Twin is related to the average firing rate νj of that neuron: ν j: = μ (N j) T w i n (10) Recall that the activity in each representative cell is modeled by: τ d x j d t = - x j + μ j + σ j η j + ∑ k g j k F (x k) (11) where F (xk) is a transfer function mapping activity to firing rate. Thus, the firing rate is: ν j = F (x j). (12) The index of each region is denoted as follows: j ∈ {1,2, 3} for the 3 OB cells, and j ∈ {4,5, 6} for the 3 PC cells, with j = 1 as the inhibitory granule OB cell and j = 4 as the inhibitory PC cell (see Fig 4). In this paper, we set σ1 = σ2 = σ3 = σOB and σ4 = σ5 = σ6 = σPC (see Table 2). In the absence of coupling (i. e. gjk = 0), any pair of activity variables, (xj, xk), are bivariate normally distributed because the equations: τ d x j d t = - x j + μ j + σ j (1 - c j k ξ j (t) + c j k ξ c (t) ) (13) τ d x k d t = - x k + μ k + σ k (1 - c j k ξ k (t) + c j k ξ c (t) ) (14) describe a multi-dimensional Ornstein-Uhlenbeck process [78]. Note that we have re-written ηj/k (t) as sums of independent white noise processes ξ (t), which is always possible for Gaussian white noise. Since x j (t) = 1 τ ∫ 0 t e - (t - u) / τ [ μ j + σ j η j (u) ] d u, we calculate marginal statistics as follows: μ (j) ≡ ⟨ x j ⟩ = μ j + 0 (15) σ 2 (j) ≡ ⟨ (x j - μ (j) ) 2 ⟩ = ⟨ σ j 2 τ 2 ∫ 0 t ∫ 0 t e - (t - u) / τ η j (u) e - (t - v) / τ η j (v) d u d v ⟩ = σ j 2 τ 2 lim t → ∞ ∫ 0 t e - 2 (t - u) / τ d u = σ j 2 2 τ A similar calculation shows that in general we have: Cov (j, k) = c j k 2 τ σ j σ k (16) Thus, (xj, xk) ~N ( (μjμk), 12τ (σj2σjσkcjkσjσkcjkσk2) ). To simplify notation, we define: ρ S N (y): = 1 2 π e - y 2 / 2, the standard normal PDF (17) ρ2D (y1, y2): =12π1−cjk2exp (−12y→T (1cjkcjk1) −1y→), bivariatestandardnormal (18) With coupling, an exact expression for a joint distribution for (x1, x2, x3, x4, x5, x6) is not explicitly known. However, we can estimate this distribution (and any derived statistics, such as means and variances) using Monte Carlo simulations. All Monte Carlo simulations of the six (6) coupled SDEs were performed using a time step of 0. 01 with a standard Euler-Maruyama method, for a time of 500 units (arbitrary, but relative to the characteristic time scale τ = 1) for each of the 3000 realizations. The activity xj was sampled at each time step after an equilibration period. Furthermore, we can approximate moments of the joint distribution under the assumption of weak coupling, as described in the next section. We will now show how to compute approximate first and second order statistics for the firing rate model with coupling; i. e. , we aim to compute the mean activity 〈xj〉, mean firing rate 〈F (xj) 〉, variance and covariances of both: 〈xj xk〉 and 〈F (xj) F (xk) 〉. For a simpler exposition, we have only included twelve synaptic connections; we have excluded self (autaptic) connections and E→E connections. An equation for each statistic can be derived by first writing Eq 11 as a low-pass filter of the right-hand-side: x j (t) = 1 τ ∫ 0 t e - (t - u) / τ [ μ j + σ j η j (u) + ∑ k g j k F (x k) ] d u (19) We then take expectations, letting t → ∞, we have: μ (j): = ⟨ x j ⟩ = μ j + ⟨ ∑ k g j k F (x k) ⟩ = μ j + ∑ k g j k ⟨ F (x k) ⟩ (20) We assume the stochastic processes are ergodic, which is generally true for these types of stochastic differential equations, so that averaging over time is equivalent to averaging over the invariant measure. We will make several assumptions for computational efficiency. First, we only account for direct connections in the formulas for the first and second order statistics, assuming the terms from the indirect connections are either small or already accounted for in the direct connections. We further make the following assumptions to simplify the calculations: ⟨ ∫ 0 t F (x k (u) ) e - (t - u) / τ d u ∫ 0 t F (x k (v) ) e - (t - v) / τ d v ⟩ ≈ τ 2 E [ F 2 (x k) ] + (τ E [ F (x k) ]) 2 (21) where E [ F n (x k) ]: = ∫ F n (σ (k) y + μ (k) ) ρ S N (y) d y (22) ⟨ ∫ 0 t σ j η j (u) e - (t - u) / τ d u ∫ 0 t F (x k (v) ) e - (t - v) / τ d v ⟩ ≈ τ 2 E [ N j F (x k) ], if j ≠ k (23) where E [ N j F (x k) ]: = σ j 2 ∫ ∫ y 1 F (σ (k) y 2 + μ (k) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (24) 〈∫0tσjηj (u) e− (t−u) /τdu∫0tF (xk (v) ) e− (t−v) /τdv〉≈τ2σk2∫yF (σ (k) y+μ (k) ) ρSN (y) dy, ifj=k (25) 〈∫0tF (xj (u) ) e− (t−u) /τdu∫0tF (xk (v) ) e− (t−v) /τdv〉≈τ2E[ F (xj) F (xk) ]+τ2E[ F (xj) ] E [F (xk) ] (26) whereE[ F (xj) F (xk) ]: =∬F (σ (j) y1+μ (j) ) F (σ (k) y2+μ (k) ) ρ2D (y1, y2) dy1dy2 (27) and Nj denotes the random variable ∫ 0 t σ j η j (u) e - (t - u) / τ d u, which is by itself normally distributed with mean 0 and variance σ j 2 τ / 2. The first assumption, Eq 21, states that time-average of F (xj (t) ) multiplied by an exponential function (low-pass filter) is equal to the expected value scaled by τ/2; the second and third, Eqs 23 and 25, address Nj and F (xk (t) ), for j ≠ k and j = k respectively (similarly for Eq 26). In all of the definitions for the expected values with ρ2D, note that the underlying correlation cjk depend on the pair of interest (j, k). Finally, we assume that the activity variables (xj, xk) are pairwise normally distributed with the subsequent statistics; this is sufficient to “close” our model and solve for the statistical quantities self-consistently. This is implicitly a weak coupling assumption because with no coupling, (xj, xk) are bivariate normal random variables. The resulting approximations for the mean activity are: μ (1) = μ 1 + ∑ k = 2,3, 5,6 g 1 k ∫ F (σ (k) y + μ (k) ) ρ S N (y) d y (28) μ (2) = μ 2 + g 21 ∫ F (σ (1) y + μ (1) ) ρ S N (y) d y (29) μ (3) = μ 3 + g 31 ∫ F (σ (1) y + μ (1) ) ρ S N (y) d y (30) μ (4) = μ 4 + ∑ k = 2,3, 5,6 g 4 k ∫ F (σ (k) y + μ (k) ) ρ S N (y) d y (31) μ (5) = μ 5 + g 54 ∫ F (σ (4) y + μ (4) ) ρ S N (y) d y (32) μ (6) = μ 6 + g 64 ∫ F (σ (4) y + μ (4) ) ρ S N (y) d y. (33) The resulting approximation to the variances of the mean activity are: τ σ 2 (1) = σ 1 2 2 + ∑ k = 2,3, 5,6 g 1 k 2 2 Var (F (σ (k) Y + μ (k) ) ) + ∑ (j, k) ∈ { (2,3); (5,6) } g 1 j g 1 k Cov (F (σ (j) Y 1 + μ (j) ), F (σ (k) Y 2 + μ (k) ) ) (34) τ σ 2 (2) = σ 2 2 2 + g 21 2 2 Var (F (σ (1) Y + μ (1) ) ) + σ 2 g 21 ∫ ∫ y 1 2 F (σ (1) y 2 + μ (1) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (35) τ σ 2 (3) = σ 3 2 2 + g 31 2 2 Var (F (σ (1) Y + μ (1) ) ) + σ 3 g 31 ∫ ∫ y 1 2 F (σ (1) y 2 + μ (1) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (36) τ σ 2 (4) = σ 4 2 2 + ∑ k = 2,3, 5,6 g 4 k 2 2 Var (F (σ (k) Y + μ (k) ) ) + ∑ (j, k) ∈ { (2,3); (5,6) } g 4 j g 4 k Cov (F (σ (j) Y 1 + μ (j) ), F (σ (k) Y 2 + μ (k) ) ) (37) τ σ 2 (5) = σ 5 2 2 + g 54 2 2 Var (F (σ (4) Y + μ (4) ) ) + σ 5 g 54 ∫ ∫ y 1 2 F (σ (4) y 2 + μ (4) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (38) τ σ 2 (6) = σ 6 2 2 + g 64 2 2 Var (F (σ (4) Y + μ (4) ) ) + σ 6 g 64 ∫ ∫ y 1 2 F (σ (4) y 2 + μ (4) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (39) In Eqs 28–39, all of the variances and covariances are computed with respect to Y ∼ N (0,1) (for Var) and (Y1, Y2) ~N ( (00), 12 (1cjkcjk1) ) (for Cov); both are easy to calculate. The value cjk depends on the pairs; for example in Eq 35, the ρ2D has cjk = cOB, the background correlation value in the olfactory bulb but in Eq 34, the Cov term is with respect to ρ2D with cjk = cPC, the background correlation value in the piriform cortex. Lastly, we state the formulas for the approximations to the covariances. Although there are 15 total covariance values, we are only concerned with 6 covariance values (3 within OB and 3 within PC); we neglect all covariances between regions. First, our experimental data set shows that these covariance (and correlation) values are small (see Fig S9 in S2 Text). Second, because there is no background correlation (i. e. , common input) between PC and OB in our model, any nonzero covariance/correlation arises strictly via direct coupling. Thus, we cannot view OB-PC covariance from coupling as a small perturbation of the background (uncoupled) state; we do not expect our model to yield qualitatively accurate predictions for these statistics. The formulas for the covariances of interest are: τ Cov (1,2) = 1 2 c O B σ 1 σ 2 + σ 1 g 21 2 ∫ y 2 F (σ (1) y + μ (1) ) ρ S N (y) d y + σ 2 g 12 2 ∫ y 2 F (σ (2) y + μ (2) ) ρ S N (y) d y + σ 2 g 13 2 ∫ y 2 F (σ (3) y + μ (3) ) ρ S N (y) d y + 1 2 ∑ (j, k) g 1 j g 2 k C (j, k) (40) τ Cov (1,3) = 1 2 c O B σ 1 σ 3 + σ 1 g 31 2 ∫ y 2 F (σ (1) y + μ (1) ) ρ S N (y) d y + σ 3 g 12 2 ∫ y 2 F (σ (2) y + μ (2) ) ρ S N (y) d y + σ 3 g 13 2 ∫ y 2 F (σ (3) y + μ (3) ) ρ S N (y) d y + 1 2 ∑ (j, k) g 1 j g 3 k C (j, k) (41) τ Cov (2,3) = 1 2 c O B σ 2 σ 3 + g 21 g 31 2 Var (F (σ (1) Y + μ (1) ) ) + σ 3 g 21 + σ 2 g 31 2 ∫ ∫ y 1 2 F (σ (1) y 2 + μ (1) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (42) τ Cov (4,5) = 1 2 c P C σ 4 σ 5 + σ 4 g 54 2 ∫ y 2 F (σ (4) y + μ (4) ) ρ S N (y) d y + σ 5 g 45 2 ∫ y 2 F (σ (5) y + μ (5) ) ρ S N (y) d y + σ 5 g 46 2 ∫ y 2 F (σ (6) y + μ (6) ) ρ S N (y) d y + 1 2 ∑ (j, k) g 4 j g 5 k C (j, k) (43) τ Cov (4,6) = 1 2 c P C σ 4 σ 6 + σ 4 g 64 2 ∫ y 2 F (σ (4) y + μ (4) ) ρ S N (y) d y + σ 6 g 45 2 ∫ y 2 F (σ (5) y + μ (5) ) ρ S N (y) d y + σ 6 g 46 2 ∫ y 2 F (σ (6) y + μ (6) ) ρ S N (y) d y + 1 2 ∑ (j, k) g 4 j g 6 k C (j, k) (44) τ Cov (5,6) = 1 2 c P C σ 5 σ 6 + g 54 g 64 2 Var (F (σ (4) Y + μ (4) ) ) + σ 6 g 54 + σ 5 g 64 2 ∫ ∫ y 1 2 F (σ (4) y 2 + μ (4) ) ρ 2 D (y 1, y 2) d y 1 d y 2 (45) where C (j, k) = ∫ ∫ F (σ (j) y 1 + μ (j) ) F (σ (k) y 2 + μ (k) ) ρ 2 D (y 1, y 2) d y 1 d y 2 - (∫ F (σ (j) y + μ (j) ) ρ S N (y) d y) (∫ F (σ (k) y + μ (k) ) ρ S N (y) d y) (46) We use a generic spiking neural network model of leaky integrate-and-fire neurons to test the results of the theory. There were NOB = 100 total OB cells, of which we set 80% (80) to be granule (I-) cells and 20% (20) to be mitral/tufted (M/T) E-cells. There are known to be many more granule cells than M/T cells in the OB; this ratio of 4-to-1 is similar to other models of OB (see [29] who used 3-to-1). The equations for the OB cells are, indexed by k ∈ {1,2, …, NOB}: τ m d v k d t = μ O B - v k - g k, X I (t) (v k - E I) - g k, X E (t) (v k - E E) - g k, X P C (t - τ Δ, P C) (v k - E E) + σ O B (1 - c ˜ O B η k (t) + c ˜ O B ξ o (t) ) v k (t *) ≥ θ k ⇒ v k (t * + τ r e f) = 0 g k, X E (t) = γ X E p X E (0. 2 N O B) ∑ k ′ ∈ { presyn OB E-cells } G k ′ (t) g k, X I (t) = γ X I p X I (0. 8 N O B) ∑ k ′ ∈ { presyn OB I-cells } G k ′ (t) g k, X P C (t) = γ X, P C p X, P C (0. 8 N P C) ∑ j ′ ∈ { presyn PC E-cells } G j ′ (t) τ d, X d G k d t = - G k + A k τ r, X d A k d t = - A k + τ r, X α X ∑ l δ (t - t k, l). (65) The conductance values in the first equation gk, XI, gk, XE, and gk, XPC depend on the type of neuron vk (X ∈ {E, I}). The last conductance, g X, P C (t - τ Δ, P C) (v k - E E), models the excitatory presynaptic input (feedback) from the PC cells with a time delay of τΔ, PC. The conductance variables gk, XY (t) are dimensionless because this model was derived from scaling the original (raw) conductance variables by the leak conductance with the same dimension. The leak, inhibitory and excitatory reversal potentials are 0, E I, and E E, respectively with E I < 0 < E E (the voltage is scaled to be dimensionless, see Table 4). ξk (t) are uncorrelated white noise processes and ξo (t) is the common noise term to all NOB cells. The second equation describes the refractory period at spike time t*: when the neuron’s voltage crosses threshold θj (see below for distribution of thresholds), the neuron goes into a refractory period for τref, after which we set the neuron’s voltage to 0. The parameter γXY gives the relative weight of a connection from neuron type Y to neuron type X; the parameter pXY is the probability that any such connection exists (X, Y ∈ {E, I}). Gk is the synaptic variable associated with each cell, and dependent only on that cell’s spike times; its dynamics are given by the final two equations in Eq 65 and depend on whether k ∈ {E, I}. Finally, two of the parameters above can be equated with coupling parameters in the reduced model: g E P = γ E, P C; g I O = γ E I (66) which are dimensionless scale factors for the synaptic conductances. The PC cells had similar functional form but with different parameters (see Table 4 for parameter values). We modeled NPC = 100 total PC cells, of which 80% were excitatory and 20% inhibitory. The equations, indexed by j ∈ {1,2, …, NPC} are: τ m d v j d t = μ P C - v j - g j, X I (t) (v j - E I) - g j, X E (t) (v j - E E) - g j, X O B (t - τ Δ, O B) (v j - E E) + σ P C (1 - c ˜ P C η j (t) + c ˜ P C ξ p (t) ) v j (t *) ≥ θ j ⇒ v j (t * + τ r e f) = 0 g j, X E (t) = γ X E p X E (0. 8 N P C) ∑ j ′ ∈ { presyn PC E-cells } G j ′ (t) g j, X I (t) = γ X I p X I (0. 2 N P C) ∑ j ′ ∈ { presyn PC I-cells } G j ′ (t) g j, X O B (t) = γ X, O B p X, O B (0. 2 N O B) ∑ k ′ ∈ { presyn OB E-cells } G k ′ (t) τ d, X d G j d t = - G j + A j τ r, X d A j d t = - A j + τ r, X α X ∑ l δ (t - t j, l). (67) Excitatory synaptic input from the OB cells along the lateral olfactory tract is modeled by: g X, O B (t - τ Δ, O B) (v j - E E). The common noise term for the PC cells ξp (t) is independent of the common noise term for the OB cells ξo (t). Two of the parameters above can be equated with coupling parameters in the reduced model: g E O = γ E, O B; g I P = γ E I (68) The values of the parameters that were not stated in Table 4 were varied in the paper: g I O, g E O, g I P, g E P. To model two activity states, we allowed mean inputs to vary (see Table 4). In contrast to the reduced model, we increased both inputs to PC cells (from μPC = 0 in the spontaneous state to μPC = 0. 4 in the evoked state) as well as to OB cells; μOB = 0. 6 in the spontaneous state to μOB = 0. 9 in the evoked state only for M/T cells (OB granule cell input is the same for spontaneous and evoked). Finally, we model heterogeneity by setting the threshold values θj in the following way. Both OB and PC cells had the following distributions for θj: θ j ∼ e N (69) where N is normal distribution with mean - σ θ 2 / 2 and standard deviation σθ, so that {θj} has a log-normal distribution with mean 1 and variance: e σ θ 2 - 1. We set σθ = 0. 1, which results in firing rates ranges seen in the experimental data. Since the number of cells are modest with regards to sampling (NOB = 100, NPC = 100), we evenly sampled the log-normal distribution from the 5th to 95th percentiles (inclusive). We remark that the synaptic delays of τΔ, PC and τΔ, OB were set to modest values to capture the appreciable distances between OB and PC. This is a reasonable choice based on evidence that stimulation in PC elicit a response in OB 5-10 ms later [79]. In all Monte Carlo simulations of the coupled LIF network, we used a time step of 0. 1 ms, with 2 s of biology time for each of the 50,000 realizations (i. e. , over 27. 7 hours of biology time), enough simulated statistics to effectively have convergence.
Sensory processing is known to span multiple regions of the nervous system. However, electrophysiological recordings during sensory processing have traditionally been limited to a single region or brain layer. With recent advances in experimental techniques, recorded spiking activity from multiple regions simultaneously is feasible. However, other important quantities— such as inter-region connection strengths—cannot yet be measured. Here, we develop new theoretical tools to leverage data obtained by recording from two different brain regions simultaneously. We address the following questions: what are the crucial neural network attributes that enable sensory processing across different regions, and how are these attributes related to one another? With a novel theoretical framework to efficiently calculate spiking statistics, we can characterize a high dimensional parameter space that satisfies data constraints. We apply our results to the olfactory system to make specific predictions about effective network connectivity. Our framework relies on incorporating relatively easy-to-measure quantities to predict hard-to-measure interactions across multiple brain regions. Because this work is adaptable to other systems, we anticipate it will be a valuable tool for analysis of other larger scale brain recordings.
action potentials
medicine and health sciences
neural networks
membrane potential
brain
electrophysiology
random variables
neuroscience
covariance
mathematics
statistics (mathematics)
network analysis
research and analysis methods
computer and information sciences
animal cells
mathematical and statistical techniques
monte carlo method
probability theory
computer networks
cellular neuroscience
olfactory bulb
cell biology
neurophysiology
anatomy
physiology
neurons
biology and life sciences
cellular types
physical sciences
statistical methods
A theoretical framework for analyzing coupled neuronal networks: Application to the olfactory system