Differences in olfactory bulb mitral cell spiking with ortho- and retronasal stimulation revealed by data-driven models

The majority of olfaction studies focus on orthonasal stimulation where odors enter via the front nasal cavity, while retronasal olfaction, where odors enter the rear of the nasal cavity during feeding, is understudied. The coding of retronasal odors via coordinated spiking of neurons in the olfactory bulb (OB) is largely unknown despite evidence that higher level processing is different than orthonasal. To this end, we use multi-electrode array in vivo recordings of rat OB mitral cells (MC) in response to a food odor with both modes of stimulation, and find significant differences in evoked firing rates and spike count covariances (i.e., noise correlations). Differences in spiking activity often have implications for sensory coding, thus we develop a single-compartment biophysical OB model that is able to reproduce key properties of important OB cell types. Prior experiments in olfactory receptor neurons (ORN) showed retro stimulation yields slower and spatially smaller ORN inputs than with ortho, yet whether this is consequential for OB activity remains unknown. Indeed with these specifications for ORN inputs, our OB model captures the salient trends in our OB data. We also analyze how first and second order ORN input statistics dynamically transfer to MC spiking statistics with a phenomenological linear-nonlinear filter model, and find that retro inputs result in larger linear filters than ortho inputs. Finally, our models show that the temporal profile of ORN is crucial for capturing our data and is thus a distinguishing feature between ortho and retro stimulation, even at the OB. Using data-driven modeling, we detail how ORN inputs result in differences in OB dynamics and MC spiking statistics. These differences may ultimately shape how ortho and retro odors are coded.


Introduction
Olfactory processing naturally occurs in two distinct modes: orthonasal (ortho) where odors enter the front of the nasal cavity and retronasal (retro) where odors enter the rear through the throat. Orthonasal olfaction is essential for avoiding predators [1,2], social interactions, and finding food, and has been studied most extensively in olfaction research. Retronasal olfaction is far less studied, but has a critical role in eating behaviors as chewed foods generate odorants that enter the nasal cavity upon exhalation. Retronasal olfaction drives flavor perception [3][4][5] and aids in avoiding harmful foods. Moreover, studies have shown that olfactory dysfunction with food odors is directly linked to obesity [6][7][8]. Previous studies have reported differences in cortical fMRI BOLD signals for ortho versus retro stimuli [9] and recent evidence suggests that food odors are easier to recognize when delivered retronasally versus orthonasally [10]. Calcium imaging studies have shown that the input to olfactory bulb from the nose differs for ortho versus retro stimulation [11]. However, the neural mechanisms that differentiate ortho versus retro olfactory processing at the level of spiking activity in olfactory bulb remain unknown.
Odor information is primarily processed in the olfactory bulb (OB) and then subsequently relayed to cortical areas via mitral cell (MC) (and tufted cell) spiking. Thus, any differences in MC spiking between ortho and retro are related to both the efficiency and accuracy [12][13][14][15] of odor coding, but any such differences are largely unknown. Presynaptic to the OB are olfactory receptor neurons (ORNs) whose activity is known to differ for ortho versus retro stimulation, as observed in prior imaging studies with fMRI [16], calcium imaging [11], and optical imaging in transgenic mice [17]. These and other prior studies [18][19][20] suggest that ORN synaptic inputs is a key factor for differences in OB activity. The two routes of stimulation make contact on different locations of the olfactory epithelium (shown in light green in Fig 1 experimental  diagram) and thus activate different ORN receptor types within the epithelium. However, the implications of these differences in ORN activity for MC spiking have yet to be explored.
We perform in vivo recordings of rat OB mitral cells using multi-electrode arrays with a food odor (Ethyl Butyrate) stimulus, delivered by both modes of stimulation, to determine whether differences exist. We find significant differences in odor-evoked MC spiking with ortho versus retro stimulation in both firing rate (larger with retro) and spike count covariance (larger with ortho). However, understanding how retro stimulation can elicit both larger firing rates and smaller co-variablity than ortho is generally difficult in recurrent networks because of the numerous attributes that shape spike statistics [21][22][23][24]. Additionally, dissecting how components of ORN inputs alter OB spiking is difficult experimentally due to the complexity of both the recurrent circuitry in the OB [25,26] and resulting spatiotemporal ORN responses [18,20]. So we develop a single-compartment biophysical OB model that accounts for differences in ORN input to investigate how they affect MC spiking responses. Specifically, we model ORN input as a time-varying inhomogeneous Poisson Process [27], where the input rate has slower increase and decay for retro than ortho [11,17], and the ORN input correlation is smaller for retro than ortho [11,17]. With these specifications, our biophysical OB network model is able to capture the salient ortho versus retro MC spiking response trends in our experimental data.
However, our biophysical OB model is too complex to directly analyze mathematically in order to address the neural encoding problem of characterizing how MCs convert ORN input to spike responses. We use a simple linear-nonlinear (LN) model framework to assess how our biophysical OB network transfers input statistics (from ORN) to outputs (MC spike statistics). We find that the linear filter component of the LN model, i.e., convolution with ORN inputs, consistently has larger absolute values with retro than with ortho input. Thus the OB network model is more sensitive to ORN fluctuations with retro-like inputs than with ortho. Finally, we use our models to examine which key attribute(s) of ORN inputs (temporal profile, amplitude, input correlation) are most significant for capturing our data. We find that temporal profile is the critical attribute for ortho versus retronasal stimulus response.
This work provides a framework for how to analyze the sources driving different OB spiking responses to different modes of olfaction, as well as important insights that have implications for how the brain codes odors.

