Analyzing the differences in olfactory bulb mitral cell spiking with ortho- and retronasal stimulation

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 processing of retronasal odors via coordinated spiking of neurons in the olfactory bulb (OB) is largely unknown. 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). To better understand these differences, 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 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 temporal 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. Author summary Olfaction is a key sense for many cognitive and behavioral tasks, and is particularly unique because odors can naturally enter the nasal cavity from the front or rear, i.e., ortho- and retro-nasal, respectively. Yet little is known about the differences in coordinated spiking in the olfactory bulb with ortho versus retro stimulation, let alone how these different modes of olfaction may alter coding of odors. We simultaneously record many cells in rat olfactory bulb to assess the differences in spiking statistics, and develop a biophysical olfactory bulb network model to study the reasons for these differences. Using theoretical and computational methods, we find that the olfactory bulb transfers input statistics differently for retro stimulation relative to ortho stimulation. Furthermore, our models show that the temporal profile of inputs is crucial for capturing our data and is thus a distinguishing feature between ortho and retro stimulation, even at the olfactory bulb. Understanding the spiking dynamics of the olfactory bulb with both ortho and retro stimulation is a key step for ultimately understanding how the brain codes odors with different modes of olfaction.

Olfactory processing naturally occurs in two distinct modes: orthonasal (ortho) where 2 odors enter the front of the nasal cavity and retronasal (retro) where odors enter the 3 rear. Despite the importance of retronasal olfaction, i.e., retro naturally occurs during 4 feeding, it is relatively understudied. Specifically, it is unknown whether odors in higher 5 brain regions are processed differently depending on the mode of olfaction (ortho versus 6 retro). The olfactory bulb (OB) is the main area where odor information is processed 7 and subsequently transferred to cortex via mitral cell (MC) (and tufted cell) spiking. 8 The differences in MC spiking with ortho and retro stimulation have implications for 9 odor processing, but any such differences are largely unknown. 10 Studies have shown that olfactory receptor neuron (ORN) activity, which are 11 presynaptic to the OB, differ for ortho versus retro stimulation. This has been shown in 12 various ways, including with fMRI [1], calcium imaging [2], and optical imaging in 13 transgenic mice [3]. These and other prior studies [4][5][6] suggest that ORN inputs likely 14 lead to any observed differences in OB activity. However, the implications of these 15 differences in ORN for MC spiking have yet to be explored. 16 We perform in vivo recordings of rat OB mitral cells using multi-electrode arrays 17 with a food odor stimulus, delivered by both modes of stimulation, to determine 18 whether there are differences in MC population spiking. We find significant differences 19 in odor-evoked MC spiking with ortho versus retro stimulation in both firing rate 20 (larger with retro) and spike count covariance (larger with ortho). Dissecting how 21 components of ORN inputs alter OB spiking is difficult experimentally due to the 22 complexity of both the recurrent circuitry in the OB [7,8] and resulting spatiotemporal 23 ORN responses [4,6]. So we developed a single-compartment biophysical OB model to 24 investigate how ORN input statistics affect the first and second order MC spiking 25 statistics. Specifically, we model ORN input as a time-varying inhomogeneous Poisson 26 Process [9], where the input rate has slower increase and decay for retro than 27 ortho [2,3], and the ORN input correlation is smaller for retro than ortho [2,3]. With 28 these specifications, our biophysical OB network model is able to capture the ortho 29 versus retro MC spiking response trends in our experimental data. 30 Understanding how retro stimulation can elicit both larger firing rates and smaller 31 co-variablity than ortho is generally difficult in recurrent networks because of the 32 numerous attributes that shape spike statistics [10][11][12][13]. Since our biophysical OB model 33 is too complex to directly analyze, we use a linear-nonlinear (LN) model framework to 34 analyze how our realistic OB network transfers input statistics (from ORN) to outputs 35 (MC spike statistics). We find that with retro inputs, the OB network effectively filters 36 input statistics (in time) with larger absolute values than with ortho inputs. Thus the spiking responses to different modes of olfaction, as well as important insights that have 44 implications for how the brain codes odors. 45 Results 46 We performed in vivo multi-electrode array recordings of the OB in the mitral cell layer 47 of anesthetized rats (see Materials and methods: Electrophysiological 48 recordings) to capture odor evoked spiking activity of populations of putative MCs. 49 This yielded a large number of cells (94) and simultaneously recorded pairs of cells 50 (1435) with which to assess population average spiking statistics. The first and second 51 order spike statistics are summarized in Figure 1, including the firing rate (peri-stimulus 52 time histogram, PSTH, Fig 1A), the spike count variance (Fig 1B), the spike count 53 covariance (Fig 1C), Fano Factor (variance divided by mean, Fig 1D), and Pearson's 54 correlation ( Fig 1E). Spike count statistics were calculated with half-overlapping 100 ms 55 time windows. The time window 100 ms is an intermediate value between shorter 56 (membrane time constants, AMPA, GABA A , etc.) and longer time scales (NMDA,57 calcium, and other ionic currents) known to exist in the OB. 58 We find statistically significant differences between ortho and retro stimulation in 59 almost all of the first and second order MC spike count statistics. At odor onset, 60 orthonasal stimulation elicited larger firing rates with a faster rise than retronasal, after 61 which retronasal firing is larger and remains elevated longer than with orthonasal. These 62 trends are consistent with imaging studies of the glomeruli layer in OB in transgenic 63 mice (see [3], their figure 2) as well as EOG recordings of the superficial layers of the 64 OB in rats (see [5], their figure 7). More specifically, we find statistical significance 65 (α = 0.01) between ortho-and retronasal firing rate after and for the duration of the 66 odor stimulation. We also find that MC spike count covariance for orthonasal stimulus 67 was significantly larger than retronasal for the entirety of the evoked state. MC spike 68 count variance, however, was not significantly different between ortho and retronasal 69 stimulus. For detailed plots of significance in time of the first and second order spike 70 statistics, see Fig S1 in S1. Supplementary Material. Hereafter, we mainly focus on 71 understanding the differences in firing rate and spike count covariance because they 72 directly impact common coding metrics (e.g. the Fisher information) in contrast to 73 scaled measures of variability (Fano factor and Pearson's correlation) which do not 74 directly impact common coding metrics [14]. Moreover, Fano factor and correlation 75 both depend on variance, which is not statistically different with ortho and retro (but 76 see Fig S1D and  Cleland's multicompartment model [15,16]. Our model is more computationally efficient 81 than their larger multi-compartment models [15,16], requiring a fraction of the variables 82 (tens of state variables instead of thousands). Importantly, our single-compartment 83 model retains important biophysical features (Fig 2A).

84
In Fig 2A, [22]. Thus, we have condensed the OB model by using 96 far less equations than Cleland's models while retaining many of the biophysical 97 dynamics known to exist in these 3 important OB cell types.

98
Since our focus is on first and second order population-averaged spiking statistics, 99 we use a minimal OB network model with 2 glomeruli (Fig. 2B). Each glomerulus has a 100 PGC, MC and GC cell; we also include a common GC that provides shared inhibition 101 to both MCs because GCs are known to span multiple glomeruli and shape MC spike 102 correlation [8,23,24]. Within the OB network, the PGC and GC cells provide able to qualitatively capture trends seen in our data. Firing rates in Fig 3A show that 115 both the model and data exhibit larger firing rates for ortho at odor onset followed by a 116 sharper decline. After the initial increase in ortho firing rates, retro firing rates continue 117 to increase, eventually becoming larger than ortho and remaining elevated longer, 118 consistent with optical imaging experiments (see [3] their Fig 2). Although there is no 119 significant differences in the spike count variance between ortho and retro in our 120 experimental data, we show our data with model for completeness ( Fig 3B).

121
Our OB model captures the trend that ortho spike count covariance is larger than 122 retro after odor onset, Fig 3C (left). The OB model certainly does not capture the 123 magnitude of the spike count covariance in the data; recall that covariance in our 124 experimental data is the population average over all 1435 simultaneously recorded pairs 125 with significant heterogeneity while our model is homogeneous. But the relative 126 differences between retro and ortho (as measured by the ratio of retro to ortho 127 covariance in the evoked state) are similar (Fig 3C, right). Thus our OB model captures 128 the trends of the population-averaged spike count statistics. We also show comparisons 129 of Fano Factor ( Fig 3D) and Pearson's correlation ( Fig 3E) for completeness. Consistent 130 with our data, our OB model has larger Fano Factor and spike count correlation for 131 ortho than with retro. In the evoked state, the OB model matches spike count 132 correlation for both ortho and retro well. The larger ortho Fano factor in our data is 133 captured in our model, but the difference is very modest.

134
How OB network transfers ORN input statistics 135 We next sought to better understand how our OB network model operates with 136 different ORN inputs. In particular, we investigated whether the same OB network 137 model transferred ortho and retro ORN inputs to MC spike outputs differently or not. 138 We addressed this in a simple and transparent manner, using a phenomenological LN 139 model ( Fig 4A) to approximate the overall effects of the OB network on ORN inputs.

140
LN-type models have often been to circumvent the complexities in biophysical spiking 141 models (see [25][26][27] and Discussion). 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. Right: comparison of the ratio of retro covariance to ortho covariance in the evoked state, t ≥ 0 s, shows that the model captures the relative differences between ortho and retro. Comparisons of the (D) Fano factor and (E) Pearson's spike count correlation shown for 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.
produce an output R(t): MC spiking response. For simplicity, we only consider: The LN model only accounts for specific input statistics transferred to their 149 corresponding output statistics, without any mixing effects (e.g., σ 2 S (t) does not directly 150 affect PSTH(t)). Although output statistics generally depend on all input 151 statistics [28][29][30], we emphasize that our ad-hoc approach here is meant to better 152 understand how the OB model operates and is not a principled alternative model.

153
For the inputs to the LN model, we use an exact theoretical calculation for  Fig 5A), it is not a part of the temporal processing of ORN inputs. 177 Table 1. Parameter b for LN model fits to MC spiking statistics in Figure 5.

