Skip to main content


Front. Comput. Neurosci., 13 July 2020
Volume 14 - 2020 |

Top-Down Control of Inhibitory Granule Cells in the Main Olfactory Bulb Reshapes Neural Dynamics Giving Rise to a Diversity of Computations

  • 1Department of Brain and Cognitive Sciences, University of Rochester, Rochester, NY, United States
  • 2Department of Neuroscience, University of Rochester School of Medicine and Dentistry, Rochester, NY, United States

Growing evidence shows that top-down projections from excitatory neurons in piriform cortex selectively synapse onto local inhibitory granule cells in the main olfactory bulb, effectively gating their own inputs by controlling inhibition. An open question in olfaction is the role this feedback plays in shaping the dynamics of local circuits, and the resultant computational benefits it provides. Using rate models of neuronal firing in a network consisting of excitatory mitral and tufted cells, inhibitory granule cells and top-down piriform cortical neurons, we found that changes in the weight of feedback to inhibitory neurons generated diverse network dynamics and complex transitions between these dynamics. Changes in the weight of top-down feedback supported a number of computations, including both pattern separation and oscillatory synchrony. Additionally, the network could generate gamma oscillations though a mechanism we termed Top-down control of Inhibitory Neuron Gamma (TING). Collectively, these functions arose from a codimension-2 bifurcation in the dynamical system. Our results highlight a key role for this top-down feedback, gating inhibition to facilitate often diametrically different computations.


Growing evidence suggests that top-down centrifugal feedback from higher cortical areas specifically target inhibitory interneurons in primary sensory regions. In the olfactory system, axons from excitatory neurons in the piriform cortex (PCx) synapse onto the inhibitory granule cells in olfactory bulb (OB), whereby they can modulate the function of the mitral/tufted cells (M/T), the principal relays of olfactory information from the bulb to the brain (Shipley and Adamek, 1984; Boyd et al., 2012; Markopoulos et al., 2012; Oswald and Urban, 2012b; Padmanabhan et al., 2016, 2019). This circuit motif results in piriform cortical neurons receiving input from only excitatory M/T cells but exerting influence on the local circuit dynamics in the OB via inhibitory populations. Consequently, the information relayed to piriform cortical neurons comes from M/T cells, but feedback intervenes in network dynamics through the local inhibitory granule cells. Although this network motif constitutes major feature of the olfactory system, the computational role of this top-down control of inhibition remains largely unknown.

A number of studies have previously explored the dynamics of M/T cells and granule cells in OB as a two-population network of excitatory (E) and inhibitory (I) neurons (Cleland and Linster, 2005; Brea et al., 2009; Kay et al., 2009; Li and Cleland, 2017) in parallel with a broader literature on excitatory-inhibitory (E-I) networks (Wilson and Cowan, 1972; Ermentrout and Kopell, 1990; Tsodyks et al., 1997; Tiesinga and Sejnowski, 2009; Ledoux and Brunel, 2011; Franci et al., 2018). The studies in olfaction have revealed not only the mechanisms by which these dynamics emerge, but also how changes in the oscillatory power (Nusser et al., 2001) result in alterations in behavior, including in odor discrimination tasks (Abraham et al., 2010). The recent evidence that cortical feedback directly synapses onto inhibitory interneurons (Boyd et al., 2012; Markopoulos et al., 2012) suggests that the local dynamics of excitatory and inhibitory neurons in OB can be gated by centrifugal projections from olfactory cortex. How the local E-I network's activity in the bulb is changed by centrifugal input, what these changes mean more broadly for neural dynamics in the early olfactory system, and the role of these dynamics play in neural computation remains an open question.

To address this, we built a three-node network model consisting of an excitatory population of mitral/tufted cells (M), an inhibitory population of granule cells (G) and a top-down population of pyramidal/semilunar cells (P) in PCx, and studied how firing rate dynamics were influenced by top-down weights onto inhibition. Changing the weight of the top-down connections to local inhibitory neurons reshaped the dynamics of the local E-I circuit in a way that enhanced sensory discrimination as well as generated oscillatory synchrony including entraining gamma oscillations in the local circuit [Top-down control of Inhibitory Neuron Gamma, (TING)]. Finally, the mechanism underlying the dynamics, as well as the functional roles played by top-down control of inhibition occurred via a codimension-2 bifurcation in the dynamical system. By gating the weight of connections from piriform cortex to the inhibitory neurons in the bulb, a number of seemingly disparate computations could be supported by a single circuit, providing an additional framework for the diversity of inhibitory interneuron function in the olfactory bulb.

Materials and Methods

Network Model

The network model was composed of three nodes, the local excitatory population, corresponding to mitral and tufted cells (M) and inhibitory population corresponding to granule cells (G) which were reciprocally coupled, and a top-down population corresponding to the principal neurons in piriform cortex (P) that received input from the local M population and projected back to the inhibitory G population (Figure 1A). In the model, ri (t), i = 1, 2, 3 represented the firing rates of the three neuron populations, respectively, whose dynamics were determined by Wilson-Cowan equations (Wilson and Cowan, 1972) as follows:

{τ1r1.=r1+S(w11r1+w12r2+μ)τ2r2.=r2+S(w21r1+w22r2+w23r3)τ3r3.=r3+tanh(w31r1+w33r3)    (1)

where S is the sigmoid function:

S(x)=11+e-x    (2)

which described the non-linear relationship between the mean synaptic input and average firing rate (normalized to a range between 0 and 1). The parameter τi, i = 1, 2, 3 was the time constant for each population, characterizing how quickly the dynamics of each population evolved. The mitral/tufted cell population (M) received an external stimulus μ, that represented the only external input to the system. The connection weight from population j to population i was denoted by wij, i, j = 1, 2, 3, among which w11, w21, w31, w23, w33 > 0 and w12, w22 < 0. The connection weight wij, i, j = 1, 2, 3 represented the average synaptic input received by the neuron population i from the population j. Throughout this paper, we set the parameters as follows: w11 = 8.7, w12 = −10, w21 = 7.0, w22 = −13, w31 = 1.5, w33 = 0.5. The parameters in the model were chosen based on previous studies of Wilson-Cowan rate model (Ermentrout and Kopell, 1990; Ledoux and Brunel, 2011; Veltz and Sejnowski, 2015), and their relative values were adjusted according to experimentally recorded excitatory and inhibitory postsynaptic inputs of M/T and granule cells in OB (Urban and Sakmann, 2002; Egger et al., 2005; Kapoor and Urban, 2006). For instance, the value of recurrent excitation (w11) was determined by mapping recorded excitatory post-synaptic potentials (EPSP) in M/T cells and the weight of inhibition from granule cells (w12) was determined using similar mappings of whole cell recordings of inhibitory post-synaptic potentials (IPSP) in M/T cells (Urban and Sakmann, 2002). The weights between populations were determined by integrating synaptic potentials with known connectivity densities using trans-synaptic viral tracers that allowed for estimates of the number of pre-synaptic cells for each population (Willhite et al., 2006; Miyamichi et al., 2011; Padmanabhan et al., 2016) to yield the relative weights for instance, |w12| > |w21|. Estimates of synaptic weights for the projections to piriform cortex and the feedback connections were estimated from Franks and Isaacson (2006), Suzuki and Bekkers (2011), and Boyd et al. (2012). Time constants, for example the time constant τ2 for granule cells, were derived from data using both calcium imaging and whole cell recordings across the types of neurons that constituted the populations in our model (Franks and Isaacson, 2006; Kapoor and Urban, 2006; Suzuki and Bekkers, 2011). As there was no inhibitory synaptic input into the piriform pyramidal/semilunar cell population (P), the combination w31r1 + w33r3 was non-negative due to w31 > 0 and w33 > 0. Sigmoid function (2) for the third equation of system (1) would mean that the lowest r3 value that an equilibrium could reach would be 0.5. In this framework, even without an input, at least half of the pyramidal/semilunar cells would keep firing. Thus, we used the hyperbolic tangent function tanh (x) in order ensure that piriform cortical cell firing rates were low in the absence of odors (Stettler and Axel, 2009; Davison and Ehlers, 2011), and to exploit the entire range [0, 1] of r3 (Figure 1C). However, it should be noted that such choice of non-linearity did not affect our findings since all bifurcations supporting the computations of the network we identified were found in the system using the sigmoid function (2).


Figure 1. Reduced network model exhibiting complex dynamics. (A) A schematic diagram illustrating the topology of the reduced network model. Each node denotes a neuron population and the connection weights are defined by wij (excitatory: arrows; inhibitory: circles). (B) Network responses [color-coded firing rates match population nodes in (A)] to different levels of stimulus strength (left: μlow = 0.1, middle: μmed = 1.5, right: μhigh = 3.0) with the other parameters fixed. (C) A distribution of all possible steady states r3 over a wide range of parameter choices: μ ∈ [0, 5], w31 ∈ [0, 5], w23 ∈ [0, 30] shows the diversity of responses the network can generate. (D) Trajectories of the firing rate responses plotted in (B) are visualized in the phase space spanned by (r1, r2, r3). The color bar indicates stimulus strengths for three representative levels as in (B).

Definition of Period of a Limit Cycle

From the perspective of dynamical system, a limit cycle in the phase space corresponds to the oscillation of firing rates in the temporal space. Since the time constants in our model had units of millisecond, the frequency of oscillations was defined as 1, 000/T, with T denoting the period of limit cycles, which was defined as follows: if ri (t), i = 1, 2, 3 denote the firing rate of neural population in the model, then a limit cycle satisfies the periodicity:

ri(t)=ri(t+T), i=1,2,3    (3)

for some T > 0 and all t∈ ℝ. The minimal T for which the above equality holds is the period of the limit cycle.

Metric Definition

As representations of the network (ω-limit sets) could take on different forms, an equilibrium or a limit cycle, we defined two quantitative metrics: dE(Ω1,Ω2) and dS(r3Ω1,r3Ω2) which served to measure the distance between different types of ω-limit sets in responses to any given stimulus pair (μ1, μ2) where μ2 = μ1 + Δ. dE(Ω1,Ω2) denoted the average Euclidean distance between Ω1 and Ω2 in the three-dimensional (3D) phase space of firing rates (r1, r2, r3), and the spectrum distance dS(r3Ω1,r3Ω2) was a sum of the squared differences between both direct components (DC) and alternating components (AC) in the amplitude-frequency domain of the Fourier transforms to the signals r3 (t) associated with the two stimuli.

The Euclidean distance was defined as follows: supposing that Ω1 and Ω2 are two ω-limit sets composed of N1 and N2 discrete points in three-dimensional phase space (r1, r2, r3), respectively, denoted as {α1, α2, …, αN1} and {β1, β2, …, βN2} where αi and βj are three-dimensional vectors, then

dE=d(αi, βj), i=1,, N1, j=1,, N2    (4)

Where 〈·〉 denotes the average and d(x,y)=(x-y)T(x-y) which is the standard distance between any two points in the three-dimensional Euclidean space. Note that in the case of equilibria, the number of discrete points in the ω-limit set was one: N = 1, and in the case of a limit cycle, we set N = T/dt, where T was the period and dt was time bin for numerical integration.

The Euclidean distance dE worked well in measuring the distance between two equilibria in the phase space, but when one ω-limit set was a limit cycle, the averaging operation in the definition made it only a coarse and lagged estimate of the separation for equilibrium vs. limit cycle and limit cycle vs. limit cycle. In particular, the activity of the P population should be decodable with respect to the stimulus information, something that was problematic for when using dE. Therefore, we defined the spectrum distance dS to address the question of distances between representations that were sensitive to those representations being oscillations. To calculate the spectrum distance dS between two ω-limit sets Ω1 and Ω2, only the sequence of their r3 component was decomposed by Fourier transform which converted the firing rate signal in the temporal domain into a representation in the frequency domain.

The single-sided amplitude spectrum for the Fourier transform of the firing rate signal r3 (t), was used to obtain peaks around frequency values. For an equilibrium corresponding to constant firing rate r3 = A, there existed only one peak around zero frequency with its amplitude proportional to A, since the Fourier transform of a constant function is a delta-function. We referred this component as the direct component (DC) of the signal. For the case of a limit cycle, in addition to one peak around zero frequency, there existed another peak around frequency 1, 000/T where T denoted the period of the limit cycle. We referred this additional peak of a limit cycle as the alternating component (AC). Thus, the spectrum distance dS was a sum of the differences between both direct components (DC) of two limit sets and alternating components (AC) of two limit sets, which was formalized as follows: supposing that D1 and D2 were the amplitudes of the peaks at zero frequency for two ω-limit sets Ω1 and Ω2, and ai = (fi, Ai), i = 1, 2 denoted the corresponding alternating components of Ωi, where fi was the non-zero frequency and Ai was the amplitude of the peak around fi, then we have

dS=|D1-D2|+d(a1,a2)    (5)

Note that we set ai = (0, 0) if the ω-limit sets Ωi was an equilibrium. Thus, when the two limit sets were both equilibria, dS only contained the first term measuring the difference between the direct components. In this case, dS was only a linear projection (up to a constant factor) of the Euclidean distance dE.

When both Ω1 and Ω2 were equilibria, the spectrum distance dS was a projection of the Euclidean distance dE onto the r3 axis. The spectrum distance was however more sensitive to bifurcations when one of the two ω-limit sets Ω1 and Ω2 transitioned into a limit cycle as the frequency of the limit cycle started from non-zero values at the onset of a bifurcation (Figure 2C1), causing discontinuous jumps in the spectrum distance. However, for the purpose of pattern separation, the two metrics did not give qualitatively different results when assessing the distances due to changes the feedback weight w23 (Figure 4). Additionally, we found that the non-monotonic dependence of distance in both dE and dS on the feedback weight (w23max) persisted, and that an optimal value for any given pair of stimuli could be found (Figure 3D).


Figure 2. Network dynamics are controlled by top-down modulation. (A) Top: changing top-down input w23 reshapes network activity to the same stimulus (color bar indicates values of w23); bottom: firing rates at three representative values of w23 while the stimulus is held constant (μ = 1.5). (B) Trajectories in the phase space for the same w23 as in (A); left: oscillations occur around w23 = 5.5 for μ = 1.5 (same as A); right: oscillations occur around w23 = 10.5 for a different stimulus μ = 0.6, revealing that the dynamics are diverse across different combinations of stimuli and top-down input. (C) Modulation of oscillation frequency by top-down input (for μ = 1.5). (C1) Dependence of frequency on top-down input w23. Inset: time series of r3 (t) for two example values of w23 (squares). (C2) Frequency modulation by top-down input occurred over a range of feedforward drive w31. (C3) A distribution of oscillation frequencies that can be generated by the network for all possible combinations of w23 ∈ [0, 30] and μ ∈ [0, 5]. (D) Similar to (C) but for oscillation amplitudes.


Figure 3. Pattern separation via top-down control. (A) A schematic diagram illustrating the separation maximization between the response patterns to a pair of stimuli (μ1 and μ2) by changing top-down input w23 (color bar indicates values of w23). (B) Time series of r3 (t) in response to μ1 = 1.0 and μ2 = 1.1 at three representative values of w23 with the spectrum distance dS indicated. The two responses r3 (t) are close at some top-down inputs (left: w23 = 3; right: w23 = 10), but pushed apart at other top-down input (middle: w23 = 3.8). (C) Phase trajectories and network representations of two stimuli from which the Euclidean distance dE is calculated. From left to right the top-down input w23 correspond to those in (B) for the same stimulus pair (μ1 = 1.0 and μ2 = 1.1). (D) Non-monotonic dependence of both dE and dS on top-down input w23, with maximum achieved at w23 = 3.8. The squares are color coded as in (B,C).

Bifurcation Analysis

We denoted the dynamical system of Equation (1) as a parameterized form

r˙=f(r,Θ)    (6)

where r ∈ ℝ3 was the vector of firing rates and Θ∈ ℝp was the vector of parameters. The vector field f = (f1,f2,f3)T was a smooth function on some open set of ℝ3 × ℝp. The dimensionality of Θ could be up to eight dimensions maximally to include all connection weights wij and the external stimulus μ. However, since we were only interested in the top-down control, Θ was restricted to be two dimensions (p = 2) including the top-down weight w23 and the external stimulus μ. All the other parameters were fixed as constants determined based on previous experimental work (Whittington et al., 2000). As the system Equation (6) had an equilibrium at (r, Θ) = (r0, Θ0), i.e.,

f(r0, Θ0)=0    (7)

the stability of this equilibrium could be determined from the linearized vector field of Equation (6) given by

ξ˙=Drf(r0,Θ0)ξ,ξ 3    (8)

where Drffr was the Jacobian matrix of the vector field f.

If none of the eigenvalues of Drf (r0, Θ0) lied on the imaginary axis (i.e., the equilibrium was hyperbolic), the local stability of (r0, Θ0) in the non-linear system (6) could be determined by the linear system (8). The equilibrium was stable if all eigenvalues of Drf (r0, Θ0) had negative real parts. In the case of a hyperbolic equilibrium, varying slightly the parameter Θ would not change the stability as taking Equation (7) and the invertibility of Drf (r0, Θ0), there existed a unique smooth function h : ℝp → ℝ3 such that

h(Θ0)=r0 and f(h(Θ),Θ)=0    (9)

for Θ sufficiently close to Θ0. By continuity of the eigenvalues with respect to parameters, Drf (h (Θ), Θ) had no eigenvalue on the imaginary axis for Θ sufficiently close to Θ0. Therefore, the hyperbolicity of the equilibrium persisted and its stability type remained unchanged for in close vicinity of Θ0. By contrast, when some of the eigenvalues of Drf (r0, Θ0) lied on the imaginary axis, for example, a zero eigenvalue or a pair of purely imaginary eigenvalues, new topologically different dynamical behaviors occurred by a small change in Θ. Equilibria could be created or annihilated, and periodic dynamics could emerge.

The parameterized system (6) thus underwent a bifurcation at (r0, Θ0) if the Jacobian matrix Drf (r0, Θ0) has an eigenvalue of zero real part. In our model, a saddle-node bifurcation (SN) occurred when Drf (r0, Θ0) had a single zero eigenvalue (in addition to some non-degenerate conditions), and a Hopf bifurcation (H) occurred when Drf (r0, Θ0) had a pair of purely imaginary eigenvalues. The bifurcation point was found numerically by XPPAUT or the Matlab toolbox MATCONT.

The number of parameters that must be varied simultaneously to evoke a bifurcation is defined as the codimension of this bifurcation (Guckenheimer and Holmes, 2013; Kuznetsov, 2013). Considering the infinite-dimensional space H of all vector fields defined on the n-dimensional Euclidean space ℝn, a vector field f0 undergoing a bifurcation, for example, a Hopf bifurcation, corresponds to a point in the space H. All nearby vector fields with the same singularity as f0 (i.e., vector fields that are orbitally topologically equivalent to f0) form a submanifold L of co-dimension k, which is an equivalence class of the singular vector field f0. Therefore, within the space H it requires another submanifold L~ of at least k dimensions to intersect transversely with L at point f0, such that the singularity of f0 persists under small perturbations of the vector field. The submanifold L~ was obtained through a parametrized family of vector fields involving at least k parameters. The least number k is then defined as the codimension of f0. The parametrized vector field f (r, Θ) in Equation (6) can be thought of as one realization of the submanifold L~ which passes through the vector field f0f (r0, Θ0) undergoing a bifurcation with the two parameters corresponding to the top-down weight w23 and the external stimulus μ.

Simulating Sniffing With a Periodically Driven Non-autonomous System

Olfactory sensation is an active process, with sensory stimuli being sampled by sniffing on the time scale of 215Hz in animal models (Carey and Wachowiak, 2011; Wachowiak, 2011). To simulate sniffing, a vector field for the dynamical system that depended explicitly on time and was also periodic with fixed period T = 2π / ω > 0, i.e.,

r˙  =  f(r, t) and f(r, t+T)=f(r, t)    (10)

could be rewritten in the form of an autonomous system by defining the function

φ(t)  =  ωt, mod 2π    (11)

such that using Equation (11), Equation (10) became

r˙  =  f(r, φ)φ˙=ω (r, φ)3×S1    (12)

where S1 denoted a circle. To construct the Poincaré map, we defined a cross-section of Equation (12) by

Σ0φ¯={(r, φ)3×S1 |φ=φ¯0(0, 2π] }    (13)

such that a fixed point of the Poincaré map Pφ¯0:Σφ¯0Σφ¯0 corresponded to a limit cycle of the extended system Equation (12), and a limit cycle of Pφ¯0 corresponded to a two-dimensional (2D) torus of Equation (12).

Topological changes in the ω-limit sets of the extended system Equation (12) could thus be understood via bifurcations of the discrete map Pφ¯0. Specifically, the bifurcation analysis we performed for autonomous system (6) also applied to the Poincaré map Pφ¯0. Hopf bifurcations undergone in autonomous system (6) which gave rise to limit cycles in 3D phase space corresponded to Neimark-Sacker bifurcations of Pφ¯0 which gave birth to a 2D torus in the extended space. The torus oscillation thus had two periodic components: one (the toroidal direction) driven extrinsically by the frequency of sniffs and the other (the poloidal direction) governed by the intrinsic network dynamics. Therefore, the Neimark-Sacker bifurcations provided an analogous bifurcation mechanism for non-autonomous system (10) as the Hopf bifurcations did for autonomous system (6).


Reduced Network Model Generates Complex Dynamics

To understand the functional role of top-down projections onto inhibitory neurons, we built a three-node network model (Figure 1A, see Methods) that recapitulated a circuit architecture identified both structurally (Padmanabhan et al., 2019) and functionally (Boyd et al., 2012; Markopoulos et al., 2012) across a number of brain areas. For different stimuli μ, the network exhibited a variety of dynamics (Figure 1B). For instance, when the stimulus was small, the firing rates ri, i = 1, 2, 3 had a fast-transient increase followed by damping oscillations that converged to a stationary state (Figure 1B, left). A sufficiently large stimulus μ elevated the firing rates to near saturation, where they then remained at the upper bound of the non-linear sigmoid function throughout the duration of the stimulus (Figure 1B, right). For small or large stimuli, the network responses converged to a constant firing rate after transient dynamics. By contrast, for medium values of μ, more complex firing rate dynamics emerged, including oscillations (Figure 1B, middle). To visualize the collective behaviors of M, G, and P populations to these different stimuli, we turned to a three dimensional dynamical system representation of the model where the time evolution of the firing rates (i.e., state variables) was a trajectory (or an orbit) in the phase space (r1, r2, r3) and the tangent vector defining the velocity of each point along a trajectory was given by the vector field f=(f1,f2,f3)T (see Methods) of Equation (1). The firing rates over time in Figure 1B thus corresponded to trajectories in Figure 1D starting from the origin O (where all three populations were silent). For small stimuli, the trajectory made an excursion before spiraling into an equilibrium indicated by the solid dot (Figure 1D, orange). Similarly, when the stimulus μ was large, the trajectory again settled into an equilibrium, but one that was translated within the phase space to the top-right corner (Figure 1D, black). Finally, for medium stimuli, the time-varying oscillation of firing rates manifested as a periodic orbit (or a limit cycle) in the 3D phase space (Figure 1D, brown). By convention, we defined the steady-state dynamics as the ω-limit set of the system.

Top-Down Weight Reshapes Network Dynamics and Modulates Neural Oscillations

Next, to explore how top-down down projections onto the inhibitory granule cell population (G) shaped the dynamics of the network, we studied the effects of changes in the connection weight w23 on firing rate dynamics. First, we varied the top-down weight w23 (Figure 2A, top) from the piriform population (P) to the inhibitory granule cell population (G) and studied the effect of these changes on the firing rate dynamics of the network. For a fixed stimulus (μ = 1.5) the dynamics of the firing rates ri (t), i = 1, 2, 3 were sensitive to different values of w23 (Figure 2A, bottom). When the top down weight was small (w23 = 4), firing rates approached the equilibrium exponentially (Figure 2A, bottom and Figure 2B, left, black traces). Conversely, when the top-down weight was large (w23 = 10.5), the firing rate of excitatory cells (r1) increased initially, but was suppressed as inhibition reduced the activity, until the firing rates ultimately settled to an equilibrium (Figure 2A, top and Figure 2B, left, light magenta traces). When the magnitude of top-down weight was changed to an intermediate value (w23 = 5.5), the same stimulus generated oscillatory activity in the network, with the steady-state dynamics transitioning to a periodic orbit (a limit cycle). Changing the weight of top-down projections onto the local inhibitory population for a single stimulus produced the same diversity of firing rate dynamics that occurred from changes in the stimulus. Furthermore, for a given top-down weight (w23), the effects on the network dynamics stimulus was unique to that stimulus (Figure 2B, right vs. left).

In regimes where specific weights of top-down weight generated sustained oscillatory activity for a given stimulus μ, we characterized the frequency and amplitude of these oscillations (Figure 2B, left w23 = 5.5, right w23 = 10.5) as changes in both have been tied to circuit function and behavior (Buzsaki and Draguhn, 2004; Kay et al., 2009). For a given stimulus, oscillations emerged between two values of w23, with the frequency of the oscillation varying monotonically (Figure 2C1). By contrast, while the amplitudes of the oscillations started from zero at the two critical values of w23, they reached a maximum in between (Figure 2D1). The control of both the frequency and amplitude via changes in w23 occurred across an array of weights (w31) associated with the feedforward drive from the mitral/tufted population (M) to the piriform population (P) (Figures 2C2,D2). Furthermore, the magnitude of synaptic weights from mitral/tufted cells to piriform cortical neurons established the dynamic range within which changes in top-down weights (w23) influenced the frequencies (Figure 2C3, 0 − 45Hz) and amplitudes (Figure 2D3, 0 − 1 A.U.) of network oscillations, spanning frequencies in the alpha, beta and gamma bands.

Top-Down Weight Contributes to Pattern Separation

As changing the top-down weight onto inhibitory neurons could generate complex activity patterns we next asked what computations could be performed by this control. For example, both behavioral and neurophysiological measures show that as the representations of two stimuli by neuronal circuits become different, distinguishing between them becomes easier (Friedrich and Laurent, 2001; Leutgeb et al., 2007; Yassa and Stark, 2011) Control of inhibition, via top-down centrifugal projections, may be one way that such stimulus discrimination is implemented by the circuit.

To test this hypothesis, we presented our network with a pair of stimuli, denoted by μ1 and μ2 (corresponding to stimuli arranged along a one-dimensional axis) and studied how control of inhibition altered the representations of the two stimuli (Figure 3A). Conceptually, these two stimuli could be two different concentrations of an odor or two odors that share a similar physiochemical feature (two odors with different carbon chain lengths). For a set of stimuli μi, i = 1, 2, we defined the steady-state representation of network activity as the ω-limit set Ωi, i = 1, 2. The distance between the two stimuli μ1 and μ2 in the stimulus space was defined as Δμ, and the resultant distance in the firing rate phase space between the two ω-limit sets (Ω1 and Ω2) we defined as a metric d (see Methods). The smaller the Δμ, the more similar the two stimuli were. We hypothesized that changes in the weight of feedback onto the inhibitory neuron population (w23) could increase the value of d, making the representations of those stimuli more distinct (Figure 3A).

In a representative example where μ1 = 1.0 and μ2 = 1.1, when the top-down weight was low (w23 = 3), the representations of the two stimuli were close (Figures 3B,C, left). As w23 was increased, the representations of the two stimuli were pushed apart making them more separable (w23 = 3.8, Figures 3B,C). Interestingly, as w23 was increased further (w23 = 10), the representations of the two stimuli became close to one another again (Figures 3B,C, right). As representations could be either oscillations in the state space, or equilibria, we compared how the distances of these representations changed across different measures (see methods). Interestingly, although the absolute values given by the Euclidean distance dE (Figure 3C) and the spectrum distance dS (Figure 3B) were different, they both occurred at the same feedback weight (Figure 3D). We visualized the distance landscapes defined by dE and dS over all combinations of μ1 and Δμ as a function of a change in the weight of the top-down weight (Figures 4A,B). Irrespective of which distance definition was exploited for measurement, we found an optimal value of w23max that maximized the distance between the two resultant representations for any given pair of stimuli.


Figure 4. Optimal top-down input for pattern separation. (A) Landscapes of the Euclidean distance dE over all stimulus pairs at three representative values of top-down input. (B) Landscapes of the spectrum distance dS unfolds as in (A). (C) The matrix w23max which maximizes d between network representations in response to all combinations of stimuli (μ1 and μ2) presented to the network. Inset: the same matrix of w23max organized by one stimulus μ1 vs. stimulus difference Δμ. (D) Dependence of the stationary firing rate of inhibitory population r2 on feedback w23 at different levels of input strength (indicated by color bar). (E) The correlation between w23max obtained from (C) for all pairs of stimuli μ1 and μ2 and the w23mid corresponding to the mid-point of the two inhibitory firing rate maxima associated with the same pair of stimuli (upper left inset) reveals that top-down input optimizes pattern separation by gating the G-M inhibition as well as recurrent G inhibition. The gray line denotes the utility line. Lower right inset shows the correlation coefficient and the slope of the linear regression.

A landscape of the optimal w23max across all pairs of stimuli (μ1, μ2) was shown in Figure 4C. Thus, changing the weight of top-down projections onto the inhibitory neuron population could be used to facilitate stimulus separation dynamically. To understand why, we examined the effect that varying the top-down weight had on the firing rate responses of both local excitatory mitral/tufted cell and inhibitory granule cell populations (r1 and r2). For a given stimulus, an increase in w23 led to a monotonic decrease in r1, suggesting persistent suppression onto the local M population. By contrast, the response of the local G population r2 was elevated first with increasing w23 until reaching the maximum r2max, after which r2 dropped significantly (Figure 4D). Across different stimuli μ, the shape of the firing rate r2 as a function of w23 remained the same but shifted vertically. To determine if these differences in the firing rate of inhibitory neurons (r2) were related to the values of top-down weights that maximally separated the distance between two stimulus representations, we plotted the w23max (abscissa) obtained from the landscape in Figure 4C vs. the midpoint w23mid (ordinate) between the r2max of the same stimulus pair (inset of Figure 4E). The response r2 to one stimulus on the left of the midpoint dropped significantly, while the response to a similar stimulus on the right still had a high inhibitory firing rate. The optimal w23max was correlated to w23mid (R2 = 0.98), the value at which inhibitory neuron activity from one stimulus was suppressed while activity from the other similar stimulus remained persistently high. Consequently, stimulus separation arose from the differential sensitivity of inhibitory neurons to the balance between top-down feedback and recurrent inhibition; an imbalance occurred between the top-down feedback and the recurrent inhibition for one stimulus while that balance was preserved for the second stimulus.

Top-Down Weight Contributes to Oscillation Synchrony

Stimulus-evoked oscillations also appeared in our model, and were modulated by the top-down weights (Figure 2) covering a wide range of frequency and amplitude. This suggested that oscillatory responses to different stimuli could be synchronized by tuning w23. To explore this, we first examined the oscillations in the firing rate generated by two different stimuli μ1 and μ2. At a given value of top-down weight (w23 = 12.0), one stimulus (μ1 = 1.15) generated oscillations (f1 = 33.5 Hz, Figure 5A1, before) in the piriform population's firing rate that were different in both frequency and amplitude from the oscillations (f2 = 37.Hz, Figure 5A1, before) in response to a second stimulus (μ2 = 0.6). However, a change in the top-down weight (w23 = 8), resulted in firing rate oscillations becoming more similar for the same two stimuli (Figure 5A1, after, f1 = 30.4 Hz, f2 = 30.8 Hz). This increase in the firing rate synchrony was also apparent when visualized in the 3D phase space (Figure 5A2). To quantify the synchrony between the oscillations responses to μ1 and μ2, we calculated the spectrum distance dS (see Methods) between the network representations for the two stimuli before and after changes in top-down weight (Figure 5B1). Changes in the feedback to inhibitory neurons w23 synchronized activity in the network stimuli (Figure 5B2), and while the effect was greatest when stimuli were similar, we found examples for stimuli that were initially as far apart as 20 Hz. As with stimulus discrimination, a systematic relationship emerged corresponding to the optimal top-down weight w23min across combinations of stimuli (μ1 vs. μ2) that was most effective at generating synchronous oscillations (Figure 5C).


Figure 5. Oscillation synchrony via top-down control. (A) Oscillatory responses to two example stimuli μ1 = 1.15 and μ2 = 0.6 become synchronized promptly after changing the top-down input w23. (A1) Time series of r3 (t) before and after changing the top-down input. (A2) Limit cycles corresponding to the oscillatory responses in (A1) are plotted in phase space (transitions not shown). (B) Changing w23 can make both the frequency and amplitude of two oscillations closer to each other. (B1) Frequency and amplitude components of the two oscillations shown in (A). (B2) Changing top-down input reduces the frequency differences of responses to two distinct stimuli, effectively using frequency to synchronize the representations in the phase space. (C) The matrix w23min which minimizes the distance d between oscillatory responses to all combinations of stimuli μ1 and μ2. (D) Schematic diagram illustrating that the same value of w23 which minimizes the distance between oscillations responding to stimuli μ1 and μ2 can maximize the distance between responses to stimuli μ1 and μ3. (E) Scatter plot in μ1- μ2- μ3 space where each sphere denotes a top-down input w23 as illustrated in (D) and is coded by color and size. (F) Correlation between the differences of those stimuli of which the response distances are simultaneously minimized and maximized. Inset: correlation coefficient and the slope of linear regression.

Although we have thus far treated stimulus discrimination and synchrony separately, neural circuits perform both operations simultaneously, bringing the network representation of one stimulus closer to another, while simultaneously pushing the representation of that stimulus farther from a third. We therefore tested if a single change in the top-down weight w23 accomplish both of these operations; minimize the distance between the responses to one pair of stimuli (μ1 vs. μ2) while also maximize the response distance to another other pair of stimuli (μ1 vs. μ3, Figure 5D). To do this, we generated a 3D scatter plot of values of w23 that were optimal for synchrony between oscillations generated by stimulus μ1 and μ2 (Figure 5C) and also produced a maximum separation between the representations of stimulus μ1 and μ3 (Figure 5E). The values of top-down weight w23 for each point that fulfilled these diametrically distinct functions were coded by color and size (Figure 5E). We found that the top-down weight corresponding to both operations scaled with the stimuli, such that when μ1, μ2 and μ3 were small, the top-down weight was also small, but as the three stimuli increased in magnitude, the top-down weight needed to synchronize one pair and separate the other pair also increased. Finally, we found a strong correlation between the values of stimulus differences: |μ1 − μ2| and |μ1 − μ3| (Figure 5F) at which an w23 weight was optimal for stimulus separation and oscillatory synchrony.

Generalization to Oscillatory Stimulus Driven by Sniffs

Although we used a constant stimulus μ to represent the average input to mitral/tufted cells, in mammals sniffing brings odors into the nasal epithelium in a periodic fashion (Wachowiak, 2011). Sniff cycles carry different amounts of information about odor identity and concentration (Miura et al., 2012) and a single sniff cycle is sufficient for animals to discriminate accurately between two odors (Uchida and Mainen, 2003; Wesson et al., 2008). To explore how changing top-down weights can reshape network responses to oscillatory stimuli, we modeled our stimulus μ as a sinusoidal function μ~(t)=μcos(ωst+φ0), where the different odors had different amplitudes μ, the sniffing frequency ωs was set to ~4 Hz and φ0 characterized the initial phase of sniffing (Carey and Wachowiak, 2011; Shusterman et al., 2011) (Figure 6A).


Figure 6. Pattern separation and oscillation synchrony for sniff-modulated oscillatory stimulus. (A) Schematic diagram illustrating that varying the top-down weight of the network model (middle) can accomplish both pattern separation and oscillation synchrony (bottom) for a pair of oscillatory stimuli μ1 and μ2 modulated by sniffs (top). (B) Pattern separation for two oscillatory stimuli with closely related amplitudes: μ1 = 1.4 and μ2 = 1.5. (B1) Firing rate of the piriform population r3 (t) before and after changing the top-down weight: before, w23 = 26; after, w23 = 21. (B2) The limit cycles in the phase space corresponding to the firing rate in (B1). (C) Oscillation synchrony for two oscillatory stimuli with distinct amplitudes: μ1 = 1.15 and μ2 = 0.6. (C1) Changing the top-down weight from w23 = 6.1 to w23 = 9.2 synchronizes the intrinsic oscillations of the firing rate of the piriform population r3 (t). (C2) Frequency and amplitude components of the two intrinsic oscillations shown in (C1).

For a pair of oscillatory stimuli μ~(t) with two similar amplitudes μ1 and μ2, the firing rate responses were also similar for the piriform population (P) (Figure 6B1, before top-down change) and the entire network in the phase space (Figure 6B2, left). If we changed the top-down weight w23, as we had done for a fixed stimulus, both the piriform population firing (Figure 6B1, after top-down change) and the network representations became more distinct (Figure 6B2, right). The conch-shaped limit cycle in Figure 6B2 (right) arose from oscillations occurring at two different time scales (see methods): a slower oscillation on the time scale of sniff cycles and a faster oscillation governed by the intrinsic dynamics of the network. As a consequence, the same pattern separation achieved by changing top-down weight for constant stimuli could also be accomplished for oscillatory stimuli.

To explore if network representations of oscillatory stimuli could be made synchronous by changing top-down weights, we presented two stimuli with distinct amplitudes (μ1, μ2) to the network (Figure 6A). Following a change in the top-down weight w23, the network representations became synchronous (Figure 6C1), with the oscillations of firing rates occurring at the same frequency (Figure 6C2, right). Importantly, these high frequency oscillations occurred at the gamma band, and rode on top of the slower oscillations corresponding to sniff cycles, further revealing the computational decoupling of sniffing and inhibitory dynamics across two different time scales. Taken together, the mechanisms giving rise to both pattern separation and oscillatory synchrony were general to constant and oscillatory inputs.

Bifurcation Mechanism for Top-Down Control of Inhibition

Finally, to understand mathematically how such operations emerged from changes in the top-down weight to inhibition, we studied the structure of the transitions in network firing rates dynamics (Figure 2A). These transitions were associated with qualitative or topological changes in the ω-limit sets of the system, indicative of the occurrence of bifurcations in the system.

To explore this further, we first examined the how the ω-limit sets of the system receiving constant stimuli changed with different top-down weights. An equilibrium corresponding to constant firing rates in the network arose from the intersection of three nullclines (Figure 7, yellow thick lines), each resulting from pairwise intersections of three “nullplanes” that characterized the geometric surface on which the firing rate derivatives of one node equaled to zero (Figure 7, transparent surfaces). Global phase structures for two representative values of top-down weights (Figure 7A: ω23 = 4, Figure 7C: ω23 = 15) illustrated how these equilibria varied within the firing rate phase space. In these two examples, both equilibria were stable and attractive, with all nearby trajectories (Figures 7A,C, black thin lines) moving toward them. This was, however, not true for all values of w23. At some critical values of w23, the equilibrium lost stability, and a small-amplitude limit cycle branched from that unstable equilibrium, resulting in the oscillations observed in the dynamics (Figure 7B). This transition signified a Hopf bifurcation of the system (see Methods), which arose when the top-down weight w23 was within a specific regime. Therefore, across all combinations of external stimuli μ and top-down weight w23, we obtained a smooth manifold in the phase space (Figure 8A, left), corresponding to a family of ω-limit sets on which the network dynamics settled from any set of initial conditions. Sustained oscillations corresponded to the red region-LC (LC: limit cycle) where each equilibrium (unstable) was paired with exactly one limit cycle born simultaneously via a Hopf bifurcation (the purple empty square vs. the dot-dashed curve). Constant firing rates corresponded to the gray region-EE (EE: exponential equilibrium) and green region-SE (SE: spiraling equilibrium), where the equilibria were stable, approached either exponentially (region-EE) or via damping-oscillations (region-SE). Finally, the two blue regions on the manifold were bounded by saddle-node bifurcations near Bogdanov–Takens (BT), an example global phase structure of which was shown in Figure 7D. The equilibrium manifold thus defined the entire family of network representations for all possible combinations of stimuli μ and top-down weigh tw23.


Figure 7. Global phase structure changes with the top-down input. (A–C) Global phase structures showing the nullclines (yellow thick curves), nullplanes (transparent surfaces with the same color code as M, G, P population) and several representative trajectories (black thin curves) for the same stimulus μ = 1.5 and three different values of w23. A, w23 = 4; B, w23 = 6; C, w23 = 15. Varying the top-down input tilts the nullplanes, thus changing the position of the equilibrium as well as its stability. (D) An example phase structure where three equilibria were present simultaneously (two stable and one unstable), corresponding one of the blue regions of the manifold in Figure 8A.


Figure 8. Bifurcation mechanism of top-down control to support both pattern separation and oscillation synchrony. (A) Illustration of the bifurcation mechanism and the transition boundary of dynamics in phase space and parameter space. (A1) Pattern separation. Left: the equilibrium manifold in the firing rate phase space divided by the separatrix emitting from multiple codimension-2 bifurcation points: BT (Bogdanov–Takens bifurcation) into several regions. Region-LC: each equilibrium was unstable and had exactly one corresponding stable limit cycle (dot dashed cycle) arising from a Hopf bifurcation. Region-SE: each equilibrium was stable all trajectories spiraled into it. Region-EE: each equilibrium was stable all trajectories approached it exponentially. Blue regions: regions where multiple equilibria coexisted. For two example stimuli μ1 = 2.0 and μ2 = 3.0 given in the middle of (A1), two paths of equilibria were induced on the equilibrium manifold and traversed across different regions as changing top-down input w23. Middle: different regions on the equilibrium manifold corresponded to different regimes in the parameter space of μ and w23 in the same color scheme [parameters for blue regions in Left were largely beyond the range thus not shown]. The transition boundary ΓH specified the pair (μ, w23) at which the network underwent a Hopf bifurcation and corresponded to the separatrix enclosing the region-LC in Left. Two given stimuli were denoted by two vertical lines and three example values of w23 corresponded to three horizontal dashed lines, giving rise to a pair of junctions for each. These junctions were also plotted as squares in the left of (A1) denoting the corresponding ω-limit sets in the same color (solid square: stable equilibrium; empty square: unstable equilibrium). Right: the distance between the ω-limit sets to represent the two given stimuli. The maximal distance was achieved when the two junctions were on opposite sides of ΓH. (A2) Same as (A1) but for oscillation synchrony occurring when the two junctions were both inside ΓH. (B) Comparisons between the translated transition boundary ΓΔ (dashed curve) depending on μ (left: Δμ = 0.1, middle: Δμ = 0.5) and the sliced section of the w23max at the same Γμ (solid curve). Right: a series of translated ΓΔ for three representative values of Γμ.

Within the manifold of the stimulus μ and top-down weight w23, we identified a transition boundary ΓH (black solid curves, Figure 8A, middle) corresponding to the separatrix enclosing the region-LC. ΓH specified the parameter pairs (μ, w23) at which a Hopf bifurcation occurred, thereby dividing the parameter space into regimes with different dynamics (same color coded as Figure 8A1, left). For a given pair of stimuli (for example, μ1 = 2.0, μ2 = 3.0), changing w23 corresponded to shifting the horizontal dashed line vertically (three representatives were shown in Figure 8A1, middle), thereby shifting the junctions with the two stimuli (vertical solid lines, Figure 8A1, middle) across different regimes in the parameter space. In the firing rate phase space (Figure 8A, left), these changes in w23 for one stimulus moved the equilibrium through different regions of the manifold: EE-LC-SE, while for another stimulus, a parallel curve on the manifold could also be traced. When the two junctions in Figure 8A1 (middle) were on different sides of the transition boundary, with one equilibrium in region-LC and the other in region-SE (Figure 8A1, left), the two network representations became topologically different from each other; the former a limit cycle, and the latter an equilibrium point. Thus for a combination of stimulus pairs, the optimal w23max for pattern separation was then achieved when the ω-limit sets were on different sides of transition boundary (Figure 8A1, right).

Furthermore, when the junction of feedback weight and stimulus pair were both inside the transition boundary (Figure 8A2, middle) two limit cycles emerged (one for each stimulus, for example, μ1 = 3.45, μ2 = 3.95, Figure 8A2) synchronize the network representations. Changes in top-down weight moved the junctions for pairs of stimuli within the parameter space, revealing a shared mechanism supported both stimulus separation and oscillation synchrony, depending on the relative positions of the junctions with respect to the transition boundary.

Finally we determined if the transition boundary identified via analysis of the dynamical system corresponded to the w23max matrix found in Figure 4Ci. To do this, we considered a set of initial stimuli μ1, and a set of distances to a second set of stimuli Δμ, wherein each value was an array that defined a set of stimulus pairs {(μ1, μ1 + Δμ)| μ1 ∈ [0, 4]}. For a given Δμ > 0, distinguishing the pair (μ1, μ1 + Δμ) was the same as distinguishing (μ1 + μ, Δμ1) in terms of pattern separation. In this framework two different distances, for instance, Δμ = 0.1 or μ = 0.5 (Figure 8B, left and middle), the set of stimulus pairs had a unique transition boundary ΓΔ (Figure 8B, right). The section of the w23max matrix in Figure 4 Ci for the set of stimulus pairs {(μ1, μ1 + Δμ)| μ1 ∈ [0, 4], Δμ is given} was correlated with the transition boundary ΓΔ of the same value μ (Figure 8B). For small Δμ = 0.1, the slice of the w23max followed closely with the transition boundary ΓΔ (Figure 8B, left). As the stimulus difference increased (Figure 8B, middle), w23max deviated from the boundary ΓΔ. Larger Δμ's increased the bifurcation lag between two stimuli such that the stimulus that caused a bifurcation first had more parameter space to develop before the bifurcation of the other stimulus. Conversely, stimulus discrimination was harder as μ decreased because the range of top-down weights w23 that separate two stimuli shrank significantly around a close vicinity of the transition boundary. Thus, subtle adjustments of top-down weight around the transition boundary were required to separate similar stimuli from each other. The same analysis could also be performed for the non-autonomous system receiving oscillatory stimuli μ~(t) by investigating bifurcations of fixed points of the constructed Poincaré map on a given cross section (see Methods) with the same computational mechanism arising via the discrete version of Hopf bifurcation, i.e., a Neimark-Sacker bifurcation (Kuznetsov, 2013). Taken together, these results provide a bridge linking the mechanisms that give rise to the dynamics of the neural circuit with the computations performed by the circuit.


Using a three-node model, which included top-down projections from piriform cortical cells onto inhibitory granule cells in the main olfactory bulb, we identified a network capable of complex dynamic behaviors, ranging from an attractor to stable oscillations across a range of frequencies and amplitudes. By changing the weight of these top-down projections, the network could either facilitate pattern separation between two similar stimuli, or synchronize the oscillatory activity produced by two different stimuli. A bifurcation analysis of the dynamical system revealed that both mechanisms emerged from the transition boundary of Hopf bifurcations which branched from co-dimensional two bifurcation points (i.e., the Bogdanov-Takens bifurcation). Furthermore, these computations could be accomplished even when the stimuli were periodic, fluctuating at the frequency of sniffing (Neimark-Sacker bifurcation), suggesting that these findings are a general feature of this network. Our results provide both a mathematical framework for how top-down control of inhibition shapes the dynamics of a network, and a link between such dynamics and the computations that neural circuits can perform.

An important point to consider is how changes in top-down weights may be instantiated biologically? This point depends on the timescale of weight changes. On short time scales, changes in inhibitory drive to granule cells can facilitate olfactory discrimination (Abraham et al., 2010; Nunes and Kuner, 2015) and generate synchronous oscillatory activity among mitral cells in the bulb (Galan et al., 2006). Neuromodulators such as serotonin (Petzold et al., 2009; Kapoor et al., 2016) can act on fast sub second time scales to support both oscillatory synchrony and stimulus discrimination, providing one biological mechanism by which weights can be changed dynamically. By contrast, long-term changes in the bulb may be instantiated by classis synaptic plasticity mechanisms such as LTP (Cauthron and Stripling, 2014), or via the remodeling of synaptic connectivity (Arenkiel et al., 2011; Deshpande et al., 2013), for instance due to adult neurogenesis (Lledo et al., 2006). In these examples, the changes in feedback weight likely reflect slow alterations in network structure that result in stable changes in neural representations, possibly corresponding to learning.

While the biological mechanisms by which the top-down synaptic weights change onto inhibitory neurons may be diverse depending on timescale, we find that such alterations give rise to functionally equivalent changes supporting an array of computations. For instance, changes in the top-down weight would render two stimuli more distinct at the level of firing rates in the population, a process referred to as pattern separation (Cayco-Gajic and Silver, 2019) or decorrelation (Friedrich and Laurent, 2001). Our model predicts that pattern separation arises from the non-monotonic change in firing in granule cells (at the balance between op-down excitation and recurrent inhibition). The top-down weight onto inhibitory neurons sets a gate, allowing some stimuli to cross a threshold of recurrent inhibition, while others do not.

In parallel, changing top-down weights onto inhibitory neurons can increase the synchrony between two stimuli that were initially asynchronous. A number of experimental and theoretical studies have explored the privileged role that inhibitory interneurons play in generating gamma oscillations (Whittington et al., 1995; Hasenstaub et al., 2005; Cardin et al., 2009; Sohal et al., 2009; Tiesinga and Sejnowski, 2009). Among these, the two most common models are when gamma arises from reciprocal coupling between pyramidal cells and inhibitory interneurons (PING), and recurrent connections among inhibitory interneurons (ING) (Whittington et al., 2000; Tiesinga and Sejnowski, 2009). In both, oscillatory activity arises from the structure of local connectivity. In our work, we identified another motif by which gamma oscillations can arise—Top-down control of Inhibitory Neuron Gamma (TING). Local excitatory mitral and tufted cells broadcast activity patterns to a pyramidal/semilunar cell population in piriform cortex, that then synapses back onto inhibitory granule cells.

Studies on dynamics of local excitatory and inhibitory neurons in the olfactory system both experimentally and mathematically are extensive (Wilson and Cowan, 1972; Ermentrout and Kopell, 1990; Kay et al., 2009; Li and Cleland, 2017). To these models we add a description of how an external (in this case, top-down input from piriform cortex) source controlling the inhibitory neuron population can influence dynamics. In studying the dynamical system defined by this network, we found that the bifurcations largely result from the singularity (linearized Jacobian matrix is non-hyperbolic) inherently embedded in the system itself. Thus, although the exact parameter values (defined by the weights of connections) influence when the dynamics of the network undergoes a bifurcation, the types of bifurcations that arise are determined by the normal form of the system (Guckenheimer and Holmes, 2013; Kuznetsov, 2013); revealing that the behaviors observed in this three node population are a fundamental feature of the network architecture. For Wilson-Cowan equations we used, the Bogdanov–Takens bifurcation is the inherent codimension-2 singularity (Cowan et al., 2016), meaning that the diversity of dynamics exists for a broad range of parameter settings, and that the unfolding of these dynamics can be implemented by modulating the top-down connection weight. Our model address this in the context of olfaction (Oswald and Urban, 2012a), but it may be applicable to a number of other sensory systems that share a similar architecture. For instance, the axonal projections from the cingulate of frontal cortex to GABAergic inhibitory neurons in V1 of the mouse visual system are organized (Zhang et al., 2014), and may therefore serve an analogous function as piriform projections to granule cells. Consequently, we identified a generalized principle by which control of inhibition via top-down weights can support a number of computations essential for neural circuit function.

Finally, we found that the firing rate representations of mitral/tufted cells, granule cells, and piriform neurons resided within distinct domains on a manifold defined by the stimulus and the weight of feedback. These domains corresponded to transitions in the dynamics of the system. Changes in the top-down weights moved a transition boundary that delineating these domains across different stimuli. When two stimuli were on opposites sides of this transition boundary, their dynamics operated under two different regimes, and their representations were pushed further apart. By contrast, when the stimuli were both on the same side of the transition boundary, within regimes corresponding to similar dynamics, their activity became more synchronous; effectively binding those stimuli together. Changes in top-down weight were therefore changes in the location of the transition boundary that could either marshal the representations of two stimuli together or push them apart. In conclusion, we identified a model that links the dynamics of neural systems with the computations they are hypothesized to perform and may be used as a generalized framework to study the diverse effects of feedback onto inhibitory populations.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

KP conceived and supervised the project. ZC performed all the experiments and analysis. ZC and KP made the figures and wrote the manuscript. All authors contributed to the article and approved the submitted version.


This study was supported by funding from the National Institutes of Health (NIH) and the National Science Foundation (NSF). KP was funded by NIH R01 MH113924, NSF CAREER 1749772, the Cystinosis Research Foundation, and the Kilian J. and Caroline F. Schmitt Foundation. This manuscript has been released as a pre-print at (Chen and Padmanabhan, 2020).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Abraham, N. M., Egger, V., Shimshek, D. R., Renden, R., Fukunaga, I., Sprengel, R., et al. (2010). Synaptic inhibition in the olfactory bulb accelerates odor discrimination in mice. Neuron 65, 399–411. doi: 10.1016/j.neuron.2010.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Arenkiel, B. R., Hasegawa, H., Yi, J. J., Larsen, R. S., Wallace, M. L., Philpot, B. D., et al. (2011). Activity-induced remodeling of olfactory bulb microcircuits revealed by monosynaptic tracing. PLoS ONE 6:e29423. doi: 10.1371/journal.pone.0029423

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyd, A. M., Sturgill, J. F., Poo, C., and Isaacson, J. S. (2012). Cortical feedback control of olfactory bulb circuits. Neuron 76, 1161–1174. doi: 10.1016/j.neuron.2012.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Brea, J. N., Kay, L. M., and Kopell, N. J. (2009). Biophysical model for gamma rhythms in the olfactory bulb via subthreshold oscillations. Proc. Natl. Acad. Sci. U.S.A. 106, 21954–21959. doi: 10.1073/pnas.0910964106

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardin, J. A., Carlén, M., Meletis, K., Knoblich, U., Zhang, F., Deisseroth, K., et al. (2009). Driving fast-spiking cells induces gamma rhythm and controls sensory responses. Nature 459, 663–667. doi: 10.1038/nature08002

PubMed Abstract | CrossRef Full Text | Google Scholar

Carey, R. M., and Wachowiak, M. (2011). Effect of sniffing on the temporal structure of mitral/tufted cell output from the olfactory bulb. J. Neurosci. 31, 10615–10626. doi: 10.1523/JNEUROSCI.1805-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cauthron, J. L., and Stripling, J. S. (2014). Long-term plasticity in the regulation of olfactory bulb activity by centrifugal fibers from piriform cortex. J. Neurosci. 34, 9677–9687. doi: 10.1523/JNEUROSCI.1314-14.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Cayco-Gajic, N. A., and Silver, R. A. (2019). Re-evaluating circuit mechanisms underlying pattern separation. Neuron 101, 584–602. doi: 10.1016/j.neuron.2019.01.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., and Padmanabhan, K. (2020). Top-down control of inhibition reshapes neural dynamics giving rise to a diversity of computations. bioRxiv [Preprint]. doi: 10.1101/2020.02.25.964965

CrossRef Full Text | Google Scholar

Cleland, T. A., and Linster, C. (2005). Computation in the olfactory system. Chem. Sens. 30, 801–813. doi: 10.1093/chemse/bji072

PubMed Abstract | CrossRef Full Text | Google Scholar

Cowan, J. D., Neuman, J., and van Drongelen, W. (2016). Wilson–cowan equations for neocortical dynamics. J. Math. Neurosci. 6, 1–24. doi: 10.1186/s13408-015-0034-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Davison, I. G., and Ehlers, M. D. (2011). Neural circuit mechanisms for pattern detection and feature combination in olfactory cortex. Neuron 70, 82–94. doi: 10.1016/j.neuron.2011.02.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Deshpande, A., Bergami, M., Ghanem, A., Conzelmann, K.-K., Lepier, A., Götz, M., et al. (2013). Retrograde monosynaptic tracing reveals the temporal evolution of inputs onto new neurons in the adult dentate gyrus and olfactory bulb. Proc. Natl. Acad. Sci. U.S.A. 110, E1152–E1161. doi: 10.1073/pnas.1218991110

PubMed Abstract | CrossRef Full Text | Google Scholar

Egger, V., Svoboda, K., and Mainen, Z. F. (2005). Dendrodendritic synaptic signals in olfactory bulb granule cells: local spine boost and global low-threshold spike. J. Neurosci. 25, 3521–3530. doi: 10.1523/JNEUROSCI.4746-04.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1990). Oscillator death in systems of coupled neural oscillators. SIAM J. Appl. Math. 50, 125–146. doi: 10.1137/0150009

