An exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks

TL;DR

This paper introduces an exact firing rate model for heterogeneous QIF neuron networks with chemical and electrical synapses, revealing how electrical coupling promotes synchrony and alters bifurcation scenarios.

Contribution

The paper presents a novel firing rate model that exactly captures the mean field dynamics of networks with both synapse types, including electrical coupling effects.

Findings

Electrical coupling induces neuronal synchrony via a supercritical Hopf bifurcation.

Electrical coupling transforms bifurcation scenarios, reducing persistent states.

Model predictions align with numerical simulations of large spiking neuron networks.

Abstract

Chemical and electrical synapses shape the dynamics of neuronal networks. Numerous theoretical studies have investigated how each of these types of synapses contributes to the generation of neuronal oscillations, but their combined effect is less understood. This limitation is further magnified by the impossibility of traditional neuronal mean-field models ---also known as firing rate models, or firing rate equations--- to account for electrical synapses. Here we introduce a novel firing rate model that exactly describes the mean field dynamics of heterogeneous populations of quadratic integrate-and-fire (QIF) neurons with both chemical and electrical synapses. The mathematical analysis of the firing rate model reveals a well-established bifurcation scenario for networks with chemical synapses, characterized by a codimension-2 Cusp point and persistent states for strong recurrent excitatory coupling. The inclusion of electrical coupling generally implies neuronal synchrony by virtue of a supercritical Hopf bifurcation. This transforms the Cusp scenario into a bifurcation scenario characterized by three codimension-2 points (Cusp, Takens-Bogdanov, and Saddle-Node Separatrix Loop), which greatly reduces the possibility for persistent states. This is generic for heterogeneous QIF networks with both chemical and electrical coupling. Our results agree with several numerical studies on the dynamics of large networks of heterogeneous spiking neurons with electrical and chemical coupling.

Equations

$$\tau {\dot{V}}_j=V_j^2+{\eta }_j+g(v-V_j)+J\tau s,$$

$$v=\frac{1}{N}\sum _{k=1}^N{V}_k,$$

$$s(t)=\frac{1}{N}\sum _{j=1}^N\frac{1}{{\tau }_s}\underset{t-{\tau }_s}{\overset{t}{\int }}\sum _k\delta \left(t^{\prime }-t_j^k\right)d{t}^{\prime },$$

$$f_i=\frac{1}{\tau \pi }\sqrt{{\eta }_i}.$$

$$\bar{\eta }=\frac{{\eta }_1+{\eta }_2}{2}.$$

$$L_{\Delta ,\bar{\eta }}(\eta )=\frac{1}{\pi }\frac{\Delta }{(\eta -\bar{\eta }{)}^2+{\Delta }^2},$$

$$\tau \dot{r}$$

$$\tau \dot{v}$$

$$\rho (V,t)=\frac{1}{\pi }\frac{\pi \tau r(t)}{{\left[V-v(t)\right]}^2+{\left[\pi \tau r(t)\right]}^2}.$$

$$\pi \tau r=\sqrt{v^2+\bar{\eta }}.$$

$${\left(\frac{\bar{\eta }}{\Delta }\right)}_{\text{H}}=\frac{4\Delta }{g^2}-\frac{g^2}{16\Delta }-\frac{2J}{\pi g},$$

$$f_{\text{H}}=\frac{1}{\pi \tau }\sqrt{\bar{\eta }+\frac{\Delta }{\pi }\frac{J}{g}}.$$

$$g_H\approx \frac{2\Delta }{\sqrt{\bar{\eta }}}.$$

$$\bar{f}=\frac{\sqrt{\bar{\eta }}}{\pi \tau }.$$

$$r_{\ast }=\Phi (\bar{\eta }+J\tau r_{\ast }),$$

$$\tau {\partial }_t\rho +{\partial }_V\left[(V^2+\eta +g(v-V)+J\tau r)\rho \right]=0,$$

$$v(t)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\rho (V|\eta ,t)V{L}_{\Delta ,\bar{\eta }}(\eta )dVd\eta .$$

$$\rho (V|\eta ,t)=\frac{1}{\pi }\frac{x(\eta ,t)}{{\left[V-y(\eta ,t)\right]}^2+x(\eta ,t{)}^2},$$

$$\tau \dot{x}(\eta ,t)$$

$$\tau \dot{y}(\eta ,t)$$

$$\tau\partial_{t}w(\eta,t)=i\big{[}\eta-w(\eta,t)^{2}+J\tau r\big{]}+g\big{[}iv-w(\eta,t)\big{]}$$

$$x(\eta ,t)=\pi \tau r(\eta ,t).$$

$$r(t)=\frac{1}{\tau \pi }{\int }_{-\infty }^{\infty }x(\eta ,t)L_{\Delta ,\bar{\eta }}(\eta )d\eta .$$

$$y(\eta ,t)=\text{p.v.}{\int }_{-\infty }^{\infty }\rho (V|\eta ,t)VdV.$$

$$v(t)={\int }_{-\infty }^{\infty }y(\eta ,t)L_{\Delta ,\bar{\eta }}(\eta )d\eta .$$

$$\pi \tau r(t)+iv(t)=w(\bar{\eta }-i\Delta ,t),$$

$$\tau \dot{r}$$

$$\tau \dot{v}$$

$$\rho (V,t)=\frac{1}{\pi }\frac{\pi \tau r(t)}{[V-v(t){]}^2+{\pi }^2{\tau }^{2}r(t{)}^2},$$

$$g[Q(t)-v(t)],$$

$$Q(t)=\frac{i}{2}\sum _{m=1}^{\infty }\frac{{\rho }^{m+1}-{\rho }^{m-1}}{\rho +1+ϵ}\left[z^m-{\bar{z}}^m\right]$$

$$\pi r+iv=\frac{1-\bar{z}}{1+\bar{z}}.$$

Full text

I Introduction

Collective oscillations and synchrony are prominent features of neuronal circuits, and are fundamental for the well-timed coordination of neuronal activity. Such oscillations are profoundly shaped by the presence of chemical synapses Wang (2010). An increasing number of experimental studies indicate both the prevalence and functional importance of electrical synapses (formed by gap junctions between neurons) in many diverse regions of central nervous systems, especially in inhibitory interneurons Nagy et al. (2018); Connors (2017); Traub et al. (2018). Electrical synapses participate in mediating synchronization of neuronal network activity Bennett and Zukin (2004); Connors and Long (2004), suggesting that electrical interaction may be interrelated with the generation of oscillations via chemical transmission.

