Bottom-up approach to torus bifurcation in neuron models

Deterministic mathematical modeling of oscillatory cells has been originally proposed within a framework of slow-fast dynamical systems written in a generic form,

where $x∈Rn$⁠, $n≥2$⁠, and $y∈R1$ represent fast and slow variables, respectively, $|μ|≪1$⁠, and $F,G$ are some smooth functions with $α$ being a parameter vector. For the sake of simplicity and convenience, the functions in mathematical models of cells are often chosen so that $G(x,y)=[f(x)−y]$ to have a Jordan block in the linearization matrix, with $G$ being bi-linear in its arguments. A geometric approach to the study of three-dimensional (3D) models of bursting neurons was proposed and developed by Rinzel.¹ The essence of the geometric understanding of neuronal dynamics is in the topology of the so-called slow-motion or critical manifolds in the phase space of the corresponding slow-fast model. The slow-fast dissection in the singular limit, $μ=0$⁠, freezes the slowest variable, $y$⁠, which is then treated as the control parameter to detect and parametrically continue critical, slow-motion manifolds: one-dimensional (1D) (given by $F=0$⁠) and two-dimensional ones made, respectively, of equilibria (eq.) and periodic orbits (PO) of the fast subsystem.

These manifolds, called quiescent M$eq$ and tonic-spiking M$PO$ in the neuroscience context, as some backbones, shape and determine the kinds of bursting activity in typical models of individual neurons both formal and derived through the Hodgkin-Huxley formalism. Moreover, one can loosely describe the particular bursting using the bifurcations that initiate and terminate these quiescent M$eq$ and tonic-spiking M$PO$ manifolds.² For example, “square-wave” bursting shown in Fig. FIG. 1,FIG. 2,FIG. 3,FIG. 4,FIG. 5,FIG. 6,FIG. 7,FIG. 8,9 can be alternatively called “fold-homoclinic,” while “square-wave” bursting depicted in Fig. FIG. 10,FIG. 11,FIG. 12,13 is due to two saddle-node bifurcations: one for equilibria and the other for periodic orbits that occur on the folds of the corresponding manifolds, M$eq$ and M$PO$⁠.

The dissection method, allowing a clear geometric description and understanding of the dynamics of low- and even high-dimensional models, is also helpful to predict possible bifurcation transitions between activity types, e.g., tonic-spiking and bursting. While the slow-fast dissection paradigm works well for most models with distinct timescales, it provides less or even misleading insights to understand transition mechanisms between activity patterns, which are due to nonlinear reciprocal interactions of slow and fast dynamics in both subsystems even for small but finite $μ$-values, which is typical for models of individual neurons discussed below. Such reciprocal interactions become not only non-negligible but pivotal for a comprehensive understanding of both formal and biologically plausible models in question when dynamics of the fast subsystem slows down to the time scale of the slow dynamics. This occurs universally when the fast subsystem undergoes an Andronov-Hopf (AH) bifurcation so that the rate of convergence to a critical equilibrium state is no longer exponentially fast but logarithmically slow, or when its equilibria undergo a saddle-node or fold bifurcation with a long dwelling time due to a bottleneck effect near the ghost of vanished equilibria. More complex yet very interesting cases include non-local bifurcations through which the fast subsystem possesses a stable orbit of long period that is about to become a homoclinic orbit to a saddle or to a saddle-node (SNIC), or when a pair of periodic orbits, stable and unstable, merges and vanishes in the phase space through the saddle-node bifurcation. One of the well-known examples of such slow reciprocal interactions is related to the phenomenon of canard oscillations of small amplitudes in a relaxation oscillator model (also known as the FitzHugh-Nagumo model in computational neuroscience³). The canard oscillations emerge through an AH bifurcation of an equilibrium state with characteristic exponents close to zero in the norm (⁠$ω∼μ$⁠). The shape of these small oscillations is defined by the fold on the characteristic $S$-shaped branch (called the fast nullcline and given by $F=0$⁠) of equilibria of the fast subsystem [see Figs. FIG. 14^,FIG. 15^,FIG. 16^,FIG. 17^,FIG. 18^,FIG. 19^,20(a) and 20(b)]. This equilibrium state of the full system is located at the transverse intersection of the fast nullcline with the slow (⁠$|μ|≪1$⁠) nullcline given by $G=0$⁠. [Note that if the intersection of the slow and fast nullclines is not transverse, i.e., when $∇F||∇G$⁠, then the equilibrium states in the full system bifurcate through a saddle-node]. The stability of the canard limit cycle emerging from the persisting equilibrium state through either a sub-critical or super-critical AH bifurcation is determined solely by the function in the right-hand side of the fast subsystem in the singular limit, $μ=0$⁠. Specifically, the type of AH bifurcation is determined by the sign of the first Lyapunov coefficient $l1$⁠,⁴ which is related to the smoothness of the function $f(x)$⁠. More specifically, given that $f′=0$ at the fold, and the constant concavity of the graph of $f(x)$ at the equilibrium state in question, then the sign of $l1$ is related to the sign of the coefficient of the cubic term in the Taylor expansion of $f(x)$ at the bifurcating equilibrium state. For example, depending on the choice of the function, say $x2$⁠, $x3$⁠, $1/(x2+1)$⁠, etc., whose all graphs have such a fold locally, the AH bifurcation can be either subcritical, $l1>0$⁠, or supercritical if $l1<0$⁠. As an illustration, the AH-bifurcation in the fast subsystem of Eq. (6) is subcritical with the term $v3/3$ but will become supercritical if it is replaced with $v3$ in the first equation. The same holds true when one considers the difference equations with discrete time variable, $t∈Z+$⁠, rather than differential ones, written in a form similar to Eq. (1) to analyze the local stability of invariant circle-canards emerging from a stable fixed point in the fold of the crucial manifold.^5,6 It happens often too that the sign of $l1$ derived through the fast subsystem in the singular limit $μ=0$ changes to the opposite in the full system even at small but finite $μ$ values. This indicates that the first Lyapunov value $l1$ is near zero and that the system is close to the codimension-2 bifurcation, occurring at the so-called Bautin point given by $l1=0$⁠, in the parameter space. This implies also that the slow-fast dissection in the singular limit should be only considered as a first-order approximation for comprehensive understanding of the dynamics of the full neuronal model with less desperate time scales.

The feature of the Bautin bifurcation is that its unfolding includes a saddle-node bifurcation curve corresponding to a double periodic orbit that either disappears or de-couples into stable and unstable orbits (repelling when $x∈R2$ and saddle in higher dimensions). As a result, the tonic-spiking manifold M$PO$ in the phase space of the system has a distinct fold (like the ones shown in Figs. 9, 12–15, f21,22, and 23) at a merger of stable, $MPOS$⁠, and unstable, $MPOU$⁠, branches of the tonic spiking manifold, M$PO$⁠. In slow-fast systems with a single slow variable, the existence of such a fold is a prerequisite, but not a guarantee, for the torus bifurcation at the transition between tonic-spiking and bursting activity.

A canard-torus emerges from a periodic orbit with a pair of multipliers $e±iϕ$ (⁠$ϕ∼μ$⁠) on a unit circle. However, the question about its stability is more challenging compared to that of a limit cycle canard emerging from an equilibrium state with a pair of complex conjugate exponents $±iω$ (⁠$ω∼μ$⁠), because the torus stability cannot be determined a priori by a local stability analysis. En route to bursting the stable periodic orbit, associated with tonic spiking activity, may lose its stability via a period doubling bifurcation (a single multiplier $eiπ=−1$⁠) initiating a cascade of such consecutive bifurcations or through a torus bifurcation which is unknown a priori whether it will be super- or subcritical.^7,8

