Network mechanism for insect olfaction

Abstract

Early olfactory pathway responses to the presentation of an odor exhibit remarkably similar dynamical behavior across phyla from insects to mammals, and frequently involve transitions among quiescence, collective network oscillations, and asynchronous firing. We hypothesize that the time scales of fast excitation and fast and slow inhibition present in these networks may be the essential element underlying this similar behavior, and design an idealized, conductance-based integrate-and-fire model to verify this hypothesis via numerical simulations. To better understand the mathematical structure underlying the common dynamical behavior across species, we derive a firing-rate model and use it to extract a slow passage through a saddle-node-on-an-invariant-circle bifurcation structure. We expect this bifurcation structure to provide new insights into the understanding of the dynamical behavior of neuronal assemblies and that a similar structure can be found in other sensory systems.

Appendices

Appendix

Appendix 1: Firing-rate model derivation

In this appendix, we present the details of obtaining the FR model, Eqs. (4) and (5) and the inequality (7) presented in “Firing-Rate Model” section, starting from the I&F model in Eqs. (1) and (2) driven by a “white-odor” stimulus. We, in fact, derive the FR model for a slightly more general connectivity architecture than the all-to-all coupled network, discussed in “Firing-Rate Model” section, namely, unstructured, uniform random connectivity as assumed for the most general version of the I&F model in “Integrate-and-Fire Model” section.

We heuristically derive the FR model by employing the following assumptions:

• (i) We assume that the excitatory and fast inhibitory conductances rise instantaneously compared to the slow inhibitory conductance.

• (ii) We consider the mean-driven regime, in which both network-generated and external-drive spikes are small and arrive at high rates. Sums over incoming spike trains in Eqs. (2) are thus replaced by continuous functions representing the incoming firing rates, since they produce the same conductance dynamics.

• (iii) Individual synaptic connections are replaced by the corresponding connection probabilities, $p^E_i$, $p^F_i$, and $p^S_i$, where each probability takes one of the corresponding two values, $p^P_Q$, $P=E$, F, and S, $Q=E$ or I, depending on the $i{{\text {th}}}$ neuron’s type, as described in the second paragraph below Eqs. (2). In this way, firing rates are scaled by the average number of synaptic connections. Together with the coupling strengths, $S^E_i$, $S^F_i$, and $S^S_i$, where each coupling strength likewise takes one of the corresponding two values, $S^P_Q$, $P=E$, F, and S, $Q=E$ or I, the network-drive terms in Eqs. (2) become $S^E_ip^E_im_E(t)$, $S^F_ip^F_im_I(t)$, and $S^S_ip^S_im_I(t)$, respectively.

• (iv) We assume that the per-neuron firing rates, $m_E(t)$ and $m_I(t)$, vary slowly in comparison to changes in the membrane potential and further assume that the total conductance is high, allowing us to treat the voltage in Eq. (1) as a constant coefficient differential equation. We can then directly solve for $m_E(t)$ and $m_I(t)$ as the multiplicative inverses of the time for the corresponding voltage solutions to reach threshold from reset.

• (v) We neglect the refractory period.

Using Assumption (i), we can adiabatically eliminate the auxiliary variables $h_i^E(t)$ and $h_i^F(t)$ in Eqs. (2b) and (2c), allowing the incoming spikes to cause instantaneous jumps in the fast conductances, $g_i^E(t)$ and $g_i^F(t)$, respectively. In other words, we describe the dynamics of the fast conductance variables themselves using only first-order kinetics in contrast to Eqs. (2) in the I&F model which use second-order kinetics. We thus replace Eqs. (2) with the equations

$$\begin{aligned} \sigma _E \frac{dg_i^E(t)}{dt}&= -g_i^E(t) +f_{i}^{{\text {odor}}} \sum _{l} \delta (t - \tau ^i_l)+ f_{i} \sum _{k} \delta (t - \gamma ^i_k) \nonumber \\&\quad + \frac{S^E_{i}}{N_E}\sum _{j \ne i}p^E_{ji}\sum _{\mu }\delta (t-t^j_{\mu }), \end{aligned}$$

(9a)

$$\begin{aligned} \sigma _F \frac{dg_i^F(t)}{dt}&= -g_i^F(t) + \frac{S^F_{i}}{N_I}\sum _{j \ne i}p^F_{ji}\sum _{\mu }\delta (t-t^j_{\mu }), \end{aligned}$$

(9b)

$$\begin{aligned} \sigma _S \frac{dg_i^S(t)}{dt}&= -g_i^S(t) +h^S_i(t), \end{aligned}$$

(9c)

$$\begin{aligned} \rho _S \frac{dh_i^S(t)}{dt}&= -h_i^S(t) + \frac{S^S_{i}}{N_I}\sum _{j \ne i}p^S_{ji}\sum _{\mu }\delta (t-t^j_{\mu }), \end{aligned}$$

(9d)

where all variables and parameters are defined for Eqs. (2) in “Integrate-and-Fire Model” section. Consequently, we must readjust the values of the excitatory and fast inhibitory conductance decay rates, $\sigma _E$ and $\sigma _F$, respectively, so that the conductance responses have the appropriate duration lengths.

The next step is to consider the mean-driven limit of the incoming spikes in the model driven by a “white odor” stimulus. To mimic the “white odor” stimulus, the external spikes from the background and odor are combined to form a single source of external input with Poisson rate $\nu$, and spike strengths $f_E$ and $f_I$ for the excitatory and inhibitory neurons, respectively. We then further assume that the external drive operates in the mean-driven regime, in which each individual spike is small but the spikes arrive at high rates (i.e., $f_i\rightarrow 0,\,\nu \rightarrow \infty$ with $f_i\nu$ held constant), replacing the sums over the incoming spike trains in Eq. (9a), $f_{i}^{{\text {odor}}} \sum _{l} \delta (t - \tau ^i_l)+ f_{i} \sum _{k} \delta (t - \gamma ^i_k)$, with their statistical averages, the constant functions $f_i\nu$, where $i = E$ or I depending on whether the $i^{\text {th}}$ neuron belongs to the excitatory or inhibitory population, respectively.