The mechanisms by which chemical synapses mediate large-scale synchronous activity have been extensively investigated, see e.g. Wang (2010); Whittington et al. (2000). However, only a few studies addressed the synchronization of large networks in which neurons are not only interacting via excitation and/or inhibition, but also via electrical synapses  Kopell and Ermentrout (2004); Pfeuty et al. (2007); Ermentrout (2006); Viriyopase et al. (2016); Holzbecher and Kempter (2018); Guo et al. (2012); Tchumatchenko and Clopath (2014); Lewis and Rinzel (2003); Chow and Kopell (2000); Ostojic et al. (2009); Pfeuty et al. (2003); Coombes (2008); Mancilla et al. (2007). This limited theoretical progress for networks of electrically coupled neurons, compared to chemically coupled networks, is magnified due to the technical challenges faced when developing simplified mean field models —often called firing rate models, or firing rate equations (FRE)— for networks involving electrical synapses. While firing rate models turned out to be very useful to explain key aspects of the dynamics of spiking neuron networks with chemical synapses  Wilson and Cowan (1972); Dayan and Abbott (2001); Ermentrout and Terman (2010); Hopfield (1984); Mongillo et al. (2008); Ben-Yishai et al. (1995); Hansel and Sompolinsky (1998); Tabak et al. (2000); Rankin et al. (2015); Martí and Rinzel (2013); Tsodyks M. (1998); Roxin et al. (2005); Roxin and Montbrió (2011), it remains an open question whether there are similar simplified mean field theories for networks involving electrical interactions.

Recently, a novel method has been found to exactly derive FRE for populations of heterogeneous quadratic integrate-and-fire (QIF) neurons with chemical coupling Montbrió et al. (2015). The method, related to the so-called Ott-Antonsen ansatz Ott and Antonsen (2008, 2009); Ott et al. (2011); Pikovsky and Rosenblum (2008); Marvel et al. (2009); Pikovsky and Rosenblum (2011); Pietras and Daffertshofer (2016), allows to obtain exact, low-dimensional firing rate equations for ensembles of QIF neurons, see also Luke et al. (2013); So et al. (2014); Laing (2014). The FRE for QIF neurons have been used to investigate numerous problems regarding the dynamics of networks of chemically coupled QIF neurons  Pazó and Montbrió (2016); Ratas and Pyragas (2016); Roulet and Mindlin (2016); Devalle et al. (2018, 2017); Esnaola-Acebes et al. (2017); Ratas and Pyragas (2017); Dumont et al. (2017); Byrne et al. (2017); Laing (2018); Schmidt et al. (2018); Ratas and Pyragas (2018); di Volo and Torcini (2018); Dumont and Gutkin (2019); Akao et al. (2019); Bi et al. (2019); Coombes and Byrne (2019); Keeley et al. (2019); Byrne et al. (2019); Boari et al. (2019). Remarkably, previous work has also sought to apply this approach to networks with both chemical and electrical coupling Laing (2015). However, in Laing (2015), the electrical coupling has been treated by making use of an approximation which renders the resulting FRE analytically intractable. We build on this previous work and derive a set of FRE for networks with chemical and electrical coupling, but without the need for any approximation. The resulting system is not only analytically tractable but also allows, in a unified framework, for carrying out a complete analysis of the possible dynamics and bifurcations of networks with mixed chemical and electrical synapses. In Appendix B we show that our exact FRE are recovered by appropriately relaxing the approximation invoked in Laing (2015).

The structure of the paper is as follows: In Section II, we describe the spiking neuron network under investigation, and briefly illustrate the impact of electrical coupling in the dynamics of two nonidentical QIF neurons. In Section III, we introduce the FRE corresponding to the thermodynamic limit of the QIF network. The detailed derivation is performed in Appendix A. In Section IV, we perform a comparative analysis of the fixed points and their bifurcations in networks with electrical coupling vs. networks with chemical coupling. Finally, we investigate the dynamics of a QIF network with both electrical and chemical synapses and demonstrate that the presence of electrical coupling critically determines the bifurcation scenario of the neuronal network. Finally, we discuss our results in Section V.

II Quadratic integrate-and-fire neurons with electrical and chemical synapses

We consider a large population of globally electrically and chemically-coupled QIF neurons, with membrane potentials $\{V_{j}\}_{j=1,\dots,N}$ and $N\!\gg\!1$. Their dynamics reads

[TABLE]

where $\tau$ denotes the cells’ common membrane time constant, and parameter $\eta_{j}$ represents an external input current flowing into cell $j$. To model the action potential, the continuous dynamics Eq. (1) is supplemented by a discrete resetting rule. Here, we assume that if $V_{j}$ reaches infinity, neuron $j$ emits a spike and its membrane potential is reset to minus infinity 111In our numerical simulations (Euler scheme, $\delta t=10^{-4}$), the resetting rule was applied as follows: When $V_{j}\geq 100$, the membrane voltage is held at $V_{j}$ for a time interval $\tau/V_{j}$. Then, a spike is emitted, and the voltage is reset and kept at $-V_{j}$ for a subsequent interval $\tau/V_{j}$. In Figs. 4 and 5, to numerically evaluate the mean membrane potential $v$, the population average is computed discarding those neurons with $|V_{j}|\geq 100$. . The mean membrane voltage

[TABLE]

to which all cells are diffusively coupled with strength $g\geq 0$, mediates the electrical coupling. The constant $J$ quantifies the coupling strength of chemical synapses. The coupling via chemical synapses is mediated by the mean synaptic activation function

[TABLE]

where $t_{j}^{k}$ denotes the time of the $k$-th spike of the $j$-th neuron, $\delta(t)$ is the Dirac delta function, and $\tau_{s}$ is a synaptic time constant  222We used $\tau_{s}=10^{-2}$ms, and $\tau=10$ ms. The instantaneous firing rates shown in Fig. 4 were obtained by binning time and counting spikes within a sliding time window of size $\delta t=2.5\times 10^{-2}$ ms.. The synaptic weight $J$ can be positive or negative depending on whether the chemical synapses are excitatory or inhibitory, respectively

In the absence of coupling, $J=g=0$, the QIF neurons are either quiescent ($\eta_{i}<0$), or oscillatory ($\eta_{i}>0$) with frequency

[TABLE]

