Odor modality is transmitted to cortical brain regions from the olfactory bulb

Odor perception is the impetus for important animal behaviors with two predominate modes of processing: odors pass through the front of the nose (orthonasal) while inhaling and snifﬁng, or through the rear (retronasal) during exhalation and while eating. Despite the importance of ol-faction for an animal’s well-being and that ortho and retro naturally occur, it is unknown how the modality (ortho versus retro) is even transmitted to cortical brain regions, which could signiﬁ-cantly affect how odors are processed and perceived. Using multi-electrode array recordings in tra-cheotomized anesthetized rats, which decouples ortho-retro modality from breathing, we show that mitral cells in rat olfactory bulb can reliably and directly transmit ortho versus retronasal modality with ethyl butyrate, a common food odor. Drug manipulations affecting synaptic inhibition via GABA A lead to worse decoding of ortho versus retro, independent of whether overall inhibition increases or decreases, suggesting that the olfactory bulb circuit may naturally favor encoding this important aspect of odors. Detailed data analysis paired with a ﬁring rate model that captures population trends in spiking statistics shows how this circuit can encode odor modality. We have not only demonstrated that ortho/retro information is encoded to downstream brain regions, but also use modeling to demonstrate a plausible mechanism for this encoding: due to synaptic adaptation, it is the slower time course of the retronasal stimulation that causes retronasal responses to be stronger and less sensitive to inhibitory drug manipulations than orthonasal responses


Introduction
Olfaction is driven by odor molecules entering the nasal cavity, inducing a cascade of action potentials in the nervous system that transmit and process odor information.There are two routes by which odor molecules can enter the nasal cavity: through the nostrils during inhalation and sniffing (ortho), or from the throat during exhalation and while eating (retro).Orthonasal stimulation is by far the most commonly studied modality, despite the fact that retonasal delivery naturally occurs during eating.Prior imaging studies have shown differences in the activation of the regions of the nasal epithelium with ortho and retro stimulation (Sanganahalli et al., 2020;Gautam and Verhagen, 2012b;Furudono et al., 2013).Moreover, humans are able to discriminate whether food odors are delivered ortho or retronasally without being told which modality and without active sniffing (Frasnelli et al., 2008).Imaging studies have shown that at least with food odors, ortho/retro are two distinct modalities, rather than 2 different routes to the same modality, independent of odor intensity (Bender et al., 2009).Importantly, there is evidence that cortical processing depends on whether odors are delivered ortho-versus retro-nasally (Blankenship et al., 2019;Gautam and Verhagen, 2012a), and humans report perceiving the same smells differently depending on ortho or retro modality (Hannum et al., 2021;Small et al., 2005).Together these studies convincingly show that ortho versus retro information is transmitted to cortical brain regions, but the details of how this occurs and the core network mechanisms that facilitate this encoding are largely unknown.
There are several ways that ortho/retro information could be made available to higher brain regions.Since ortho and retro stimuli are naturally dictated by inhalation and exhalation, the motor signals that control breathing could provide efferent copies of this information to various brain regions.Moreover, even without such motor signals, the mechanical stimulation of neurons in the nasal cavity by breath-driven airflow generates strong breath related signals throughout the olfactory system (Iwata et al., 2017;Shusterman et al., 2011).However, previous experiments have delivered stimuli independently of the breath cycle, and still found significant differences in or-tho versus retro response in cortex (Small et al., 2005;Heilmann and Hummel, 2004) and ORN input to OB (Gautam and Verhagen, 2012b).Thus, there are important breath-independent aspects of ortho/retro coding.Although these different mechanisms for odor modality encoding are all plausible, here we designed our experiments to study the 'bottom up', breath-independent mechanisms.We specifically performed a double tracheotomy which redirected the rats' breath to bypass the nasal cavity entirely and allowed us to deliver ortho and retro stimuli at times that were independent of inhale and exhale times.So the mechanisms of ortho and retro coding that we study here are independent of any motor or mechanosensory signals caused by breathing.
Before odor information reaches the brain, it is processed in the olfactory bulb (OB) and relayed to higher brain regions via excitatory mitral cells (MCs) (and tufted cells).Thus, the OB is critical for determining whether ortho versus retronasal odors are encoded before being conveyed to the brain for processing and perception.More broadly, OB activity is tied to odor perception (Mandairon and Linster, 2009).While the role of inhibitory cells in influencing OB processing of ortho versus retro odors is unknown, inhibition is known to play a key role in many other aspects of OB processing.OB inhibition is known to alter activity patterns that represent odors (Wilson and Laurent, 2005), granule cells that provide inhibition reflect changes in odor concentration (Tan et al., 2010), and OB inhibition levels alter odor discrimination dynamics (Abraham et al., 2010).These facts motivated us to both decrease and increase the levels of GABA A synaptic inhibition via drug manipulations in our experiments, which consist of in vivo recordings of MCs in the OB on double tracheostomized rats.
We show that the mode of olfaction (ortho versus retro) is indeed transmitted to cortical brain regions with a common food odor, ethyl butyrate (EB).We find that encoding is generally good and well above chance level for the intact circuit and even with altered levels of inhibition.Thus downstream brain regions may readily have access to ortho versus retronasal odor modality encoded via firing rate from the OB.Importantly, the altered circuits with both increases (via Muscimol application, a GABA A agonist) and decreases in inhibition (via Bicuculline, a GABA A antagonist) result in worse individual encoding on average.This suggests that the intact OB circuit may be advantageous for encoding ortho versus retro modality.
Our data enables investigation of the OB circuit mechanisms that promote this encoding with a computational model.We show how a simple excitatory-inhibitory (E-I) reciprocally coupled model captures the modality dependent differences in the population firing rate.Our data shows that the GABA A (ant-)agonists do not significantly alter the odor-evoked population firing rate with slower longer-lasting retronasal stimulation, in contrast to ortho stimulation.Our model captures these aspects of our data, exploiting the temporal differences, i.e., the slow retro stimulation leads to short-term synaptic depression of the excitatory synapse to diminish the effects of the inhibition.Furthermore, we model trial-to-trial spike rate variability to be largely consistent with our data, so that the model can encode orth-/retro-like stimuli.We find that average encoding is degraded when inhibition is either increased or decreased, as observed in our data.Given the commonality of reciprocally coupled excitatory-inhibitory circuits, these results may apply outside of the OB, specifically to where stimuli have temporal differences.

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.