In the case of the torus bifurcation, the emerging stable or unstable (repelling in 3D or saddle in 4D+ phase spaces) torus can be resonant or ergodic, depending on the angle $ϕ$ of its multipliers on a unit circle. For example, strong resonances such as 1:4 and 1:3 occur when $ϕ=π/2$ and $ϕ=2π/3$⁠, respectively.^4,9 The torus is born ergodic when $ϕ$ is not commensurable with $π$⁠. The example of a resonant 3:1 torus with a rational Poincaré winding number 3/1 is shown in Fig. 1. In simple terms, this means that there is a pair of periodic orbits, stable and unstable, on this 2D torus that correspond to a pair of period-3 points on the stable closed invariant curve, or invariant circle (IC) for short, on a 2D local cross section transverse to the torus. The winding number remains rational and constant within a synchronization zone (see Fig. 2), also known as an Arnold’s tongue, and takes varying irrational values outside of the synchronization zone. In this example of the resonant torus, the frequency locking ratio 3/1 means that the “horizontal” oscillations are three times faster than the “vertical” slow-components due to the timescales of $x$ and $y$-variables of Eq. (1).

Andronov and Vitt¹⁰ were first to study synchronization or phase-locking in the periodically forced van der Pol equation.¹¹ They showed that each wedge-shaped resonant zone is bounded by the homoclinic saddle-node bifurcation (SNIC) curves in the parameter plane so that having crossed one, the resonant torus becomes ergodic with everywhere dense covering (trajectory) on it. Later, Afraimovich and Shilnikov^12–15 proposed a phenomenological scenario of torus breakdown. It is illustrated in Fig. 2, which depicts pivotal evolution stages of the invariant circle of a 2D Poincaré return map and fixed points on it inside and outside of the resonant zone in the parameter plane. Here, $μ2=0$ corresponds to the torus bifurcation, and $μ1$ controls the angle of the multipliers, $e±iϕ$ of the bifurcating periodic orbit within the range from zero through $π$⁠. Outside the resonant zone, the invariant circle (IC) (read the 2D torus) is initially smooth and rounded at small $μ2$ values. It possesses a saddle-node fixed point (FP) on the SN-curves such that the closure of its 1D unstable set constitutes the IC. Inside the resonant zone, there are two FPs, stable and saddle, on the IC. As $μ2$ is increased, the IC becomes “non-smooth” and the stable FP undergoes a period-doubling bifurcation (on the PD curves) that breaks the torus down. Next, the unstable sets of the saddle FPs touch (on the HT-curves) start crossing their stable sets thereby creating homoclinic tangles which give rise to the onset of complex chaotic dynamics in the system. The reader is welcome to consult with the review¹¹ and the reference texts⁴ for more details about the torus breakdown.

In-depth understanding of the universal mechanisms of transitions between oscillatory activity patterns in single neuron models is a challenge for the theory of bifurcations. It requires the use of wide range of advanced apparatus of the bifurcation theory, ranging from the blue sky catastrophe to various homoclinic bifurcations of saddle orbits and equilibria.^{4,7,8,16–26}

In what follows, we will present several examples of biologically plausible and formal mathematical slow-fast models to demonstrate the universality and complexity of the torus bifurcations near the distinct fold on the tonic-spiking manifold. We use the bottom-up approach, starting with highly detailed models and ending with simple toy models, all featuring similar and generic properties. These bifurcations can result in the onset of stable, repelling, and saddle tori as well as their coexistence in the phase space leading to bi-stability of tonic spiking and bursting activities, as well as the onset of complex dynamics in nonlinear systems. We will also discuss the mechanisms of torus breakdowns through various resonances and conditions for period-doubling bifurcations that may alternatively occur near the same fold. In our first example, we perform an original analysis of period-doubling and torus bifurcations in a 12D hair cell model of the Hodgkin-Huxley type (HH-type). Using various computational toolkits, we demonstrate how stable tori break down after they become resonant. In addition, we show how the stability of the emergent torus can be changed by scaling down the changing rate of the single slowest variable of the model that gives rise to the bistability of coexisting tonic-spiking and bursting activities at the transition. The second occurrence of the stable torus bifurcation takes place in another highly detailed, 14D HH-type model describing complex dynamics in the pyramidal cells.²⁷ We argue that one can predict the torus-bifurcation once the topology of the slow-motion manifolds is revealed in the phase space by applying the parameter continuation technique to the full model in question. The third example is the 5D HH-type model of the Purkinje cells that features the saddle (not stable) torus that makes the co-existence of stable tonic-spiking and bursting activities possible in the model. We also show that the model can exhibit a stable torus bifurcation as was originally reported in Ref. 28. The analysis of the 4D HH-type model of parabolic bursters²⁹ reveals that the saddle-torus canard may unexpectedly occur on the emergent fold near the homoclinic saddle-node bifurcation (SNIC) of periodic orbits in the dissected fast-subsystem. Next, we consider a phenomenological slow-fast system that can exhibit a saddle-node bifurcation of a pair of stable and repelling tori coexisting in its 3D phase space. Using this model, we demonstrate how the torus-canard can evolve through a series of strong resonances up to the final 1:2-resonance that is further followed by the period-doubling bifurcation similarly to the hair cell model. The paper is concluded by the consideration of several 3D slow-fast normal forms derived to exhibit various torus bifurcations, stable and repelling, on a route from tonic-spiking to bursting. Our novel findings and computational approaches should be helpful for generalization onto other systems, especially in the context of life sciences, as well as should foster better understanding of singular torus-canards occurring at transitions between tonic-spiking and bursting in various slow-fast models of individual neurons in mathematical neuroscience, see Refs. 22,28,30–35, and references therein.

Our first example is a detailed Hodgkin-Huxley type model, whose 12 ODEs describe spontaneous dynamics of the membrane potential of a hair cell. Hair cells are peripheral receptors which transform mechanical stimuli into electrical signals in the senses of hearing and balance of vertebrates. These sensors rely on nonlinear active processes to achieve astounding sensitivity, selectivity, and dynamical range.³⁶ In particular, it was proposed that ciliary bundles of hair cells may self-tune to operate on the verge of (degenerate) AH bifurcations causing a sharpened selectivity to weak periodic mechanical forcing.^37,38 In amphibians, some hair cells demonstrate various spontaneous oscillatory activities. In particular, several experimental studies have documented various oscillation regimes of the membrane potential in the frog saccular hair cells.^39,40 In humans, the sacculus is a part of the peripheral vestibular system, encoding linear acceleration and head tilts in the vertical plane. The frog sacculus detects low-frequency seismic vibrations and is a part of the auditory system.^41,42 Although physiological implications of the membrane potential oscillations on sensory function of hair cells are not yet clear, experimental study⁴⁰ showed that spontaneous voltage oscillations drive sensory neurons resulting in periodic firing. Furthermore, the membrane potential of the hair cell affects the dynamics of its apical hair bundle compartment, as documented experimentally⁴³ and in modeling work.^44,45 Modeling work^44,46,47 predicted that oscillatory dynamics of the membrane potential may lead to the improved sensory performance of the hair cell.

The model is based on several experimental studies of basolateral ionic currents in saccular hair cells in bullfrog^39,40,48,49 and is a modification of a model developed in Ref. 49. It includes 6 ion currents plus a leak current with corresponding equations for the kinetics of ion channels and for the dynamics of Ca$2+$ concentration, resulting in a Hodgkin-Huxley type system of 12 coupled nonlinear ordinary differential equations. A detailed description of the model is available in an open access publication,³⁴ and a code for the right-hand sides of corresponding ODEs is provided in the supplementary material.

The membrane potential dynamics in the hair cell model is described by this equation:

The six ionic currents in Eq. (2) can be sub-grouped according to their activation patterns. The K$+$ inward-rectifier (⁠$IK1$⁠) and the cation h-type (⁠$Ih$⁠) currents are hyperpolarization activated. The rest four currents are depolarization activated: the voltage-activated K$+$ direct-rectifier (⁠$IDRK$⁠) and Ca$2+$ (⁠$ICa$⁠) currents; the calcium-activated K$+$ steady (⁠$IBKS$⁠) and transient (⁠$IBKT$⁠) currents. Finally, the leak current (⁠$IL$⁠) stands for the mechano-electrical transduction from the hair bundle compartment of the cell. The previous computational study³⁴ showed that a convenient choice of control parameters is the maximal conductance of K1-current, $gK1$⁠, and the strength of the Ca$2+$-activated potassium currents, $b$⁠.

Experimental work showed that depolarization-activated currents are responsible for the so-called phenomenon of electrical resonance, whereby the hair cell shows fast oscillatory responses when knocked by an external current pulse.^48,50–53 Furthermore, self-sustained voltage oscillations were observed experimentally and in a model which contained the voltage-activated Ca$2+$ and the Ca$2+$-activated potassium currents.⁵⁴ On the other hand, hyperpolarization-activated currents give rise to slow (3–4 Hz) large-amplitude oscillations of the membrane potential, owing to slow kinetics of $Ih$ current.⁴⁹ Taken together, the membrane potential of the frog hair cells shows diverse patterns of oscillatory activity documented experimentally⁴⁰ and well captured by the model.³⁴ In particular, the two distinct oscillatory mechanisms mentioned above may lead to quasi-periodic oscillations, with a torus bifurcation and torus breakdown,³⁴ which will be studied here in details.

To unravel the dynamical mechanisms behind the transition to quasiperiodic oscillations, we built a bifurcation diagram of the model for two control parameters, $gK1$ and $b$⁠. The diagram shown in Fig. 3 was constructed using the combined techniques of parameter sweeping and continuation.

In Fig. 3(a), the AH bifurcation (red curve) demarcates the region of the quiescence state from the oscillatory dynamics. The solid and dashed segments of the AH line correspond to super- and sub-critical AH bifurcation. The criticality of the bifurcation means that either a stable or a saddle periodic orbit emerges and collapses into the stable equilibrium state after the corresponding segment of the AH-curve is crossed. The point, labeled BP, is where the criticality of the AH bifurcation changes. This bifurcation was first studied by N. Bautin⁴ who found that its unfolding includes another curve corresponding to a double, saddle-node periodic orbit around the equilibrium state. The orbit results from a merger of stable and saddle ones, which earlier originate from a stable depolarized equilibrium state through super- and sub-critical AH bifurcations, respectively, upon crossing the solid and dashed branches of the (red) curve, AH, in the diagram in Fig. 3(a).

This saddle-node bifurcation is associated with the characteristic fold where two branches, foliated by stable and saddle periodic orbits, merge, on the tonic spiking manifold, $MPO$⁠, in the phase space of the model. The shape of the manifold can be revealed by the parameter continuation of the periodic orbit in the full system.⁸ The notion of the fold is imperative here, as it is where a tonic-spiking periodic orbit loses stability that gives rise to the onset of bursting activity through either a cascade of period doubling bifurcations or through a torus bifurcation. These scenarios, first described in Ref. 7, have happened to occur typically in various slow-fast systems with such folds,⁵⁵ particularly in neuronal models,^30,32 including the reduced model of Purkinje cells²⁸ below.

The stable periodic orbit born through the supercritical AH bifurcation loses its stability either through a torus bifurcation [TB, dashed black line in Fig. 3(a)] or through a period-doubling bifurcation [PD, solid green line in Fig. 3(a)]. A region bounded by these lines is characterized by various bursting patterns and by period-adding bifurcations. The colormap presented in Fig. 3(a) and showing the number of spikes per burst indicates that the number of spikes increases toward the torus bifurcation. To the right of the PD-curve, i.e., for $0.02<b<1$⁠, the model produces robustly tonic oscillations,³⁴ which underlines the stabilizing physiological role of the Ca$2+$-activated potassium currents (⁠$IBKS,BKT$⁠).⁴⁰ 

Next, we investigate the torus bifurcation and the torus breakdown using the slow-fast manifolds and Poincaré return maps on a route from tonic spiking to bursting. We also examine the transition to bursting on another route featuring a cascade of period-doubling bifurcations and contrast the corresponding Poincaré return maps in both cases. The Poincaré maps are well suited for exploration of complex oscillations to reveal the mechanisms of underlying bifurcations and their precursors. Since the voltage is the only observable variable in many experimental setups, we will employ the Poincaré return map $T$ defined on consecutive minima of the membrane voltage trace as follows: $T:Vmin(n)→Vmin(n+1)$⁠. In the case of a sparse map (insufficiency of various points $Vmin$⁠), e.g., for periodic or weak chaotic bursting, a small random perturbation can be added to the model, resulting in variability sufficient to obtain densely populated iterates, revealing some key dynamical features. We point out that the sensitivity of solutions to external perturbations is a common feature of slow-fast neuronal models.⁵⁶ 

The described two scenarios of the transition from tonic-spiking to bursting (through the torus or period-doubling bifurcations) can be differentiated with the aid of the Poincaré maps for consecutive minimal values of the membrane voltage traces shown in Fig. 4.

Let us first examine the torus bifurcation scenario, as $gK1$ increases with a fixed value of $b=0.01$⁠. Typical evolution of the voltage traces is shown in Fig. 3(b). As the parameter $gK1$ increases, small-amplitude periodic oscillations become quasi-periodic (slow amplitude modulation) through the supercritical (stable) torus bifurcation. When the torus breaks down, it gives rise to large-amplitude bursting oscillations which may be chaotic or regular in alternation. This is illustrated in the bifurcation diagram of Fig. 4(a), which was built using the Poincaré return map of consecutive voltage minima. In this diagram, a single point for a given parameter value corresponds to a stable periodic orbit. Several points per a single parameter value may correspond to different regimes. Dense points above $−70$ mV correspond to ergodic tori, whereas multiple extended branches correspond to stable periodic orbit(s) on the resonant torus exiting within its resonant zone. Points below $−70$ mV refer to bursting orbits where the cell becomes strongly hyperpolarized. Bursting orbits are stable when their several branches in the bifurcation diagram are continued without interruption within a stability window, followed by regions of chaotic bursting.

One can see from the bifurcation diagram in Fig. 4(a) that the stable periodic orbit, earlier emerging through the supercritical AH bifurcation, loses its stability when $gK1$ is increased through $29.19$ nS. This stability loss gives rise to the emergence of a 2D stable torus, first appearing like ergodic or very weakly resonant, as a widening parabola-shaped solid region as the diagram suggests. In what follows, the breakdown of the torus comes along with the generic scenario lines proposed by Afraimovich and Shilnikov.⁵⁷ The torus becomes resonant as $gK1$ increases through a synchronization boundary (of an Arnold’s tongue).^11,58 This assertion is supported by the occurrence of several branches corresponding to the local minima of a stable periodic orbit coexisting with a saddle orbit on the torus after a saddle-node bifurcation at $gK1≃29.199$ nS. Further increase of the control parameter causes the torus breakdown, which leads to bursting oscillations with modulatory (quasi-periodic) tonic spiking episodes interrupted by slow hyperpolarized quiescent periods [Fig. 3(b) at $gK1=29.213$ nS]. Figure 4(a) also reveals a number of spike adding bifurcations of bursting. The cascade of spike adding ends on the upper branch of the green bifurcation curve, which corresponds to a reverse period-doubling bifurcation [Fig. 3(a)]. To the right of the curve, the hair cell model produces stably large-amplitude voltage oscillations [Fig. 3(b) at $gK1=40$ nS]. The black dashed curve in Fig. 3(a) corresponds to the torus formation and therefore quasi-periodic tonic spiking oscillations in the model.