CrossRef Full Text | Google Scholar

Franci, A., Herrera-Valdez, M. A., Lara-Aparicio, M., and Padilla-Longoria, P. (2018). Synchronization, oscillator death, and frequency modulation in a class of biologically inspired coupled oscillators. Front. Appl. Math. Stat. 4:51. doi: 10.3389/fams.2018.00051

CrossRef Full Text | Google Scholar

Franks, K. M., and Isaacson, J. S. (2006). Strong single-fiber sensory inputs to olfactory cortex: implications for olfactory coding. Neuron 49, 357–363. doi: 10.1016/j.neuron.2005.12.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedrich, R. W., and Laurent, G. (2001). Dynamic optimization of odor representations by slow temporal patterning of mitral cell activity. Science 291, 889–894. doi: 10.1126/science.291.5505.889

PubMed Abstract | CrossRef Full Text | Google Scholar

Galan, R. F., Fourcaud-Trocme, N., Ermentrout, G. B., and Urban, N. N. (2006). Correlation-induced synchronization of oscillations in olfactory bulb neurons. J. Neurosci. 26, 3646–3655. doi: 10.1523/JNEUROSCI.4605-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Guckenheimer, J., and Holmes, P. (2013). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York, NY: Springer Science & Business Media.

Google Scholar