Electrophysiological recordings
Data were collected from 8 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 (Gautam and Verhagen, 2012b).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 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 to deliver orthonasal stimuli, 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 second in duration at 1 minute intervals, with a flow rate of 250 ml/min and 1% of saturated vapor.Two odors were used: ethyl butyrate (EB, a food odor) and 1-Hexanol (Hex, a non-food odor).Here we limited our analysis to EB trials because food odors are perceived ortho-and retro-nasally (Small et al., 2005); non-food odors do not naturally occur retronasally.
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).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 band-pass filtered between 300 and 3000 Hz and semiautomatic spike sorting was performed using Klustakwik software, which is well suited to the type of electrode arrays used here (Rossant et al., 2016).
Pharmacology.We topically applied GABA antagonists (bicuculline 20 µM) or agonists (muscimol 20 µM) dissolved in artificial cerebrospinal fluid (ACSF).Gel foam pieces soaked in the ACSF+drug solutions were placed on the surface of OB surrounding the recording electrodes.
Note that the concentration of (ant-)agonists can vary over orders of magnitude; some labs use 10 µM (Li et al., 2019;Sharp and Finger, 2002) and others use up to 10 to 10 4 higher concentration levels (Wachowiak and Cohen, 1999;Kollo et al., 2014).The gel foams were placed 10 minutes before beginning a recording.

Data analysis
Data was collected in vivo from the mitral cell (MC) layer in the olfactory bulb (OB) of multiple anesthetised rats using a multi-electrode array recording.The data consisted of spike recordings of multiple MC spike responses to EB and Hex presented by the two routes of stimulation, orthonasally and retronasally, for a total of 20 trials, 10 for each ortho and retro, with the order alternating (10 ortho before 10 retro, then 10 retro before 10 ortho trials).The spike counts were calculated using 100 ms overlapping time windows.Three separate drug preparations were used in order to analyze inhibitory effects on MC spiking responses: no drug (control), Bicuculline (GABA A antagonist, i.e., decreasing inhibition), and Muscimol (GABA A agonist, i.e., increasing inhibition).Data was obtained from 8 viable rats in total.After spike sorting, we identified and eliminated units (cells) that had firing rates below 0.008 Hz or more than 49 Hz, calculated over the entire recording duration.We took this conservative approach to eliminate potential artifacts yielding units that had unrealistically high or low firing; see Table 2 for the reported number of total rats and individual cells for each drug preparation.Note that some cells around the odor onset (t = 0) will have much larger (> 49 Hz) or smaller firing rates than our criteria to eliminate units, but the entire recording duration is much longer with approximately a minute between trials.

Defining single cell decoding
We used a form of linear discriminant analysis (LDA) for decoding accuracy.Decoding accuracy is defined as the fraction of trials correctly classified by a threshold that maximizes decoding accuracy.Letting s k (t) be the spike train of the k th trial (k = 1, 2, ..., N) that at each time t is either = 0 (no spike) or = δ(t − t j,k ), where t j,k denotes the j th spike on the k th trial .The entries of the spiking rate x are: where t = 0 is the time the odor is presented.There are N = 20 trials total (10 ortho, 10 retro) for each individual cell.We consider N different threshold values θ k : and use both classification schemes (decode stimulus as ortho if x k < θ and retro if x k > θ, or as ortho if x k > θ).From this comprehensive list, we select an optimal threshold θ * and classification scheme that results in the best decoding accuracy.Here ǫ ≈ 10 −16 is machine precision and is there to insure that each trial is categorized as ortho or retro.Since decoding accuracy is the percent of correctly classified trials (ortho or retro) and N = 20, it varies between 0.5 (chance) and up to 1 in increments of 0.05 for each cell.
To systematically compare differences in decoding accuracy, we select a fixed time window for all drug preparations and modalities based on two statistical tests (two sample t−test assuming unequal variances and Wilcoxon rank sum test) and their p−values.We evaluated decoding accuracy from T ev =100 ms to T ev =1 s in 100 ms increments in the evoked state for each drug preparation.
We analyzed differences in mean between no drug preparation and Bicuculline, as well as no drug