Figures 5 and 6 represent the Poincaré return maps depicting the characteristic stages of the torus formation and breakdown. In Fig. 6(a), a family of invariant circles (ICs) corresponds to the growing of the first smooth and later non-smooth ergodic stable IC emerging from a tonic-spiking periodic orbit, as the latter has lost “its skin” to the 2D torus [Figs. 5(a) and 5(b)]. As $gK1$ increases, the torus becomes resonant with a stable periodic orbit on it. The degree of the resonance can be evaluated by the number (here 7) of periodic points of the stable orbit, which equals the number of the “sleeves” originating from an unstable FP toward the non-smooth invariant curve. This agrees well with the theory of torus breakdown that causes the onset of complex chaotic dynamics. In essence, the torus breakdown begins with the emergence of two periodic orbits, stable and saddle, on the torus when it becomes resonant through a saddle-node bifurcation. Later, the torus becomes non-smooth after the stable and unstable manifolds (sets) of the saddle orbit start making wiggles transitioning to homoclinic tangles with the stability loss of the counterpart periodic orbit through a period-doubling bifurcation. This bifurcation can be identified by end-branching on the resonant traces left by the periodic-point as the control parameter $gK1$ is increased. One can observe from Fig. 5(d) that the stable period-7 orbit has a leading pair of complex-conjugate multipliers, which indicates the proximity of forthcoming period-doubling bifurcations (as documented by the green branches indicated in Fig. 7).

We note that regular bursting can already co-exist with the stable torus in the hair cell model. Bursting becomes chaotic with the torus breakdown starting at $gK1=29.213$ nS. The corresponding 1D Poincare return map $T:Vmin(n)→Vmin(n+1)$ is shown in Fig. 6(b). It has a distinct shape that is characteristic for the square wave bursters.^25,56 The map has three components: the transient “toroidal” section revealing a ghost of the disappeared non-smooth IC around $Vmin=−65$ mV, a dropping-down middle section, and a flat (contracting) section corresponding to the slow hyperpolarized quiescence phase of bursting reaching $Vmin=−85$ mV. The voltage trace generating this map at $gK1$ is depicted in Fig. 3(b) at $gK1=29.213$ nS; it demonstrates slowly modulated chaotic tonic-spiking phases alternating with hyperpolarized voltage sags.

The resonance of a torus is determined by its winding number, calculated as the ratio of 2$π$ and the angle $ϕ$ of the Floquet multipliers $e±iϕ$ on a unit circle.⁴ For example, when the angle is $π/2$ [Fig. 7^,$(a2)$], the winding number is 4, meaning that the resonance of the torus is 1:4. The irrational and rational winding numbers correspond to ergodic (quasi-periodic) and resonant (periodic) torus, respectively. The cases 1:1, 1:2, 1:3, and 1:4 are considered as strong resonances and others as weak.^9,59 Figure 7 shows the emergence, evolution, and breakdown of the 2D tori through three strong resonances, namely, 1:2, 1:3, 1:4, and one weak case. When the torus becomes resonant, the corresponding IC (gray color) of the Poincaré return map possesses a stable orbit [made of a specific number of periodic points, like the ones shown as the blue branches in Figs. 7,$(a1)$ and 7,$(b1)$]. The resonant ICs become non-smooth and then starts breaking down when the leading multipliers of the stable periodic orbit become complex conjugate toward forthcoming PD bifurcations (shown as green and orange dots in the return maps in Fig. 7) with the increase of the parameter $gK1$ at the indicated values of parameter $b$⁠, so that the unstable sets of the saddle periodic orbits no longer converge monotonically but start spiraling onto the stable periodic points of the IC [see Fig. 5(d)].

Next, we consider the pathway at $b=0.015$ passing through the period-doubling curve (green line) in the 2D parameter plane shown in Fig. 3(a). Having crossed this curve, the stable periodic orbit (representing tonic-spiking oscillations) loses its stability to the one with the doubled period. This bifurcation and the new periodic orbit correspond to the branching in the diagram shown in Fig. 4(b). A cascade of forthcoming period-doubling bifurcations develops rapidly (this is a typical, like canard, phenomenon of slow-fast systems) as $gK1$ increases and terminates with the onset of regular bursting in the hair cell model. Spike adding causes the onset of chaotic bursting in the transition between stable 5-spike and 4-spike bursting as this bifurcation diagram reveals.

The Poincaré return map of the chaotic bursting shown in Fig. 8(a) has a wider drop-shaped section, compared to chaotic bursting emerging through the torus scenario [Fig. 6(b)]. Quantitatively, the width of the drop indicates how fast the model switches back and forth between tonic spiking and hyperpolarized quiescence in bursting. The wider the drop in the map is, the less rigid the switching is, or, alternatively, the less contrast between fast and slow dynamics is. Since the switches are quite loose, not constrained, there is a wide variability in transitions that gives rise to multiple branches of the return map depicted in Fig. 8(a). The multiple branches look like as they represent snapshots of some return maps at several parameter values in progression. Each branch reveals a generic unimodal (non-sharp) shape of the return map that is needed for a cascade of period-doubling bifurcations to occur. We note that other Hodgkin-Huxley type models with the fold on the tonic spiking manifold can exhibit a period-doubling cascade leading to a chaotic tonic-spiking attractor in the phase space.³⁵ The progression of the branches also suggests that the stable periodic orbit emerges through a saddle-node bifurcation when the graph of the map touches the $45°$ line. Recall that the corresponding bifurcation comes out of the codimension-two Bautin point in the parameter plane in Fig. 3(a). Next, as the map graph lowers further down and becomes steeper, the stable fixed point becomes unstable, thus giving rise to a period-two orbit; its points are indicated by two green solid circles in Fig. 8(b). The bottom four branches of the map correspond to bursting where each hyperpolarized voltage drop is followed by chaotic tonic-spiking phases developed through all steps of the period doubling cascade.

Precursors of bursting can be revealed by adding small random perturbations to the model. Gaussian white noise $ξ(t)$ is injected into the right-hand side of Eq. (2) with the term $σξ(t)$⁠, where $σ$ is the standard deviation of random perturbation. Figure 8(b) demonstrates that weak random perturbations can effectively reveal precursors of transition to bursting. For the value of $gK1$⁠, the model possesses a stable period-2 orbit, which contains only two points in the Poincaré return map. Weak noise induces bursting, resulting in the map which has a similar shape as the deterministic chaotic bursting.

We stop here to presume that without 1D Poincaré return maps generated through voltage oscillations, it would be virtually impossible to characterize the exact formation mechanisms underlying transitions from simple tonic spiking to bursting and back in this model as well as in other Hodgkin-Huxley type neuronal models.^25,56

Bifurcations underlying the transitions between different oscillatory patterns in the hair cell model can be explained geometrically using the slow-fast dissection as illustrated in Fig. 9. The essence of the approach is based on dissecting the original dynamics of the model into slow and fast components. In the model, the activation of the $h$-current, $mh$⁠, is the slowest dynamical variable, while all other variables are fast compared to it. In particular, we choose the membrane potential, $V$⁠, and the activation of the K$+$ DRK-current, $mDRK$⁠, as representatives of the fast subsystem. The feature of a slow-fast system that often reduces the complexity of its dynamics is that the phase trajectory stays close to the slow-motion manifolds of low dimensions: the 2D tonic-spiking manifold, $MPO$⁠, foliated by small-amplitude periodic orbits corresponding to tonic spiking activity, and the 1D $S$-shaped quiescent manifold, $Meq$⁠, is comprised of hyper- and depolarized equilibrium states of the model. This approach implies that bursting, being a multiple-time-scale regime, is interpreted as repetitive fast switching between these manifolds at their terminal phases, followed by slow passages of trajectories turning around $MPO$ and then shifting along $Meq$⁠, corresponding to alternating episodes of tonic-spiking and slow hyperpolarized quiescence, respectively. Figure 9 depicts the topology of these slow-motion manifolds for the parameters close to the torus (dark blue trajectory) bifurcation, as well as shows several bursting orbits superimposed (green, light blue, and brown trajectories) in the phase-space projection. These spiking and quiescent manifolds are obtained by the parameter continuation of periodic orbits and equilibria of the model with the auxiliary parameter shifting the slow nullcline $mh′=0$ up and down.