Hasenstaub, A., Shu, Y., Haider, B., Kraushaar, U., Duque, A., and McCormick, D. A. (2005). Inhibitory postsynaptic potentials carry synchronized frequency information in active cortical networks. Neuron 47, 423–435. doi: 10.1016/j.neuron.2005.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapoor, V., Provost, A. C., Agarwal, P., and Murthy, V. N. (2016). Activation of raphe nuclei triggers rapid and distinct effects on parallel olfactory bulb output channels. Nat. Neurosci. 19, 271–282. doi: 10.1038/nn.4219

PubMed Abstract | CrossRef Full Text | Google Scholar

Kapoor, V., and Urban, N. N. (2006). Glomerulus-specific, long-latency activity in the olfactory bulb granule cell network. J. Neurosci. 26, 11709–11719. doi: 10.1523/JNEUROSCI.3371-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kay, L. M., Beshel, J., Brea, J., Martin, C., Rojas-Libano, D., and Kopell, N. (2009). Olfactory oscillations: the what, how and what for. Trends Neurosci. 32, 207–214. doi: 10.1016/j.tins.2008.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuznetsov, Y. A. (2013). Elements of Applied Bifurcation Theory. New York, NY: Springer Science & Business Media.

Google Scholar

Ledoux, E., and Brunel, N. (2011). Dynamics of networks of excitatory and inhibitory neurons in response to time-dependent inputs. Front. Comput. Neurosci. 5:25. doi: 10.3389/fncom.2011.00025