ORN input signatures for ortho/retro
178 Despite retro eliciting larger firing rates than ortho, the spike count covariance (as well 179 as correlation and Fano factor) with retro stimulation is smaller than with ortho. It has 180 June 1, 2021 9/23 long been known theoretically and experimentally that in uncoupled cells, the spike 181 correlation increases with firing rate (at least with moderate to larger window sizes) [31], 182 in contrast with our data. In coupled (recurrent) networks, the change in correlation 183 with firing rate is complicated and depends on numerous factors [10][11][12][13]. Thus, the 184 components of ORN inputs that result in these differences (higher firing and less 185 covariance for retro than with ortho) in the same OB network are not obvious.

186
So we use our computational framework to uncover the important feature(s) of ORN 187 input that: i) results in MC spike statistics consistent with our data trends, and ii) 188 filters with larger absolute values with retro than with ortho input. Here we disregard 189 the biological differences in ortho and retro ORN inputs to consider 3 core attributes of 190 ORN inputs that influence how the OB model operates:  Table 2 (1)) of the various parameters for temporal profile, amplitude, and input correlation. Amplitude and ORN input correlation profiles as defined for Figs 3-5 and associated values previously listed in Table 1 are noted in bold. Figure 7 shows all 8 OB model results for each spike statistic. For all first and second 208 order statistics, including scaled measures of variability, the most distinct attribute that 209 distinguishes our model results is the temporal profile of input. Importantly, the 210 temporal profile is the key attribute to best capture the differences in ortho and retro 211 our experimental data (see Fig. 3). The slow increase and decay in input rate     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 differences 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 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.