Statistical significance
We use three tests to assess whether a given effect is statistical significant: 1) two-sample t−test with unequal variance, 2) Wilcoxon rank sum test, 3) one-way ANOVA.Each of these tests are used to rule out the null hypothesis that the population means are the same in two different categories (Cohen, 2013).Each test is different and has various underlying assumptions: 1) normal samples, 2) non-parametric test assuming independent groups and equal variance, 3) including sample variance to assess differences in means assuming equal variance and normally distributed residuals.We consider three tests to provide a more complete picture of the results and to demonstrate robustness, or lack thereof, for a given effect.
In addition to the p−values that indicate the probability of the null hypothesis holding (no effect, same population means), we also report the effect sizes.We use the conventions outlined in Cohen (2013); Tomczak and Tomczak (2014) that provide formulas for effect sizes, and a qualitative 'rule of thumb' for when effect sizes are: small, medium or large.For 2 groups of length n 1 and n 2 respectively, the effect sizes we use: For t − test, use: For Wilcoxon rank sum test, use: , where ζ is the z-score. (4) For one-way ANOVA, use: SS group SS total , ratio of particular sum of squares. (5) The variables µ j and S 2 j denote the sample mean and variances of the j th group.Effect sizes for t−test and Wilcoxon rank sum test with values: (0, 0.2] are considered small, (0.2, 0.5] are medium, (0.5, 0.8] and above are large (Cohen, 2013;Tomczak and Tomczak, 2014).Effect sizes for one-way ANOVA with values (0, 0.01] are considered small, (0.01, 0.06] are medium, (0.06, 0.14] and above are large (Cohen, 2013).

Population decoding
The two methods are applied to each population of simultaneously recorded MCs.
PCA+LDA: For each recording, Principle Components Analysis (PCA) was applied to the concatenated matrix X of size 20 × M, where M is the number of MCs (Fig 6A ), rows 1-10 correspond to the spiking rate with ortho and rows 11-20 correspond to retro.We only use the first two principal components and apply Linear Discriminate Analysis (LDA) to find a line in the 2-dimensional plane that maximizes classification accuracy (Cunningham and Yu, 2014).We used built-in MATLAB routines pca, fitcdiscr, predict (see GitHub for code).
SVM: We use support vector machine with supervised 'learning' to classify data into ortho or retro.We use a nonlinear classification method with Gaussian kernels, and Bayesian optimization to find the best kernel scale and box constraints for the SVM (we insured that for all recordings and drug preparations, the objective was minimized within the allowable function evaluations).After the SVM is fit to the data, the decoding accuracy is the correct classification rate obtained from the MATLAB function kfoldLoss using 10-fold cross-validation.We used built-in MATLAB routines (see GitHub link for code).

Firing rate model
We use a space-clamped Wilson-Cowan rate model of coupled E-I cells: E (MC, excitatory) and I (inhibitory from PGC, GC, etc).The models of the respective cell firing rate, A j (t) j ∈ (M, P ), are represented by the following ordinary differential equations: where the synaptic term, S j j ∈ {E, I}, is defined by synapses with rise and decay time scales (same τ rise , τ decay for E/I) as follows: In Eqs 6 and 7, I j (t) is the sum of external currents (e.g., stimulus input and other currents) that varies over time to account for stimulus input, w orn,j are the coupling strengths of olfactory receptor neuron (ORN) input to j, w jk are the coupling strengths from cell k to j, and the transfer function F is a threshold linear function: with initial value w IE (0) set to its steady-state value at t = −0.5 s in Eqs( 6)-( 12).
The trial-to-trial variability of spiking rate is modeled as a negative binomial random variable: Given output from the Wilson-Cowan model A E (t), we set: by manually determining the ρ parameter in the negative binomial distribution (see Table 1) for each drug preparation, and setting the parameter r = µ N B ρ/(1 − ρ).

Results
We collected spike data with multi-electrode array recordings of urethane anesthetized rats from the OB mitral cell layer, where each cell (MC) was subject to the same odor delivered with ortho and retro stimulation (Fig 1Ai), repeated for 10 trials for each modality.In addition, two drugs were applied to alter the circuit: a GABA A antagonist (Bicuculline) and a GABA A agonist (Muscimol), see Table 2.For reference, Fig 1Aii shows the population-averaged firing rate (trial-average and population-averaged firing rate) by modality for a given drug preparation, with t = 0 denoting time of odor presentation, held for 1 s.This is commonly referred to as the peri-stimulus time histogram (PSTH).

Individual MCs encode ortho versus retro
To remain agnostic to how higher brain regions decode population activity, we predominately consider individual MCs odor modality encoding, a logical first step (Rolls and Treves, 2011).We use odor-evoked firing rate in a given trial as the source of encoding.Specifically, the firing rate in the k th trial, x k , is the sum of the evoked spike counts normalized by time window T ev = 0.9 s: Ai) Experiment setup to test whether OB MCs encode ortho versus retro stimulation, with example raster plot for 1 MC.Aii) Population-averaged firing rate (trialand population-averaged), i.e., PSTH for each modality and drug preparation; t = 0 denotes when ethyl butyrate (food odor) is presented.Aiii) The proportion of population that fire more for ortho (or retro), i.e., ortho selective ⇔ higher trial averaged firing rate for ortho than retro.B) Schematic of how decoding accuracies are calculated for each MC based on trial firing rates: µ R/O are the trial-averaged firing rates.Shaded regions correspond to incorrect decoding, unshaded regions to correct decoding.C) Distribution of decoding accuracies of MC with 3 different preparations: intact no drug (gray), less inhibition via Bic (green) and more inhibition via Mus (purple).Mean decoding accuracy for no drug is 0.724; with Bic and Mus the means are: 0.680, 0.690 (resp.).Differences are all statistically significant (α = 0.01 with twosample t−test assuming unequal var, Wilcoxon rank sum test, and one-way ANOVA) with mostly 'medium' effect size, see Table 3 A1 for analysis of Off-Response cells).
We define decoding accuracy as the fraction of trials correctly classified by a threshold that maximizes decoding accuracy; this is a measure of information accessible by an ideal observer from each MC individually (Fig 1B).We only consider EB, a common food odor, because typically only food odors are perceived retronasally (Small et al., 2005).It has been shown that nonfood odors artificially delivered retronasally can be perceived as well as food odors (Small et al., 2005).However, numerous studies have shown that human perception of non-food odors delivered retronasally is degraded compared to food odors (Hannum et al., 2021;Bender et al., 2009;Frasnelli et al., 2008).
The OB circuit encodes ortho versus retro very well.Figure 1C shows the spiking activity of individual MCs encode (measured by decoding accuracy) ortho versus retro for the vast majority of cells.MCs have a wide range of decoding accuracies varying between 0.5 (chance) and 1 (perfect decoding).With no drugs, the mean decoding accuracy (0.724) is better than both Bic and Mus, which have lower average decoding accuracies of 0.680, 0.690, respectively.These averages are numerically small in difference, but note that the range of decoding accuracies is only 0.5.Furthermore, the differences in the mean decoding accuracy are all statistically significant (α = 0.01) using three tests: two-sample t−test, Wilcoxon rank sum test, and one-way ANOVA.
Table 3 demonstrates that the p−values are all small, and that the effect size is generally 'medium', with some 'small' (see Statistical significance section in Materials and methods).
Figure 2 shows a more complete picture of the entire (pooled) population's spike rate statistics across drug preparations and modality, beyond the PSTH in Figure 1Aii.There are very strong positive correlations between the mean firing rates and standard deviations (across trials) in all drug preparations and both modalities (see Figure A3, and Table A2).For further discussion on the modulation of evoked spiking variability in the OB, see Ly et al. (2021).