Furthermore, as we only consider networks with unstructured, uniform random connectivity, we can assume that probability of a synaptic connection between any pair of neurons of given types will be statistically equivalent to that of any other synaptic connections between any other pair of the same types of neurons. Thus, statistically, we can replace the individual synaptic connection weights, $p^E_{ij}$, $p^F_{ij}$, and $p^S_{ij}$, by the corresponding connection probabilities between the corresponding populations, $p_Q^E$, $p_Q^F$, and $p_Q^S$, respectively.

From the coarse-graining process carried out in the previous two paragraphs, we can conclude that the resulting conductance equations describe the conductance dynamics of a typical (excitatory or inhibitory) neuron, i.e., one that is statistically equivalent to all other neurons of the same type. Thus, we replace the individual neuronal conductance equations in Eqs. (9) by equations for the conductances of the corresponding neuronal populations. This replacement eliminates the need for the index i and lets us introduce the subscript notation of Q in its place, where Q is either E or I to represent quantities associated with the typical excitatory or inhibitory neuron, or, equivalently, the corresponding populations, respectively.

We now apply the same mean-driven limit to the network-generated spikes as we did to the external spikes. Due to the uniform, unstructured random connectivity of the network, we can assume that each excitatory and inhibitory neuron, respectively, is driven by the average per-neuron firing rate of the appropriate type, $m_E(t)$ and $m_I(t)$. We scale these by the number of neurons of each type in the network and their synaptic connection probabilities to find the population-averaged firing rates, $N_Ep^E_Qm_E(t)$, $N_Ip^F_Qm_I(t)$, and $N_Ip^S_Qm_I(t)$, arriving at the acetycholine, $\hbox {GABA}_A$, and slow inhibitory receptors, respectively, on the postsynaptic excitatory ($Q=E$) and inhibitory ($Q=I$) neurons. Thus, the sums over the network spikes, $\frac{S^E_{i}}{N_E}\sum _{j \ne i}p^E_{ji}\sum _{\mu }\delta (t-t^j_{\mu })$, $\frac{S^F_{i}}{N_I}\sum _{j \ne i}p^F_{ji}\sum _{\mu }\delta (t-t^j_{\mu })$, and $\frac{S^S_{i}}{N_I}\sum _{j \ne i}p^S_{ji}\sum _{\mu }\delta (t-t^j_{\mu })$, in Eqs. (9a), (9b), and (9d), are replaced by the average network drive terms $S^E_Qp^E_Qm_E(t)$, $S^F_Qp^F_Qm_I(t)$, and $S^S_Qp^S_Qm_I(t)$, respectively.

In this way, since each excitatory or inhibitory neuron is statistically equivalent to all other neurons of the same type, the set of 4N equations for the neurons’ conductances in Eqs. (9) can be replaced by the following eight representative equations for the two neuronal populations:

$$\begin{aligned} \sigma _E\frac{dg_Q^E}{dt}&= -g_Q^E+ f_Q\nu +S^E_Qp^E_Qm_E(t), \end{aligned}$$

(10a)

$$\begin{aligned} \sigma _F\frac{dg_Q^F}{dt}&= -g_Q^F+S^F_Qp^F_Qm_I(t), \end{aligned}$$

(10b)

$$\begin{aligned} \sigma _S\frac{dg^S_Q}{dt}&= -g^S_Q + h_Q, \end{aligned}$$

(10c)

$$\begin{aligned} \rho _S\frac{dh_Q}{dt}&= -h_Q +S^S_Qp^S_Qm_I(t), \end{aligned}$$

(10d)

with $Q=E$ or I for the excitatory or inhibitory population, respectively.

To solve for the firing rates, we treat the voltage equation in Eq. (1) as a constant-coefficient differential equation; we have one such equation for each of the two populations. This is justified by assuming that the per-neuron firing rates, $m_E(t)$ and $m_I(t)$, and therefore also the associated conductances, vary slowly compared to the voltages, while also assuming that the conductances are relatively high so that the voltages vary fast. This means that the (conductance-induced) time scale of the voltage is shorter than the shortest conductance time scale (see ref. Shelley et al. 2002). For constant conductance values $g^E_Q,\,g^F_Q,$ and $g^S_Q$, the first-order linear equation in Eq. (1) becomes a constant coefficient equation, and takes the form

$$\begin{aligned} \begin{aligned} \tau \frac{dv_Q(t)}{dt} =&-\left( 1+g^E_Q+ g^F_Q + g^S_Q\right) v_Q(t) + \left( \varepsilon _R + g^E_Q \varepsilon _E + g^F_Q\varepsilon _F + g^S_Q\varepsilon _S\right) , \end{aligned} \end{aligned}$$

(11)

which has the solution

$$\begin{aligned} \begin{aligned} v_Q(t)&= Ce^{-\left( 1+g^E_Q+ g^F_Q + g^S_Q\right) t/\tau } + \frac{\varepsilon _R + g^E_Q \varepsilon _E + g^F_Q\varepsilon _F + g^S_Q\varepsilon _S}{1+g^E_Q+ g^F_Q + g^S_Q}, \end{aligned} \end{aligned}$$

(12)

with arbitrary constant of integration, C. This constant is determined by the initial condition that the membrane potential starts at its resting potential value, $\varepsilon _R$, just after spiking,

$$\begin{aligned} v_Q(0) = \varepsilon _R, \end{aligned}$$