217
We have investigated how odors processed via the orthonasal and retronasal routes 218 result in different OB spike statistics, analyzing in detail how ORN inputs transfer to 219 MC spike outputs. Motivated by our in vivo rat recordings that show significant 220 differences in first and second order spiking statistics of MC, we developed a realistic 221 OB network model to investigate the dynamics of stimulus-evoked spike statistic 222 modulation (higher firing and lower covariance/correlation with retro than with ortho). 223 Our OB model balances biophysical attributes [15,16] with computational efficiency.

224
The OB model is able to capture trends in our data with both ortho and retro 225 stimulation, and should be useful for future studies of OB. We successfully used the 226 biophysical OB model, paired with a phenomenological LN model, to analyze how 227 different ORN inputs lead to different dynamic transfer of input statistics. We also 228 showed that the temporal profile of ORN inputs is a key determinant of ortho versus 229 retro input via model matching our data. The output spike statistics are crucial because 230 the OB relays odor information to higher cortical regions, and thus our work may have 231 implications for odor processing with different modes of olfaction.

232
To the best of our knowledge, our experiments detail the differences in MC spiking 233 with ortho and retro stimuli for the first time. However, the work of Scott et al. [5] is 234 related; they used 4 electrodes to record OB spiking activity in the superficial layers of 235 OB in rats. Their results are difficult to directly compare to ours as they focus on 236 superficial OB in the epithelium rather than the mitral cell layer, but at least the 237 trial-averaged firing rates in their data appear to be consistent with our data. Moreover, 238 our multi-electrode array recordings enable us to consider trial-to-trial covariance of 239 spiking.