These two dynamical regimes of individual neurons are connected by a saddle-node on the invariant circle (SNIC) bifurcation, which occurs when $\eta_{i}=0$, with $f_{i}=0$.

Electrical coupling tends to equalize the membrane potentials of the neurons they connect and may favor synchrony. Yet, if a large fraction of cells in the network is quiescent, gap junctions may suppress oscillations and neural synchrony. Next we illustrate this phenomenon for two nonidentical QIF neurons that are coupled via a gap junction 333This phenomenon has been termed ‘Aging transition’ in the literature Daido and Nakanishi (2004); Pazó and Montbrió (2006).. The results of this analysis will later be useful to understand some aspects of the dynamics of a large network of electrically coupled QIF neurons.

II.1 Strong coupling limit of two electrically coupled QIF neurons

We consider a network of $N=2$ nonidentical QIF neurons with dynamics Eq. (1). The neurons are coupled via gap junctions only, i.e. $J=0$ but $g>0$. We are interested in the strong coupling limit $g\gg 0$ when $\eta_{1}>0$ and $\eta_{2}<0$. In Fig. 1 we depict the corresponding time series of cell 1 (panel a) and cell 2 (panel b). Black thin curves correspond to the dynamics of the uncoupled ($g=0$) cells: cell 1 fires periodically, while cell 2 remains quiescent. When the neurons are electrically coupled (red curves), the membrane voltage of cell 2 displays a series of so-called ‘spikelets’ 444 Note that these spikelets depend on the shape of the presynaptic spike, and thus on the particular neuron model considered.. Moreover, the electrical interaction brings cell 1 closer to its firing threshold and, hence, its frequency $f_{1}$ is reduced. When $g$ is increased further, cell 1 becomes quiescent (blue thick curves).

Although analyzing the dynamics of the two cells for arbitrary coupling strength $g$ is a challenge, there exists a simple and general result valid in the large $g$ limit, and of relevance for the large-$N$ analysis carried out below. Indeed, for large $g$, the dynamics of the $N=2$ network simply depends on the sign of the mean current Pazó and Montbrió (2006)

[TABLE]

For $\bar{\eta}>0$, the quiescent cell eventually becomes self-oscillatory as $g$ is increased from zero. By contrast, for $\bar{\eta}<0$, the oscillatory cell eventually turns quiescent in the strong coupling limit; see the blue lines in Fig. 1 555The value $\bar{\eta}=0$ determines a boundary separating network oscillations and quiescence, cf. Eq (6) with $p_{c}^{\infty}=1/2$ in Pazó and Montbrió (2006)..

III Firing Rate Model

In the following, we introduce the FRE corresponding to the thermodynamic limit of Eqs. (1). The detailed derivation of the model closely follows the lines of Montbrió et al. (2015) and is given in Appendix A.

For $N\to\infty$, one can drop the indices in Eq. (1) and define a density function $\rho$ such that $\rho(V|\eta,t)~{}dV$ denotes the fraction of neurons with membrane potentials between $V$ and $V\!+\!dV$ and parameter $\eta$ at time $t$. In the limit of instantaneous synaptic processing, i.e. for $\tau_{s}\to 0$, Eq. (2) reduces to $s(t)=r(t)$ with $r(t)$ being the population-mean firing rate. If the external currents are distributed according to a Lorentzian distribution centered around $\eta=\bar{\eta}$ with half-width $\Delta$,

[TABLE]

we find that the asymptotic mean-field dynamics evolves according to the following FRE 666See Appendix B for the comparison of the firing rate model Eqs.(5) with the FRE derived in Laing (2015).

[TABLE]

The variables $r$ and $v$ are the mean firing rate and mean membrane potential, respectively. They determine the total voltage density for the network Eq. (1), which turns out to be a Lorentzian distribution centered at $v(t)$ and of half-width $\pi\tau r(t)$,

[TABLE]

The structure of the FRE Eqs. (5) reveals an interesting feature: Electrical coupling is solely mediated by the firing rate through the negative feedback term $-gr$ in the $r$-dynamics Eq. (5a), and not by membrane potential differences  777This might be understood as follows: The evolution equation for the mean membrane potential $v$ is obtained summing up the $N$ differential equations Eq. (1): $\tfrac{1}{N}\sum_{i=1}^{N}\dot{V}_{i}=\dot{v}=\tfrac{1}{N}\sum_{i=1}^{N}(V_{i}^{2}+\eta_{i})+Js+g(v-\tfrac{1}{N}\sum_{i=1}^{N}V_{i})$. Although this is not a closed equation for $v$ and $r$, one finds that the last term —corresponding to diffusive coupling— cancels to zero.. That is, electrical coupling leads to a narrowing of the voltage distribution Eq. (6), i.e. a decrease in firing rate. This confirms our initial sketch that electrical coupling tends to equalize the neurons’ membrane potentials and, under suitable conditions, this may promote synchrony. By contrast, chemical coupling shifts the center of the distribution Eq. (6) of voltages via the feedback term $Jr$ in the $v$-dynamics Eq. (5b). The following phase plane and bifurcation analysis of the FRE (5) allows for understanding the collective dynamics of the QIF network.

IV Analysis of the Firing Rate Equations

IV.1 Electrical vs. chemical coupling

In the absence of chemical coupling, our previous discussion of the case $N\!=\!2$ hints at two distinct dynamical regimes for positive and negative values of $\bar{\eta}$. With respect to the fixed points $(r_{*},v_{*})$ of the FRE (5) for $J=0$, we find the $v$-nullcline to be

[TABLE]

Note that if $\bar{\eta}$ is negative, there exists a range of ‘forbidden’ values of $v$. Fig. 2(a) shows the nullclines for $\bar{\eta}<0$ and for different values of the ratio $g/\sqrt{\Delta}$. Since the majority of the neurons are quiescent, an increase in coupling strength $g$ causes active neurons to reduce firing, which leads to a progressive decrease of the firing rate $r_{*}$. By contrast, in Fig. 2(b) the majority of the cells are self-oscillatory, $\bar{\eta}>0$, and strong electrical coupling forces quiescent neurons to fire. This yields an increase of $v_{*}$. Interestingly, the firing rate $r_{*}$ is a non-monotonic function of $g/\sqrt{\Delta}$: While $v_{*}$ remains negative, the voltages are pushed to subthreshold values, decreasing the firing rate. This behavior is reverted when $v_{*}$ becomes positive and all voltages are pushed towards values above the firing threshold. The different behaviors of Eqs. (5) with electrical coupling for positive and negative values of $\bar{\eta}$ are clearly revealed in the corresponding bifurcation diagrams shown in Figs. 3(a,c).