Drug effects on average population spiking
To better understand the dynamics of MC decoding accuracy, we consider how GABA A (ant- )agonists effect population activity.Although population activity is a coarse measure that neglects individual cell heterogeneity, it can still be insightful for determining average trends.A characteristic of an individual MC that might be indicative of encoding fidelity is the distance between the scaled trial-averaged firing rate of ortho and retro stimulation . We expect larger distances to coincide with better decoding accuracy than smaller distances, if all other factors are equal.Firing Rate (Hz) No Drugs

Time (s)
Firing Rate (Hz) No Drugs

Drug effects on trial-to-trial variability
The trial-to-trial variability plays a key role in decoding accuracy differences.We analyze trial variance for each MC to test whether the population average is different across drug preparations.
We find with Mus application, the average trial variability is smallest, followed by Bic, then no drug (i.e., Mus < Bic < ND), for both ortho and retro.The significance measured by p−values and effect size are shown in Tables 4 and 5 for ortho and retro trials, respectively, using the same three statistical tests as before.For ortho stimulation, Mus < ND is significant for two of the three tests (α = 0.01), while Bic < ND and Mus < Bic are only significant with Wilcoxon rank sum test.
For retro stimulation, again only the Wilcoxon rank sum test shows this relationship is significant (all with α = 0.01 except for Bic < ND for retro where p = 0.038).Overall, the Bic < ND relationship is not as strong as the others, and the trend Mus < Bic < ND is stronger with ortho stimulation.Note that the effect size is considered 'small' in all cases, see Tables 4 and 5.