In what follows, we outline how the manifolds are detected in the 12D phase space of the hair cell model. We use a new and practical way for localization of both slow motion manifolds, 1D quiescent, $Meq$⁠, and 2D tonic-spiking, $MPO$⁠, composed of equilibria and periodic orbits of the full system. The method capitalizes on the slow-fast dissection as well as on the parameter continuation technique.⁶⁰ The method supplements the conventional slow-fast dissection in the singular limit with the slow variable frozen as a parameter. The main essence is the parameter continuation technique applied to the entire set of the model equations. The advantage of our approach is that it yields the sought manifolds themselves in the phase space of the intact model rather than the manifolds of its fast subsystem without taking into account slow, yet finite component of the overall dynamics. This approach is especially valuable for complex models of higher dimensions where slow-fast dissections could be problematic because of multiple time-scales of various ionic currents involved in complex dynamics especially at transitions. Furthermore, this approach can detect and explain bifurcations, often non-local due to reciprocal interactions of slow and fast dynamical variables that underlie these transitions in the model in question.

Specifically, we introduce a faux control parameter, which is a voltage offset that affects only the slowest dynamics of the model; this offset shifts the voltage threshold at which the slow gating variable/channel is half-open. Geometrically, its variations lift or lower the slow nullcline, $mh′=0$⁠, in the 3D phase-space projection depicted in Fig. 9. As the voltage offset changes, the intersection point of the slow nullcline (sigmoidal surface) with $Meq$ (the equilibrium state of the model) is shifted along $Meq$⁠, thus tracing $Meq$ down and explicitly revealing its position and shape in the phase space. The same is true for the other manifold $MPO$⁠. More details on this approach can be found in Ref. 35.

Like many neuron models, the quiescent manifold in the hair cell system has a characteristic $S$-shape in the projection onto the slow-fast (⁠$mh,V$⁠)-plane. The shape of $Meq$ implicitly endows the model with a hysteresis required for dynamic switching between hyper- and depolarized phases in large-amplitude bursting. Reaching the low, hyperpolarized fold on $Meq$ indicates the beginning of a burst. In addition, the two real bifurcation parameters, $gK1$ and $b$⁠, of the model can shift the position of the manifold $Meq$ relative to that of the slow nullcline, $mh′=0$ in the phase space, so that the intersection point can move onto the stable upper or low section of $Meq$⁠. In those cases, the hair cell rests on the depolarized or hyperpolarized steady state, respectively.

The hair cell model oscillates tonically when there is a stable periodic orbit on $MPO$⁠. This orbit emerges from the depolarized equilibrium state through the supercritical AH bifurcation (discussed above). Moving the slow nullcline by varying the faux parameter makes the periodic orbit slide along $MPO$⁠, thus revealing its shape. The manifold terminates through a homoclinic bifurcation occurring after the periodic orbit of long period touches the middle, saddle branch of $Meq$ to become the homoclinic orbit of the saddle. This is a generic configuration for square-wave bursters.

In the case where the stable branches of the tonic spiking $MPO$ and quiescent $Meq$ manifold are not cut through by the slow nullcline, $mh′=0$⁠, i.e., both manifolds are freely transient for trajectories passing by, and the hair cell model exhibits bursting. Bursting is represented by solutions repeatedly switching between the low, hyperpolarized branch of $Meq$ and the tonic spiking manifold $MPO$⁠. The number of complete revolutions of the solution around $MPO$ is that of spikes per burst in the voltage traces. This number is used to classify the bursting activity. The larger the number is, the longer the burst duration lasts. The relative position of the torus to the fold on $MPO$ can be changed by varying the rate of the slow variable $mh$⁠. The slower the rate is, the closer the stable torus emerges to the fold as shown in Figs. 10(a) and 10(b), as well as the closer the torus bifurcation (TB) is detected near the fold by the bifurcation package MATCONT.^60,61 More interestingly, when the rate of $mh$ is about two times slower than the original one, the torus bifurcation becomes sub-critical and generates a saddle torus at the fold instead [see Fig. 10(c)]. This gives rise to bistability in the system, where both tonic spiking and bursting can be produced by the model depending on the initial conditions, as illustrated in Fig. 11. A small basin of attraction of the stable tonic-spiking orbit is “bounded” by the 2D saddle torus in the 12D phase space from that of the stable bursting orbit. We note that similar bistability is also found in the Purkinje cell model, parabolic burster model, and FitzHugh-Nagumo-Rinzel (FNR) model (see Secs. IV–VI).

In a system with drastically distinct slow and fast time scales, the tonic-spiking periodic orbit would lose its stability to a 2D torus or through a period-doubling bifurcation cascade right at the location of the fold.³² The type of the bifurcation depends on nonlocal properties of the flow near the fold. Specifically, depending on whether the phase space volume contracts or not, the tonic spiking orbit undergoes period-doubling or torus bifurcation in the system under consideration. While the verification of this condition seems doable in 3D slow-fast systems (see the FNR section below), the high-order models may present a challenge as one has to identify the set of the variables (from a 3D subspace containing the central manifold in which the bifurcation takes place) that are involved directly in the ongoing bifurcation.

Our next example is another highly detailed 14D Hodgkin-Huxley type two-compartmental (for soma and dendrite) model describing interactions of pyramidal cells and fast-spiking inhibitory interneurons.⁶² The reader can find details of all currents constituting the model in Ref. 27, which describes how electrogenic properties of the $Na+/K+$ ATPase control transitions between normal and pathological brain states, specifically during epileptic seizures. The equations describing the dendritic and somatic voltage dynamics are

where subscripts s and d stand for soma and dendrite, respectively, $Cm$ is the capacitance of the dendritic compartment, $V$ is the membrane potential, $Is$ are various ionic currents including $Na+$/$K+$ pump current, variable s stands for the surface area of the corresponding compartment, $gc$ is the coupling conductance of the soma and the dendrite, and $α$ is a scaling coefficient. The parameters and a code of the right-hand sides of corresponding ODEs are provided in the supplementary material.

Unlike in the hair cell model, where the dynamics of intracellular calcium concentration, $[Ca2+]i$⁠, is relatively fast, in the given pyramidal cell model $[Ca2+]i$ is the slowest dynamical variable. As before, we will use the 3D phase projection to disclose the topology of the slow-motion manifolds that determine the shape of bursting oscillations in the model. Figure 12 depicts a similarly shaped 1D quiescent manifold $Meq$⁠. Its lower, hyper-polarizing knee corresponds to the homoclinic saddle-node bifurcation giving rise to the stable (green) outer section of the 2D tonic-spiking manifold $MPO$⁠. This manifold wraps back inwardly so that its unstable (pink) shrinking section terminates on the depolarizing branch of $Meq$ through a sub-critical AH bifurcation. The surface $[Ca2+]i′=0$ (yellow) is the slow nullcline: $[Ca2+]i$ increases/decreases above/below it in the phase space. When this nullcline is slightly below its position depicted in Fig. 12, there would be a stable periodic orbit on the stable branch of $MPO$ that corresponds to tonic-spiking activity. As the nullcline $[Ca2+]i′=0$ is shifted up, the stable orbit moves closer to the fold where it loses its stability through a supercritical torus bifurcation. The phase point that moves along the orbit on the surface of the stable newborn torus oscillates back and forth between the stable and unstable branches of $MPO$ near the fold. The torus is stable because it lingers longer on the stable branch than on the unstable one on average. Otherwise, it would be of the saddle type, i.e., repelling in the $mNa$ and $[Ca2+]i$ coordinates and stable in all the other dimensions. One can observe from Fig. 12(a) that gaps between the successive points filling in the stable invariant curve are large enough to assume that dynamics of $[Ca2+]i$ is not that slow compared to the torus-canard observed Purkinje cell model discussed next. After the torus breaks down as the slow nullcline is further shifted up, the pyramidal cell model demonstrates bursting activity; see Ref. 27.