The case of networks with only chemical coupling, $g=0$, is simpler Montbrió et al. (2015). The bifurcation diagram depicted in Fig. 3(d) shows that $v_{*}$ remains always negative and converges asymptotically to zero as $J$ increases. The firing rate $r_{*}$, depicted in Fig. 3(b), also increases with $J$. For $\bar{\eta}<0$ and strong recurrent excitatory coupling, the system undergoes a cusp bifurcation and two saddle-node (SN) bifurcations are created. This implies the existence of a parameter regime where a persistent, high-activity state (stable focus) coexists with a low-activity state (stable node) —see Fig. 7(a), and  Montbrió et al. (2015). This coexistence between persistent and low-activity states also occurs in networks with electrical synapses, but it is located in a very small region of parameters as we show below, see Fig. 6(b).

We next explore the linear stability of the fixed points of Eqs. (5), see also Appendix C. We find that a Hopf bifurcation occurs along the boundary

[TABLE]

with frequency

[TABLE]

The Hopf boundary Eq. (7) is depicted in red in the phase diagrams of Figs. 6, 7. Note that $\bar{\eta}/\Delta\to+\infty$ as $g\to 0$ according to Eq. (7), which indicates that electrical coupling is a necessary ingredient for the Hopf bifurcation to exist  888For non-instantaneous inhibitory synapses, oscillations emerging through a Hopf bifurcation are also encountered for weak coupling and weak heterogeneity Devalle et al. (2017, 2018). These oscillations are often referred to as ‘Interneuronal Gamma Oscillations’, ING Whittington et al. (2000). To keep our analysis simple, here we do not consider ING oscillations, since Eqs. (5) become higher-dimensional and phase plane analysis is no longer possible in this case..

To confirm the presence of collective oscillations in the original network of electrically coupled QIF neurons with dynamics Eq. (1), we carried out numerical simulations and compared them with those of the FRE (5). Fig. 4 shows the time series of the firing rate in the full and in the reduced system, which display a very good agreement. In panel (a) we considered a network with electrical coupling only. The frequency of the oscillations, $f\approx 30.1$ Hz, is close to the theoretical value at criticality, given by Eq. (8): $f_{\text{\scriptsize{H}}}=100/\pi\approx 31.8$ Hz. Therefore, in absence of chemical coupling and near the Hopf bifurcation, the frequency of the oscillations is almost independent of the coupling strength $g$ and closely follows Eq. (8). To further test the validity of Eq. (8) far from criticality, we numerically evaluated the frequency of the limit cycle of the FRE (5) (black solid line, Fig. 5) as the the coupling strength $g$ is increased from the Hopf bifurcation (at $g_{H}\approx 1.8$). The black dotted line corresponds to the Hopf frequency Eq. (8). We find that the frequency of the limit cycle remains close to this for a broad range of $g$ values.

Hopf instability in networks of electrically coupled QIF neurons occurs like the transition to synchronization in the Kuramoto model of coupled phase oscillators Kuramoto (1975). Considering $J=0$, we find the main features of the Kuramoto transition to collective synchronization: ($i$) In the limit of weak electrical coupling $g$, the Hopf boundary Eq. (7) can be written as

[TABLE]

For $\bar{\eta}=1$, Eq. (9) coincides with Kuramoto’s critical coupling for synchrony. ($ii$) As previously discussed, macroscopic oscillations emerge with a frequency determined by the most likely value of the natural frequencies in the network, see Eq. (3). For the case of the Lorentzian distribution of currents, Eq. (4), the most likely value of the frequency is

[TABLE]

($iii$) The Hopf bifurcation is always supercritical; cf. Appendix D. Taken together, for $\bar{\eta}>0$ and given a certain level of heterogeneity $\Delta$, synchronization occurs —at a critical coupling approximately given by Eq. (9)— with the nucleation of a small cluster of oscillators with natural frequencies Eq. (3) near $\bar{f}$. As electrical coupling $g$ is further increased, more and more oscillators become entrained to the frequency $\bar{f}$, resulting in a continuous and monotonous increase in the amplitude of the oscillations. This transition is in contrast to that of networks with inhibitory coupling and synaptic kinetics and/or delays, where synchrony is only achieved for weak heterogeneity and weak coupling, see, e.g., Devalle et al. (2017, 2018).

The phase diagram depicted in Fig. 6 characterizes the dynamics of the firing rate model Eq. (5) with only electrical coupling. The red curve corresponds to the Hopf bifurcation line given by Eq. (7). According to Eq. (8), the frequency of the collective oscillations approaches zero as $\bar{\eta}\to 0$. This indicates that the Hopf line ends in a Takens-Bogdanov (TB) bifurcation at $\bar{\eta}=0$, see Fig. 6(b). At this codimension-2 point, the Hopf boundary tangentially meets a SN bifurcation and a homoclinic bifurcation. The homoclinic line moves parallel to the Hopf line for a while, it makes a sharp backward turn and then tangentially joins onto the upper branch of the SN bifurcation curve (two branches of SN bifurcations are created at the Cusp point), at a saddle-node-separatrix-loop (SNSL) point. At this point the SN boundary becomes a SNIC boundary that, together with the Hopf and homoclinic lines, encloses the region of synchronization (Sync) featuring collective oscillations. Note that in Fig. 6(b) we encounter a very small region of bistability between a Low-Activity State (LAS, node) and a persistent state (focus). Electrical coupling destabilizes the persistent state almost immediately after the SN line, leading to another small region of bistability between LAS and a small amplitude limit cycle (Sync) —which disappears in the homoclinic bifurcation.

Finally, the SNIC curve asymptotically approaches $\bar{\eta}=0$ as $g/\sqrt{\Delta}\to\infty$ (as suggested by the $N=2$ analysis in Section II.1). In this limit, all neurons are strongly coupled ($g\to\infty$) and/or are nearly identical ($\Delta\to 0$) so that they behave as a single QIF neuron with input current $\bar{\eta}$ 999 Exactly the same scenario is found in systems of globally, sine coupled ‘active rotators’. Such systems are, in fact, closely related to the dynamics (1) with $J=0$ Sakaguchi et al. (1988); Zaks et al. (2003); M.Childs and Strogatz (2008); Lafuerza et al. (2010)..

