A Theoretical Framework for Analyzing Coupled Neuronal Networks: Application to the Olfactory System

TL;DR

This paper introduces a new theoretical framework for analyzing interactions between coupled neuronal networks using simultaneous spiking data, validated on olfactory system recordings, and applicable to other neural systems.

Contribution

The study develops a novel method to infer and validate neural interaction parameters from experimental spiking data, bridging electrophysiological recordings and computational modeling.

Findings

Inhibition within OB is weaker than in PC.

Excitation from PC to OB exceeds that from OB to PC.

Strong excitation from PC to OB and inhibition within PC are necessary for data consistency.

Abstract

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.

Tables

Table 1. Table 1: The 12 relationships (constraints) that hold in the experimental data across all odors.

	Spont.	Evoked	Spon. to Evoked
Firing Rate
Variability	$F{F}_{PC}>F{F}_{OB}$		$F{F}_{PC}^{Sp}>F{F}_{PC}^{Ev}$
Co-variability	${\rho }_{PC}>{\rho }_{OB}$		${\rho }_{PC}^{Sp}>{\rho }_{PC}^{Ev}$

Table 2. Table 2: Population firing rate statistics from an LIF model of the OB–PC pathway.

	Mean Firing Rate (Hz)	Std. Dev. (Hz)
${\nu }_{OB}^{Sp}$	5.5	4.6
${\nu }_{OB}^{Ev}$	6.2	4.8
${\nu }_{PC}^{Sp}$	2.1	2.6
${\nu }_{PC}^{Ev}$	4.1	5.8

Table 3. Table 3: Parameters of the rate model (Eq  1 ). The only difference between the spontaneous and evoked states, is that the mean input to OB increased in the evoked state. We set τ = 1 𝜏 1 \tau=1 throughout.

	Parameter	Definition	Spontaneous Value	Evoked Value
Olfactory Bulb	${\mu }_1$	Mean Input	13/60	26/60
	${\mu }_2$		9/60	18/60
	${\mu }_3$		7/60	14/60
	${\sigma }_{OB}$	Background Noise Level	1.4	1.4
	$c_{OB}$	OB Background Correlation	0.3	0.3
Piriform Cortex	${\mu }_4$	Mean Input	9/60	9/60
	${\mu }_5$		5/60	5/60
	${\mu }_6$		3/60	3/60
	${\sigma }_{PC}$	Background Noise Level	2	2
	$c_{PC}$	PC Background Correlation	0.35	0.35

Table 4. Table 4: Fixed parameters for the LIF OB–PC model, see Eqs  Leaky Integrate-and-Fire Model of the generic Cortical–Cortical Circuit – Leaky Integrate-and-Fire Model of the generic Cortical–Cortical Circuit .

Same for both OB and PC
Parameter	${\tau }_m$	${\tau }_{ref}$	${ℰ}_I$	${ℰ}_E$	${\tau }_{d,I}$	${\tau }_{r,I}$	${\tau }_{d,E}$	${\tau }_{r,E}$	${\alpha }_I$	${\alpha }_E$
	20 ms	2 ms	-2.5	6.5	10 ms	2 ms	5 ms	1 ms	2 Hz	1 Hz
Parameter	$N$	Spont. $\mu$	Evoked $\mu$	$\sigma$	$\tilde{c}$	${\gamma }_{EE}$	${\gamma }_{IE}$	${\gamma }_{II}$	${\tau }_{\Delta ,PC/OB}$
OB	100	0.6	0.9∗	0.05	0.5	2	4	2	10 ms
PC	100	0	0.4	0.1	0.8	5	8	6	5 ms

Table 5. Table S1: Average population firing rate (Hz) by odor and activity state.

Mean across population	All Odors	Odor 1	Odor 2	Odor 3	Odor 4
${\nu }_{OB}^{Sp}$	1.97	1.44	1.8	2.62	2.24
${\nu }_{OB}^{Ev}$	4.66	4.91	3.28	5.41	5.34
${\nu }_{PC}^{Sp}$	0.75	0.56	0.91	0.74	0.79
${\nu }_{PC}^{Ev}$	1.45	1.6	1.26	1.7	1.23

Table 6. Table S2: Standard deviation of firing rate across the population (Hz) by odor and activity state.

Std. Dev. across population	All Odors	Odor 1	Odor 2	Odor 3	Odor 4
${\nu }_{OB}^{Sp}$	3.28	2.34	3.07	4.32	3.58
${\nu }_{OB}^{Ev}$	7.14	7.55	5.55	8	8.04
${\nu }_{PC}^{Sp}$	0.93	0.83	1.08	0.96	0.95
${\nu }_{PC}^{Ev}$	1.58	2.09	1.45	1.93	1.18

Table 7. Table S3: Fixed parameters for the LIF Cortical–Cortical model.

Same for both C1 and C2
Parameter	${\tau }_m$	${\tau }_{ref}$	${ℰ}_I$	${ℰ}_E$	${\tau }_{d,I}$	${\tau }_{r,I}$	${\tau }_{d,E}$	${\tau }_{r,E}$	${\alpha }_I$	${\alpha }_E$
	20 ms	2 ms	-2.5	6.5	10 ms	2 ms	5 ms	1 ms	2 Hz	1 Hz
Parameter	$N$	Spont. $\mu$	Evoked $\mu$	$\sigma$	$\tilde{c}$	${\gamma }_{EE}$	${\gamma }_{IE}$	${\gamma }_{II}$	${\gamma }_{E,C2/C1}$	${\tau }_{\Delta ,C2/C1}$
C1	60	0.6	0.9	0.05	0.5	2	4	6	1	10 ms
C2	100	0	0.4	0.1	0.8	2	4	6	1	5 ms

Equations

$$\tau \frac{d{x}_j}{dt}$$

$${\nu }_j=F(x_j).$$