240
The key attribute(s) of ORN inputs that can result in different ortho and retro 241 statistics consistent with our data are not obvious. Indeed, retro stimulation resulted in 242 larger firing rates than ortho, and the spike count covariance (as well as correlation and 243 Fano factor) with retro stimulation is smaller than with ortho, in contrast to uncoupled 244 cells where correlation increases with firing rate [31]. Using various models, we were 245 able to consider how three components of ORN inputs (temporal profile, amplitude, and 246 input correlation) result in different OB dynamics with regards to transferring input 247 statistics to outputs. Prior experiments [1][2][3] have shown these input components can 248 differ with ortho and retro inputs. We found that the temporal profile (fast versus slow) 249 plays a critical role for both capturing our data and for shaping the transfer of inputs to 250 outputs, i.e., retro inputs consistently resulted in larger temporal filter values, so the 251 OB network is more sensitive to fluctuations of retro input statistics than ortho. The 252 slower input rate (rise and decay) is a key signature of retronasal stimulation [2,3,5] to 253 capture the trends in our data with retro stimulation, while faster rise and decay is 254 similarly a key signature of orthonasal stimulus.

255
The temporal differences between ortho versus retro have previously been thought to 256 play a role in distinguishing ortho/retro stimulation at the ORN [1-3, 5, 6], but whether 257 this carried over to the OB and if this held at the level of spiking was unknown. Here 258 we demonstrate the importance of different temporal input to OB for ortho versus retro. 259 We used an ad-hoc LN model framework because many of biological complexities are 260 removed yet important features are retained. That is, neurons are known to linearly 261 filter inputs, i.e., finding linear filters of neurons is not new [27], and they are related to 262 the spike-triggered average [32], and spike generation is inherently nonlinear. Thus, 263 LN-type models have been used in many contexts, and often to circumvent biophysical 264 modeling, most notably with generalized linear models [33,34] (also see [26]) where 265 various filters and model components are fit to data using maximum likelihood.

266
Connecting the large gap between biophysical models and LN models is daunting, but 267 see Ostojic and Brunel [25] who relate stochastic integrate-and-fire type models to LN. 268 Our approach here is much simpler than the aforementioned works because we simply 269 wanted to assess how a particular statistic (mean, variance or covariance) mapped via 270 the OB network in a simple and transparent manner.

271
With a combination of experiments and different scales of neural network modeling, 272 we provide a basis for understanding how differences in OB spiking statistics arise with 273 these 2 natural modes of olfaction. More generally, our model framework provides a 274 road map for how to analyze attributes responsible for different OB spiking when driven 275 by differences in ORN inputs.  Cleland Lab [15,16]. We consider two glomeruli each with a representative MC, PGC, 283 GC (see Fig 2B). Each cell is a conductance-based model with intrinsic ionic currents. 284 The voltage responses of all three cell types, measured in experiments and in a 285 multi-compartment model [15,16], are generally captured in our single-compartmental 286 model, see Fig 2A. Here we describe all of the pertinent model details thoroughly; for 287 other extraneous details and implementation, please refer to provided code on GitHub. 288