IV.2 Networks with both chemical and electrical coupling

We finally analyze the dynamics of a population of QIF neurons with mixed, chemical and electrical synapses. Fig. 7(a) presents the possible dynamical regimes of a population with chemical synapses only, $g=0$. In contrast to networks with pure electrical coupling, where the bifurcation scenario is determined by the presence of three codimension-2 points, cf. Fig. 6, here there is only a Cusp point, see also Montbrió et al. (2015). This entails the presence of a persistent state (focus) coexisting with an asynchronous, LAS (node) within the cusp-shaped region in the top-left corner of Fig. 7(a). Additionally, the dashed line indicates that the asynchronous state is of focus type in a vast region of parameters for excitatory coupling, and always for inhibitory coupling.

Including electrical coupling, $g>0$, yields the Hopf bifurcation given by Eq. (7), which joins onto the lower branch of the SN bifurcation curve at a TB point, see Figs. 7(b,c). Hence, the bifurcation scenario for networks with electrical and chemical synapses matches that for networks with electrical synapses only: Similar to Fig. 6, the Hopf line cuts through the cusp-shaped region —the TB bifurcation demarcates the point where the Hopf boundary and the lower SN line intersect. Then, due to the presence of electrical coupling, the persistent state becomes only stable in a small parameter region confined between the Hopf and the lower SN line, see Fig. 7(b). As electrical coupling is increased, the TB point approaches the Cusp bifurcation, which results in an even smaller range of parameters for which the persistent state is stable. This agrees with numerical results using large networks of noisy, conductance-based and QIF neurons, and has been hypothesized to be a possible reason why electrical synapses are rarely found between excitatory neurons Ermentrout (2006).

Returning to the analysis of the FRE (5), we find for low values of $g$ that synchronization emerges predominantly for excitatory coupling, $J>0$, see Fig. 7(b). As electrical coupling is increased, the Sync region extends to the inhibitory region, $J<0$, and to larger values of $\bar{\eta}$ —note that in this coupling regime the emergence of collective oscillations mainly occurs via a SNIC bifurcation for excitation and via a Hopf bifurcation for inhibition, see Fig. 7(c). For even larger electrical coupling, the TB point moves further into the inhibitory region. That is, for strong electrical coupling the $J$-coordinate of the TB bifurcation rapidly decreases towards minus infinity whereas the other coordinate stays relatively close to the $\bar{\eta}=0$ axis 101010 It can be shown that the TB point behaves proportional to $(\tilde{\eta},\tilde{J})_{\text{\scriptsize{TB}}}\propto(g^{2},-g^{3})$ in the limit $g\to\infty$. . The SNIC bifurcation tilts towards a vertical line close to the $\bar{\eta}=0$ axis, because strong electrical coupling coerces all neurons to behave as a single QIF neuron with common input $\eta=\bar{\eta}$. Then, the SNIC bifurcation becomes the only transition between the two possible dynamical regimes, asynchrony or synchrony, see Fig. 7(d).

Fig. 4 shows how the addition of inhibitory coupling into a network with only electrical synapses degrades synchrony —parameters used in Fig. 4 correspond to the symbols shown in Fig. 7(c). The presence of inhibition clearly slows down the oscillations, as predicted by Eq. (8).

Although Eq. (8) is strictly valid only at the Hopf bifurcation, it is a good estimate of the frequency of the oscillations of the FRE (5). Fig. 5 depicts the comparison between Eq. (8) as a function of $g$ (dotted lines), with the actual frequencies numerically obtained using the FRE (5) (solid lines) and the QIF network Eq. (1) (symbols). In excitatory networks, the oscillations already emerge for weak values of $g$. In contrast, synchronizing inhibitory networks requires a much larger value of $g$, i.e. inhibition does not promote synchronization. Remarkably, only the presence of chemical coupling allows the frequency of the oscillations to deviate from $\bar{f}$, see Eqs. (8,10). Oscillations emerge with $f>\bar{f}$ for excitation and with $f<\bar{f}$ for inhibition, while they remain $f\approx\bar{f}$ for networks with only electrical coupling. As $g$ increases, the effects of chemical coupling are gradually washed out since an increasing number of neurons are entrained by electrical coupling to the most-likely frequency of the uncoupled network, $\bar{f}$. This dependence is well described by Eq. (8). Finally, since the level of heterogeneity $\Delta$ degrades synchrony, in Eq. (8) this term favors the deviation of the frequency from $\bar{f}$, and compensates for the homogenizing effect of electrical coupling. In the limit of identical neurons $\Delta\to 0$, the effects of instantaneous chemical coupling on the frequency vanish, since neurons synchronize in-phase and, at the instant of firing, all neurons become refractory.

V Conclusion and discussion

Firing rate models are very useful tools for investigating the dynamics of large networks of spiking neurons that interact via excitation and inhibition, see e.g. Wilson and Cowan (1972); Dayan and Abbott (2001); Ermentrout and Terman (2010); Hopfield (1984); Mongillo et al. (2008); Ben-Yishai et al. (1995); Hansel and Sompolinsky (1998); Tabak et al. (2000); Rankin et al. (2015); Martí and Rinzel (2013); Tsodyks M. (1998); Roxin et al. (2005); Roxin and Montbrió (2011). Remarkably, using a recently proposed approach to derive exact firing rate equations for networks of excitatory and/or inhibitory QIF neurons Montbrió et al. (2015); Luke et al. (2013), Laing has found that electrical synapses can also be incorporated in the framework of firing rate models Laing (2015). Yet, the FRE in Laing (2015) are not exact, and their mathematical form makes the analysis intractable.

Here we showed that the FRE corresponding to a network of QIF neurons with both chemical and electrical synapses can be exactly obtained, without the need for any approximation. Much in the spirit of firing rate models, the resulting FRE (5) are simple in form and highly amenable to analysis. Moreover, in Appendix B, we demonstrate that relaxing the approximation invoked in Laing (2015) the FRE derived by Laing simplify to our Eqs. (5).

At first glance, the mathematical form of Eqs. (5) already unveils two interesting features of electrical and chemical coupling, see also Eq. (6): ($i$) Chemical coupling tends to shift the center of the distribution of membrane potentials (given by $v$), while electrical coupling tends to reduce the width of the distribution (given by $r$), potentially promoting the emergence of synchronization; ($ii$) While in the original network of QIF neurons electrical coupling is mediated by membrane potential differences, at the mean-field level the electrical interaction is solely mediated by the mean firing rate $r$.