Results
We performed in vivo multi-electrode array recordings of the OB in the mitral cell layer of anesthetized rats (see Materials and methods: Electrophysiological recordings) to capture odor evoked spiking activity of populations of putative MCs. This yielded a large number of cells (94) and simultaneously recorded pairs of cells (1435) with which to assess population average spiking statistics. The spike statistics are trial-averaged responses of a rat to a single odorant, Ethyl Butyrate (food odor). We focus on a food odor because they dominate retronasal smells, and a recent study showed that humans can more accurately detect food odors (vs. non-food odors) delivered retronasally [10]. In addition, an fMRI study showed different cortical activity [9] in humans for ortho versus retronasal stimulus, specifically with food odors.
The first and second order spike statistics are summarized in Fig 1, including the firing rate (peri-stimulus time histogram, PSTH, Fig 1A), the spike count variance (Fig 1B), the spike count covariance (Fig 1C), Fano Factor (variance divided by mean, Fig 1D), and Pearson's correlation ( Fig 1E). For each cell and simultaneously record pair of cells, we calculated the trialaveraged spike statistics with half-overlapping 100 ms time windows. The time window 100 ms is an intermediate value between shorter (membrane time constants, AMPA, GABA A , etc.) and longer time scales (NMDA, calcium, and other ionic currents) known to exist in the OB.
We find statistically significant differences between ortho and retro stimulation in almost all of the first and second order MC spike count statistics. At odor onset, orthonasal stimulation elicits larger firing rates with a faster rise than retronasal, after which retronasal firing is larger and remains elevated longer than with orthonasal. These trends are consistent with imaging studies of the glomeruli layer in OB in transgenic mice (see [17], their Fig 2) as well as EOG recordings of the superficial layers of the OB in rats (see [19], their Fig 7). More specifically, we find statistical significance (α = 0.01) between ortho-and retronasal firing rate after and for the duration of the odor stimulation. We also find that MC spike count covariance for ortho is significantly larger than retro for the entirety of the evoked state. MC spike count variance, however, is not significantly different. Note that we specifically tested whether the population averages (averaged over all cells for PSTH and spike variance, over all simultaneously recorded pairs for spike covariance) are significantly different between ortho and retro, via a two-sample t-test assuming unequal variances with the null hypothesis of equal population averages (see S1 Fig). Further, we calculated Cohen's d value to measure effect size [28], and find medium effect size for the majority of the evoked time period for all spike statistics considered except for variance, see S2 Fig. Hereafter, we mainly focus on understanding the differences in firing rate and spike count covariance because they directly impact common coding metrics (e.g. the Fisher information) in contrast to scaled measures of variability (Fano factor and Pearson's correlation) that are nonlinear functions of the entities that impact coding [12]. Moreover, Fano factor and correlation both depend on variance, which is not statistically different with ortho and retro (but see S1(D) and S1(E) Fig for completeness).

OB network model captures data trends
To better understand how differences in MC spiking with ortho and retro stimulation come about, we developed a single-compartment OB network model based on Li & Cleland's multicompartment model [29,30]. Our model is more computationally efficient than their larger multi-compartment models [29,30], requiring a fraction of the variables (tens of state variables instead of thousands). Importantly, our single-compartment model retains important biophysical features (Fig 2A).
In Fig 2A, we see in both models of MC (uncoupled) that the time to spiking decreases with increasing current values, and the number of spikes in a cluster increases with current values consistent with prior electrophysiological experiments [31][32][33]. The spacing between spike clusters and number of spikes in a cluster in our model (right) qualitatively match the Li & Cleland model (left). The sub-threshold oscillations are not as prominent as in Li & Cleland, but still apparent. In the uncoupled GC models, both ours and Li & Cleland's models exhibit a delay to the first spike with weak current step [34] (Fig 2A, bottom) and tonic firing without appreciable delay for higher current injections [35] (Fig 2A, middle and top). In the uncoupled PGC models, we do not observe repetitive firing in either models (Fig 2A, top and middle). Also, releasing from a hyperpolarizing current injection (bottom) can illicit spiking in both models, as observed by McQuiston & Katz [36]. Thus, we have a condensed OB model by using far less equations than Cleland's models while retaining many of the biophysical dynamics known to exist in these 3 important OB cell types.
Since our focus is on first and second order population-averaged spiking statistics, we use a minimal OB network model with 2 glomeruli (Fig 2B). Each glomerulus has a PGC, MC and GC; we also include a common GC that provides shared inhibition to both MCs because GCs are known to span multiple glomeruli and shape MC spike correlation [26,37,38]. Within the OB network, the PGC and GC cells provide presynaptic GABA A inhibition to the MCs they are coupled to, while MC provide both AMPA and NMDA excitation to PGC and GCs (see Materials and methods: Single-Compartment Biophysical Model for further details). The ORN synaptic inputs are an important component of this coupled OB network; they are driven by correlated inhomogeneous Poisson Process with increases in rate and correlation at odor onset. The specific time-varying input rate and correlation we use are shown in Fig 2Ci and 2Cii, respectively. The differences in ortho versus retro (Fig 2Ci and 2Cii) are based on prior studies of ORN input to the OB in response to both ortho and retro stimulation [11,17]. We fixed all model components and manually varied the ORN input rate λ O/R (t), see Materials and methods: Fitting biophysical network model to data for further details.
A comparison of first and second-order statistics between our OB model and in vivo data is shown in Fig 3. With the ORN activity specified in Fig 2C, our OB model is able to qualitatively capture trends seen in our data for firing rate and spike count covariance. Firing rates in Fig  3A show that both the model and data exhibit larger firing rates for ortho at odor onset followed by a sharper decline. After the initial increase in ortho firing rates, retro firing rates  Comparison of all first and second-order statistics of coupled OB network model to our data. A) Firing rate of ortho increases and decays faster than retro in both data and model. B) Variance of spike counts for ortho and retro shown for completeness, but recall that in experimental data that they are not statistically different. C) Covariance of spike counts is larger for ortho than retro in both data and model (left), but the magnitudes of data and model differ. Comparison of the ratio of retro covariance to ortho covariance in the evoked state (right) shows that the model captures the relative differences between ortho and retro-here μ (resp. σ) is the average (resp. standard deviation) ratio over 20 time bins in the evoked state. For A-C, top shaded error regions of data (retro in A, ortho in B,C) are cut-off to better compare model and data. Comparisons of the (D) Fano factor and (E) Pearson's spike count correlation shown for continue to increase, eventually becoming larger than ortho and remaining elevated longer, consistent with optical imaging experiments (see [17] their Fig 2). Although there is no significant difference in the spike count variance between ortho and retro in our experimental data, we show our data with model for completeness (Fig 3B).
Our OB model captures the trend that ortho spike count covariance is larger than retro after odor onset, Fig 3C (left). The OB model certainly does not capture the magnitude of the spike count covariance in the data; recall that covariance in our experimental data is the population average over all 1435 simultaneously recorded pairs with significant heterogeneity while our model is homogeneous. But the relative differences between retro and ortho (as measured by the ratio of retro to ortho covariance in the evoked state) are similar (Fig 3C, right). Thus our OB model captures the salient trends of the population-averaged spike count statistics. We also show comparisons of Fano Factor ( Fig 3D) and Pearson's correlation ( Fig 3E) for completeness. Consistent with our data, our OB model has larger Fano Factor and spike count correlation for ortho than with retro. In the evoked state, the OB model matches spike count correlation for both ortho and retro well. The larger ortho Fano factor in our data is captured in our model, but the difference is very modest.