(13)

and thus,

$$\begin{aligned} C = \frac{g^E_Q(\varepsilon _R - \varepsilon _E)+g^F_Q(\varepsilon _R -\varepsilon _F)+g^S_Q(\varepsilon _R -\varepsilon _S)}{1+g^E_Q+g^F_Q+ g^S_Q}. \end{aligned}$$

(14)

At the time $t^{\text {spike}}_Q$, the membrane potential in Eq. (12) with C in Eq. (14), reaches the threshold,

$$\begin{aligned} v_Q(t^{\text {spike}}_Q) = V_T. \end{aligned}$$

(15)

Therefore, $t^{\text {spike}}_Q$ is the amount of time between a neuron’s consecutive spikes given the particular values of the conductances. The reciprocal of $t^{\text {spike}}_Q$ is thus the firing rate of the typical neuron in the population, $m_Q$, at those conductance values.

The membrane potential, $v_Q$, only crosses the firing threshold when the corresponding slaving potential,

$$\begin{aligned} V_{s,Q} = \frac{ \varepsilon _R+ g^E_Q \varepsilon _E+g^F_Q\varepsilon _F+g^S_Q\varepsilon _S}{1+g^E_Q+g^F_Q+g^S_Q}, \end{aligned}$$

(16)

is larger than the spiking threshold, $V_T$, as the membrane potential is always being drawn toward $V_{s,Q}$. The condition $V_{s,Q} > V_T$ can be rewritten as a condition on the effective excitatory conductance $g^E_Q$ such that

$$\begin{aligned} g^E_Q > \frac{V_T - \varepsilon _R}{\varepsilon _E - V_T} + g^F_Q\left( \frac{V_T - \varepsilon _F}{\varepsilon _Q - V_T}\right) +g^S_Q\left( \frac{V_T - \varepsilon _S}{\varepsilon _Q - V_T}\right) . \end{aligned}$$

(17)

Using Eqs. (12) through (15), we can calculate the spike time $t^{\text {spike}}_Q$ for the effective conductance region where these conditions hold, arriving at the expression

$$\begin{aligned} t^{\text {spike}}_Q=&\frac{\tau }{1+g^E_Q+g^F_Q+g^S_Q} \ln \left[ \tfrac{g^E_Q(\varepsilon _E - \varepsilon _R)+g^F_Q(\varepsilon _F - \varepsilon _R)+g^S_Q(\varepsilon _S- \varepsilon _R)}{\left\{ \varDelta \varepsilon _R+ g^E_Q\varDelta \varepsilon _E+g^F_Q\varDelta \varepsilon _F+g^S_Q\varDelta \varepsilon _S\right\} ^+}\right] , \end{aligned}$$

(18)

where $\varDelta \varepsilon _Z$ is defined in Eq. (6) and $\{ \cdot \}^+$ is defined as after Eqs. (5).

As stated above, the reciprocal of $t^{\text {spike}}_Q$ in Eq. (18) is the per-neuron firing rate $m_Q$, given by the equation

$$\begin{aligned} m_Q(t) = \frac{1+g_Q^E+ g^F_Q+ g^S_Q}{\tau \ln \left[ \frac{g_Q^E(\varepsilon _E - \varepsilon _R)+ g^F_Q(\varepsilon _F-\varepsilon _R)+g^S_Q(\varepsilon _S-\varepsilon _R)}{ \left\{ \varDelta \varepsilon _R + g_Q^E \varDelta \varepsilon _E + g^F_Q\varDelta \varepsilon _F+g^S_Q\varDelta \varepsilon _S\right\} ^+ }\right] }, \end{aligned}$$

(19)

where, again, $Q=E$ or I. Note that, in the case when $V_{s,Q} < V_T$, no spiking occurs and the average time to spike, $t^{\text {spike}}_Q$, would be taken as infinite, so that the per-neuron firing rate vanishes. The firing rate, $m_Q$ in Eq. (19), is inserted into the conductance equations in Eqs. (10) to close the FR model.

As stated in the previous paragraph, when the denominator inside the logarithm in Eq. (19) vanishes, we define the corresponding firing rate, $m_Q$, $Q \in \left\{ E,I\right\}$, to also vanish. Thus, we conclude that the firing rate $m_Q(t)$ is a piecewise-defined, continuous function with $m_Q(t)=0$ when

$$\begin{aligned} \varDelta \varepsilon _R + g_Q^E \varDelta \varepsilon _E + g^F_Q\varDelta \varepsilon _F+g^S_Q\varDelta \varepsilon _S \le 0, \end{aligned}$$

(20)

and is given by Eq. (19) otherwise.

Using the relations

$$\begin{aligned} g^E_Q&= S^E_Qp^E_Qg^E+f_Q\nu , \end{aligned}$$

(21a)

$$\begin{aligned} g^F_Q&= S^F_Qp^F_Qg^F, \end{aligned}$$

(21b)

$$\begin{aligned} g^S_Q&= S^S_Qp^S_Qg^S, \end{aligned}$$

(21c)

$$\begin{aligned} h_Q&= S^S_Qp^S_Qh, \end{aligned}$$

(21d)

pairs of equations in Eqs. (10) for each type of conductance become redundant and the FR system, defined by Eqs. (10) and (19) and the inequalities in Eqs. (20) can be reduced to the four-equation model for the effective conductance variables $g^E, g^F, g^S$, and h, given by the conductance equations in Eqs. (4) presented in “Firing-Rate Model” section, together with the firing-rate equations

$$\begin{aligned} m_Q(t) =\frac{1+S^E_Qp^E_Qg^E+f_Q\nu +S^F_Qp^F_Qg^F+S^S_Qp^S_Qg^S}{\tau \ln \left[ M_Q\right] } , \end{aligned}$$