PubMed Abstract | CrossRef Full Text | Google Scholar

Leutgeb, J. K., Leutgeb, S., Moser, M.-B., and Moser, E. I. (2007). Pattern separation in the dentate gyrus and CA3 of the hippocampus. Science 315, 961–966. doi: 10.1126/science.1135801

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G., and Cleland, T. A. (2017). A coupled-oscillator model of olfactory bulb gamma oscillations. PLoS Comput. Biol. 13:e1005760. doi: 10.1371/journal.pcbi.1005760

PubMed Abstract | CrossRef Full Text | Google Scholar

Lledo, P.-M., Alonso, M., and Grubb, M. S. (2006). Adult neurogenesis and functional plasticity in neuronal circuits. Nat. Rev. Neurosci. 7, 179–193. doi: 10.1038/nrn1867

PubMed Abstract | CrossRef Full Text | Google Scholar

Markopoulos, F., Rokni, D., Gire, D. H., and Murthy, V. N. (2012). Functional properties of cortical feedback projections to the olfactory bulb. Neuron 76, 1175–1188. doi: 10.1016/j.neuron.2012.10.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Miura, K., Mainen, Z. F., and Uchida, N. (2012). Odor representations in olfactory cortex: distributed rate coding and decorrelated population activity. Neuron 74, 1087–1098. doi: 10.1016/j.neuron.2012.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyamichi, K., Amat, F., Moussavi, F., Wang, C., Wickersham, I., Wall, N. R., et al. (2011). Cortical representations of olfactory input by trans-synaptic tracing. Nature 472, 191–196. doi: 10.1038/nature09714

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunes, D., and Kuner, T. (2015). Disinhibition of olfactory bulb granule cells accelerates odour discrimination in mice. Nat. Commun. 6:8950. doi: 10.1038/ncomms9950

