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

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.


INTRODUCTION
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; 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 excitatoryinhibitory (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.

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, r i (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:    τ 1ṙ1 = −r 1 + S (w 11 r 1 + w 12 r 2 + µ) τ 2ṙ2 = −r 2 + S (w 21 r 1 + w 22 r 2 + w 23 r 3 ) τ 3ṙ3 = −r 3 + tanh (w 31 r 1 + w 33 r 3 ) (1) where S is the sigmoid function: 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 w ij , i, j = 1, 2, 3, among which w 11 , w 21 , w 31 , w 23 , w 33 > 0 and w 12 , w 22 < 0. The connection weight w ij , 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: w 11 = 8.7, w 12 = −10, w 21 = 7.0, w 22 = −13, w 31 = 1.5, w 33 = 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 (w 11 ) was determined by mapping recorded excitatory postsynaptic potentials (EPSP) in M/T cells and the weight of inhibition from granule cells (w 12 ) 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, |w 12 | > |w 21 |. Estimates of synaptic weights for the projections to piriform cortex and the feedback connections were estimated from Franks and Isaacson (2006), Suzuki andBekkers (2011), andBoyd 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 w 31 r 1 + w 33 r 3 was nonnegative due to w 31 > 0 and w 33 > 0. Sigmoid function (2) for the third equation of system (1) would mean that the lowest r 3 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 r 3 ( 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).

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 r i (t) , i = 1, 2, 3 denote the firing rate of neural population in the model, then a limit cycle satisfies the periodicity: for some T > 0 and all t∈ R. 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: d E ( 1 , 2 ) and d S r 3 | 1 , r 3 | 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 + . d E ( 1 , 2 ) denoted the average Euclidean distance between 1 and 2 in the three-dimensional (3D) phase space of firing rates (r 1 , r 2 , r 3 ), and the spectrum distance  d S r 3 | 1 , r 3 | 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 r 3 (t) associated with the two stimuli. The Euclidean distance was defined as follows: supposing that 1 and 2 are two ω-limit sets composed of N 1 and N 2 discrete points in three-dimensional phase space (r 1 , r 2 , r 3 ) , respectively, denoted as α 1 , α 2 , . . . , α N 1 and β 1 , β 2 , . . . , β N 2 where α i and β j are three-dimensional vectors, then 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 d E 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 d E . Therefore, we defined the spectrum distance d S to address the question of distances between representations that were sensitive to those representations being oscillations. To calculate the spectrum distance d S between two ω-limit sets 1 and 2 , only the sequence of their r 3 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 r 3 (t), was used to obtain peaks around frequency values. For an equilibrium corresponding to constant firing rate r 3 = 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 d S 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 D 1 and D 2 were the amplitudes of the peaks at zero frequency for two ω-limit sets 1 and 2 , and a i = f i , A i , i = 1, 2 denoted the corresponding alternating components of i , where f i was the non-zero frequency and A i was the amplitude of the peak around f i , then we have 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 w min 23 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 w 23 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 w 23 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.
Note that we set a i = (0, 0) if the ω-limit sets i was an equilibrium. Thus, when the two limit sets were both equilibria, d S only contained the first term measuring the difference between the direct components. In this case, d S was only a linear projection (up to a constant factor) of the Euclidean distance d E .
When both 1 and 2 were equilibria, the spectrum distance d S was a projection of the Euclidean distance d E onto the r 3 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 w 23 (Figure 4). Additionally, we found that the non-monotonic dependence of distance in both d E and d S on the feedback weight (w max 23 ) persisted, and that an optimal value for any given pair of stimuli could be found ( Figure 3D).

Bifurcation Analysis
We denoted the dynamical system of Equation (1) where r ∈ R 3 was the vector of firing rates and ∈ R p was the vector of parameters. The vector field f = f 1 , f 2 , f 3 T was a smooth function on some open set of R 3 ×R p . The dimensionality of could be up to eight dimensions maximally to include all connection weights w ij 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 w 23 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, ) = (r 0 , 0 ), i.e., the stability of this equilibrium could be determined from the linearized vector field of Equation (6) given bẏ where D r f ∂f ∂r was the Jacobian matrix of the vector field f.
If none of the eigenvalues of D r f (r 0 , 0 ) lied on the imaginary axis (i.e., the equilibrium was hyperbolic), the local stability of (r 0 , 0 ) in the non-linear system (6) could be determined by the linear system (8). The equilibrium was stable if all eigenvalues of D r f (r 0 , 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 D r f (r 0 , 0 ), there existed a unique smooth function h : R p →R 3 such that for sufficiently close to 0 . By continuity of the eigenvalues with respect to parameters, D r f (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 D r f (r 0 , 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 (r 0 , 0 ) if the Jacobian matrix D r f (r 0 , 0 ) has an eigenvalue of zero real part. In our model, a saddle-node bifurcation (SN) occurred when D r f (r 0 , 0 ) had a single zero eigenvalue (in addition to some non-degenerate conditions), and a Hopf bifurcation (H) occurred when D r f (r 0 , 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 R n , a vector field f 0 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 f 0 (i.e., vector fields that are orbitally topologically equivalent to f 0 ) form a submanifold L of codimension k, which is an equivalence class of the singular vector field f 0 . Therefore, within the space H it requires another submanifoldL of at least k dimensions to intersect transversely with L at point f 0 , such that the singularity of f 0 persists under small perturbations of the vector field. The submanifoldL was obtained through a parametrized family of vector fields involving at least k parameters. The least number k is then defined as A1 A2 B 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 w 23 . Middle: different regions on the equilibrium manifold corresponded to different regimes in the parameter space of µ and w 23 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 w 23 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 w max 23 at the same Ŵµ (solid curve). Right: a series of translated Ŵ for three representative values of Ŵµ.
the codimension of f 0 . The parametrized vector field f (r, ) in Equation (6) can be thought of as one realization of the submanifoldL which passes through the vector field f 0 f (r 0 , 0 ) undergoing a bifurcation with the two parameters corresponding to the top-down weight w 23 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) could be rewritten in the form of an autonomous system by defining the function such that using Equation (11), Equation (10) becamė where S 1 denoted a circle. To construct the Poincaré map, we defined a cross-section of Equation (12) 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 r i , 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 nonlinear 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 (r 1 , r 2 , r 3 ) and the tangent vector defining the velocity of each point along a trajectory was given by the vector field f = f 1 , f 2 , f 3 T (see 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 w 23 on firing rate dynamics. First, we varied the topdown weight w 23 (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 r i (t) , i = 1, 2, 3 were sensitive to different values of w 23 (Figure 2A, bottom). When the top down weight was small (w 23 = 4), firing rates approached the equilibrium exponentially (Figure 2A, bottom and Figure 2B, left, black traces). Conversely, when the top-down weight was large (w 23 = 10.5), the firing rate of excitatory cells (r 1 ) 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 (w 23 = 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 topdown 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 (w 23 ), 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 w 23 = 5.5, right w 23 = 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 w 23 , 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 w 23 , they reached a maximum in between ( Figure 2D1). The control of both the frequency and amplitude via changes in w 23 occurred across an array of weights (w 31 ) 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 topdown weights (w 23 ) 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 (w 23 ) 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 (w 23 = 3), the representations of the two stimuli were close (Figures 3B,C,  left). As w 23 was increased, the representations of the two stimuli were pushed apart making them more separable (w 23 = 3.8, Figures 3B,C). Interestingly, as w 23 was increased further (w 23 = 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 d E (Figure 3C) and the spectrum distance d S (Figure 3B) were different, they both occurred at the same feedback weight (Figure 3D). We visualized the distance landscapes defined by d E and d S 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 w max 23 that maximized the distance between the two resultant representations for any given pair of stimuli.
A landscape of the optimal w max 23 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 (r 1 and r 2 ). For a given stimulus, an increase in w 23 led to a monotonic decrease in r 1 , suggesting persistent suppression onto the local M population. By contrast, the response of the local G population r 2 was elevated first with increasing w 23 until reaching the maximum r max 2 , after which r 2 dropped significantly ( Figure 4D). Across different stimuli µ, the shape of the firing rate r 2 as a function of w 23 remained the same but shifted vertically. To determine if these differences in the firing rate of inhibitory neurons (r 2 ) were related to the values of top-down weights that maximally separated the distance between two stimulus representations, we plotted the w max 23 (abscissa) obtained from the landscape in Figure 4C vs. the midpoint w mid 23 (ordinate) between the r max 2 of the same stimulus pair (inset of Figure 4E). The response r 2 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 w max 23 was correlated to w mid 23 (R 2 = 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 w 23 . 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 (w 23 = 12.0), one stimulus (µ 1 = 1.15) generated oscillations (f 1 = 33.5 Hz, Figure 5A1, before) in the piriform population's firing rate that were different in both frequency and amplitude from the oscillations (f 2 = 37.Hz, Figure 5A1, before) in response to a second stimulus (µ 2 = 0.6). However, a change in the topdown weight (w 23 = 8), resulted in firing rate oscillations becoming more similar for the same two stimuli (Figure 5A1, after, f 1 = 30.4 Hz, f 2 = 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 d S (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 w 23 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 w min 23 across combinations of stimuli (µ 1 vs. µ 2 ) that was most effective at generating synchronous oscillations ( Figure 5C).
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 w 23 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 w 23 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 w 23 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 w 23 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 (ω s t + ϕ 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).
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 w 23 , 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 w 23 , 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 topdown 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 w 23 . At some critical values of w 23 , 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 w 23 was within a specific regime. Therefore, across all combinations of external stimuli µ and top-down weight w 23 , 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 tw 23 .
Within the manifold of the stimulus µ and top-down weight w 23 , 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 (µ, w 23 ) 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 w 23 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 w 23 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 w 23 max 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 w max 23 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 w max 23 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 w max 23 followed closely with the transition boundary Ŵ (Figure 8B, left). As the stimulus difference increased (Figure 8B, middle), w max 23 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 w 23 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.

DISCUSSION
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, longterm 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 topdown 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.

FUNDING
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 https://www.biorxiv.org/content/10.1101/2020.02. 25.964965v1 (Chen and Padmanabhan, 2020).