The mathematical analysis of the FRE (5) unravels how chemical and electrical coupling shape the dynamics of globally coupled populations of QIF neurons with Lorentzian heterogeneity. Some of our results were already reported in previous work, and confirm the value and the validity of the FRE (5). An important conclusion of our study is that the presence of electrical coupling, $g\neq 0$, generally implies the appearance of a supercritical Hopf bifurcation, see Eq. (7). This Hopf bifurcation meets a SN bifurcating line in a codimension-2 TB point, causing a drastic reduction of the region of bistability between low-activity and persistent asynchronous states, see Figs. 6 and 7. The Hopf bifurcation destabilizes the persistent state producing synchronous oscillations, which then are abolished via a homoclinic bifurcation. Previous studies of networks of excitatory and inhibitory neurons showed that synchrony often destroys persistent states Hansel and Mato (2001, 2003); Gutkin et al. (2001); Laing and Chow (2001); Kanamaru and Sekine (2003). Moreover, the generality of the bifurcation scenario of Figs. 6 and 7 —characterized by three codimension-2 points, TB, Cusp and SNSL—, is confirmed in previous studies analyzing closely related systems Sakaguchi et al. (1988); Zaks et al. (2003); M.Childs and Strogatz (2008); Lafuerza et al. (2010).

Networks of spiking neurons with strong excitatory coupling display robust persistent states. These states emerge at a Cusp bifurcation, see Fig. 7(a). Of particular relevance to our study is the work by Ermentrout Ermentrout (2006). He found that electrical coupling tends to synchronize neurons, and that this anihilates persistent states via the bifurcation scenario described in Figs. 6 and  7. Persistent activity may underlie important cognitive functions such as working memory, and has been suggested as a possible reason for the lack of electrical coupling between excitatory neurons Ermentrout (2006). According to Ermentrout (2006), ‘the main role for gap junctions is to encourage synchronization during rhythmic behavior. Synchrony, because it leads to a shared refractory period between neurons can lead to the extinction of persistent activity’.

The Hopf bifurcation is always supercritical. In Ostojic et al. (2009), Ostojic et al. analyzed the super- or sub-critical character of the Hopf bifurcation in networks of electrically coupled leaky integrate and fire (LIF) neurons. At variance with QIF neurons, LIF neurons do not have spikes and, hence, modeling electrical coupling requires an additional parameter Lewis and Rinzel (2003). This parameter enables one to adjust the shape of the spikelet elicited in the postsynaptic cell due to an action potential in the presynaptic cell. Ostojic et al. Ostojic et al. (2009) found that the Hopf bifurcation is supercritical when the spikelets are effectively excitatory, while inhibitory spikelets lead to subcritical Hopf bifurcations. For the QIF model, the spikelet elicited in a postsynaptic cell by the transmission of a presynaptic spike has a net excitatory effect —see Fig. 1(b)—, and hence our result that the Hopf bifurcation is supercritical is in agreement with the results in Ostojic et al. (2009). Yet, we note that our result also includes networks with chemical synapses, and not only networks with electrical synapses, as in Ostojic et al. (2009).

Another important result by Ostojic et al. Ostojic et al. (2009) is that electrical coupling can lead to oscillations even in the presence of strong heterogeneity. Our Eq. (7) is consistent with this. Kopell and Ermentrout Kopell and Ermentrout (2004) also investigated the robustness of synchrony against current heterogeneities in networks with both electrical and inhibitory synapses. They found that a small amount of electrical coupling, added to an already significant inhibitory coupling, can increase synchronization more than a very large increase in the inhibitory coupling. In Fig. 4 we show that increasing inhibition reduces the amplitude of the oscillations in a network with $g\neq 0$. In addition, Fig. 7 shows that, for a given value of $g$, increasing inhibition leads to asynchrony. This level of inhibition increases with electrical coupling, in line with the results in Kopell and Ermentrout (2004). Two studies Ostojic et al. (2009); Viriyopase et al. (2016) also investigated the frequency of the emerging oscillations in networks with electrical synapses. This frequency remains tied to the mean firing rate $f_{i}$ in the network (i.e. near $\bar{f}$), as our Eq. (8) suggests. In Fig. 5 we confirm that, in networks with only electrical synapses, the frequency of the oscillations remains near the most likely $f_{i}$-value: $\bar{f}$.

The result that the Hopf bifurcation is always supercritical, and that the frequency of the emerging oscillations is given by $\bar{f}$ evoke the paradigmatic synchronization transition in the Kuramoto model Kuramoto (1975). For weakly electrically-coupled networks, we find that the onset of oscillations occurs at the Kuramoto’s critical coupling for synchrony, Eq. (9). When considering chemical coupling, the frequency of the oscillations deviates from $\bar{f}$, increasing/decreasing for excitatory/inhibitory coupling. The intensity of this deviation depends on the ratio of chemical to electrical coupling, as Eq. (8) suggests. Fig. 5 confirms that strong electrical coupling overcomes the effect of excitation/inhibition onto the frequency of the oscillations, approximately as dictated by Eq. (8).

Together with the firing rate model derived in Laing (2015), the FRE (5) constitute a unique example of a firing rate model with both electrical and chemical coupling. The numerical simulations of the original QIF network Eq. (1) are in agreement with the FRE (5) —see Figs. 4, 5—, underlining the validity of the reduction method applied. Interestingly, the fixed points of the FRE (5) with chemical synapses ($g=0$, $J\neq 0$) can be cast in the form of a traditional firing rate model

[TABLE]

where $\Phi(x)=\sqrt{x+\sqrt{x^{2}+\Delta^{2}}}/(\sqrt{2}\pi\tau)$ is the so-called transfer function of the heterogeneous QIF network  Esnaola-Acebes et al. (2017); Devalle et al. (2017, 2018). The FRE (5) with electrical synapses ($g\neq 0$), however, cannot be written in the form of Eq. (11). Therefore, the link pointed by Eq. (11) between traditional firing rate models and Eqs. (5) is lost when electrical coupling is considered.

Acknowledgements

This work was supported by ITN COSMOS funded by the EU Horizon 2020 Research and Innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 642563. EM acknowledges support by the Spanish Ministry of Economy and Competitiveness under Project No. PSI2016-75688-P.

Appendix A: Derivation of the Firing Rate Equations