How OB network transfers ORN input statistics
We next sought to better understand how our OB network model operates with different ORN inputs. In particular, we investigated whether the same OB network model transfers ortho and retro ORN inputs to MC spike outputs differently or not. We addressed this in a simple and transparent manner with a phenomenological LN model ( Fig 4A) to approximate the overall effects of the OB network on ORN inputs. LN-type models have often been used to circumvent the complexities in biophysical spiking models (see [39][40][41] and Discussion).
Description of the LN model. The LN model first applies a linear filter to the input, X(t), i.e., a convolution with a fixed temporal linear filter k, shifts the result by b, followed by a static non-linearity (exponential function) to produce an output Y(t), see Fig 4A: For our purposes, XðtÞ 2 fm S ðtÞ; s 2 S ðtÞ; Cov ðS 1 ðtÞ; S 2 ðtÞÞg are the statistics of ORN input synapses to the MCs, and Y(t) is an approximation to the statistics of MC spiking response: fPSTHðtÞ; s 2 R ðtÞ; Cov ðR 1 ðtÞ; R 2 ðtÞÞg. We calculate Y(t) (Eq (1)) by minimizing the L 2 -norm of the difference between Y (t) and the simulated MC spike statistic from the biophysical OB model. The LN model is applied separately to each statistic (further details to follow, see Eqs (also see Materials and methods: Linear-Nonlinear (LN) model: numerical details). This completeness despite both measures depending on spike count variance. D) The model has slightly larger Fano factor with ortho, consistent with the data. E) The model does qualitatively capture the spike count correlation for both ortho and retro, at least in the evoked state. That is, we consider different, separate LN models for each statistic, without any mixing effects (e.g., s 2 S ðtÞ does not directly affect PSTH(t)). Although output statistics generally depend on all input statistics [42][43][44], we emphasize that our ad-hoc approach here is meant to better understand how the OB model operates on each statistic and is not a principled alternative model.
By construction, in the biophysical OB model, both the inputs to each MC and the spike output of each MC have identical marginal statistics, so we are using the LN model to assess how univariate input statistics (mean/var) are mapped to univariate output statistics (mean/ var). The covariances depend on 2 variables (bivariate: (S 1 , S 2 ) for input and (R 1 , R 2 ) for output), but the LN model is only for assessing how covariance of inputs maps to covariance of outputs without directly modeling multiple random variables.
For the inputs to the LN model, we use an exact theoretical calculation for m S ðtÞ; s 2 S ðtÞ; Cov ðS 1 ðtÞ; S 2 ðtÞÞ rather than relying on Monte Carlo simulations. The ORN input synapses are driven by correlated time-varying inhomogenous Poisson processes yet we are still able to calculate the first and second order statistics of the ORN inputs in the limit of infinite number of realizations; detailed in Materials and methods: Calculating time-varying ORN input synapses, Eqs (13), (18) and (22). A comparison of Monte Carlo simulations of the actual ORN inputs used in our OB model results (Eqs (11) and (12)) to the theoretical calculation (Eqs (13), (18) and (22)) is shown in Fig 4B. We clearly see that the calculations (labeled 'Theory') matches all three ORN input statistics with smooth curves, properly accounting for both time-varying ORN input and time-varying input correlation. These calculations do not depend on any asymptotic assumptions; see S4 Applying LN models to biophysical OB model results. The LN model is able to fit well to the biophysical OB model output MC spike statistics for both ortho and retro stimuli, shown in Fig 5A. For this reason, we can assume that the LN model provides a decent  Table 1; they represent an absolute shift independent of the temporal dynamics. The b values are similar for ortho and retro for all statistics except spike count covariance. Although b is important for the resulting LN curves (dot-black in Fig 5A), it is not a part of the temporal processing of ORN inputs.