Model connects network dynamics and decoding accuracy
To reveal the neural network dynamics that explain our experimental data results, we use a simple firing rate model.Our rich data set with 3 total drug preparations enables a model framework that is highly constrained by data (Xiao et al., 2021;Ly et al., 2021;Barreiro et al., 2017).In addition, we incorporate known differences in temporal dynamics of ORN inputs to MCs with ortho versus retro stimulation (Sanganahalli et al., 2020;Gautam and Verhagen, 2012b;Furudono et al., 2013;Scott et al., 2007;Carey et al., 2009).Ortho stimulation results in fast increase and fast decrease of ORN inputs to MCs while retro input results in slower increase and slower decrease ORN inputs than ortho (Fig 4B  2021)) that is also computationally expensive.
We first set out to implement an OB model that accounts for the large changes in population PSTH with GABA A (ant-)agonists with ortho stimulation only, and small (or no significant) changes in population PSTH with retro (Fig 3B).We use a reciprocally coupled 2 cell E-I network (Fig 4A ), consistent with OB networks that have many such E-I pairs (granule and periglomerular cells are reciprocally connected to MCs with fast dendrodendritic synapses (Rall et al., 1966)).We model the drug effects Bic (Mus) by decreasing (increasing) the coupling strength from I to E cell by 25% of the baseline value (w EI = 0.2 with no drugs, arbitrary units, see Eqs ( 6)-( 7) and Table 1).
The data shows that affects of drugs altering inhibition are largely absent with slower longerlasting retro stimulation (Fig 3B ), thus finding a neural network principle that effectively removes the inhibitory cell input is imperative, and might be indicative of the core mechanisms.We hypothesize that short-term plasticity on the excitatory synapse (E → I) could be a factor in explaining the differences of population firing between ortho versus retro with GABA A (ant-)agonists because of known temporal differences in ORN inputs.There is indeed evidence of plasticity in the OB circuit (Cang and Isaacson, 2003), specifically for short-term synaptic depression in the main OB (Wang et al., 2012) and from MC to granule cells (Dietz and Murthy, 2005).Since the relative differences in firing rate with ortho are already large without plasticity (1st column in Table 6), honing in on decreasing the firing rate differences in retro (3rd column in Table 6) would make the model more consistent with the data (Fig 3B).We use short-term synaptic depression as a means to nonlinearly change the retro firing rates to diminish the differences because retro has longer lasting stimuli.Although Dietz and Murthy (2005) showed that E → I synapses in OB have both facilitation and depression with facilitation appearing more prominent, we found that short-term facilitation alone made our model much less consistent with our data.(Sanganahalli et al., 2020;Gautam and Verhagen, 2012b;Furudono et al., 2013;Scott et al., 2007;Carey et al., 2009) (see main text).C) Wilson-Cowan firing rate model (Eqs ( 6)-( 12)) output for ortho (top row) and retro (bottom row) stimulation, compared when synaptic depression is included (solid lines) or not (dashed lines).Left column: E cell (MC) output; center column: I cell output; right column: synaptic strengths (W IE ).Data from experiments also shown for comparison (dot-lines) in left column.
Figure 4C shows comparisons of various model outputs (Eqs ( 6)-( 12)).The model indeed follows the firing rate modulation with GABA A (ant-)agonists with both odor modalities, whether or  4).The effects of synaptic depression significantly narrows the differences in retro population PSTH among all 3 drug applications, while with ortho the differences across drug preps remain large (see Table 6 for relative differences).Taken together, synaptic depression is crucial for capturing the data (Fig 3B).With retro stimuli and all 3 drug preparations, the E → I synaptic connections are depressed for longer so that the I-cell fires less, giving a boost to E-cell firing.In both modalities, there is overall higher E (MC) firing with short-term synaptic depression because of less inhibition (for all 3 drug preparations).Synaptic depression effectively diminishes the reciprocal loop so that the inhibitory drug effects are weakened.In principle Mus should have the opposite effect as Bic, but: the coupling strength w EI has moderate variation: 0.15 (Bic), 0. 2 (ND), 0.25 (Mus).Also, we modeled short-term depression so that inhibition is almost completely turned off with Retro, resulting in MC being effectively decoupled from the I cell.Our model does not have realistic levels of synaptic depression, Cang and Isaacson (2003) report moderate levels of synaptic depression.If synaptic depression was limited in our model so that the minimum w IE was 0.5, the ortho firing would be similar but the retro firing would strongly modulate with drug effects, which is inconsistent with our data.And anesthetized rats have less active granule cells (Cang and Isaacson, 2003) than awake, so perhaps other interneurons would have to be active for our model to be more consistent.The exact type of inhibitory neurons in the OB that would have to be involved and their corresponding dynamics with excitation are unknown to us.
With the neural dynamics captured in our model, we turn to connecting the model to decoding accuracy.The firing rate model does not explicitly have trial-to-trial variability, an important piece for decoding accuracy.A natural choice would be to assume Poisson statistics for the spike counts, but this single parameter distribution results in trial variance ordered by the firing rate: Mus<ND<Bic, which is inconsistent with our data and has perfect (100%) decoding accuracy with these firing rate models.Capturing the trial variability trends in our data seems to require tunable parameters that are not directly tied to the (mean) firing rate.We use a simple two parameter model that makes simulation easy and transparent: the negative binomial distribution.Note that there may be other models of trial variability that are more realistic relation to MCs, but that is unlikely to be as pragmatic for capturing firing rate and decoding accuracy trends.Specifically, for each drug (ND/Bic/Mus) and each input (ortho/retro), we fit a negative binomial distribution to the average decoding accuracy.The two distribution parameters were chosen so that: 1) the mean of distribution coincided with the model spiking rate column), and 2) the decoding accuracy matches the averages in the experimental data.
The negative binomial distribution has distribution: x with parameters: ρ ∈ (0, 1) and r > 0 that are not directly related to the physiology.We manually (by trial and error) determine the parameter ρ for each modality and drug preparation, which are all the same ρ except for retro with Bic and Mus (see Table 1), after which the r parameter is determined by: since the mean of the negative binomial is r(1 − ρ) ρ .We chose ρ for each regime so that the resulting decoding accuracy is close to the averages from the experiments: 0.724 (no drug), 0.680 (Bic), and 0.690 (Mus).
The results are summarized in Fig 5A .We simulated 50,000 trials and found the optimal threshold to get decoding accuracies that match the data.We check whether our model results are consistent with the trial-to-trial variability of our data by calculating the variances of firing rate over trials, which is simply r(1 − ρ) ρ 2 for the negative binomial distribution.The values listed in Table 7 show that trial variance is consistent with the data: Mus < ND < Bic, except for retro Mus firing that has the largest variance.Recall in our data that the trends for variance of firing rate was weakest with retro stimulation with larger p−values (Table 5).
equally spaced points and M R ∈ {−6, . . ., 9} Hz on 31 equally spaced points.Once a mean is specified, we use the same ρ (Table 1) and set r * to correspond to the mean µ * : we simulate 50,000 trials to get a decoding accuracy value for each set of parameters.We get

Population decoding results
Thus far we have focused on decoding odor modality using the firing rate of individual MCs.
While this approach avoids making assumptions about how higher level cortical regions might use the population as a whole (Rolls and Treves, 2011), population coding is known to be important for complicated tasks (Saxena and Cunningham, 2019).We apply standard population coding metrics to assess how well odor modality can be decoded from populations of MCs.We show that modality is well-encoded, but that the differences with GABA A maniuplations are minor.We applied two approaches to population decoding.First, we used principal component analysis (PCA) to project the population response onto a two-dimensional subspace that captures the most variance (first 2 PCs), and then applied linear discriminant analysis (LDA) to find an optimal decoder in that subspace (PCA+LDA).Second, we used a support vector machine (SVM), a common supervised learning algorithm to (nonlinearly) classify data.We use the first 2 PCs in PCA+LDA because it is a natural extension of 1D decoding with LDA that was the primary focus of prior sections.The SVM is at another extreme using all dimensions possible; together, these two provide a glimpse into the various accuracies of population decoding.