In the thermodynamic limit, $N\to\infty$, we drop the indices for the individual neuronal dynamics Eq. (1), and denote $\rho(V|\eta,t)dV$ as the fraction of neurons with membrane potentials between $V$ and $V+dV$, and parameter $\eta$ at time $t$. Accordingly, the parameter $\eta$ becomes a continuous random variable that is distributed according to a probability distribution function, which here is considered to be a Lorentzian $L_{\Delta,\bar{\eta}}(\eta)$ of half-width $\Delta$ and centered at $\bar{\eta}$, see Eq. (4). The conservation of the number of neurons leads to the continuity equation

[TABLE]

where we explicitly included the velocity given by the continuous equivalent of Eq. (1). We also defined the mean value of the membrane potential as

[TABLE]

Next, we consider the family of conditional density functions Montbrió et al. (2015)

[TABLE]

which are Lorentzian functions with time-dependent half-width $x(\eta,t)$, centered at $y(\eta,t)$. Substituting (A3) into the continuity equation (A1), we find that, for each value of $\eta$, variables $x$ and $y$ must obey two coupled equations,

[TABLE]

that can be written in complex form as

[TABLE]

where $w(\eta,t)\equiv x(\eta,t)+iy(\eta,t)$. For a particular value of $\eta$, the firing rate $r$ of the population of QIF neurons is related to the width $x$ of the Lorentzian ansatz (A3). Specifically, the firing rate $r(\eta,t)$ for each $\eta$ value at time $t$ is the probability flux at infinity: $r(\eta,t)=\rho(V\to\infty|\eta,t)\dot{V}(V\to\infty|\eta,t)$, which yields the identity

[TABLE]

Hence, integrating this quantity over the distributions of currents Eq. (4) provides the mean firing rate

[TABLE]

Likewise, we can link the center $y(\eta,t)$ of the Lorentzian ansatz Eq. (A3) with the mean of the (conditional) membrane potential via

[TABLE]

Note that the Lorentzian distribution does not have finite moments so that the integral in Eq. (A8) needs to be taken as the Cauchy principal value (i.e. $\text{p.v.}\int_{-\infty}^{\infty}\rho VdV=\lim_{R\to\infty}\int_{-R}^{R}\rho VdV$). Then, Eq. (A2) becomes

[TABLE]

The integrals in (A7,A9) can be evaluated closing the integral contour in the complex $\eta$-plane and using Cauchy’s residue theorem. The integrals must however be performed carefully, so that the variable $x(\eta,t)$ remains non-negative. To make the analytic continuation of $w(\eta,t)$ from real to complex-valued $\eta$, we define $\eta\equiv\eta_{r}+i\eta_{i}$. This continuation is possible into the lower half-plane $\eta_{i}<0$, since this guarantees the half-width $x(\eta,t)$ to remain non-negative: $\partial_{t}x(\eta,t)=-\eta_{i}>0$, at $x=0$. Therefore, we perform contour integration in Eq. (A7) and Eq. (A9) along the arc $|\eta|e^{i\vartheta}$ with $|\eta|\to\infty$ and $\vartheta\in(-\pi,0)$. This contour encloses one pole of the Lorentzian distribution Eq. (4). Then, we find that the firing rate and the mean membrane potential depend only on the value of $w$ at the pole of $L_{\Delta,\bar{\eta}}(\eta)$ in the lower half $\eta$-plane:

[TABLE]

As a result, we only need to evaluate Eq. (A5) at $\eta=\bar{\eta}-i\Delta$, and obtain a system of FRE composed of two ordinary differential equations as given in Eq. (5),

[TABLE]

in terms of the population-mean firing rate $r$ and the population-mean membrane potential $v$. Multiplying the Lorentzian ansatz Eq. (A3) by $L_{\Delta,\bar{\eta}}(\eta)$ and integrating over $\eta$, we finally obtain the total density of neurons Eq. (6) as

[TABLE]

where we again applied Cauchy’s residue theorem by using that the ansatz Eq. (A3) is analytic in the lower $\eta$-complex plane. Hence, the total density of the population of QIF neurons is a Lorentzian distribution centered at $v(t)$ and half-width $\pi\tau r(t)$, which evolves according to the FRE (5).

Appendix B: Connection between the FRE in Laing (2015) and Eq. (5)

The derivation of the FRE (5) is exact in the thermodynamic limit, and does not rely on any approximation. Here we show that the Eqs. (2.35&2.36) in Laing (2015) reduce to our Eq. (5) after adopting a limit in which the derivation performed in Laing (2015) becomes exact.

In contrast to our Eq. (5b), note that Eq. (2.36) in Laing (2015) contains a diffusive term,

[TABLE]

where the function $Q(t)$ is defined as

[TABLE]

with $0<\epsilon\ll 1$, and $\rho=\sqrt{2\epsilon+\epsilon^{2}}-1-\epsilon$. The variable $z$ in Eq. (B2) is the complex Kuramoto order parameter (the bar denotes complex conjugation), which is related to the variables $r$ and $v$ in the FRE (5) via the change of variables Montbrió et al. (2015)

[TABLE]

The parameter $\epsilon$, defined in Eq. (2.7) in Laing (2015), was used to approximate the mean voltage $v$, see also Ermentrout (2006). In the limit $\epsilon\to 0$ this approximation becomes exact, but this limit was not considered in Laing (2015). In consequence, to use the Eqs. (2.35&2.36) in Laing (2015), the infinite series Eq. (B2) was truncated after 100 terms, and the bifurcation analysis of the mean-field model could only be performed numerically.

Using the geometric series formula ($|z|<1$) and the transformation of variables Eq. (B3), we find

[TABLE]

Hence we have showed that the diffusive term Eq. (B1) identically vanishes when the mean-field reduction becomes exact (i.e. in the limit $\epsilon\to 0$), and the FRE in Laing (2015) reduce to Eqs. (5).

Appendix C: Bifurcation analysis of the Firing Rate Equations

The FRE (5) have five free parameters. The number of effective parameters can be reduced to three through non-dimensionalization, defining

[TABLE]

and rescaling variables as

[TABLE]

Then, the firing rate model becomes

[TABLE]

The fixed points $(\tilde{r}_{*},\tilde{v}_{*})$ of Eq. (C1) satisfy

[TABLE]

Linearization about the fixed points Eq. (C2) gives the eigenvalues

[TABLE]