Individual cell model
The voltages of all model cells are governed by a Hodgkin-Huxley type current  Table 3 and 4 for units and numerical values,   293 respectively. For our modeling purposes, the ionic currents and the ORN inputs are 294 modified from [15,16] and described below.

295
Ionic currents The ionic currents are defined by Eq (6) above (for specific ion type i) and account 297 for maximal conductance (g), activation variable (m) with exponent (p), inactivation 298 variable (h) with exponent (q), time-varying voltage (V assumed to be isopotential), 299 and reversal potential (E i ). All parameters and function for intrinsic ionic currents and 300 their gating variables are the same as in [15,16] with the exception of maximal 301 conductance. We chose to condense the model as defined in [15,16] by collapsing all 302 compartments to a single-compartment, and we set the maximal conductance as the 303 sum of all maximal conductance values (e.g., in PGC, I Na has maximal conductance 304 g Na = 70 mS/cm 2 because [15] set g Na = 50 mS/cm 2 in the soma and g Na = 20 mS/cm 2 305 in the spine). All summed maximal conductance values used are listed for reference in 306  Table 4. The calcium dynamics used to define the calcium-related ionic currents are the 307 same as in [15,16].
Eqs (7) and (8)   The channel opening rate constants (α and β) are normalized sigmoidal function of 314 presynaptic membrane potential (F (V pre ) in Eq (8)), the same as in [15,16]. We also 315 define the conductance parameter (g syn ) and reversal potentials (E syn ) as [15,16] have, 316 with g GABA = 1.5 nS for GC→MC synapses, g GABA = 2 nS for PGC→MC synapses, 317 Table 4. Parameter values for each cell type. Each of these values are the same as defined by [16] with the exception of maximal conductance values which are the sum of all cell compartments (soma, dendrite, axon, etc.) as defined by [16]. Additionally, any conductance value denoted by − implies that this ionic current is not included in the associated cell.  In order to capture different effects of coupling, we include coupling strengths w. The 320 synaptic coupling strengths are fixed and set to: w M ←G = 3 (independent inhibition), 321 w M ←Gc = 0.3 (common inhibition to MC), w G←M = 1 (same for both AMPA, NMDA), 322 w Gc←M = 0.5 (inhibition to common GC), w P ←M = 1 and w M ←P = 2 (same for both 323 AMPA, NMDA). These coupling strengths were established based on preliminary results 324 by Ly et. al [35] who use a related biophysical OB network model to evaluate regions of 325 parameter spaces that provide best model fits to our experimental data. Similar to as 326 seen in Ly et. al [35] (see their Figs 2 and 3) we define independent inhibition defined 327 for GC to be greater than excitation from MC and common GC inhibition to be less 328 than both GC inhibition and MC excitation (i.e., w Gc←M ≤ w G←M ≤ w M ←G ).

ORN input
The ORN inputs for each cell consist of both excitatory and inhibitory inputs as 331 specified in Eqs (9) and (10)  Poisson-like in the spontaneous state [9]. Thus, we extend the notion that ORN spiking 337 would be Poisson-like in the evoked state with increased rate λ X (t) varying in time.