$$\displaystyle c_{jk}=\left\{\begin{array}[]{ll}0,&\hbox{if }j\in\{1,2,3\};k\in\{4,5,6\}\\ 1,&\hbox{if }j=k\\ c_{OB}&\hbox{if }j\neq k;j,k\in\{1,2,3\}\\ c_{PC}&\hbox{if }j\neq k;j,k\in\{4,5,6\}\\ \end{array}\right.$$

$$c_{PC}>c_{OB}.$$

$$A(j,:)=[gIO(j),gEO(j),gIP(j),gEP(j)]$$

$$A=U\Sigma V^T,$$

$$|gIO|<gEO<gEP<|gIP|.$$

$$gIO$$

$$gEO$$

$$gIP$$

$$gEP$$

$$N_j=\sum _k{\int }_t^{t+T_{win}}\delta (t-t_k)dt$$

$${\rho }_T=\frac{\text{Cov}(N_j,N_k)}{\sqrt{\text{Var}(N_j)\text{Var}(N_k)}},$$

$$\text{Cov}(N_j,N_k)=\frac{1}{n-1}\sum \left(N_j-\mu (N_j)\right)\left(N_k-\mu (N_k)\right).$$

$$F{F}_k=\frac{\text{Var}(N_k)}{\mu (N_k)}.$$

$${\nu }_j:=\frac{\mu (N_j)}{T_{win}}$$

$$\tau \frac{d{x}_j}{dt}=-x_j+{\mu }_j+{\sigma }_j{\eta }_j+\sum _k{g}_{jk}F(x_k)$$

$${\nu }_j=F(x_j).$$

$$\tau \frac{d{x}_j}{dt}$$

$$\tau \frac{d{x}_k}{dt}$$

$$\mu (j)\equiv ⟨x_j⟩={\mu }_j+0$$

$${\sigma }^2(j)$$

$$\text{Cov}(j,k)=\frac{c_{jk}}{2\tau }{\sigma }_j{\sigma }_k$$

$${\rho }_{SN}(y)$$

$${\rho }_{2D}(y_1,y_2)$$

$$x_j(t)$$

$$\mu (j):=⟨x_j⟩$$

$$⟨{\int }_0^{t}F(x_k(u))e^{-(t-u)/\tau }du{\int }_0^{t}F(x_k(v))e^{-(t-v)/\tau }dv⟩\approx \frac{\tau }{2}\mathbb{E}\left[F^2(x_k)\right]$$

$$where\mathbb{E}\left[F^2(x_k)\right]:=\int F^2(\sigma (k)y+\mu (k)){\rho }_{SN}(y)dy$$

$$⟨{\int }_0^t{\sigma }_j{\eta }_j(u)e^{-(t-u)/\tau }du{\int }_0^{t}F(x_k(v))e^{-(t-v)/\tau }dv⟩\approx \frac{\tau }{2}\mathbb{E}\left[N_{j}F(x_k)\right],ifj\ne k$$

$$where\mathbb{E}\left[N_{j}F(x_k)\right]:=\frac{{\sigma }_j}{\sqrt{2}}\iint y_{1}F(\sigma (k)y_2+\mu (k)){\rho }_{2D}(y_1,y_2)d{y}_{1}d{y}_2$$

$$⟨{\int }_0^t{\sigma }_j{\eta }_j(u)e^{-(t-u)/\tau }du{\int }_0^{t}F(x_k(v))e^{-(t-v)/\tau }dv⟩$$

$$\approx \frac{\tau }{2}\frac{{\sigma }_k}{\sqrt{2}}\int yF(\sigma (k)y+\mu (k)){\rho }_{SN}(y)dy,ifj=k$$

$$⟨{\int }_0^{t}F(x_j(u))e^{-(t-u)/\tau }du{\int }_0^{t}F(x_k(v))e^{-(t-v)/\tau }dv⟩\approx \frac{\tau }{2}\mathbb{E}\left[F(x_j)F(x_k)\right]$$

$$where\mathbb{E}\left[F(x_j)F(x_k)\right]:=$$

$$\iint F(\sigma (j)y_1+\mu (j))F(\sigma (k)y_2+\mu (k)){\rho }_{2D}(y_1,y_2)d{y}_{1}d{y}_2$$

$$\mu (1)={\mu }_1+\sum _{k=2,3,5,6}{g}_{1k}\int F(\sigma (k)y+\mu (k)){\rho }_{SN}(y)dy$$

Full text

A Theoretical Framework for Analyzing Coupled Neuronal Networks: Application to the Olfactory System

Andrea K. Barreiro 1, Shree Hari Gautam 2, Woodrow L. Shew 2, Cheng Ly 3*,

1 Department of Mathematics, Southern Methodist University, Dallas, TX 75275 U.S.A.

2 Department of Physics, University of Arkansas, Fayetteville, AR 72701 U.S.A.

3 Department of Statistical Sciences and Operations Research, Virginia Commonwealth University, Richmond, VA 23284 U.S.A.

• CLy@vcu.edu

Abstract

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.

Author Summary

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.

Introduction

As experimental tools advance, measuring whole-brain dynamics with single-neuron resolution becomes closer to reality [1, 2, 3, 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, 6, 7, 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, 10, 11, 12, 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.

Results

Our main result is the development of a theoretical framework to infer hard-to-measure connection strengths in a minimal firing rate model, constrained by spike count statistics from simultaneous array recordings.

We performed simultaneous dual micro-electrode recordings in the olfactory bulb (OB) and the anterior piriform cortex (PC) (see Materials and Methods). First, we use the experimental data to compute population-averaged (across cells or cell pairs) first and second order spike count statistics, comparing across regions (OB or PC) and activity states (spontaneous or stimulus-evoked). We use these statistics to constrain a minimal firing rate model of the coupled OB-PC system, aided by a quick and efficient method for calculating firing statistics without Monte Carlo simulations.

As a test case for our methods, we investigate the structure of four important parameters: within-region inhibitory connection strengths and between-region excitatory connection strengths. We find several relationships that must hold, in order to satisfy all constraints from the experimental data. These results are then validated with a large spiking network of leaky integrate-and-fire (LIF) model neurons.

Consistent Trends in the Experimental Data

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–B). 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 $\nu_{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 $T_{win}$ 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 $\nu_{OB}^{Ev}>\nu_{OB}^{Sp}$ (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 Fig 1C–D, Fig 2, Fig 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–D). 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–D).

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–D 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.

A Minimal Firing Rate Model to Capture Data Constraints

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 [15, 7] (also see [17] who similarly used these same cell types to analyze spiking correlations in larger spiking network models).

We use $j\in\{1,2,3\}$ to denote three OB “cells” and $j\in\{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:

[TABLE]

where $F(x_{k})$ is a transfer function mapping activity to firing rate. Thus, the firing rate is:

[TABLE]

We set the transfer function to $F(X)=\frac{1}{2}\left(1+\tanh((X-0.5)/0.1)\right)$, a commonly used sigmoidal function [14] for all cells; experimental recordings of this function demonstrate it can be sigmoidal [18, 19, 20]. All cells receive noise $\eta_{j}$, the increment of a Weiner process, uncorrelated in time but correlated within a region: i.e. $\langle\eta_{j}(t)\rangle=0$, $\langle\eta_{j}(t)\eta_{j}(t+s)\rangle=\delta(s)$, and $\langle\eta_{j}(t)\eta_{k}(t+s)\rangle=c_{jk}\delta(s)$. We set $c_{jk}$ to:

[TABLE]

The parameters $\mu_{j}$ and $\sigma_{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.,

[TABLE]

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 $\mu_{1},\mu_{2},\mu_{3}$, representing mean input to the olfactory bulb (values given in Table 3). All other parameters were the same for both states. While increasing the input to the I cells in OB in the evoked state ($\mu_{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 Eq 32–49, 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, $g_{jk}$, 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:

• •

No cross-region inhibitory projections: $g_{41}=g_{51}=g_{61}=g_{14}=g_{24}=g_{34}=0$.

• •

Excitatory $PC\rightarrow OB$ output will synapse only onto the inhibitory population: $g_{25}=g_{26}=g_{35}=g_{36}=0$ (see [23]). This reflects experimental evidence that the feedback projections from PC to OB are dominated by inhibition [24, 25].

• •

Excitatory $OB\rightarrow PC$ output will synapse only onto the inhibitory population: $g_{52}=g_{62}=g_{53}=g_{63}=0$. Although this is not anatomically accurate because the mitral/tufted cells also project to the E cells in PC, our goal is to (heuristically) model the prominent role of I cells in PC. Recent work has shown that within PC the recurrent activity is dominated by inhibition [26]. Previous work has also shown that inhibitory synaptic events are much more common in PC and are much easier to elicit [27]. Thus, the connections from excitatory OB to inhibitory PC (Fig 4) should be thought of as the net effect of OB-to-PC connections.

Within OB, there is also excitatory (M/T) input to the inhibitory (granule) cells: $g_{12}=g_{13}=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: $g_{45}=g_{46}=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:

• •

Feedforward inhibitory connections within a population were identical: $gIO\equiv g_{21}=g_{31}$ and $gIP\equiv g_{54}=g_{64}$.

• •

Excitatory connections projecting outward from each region to the other region were identical: $gEO\equiv g_{42}=g_{43}$ and $gEP\equiv g_{15}=g_{16}$.

• •

No within-region excitatory connections; $g_{23}=g_{32}=g_{56}=g_{65}=0$.

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 $20^{4}=1.6\times 10^{5}$ total models. Given each choice of 4-tuple $\{gIO,gIP,gEO,gEP\}$, we computed first- and second-order statistics of both activity $x_{k}$ and firing rates $F(x_{k})$ using the formulas given in Eq 32–49, 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 magnitude of the inhibition within PC, $|gIP|$, is greater than the magnitude of the inhibition within OB, $|gIO|$ (Fig 4A: all dots are above the diagonal line).

• •

The excitation from PC to OB, $gEP$, is generally larger than the excitation from OB to PC, $gEO$ (Fig 4B).

• •

The inhibition within OB is generally weak (Fig 4C: dots are to the left of the vertical line).

• •

The inhibition within PC is generally strong (Fig 4D: dots are to the right of the vertical line).

The statistics computed in Eq 32–49 rely on the assumption that the activity distributions $x_{k}$ 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–D (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.

Admissible firing rate model parameters

How do each of the 12 data constraints (Table 1) restrict the set of possible model parameters? Figure 5 addresses this question in two ways. In Fig 5A, we show, for each constraint, the fraction (as a percent) of all $20^{4}$ 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 $20^{4}$ 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 $\rho_{OB}^{Sp}<\rho_{PC}^{Sp}$ 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 $\rho$ calculation.

Another way to succinctly examine the structure of the four neural attributes: $gIO,gEO,gIP,gEP$ is to consider a matrix:

[TABLE]

where the $j^{th}$ 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,

[TABLE]

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:

[TABLE]

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.

Generality of Firing Rate Model Predictions

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:

• •

Include E-to-E connections from OB to PC (lateral olfactory tract). Also include strong E-to-I drive within PC because input from OB results in balanced excitation and inhibition in PC [26];

• •

Remove the E-to-I connections from OB to PC ($gEO$ in the firing rate model) so that the recurrent activity in PC is driven by E inputs along the lateral olfactory tract;

• •

Remove the direct sensory input to I cells in OB since granule cells do not receive direct sensory input [23];

• •

Include substantial recurrent E-to-E connections within PC (see Table 4 for strength relative to other connections).

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.).

Results are Validated in a Spiking LIF Network

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:

[TABLE]

See Fig 6 or Eq Leaky Integrate-and-Fire Model of the generic Cortical–Cortical Circuit–Leaky Integrate-and-Fire Model of the generic Cortical–Cortical Circuit 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 9, and other standard parameter values (see Table 4) in a typical LIF model, we were able to easily satisfy all 12 constraints: see Table 2 and Fig 6.

While the firing rates in the LIF network (Table 2) 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].

Results of Violating Derived Relationships Between Conductance Strengths

What happens in the full LIF spiking network when the derived relationships between the conductance strengths are violated? Since the minimal firing rate model is different than the detailed spiking model in many ways, we do not expect the relationships between the conductance strengths to hold precisely. However, the minimal firing rate model is still useful in providing intuition for what would otherwise be a complicated network with a high-dimensional parameter space. We now demonstrate that when the relationships derived in firing rate model are violated, a subset of the constraints in the experimental data (Table 1) will no longer be satisfied in the large spiking network.

Because our network is heterogeneous, our ability to subsample cell pairs is limited, relative to a homogeneous network of the same size. Also, computation for even a single parameter set in the spiking network require enormous computing resources. Thus, we cannot exhaustively explore the parameter space; indeed, the purpose of the reduction method of the firing rate model is to probe large dimensions quickly. Instead, we perform three tests that violate the firing rate model results:

Make $gIO>gIP$ by setting $gIO=20$ and $gIP=7$. 2. 2.

Make $gEO>gEP$ by setting $gEO=15$ and $gEP=1$ 3. 3.

Make $gEP$ and $gIP$ relatively smaller by setting $gEP=10$ and $gIP=10$

The original values (used in Fig 6) for these parameters were given in Eq 9.

The result of Test 1 is that 2 of the 12 constraints are violated (see Fig S14 in S2 Text.); most importantly stimulus-induced decorrelation of the PC cells, which is particularly important in the context of coding, was not present. In addition, the evoked PC correlation is larger than evoked OB correlation, violating another constraint.

The result of Test 2 is that 3 of the 12 constraints are violated (see Fig S15 in S2 Text.). The evoked PC correlation is larger than evoked OB correlation, and both the variance and covariance in PC are larger than the corresponding quantities in OB in the evoked state, which is not consistent with our data.

The result of Test 3 is that 3 of the 12 constraints are violated: they are the same constraints that are violated in Test 2, despite quantitative differences in the statistics (see Fig S16 in S2 Text.). The stimulus-induced decorrelation of the PC cells does not hold for small windows, but this is also observed in our data (Fig 2A), so we do not formally count this as a clear violation of data constraints. However, Test 1 and Test 3 show that strong PC inhibition is key for stimulus-induced decorrelation [7, 31, 32, 17, 33, 34, 22].

Discussion

As electrophysiological recording technology advances, there will be more datasets with simultaneous recordings of neurons, spanning larger regions of the nervous system. Such networks are inherently high-dimensional, making mechanistic analyses generally intractable without fast and reasonably accurate approximation methods. We have developed a computational reduction method for a multi-population firing rate model [14] that enables analysis of the spiking statistics. Our work specifically enables theoretical characterizations of an important, yet hard-to-measure quantity – synaptic connection strength – using easy-to-measure spiking statistics. The method is computationally efficient, is validated with Monte Carlo simulations of spiking neural networks, and can provide insight into network structure.

We applied our computational methods to simultaneous dual-array recordings in two distinct regions of the olfactory system: the olfactory bulb (OB) and anterior piriform cortex (PC). Our unique experimental dataset enables a detailed analysis of the first- and second-order spike count statistics in two activity states, and a comparison of how these statistics are related between OB and PC cells. We found twelve (12) consistent trends that held across four odors in the dataset (Table 1), and sought to identify what neural network attributes would account for these trends. We focused on four important network attributes, specifically the conductance strengths in the following connections: feedforward inhibition within OB and within PC, excitatory projections from OB to PC neurons, and finally excitatory projections from PC to OB. Our reduced firing rate model predicts several relationships that are then verified with a more detailed spiking network model, specifically: i) inhibition within the 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) connections that originate within PC have to relatively strong compared to connections that originate within OB. These results make a strong prediction that to the best of our knowledge is new and might be testable with simultaneous patch-clamp recordings.

In principle our theory could be used to study the structure of other network features such as background correlation, noise level, transfer function, etc.. It is straightforward mathematically to incorporate other desired neural attributes (with the caveat of perhaps increasing the overall number of equations and terms in the approximations) without changing the basic structure of the framework. Here we have focused on the role of the strength of synaptic coupling; of course, other neural attributes can affect spike statistics (in particular, spike count correlation [35, 7]), some of which can conceivably change with stimuli. Spike count correlations can depend on intrinsic neural properties [36, 37, 38, 39, 40, 41], network architecture [42, 30, 43] and synaptic inputs [22, 31, 34, 33, 32, 17] (or combinations of these [44, 45, 46, 8]), plasticity [47], as well as top-down mechanisms [48, 49, 50]. Thus, correlation modulation is a rich and deep field of study, and we do not presume our result is the only plausible explanation for spike statistics modulation.

Although the minimal firing rate model did not include certain anatomical connections that are known to exist (e.g., recurrent excitation in the PC), the model is meant for deriving qualitative principles rather than precise quantitative modeling of the pathway. We based our simplifications on insights from recent experimental work: recent slice physiology work has shown that within PC, recurrent activity is dominated by inhibition [26]; previous work has also shown that inhibitory synaptic events are much more common (than excitatory synaptic events) in PC and are much easier to elicit [27]. Thus, the connection from excitatory OB cells to inhibitory PC cells ($gEO$ in Fig 4) should be thought of as the net effect of these connections along the lateral olfactory tract. Other theoretical analyses of effective feedforward inhibitory networks have also neglected anatomical E-to-E connections [32, 17]. Furthermore, this minimal model was validated with a more realistic, recurrently coupled spiking network, which did include within-region excitatory connections (see Fig 6 and Fig S14–S16 in S2 Text., as well as S3 Text.).

We have only focused on first- and second-order firing statistics, even though in principle other, higher-order statistics may be important [51, 52, 53]. If downstream neurons use a linear decoding scheme, then first- and second-order spiking statistics are sufficient in quantitative measures of neural coding [54, 55]. It is currently unknown whether downstream neurons decode olfactory signals with a nonlinear decoder, but there is evidence in other sensory systems that second-order statistics are sufficient [13]. Recent work has shown conflicting results for coding in olfactory bulb; one study found that decoding an odor in the presence of other odors might be more efficient using nonlinear decoding [29], but another has shown that linear decoding is still plausible [56].

A second reason to neglect higher-order statistics is suggested by Fig. 5, where we show how the various data constraints narrow the scope of plausible models. Here, we saw that even with first and second- order statistics, only 1% of the parameter sets satisfy the data constraints; including more constraints would limit the space further. In order to usefully include higher-order constraints, we would need to use a more detailed model and/or larger parameter spaces.

As a test case for our method, we used recordings from anesthetized animals. The absence of breathing in tracheotomized rats in these experiments is only an approximation to olfactory processing in awake animals. However, there is a benefit to tracheotomized animals: the complex temporal firing patterns are removed, so that firing statistics are closer to stationarity. In principle, we can incorporate breathing dynamics into our framework by including an oscillatory forcing term in Eq 1; this will be the subject of future work. In support of this simplification, we note that there is evidence that in the anterior piriform cortex, spike count — rather than the timing — is most consequential for odor discrimination [16]. However, other studies have reported that timing of the stimuli in the olfactory bulb is important: [57, 58, 29] showed decoding performance is best at the onset of odors in mammals and worsens as time proceeds, whereas [59] found that decoding performance improved with time in zebrafish. These important issues are beyond the scope of this current study.

Relationship to Other Reduction Methods

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, 64, 62, 65, 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.

Conclusion

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.

Materials and Methods

Electrophysiological Recordings

Subjects. 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). Experimental data was obtained from one adult male rat (289 g ; Rattus Norvegicus, Sprague-Dawley outbred, Harlan Laboratories, TX, USA) housed in an environment of controlled humidity (60%) and temperature (23*∘*C) with 12 h light-dark cycles. The experiments were performed in the light phase.

Anesthesia. Anesthesia was induced with isoflurane inhalation and maintained with urethane (1.5 g/kg body weight (bw) dissolved in saline, intraperitoneal injection (ip)). Dexamethasone (2 mg/kg bw, ip) and atropine sulphate (0.4 mg/kg bw, ip) were administered before performing surgical procedures.

Double tracheotomy surgery. To facilitate ortho- and retronasal delivery of the odorants a double tracheotomy surgery was performed as described previously [76]. This allowed for the rat to sniff artificially while breathing naturally through the trachea bypassing the nose. A Teflon tube (OD 2.1 mm, upper tracheotomy tube) was inserted 10 mm into the nasopharynx through the rostral end of the tracheal cut. Another Teflon tube (OD 2.3 mm, lower tracheotomy tube) was inserted in to the caudal end of the tracheal cut to allow breathing. Both tubes were fixed and sealed to the tissues using surgical thread. Local anesthetic (2% Lidocaine) was applied at all pressure points and incisions. Throughout the surgery and electrophysiological recordings rats’ core body temperature was maintained at 37*∘*C with a thermostatically controlled heating pad.

Craniotomy surgery. Subsequently, a craniotomy surgery was performed on the dorsal surface of the skull at two locations, one over the right Olfactory Bulb (2 mm $\times$ 2 mm, centered 8.5 mm rostral to bregma and 1.5 mm lateral from midline) and the other over the right anterior Pyriform Cortex (2 mm $\times$ 2 mm, centered 1.5ṁm caudal to bregma and 5.5 mm lateral from midline).

Presentation of ortho- and retronasal odorants. The bidirectional artificial sniffing paradigm previously used for the presentation of ortho- and retronasal odorants [76] were slightly modified such that instead of a nose mask a Teflon tube was inserted into the right nostril and the left nostril was sealed by suturing. The upper tracheotomy tube inserted into the nasopharynx was used to deliver odor stimuli retronasally (Fig 1. We used two different odorants, Hexanal (Hexa) and Ethyl Butyrate (EB) by both ortho- and retronasal routes, there by constituting 4 different odor stimuli. Each trial consisted of 10 one-second pulse presentations of an odor with 30 second interval in between two pulses, and 2-3 min in between two trials.

Electrophysiology. Extracellular voltage was recorded simultaneously from OB and aPC using two different sets of 32-channel microelectrode arrays (MEAs). (OB: A4x2tet, 4 shanks x 2 iridium tetrodes per shank, inserted 400 $\mu$m deep from dorsal surface; aPC: Buzsaki 32L, 4 shanks x 8 iridium electrode sites per shank, 6.5 mm deep from dorsal surface; NeuroNexus, MI, USA). Voltages were measured with respect to an AgCl ground pellet placed in the saline-soaked gel foams covering the exposed brain surface around the inserted MEAs. Voltages were digitized with 30 kHz sample rate as described previously [77] using Cereplex + Cerebus, Blackrock Microsystems (UT, USA).

Recordings were filtered between 300 and 3000 Hz and semiautomatic spike sorting was performed using Klustakwik software, which is optimized for the types of electode arrays used here [78]. After automatic sorting, each unit was visually inspected to ensure quality of sorting.

Data processing

After the array recordings were spike sorted to identify activity from distinct cells, we further processed the data as follows:

• •

We computed average firing rate for each cell, where the average was taken over all trials and over the entire trial length (i.e., not distinguishing between spontaneous and evoked periods); units with firing rates below 0.008 Hz and above 49 Hz were excluded.

• •

When spike times from the same unit were within 0.1 ms of each other, only the first (smaller) of the spike time was used and the subsequent spike times were discarded

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 $T_{win}$. For a particular time window, the $j^{th}$ neuron has a spike count instance $N_{j}$ in the time interval $[t,t+T_{win})$:

[TABLE]

The spike count correlation between cells $j$ and $k$ is given by:

[TABLE]

where the covariance of spike counts is:

[TABLE]

Here $n$ is the total number of observations of $N_{j/k}$, and $\mu(N_{j}):=\frac{1}{n}\sum N_{j}$ is the mean spike count across $T_{win}$-windows and trials. The correlation $\rho_{T}$ is a normalized measure of the the trial-to-trial variability (i.e., noise correlation), satisfying $\rho_{T}\in[-1,1]$; it is also referred to as the Pearson’s correlation coefficient. For each cell pair, the covariance $\text{Cov}(N_{j},N_{k})$ and variance $\text{Var}(N_{j})$ 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:

[TABLE]

In principle, any of the statistics defined here might depend on the time $t$ as well as time window size $T_{win}$; here, we assume that Var, Cov, $FF$, and $\rho_{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 windows111an exception is when $T_{win}=2\,$s; in the evoked state, there is only 1 window per trial; the exact number depends on the state, the value of $T_{win}$, and whether disjoint or overlapping windows are used. In this paper we use overlapping windows by half the length of $T_{win}$222e.g. if the trial length is 2 s and $T_{win}=1\,$s, then there are 3 total windows per trial: [0 s, 1 s], [0.5 s, 1.5 s], and [1 s, 2 s] to calculate the spiking statistics. The results are qualitatively similar for disjoint windows and importantly the relationships/constraints are the same with disjoint windows. We limit the size of $T_{win}\leq 2\,$s because this is the maximum duration of the evoked state, within each trial.

The average spike count $\mu(N_{j})$ of the $j^{th}$ neuron with a particular time window $T_{win}$ is related to the average firing rate $\nu_{j}$ of that neuron:

[TABLE]

Firing Rate Model

Recall that the activity in each representative cell is modeled by:

[TABLE]

where $F(x_{k})$ is a transfer function mapping activity to firing rate. Thus, the firing rate is:

[TABLE]

The index of each region is denoted as follows: $j\in\{1,2,3\}$ for the 3 OB cells, and $j\in\{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 $\sigma_{1}=\sigma_{2}=\sigma_{3}=\sigma_{OB}$ and $\sigma_{4}=\sigma_{5}=\sigma_{6}=\sigma_{PC}$ (see Table 3).

In the absence of coupling (i.e. $g_{jk}=0$), any pair of activity variables, $(x_{j},x_{k})$, are bivariate normally distributed because the equations:

[TABLE]

describe a multi-dimensional Ornstein-Uhlenbeck process [79]. Note that we have re-written $\eta_{j/k}(t)$ as sums of independent white noise processes $\xi(t)$, which is always possible for Gaussian white noise. Since x_{j}(t)=\frac{1}{\tau}\int_{0}^{t}e^{-(t-u)/\tau}\Big{[}\mu_{j}+\sigma_{j}\eta_{j}(u)\Big{]}\,du, we calculate marginal statistics as follows:

[TABLE]

[TABLE]

A similar calculation shows that in general we have:

[TABLE]

Thus, $(x_{j},x_{k})\sim\mathcal{N}\left(\left(\begin{smallmatrix}\mu_{j}\\ \mu_{k}\end{smallmatrix}\right),\frac{1}{2\tau}\left(\begin{smallmatrix}\sigma^{2}_{j}&\sigma_{j}\sigma_{k}c_{jk}\\ \sigma_{j}\sigma_{k}c_{jk}&\sigma^{2}_{k}\end{smallmatrix}\right)\right)$.

To simplify notation, we define:

[TABLE]

With coupling, an exact expression for a joint distribution for $(x_{1},x_{2},x_{3},x_{4},x_{5},x_{6})$ 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 $\tau=1$) for each of the 3000 realizations. The activity $x_{j}$ 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.

Approximation of Firing Statistics in the Firing Rate Model

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 $\langle x_{j}\rangle$, mean firing rate $\langle F(x_{j})\rangle$, variance and covariances of both: $\langle x_{j}x_{k}\rangle$ and $\langle F(x_{j})F(x_{k})\rangle$. For a simpler exposition, we have only included twelve synaptic connections; we have excluded self (autaptic) connections and E$\to$E connections.

An equation for each statistic can be derived by first writing Eq 15 as a low-pass filter of the right-hand-side:

[TABLE]

We then take expectations, letting $t\rightarrow\infty$, we have:

[TABLE]

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:

[TABLE]

and $N_{j}$ denotes the random variable $\int_{0}^{t}\sigma_{j}\eta_{j}(u)e^{-(t-u)/\tau}\,du$, which is by itself normally distributed with mean 0 and variance $\sigma_{j}^{2}\tau/2$.

The first assumption, Eq 25, states that time-average of $F(x_{j}(t))$ multiplied by an exponential function (low-pass filter) is equal to the expected value scaled by $\tau/2$; the second and third, Eq 27 and Eq 29, address $N_{j}$ and $F(x_{k}(t))$, for $j\not=k$ and $j=k$ respectively (similarly for Eq 30).

In all of the definitions for the expected values with $\rho_{2D}$, note that the underlying correlation correlation $c_{jk}$ depend on the pair of interest $(j,k)$. Finally, we assume that the activity variables $(x_{j},x_{k})$ 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, $(x_{j},x_{k})$ are bivariate normal random variables.

The resulting approximations for the mean activity are:

[TABLE]

The resulting approximation to the variances of the mean activity are:

[TABLE]

In Eq 32–43, all of the Var and Cov are with respect to $Y\sim\mathcal{N}(0,1)$ (for Var) and $(Y_{1},Y_{2})\sim\mathcal{N}\left(\left(\begin{smallmatrix}0\\ 0\end{smallmatrix}\right),\frac{1}{2}\left(\begin{smallmatrix}1&c_{jk}\\ c_{jk}&1\end{smallmatrix}\right)\right)$ (for Cov); both are easy to calculate. The value $c_{jk}$ depends on the pairs; for example in Eq 39, the $\rho_{2D}$ has $c_{jk}=c_{OB}$, the background correlation value in the olfactory bulb but in Eq 38, the Cov term is with respect to $\rho_{2D}$ with $c_{jk}=c_{PC}$, 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 Cov of interest are:

[TABLE]

where

[TABLE]

Iteration procedure to solve for the approximate statistics self-consistently

Based on the approximations and resulting equations described in the previous section, our objective is to solve for the statistics of $x_{j}$ self-consistently. Once these are determined, the statistics of the firing rates $F(x_{j})$ are approximated with the same pairwise normal assumption on $(x_{j},x_{k})$; we are not assuming that $(F(x_{j}),F(x_{k}))$ are bivariate normal random variables.

We use a simple iterative procedure to solve the system of coupled algebraic expression for the statistics of $x_{j}$. We first solve the system in the absence of coupling (i.e. Eq 19, 20), and use these values to start the iteration; at each step, the formulas for the means (Eq 32–37), variances (Eq 38–43), and covariances (Eq 44–49) are recalculated numerically, using the results of the previous step. The iteration stops once all 18 statistical quantities of the activity match up to a relative tolerance of $10^{-6}$ (convergence), or after 50 total iterations (non-convergence). The result with a given parameter set can either be: i) convergence, ii) non-convergence, iii) a pair of statistics with invalid covariance (non-positive definite covariance matrix), which is checked after i) and ii). We only consider parameter sets where the iteration has converged and all of the covariances are valid, after which we determine whether the constraints are satisfied.