PubMed Abstract | CrossRef Full Text

Nusser, Z., Kay, L. M., Laurent, G., Homanics, G. E., and Mody, I. (2001). Disruption of GABAA receptors on GABAergic interneurons leads to increased oscillatory power in the olfactory bulb network. J. Neurophysiol. 86, 2823–2833. doi: 10.1152/jn.2001.86.6.2823

PubMed Abstract | CrossRef Full Text | Google Scholar

Oswald, A. M., and Urban, N. N. (2012a). There and back again: the corticobulbar loop. Neuron 76, 1045–1047. doi: 10.1016/j.neuron.2012.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Oswald, A. M. M., and Urban, N. N. (2012b). Interactions between behaviorally relevant rhythms and synaptic plasticity alter coding in the piriform cortex. J. Neurosci. 32, 6092–6104. doi: 10.1523/JNEUROSCI.6285-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Padmanabhan, K., Osakada, F., Tarabrina, A., Kizer, E., Callaway, E. M., Gage, F. H., et al. (2016). Diverse representations of olfactory information in centrifugal feedback projections. J. Neurosci. 36, 7535–7545. doi: 10.1523/JNEUROSCI.3358-15.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Padmanabhan, K., Osakada, F., Tarabrina, A., Kizer, E., Callaway, E. M., Gage, F. H., et al. (2019). Centrifugal inputs to the main olfactory bulb revealed through whole brain circuit-mapping. Front. Neuroanat. 12:115. doi: 10.3389/fnana.2018.00115