Decoding Accuracy Fraction of Recordings
With both methods, we find odor modality is still encoded well above chance even when inhibition is altered (Fig 6B ,C).Not only is the average decoding accuracy (across separate recordings) higher than before (Fig 1C) but there are no longer decoding accuracies near chance level (0.5).
This is perhaps expected because these metrics (PCA+LDA, SVM) exploit the higher dimensional population activity to perform optimally in a binary classification problem when in reality the MCs encode other aspects of stimuli.Note that nearly 35% of recordings with no drug and Mus have perfect decoding accuracy (Fig 6B ,C).
In contrast to the case of a single MC's encoding (Fig 1C), we do not find any statistically significant differences between the mean decoding accuracy across drug preparations.In fact, the minimum p−value for all relationships and all 3 statistical tests is p ≥ 0.32 for PCA+LDA, p ≥ 0.13 for SVM.We expect population coding to perform better than individual cell coding, especially given the large dimensionality of population activity and that these metrics (PCA+LDA, SVM) were optimized for a binary classification task.

Discussion
We have shown that MCs in the OB encode odor modality with our experimental data, and use an OB model inspired by many prior studies and constrained by our rich data to analyze plausible network components that support these dynamics.Our study uncovers further details for how ortho versus retro information is encoded to early cortical regions via the OB, which is consistent with reports that humans are able to discriminate whether food odors are delivered ortho or retronasally without being told which modality (Frasnelli et al., 2008).Imaging studies of the human mouth have shown that with food odors, ortho stimuli delivery versus retro are two distinct modalities, rather than 2 different stimulus routes to the same modality, independent of odor intensity (Bender et al., 2009).Based on in vivo rat data, we found significant differences in decoding accuracies (our proxy for encoding assuming an ideal observer) for classifying ortho and retronasal odors across different drug preparations.The intact (no drug) circuit had the best (average) decoding accuracy, suggesting that the OB might be structured to encode this important aspect of odor input.
Our model framework used prior results, data constraints (Xiao et al., 2021;Ly et al., 2021;Barreiro et al., 2017), and simplicity to account for the results in our data.Many prior studies (Sanganahalli et al., 2020;Gautam and Verhagen, 2012b;Furudono et al., 2013;Scott et al., 2007;Carey et al., 2009) showed that there are differences in the temporal dynamics of ORN response to ortho versus retro stimulation.We used this to investigate the OB circuit components that promote efficient coding of individual mitral cells with drug manipulations of inhibitory synapses.Further data analysis revealed that inhibitory drugs have a stronger effect on the population firing with ortho stimulation than retro.Using these insights, we constructed a standard firing rate model that captures the various drug effects for ortho versus retro stimulus, using only a pair of E − I cells.Including synaptic depression was crucial for capturing the different drug effects of firing rate.We did not explore varying the many parameters of the model to find a 'best fit' to data, There have been prior studies that implemented an OB circuit model with varying levels of realism, but none that we are aware of that uses an OB model to account for coding of odor modality and firing rate dynamics.We previously developed a realistic biophysical OB network model (Craft et al., 2021) and used it to account for different neural network dynamics (trial-averaged) with ortho/retro inputs, but did not focus on decoding accuracy for a given trial (which is more computationally taxing).In our prior work (Craft et al., 2021), the PSTH for ortho versus retro are different than here (Figs 1Aii, 3B) because we previously used only 1 rat versus 8 rats here.There are biophysical models of particular OB cells (Viertel and Borisyuk, 2019;Li and Cleland, 2013) and of the OB network (Bathellier et al., 2006;Li and Cleland, 2017), and other have used models (firing rate or Bayesian) to investigate coding of mixed odors (Grabska-Barwińska et al., 2017) and learning new odors (and their concentration levels) (Hiratani and Latham, 2020), but none that we know of that focus on coding of odor modality.Note that muscimol increases tonic E cell inhibition independent of I cell firing; our simple model does not disentangle these two components.
A limitation of this study is the use of one odor (EB, a food odor) with one concentration.We collected data from a non-food odor (1-Hexanol), but since typically only food odors experienced retronasally (Small et al., 2005), we did not include 1-Hexanol in our study.It is well-known that non-food odors delivered retronasally are not perceived as well as retro food odors by humans (Hannum et al., 2021;Bender et al., 2009;Frasnelli et al., 2008).Other food odors have different physicochemical properties that impact differences in ortho versus retro responses (Scott et al., 2007), so the generalization of our results merits further investigation.Firing rates vary with odor concentration (Tan et al., 2010) and is thus a confounding factor in distinguishing ortho vs retro from firing rates alone.Whether higher brain regions rely on other signals to account for concentration or different patterns of activation in OB to reliably estimate ortho/retro is important but beyond the scope of this study.Our data was collected from anesthetized rats with forced air to model ortho/retro (Gautam and Verhagen, 2012b,a;Gautam et al., 2014) which has clear advantages: control and a fair decoding 'task' so that ortho and retro are mechanically similar with same stimulus duration.In addition to the intact circuit, our data also had both increased and decreased GABA A synaptic inhibition with direct manipulations in the OB.Via pooling data from many rats, we had a large number MCs.Another limitation is that the rats were not awake and were not eating.We are currently unaware of any data with awake rodents that record MC spiking activity for the purposes of comparing ortho and retro (food) odors during feeding.There are freely available MC spike data with awake rodents (mice in Bolding andFranks (2018, 2021)) but not during feeding and not with direct manipulation of OB inhibitory synapses.
Since the temporal dynamics of odors can be crucial in OB processing, we systematically varied the evoked time window and even considered removing the first 300 ms of the evoked state (see Fig A2).We found that the population-average (and standard deviation) decoding accuracy varied slightly (Fig A2B ,D), and indeed with the first part of the evoked state removed the average decoding accuracy increased.In the context of our study, the temporal patterning differences between ortho and retro have minor effects on decoding.In awake, intact animals, the timing of the breath cycle determines the timing of ortho and retro stimulation (Shusterman et al., 2011).
In principle, brain areas downstream from OB could use a breath-timing-based decoding scheme to entirely distinguish ortho from retro.However, our experiments show that timing is not the only mechanism that impacts ortho versus retro coding.In our experiments, the breath of the anesthetized rats did not pass through the nose at all; it was redirected through a tracheotomy tube.Thus, the timing of ortho and retro signals in our experiments was decoupled from breath timing and we kept the stimulation parameters (concentration and flow rate) the same, changing only the direction of airflow.Yet, we found differences in ortho and retro response in OB.Also, when instructed not to breath and actively sniff, humans can correctly report ortho versus retro odor modality with food odor without being told the modality (Frasnelli et al., 2008).Thus, we conclude that the OB encodes ortho and retro stimuli not just with the timing of the breath, but also with other mechanisms.Nonetheless, it is likely that intact, awake animals would employ both timing (inhale vs exhale) and the mechanisms we identify here.Disentangling the relative importance of these different possible contributions to ortho versus retro coding will require further experiments.
Our results predominately center on population averages of individual MC decoding accuracy because it is a logical first step (Rolls and Treves, 2011).Given the importance of population coding for many complicated tasks (Saxena and Cunningham, 2019), we also considered 2 standard population decoding accuracies measures: PCA+LDA and SVM and found in both cases that the intact circuit was not statistically significant better (on average).However, decoding accuracy was again much higher than chance levels and better compared to average individual MC decoding accuracy.Developing a heterogeneous OB model that accounts for the population coding results and dynamics is a potential next step but beyond the scope of this current study.shows that the observation that ortho stim is more effected with GABA A manipulation than retro has consequences for decoding accuracy modulation, even if the predictions are imperfect.