ORN input signatures for ortho/retro
Despite retro eliciting larger firing rates than ortho, the spike count covariance (as well as correlation and Fano factor) with retro stimulation is smaller than with ortho. It has long been known theoretically and experimentally that in uncoupled cells, the spike correlation increases with firing rate (at least with moderate to larger window sizes) [45], in contrast with our data. In coupled networks, the change in correlation with firing rate is complicated and depends on numerous factors [21][22][23][24]. Thus, the components of ORN inputs that result in these differences (higher firing and less covariance for retro than with ortho) in the same OB network are not obvious.
So we use our computational framework to uncover the important feature(s) of ORN input that: i) results in MC spike statistics consistent with our salient data trends, and ii) linearly filters ORN inputs with larger values with retro than with ortho input. Here we disregard the biological differences in ortho and retro ORN inputs to consider 3 core attributes of ORN inputs that influence how the biophysical OB model operates: • Temporal (faster increase and decay, or slower increase and decay; see • Input correlation (lower or higher, black and gray curves respectively, in Fig 6A, right) We consider a total of 8 different ORN input profiles consisting of various combinations of amplitude, input correlation, temporal profiles. The LN model fit to the OB model (i.e., MC spike statistics) for these 8 different ORN input profiles are all similar, well approximating how the OB coupled network transfers input statistics (see Fig 5A and S5 Fig). Fig 6B clearly shows that the slower increase and decay in input rate (redish/lighter) consistently results in linear filters k(t) with larger absolute values than with faster increase/decay (bluish/darker). The larger filter values holds with all 3 statistics, and with all variations of amplitude and input correlation. Thus, the OB network consistently has filters with larger absolute values when the input profile is slower (i.e., retronasal-like). The resulting LN model b values are listed in Table 2 for The parameter b for the LN model fits (Eq (1)) between orthonasal and retronasal are similar for a given statistic, except for spike count covariance. https://doi.org/10.1371/journal.pcbi.1009169.t001

PLOS COMPUTATIONAL BIOLOGY
Ortho-and retronasal olfactory bulb spiking differences reference, although these values represent an absolute scaling independent of the temporal dynamics. Fig 7 shows all 8 OB model results for each spike statistic. For all first and second order statistics, including scaled measures of variability, the most distinct attribute that distinguishes our model results is the temporal profile of input. Importantly, the temporal profile is the key attribute to best capture the differences in ortho and retro our experimental data (see Fig 3). The slow increase and decay in input rate consistently results in retro-like spiking trends while the fast increase and decay in input rate results in ortho-like spiking trends. Thus, our models show that the temporal profile is a signature of retro and ortho stimulation, and emphasizes the critical role of ORN inputs for shaping how the same OB network modulates ortho and retro stimuli.

Discussion
We investigated how odors processed via the orthonasal and retronasal routes result in different OB spike statistics, analyzing in detail how ORN inputs transfer to MC spike outputs. Motivated by our in vivo rat recordings that showed significant differences in first and second order spiking statistics of MC, we developed a realistic OB network model to investigate the dynamics of stimulus-evoked spike statistic modulation (higher firing and lower covariance/ correlation with retro than with ortho). Our OB model balances biophysical attributes [29,30] with computational efficiency. The OB model is able to capture salient trends in our data with both ortho and retro stimulation, and should be useful for future studies of OB. We successfully used the biophysical OB model, paired with a phenomenological LN model, to analyze how different ORN inputs lead to different dynamic transfer of input statistics. We also showed that the temporal profile of ORN inputs is a key determinant of ortho versus retro input via model matching our data. The output spike statistics are crucial because the OB relays odor information to higher cortical regions, and thus our work may have implications for odor processing with different modes of olfaction [9][10][11].
To the best of our knowledge, our experiments detail the differences in MC spiking with ortho and retro stimuli for the first time. However, the work of Scott et al. [19] is related; they used 4 electrodes to record OB spiking activity in the superficial layers of OB in rats. Their results are difficult to directly compare to ours as they focus on superficial OB in the epithelium rather than the mitral cell layer, but at least the trial-averaged firing rates in their data appear to be consistent with our data. Moreover, our multi-electrode array recordings enable us to consider trial-to-trial covariance of spiking.  Fig 6A). Different temporal profiles is key for both having different model spike statistics and for best matching qualitative differences in our data (see Fig 3). A) Firing rate in Hz (left) is slightly lower with low input rate amplitude, but no significant difference with different input correlations. B) Spike count variance, similar to firing rate, has only slightly lower values with low input rate amplitude. C) Spike count covariance is lower with lower input correlation for all of ortho evoked state (not surprisingly). However, retro (fast) with lower amplitude steadily increases above higher amplitude after about 1 s in the evoked state. D) Fano Factor model results only change modestly. E) Pearson's spike The key attribute(s) of ORN inputs that can result in different ortho and retro MC spike statistics consistent with our data are not obvious. Indeed, retro stimulation resulted in larger firing rates than ortho, and the spike count covariance (as well as correlation and Fano factor) with retro stimulation is smaller than with ortho, in contrast to uncoupled cells where correlation increases with firing rate [45]. Using various models, we were able to consider how three components of ORN inputs (temporal profile, amplitude, and input correlation) result in different OB dynamics with regards to transferring input statistics to outputs. Prior experiments [11,16,17] have shown these input components can differ with ortho and retro inputs. We found that the temporal profile (fast versus slow) plays a critical role for both capturing our data and for shaping the transfer of inputs to outputs, i.e., retro inputs consistently resulted in larger temporal filter values, so the OB network is more sensitive to fluctuations of retro input statistics than ortho. To capture the salient trends in our data, we find slower input rate (rise and decay) is a key signature of retronasal stimulation, while faster rise and decay [11,17,19] is similarly a key signature of orthonasal stimulus.
The temporal differences between ortho versus retro have previously been thought to play a role in distinguishing ortho/retro stimulation at the ORN [11,16,17,19,20], but whether this carried over to the OB and if this held at the level of spiking was unknown. Here we demonstrate the importance of different temporal input to OB for ortho versus retro.
We used an ad-hoc LN model framework because many of biological complexities are removed yet important features are retained. That is, neurons are known to linearly filter inputs, and spike generation is inherently nonlinear, i.e., finding linear filters of neurons is not new [41], and they are related to the spike-triggered average [46]. Thus, LN-type models have been used in many contexts, often to circumvent biophysical modeling, and most notably with generalized linear models [47,48] (also see [40]) where various filters (stimulus, post-spike) and model components are fit to data using maximum likelihood. Connecting the large gap between biophysical models and LN models is daunting, but see Ostojic and Brunel [39] who relate stochastic integrate-and-fire type models to LN. Our approach here is much simpler than the aforementioned works because we simply wanted to assess how a particular statistic (mean, variance or covariance) mapped via the OB network model in a simple and transparent manner. An enhanced data-driven approach to fitting an LN-type model that relies on experimental data of both ORN inputs and MC spike outputs with many trials might better reveal differences in how the OB operates with ortho versus retro. However, we currently do not know if such a dataset exists.
Here we list some limitations of our study. We only considered the MC response to a single food odor despite a large variety of food (and non-food) odors animals encounter. Different odors activate different olfactory receptors that could lead to qualitatively different population MC spiking activity than what we report here. Retronasal odors are predominately food odors, and studies have shown that humans can more accurately detect food odors (vs. non-food odors) delivered retronasally [10]. Frasnelli et al. [49] showed that food versus non-food odors illicit varying neural responses in humans when introduced ortho-versus retronasally. An fMRI BOLD study showed that cortical activity in humans differed when odors were introduced via the ortho or retro routes, specifically with food odors [9]. Thus our choice of a food odor is a natural first step for investigating retronasal MC responses. Also, we attributed the differences in ortho/retro MC responses solely to ORN inputs when in fact many regions synapse to OB [50]. For example, optogenetic studies [51,52] have shown that feedback from count correlation, similar to spike count covariance, is lower with lower input correlation and similarly for retro (fast), there is an increase with higher input correlation.
https://doi.org/10.1371/journal.pcbi.1009169.g007 olfactory cortex to OB is relatively strong and inhibition dominated. Whether this cortical feedback (or other external modulation) differs for ortho and retro stimulation is currently unknown. Moreover, odor-specific cortical feedback to OB [53] could alter the OB spike correlation, a factor that our modeling study did not account for. Finally, our data is from anaesthetized rats that enabled control of odor delivery and excluded confounding factors such as the breath cycle and sniff rate [54,55]. However, the MC spike activity in awake rodents can be quite different than in anesthetized [56], so whether our reported differences in ortho versus retro MC spiking hold in awake rodents is an open question. We hope our work here inspires more research into the differences between ortho versus retro olfaction, in particular in downstream olfactory circuits and with other experimental preparations.
With a combination of experiments and different scales of neural network modeling, we provide a basis for understanding how differences in OB spiking statistics arise with these 2 natural modes of olfaction. More generally, our model framework provides a road map for how to analyze attributes responsible for different OB spiking when driven by differences in ORN inputs.