One subtle point is that we did not use any of the numerically calculated Cov values in the bivariate normal distributions $\rho_{2D}$; rather, the correlation value is always $c_{jk}$ which is either 0, $c_{OB}$, or $c_{PC}$ depending on the pair. In principle, one can use a fully iterative procedure where the formulas for the Cov (Eq 44–49) are used in $\rho_{2D}$; however, we found that the resulting covariance matrices (for $\rho_{2D}$) can fail to be positive semi-definite. Handling this case requires additional code in the program and slower calculations for each parameter set, which detracts from the purpose of our method. We checked some parameter sets comparing the results of the two procedures, and the results are quantitatively similar.

The standard normal $\rho_{SN}$ and bivariate $\rho_{2D}$ PDFs have state variable(s) $y_{1,2}$ discretized from -3 to 3 with a mesh size of 0.01; integrals in Eq 32–49 are computed using the trapezoidal rule.

Simplified network with four coupling parameters

To further simplify the network, we:

• •

set $\tau=1$,

• •

assume feedforward inhibitory connections within a region have the same strength: $g_{21}=g_{31}=:gIO$ and $g_{54}=g_{64}=:gIP$,

• •

assume cross-region excitatory connections are equal from the presynaptic cell, i.e., $g_{15}=g_{16}=:gEP$ and $g_{42}=g_{43}=:gEO$.