For networks with only chemical synapses (i.e. $g=0$), the real part of the eigenvalues remains always negative (since $v_{*}<0$), and a Hopf bifurcation is not possible. However, chemical coupling has a direct influence on the real part of the eigenvalues Eq. (C3), and may produce oscillatory instabilities if the argument of the square root is a real number.

V.1 Hopf boundaries and Takens-Bogdanov point

The Hopf boundaries can be obtained when imposing Re($\lambda_{\pm})=0$ in Eq. (C3), which gives $\tilde{g}=4\tilde{v}_{*}$. Then, using Eq. (C2), we find

[TABLE]

Substituting Eq. (C4) in the $v$-fixed point equation Eq. (C1b), and solving for $\tilde{\eta}$, we obtain

[TABLE]

Solving Eq. (C4) for $\tilde{r}_{*}$ and substituting it into Eq. (C5) we obtain the Hopf boundaries in explicit form

[TABLE]

The frequency of the oscillations is given by the imaginary part of the eigenvalues Eq. (C3) at criticality that, using the fixed points of Eqs. (C1) and Eq. (C4), reduces to the explicit formula

[TABLE]

The frequency becomes zero at a Takens-Bogdanov (TB) point, when $\tilde{\eta}=-\tilde{J}/\tilde{g}$. Inserting this condition into Eq. (C6) we obtain the coordinates of the TB point

[TABLE]

see also Fig. 7. For $\tilde{J}=0$, the TB point is located at

[TABLE]

in the phase diagram Fig. 6.

V.2 Saddle-node boundaries

The boundaries of the saddle-node bifurcations are obtained by setting $\lambda_{\pm}=0$ in Eq. (C3), using Eq. (C2), and solving for $\tilde{g}$:

[TABLE]

Substituting (C10) in the $v$-fixed point Eq. (C1b), and solving for $\tilde{\eta}$, we obtain

[TABLE]

The saddle node boundaries are plotted in the $(\tilde{\eta},\tilde{g})$ phase diagram in Fig. 6. The same boundaries can be represented in the $(\tilde{\eta},\tilde{J})$ phase diagram when solving Eq. (C10) for $\tilde{J}$

[TABLE]

and replacing $\tilde{J}$ by Eq. (C12) in Eq. (C11) gives

[TABLE]

These saddle-node boundaries are shown in black for different values of $\tilde{g}$ in Fig. 7.

V.3 Focus-Node boundaries

The boundaries in the phase diagram Fig. 7 in which the stable asynchronous state changes from Focus to Node can be obtained in parametric form equating the square root in Eq. (C3) to zero. This gives

[TABLE]

Substituting $\tilde{J}_{FN}$ into the $v$-fixed point Eq. (C1b), and using Eq. (C2) we find

[TABLE]

For networks without chemical coupling, $\tilde{J}=0$, the Focus-Node boundary can be obtained in explicit form

[TABLE]

This is the dashed boundary depicted in Fig. 6.

Appendix D: Small-amplitude equation

near the Hopf bifurcation

In this Appendix we derive the small amplitude equation near the Hopf bifurcation, and show that the Hopf bifurcation is always supercritical. The derivation is performed using multiple-scales analysis, see e.g. Kuramoto (1984). We first expand the solution

[TABLE]

in powers of a small parameter $\epsilon\ll 1$, about a fixed point $(r_{0},v_{0})$ of Eqs. (C1) at the Hopf bifurcation. In addition, we introduce the deviation from the Hopf bifurcation Eq. (C6) of parameter $\tilde{\eta}$ as

[TABLE]

where $\chi$ determines the sign of the deviation. Finally, we define the slow time

[TABLE]

Then, the time differentiation is transformed as

[TABLE]

Plugging Eqs. (D1,D2,D4) into Eq. (C1) gives

[TABLE]

where

[TABLE]

and

[TABLE]

Critical eigenvectors

Using Eqs. (C3,C4), we define the critical frequency as

[TABLE]

so that the matrix Eq. (D6) can be written as

[TABLE]

The right-eigenvector of Eq. (D9) is

[TABLE]

Imposing the condition $\mathbf{u}_{L}\mathbf{u}_{R}=1$, the left-eigenvector of $L_{0}$ is

[TABLE]

Analysis of multiple scales

At order $\epsilon$, Eq. (D5) is

[TABLE]

This system of differential equations has a general solution

[TABLE]

which is the so-called neutral solution.

At order $\epsilon^{2}$, Eq. (D5) is

[TABLE]

Substituting the neutral solution Eq. (D13) into Eq. (D7) we find

[TABLE]

Next we use the following ansatz

[TABLE]

and substitute it into Eq. (D14). We find

[TABLE]

Substituting Eqs. (D13,D15) into the cubic term $N_{3}$ in Eq. (D7), we find

[TABLE]

The solvability condition at order $\epsilon^{3}$ is

[TABLE]

Substituting Eqs. (D11,D13,D16) into Eq. (D17), we find the amplitude equation

[TABLE]

with the coefficients

[TABLE]

where we used Eq. (D8) to express $a-d$ in terms of $\tilde{J}$ and $\tilde{g}$. Defining the amplitude $R$ and the phase $\Psi$ via

[TABLE]

one may alternatively write Eq. (D18) as

[TABLE]

where primes refer to differentiation with respect to $T$. An oscillatory solution with amplitude $R=R_{s}$ and phase $\psi=\omega T+\psi_{0}$, with

[TABLE]

appears in the supercritical ($\chi>0$) region for $c>0$, and in the subcritical region for $c<0$.

Remarkably, the coefficient $c$ is always positive and hence the Hopf bifurcation is always supercritical. This can be seen as follows: First, note that the denominator of $c$ remains always positive along the Hopf boundary Eq. (C6), and becomes zero at the Takens-Bogdanov point Eq. (C8). Second, note that the numerator of $c$ is positive for $\tilde{J}=0$ and may potentially change sign at $\tilde{J}_{c}=8/\tilde{g}$. However, this change of sign always occurs after the TB point ($\tilde{J}_{c}>\tilde{J}_{\text{{\scriptsize{TB}}}}$) in which the Hopf bifurcation ends, see Eq. (C8).

Finally, the approximate solution in terms of the original variables reads

[TABLE]

which describes an oscillatory motion in the critical eigenplane, with a small amplitude firing rate (for $\tilde{\eta}\geq\tilde{\eta}_{H}$)

[TABLE]