(22a)

where

$$\begin{aligned} M_Q = \frac{ \sum _{P\in \{E,F,S\}} S^P_Qp^P_Qg^P(\varepsilon _P - \varepsilon _R)+f_Q\nu (\varepsilon _E - \varepsilon _R)}{\left\{ \varDelta \varepsilon _R +f_Q\nu \varDelta \varepsilon _E+\sum _{P\in \{E,F,S\}}S^P_Qp^P_Qg^P\varDelta \varepsilon _P \right\} ^+ }. \end{aligned}$$

(22b)

The conditions for which the firing rates in Eqs. (22) do not vanish become

$$\begin{aligned} \varDelta \varepsilon _R +f_Q\nu \varDelta \varepsilon _E+\sum _{P\in \{E,F,S\}}S^P_Qp^P_Qg^P\varDelta \varepsilon _P > 0. \end{aligned}$$

(23)

For an all-to-all coupled network, we take the connection probabilities $p_Q^E=p_Q^F=p_Q^S=1$, and the above reduced model becomes that given by the equations in Eqs. (4) and (5) and the inequalities in Eqs. (7) in “Firing-Rate Model” section.

We remark on the above assumption (iv), which states that the dynamically-induced membrane-potential time scales are taken to be faster than the conductance time scales. This is rarely the case. In particular, the voltage time scale almost never falls below 50% of the conductance time scale; nevertheless, the frozen-conductance approximation in the voltage equation is quite accurate (Shelley et al. 2002). Therefore, despite this accurate agreement, the derivation of the FR model presented here is not systematic in terms of any small parameter, and thus the FR model cannot be considered as a limit of the I&F model. While we can expect qualitative agreement between the dynamics of the I&F model and the FR model, we can only expect to achieve a more quantitative agreement after possibly tuning some of the FR model parameters.

Appendix 2: Robustness of the binding-index threshold

In this appendix, we discuss the robustness of our choice for the value of the binding index threshold b, defined in “Binding index (BI)” section. As mentioned in “Odor Response and Discriminability” section, the threshold value $b = 0.65$ we use for the binding index is the same as that used in ref. Patel et al. (2013). There, it is shown that this choice of b is suitable to identify a distinct, clearly defined subset of the PNs that are temporally-bound given each odor, and only indirectly driven by the corresponding excited ORNs, which can be used to distinguish among different odors. This is also true in the results of our simulations, as we have discussed in “Odor Response and Discriminability” section and displayed in Fig. 3 and Table 2.

In particular, in Fig. 9 we use the results of the same numerical simulations as those used in Fig. 3, but we now analyze them using three different values of the binding-index threshold: $b=0.2$, $b = 0.65$, and $b=0.9$. (Note that the middle column of Fig. 9 is identical to Fig. 3.) We see that only a few neurons that are bound at $b = 0.65$ change their bound/unbound property between the values of $b=0.2$ and $b=0.9$, and most of those neurons are driven directly by the odor stimulus arriving from the ORNs. For Odor 1 in the top row of Fig. 9 (spiking activity shown in blue), we found that the only indirectly-driven neuron that changes from a bound to unbound classification with the increase in threshold from 0.65 to 0.9 is neuron 69. For Odor 3 in the bottom row of Fig. 9 (spiking activity shown in green), we found that only the indirectly-driven neurons 37 and 51 change their classification. Finally, for Odor 2 in the middle row of Fig. 9 (spiking activity shown in red), no indirectly-driven neurons change their bound/unbound classification between threshold values of 0.65 and 0.9. Additionally, when comparing spiking behavior at a threshold of $b = 0.2$ to that at $b = 0.65$ we observe that all indirectly-driven bound neurons remain bound at both thresholds for all three model odors. Therefore, Table 2 also provides the data for $b = 0.2$. Overall, the largest change for the odors shown occurs for Odor 1 from $b = 0.65$ to $b = 0.9$ where only 8% of the PNs (directly- and indirectly-driven combined) change classification from bound to unbound. Thus, the raster plots of bound neurons between the threshold values $b=0.2$ and $b=0.9$ can still clearly discriminate among different odors.

From the discussion in the previous paragraph, we conclude that when we consider the portion of indirectly-driven, temporally-bound neurons that are uniquely assocaited to a single odor in a pair-wise comparison between two odors for threshold values of $b = 0.65$ and $b = 0.9$, we find that there are overall fewer indirectly-driven, bound neurons at the higher threshold but that a greater than or equal number and a larger percentage of them are uniquely associated to that odor. This holds across all three pair-wise comparisons and is supported by the data provided in Tables 2 and 4.

We identified this interval from $b = 0.2$ to $b = 0.9$ somewhat arbitrarily to be near, but not exactly equal to the threshold values of $b = 0$ and $b = 1$, respectively. Choosing these extreme values would be prohibitively inclusive at $b = 0$ (neurons in all triplets would be classified as temporally-bound) and exclusive at $b = 1$ (only neurons in triplets with identical spiking behavior would be classified as bound).

The discussion in this appendix thus shows that there exists a large interval of threshold values of the binding index, $b$, for which the odor-discriminability data quantitatively changes little but the ability to discriminate among odors qualitatively changes none at all.

Fig. 9

Raster plots of the excitatory neurons in a network to which three different odors are presented, generated by stimulating a different 1/3 of the neuron in the network. Each column corresponds to a different binding index b, for which temporally-bound neurons are those with binding index $BI_{i,j,k} \ge b$. In all panels, black dots represent unbound neurons, stars represent temporally-bound neurons directly driven by the odor, and open symbols represent temporally-bound neurons not directly driven by the odor stimulus. The figure illustrates that many bound neurons are stable with respect to binding index threshold b. Parameters found in Table 1