• •

assume $\sigma_{1}=\sigma_{2}=\sigma_{3}=:\sigma_{OB}$ and $\sigma_{4}=\sigma_{5}=\sigma_{6}=:\sigma_{PC}$

• •

assume $g_{12}=g_{13}=g_{45}=g_{46}=:g_{\epsilon}=0.1$

Now there are only 4 variable coupling parameters: $gIO$, $gEO$, $gIP$, $gEP$.

The above formulas for the statistics of $x_{j}$ reduce to:

[TABLE]

the variances are:

[TABLE]

the covariances are:

[TABLE]

See Eq 50 for the definition of $\mathcal{C}$.

Leaky Integrate-and-Fire Model of the OB–PC Circuit

We use a generic spiking neural network model of leaky integrate-and-fire neurons to test the results of the theory. There were $N_{OB}=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\in\{1,2,\dots,N_{OB}\}$:

[TABLE]

The conductance values in the first equation $g_{k,XI}$, $g_{k,XE}$, and $g_{k,XPC}$ depend on the type of neuron $v_{k}$ ($X\in\{E,I\}$). The last conductance, $g_{X,PC}(t-\tau_{\Delta,PC})(v_{k}-\mathcal{E}_{E})$, models the excitatory presynaptic input (feedback) from the PC cells with a time delay of $\tau_{\Delta,PC}$. The conductance variables $g_{k,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, $\mathcal{E}_{I}$, and $\mathcal{E}_{E}$, respectively with $\mathcal{E}_{I}<0<\mathcal{E}_{E}$ (the voltage is scaled to be dimensionless, see Table 4). $\xi_{k}(t)$ are uncorrelated white noise processes and $\xi_{o}(t)$ is the common noise term to all $N_{OB}$ cells.

The second equation describes the refractory period at spike time $t^{*}$: when the neuron’s voltage crosses threshold $\theta_{j}$ (see below for distribution of thresholds), the neuron goes into a refractory period for $\tau_{ref}$, after which we set the neuron’s voltage to 0.

The parameter $\gamma_{XY}$ gives the relative weight of a connection from neuron type $Y$ to neuron type $X$; the parameter $p_{XY}$ is the probability that any such connection exists ($X,Y\in\{E,I\}$). $G_{k}$ 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 74 and depend on whether $k\in\{E,I\}$.

Finally, two of the parameters above can be equated with coupling parameters in the reduced model:

[TABLE]

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 $N_{PC}=100$ total PC cells, of which 80% were excitatory and 20% inhibitory. The equations, indexed by $j\in\{1,2,\dots,N_{PC}\}$ are:

[TABLE]

Excitatory synaptic input from the OB cells along the lateral olfactory tract is modeled by: $g_{X,OB}(t-\tau_{\Delta,OB})(v_{j}-\mathcal{E}_{E})$. The common noise term for the PC cells $\xi_{p}(t)$ is independent of the common noise term for the OB cells $\xi_{o}(t)$. Two of the parameters above can be equated with coupling parameters in the reduced model:

[TABLE]

The values of the parameters that were not stated in Table 4 were varied in the paper:

[TABLE]

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 $\mu_{PC}=0$ in the spontaneous state to $\mu_{PC}=0.4$ in the evoked state) as well as to OB cells; $\mu_{OB}=0.6$ in the spontaneous state to $\mu_{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 $\theta_{j}$ in the following way. Both OB and PC cells had the following distributions for $\theta_{j}$:

[TABLE]

where $\mathcal{N}$ is normal distribution with mean $-\sigma^{2}_{\theta}/2$ and standard deviation $\sigma_{\theta}$, so that $\{\theta_{j}\}$ has a log-normal distribution with mean 1 and variance: $e^{\sigma_{\theta}^{2}}-1$. We set $\sigma_{\theta}=0.1$, which results in firing rates ranges seen in the experimental data. Since the number of cells are modest with regards to sampling ($N_{OB}=100$, $N_{PC}=100$), we evenly sampled the log-normal distribution from the 5th to 95th percentiles (inclusive).

We remark that the synaptic delays of $\tau_{\Delta,PC}$ and $\tau_{\Delta,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 [80].

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.

Supporting Information

S1 Text.

Experimental Data Statistics by Odor. This file shows the trial-averaged spiking statistics of the experimental data dissected by a specific odor. Contains Figs. S1-S8, and Tables S1-S2.

S2 Text.

Supplementary Figures for the Main Modeling. This file contains supplemental figures from modeling and analysis. Contains Figs. S9-S16.

S3 Text.

Supplementary Material: Cortical-Cortical Network. This file contains supplemental modeling results on a generic cortical-cortical coupled network. Contains Figs. S17-S21, and Table S3.

S1 Table.

Average population firing rate by odor and activity state.

S2 Table.

Standard deviation of population firing rate by odor and activity state.

S3 Table.

Fixed parameters for the LIF Cortical-cortical model.

S1 Figure.

Experimental statistics by odor and activity state: Fano Factor.

S2 Figure.

Experimental statistics by odor and activity state: spike count variance.

S3 Figure.

Experimental statistics by odor and region: Fano Factor.

S4 Figure.

Experimental statistics by odor and region: spike count variance.

S5 Figure.

Experimental statistics by odor and activity state: spike count correlation.

S6 Figure.

Experimental statistics by odor and activity state: spike count covariance.

S7 Figure.

Experimental statistics by odor and region: spike count correlation.

S8 Figure.

Experimental statistics by odor and region: spike count covariance.

S9 Figure.

Cross-region correlations are smaller than within-region correlations.

S10 Figure.

Fast analytic approximation accurately captures statistics of a multi-population firing rate model.

S11 Figure.

Experimental observations constrain conductance parameters in analytic model.

S12 Figure.

Analytic approximation results are robust to choice of transfer function.

S13 Figure.

Mean input to PC must increase in the evoked state.

S14 Figure.

Violating derived relationship $gIO<gIP$ results in statistics that are inconsistent with experimental observations.

S15 Figure.

Violating derived relationship $gEP>gEO$ results in statistics that are inconsistent with experimental observations.

S16 Figure.

Violating derived relationship $gEP,gIP\gg gEO,gIO$ results in statistics that are inconsistent with experimental observations.

S17 Figure.

Minimal firing rate model to analyze synaptic conductance strengths.

S18 Figure.

Detailed spiking LIF model confirms the results from analytic rate model.

S19 Figure.

Violating derived relationship $|gI1|<|gI2|$ results in statistics that are inconsistent with experimental observations.

S20 Figure.

Violating derived relationship $gE2,gI2\gg gE1,gI1$ results in statistics that are inconsistent with experimental observations.

S21 Figure.

iolating derived relationship $gE2>gE1$ results in statistics that are inconsistent with experimental observations.

Author Contributions

Conceived and designed research: AKB CL. Derived expressions in theoretical methods: CL. Analyzed the data: AKB SHG WLS CL. Conceived and designed electrophysiological experiments: SHG WLS. Wrote the paper: AKB SHG WLS CL.

Minimal Firing Rate Model

The minimal firing rate model results in Fig S17 have the same parameters and configuration as in the main text except the E to I connections within a region are omitted. The derived relationships are qualitatively the same as in the main text:

[TABLE]

Leaky Integrate-and-Fire Model of the generic Cortical–Cortical Circuit

We use a generic spiking neural network model of leaky integrate-and-fire neurons to test the results of the theory again. The following model is very much like the LIF model in the main text, with the main differences being in the network connection strengths and the size of C1 (60 here instead of 100 for OB in the main text). There were $N_{C1}=60$ total C1 cells, of which we set 80% (48) to be excitatory and 20% (12) inhibitory. The equations for the C1 cells are, indexed by $k\in\{1,2,\dots,N_{C1}\}$:

[TABLE]

The conductance values in the first equation $g_{k,XI}$, $g_{k,XE}$, and $g_{k,XC2}$ depend on the type of neuron $v_{k}$ ($X\in\{E,I\}$). The last conductance, $g_{X,C2}(t-\tau_{\Delta,C2})(v_{k}-\mathcal{E}_{E})$, models the excitatory presynaptic input (feedback) from the C2 cells with a time delay of $\tau_{\Delta,C2}$. The conductance variables $g_{k,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, $\mathcal{E}_{I}$, and $\mathcal{E}_{E}$, respectively with $\mathcal{E}_{I}<0<\mathcal{E}_{E}$ (the voltage is scaled to be dimensionless, see Table S3). $\xi_{k}(t)$ are uncorrelated white noise processes and $\xi_{o}(t)$ is the common noise term to all $N_{C1}$ cells.

The second equation describes the refractory period at spike time $t^{*}$: when the neuron’s voltage crosses threshold $\theta_{j}$ (see below for distribution of thresholds), the neuron goes into a refractory period for $\tau_{ref}$, after which we set the neuron’s voltage to 0.

The parameter $\gamma_{XY}$ gives the relative weight of a connection from neuron type $Y$ to neuron type $X$; the parameter $p_{XY}$ is probability that any such connection exists ($X,Y\in\{E,I\}$). $G_{k}$ 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 74 and depend on whether $k\in\{E,I\}$.

Finally, two of the parameters above can be equated with coupling parameters in the reduced model:

[TABLE]

which are dimensionless scale factors for the synaptic conductances.

The C2 cells had similar functional form but with different parameters (see Table S3 for parameter values). We modeled $N_{C2}=100$ total C2 cells, of which 80% were excitatory and 20% inhibitory. The equations, indexed by $j\in\{1,2,\dots,N_{C2}\}$ are:

[TABLE]

Excitatory synaptic input from the C1 cells along the lateral olfactory tract is modeled by: $g_{X,C1}(t-\tau_{\Delta,C1})(v_{j}-\mathcal{E}_{E})$. The common noise term for the C2 cells $\xi_{p}(t)$ is independent of the common noise term for the C1 cells $\xi_{o}(t)$. Two of the parameters above can be equated with coupling parameters in the reduced model:

[TABLE]

The values of the parameters that were not stated in Table S3 were varied:

[TABLE]

To model two activity states, we allowed mean inputs to vary (see Table S3). In contrast to the reduced model, we increased both inputs to C2 cells (from $\mu_{C2}=0$ in the spontaneous state to $\mu_{C2}=0.4$ in the evoked state) as well as to C1 cells (from $\mu_{C1}=0.6$ in the spontaneous state to $\mu_{C1}=0.9$ in the evoked state).

Finally, we model heterogeneity by setting the threshold values $\theta_{j}$ in the following way. Both C1 and C2 cells had the following distributions for $\theta_{j}$:

[TABLE]

where $\mathcal{N}$ is normal distribution with mean $-\sigma^{2}_{\theta}/2$ and standard deviation $\sigma_{\theta}$, so that $\{\theta_{j}\}$ has a log-normal distribution with mean 1 and variance: $e^{\sigma_{\theta}^{2}}-1$. We set $\sigma_{\theta}=0.1$, which results in firing rates ranges seen in the experimental data. Since the number of cells are modest with regards to sampling ($N_{C1}=60$, $N_{C2}=100$), we evenly sampled the log-normal distribution from the 5th to 95th percentiles (inclusive).

Violating Derived Relationships Between Conductance Strengths

Similar to the main text, we demonstrate here that violating the relationships derived in the main text results in a subset of the 12 constraints in the experimental data no longer being satisfied in the full spiking network.

Due to the large amount of computing resources required, we cannot exhaustively explore the parameter space; recall that the purpose of the method we developed in the minimal firing rate model is for faster computation. Instead, we distill results into three tests that are exactly the same as in the main text:

Make $gI1>gI2$ by setting $gI1=20$ and $gI2=7$. 2. 2.

Make $gE1>gE2$ by setting $gE1=15$ and $gE2=1$ 3. 3.

Make $gE2$ and $gI2$ relatively smaller by setting $gE2=10$ and $gI2=10$

The result of Test 1 is that 7 of the 12 constraints are violated (see Fig S19); most importantly stimulus-induced decorrelation of the C2 cells, which is particularly important in the context of coding, was not present. In addition, the C2 firing rates are larger than the C1 firing rates in both states, the evoked C2 correlation is larger than evoked C1 correlation, the spontaneous C2 Fano Factor is larger than spontaneous C1 Fano Factor, and both the variance and covariance of C2 is larger than C1 in the evoked state (all of which violate the constraints from our data).

The result of Test 3 is that 4 of the 12 constraints are violated (see Fig S20), including again stimulus-induced decorrelation of the C2 cells. The evoked C2 correlation is larger than evoked C1 correlation, and both the variance and covariance of C2 are larger than the corresponding quantities in C1 in the evoked state. Both these two tests (1 and 3) indicate that these two qualitative relationships (stronger effective inhibition within C2 and stronger effective presynaptic inputs from C2) are robust with respect to both the detailed LIF spiking model and the minimal firing rate model.

The result of Test 2 is not as straightforward as the others. We did not exhaustively search parameter space due to the vast computational resources this would require, but in several parameter sets with $gE1>gE2$, we found the resulting network statistics could still satisfy all of the constraints (e.g., with $gE1=15$ and $gE2=1$, as well as with $gE1=20$ and $gE2=1$). The reason for this may be that in the two coupled recurrent networks we chose very different $gI1$ and $gI2$ values to begin with (7 and 20, respectively), and would thus require $gE1$ and $gE2$ to be significantly different to counter-balance this. Also, notice in the minimal firing rate model results in Fig 4B that there are a significant number of red dots below the diagonal, indicating that the relationship $gE2>gE1$ does not have to strictly hold. However, we did find a condition where this test demonstrates the value of the minimal firing rate model; we changed $\tilde{c}_{C1}$ from $0.5$ to $0.6$ (recall $\tilde{c}_{C2}=0.8$). (Note that in the minimal firing model that $c_{C1}=0.3$ and $c_{C2}=0.35$, relatively close in value.) The result of Test 2 ($gE1=15$ and $gE2=1$) with $\tilde{c}_{C1}=0.6$ is that one constraint is violated: $\rho^{Sp}_{C2}$ is no longer less than $\rho^{Sp}_{C1}$ (see Fig S21). This suggests that the relationship that $gE1>gE2$ is not as robust as the others and can be violated.