PubMed Abstract | CrossRef Full Text | Google Scholar

Petzold, G. C., Hagiwara, A., and Murthy, V. N. (2009). Serotonergic modulation of odor input to the mammalian olfactory bulb. Nat. Neurosci. 12, 784–791. doi: 10.1038/nn.2335

PubMed Abstract | CrossRef Full Text | Google Scholar

Shipley, M. T., and Adamek, G. D. (1984). The connections of the mouse olfactory bulb: a study using orthograde and retrograde transport of wheat germ agglutinin conjugated to horseradish peroxidase. Brain Res. Bull. 12, 669–688. doi: 10.1016/0361-9230(84)90148-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Shusterman, R., Smear, M. C., Koulakov, A. A., and Rinberg, D. (2011). Precise olfactory responses tile the sniff cycle. Nat. Neurosci. 14:1039. doi: 10.1038/nn.2877

PubMed Abstract | CrossRef Full Text | Google Scholar

Sohal, V. S., Zhang, F., Yizhar, O., and Deisseroth, K. (2009). Parvalbumin neurons and gamma rhythms enhance cortical circuit performance. Nature 459, 698–702. doi: 10.1038/nature07991

PubMed Abstract | CrossRef Full Text | Google Scholar

Stettler, D. D., and Axel, R. (2009). Representations of odor in the piriform cortex. Neuron 63, 854–864. doi: 10.1016/j.neuron.2009.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, N., and Bekkers, J. M. (2011). Two layers of synaptic processing by principal neurons in piriform cortex. J. Neurosci. 31, 2156–2166. doi: 10.1523/JNEUROSCI.5430-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiesinga, P., and Sejnowski, T. J. (2009). Cortical enlightenment: are attentional gamma oscillations driven by ING or PING? Neuron 63, 727–732. doi: 10.1016/j.neuron.2009.09.009