Appendix 3: Limit cycle in the linearized firing-rate model

In this appendix, we study a linearization of the FR model in Eqs. (4) and (22), with an eye on providing a demonstration that a unique, attracting limit cycle exists in the two-dimensional FR model for the fast effective conductances, $g^E$ and $g^F$, with frozen slow effective conductance, $g^S$. Thus, in  “Appendix 3.1”, we derive this linearization, and in  “Appendix 3.2” we semi-numerically verify the existence of a unique limit cycle in the resulting planar piecewise-linear system for $g^E$ and $g^F$ with frozen $g^S$. In particular, we show that we can find the trajectories explicitly in the four effective-conductance regions in which this system is linear, and splice them together using a numerical solution of a transcendental equation at each boundary between two such regions. We use these spliced trajectories to construct a Poincaré map, and demonstrate the existence of the unique limit cycle by finding a unique fixed point of this map.

Table 4 Quantities related to temporally-bound, indirectly-driven neurons during initial network oscillations characterized by a Binding Index threshold of $b = 0.9$

Appendix 3.1: Derivation of the linearized model

In this section, we present the details of linearizing the FR model in Eqs. (4) and (22), which results in a piecewise-linear model.

Our procedure is to first note that if a firing rate $m_Q$, $Q =E$ or I, vanishes, the corresponding equation(s) in Eqs. (4) is (are) already linear. If however, a firing rate is nonzero, we linearize it as a function of the conductances in order to linearize the corresponding equation(s) (4a), (4b), or (4d). This procedure creates a piecewise-linear function for the derivative of each effective conductance variable: $g^E$, $g^F$, and h; the firing-rate term does not appear in the differential equation for the slow effective conductance $g^S$, and thus Eq. (4c) does not take a piecewise form. A crucial step involves determining the boundaries between the four regions created by the various combinations of parts of the piecewise-defined firing rates: $m_E = m_I = 0$, $m_E=0$ and $m_I\ne 0$, $m_E\ne 0$ and $m_I=0$, $m_E, m_I \ne 0$, along which the system switches among the differential equations that govern the trajectory.

When a slaving voltage $V_{s,Q}$ in Eq. (16) is subthreshold, $V_{s,Q} <V_T$, the inequality in Eq. (17) is not satisfied, and the corresponding firing rate vanishes, $m_Q(t) = 0$, $Q=E$ or I. It would follow then that each equation in Eqs. (4) would take a linear form. In the superthreshold case when $V_{s,Q} > V_T$, $m_Q(t)$ is nonlinear for each $Q =E$ or I. We linearize about large values of the effective excitatory conductance $g^E$. (See more discussion at the end of this section). In order for our linearization to proceed more systematically, we first expand the firing rates $m_Q(t)$ about large values of the conductances $g^E_Q$, where $Q =E$ or I, and only then express $g^E_Q$ in terms of the effective conductance $g^E$ using formulas in Eqs. (21).

By linearizing the firing rates in the superthreshold case, we can produce a system of piecewise-linear equations for the derivatives of the effective conductances. We then use the reduced two-dimensional fast model, in which $g^S$ is treated as a frozen parameter, to demonstrate the existence of a limit cycle analytically, save for four root-finding calculations. The details of this analysis are given in “Appendix 3.2”, below.

To more easily identify the large parameter $g^E_Q$, it is helpful to rearrange the fraction inside of the logarithm in Eq. (19) as follows:

$$\begin{aligned} \begin{aligned}&\frac{g^E_Q(\varepsilon _E - \varepsilon _R)+ g^F_Q(\varepsilon _F -\varepsilon _R)+g^S_Q(\varepsilon _S -\varepsilon _R)}{\varDelta \varepsilon _R+ g^E_Q \varDelta \varepsilon _E+g^F_Q\varDelta \varepsilon _F+g^S_Q\varDelta \varepsilon _S } \\&\qquad =\frac{A\left[ 1-\left( Bg^F_Q + Xg^S_Q\right) /g^E_Q\right] }{ 1-\left[ (A-1)- \frac{\varDelta \varepsilon _F}{\varDelta \varepsilon _E}g^F_Q- \frac{\varDelta \varepsilon _S}{\varDelta \varepsilon _E}g^S_Q\right] /g^E_Q}, \end{aligned} \end{aligned}$$

(24)

with A, B, and X defined as

$$\begin{aligned}&A = \frac{\varepsilon _E - \varepsilon _R}{\varepsilon _E - V_T}, \quad B = \frac{\varepsilon _R -\varepsilon _F}{\varepsilon _E - \varepsilon _R}, \quad X = \frac{\varepsilon _R -\varepsilon _S}{\varepsilon _E - \varepsilon _R}. \end{aligned}$$

The logarithm of Eq. (24) can be rearranged so that the resulting expression may be expanded for large $g^E_Q$,

$$\begin{aligned}&\log \left\{ \frac{A\left[ 1-\left( Bg^F_Q+ Xg^S_Q\right) /g^E_Q\right] }{ 1-\left[ (A-1)- \frac{\varDelta \varepsilon _F}{\varDelta \varepsilon _E}g^F_Q- \frac{\varDelta \varepsilon _S}{\varDelta \varepsilon _E}g^S_Q\right] /g^E_Q}\right\} = \log A + \log \left[ 1-\left( Bg^F_Q + Xg^S_Q\right) /g^E_Q\right] \nonumber \\&\quad - \log \left\{ 1-\left[ (A-1)- \frac{\varDelta \varepsilon _F}{\varDelta \varepsilon _E}g^F_Q- \frac{\varDelta \varepsilon _S}{\varDelta \varepsilon _E}g^S_Q\right] /g^E_Q\right\} \nonumber \\&\quad = \left[ (A-1) \left( 1 + (B+1)g^F_Q+(X+1)g^S_Q\right) \right] /g^E_Q + \log A + \mathcal {O}\left( {g^{E}_Q}^{-2}\right) . \end{aligned}$$