Ethics statement
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.

Code availability
See https://github.com/michellecraft64/OB for MATLAB code implementing the single-compartment biophysical model, the equations for synaptic input statistics, and the linear-nonlinear (LN) model.

Single-compartment biophysical OB model
Models of all three cell types (MC, PGC, GC) are based on models developed by the Cleland Lab [29,30]. We consider two glomeruli each with a representative MC, PGC, GC (see Fig  2B). Each cell is a conductance-based model with intrinsic ionic currents. The voltage responses of all three cell types, measured in experiments and in a multi-compartment model [29,30], are generally captured in our single-compartmental model, see Fig 2A. Here we describe all of the pertinent model details thoroughly; for other extraneous details and implementation, please refer to provided code on GitHub.
Individual cell model.
The voltages of all model cells are governed by a Hodgkin-Huxley type current balance equation (Eq (6) above for the j th cell) consisting of voltage (V), membrane capacitance (C), applied current (I App ), ionic currents (I Ion ), synaptic currents (I Synapse ), and ORN inputs (I ORN ); see Tables 3 and 4 for units and numerical values, respectively. For our modeling purposes, the ionic currents and the ORN inputs are modified from [29,30] and described below.

Ionic currents.
The ionic currents are defined by Eq (7) above (for specific ion type i) and account for maximal conductance (g), activation variable (m) with exponent (p), inactivation variable (h) with exponent (q), time-varying voltage (V assumed to be isopotential), and reversal potential (E i ). All parameters and function for intrinsic ionic currents and their gating variables are the same as in [29,30] with the exception of maximal conductance. We chose to condense the model as defined in [29,30] by collapsing all compartments to a single-compartment, and we set the maximal conductance as the sum of all maximal conductance values (e.g., in PGC, I Na has maximal conductance g Na = 70 mS/cm 2 because [29] set g Na = 50 mS/cm 2 in the soma and g Na = 20 mS/cm 2 in the spine). All summed maximal conductance values used are listed for reference in Table 4. The calcium dynamics used to define the calcium-related ionic currents are the same as in [29,30]. Synaptic currents.
Eqs (8) and (9) are the equations for the synaptic variables, where all presynaptic GCs and PGCs provide GABA A inputs, and all presynaptic MCs provide both AMPA and NMDA Table 4. Parameter values for each cell type. Each of these values are the same as defined by [30] with the exception of maximal conductance values which are the sum of all cell compartments (soma, dendrite, axon, etc.) as defined by [30]. Additionally, any conductance value denoted by − implies that this ionic current is not included in the associated cell. inputs. B(V) in Eq (8) is the NMDA-specific magnesium block function (B(V) = 1 for all other synapses), and s(t) is the fraction of open synaptic channels. The channel opening rate constants (α and β) are normalized sigmoidal function of presynaptic membrane potential (F(V pre ) in Eq (9)), the same as in [29,30]. We also define the conductance parameter (g syn ) and reversal potentials (E syn ) as [29,30] have, with g GABA = 1.5 nS for GC!MC synapses, g GABA = 2 nS for PGC!MC synapses, g AMPA = 2 nS and g NMDA = 1 nS for both MC!PGC and MC!GC synapses; E syn = 0 mV for AMPA and NMDA currents, and E syn = −80 mV for GABA A currents. ORN input.