Purkinje cells are the main output of the cerebellum and are involved in motor learning/conditioning.⁶³ To study these cells, several highly detailed models were constructed based on experimental data.^64,65 A reduced yet quantitatively-similar 5D Hodgkin-Huxley type model was proposed and studied in Ref. 28. Its generic representation reads as follows

where the vector $x$ represents the gating variables: fast n, h, c and slow M, for the following four ionic currents: a delayed rectifier potassium current, a transient inactivating sodium current, a high-threshold non-inactivating calcium current, and a muscarinic receptor suppressed potassium current; here, the external current $Iapp$ is the bifurcation parameter of the model. The reader is welcome to consult with the original paper²⁸ for a more detailed description of the model and the supplementary material for its MATLAB code.

In this 5D model,the gating variable $M$ is the slowest one. A faux parameter, $Δ$⁠, is introduced in the right-hand side of the equation [see Eq. (4)] governing $M′$ to perform the parameter continuation letting us sweep the slow-motion manifolds using the package MATCONT.⁶⁰ 

Figure 13(a) shows that the tonic-spiking manifold with the distinct fold consists of two components: the inner stable (green) cylinder-shaped M$POS$ and the outer unstable (red) part M$POU$⁠. Note that M$POU$ vanishes via the homoclinic bifurcation (HB) when it touches the middle, saddle branch of the quiescent manifold M$eq$ (⁠$S$-shaped black line); see also 2D phase-space projections in Figs. 14 and 15.

It was first reported in the original paper²⁸ that during the transition from tonic spiking to bursting, the Purkinje cell model shows slow amplitude-modulation of tonic-spiking activity in the voltage traces [Fig. 13(c)]. This modulation corresponds to a stable torus-canard that emerges around the fold of the tonic spiking manifold M$PO$ in the phase space of the model [see Figs. 13(a) and 13(b)]. This torus-canard exists in a rather narrow parameter window beyond which the model quickly transforms into a square-wave buster.

To examine the stability of the torus in this model, we tracked the stable state of the model by varying the applied current $Iapp$ in small incremental steps (⁠$∼1×10−4$⁠) in both directions—increasing and decreasing. Simulations for each applied current value were long enough for the model to reach some steady states. We then represented the stable states by minimal values, $Vmin$⁠, of their voltage values as shown in Fig. 14(a).

Interestingly, for $Iapp∈[−29.4871,−29.4762]$ [denoted by the white bar in Fig. 14(a)], the stable state in $Iapp$-increasing direction is tonic spiking [green dots in Fig. 14(a)], while in $Iapp$-decreasing direction is bursting (blue dots). That is, the system is bistable within a small overlapping window of $Iapp$⁠. What “separates” the two stable states is a saddle torus-canard detected at $Iapp=−29.487$ (red box), whose corresponding voltage trace shows slow amplitude-modulation [Fig. 13(c)]. The phase-space trajectories corresponding to the two stable states and the separating torus-canard are shown in Fig. 14(b), and the corresponding voltage traces are shown in Figs. 14(c) and 14(d). The vertical red line [Fig. 14(a)] indicates the I$app$ value where MATCONT detects torus bifurcation.

As for most slow-fast systems, the location of the tonic-spiking periodic orbit on the manifold in the phase space in the Purkinje cell model at various parameters can be accurately predicted by the averaging technique; see Ref. 35 for details. The approach is illustrated in Fig. 15(a). Its core is the use of the parameter continuation to detect periodic orbit(s) on M$PO$ where the average derivative of the slow variable vanishes as well as to find the corresponding function in the right-hand side of the average slow equation to predict forthcoming bifurcations like the blue sky catastrophe.³⁵ Specifically, for the Purkinje cell model, first, $⟨M′⟩$ and $⟨M⟩$ are calculated for each periodic orbit on the tonic-spiking manifold based on the continuation data; then $⟨M′⟩$ is plotted against $⟨M⟩$ (which represents each periodic orbit on the tonic spiking manifold) to pin down where the stable tonic spiking occurs [Fig. 15(⁠$a2$⁠)]. The prediction from the averaging technique showing where $⟨M′⟩=0$ [dark green dot in Fig. 15(⁠$a2$⁠)] is quite precise, as it matches very well the result of the simulations [dark green PO in Fig. 15(⁠$a1$⁠)].

It is obvious though that the averaging technique cannot explicitly establish the occurrence of the torus bifurcation in a system with a single slow equation. Note, however, that the averaging approach illustrated in Fig. 15(⁠$a2$⁠) is done at the parameter values (⁠$Iapp$ and $Δ$⁠) when the saddle-torus separates two stable oscillatory states: spiking and bursting. Although the averaging approach is only applicable to detect tonic-spiking state orbits, it can be yet helpful to locate the torus on the spiking manifold [Fig. 15(⁠$b1$⁠)]. Actually, the averaging applied to the torus orbit in the phase space does reveal the corresponding invariant circle in the $⟨M′⟩−⟨M⟩$ projection as shown in Fig. 15(⁠$b2$⁠).

The same kind of bistability due to a saddle torus at the fold, separating the stable tonic spiking and bursting, is also observed in a simpler 4D model of a parabolic burster. This Hodgkin-Huxley type model is given by

see its complete set of equations in the supplementary material. This model for an endogenously parabolic burster, i.e., with the inter-spike frequency being maximized in the middle of each burst, was proposed and examined in Ref. 29. Its key feature is the homoclinic saddle-node bifurcation (SNIC) in the fast 2D subsystem that should demarcate the entry and exit points of the corresponding tonic-spiking and quiescent phases of bursting activity in the full system. Figure 16 shows the topology of the slow-motion manifolds in the phase-space projection of the model. Here, the 2D cylinder-shaped manifold $MPO$ originates from the depolarized branch of the 1D quiescent manifold $Meq$ through a subcritical Andronov-Hopf bifurcation. Next, the unstable (purple) section $MPOU$ merges with the stable (green) section $MPOS$ through the saddle-node bifurcation at the first fold. The second fold on $MPO$ takes place at the higher values of the $[Ca2+]$-variable, where the torus bifurcation occurs in the full system. Here, the saddle torus (red) bounds the attraction basin of the tonic-spiking orbit (hidden inside the saddle torus) from that of the stable bursting orbit (black). Inset 16(a) depicts a local transversal cross section with a repelling IC, representing the saddle-torus, that surrounds a light blue fixed point corresponding to the stable tonic-spiking orbit next to the right fold of the spiking manifold $MPO$⁠.

The FitzHugh-Nagumo-Rinzel (FNR) system is a formal mathematical neuron model of an elliptic buster³² 

where $v$ and $w$ represent the fast “voltage” and “gating” variables, while $y$ is the slow variable due to the smallness of $μ=0.002$⁠. Here, we fixed $Iext=0.3125$ and use the parameter $c$ in the slow equation to continue and reveal the topology of the slow-motion manifolds of this system; the other bifurcation parameter $δ$ is used to demonstrate how time scales of the variables in the model can determine the bifurcations (torus and period-doubling) of periodic orbits of this model. Figure 17(a) shows the tonic-spiking manifold of the FNR model consisting of an unstable inner layer M$POU$ (purple surface) emerging through a subcritical Andronov-Hopf bifurcation from the 1D quiescent manifold, M$eq$⁠, and a stable outer layer M$POS$ (green surface), both merging at a fold.