(25)

When the expression in Eq. (25) replaces the logarithm in the denominator of Eq. (19), the resulting version of the firing rate, $m_Q(t)$, can be expanded once again for large $g^E_Q$ to produce the expression

$$\begin{aligned}&\frac{1+g^E_Q+g^F_Q+g^S_Q}{\tau \log A} \bigg [1 - \frac{(A-1)}{\log A}\left[ 1 + (B+1)g^F_Q +(X+1)g^S_Q\right] \frac{1}{g^E_Q} + \mathcal {O}\left( g^{E-2}_Q\right) \bigg ]. \end{aligned}$$

(26)

To complete the linearization, we neglect higher order terms in Eq. (26) and thus arrive at a linear firing-rate equation for the superthreshold case ($V_{s,Q}>V_T$). Considering again that the firing rate is taken as vanishing in the subthreshold case, when linearized about large $g^E_Q$, Eq. (19) turns into the equation

$$\begin{aligned} \tilde{m}_Q&= \left\{ \frac{1}{\tau \log A}\left[ \left( 1- \frac{A-1}{\log A}\right) + g^E_Q \right. \right. \end{aligned}$$

$$\begin{aligned}&\quad \left. \left. +\left( 1- \frac{(A-1)(B+1)}{\log A}\right) g^F_Q +\left( 1- \frac{(A-1)(X+1)}{\log A}\right) g^S_Q\right] \right\} ^+, \end{aligned}$$

(27)

where, again, $\{x \}^+ = x$ if $x>0$, and zero otherwise. The inequality condition for $\tilde{m}_Q$ not to vanish, similarly to the condition for the nonlinear firing rate in Eq. (23), is given by the inequality

$$\begin{aligned}&g^E_Q +\left( 1- \frac{(A-1)(B+1)}{\log A}\right) g^F_Q +\left( 1- \frac{(A-1)(X+1)}{\log A}\right) g^S_Q+ \left( 1- \frac{A-1}{\log A}\right) > 0. \end{aligned}$$

(28)

Using the relations in Eqs. (21), the piecewise-linear, per-neuron firing rate for each excitatory neuron (taking $Q = E$) becomes a function of the effective conductance variables, $g^E$, $g^F$, and $g^S$:

$$\begin{aligned} \tilde{m}_E&=\left\{ \frac{1}{\tau \log A}\left[ \left( S_E^Ep_E^Eg^E+f_E\nu \right) + 1- \frac{A-1}{\log A} +\left( 1- \frac{(A-1)(B+1)}{\log A}\right) S_{E}^Fp_{E}^Fg^F\right. \right. \nonumber \\&\quad \left. \left. +\left( 1- \frac{(A-1)(X+1)}{\log A}\right) S_{E}^Sp_{E}^Sg^S \right] \right\} ^+. \end{aligned}$$

(29)

Likewise, taking $Q = I$, the piecewise-linear, per-neuron, inhibitory firing rate is defined as

$$\begin{aligned} \tilde{m}_I&= \left\{ \frac{1}{\tau \log A}\left[ \left( 1-\frac{(A-1)(B+1)}{\log A}\right) S_{I}^Fp_{I}^Fg^F+ \left( 1-\frac{(A-1)(X+1)}{\log A}\right) S_{I}^Sp_{I}^Sg^S \right. \right. \nonumber \\&\quad + \left. \left. S_I^Ep_I^Eg^E+f_I\nu +1 - \frac{A-1}{\log A} \right] \right\} ^+. \end{aligned}$$

(30)

The linearized, four-dimensional system for the derivatives of the effective conductance variables, incorporating Eqs. (29) and (30), is

$$\begin{aligned} \sigma _E\frac{dg^E}{dt}&= -g^E + \tilde{m}_E(t), \end{aligned}$$

(31a)

$$\begin{aligned} \sigma _F\frac{dg^F}{dt}&= -g^F + \tilde{m}_I(t), \end{aligned}$$

(31b)

$$\begin{aligned} \sigma _S\frac{dg^S}{dt}&= -g^S + h, \end{aligned}$$

(31c)

$$\begin{aligned} \rho _{S}\frac{dh}{dt}&= -h + \tilde{m}_I(t). \end{aligned}$$

(31d)

As the only nonlinearities in Eqs. (4) appear in the firing-rate terms, Eqs. (31) are piecewise-linear due to the piecewise-linear nature of the linearized firing rates $\tilde{m}_E$ and $\tilde{m}_I$. The conditions upon which different parts of the piecewise-linear differential equations apply are thus the conditions required for $\tilde{m}_Q$ not to vanish, given in Eq. (28) with the relations in Eqs. (21). The three-dimensional $g^Eg^Fg^S$-space is thus divided into four regions by the two interesting planes spanned by $g^E$, $g^F$, and $g^S$, and defined by setting the conditions in Eqs. (28) to zero. As solutions of Eqs. (31) are computed, the region within which the trajectory moves dictates which parts of the piecewise-defined derivatives describe the behavior in that region.

We now comment on why we linearize the FR model by expanding about large values of the excitatory conductances. This is because, in that regime, firing rates are also large. At large firing rates, in turn, the response of a neuronal network (such as an I&F model) asymptotically becomes independent of the amount of fluctuations in the dynamics; see ref. Kovačič et al. (2009), especially the results of numerical simulations shown in Fig. 8. Additionally, the response of the linearized FR model discussed in this section turns out to resemble that of a neuronal network in the fluctuation-driven regime. Therefore, the linearized FR model has a convenient physiological interpretation: it provides a FR model corresponding to an I&F model in the fluctuation-driven regime.