Resistance and
The ORN inputs for each cell consist of both excitatory and inhibitory inputs as specified in Eqs (10) and (11) where X 2 {E, I}. The reversal potential value (E X ) is much larger for excitatory inputs and smaller for inhibitory. The function Sðt þ k Þ ¼ Sðt À k Þ þ a X accounts for the random times (t k ) when S instantaneously increases by a X . The random times, t k , are governed by an inhomogeneous Poisson process with rate λ X (t). This aligns with experimental evidence that ORN spiking is Poisson-like in the spontaneous state [27]. Thus, we extend the notion that ORN spiking would be Poisson-like in the evoked state with increased rate λ X (t) varying in time. Finally, we set the synaptic rise and decay time constants (τ X ) to be 5.5 ms for PGCs and GCs, 10 ms for MCs, as in [29,30].
The ORN input rates can be pairwise correlated, which is achieved by the parameter c j,k 2 [0, 1], for cells j and k detailed by Eq (12) below: l j ðtÞ ¼l j ðtÞ À � lðtÞc j;k ðtÞ: ð12Þ wherel j ðtÞ andl k ðtÞ are the individually defined ORN input rates for cells j and k, and � lðtÞ ≔ minðl j ðtÞ;l k ðtÞÞ.

Fitting biophysical network model to data
The biophysical OB model described thus far was adopted directly from Li & Cleland, aside from our single-compartment simplification where we lumped all ionic currents into one compartment and used the sum of the (maximal) conductances from all compartments. Here we describe how the network model was tuned to capture the salient features of our experimental data. We did not systematically consider large volumes of parameter space due to the enormous computational resources required for 50,000 simulations of the model for each parameter set to accurately simulate the spike count statistics. After model parameters were set, the only manual tuning we did was to consider several Poisson input rates λ O/R (t) (see Eqs (10)-(12)) for the ORN input synapses (see S3 Fig)-even the ORN input correlations c j,k (t) that we arbitrarily chose were fixed throughout. Note that we did not further tune the intrinsic properties of the individual cells; the PGC, MC, and GC parameters are as stated above with behavior shown in Fig 2A. Specifying coupling strengths. We used the same equations for the synaptic variables as Cleland [30], but set the coupling strengths w (see Eq (8)) to: w M G = 3 (independent inhibition), w M Gc = 0.3 (common inhibition to MC), w G M = 1 (same for both AMPA, NMDA), w Gc M = 0.5 (inhibition to common GC), w P M = 1 and w M P = 2 (same for both AMPA, NMDA). These coupling strengths were chosen in part from results in Ly et. al [57] who used a related/simpler OB network model with the same 2 glomeruli architecture to find regions of parameter space that best fit orthonasal experimental data. Similar to Ly et. al [57] (see their Figs 2, 3, and 6) we set independent inhibition from GC to MC to be greater than excitation from MC to GC, and shared GC inhibition to MC to be relatively weak (i.e., w Gc M � w G M � w M G ). The coupling strengths were never tuned, they were fixed throughout.
Specifying ORN input. The ORN inputs (Eqs (10)-(12)) consists of a Poisson input rate, and input correlation between pairs of cells. We set the correlation (c j,k ) between the following cell pairs: MC and PGC pair within a glomerulus have c j,k = 0.3 because they receive inputs from the same ORN cells; the two MCs have correlated ORN input [58] (c j,k (t) time-varying as in Fig 2Cii); and between all 3 GCs because they are known to synchronize [30,59] (c j,k = 0.3 for all 3 different pairs of GCs). All other pairs of cells have no ORN input correlation. Note that input correlation for the 2 MCs increased with odor to mimic increased correlation of glomeruli activity. In Fig 2Cii, input correlation for the 2 MCs are constrained such that c R (t) < c O (t). This is based on prior imaging studies that show retronasal stimulation activates spatially smaller regions of glomeruli inputs than orthonasal, and that the activation regions from retro are subsets of ortho [11,17]. For specific algebraic formula of c R/O (t), please refer to code on GitHub.
We considered several different λ O/R (t), the inhomogeneous Poisson input rate of t k in Eq (11) (with constraints described below) and chose the ones that best matched the time-varying firing rates (Fig 7A). The ortho-vs. retronasal odor input rates, λ O/R (t), are constrained such that λ O (t) increases faster and more abruptly than λ R (t) with odor, and λ R (t) decays slower than λ O (t); this is all based on ORN imaging studies [11,17]. Inputs consist of both excitatory synapses (with rate λ O/R (t)) and inhibitory synapses (with rate 0.75λ O/R (t)) to capture other unmodeled inhibitory effects.
To first understand how MC firing rates depends on λ(t) without any consideration for ortho or retro, we used a simple alpha-function form in the evoked state: λ(t) = te −t/τ , surveying 6 different τ (see S3(A) Fig, left). The resulting MC firing rates (S3(A) Fig, right with 2,000 realizations) was informative for how to manually set the input values (spontaneous, peakedevoked, etc.). S3(B) Fig shows all of the λ O/R (t) we tried, notice that they all satisfy the constraints described above. Via trial and error with 2,000 realizations, we only looked at the resulting firing rates (PSTH), insuring the simulations matched the ortho data well. We were fortunate in fitting the retro firing rate data, trying only 2 λ R (t). The other spike statistics (e.g., covariance, Fano factor) were never accounted for in our consideration of different λ O/R (t), which is perhaps why the fit to the spike covariance data is so bad.