338
Finally, we set the synaptic rise and decay time constants (τ X ) to be 5.5 ms for PGCs 339 and GCs, 10 ms for MCs, as in [15,16].  The ORN input rates can be pairwise correlated, which is achieved by the parameter 351 c j,k ∈ [0, 1], for cells j and k detailed by Eq (11) below: 352 λ j (t) =λ j (t) −λ(t)c j,k (t). (11) whereλ j (t) andλ k (t) are the individually defined ORN input rates for cells j and k, 353 andλ(t) := min(λ j (t),λ k (t)). 354 We set the correlation (c j,k ) between the following cell pairs: MC and PGC pair 355 within a glomerulus that have inputs from the same ORN cells (c j,k = 0.3); two MCs as 356 they are known to have correlated ORN input [36] (c j,k (t) time-varying as in Fig 2Cii); 357 and between all 3 GCs because they are known to synchronize with common 358 input [16,37] (c j,k = 0.3). All other pairs of cells have no background input correlation. 359 Specifically, input correlation for the 2 MCs (c j,k (t)) varied in time to mimic increased 360 correlation of glomeruli activity with odor onset and additionally account for ortho vs. 361 retronasal odor input (c O/R (t)). As seen in Fig 2Cii, To derive the equation for variance σ 2 S (t), we multiply Eq (10) by itself. We can 373 equivalently rewrite Eq (10) as an integral: where Recall Equivalently, σ 2 S (t) satisfies the ODE: Similarly for S j (t)S k (t) correlated synapses, we have: where 379 λ(t) := min (λ j (t), λ 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 382 capturing realistic statistics because a steady-state approximation can be very 383 inaccurate, especially when the time-varying correlation and ORN spiking rate change 384 quickly relative to the time-scales (τ X ). The quasi-steady-state approximation is: Cov(S j (t), S k (t)) ≈ τ Xj τ X k τ Xj + τ X k a Xj a X k c j,k (t)λ(t) Figure S2 in S1. implemented LN-type models as an alternative to biophysical spiking models with 392 various conditions (see [25][26][27] and Discussion). Fig 4A illustrates Where we define our function f as an exponential, and we can approximate the integral 397 numerically as follows: Where n denotes the number of time points included in the linear filter, and j 399 denotes the points in time of input statistic S of size Lt. Then, we can rewrite Eq (26) 400 in matrix vector form A x = y where: A is the Toeplitz matrix of size (Lt − n + 1) × n 401 of our input statistic ( S) with an additional row of value one to account for shift; x is 402 our linear filter ( k) and shift ( b); and y is our OB network firing rate statistic to which 403 we fit our filter. Then, we solve for x using Least Squares approximation by QR Where K denotes the convolutional matrix constructed from the truncated linear filter 408 k.

409
Electrophysiological recordings 410 We decided to use recordings from a single rat, with recordings from 3 sessions. We took 411 this conservative approach to control differences in nasal cavity structure that can vary 412 across rats [38,39], which may shape differences in ortho versus retro activity [2,40]. See 413 provided GitHub code for statistical summary of experimental data.

414
All procedures were carried out in accordance with the recommendations in the intraperitoneal injection (ip)). Dexamethasone (2 mg/kg bw, ip) and atropine sulphate 424 (0.4 mg/kg bw, ip) were administered before performing surgical procedures.

425
Throughout surgery and electrophysiological recordings, core body temperature was 426 maintained at 37 • C with a thermostatically controlled heating pad. To isolate the 427 effects of olfactory stimulation from breath-related effects, we performed a double 428 tracheotomy surgery as described previously [2]. A Teflon tube (OD 2.1 mm, upper 429 tracheotomy tube) was inserted 10mm into the nasopharynx through the rostral end of 430 the tracheal cut. Another Teflon tube (OD 2.3 mm, lower tracheotomy tube) was 431 inserted into the caudal end of the tracheal cut to allow breathing, with the breath 432 bypassing the nasal cavity. Both tubes were fixed and sealed to the tissues using 433 surgical thread. Local anesthetic (2% Lidocaine) was applied at all pressure points and 434 incisions. Subsequently, a craniotomy was performed on the dorsal surface of the skull 435 over the right olfactory bulb (2 mm × 2 mm, centered 8.5 mm rostral to bregma and 436 1.5 mm lateral from midline).

437
Olfactory Stimulation. A Teflon tube was inserted into the right nostril and the 438 left nostril was sealed by suturing. The upper tracheotomy tube inserted into the 439 nasopharynx was used to deliver odor stimuli retronasally. Odorized air was delivered 440 for 1 s in duration at 1 minute intervals. The odorant was Ethyl Butyrate (EB, 441 saturated vapor). We note that the full experimental data set included additional odors, 442 but here we consider only EB. ground pellet placed in the saline-soaked gel foams covering the exposed brain surface 451 around the inserted MEAs. Voltages were digitized with 30 kHz sample rate (Cereplex 452 + Cerebus, Blackrock Microsystems, UT, USA). Recordings were band-pass filtered 453 between 300 and 3000Hz and semiautomatic spike sorting was performed using 454 Klustakwik software, which is well suited to the type of electrode arrays used here [41]. 455 Supporting Information 456 S1. Supplementary Material. Other parts of the research that support 457 the main results but are separated for a more streamlined exposition.