When the slow effective conductance $g^S$ is held fixed, we can depict the nonlinear and linearized derivatives of the excitatory and fast inhibitory effective conductance variables $g^E$ and $g^F$ given in Eqs. (4) and (31), as surfaces over the $g^E$-$g^F$-plane. In Fig. 10, we show representative slices through these surfaces, and thus compare each of the nonlinear derivatives of the effective fast conductances $g^E$ and $g^F$, shown by solid red lines, with their respective linearized derivatives, shown by dashed black lines, over an interval of $g^E$, for a fixed value of $g^F$. Notice that the two sets of curves are tangential in the limit of large excitatory effective conductance $g^E$.

Fig. 10

The curves representing the nonlinear (solid red lines) and piecewise-linear (dashed black lines) functions that define a $dg^E/dt$, in Eqs. (4a) and (31a) as $g^E$ varies, with $g^F = 4$, and b $dg^F/dt$ in Eqs. (4b) and (31b), versus $g^E$ for $g^F = 9.6$. Note the agreement between the nonlinear and piecewise-linear curves for large $g^E$ in both cases. Parameters in Table 1, except for the parameters pertaining to the slow effective conductance $g^S$. The slow effective conductance $g^S$ is fixed at $g^S=0$

Appendix 3.2: Existence of a unique limit cycle

In this section, we use an approach alternative to that presented in “Bifurcation Mechanism” section in order to verify the existence of a unique limit cycle in the piecewise-linear FR model, Eqs. (31), with the effective slow conductance variable, $g^S$, held fixed at a constant value. We achieve this verification with a combination of analytical solutions within each of the above-described regions in which the piecewise-defined linearized FR model has a linear form, and numerical solves of single-variable transcendental equations at each boundary between two such regions. Lastly, we confirm the existence of a limit cycle through an iterative process of possible initial conditions and associated solution trajectories.

With the simplification of the effective slow conductance, $g^S$, as a fixed value, the piecewise-linear system in Eq. (31) becomes a two-dimensional system of the excitatory and fast inhibitory effective conductances $g^E$ and $g^F$, respectively. The variable h is also no longer necessary when $g^S$ is held fixed. The two-dimensional system is given by the equations

$$\begin{aligned} \sigma _E\frac{dg^E}{dt}&= -g^E + \tilde{m}_E(t), \end{aligned}$$

(32a)

$$\begin{aligned} \sigma _F\frac{dg^F}{dt}&= -g^F + \tilde{m}_I(t), \end{aligned}$$

(32b)

with $\tilde{m}_E(t)$ and $\tilde{m}_I(t)$ given in Eqs. (29) and (30), respectively, for a fixed value of $g^S$.

Additionally, the three-dimensional boundary surfaces described in “Appendix 3.1” for the intact linearized FR model, become two intersecting lines in the $g^E$-$g^F$-plane given by the equations of the firing rates with the slow effective conductance $g^S$ fixed,

$$\begin{aligned} \tilde{m}_E&= 0, \end{aligned}$$

(33a)

$$\begin{aligned} \tilde{m}_I&= 0, \end{aligned}$$

(33b)

with the firing rates $\tilde{m}_E$ and $\tilde{m}_I$ given in Eqs. (29) and (30), respectively.

Within each of the four regions determined by the above boundary lines, the system in Eqs. (32) is described by a different subset of parts of the piecewise-defined functions on its right-hand sides. We solve each system explicitly up to arbitrary constants of integration, determined by initial conditions. As the trajectory moves through one region and meets a boundary line, new initial conditions are determined by the intersection of the trajectory and the boundary, and the trajectory is then computed for the next region. Computing the new initial conditions as the trajectory crosses each boundary requires solving for the roots of a transcendental function, for which we employ matlab’s built-in root-finding algorithm, fzero. Proceeding in this manner through all four regions, the solution trajectory can be determined analytically, except for these four solves.

Specifically, we consider the two-dimensional system in Eq. (32) with fixed slow inhibitory effective conductance $g^S$, and begin with an initial value for the fast inhibitory effective conductance, $g^F(0) = l$. For convenience, we consider that the initial condition for the solution lies on the inhibitory boundary line, $\tilde{m}_I = 0$, from Eq. (33b), and thus the excitatory effective conductance begins at the value given by

$$\begin{aligned} \begin{aligned} g^E(0)&= -\frac{1}{S_I^Ep_I^E} \left\{ f_I\nu + \left( 1-\tfrac{(A-1)(B+1)}{\log A}\right) S_{I}^Fp_{I}^Fl \right. \\&\quad \left. + \left( 1-\tfrac{(A-1)(X+1)}{\log A}\right) S_{I}^Sp_{I}^Sg^S +1 - \frac{A-1}{\log A} \right\} . \end{aligned} \end{aligned}$$

(34)

Starting at the above-described point, $(g^E(0), g^F(0))$, the trajectory is determined by the system of equations,

$$\begin{aligned} \sigma _E\frac{dg^E}{dt}&= -g^E , \end{aligned}$$

(35a)

$$\begin{aligned} \sigma _F\frac{dg^F}{dt}&= -g^F + \frac{1}{\tau \log A}\left[ \left( 1-\frac{(A-1)(B+1)}{\log A}\right) S_{I}^Fp_{I}^Fg^F \right. \nonumber \\&\quad \left. + \left( 1-\frac{(A-1)(X+1)}{\log A}\right) S_{I}^Sp_{I}^Sg^S+ S_I^Ep_I^Eg^E+f_I\nu +1 - \frac{A-1}{\log A} \right] , \end{aligned}$$

(35b)

in the first region through which it moves. The system in Eq. (35) has the explicit solution given by the equations