The $(c,δ)$-bifurcation diagram [Fig. 18(a)] of the FNR model is obtained using the parameter continuation package MATCONT.⁶⁰ It includes several bifurcation curves: AH stands for the Andronov-Hopf bifurcation of the only equilibrium state of the FNR model; TB stands for the torus bifurcation; along this curve, the multipliers, $e±iϕ$⁠, of the bifurcating periodic orbit run through strong resonances detected at points labeled as R$∼$1:4 with $ϕ=π/4$⁠, R$∼$1:3 with $ϕ=2π/3$⁠, and the final resonance R$∼$1:2 where $ϕ$ completes the arch of the unit circle and reaches $π$⁠. The last point is also included in the PD-curve corresponding to the period doubling bifurcation through which the periodic orbit loses its stability. The diagram also includes two bifurcation curves, SN, of saddle-node periodic orbits. These curves follow the PD-curve closely. They cannot be further continued in the $c$-parameter any further down the bifurcation diagram for the given slope of the slow nullcline $v=c−y$ (being a 2D plane in the 3D phase space) because it becomes tangent to the 2D spiking manifold M$PO$⁠. We hypothesize that the bifurcation curves, SNs and PD, will terminate at the corresponding resonances, 1:1 and 1:2, occurring near the point labeled as a red dot R (with $δ≃0$⁠) in the $(c,δ)$-parameter plane. Geometrically, this point corresponds to a cusp underlying the transition from the hysteresis to the monotone, increasing dependence of the V$max/min$-values and $⟨V⟩$-coordinates of the periodic orbits on the $c$-parameter in the bifurcation diagrams in Figs. 18(b) and 18(c), respectively.

Depending on the level of the $δ$-parameter, the periodic orbit can undergo two supercritical torus bifurcations, forward and backward, when it loses stability to the torus and re-gains for $δ=0.3$ (along the upper dashed gray line, observe the $∩$-shape of the TB-curve) or go through the forward torus bifurcation followed by the period-doubling bifurcation as $c$ is increased for $δ=0.08$ (lower dashed gray line).

The bifurcation diagrams in Figs. 18(b) and 18(c), representing the maximum, $vmax$⁠, and minimum values, $vmin$⁠, of the voltage plotted against the $c$-parameter, illustrate how the tonic-spiking manifold is shaped and where the periodic orbits that constitute it, bifurcate (indicated by vertical red lines) as $c$-parameter is varied for fixed values of $δ=0.08$ and $δ=0.3$⁠, respectively. One can clearly see that at a rather large value $δ=0.3$⁠, the torus bifurcations do not occur close to the fold of M$PO$ unlike the period-doubling bifurcation that occurs between two saddle-node bifurcations [indicated by the vertical lines labeled as SNs in Fig. 18(b)] corresponding to the two folds on M$PO$ at $δ=0.08$⁠.

To verify the bifurcations detected by the parameter continuation in the FNR-model, we ran straight-forward $c$-parameter sweeping simulations of the limiting solutions of the model when $δ=0.08$ and when $δ=0.3$⁠. The corresponding sweeping diagrams are represented in Figs. 18(d) and 18(e).

At $δ=0.3$⁠, after crossing AH-line in the bifurcation diagram of Fig. 18(a), trajectories of the FNR-model are attracted into a whirlpool stirred by a stable PO in the phase space shown in Fig. 18(⁠$f1$⁠). After crossing the left branch of the TB-curve, this stable PO loses stability (red circle) to a stable (ergodic and round) torus in the phase space depicted in Fig. 18(⁠$f2$⁠). The sweeping bifurcation diagram in Fig. 18(e) shows that with a further increase of $c$⁠, the torus becomes resonant [Fig. 18(⁠$f3$⁠)] and its shape becomes more complex giving rise to the so-called mixed-model oscillations (MMO) of alternating small-amplitude sub-threshold and large spiking ones. Once crossing the right branch of the TB line, the torus goes into a strong resonance 1:3 before it breaks down at larger values of the $c$-parameter giving rise to the stable period-3 orbit in the phase space; see Fig. 18(⁠$f4$⁠).

At $δ=0.08$⁠, after crossing the left branch of the TB-curve in the bifurcation diagram in Fig. 18(a), the PO loses stability [red line in Fig. 17(c)] to a stable torus (blue line) in the phase space. Interestingly, at $c=−0.94415$⁠, this small stable torus bounding the unstable PO is nested inside a larger repelling unstable torus (gray line). The coexistence of both tori is better seen in a local Poincaré cross section (an orange plane) chosen transversally to these orbits as shown in Fig. 17(d). Note that in a 3D system, a repelling torus can be detected and traced down in the backward time when it becomes attracting, which is impossible for saddle tori occurring in high-dimensional systems.

Again, we note that the tori emerge in the FNR model near the smooth cone-shaped end of M$POU$ instead of around its fold [Fig. 17(a) inset]. Figure 17(e) shows the slowly-modulated “voltage” traces corresponding to the two tori. As $c$-parameter further increases to $c=−0.94$⁠, elliptic bursting becomes dominantly stable in the phase space shown in Figs. 17(a) and 17(b). The period-doubling cascade starts at the left branch of the PD-curve and ends at its right branch in the bifurcation diagram in Figs. 18(a) and 18(d) [magenta line in Figs. 17(a) and 17(b), $c=−0.6191$]. After that, tonic spiking (PO in phase space) becomes stable [dark green line in Figs. 17(a) and 17(b), $c=−0.5$]. Note that in the FNR model, PD instead of TB occurs at the fold.

The same approach employed to detect bistability in the Purkinje cell model due to a hysteresis effect is used to visualize bistability in the FNR-model as its solutions are swept in the opposite directions of the $c$-parameter. The obtained bifurcation diagram, with $vmin$-values (as green/blue dots) plotted against increasing/decreasing $c$-values, is represented in Fig. 19. Here, a single dot for the given $c$-value corresponds to a stable periodic orbit of a single $vmin$-value, while several vertical dots correspond to co-existing busting orbits. The basins of these orbits are bounded by the repelling torus enclosing the stable tonic spiking orbit in the 3D phase space of the FNR-model.

Let us outline a possible origin and stability of the torus bifurcation in the FNR-model. It starts with a saddle-node (fold) bifurcation of two limit cycles, stable and repelling, that occurs in and due to the 2D fast $(v,w)$ sub-system of the model, i.e., whenever the slow $y$-variable is frozen and used as a parameter. Specifically, the value $y=0.0117$ is the critical one at which the stable (blue circle) $LCS$ and unstable (red circle) $LCU$ limit cycles emerge through a saddle-node (fold) bifurcation in the $(v,w)$-subspace; see Fig. 20(a). While the stable limit cycle $LCS$ increases in its size to become the large-amplitude orbit of the relaxation oscillator, the unstable limit cycle $LCU$ keeps shrinking as $y$ increases; see the corresponding $(v,w)$-phase portrait in Figs. 20(b) and 20(c). As a result, in the full model, the slow increase/decrease of the $y$-variable makes the trajectories drag around $LCU$ back and forth to form the 2D repelling torus in the 3D phase space of the model; see Fig. 20(d).

Either the torus bifurcation (TB) or the period-doubling bifurcation (PD) can underlie the stability loss of the periodic orbit existing on the attracting cylinder-shape section of the tonic-spiking manifold M$PO$ in the phase space of the FNR-model when it reaches the fold; see Fig. 17. We have already noticed that depending on the value of the $δ$-parameter, this large amplitude orbit loses stability through either the torus bifurcation followed by its resonances and breakdown for $δ≥0.3$ or a period doubling bifurcation at smaller values like $δ=0.08$⁠. This indicates wherever the time scales of the $v$- and $w$-variables become of a similar order, the torus bifurcation case prevails over the period doubling one; and wherever both $w,y$-variables become much slower than the $v$-variable, the period doubling one prevails.