CrossRef Full Text | Google Scholar

Tsodyks, M. V., Skaggs, W. E., Sejnowski, T. J., and McNaughton, B. L. (1997). Paradoxical effects of external modulation of inhibitory interneurons. J. Neurosci. 17, 4382–4388. doi: 10.1523/JNEUROSCI.17-11-04382.1997

PubMed Abstract | CrossRef Full Text | Google Scholar

Uchida, N., and Mainen, Z. F. (2003). Speed and accuracy of olfactory discrimination in the rat. Nat. Neurosci. 6, 1224–1229. doi: 10.1038/nn1142

PubMed Abstract | CrossRef Full Text | Google Scholar

Urban, N. N., and Sakmann, B. (2002). Reciprocal intraglomerular excitation and intra- and interglomerular lateral inhibition between mouse olfactory bulb mitral cells. J. Physiol. 542, 355–367. doi: 10.1113/jphysiol.2001.013491

PubMed Abstract | CrossRef Full Text | Google Scholar

Veltz, R., and Sejnowski, T. J. (2015). Periodic forcing of inhibition-stabilized networks: nonlinear resonances and phase-amplitude coupling. Neural. Comput. 27, 2477–2509. doi: 10.1162/NECO_a_00786

PubMed Abstract | CrossRef Full Text | Google Scholar

Wachowiak, M. (2011). All in a sniff: olfaction as a model for active sensing. Neuron 71, 962–973. doi: 10.1016/j.neuron.2011.08.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Wesson, D. W., Carey, R. M., Verhagen, J. V., and Wachowiak, M. (2008). Rapid encoding and perception of novel odors in the rat. PLoS Biol. 6:82. doi: 10.1371/journal.pbio.0060082

PubMed Abstract | CrossRef Full Text | Google Scholar

Whittington, M. A., Traub, R., Kopell, N., Ermentrout, B., and Buhl, E. (2000). Inhibition-based rhythms: experimental and mathematical observations on network dynamics. Int. J. Psychophysiol. 38, 315–336. doi: 10.1016/S0167-8760(00)00173-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Whittington, M. A., Traub, R. D., and Jefferys, J. G. (1995). Synchronized oscillations in interneuron networks driven by metabotropic glutamate receptor activation. Nature 373, 612–615. doi: 10.1038/373612a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Willhite, D. C., Nguyen, K. T., Masurkar, A. V., Greer, C. A., Shepherd, G. M., and Chen, W. R. (2006). Viral tracing identifies distributed columnar organization in the olfactory bulb. Proc. Natl Acad. Sci. U.S.A. 103, 12592–12597. doi: 10.1073/pnas.0602032103

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Yassa, M. A., and Stark, C. E. (2011). Pattern separation in the hippocampus. Trends Neurosci. 34, 515–525. doi: 10.1016/j.tins.2011.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Xu, M., Kamigaki, T., Do, J. P. H., Chang, W.-C., Jenvay, S., et al. (2014). Long-range and local circuits for top-down modulation of visual cortex processing. Science 345, 660–665. doi: 10.1126/science.1254126

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: dynamical system, olfactory bulb, oscillations, pattern separation, synchrony, bifurcation, feedback, top-down

Citation: Chen Z and Padmanabhan K (2020) Top-Down Control of Inhibitory Granule Cells in the Main Olfactory Bulb Reshapes Neural Dynamics Giving Rise to a Diversity of Computations. Front. Comput. Neurosci. 14:59. doi: 10.3389/fncom.2020.00059

Received: 31 March 2020; Accepted: 22 May 2020;
Published: 13 July 2020.

Edited by:

Daya Shankar Gupta, Camden County College, United States

Reviewed by:

Markus Rothermel, RWTH Aachen University, Germany
Christiane Linster, Cornell University, United States

Copyright © 2020 Chen and Padmanabhan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Krishnan Padmanabhan,

ORCID: Zhen Chen
Krishnan Padmanabhan