Calculating time-varying ORN input statistics of synapses
Here we describe a method to capture the effects of ORN input statistics of synapses to the biophysical OB model, in the limit of infinite realizations. These methods are very useful as inputs for the LN model, without which one would have to use averages from Monte Carlo simulations that contain deviations from finite size effects. Taking the expected value of Eq (11) results in an equation for the average of S(t), μ S (t): To derive the equation for variance s 2 S ðtÞ, we multiply Eq (11) by itself. We can equivalently rewrite Eq (11) as an integral: where DðtÞ ≔ X j dðt À t k Þ. So DðuÞDðvÞe À ðtÀ uÞ=t X e À ðtÀ vÞ=t X dudv ð15Þ Recall that E½DðuÞDðvÞ� ¼ lðvÞdðu À vÞ þ lðuÞlðvÞ, so we have: Equivalently, s 2 S ðtÞ satisfies the ODE: Similarly for S j (t)S k (t) correlated synapses, we have: By our model construction E½DðuÞDðvÞ� ¼ c j;k ðvÞ � lðvÞdðu À vÞ þ l j ðuÞl k ðvÞ, where � lðtÞ ≔ minðl j ðtÞ; l k ðtÞÞ, so we have: Cov(S j (t), S k (t)) equivalently satisfies the ODE: The calculations for the dynamic (time-varying) synapse statistics are important for capturing realistic statistics because a steady-state approximation can be very inaccurate, especially when the time-varying correlation and ORN spiking rate change quickly relative to the timescales (τ X ). The quasi-steady-state approximation is: CovðS j ðtÞ; S k ðtÞÞ � t X j t X k t X j þ t X k a X j a X k c j;k ðtÞ � lðtÞ  (13), (18) and (22)) and how inaccurate Eqs (23)-(25) can be.

Linear-Nonlinear (LN) model: Numerical details
We use the above described ODEs (Eqs (13), (18), and (22)) to simplify the calculations of ORN input statistics for use with the LN model framework. Previous work has implemented LN-type models as an alternative to biophysical spiking models with various conditions (see [39][40][41] and Discussion). Fig 4A illustrates a schematic of the LN model parameters, linear filter (k) and shift (b), that are used with the ORN input statistics (X 2 fm S ; s 2 S ; CovðS 1 ; S 2 Þg) in order to construct an approximation (Y) to the biophysical OB network model's output MC spike statistics (PSTH(t), s 2 R , Cov(R 1 , R 2 )). The LN model is summarized as: Where we define our function f as an exponential, and we can approximate the integral numerically as follows: Where n denotes the number of time points included in the linear filter, and j denotes the points in time of input statisticX of size Lt. Then, we can rewrite Eq (27) in matrix vector form Aw ¼ṽ where: A is the Toeplitz matrix of size (Lt − n + 1) × n of our input statistic (X) with an additional row of value one to account for shift;w is our linear filter (k) and shift (b); andṽ is our OB network MC spike statistic to which we fit our filter. Then, we solve forw using Least Squares approximation by QR decomposition. The linear filter (k) converges to 0 by construction, therefore we truncate the filter at −0.1 s and set k = 0 for the remaining time −1 � t < −0.1. Then, the LN output approximation (Y) of the MC spike statistic is calculated as follows: Where K denotes the convolutional matrix constructed from the truncated linear filter k.