To predict what bifurcation, torus or period doubling, the stable orbit will next undergo, without computing its multipliers, we can evaluate the sign of the divergence given by $∇=∂v′(t)∂v+∂w′(t)∂w+∂y′(t)∂y$ along the periodic orbit on the fold of the slow-motion manifold in the 3D phase space of the FNR model. The arguments are similar to the well-known Bendixson-Dulac criterion for the existence of a closed orbit, like a limit cycle and a separatrix loop of a saddle in a system on a 2D plane; note, the Dulac’s argument follows from the Green’s theorem.⁶⁶ Loosely speaking, it can be directly re-iterated as follows: no closed orbit exists in a domain of the 2D phase plane where the vector field divergence is of a constant sign. This implies that no 2D torus cannot emerge in the 3D phase space of a dissipative system with the divergence of some constant sign. For example, a negative divergence of the vector field would collapse phase volumes. As such, for a torus bifurcation to occur in the system in question, divergence must be phase-coordinate dependent and hence sign-alternating. This is a partial explanation why the torus bifurcation is a part of bifurcation unfolding of an equilibrium state with the characteristic exponents $(0,±iω)$ (the sum of the exponents vanishes), which is also referred to as the codimension-two Gavrilov-Guckenheimer or fold-Hopf bifurcation; see Refs. 4,59,67, and 68. Such a bifurcation may also occur in the FNR-model near the tip of the cone-shaped manifold M$PO$ crossed out by the slow nullcline $y′=0$ with $μ≪1$⁠.

Thus, if the averaged divergence, $⟨∇⟩$⁠, i.e., all $∇$-values summed up along the periodic orbit and averaged over its period, on the fold of the manifold M$PO$⁠, is close to zero, then the stable periodic orbit may undergo the torus bifurcation at the stability loss. Otherwise, if $⟨∇⟩<0$⁠, which is typical for most 3D dissipative systems, the periodic orbit will lose stability through the period-doubling bifurcation. The calculation of $⟨∇⟩$ in the FNR model at various points on the PD- and TB-curves supports our hypothesis. Figure 21 shows how $∇$ varies over the normalized period of the orbit losing the stability through the period-doubling and torus bifurcations. Its $⟨∇⟩$-value is about $∼0.56$ and close to 0 in the period-doubling and torus cases, respectively, at $δ=0.08$⁠. This de facto validates the divergence-based approach and that it can be used to predict the type of the bifurcations that will occur at the loss of stability of periodic orbits on the distinct fold of the slow-motion manifold in the phase space of 3D slow-fast systems.

The amplitude modulation phenomenon that happens in various complex and realistic neuron models inspires us to build a simple and adjustable mathematical model to generate and analyze torus bifurcations using evident geometric constraints. Since such tori emerge near the fold of the tonic-spiking manifold, the first step is to create the desired number of such folds in the 2D normal-form fast subsystem given by

with a polynomial function $Q=ε+∑n=1∞ln(x2+y2)n$⁠. When $Q≡0$⁠, Eq. (7) describes a linear harmonic oscillator that has a single equilibrium state—the center at the origin in the $(x,y)$-phase plane. When $Q=ε<0$⁠, this linear system has a stable focus at the origin, which becomes unstable when $ε>0$⁠. When $Q(x,y,ε)$ is negatively/positively defined, the nonlinear system has yet a single equilibrium state that globally attracts/repels all other spiral-shaped trajectories. A closed level curve given by $Q(x,y,ε)=0$ (other than the origin) corresponds to a limit cycle of the system. Such a single curve would correspond to a double or saddle-node limit cycle (LC). Two such nested curves would correspond to two limit cycles surrounding the origin. Specifically, if the origin is stable for $ε<0$⁠, then the inner LC is repelling, whereas the outer one is stable or vice versa for $ε>0$⁠. In the special case $ε=0$⁠, the stability of the limit cycles is determined by the sign of the corresponding Lyapunov coefficient $ln$ of the bifurcating equilibrium state at the origin.

This basically describes the algorithm for devising the desired systems. Let us first consider the model with a single fold. Its equations are given below

with $(x,y)$ being the fast variables and $z$ being a slow (as $μ$ is small) dynamic variable replacing fixed $ε$⁠. When $μ=0$⁠, the first two equations present a normal form for the Bautin bifurcation of the equilibrium state at the origin in the case when $ε=0$ and first Lyapunov coefficient, $l1=0$⁠, of the cubic terms (note that $l2$ can be scaled to equal $−1$⁠, for simplicity). The last constraint implies that the outer LC will be stable and a nested inner one will be repelling as $l1>0$⁠. The equilibrium state at the origin is stable as long as $z<0$ and becomes unstable when $z≥0$ after the repelling LC collapses onto it through a sub-critical AH bifurcation at $z=0$⁠. This ends up bi-stability in the fast subsystem and the stable LC is the only global attractor. The left boundary for bistability is given by $zsn=−[l1(x2+y2)+l2(x2+y2)2]$ that corresponds to the double, saddle-node LC in the $(x,y)$-plane. This double orbit vanishes for smaller negative $z$-values, making the origin the only global attractor in the phase space.

The slow equation in (8) is designed here as follows: the slow nullcline, where $z′=0$ is a cylinder given by $x2+y2=1+Δ$ of radius $1+Δ$ in the 3D phase space of the system. Inside this cylinder, $z′>0$ and $z′<0$ outside of it. The use of $Δ$ is twofold here: as a bifurcation and sweeping parameter to locate the 2D manifold M$PO$ in MATCONT.

Now, let us put all components together in the full system. Figure 22 shows the 3D phase space of the 2-layer model [Eq. (8)]: the 1D quiescent manifold $Meq$ (the $z$-axis) with stable (⁠$z<0$⁠) and unstable $(z>0)$ branches. The subcritical AH bifurcation gives rise to the unstable (red) branch $MPOU$ that folds back and continues further (rightwards) as stable (green) branch $MPOS$⁠. The slow nullcline is represented as a gray cylinder parallel to the $z$-axis. Let us pick an initial point to the left quite far away from the fold on $MPOU$ in the 3D phase space. Then, the trajectory spirals fast onto the stable branch (solid line) of M$eq$ inside the slow nullcline/cylinder along which it will be slowly dragged to the right toward the AH bifurcation in the fast subsystem. Passing through it, the phase point keeps following the unstable branch of $Meq$ (the so-called delayed loss of stability) before it spirals away to converge to the stable, outer branch $MPOS$⁠. Provided that the radius of $MPOS$ is larger than that of the slow nullcline, the phase point, turning around $MPOS$⁠, will slide leftward toward the fold on the 2D spiking manifold. If $Δ$ controlling the radius of the slow nullcline is large enough, then the phase point will converge to a stable and round periodic orbit on $MPOS$⁠. Varying $Δ$ makes this periodic orbit slide along M$PO$ thus revealing its shape, including the fold. This is the essence of the parameter continuation approach to localize slow-motion manifolds of the full system by varying the control (faux) parameter of the slow subsystem in the phase space. By decreasing $Δ$⁠, the stable PO slides down toward the fold below which the orbit continues as repelling both unstable in the $z$-direction and in the fast (⁠$x,y$⁠)-coordinates. As we pointed out earlier, the transverse intersection of M$PO$ with the slow nullcline/cylinder guarantees that the PO does not vanish but loses stability at the fold through a super-critical torus bifurcation that results in the emergence of the stable torus-canard (provided $μ$ is small enough) switching back and forth between the stable and unstable branches of the tonic spiking manifold near its fold; see Figure 22(⁠$b1$⁠). Further decreasing $Δ$ makes the stable torus morph into the “elliptic” buster with amplitude modulations and episodes of quiescent meta-states.

Our last 3-layer (torus-relaxation) model below