Orthonasal
preparation and Muscimol.This resulted in 4 different p−values for each time window; we found that T ev = 900 ms had the most significant p−values for both combinations of no drug/Bicuculline and no drug/Muscimol, see Fig A2A,C.Figure A2B,D shows the population decoding accuracies (mean across the population and one standard deviation in shaded region).

Figure 1 :
Figure1: Individual mitral cells encode modality.Ai) Experiment setup to test whether OB MCs encode ortho versus retro stimulation, with example raster plot for 1 MC.Aii) Population-averaged firing rate (trialand population-averaged), i.e., PSTH for each modality and drug preparation; t = 0 denotes when ethyl butyrate (food odor) is presented.Aiii) The proportion of population that fire more for ortho (or retro), i.e., ortho selective ⇔ higher trial averaged firing rate for ortho than retro.B) Schematic of how decoding accuracies are calculated for each MC based on trial firing rates: µ R/O are the trial-averaged firing rates.Shaded regions correspond to incorrect decoding, unshaded regions to correct decoding.C) Distribution of decoding accuracies of MC with 3 different preparations: intact no drug (gray), less inhibition via Bic (green) and more inhibition via Mus (purple).Mean decoding accuracy for no drug is 0.724; with Bic and Mus the means are: 0.680, 0.690 (resp.).Differences are all statistically significant (α = 0.01 with twosample t−test assuming unequal var, Wilcoxon rank sum test, and one-way ANOVA) with mostly 'medium' effect size, see Table3.Windows were determine systematically to maximize p−values in testing average decoding accuracy, see Fig A2A,C and Materials and methods.
Figure1: Individual mitral cells encode modality.Ai) Experiment setup to test whether OB MCs encode ortho versus retro stimulation, with example raster plot for 1 MC.Aii) Population-averaged firing rate (trialand population-averaged), i.e., PSTH for each modality and drug preparation; t = 0 denotes when ethyl butyrate (food odor) is presented.Aiii) The proportion of population that fire more for ortho (or retro), i.e., ortho selective ⇔ higher trial averaged firing rate for ortho than retro.B) Schematic of how decoding accuracies are calculated for each MC based on trial firing rates: µ R/O are the trial-averaged firing rates.Shaded regions correspond to incorrect decoding, unshaded regions to correct decoding.C) Distribution of decoding accuracies of MC with 3 different preparations: intact no drug (gray), less inhibition via Bic (green) and more inhibition via Mus (purple).Mean decoding accuracy for no drug is 0.724; with Bic and Mus the means are: 0.680, 0.690 (resp.).Differences are all statistically significant (α = 0.01 with twosample t−test assuming unequal var, Wilcoxon rank sum test, and one-way ANOVA) with mostly 'medium' effect size, see Table3.Windows were determine systematically to maximize p−values in testing average decoding accuracy, see Fig A2A,C and Materials and methods.