$$\begin{aligned} g^E(t)&= C_1e^{-t/\sigma _E} , \end{aligned}$$

(36a)

$$\begin{aligned} g^F(t)&= C_2e^{-\kappa t/\sigma _F}+ \frac{C_1S_I^Ep_I^Ee^{-t/\sigma _E} + f_I\nu +1 - \frac{A-1}{\tau \ln A}}{\tau \left( \kappa - \frac{\sigma _F}{\sigma _E}\right) \ln A} - \frac{ f_I\nu +1 - \frac{A-1}{\tau \ln A}}{\tau \frac{\sigma _E}{\sigma _F}\left( \kappa ^2 - \kappa \right) \ln A }, \end{aligned}$$

(36b)

where $\kappa = \left[ 1- \frac{S_{I}^Fp_{I}^F}{\tau \ln A} \left( 1-\frac{(A-1)(B+1)}{\log A}\right) \right]$, and $C_1$ and $C_2$ are constants of integration determined by the initial condition established above. The value of t at which this trajectory meets the excitatory boundary in Eq. (33a) is determined using matlab’s built-in root-finding algorithm, and the value of Eq. (36) at that point becomes the initial condition for the solution calculation in the next region. Similar systems of equations and solutions to these systems can be determined explicitly in the same fashion. The trajectory moves through the second region where the linearized firing-rate equations, $\tilde{m}_E$ and $\tilde{m}_I$ are nonzero, then through the third region where $\tilde{m}_E$ is nonzero and $\tilde{m}_I = 0$, and finally moves through the fourth region where $\tilde{m}_E= \tilde{m}_I = 0$, until the trajectory reaches the first region again. At each point where the trajectory meets a boundary line, new initial conditions are determined. Explicit systems that determine the dynamics of the solution trajectory with explicit solutions, up to constants of integration, determined later by initial conditions, for each region of the two-dimensional $g^E$-$g^F$-plane, can be written down. (See ref. Pyzza 2015 for the complete set of systems and solutions). Note that the initial condition choice on the boundary line ${\tilde{m}}_I=0$ is made for convenience; for our argument below, it could be made along any of the boundaries dividing the regions in which Eqs. (32) are linear, or, in principle, along any sufficiently long curve segment transverse to the trajctories of the system.

The process by which we calculate the trajectory for a given value, $g^F(0) = l$, may be repeated for a range of different values of l. We then consider the value of the fast inhibitory effective conductance, $g^F(t)$, when the trajectory returns to the part of the boundary line, $\tilde{m}_I=0$, from which it originated, after one excursion away from it. In this way, we generate a Poincaré map, P(l), of this boundary line into itself. Investigating $P(l)-l$ for a range of values of l, we find that the map P(l) has a unique fixed point. In other words, there is a particular value of l such that the trajectory will begin and end at the same point, which indicates the existence of a unique limit cycle for a given parameter set.

We illustrate this process in Fig. 11, in which the slow effective conductance is fixed at the value $g^S = 0$ (i.e., no slow inhibition), representing the case with the largest oscillations. In Figs. 11a, b, for $l=l_1 = 0.28$ and $l=l_2 = 0.8$, respectively, we show two trajectory segments that progress counter-clockwise in the $g^E$-$g^F$-plane through the four regions determined by the boundary lines given in Eqs. (33), until they return to the portion of the line $\tilde{m}_I=0$ above its intersection with the line $\tilde{m}_E=0$. The blue, cyan, green, and red trajectory sub-segments are described by the solutions to the systems of equations for $g^E$ and $g^F$ given by the appropriate combinations of the piecewise portions of Eqs. (32a) and (32b) in each of these four regions, respectively. Note that for the first of two trajectories $P(l_1)>l_1$, and for the second $P(l_2)<l_2$. In Fig. 11c, d, we display the mapping of the interval $l_1<l<l_2$ along one excursion of the trajectories that begin on this interval, back onto a smaller interval inside $l_1<l<l_2$. The corresponding trajectory segments of a number of equidistant initial points are plotted in Fig. 11c, where we see that these trajectories are so strongly contracted towards one another that they are indistinguishable already on the line ${\tilde{m}}_E=0$. A graph of the difference, $P(l)-l$, between the final and initial $g^F$-values along the trajectory segments emerging from the interval $l_1<l<l_2$ on the line ${\tilde{m}}_I=0$ is displayed in Fig. 11d. The initial condition for which this difference vanishes gives the fixed point of the mapping l to P(l), and thus the unique limit cycle.

Fig. 11

The solid (dashed) black line represents the inhibitory (excitatory) boundary line given by $\tilde{m}_I = 0$ ($\tilde{m}_E= 0$). The colored, dashed curves illustrate the solution trajectory-segments through the four regions separated by these lines for the linearized firing rate system with $g^S$ frozen. For each such trajectory segment, the corresponding dot indicates its initial condition, and the pink dot the point at which the trajectory returns to the initial portion of the boundary line. Two example solution trajectories, one with $P(l_1)>l_1$ and another with $P(l_2)<l_2$, are shown in panels a and b. The mapping of the interval $l_1<l<l_2$ on the line ${\tilde{m}}_I=0$ along the analogous trajectory segments is shown in c. Note that a single cyan, green, red, and pink dot is plotted for all the trajectories as they are indistinguishable at those points due to the strong contraction. In panel d, the initial values for $g^F(0) = l$ are compared with the final values of $g^F=P(l)$ when the trajectory returns to the starting part of the boundary line. Their difference, $P(l)-l$, is plotted. The value l of the fixed point occurs where $P(l)-l=0$, i.e., at the intersection between the curve and the dashed line. Parameters found in Table 5

Table 5 Parameters used for visualizing boundary lines and computing solutions to the firing-rate system in Eq. (33) with $g^S$ held constant in Fig. 11