Electrophysiological recordings
We decided to use recordings from a single rat, with recordings from 3 sessions. We took this conservative approach to control differences in nasal cavity structure that can vary across rats [60,61], which may shape differences in ortho versus retro activity [11,62]. See provided GitHub code for statistical summary of experimental data. 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). Data were collected from 11 adult male rats (240-427 g; Rattus Norvegicus, Sprague-Dawley outbred, Harlan Laboratories, TX, USA) housed in an environment of controlled humidity (60%) and temperature (23˚C) with 12h light-dark cycles. The experiments were performed in the light phase.
Surgical preparations. 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. Throughout surgery and electrophysiological recordings, core body temperature was maintained at 37˚C with a thermostatically controlled heating pad. To isolate the effects of olfactory stimulation from breath-related effects, we performed a double tracheotomy surgery as described previously [11]. A Teflon tube (OD 2.1 mm, upper tracheotomy tube) was inserted 10mm into the nasopharynx through the rostral end of the tracheal cut. Another Teflon tube (OD 2.3 mm, lower tracheotomy tube) was inserted into the caudal end of the tracheal cut to allow breathing, with the breath bypassing the nasal cavity. Both tubes were fixed and sealed to the tissues using surgical thread. Local anesthetic (2% Lidocaine) was applied at all pressure points and incisions. Subsequently, a craniotomy was performed on the dorsal surface of the skull over the right olfactory bulb (2 mm × 2 mm, centered 8.5 mm rostral to bregma and 1.5 mm lateral from midline).
Olfactory stimulation. 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. Odorized air was delivered for 1 s in duration at 1 minute intervals, with a flow rate of 250 ml/min and 1% of saturated vapor. The odorant was Ethyl Butyrate (EB). We note that the full experimental data set included additional odors, but here we consider only EB.
Electrophysiology. A 32-channel microelectrode array (MEA, A4x2tet, NeuroNexus, MI, USA) was inserted 400 μm deep from dorsal surface of OB targeting tufted and mitral cell populations. The MEA probe consisted of 4 shanks (diameter: 15 μm, inter-shank spacing: 200 μm), each with eight iridium recording sites arranged in two tetrode groups near the shank tip (inter-tetrode spacing: 150 μm, within tetrode spacing 25 μm). Simultaneous with the OB recordings, we recorded from a second MEA placed in anterior piriform cortex. Voltage was 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 (Cereplex + Cerebus, Blackrock Microsystems, UT, USA). Recordings were bandpass filtered between 300 and 3000Hz and semiautomatic spike sorting was performed using Klustakwik software, which is well suited to the type of electrode arrays used here [63].

S1 Fig. Statistically significant different spike count statistics.
We performed two-sample ttests assuming unequal variances for each point in time to assess whether the spike count statistics are significantly different with ortho and retro stimulation. We find statistical significance (α = 0.01) between ortho and retronasal firing rate (A) after and for the duration of odor stimulation (0.3 � t � 1 s with 100 ms time windows and 0.5 � t � 1.1 s with 200 ms time windows) as well as spike count covariance (C) for the entirety of the evoked state (0 � t � 2 s excluding t = 0 s with 200 ms time window). Spike count variance (B) is not found to have any statistical significant differences between ortho and retro. For completeness, significance of Fano Factor (D) and Pearson's correlation (E) are also significantly different for ortho and retro in the evoked state (0 < t � 2 s for Fano Factor and 0 � t � 2 s excluding t = 0 s with 100 ms time window for correlation). We calculated Cohen's d value for the nondirectional (two-tailed) case to measure effect size index for t-tests of means (see S1 Fig) in standard units. We find small (t = 0.3, 0.7 � t � 0.9 s with 100 ms; 0.6 < t �1 s with 200 ms) and medium (0.3 < t < 0.7 s, t = 1 s with 100 ms; 0.4 � t �0.6 s with 200 ms) effect size of statistical significance between ortho and retronasal firing rate (A) as well as small (0 � t � 2 s excluding t = 0 s with 200 ms time windows) effect size of spike count covariance (C). Spike count variance (B) does not have a measure of effect size since it is not found to have any statistical significant differences between ortho and retro. For completeness, effect size of Fano Factor (D) and Pearson's correlation (E) are also found to be small (0 < t �0.2 s [ 0.4 < t �1 [ 1.4 < t �2 s for Fano Factor, and 0 � t < 0.2 [ 0.5 < t < 0.8 excluding t = 0 s with 100 ms time windows for correlation) and medium (0.2 < t �0.4 [ 1 < t �1.4 s for Fano Factor and 0.2� t � 0.5 [ 0.8 � t � 2 s for correlation). (TIF) S3 Fig. Details of various ORN input rates we surveyed (achieved via trial and error) λ(t). A) Left: initial set of ORN inputs λ(t) (with evoked λ(t) = (t + 1)e −(t+1)/τ ) we surveyed to better understand the MC firing rate (right), calculated with 2,000 realizations. B) Fitting the ortho firing rate well enough required considering many λ O (t), and we even shifted the spontaneous input rate up slightly at some point. However, the only 2 retro inputs we tried (pink and red) were relatively accurate. (TIF) S4 Fig. ORN synapses statistic calculation is robust and accurate. Our theory for the ORN synaptic input statistics (Eqs (13), (18) and (22)) is accurate for time-varying inhomogeneous Poisson process rates and time-varying input correlation. A) Ortho-like input (fast rise and decay of Poisson rate) with same amplitude as retro (high), but with low input correlation used to capture data. Notice how the theory captures the fine structure of the covariance (double-hump). B) Retro-like input (slow rise and decay) with same amplitude as ortho (high), but with high input correlation used to capture data. C, D) Demonstrating accuracy of dynamic theory with much slower (unrealistic) time-scales: τ 1 = 50 ms and τ 2 = 100 ms and faster relative change in Poisson rate (all with low input correlation). Showing the quasi-steady-state approximation (Eqs (23)