Figure 2 :
Figure 2: Detailed spiking rate statistics of pooled MC populations.A) Colormap showing trial-averaged spiking rate of individual MC (individual PSTH), each of the 3 columns is sorted by time-averaged firing rate (4 s total, 2 s in spontaneous and evoked, but only showing 1 s of spontaneous for brevity) of ortho firing rate, so using 3 ways to sort.B) Colormap showing trial variability of spiking rate (measured by standard deviation across trials) of individual MC, using the same corresponding index in top row (ortho firing rate).To see relative differences across MC, the color scales were effectively thresholded so that the larger values did not obscure the entire population.See Fig A3 for more on firing rate and trial-to-trial variability relationship.

Figure 3 :
Figure 3: Population activity: retronasal response is more robust to changes in OB inhibition.A) Distance between scaled trial-averaged mean | µ R σ R − µ O σ O | is correlated with decoding accuracy.The 'Correlation' is Pearson's correlation coefficient, and 'Rank Correlation' is Spearmen's.B) Drug effects on evoked population PSTH are different with statistically significant effects only with ortho stimulation (Bi) but not retro (Bii).The right-panels with 3 subplots show the time series of p−values from two-sample t−test assuming unequal variances, comparing pairs of drug preparations.

Figure 4 :
Figure4: Firing rate model captures effects of drugs, with better alignment with data when including synaptic depression.A) Schematic of reciprocally coupled E-I pair.B) Fixed ORN input profiles for ortho and retro based on experiments(Sanganahalli et al., 2020;Gautam and Verhagen, 2012b;Furudono et al., 2013;Scott et al., 2007;Carey et al., 2009) (see main text).C) Wilson-Cowan firing rate model (Eqs (6)-(12)) output for ortho (top row) and retro (bottom row) stimulation, compared when synaptic depression is included (solid lines) or not (dashed lines).Left column: E cell (MC) output; center column: I cell output; right column: synaptic strengths (W IE ).Data from experiments also shown for comparison (dot-lines) in left column.

Figure 5 :
Figure 5: Model decoding accuracies with simulated trial variability.A): Histograms of simulated firing rate assuming negative binomial trial-to-trial variability with the mean set to: 1 T ev Tev Figure 5C shows how decoding accuracy varies with scaled distance | µ * R σ R − µ * O σ O |.The decoding accuracies can vary for a given drug preparation, but the ND and Bic have strong positive correlations (like the data, see Fig 3A), but Mus has strong negative correlations that is inconsistent with the data.

Figure 6 :
Figure6: Population decoding of ortho versus retro is reliable and robust to altered inhibition.A) Each drug preparation has a different number of recordings and each recording has a different number of simultaneously recorded MCs.Showing recording number (sorted) and MC population size.B) Decoding accuracies (one for each simultaneous recording) for all drug preparations using a 2D PCA dimension reduction following by LDA.Inset shows variance explained by PCA (sorted) by recording.C) Decoding accuracies using full dimensional support vector machine (SVM) classification.Average decoding accuracies by both methods (B)-C)) are statistically indistinguishable (minimum p−value for all relationships and all 3 statistical tests is: p ≥ 0.32 for PCA+LDA, p ≥ 0.13 for SVM).
rather we wanted a simple yet informative model that demonstrates a plausible mechanism.The model PSTH fits to the data are imperfect and qualitative (Fig 4C): the spontaneous firing of ortho and retro in the model are different because of differences in the ORN input strengths, the ortho model firing decays faster than the data.In fact for retro it looks like the model fits the data better without synaptic depression (although the differences with GABA A manipulations are large and inconsistent with the data).In theory the model could be better fit to the data PSTH, for example if the ORN inputs(Fig 4B)  are held fixed and an optimization routine is performed on a set of parameters, but then the decoding accuracy would have to be accounted for(Fig 5A), perhaps in the objective.We thus opted to manually determine the model parameters.Lastly, we simulated trial-to-trial variability by assuming a negative binomial distribution with mean defined by the firing rate model output.The resulting model captured both decoding accuracy trends (ND is better) and some aspects of firing rate trial variability.

Figure A1 :
Figure A1: Supplemental data analysis: extended time of PSTHs in main text.Extended timeframe for Figure 1Aii (A) and Figure 3B (B).
Without synaptic depression, we fix w IE = 1.With synaptic depression, w IE (t) is governed by:

Table 2 :
Number of rats and respective individual cells for each drug preparation subject to the food odor ethyl butyrate (EB).

Table 3 :
Population Decoding Accuracy: significance measured by p−values of average decoding accuracies over MC, shown in Figure 1C.Using 3 tests: two-sample t−test assuming unequal variance, Wilcoxon rank sum test, and one-way ANOVA.Bottom table shows effect size (see Materials and methods).

Table 4 :
Ortho stimulation only: significance (p−values) of average (over MC) trial variance differences between drugs with various statistical tests.Using same 3 tests as before; ND='no drug'; bottom table shows effect size.

Table 5 :
Retro stimulation only: significance (p−values) of average (over MC) trial variance differences between drugs with various statistical tests.Using same 3 tests as before; ND='no drug'; bottom table shows effect size.

Table 6 :
Relative differences of firing rate model between no drug and inhibitory drugs in both modalities.All percent differences (top table) are relative to no drug:A E (t) dt (Hz) , i.e., the components of the relative differences.
not synaptic depression is included.Specifically, ortho MC firing rates (top left in Fig 4) with different drug preparations are well-separated, while the retro MC firing rates are pretty close together (bottom left in Fig Table 1) in negative binomial distribution generally results in no drug having best decoding accuracy (see main text for details).C) Same model results in B) but showing relationship with distance between (scaled) means; red circles correspond to regimes in A).

Table 7 :
Simulated trial variance of the model: E (t) dt (solid curves in Fig 4C, left column) and perturb these values to get different means: µ