Bottom-up approach to torus bifurcation in neuron models
BBottom-up approach to torus bifurcation in neuron models
Huiwen Ju, Alexander B. Neiman, and Andrey L. Shilnikov
1, 3 Neuroscience Institute, Georgia State University Department of Physics and Astronomy, Ohio University Department of Mathematics and Statistics, Georgia State University (Dated: 31 May 2018)
We study the quasi-periodicity phenomena occurring at the transition between tonic spikingand bursting activities in exemplary biologically plausible Hodgkin-Huxley type models ofindividual cells and reduced phenomenological models with slow and fast dynamics. Usingthe geometric slow-fast dissection and the parameter continuation approach we show that thetransition is due to either the torus bifurcation or the period-doubling bifurcation of a stableperiodic orbit on the 2D slow-motion manifold near a characteristic fold. We examine varioustorus bifurcations including stable and saddle torus-canards, resonant tori, the co-existenceof nested tori and the torus breakdown leading to the onset of complex and bistable dynamicsin such systems.
Neurons exhibit a multiplicity of oscillatorypatterns such as periodic, quasiperiodic andchaotic tonic-spiking and bursting oscillations in-cluding their mixed modes. The correspondingmathematical models fall into the class of slow-fast systems, which are characterized by the ex-istence and interaction of diverse characteristictimescales. Study of typical transitions betweenthese oscillatory regimes, described in terms ofthe bifurcation theory, is one of the trends inmathematical neuroscience. Quasiperiodic oscil-lations, associated with frequency locking andsynchronization two close or commensurable fre-quencies, typically occur in coupled or periodi-cally forced nonlinear dynamical systems. Thephase space of the system with quasiperiodic os-cillations contains a resonant torus, whose dimen-sion is determined by the number of the interact-ing characteristic frequencies. In slow-fast sys-tems quasiperiodic oscillations emerge throughnonlinear reciprocal interactions when the fastsubsystem experiences a bottle-neck effect equat-ing its time scale to that of the slow subsystem.The torus breakdown, being one of the routesto onset of complex dynamics, turns out to alsounderlie a typical mechanism of transition fromtonic-spiking to regular or chaotic bursting in oneclass of slow-fast neuronal models. The study ofquasiperiodicity and torus bifurcation is a chal-lenging task that requires the aggregated use ofcomputational approaches based on non-local bi-furcation techniques. In this paper we follow theso-called bottom-up approach, starting with thehighly detailed Hodgkin-Huxley type models, andending with formal reduced models. Our goal isto demonstrate how various techniques: geomet-rical slow-fast dissection, parameter continuationand averaging, Poincar´e return maps, when com-bined, allow us to elaborate on all fine details ofthe theory of torus bifurcations and breakdownin the illustrative applications.We dedicate this paper to our dear colleagueand friend, Valentin Afraimovich, who passedaway suddenly in February, 2018. He made many fundamental contributions to the theory of dy-namical systems and bifurcations, including hisco-works on quasiperodicity that set the stage forthe current understanding of this phenomenon indissipative nonlinear systems.
I. INTRODUCTION
Deterministic modeling of oscillatory cells has beenoriginally proposed within a framework of slow-fast dy-namical systems written in a generic form: x (cid:48) = F ( x , y, α ) , y (cid:48) = µG ( x , y, α ) , (1)where x ∈ R n , n ≥
2, and y ∈ R represent fastand slow variables, respectively, | µ | (cid:28)
1, and
F, G aresome smooth functions with α being a parameter vec-tor. For sake of simplicity and convenience, the func-tions in mathematical models are often chosen so that G ( x, y ) = ( f ( x ) − y ) to have a Jordan block in the lin-earization matrix, and G is bi-linear in its arguments.Geometric approach to study of three-dimensional (3D)models of bursting neurons was proposed and developedby J. Rinzel ? . The essence for geometric understandingof neuronal dynamics is in the topology of the so-calledslow-motion or critical manifolds in the phase space ofthe corresponding slow-fast model. The slow-fast dis-section in the singular limit, µ = 0, freezes the slowestvariable, y , which is then treated as a control parameter.This allows detection and parameter continuation in y of the one-dimensional (1D) slow-motion manifold [givenby F = 0] of equilibria (eq), and 2D branches of periodicorbits (PO) of the fast subsystem.These manifolds, called quiescent M eq and tonic-spiking M PO in the neuroscience context, as some back-bones shape the kinds of bursting activity in typicalmodels of individual neurons both formal and derivedthrough the Hodgkin-Huxley formalism. Moreover, onecan loosely describe the type of bursting using the bifur-cations that initiate and terminate the manifolds qui-escent M eq and tonic-spiking M PO ? . For example,“square-wave” bursting shown in Fig. 9 can be alter-natively called “fold-homoclinic”, while “square-wave” a r X i v : . [ n li n . C D ] M a y FIG. 1. 3:1-resonant torus with a pair of periodic orbits, sta-ble and repelling, corresponding to stable and saddle period-3 points on an invariant circle (IC) around a repelling fixedpoint inside, of the Poincar´e return map on a 2D cross-section. bursting depicted in Fig. 13 is due to the fold/saddle-node of equilibria and the fold of periodic orbits.The dissection method, allowing a clear geometric de-scription and understanding the dynamics of low- andeven high-dimensional models, is also helpful to pre-dict possible transitions between activity types, e.g.tonic-spiking and bursting. While the slow-fast dissec-tion paradigm works well for most models with distincttimescales, it provides less or even misleading insights tounderstand transitions between activity patterns, whichare due to reciprocal interactions of slow and fast dy-namics in both subsystems even for small but finite µ -values, which is typical for models of individual neuronsdiscussed below. Such reciprocal interactions become notonly non-negligible but pivotal for full understanding ofboth formal and biologically plausible models in ques-tion when dynamics of the fast subsystem slows down tothe time scale of the slow dynamic. This occurs univer-sally when the fast subsystem undergoes an Andronov-Hopf (AH) bifurcation so that the rate of convergenceto a critical equilibrium state is no longer exponentiallyfast but logarithmically slow, or when its equilibria un-dergo through a saddle-node/fold bifurcation with a longdwelling time due to a bottle neck effect. More complexyet very interesting cases include non-local bifurcationsthrough which the fast subsystem possesses a stable orbitof large period that is about to become a homoclinic orbitto a saddle or to a saddle-node (SNIC), or when a pair ofperiodic orbits, stable and unstable, merges and vanishesin the phase space through the saddle-node bifurcation.One of the well-known examples of such slow reciprocalinteractions is related to the phenomenon of canard os-cillations of small amplitudes in a relaxation oscillatormodel (also known as the FitzHugh-Nagumo model incomputational neuroscience ? ). The canard oscillationsemerge through an AH bifurcation of an equilibrium statewith characteristic exponents close to zero in the norm( ± iω ∼ µ ). The shape of these small oscillations is dueto the fold of the characteristic S -shaped branch (calledthe fast nullcline and given by F = 0) of equilibria of thefast subsystem (see e.g. Fig. 20 A - B ). The equilibriumof the full system is located at the transverse intersec-tion of the fast nullcline with the slow ( | µ | (cid:28)
1) nullclinegiven by G = 0. Note that if the intersection of the slowand fast nullclines is not transverse, i.e. when ∇ F ||∇ G ,then the equilibrium states in the full system bifurcatethrough a saddle-node. The stability of the canard orbitemerging from a preserving equilibrium state is deter- mined solely by the fast subsystem in the singular limit, µ = 0, i.e., it depends on whether the AH bifurcationis sub- or super-critical. The kind of AH bifurcation isdetermined by the sign of the first Lyapunov coefficient l ? , which is related to the smoothness of the function f ( x ). More specifically, given that f (cid:48) = 0 at the fold, andthe constant concavity of the graph of f ( x ) at the equi-librium state in question, then the sign of l is related tothe sign of the Taylor expansion coefficient of f (cid:48)(cid:48)(cid:48) at thecritical equilibrium state at the bifurcation. The choiceof the function, for example x , x , 1 / ( x +1), etc, whosegraphs have such folds locally, will result in either a sub-or a super-critical AH bifurcation. The same holds truewhen one considers the difference equations with discretetime variable, t ∈ Z + , rather the differential one, writtenin a similar form like Eqs. (1) to analyze the local stabilityof invariant circle-canards emerging from a stable fixedpoint ? ? . It happens often that the sign of l derivedthrough the fast subsystem in the singular limit µ = 0changes to the opposite in the full system even at smallbut finite µ values. This indicates that l value is nearzero and that the the system is close to the codimension-2Bautin point, where l = 0, in the parameter space. Thisindicates also that the slow-fast dissection in the singularlimit is only a first order approximation for comprehen-sive understanding the full dynamics of a neuronal modelwith less desperate time scales.The feature of the Bautin bifurcation is that its un-folding includes a saddle-node bifurcation curve corre-sponding to a double periodic orbit that either disap-pears or de-couples into stable and unstable (repellingin when x ∈ R and saddle in higher dimensions) orbitsmerge. As the result, the tonic-spiking manifold M PO inthe phase space of the system has a distinct fold (likesones shown in Figs. 9,12-15, and 22-23 below) at a mergerof stable M SPO and unstable, M
UPO , branches of the tonicspiking manifold, M PO . In slow-fast systems with a singleslow variable, the existence of such a fold is a prerequi-site, but not a guarantee, for the torus bifurcation at thetransition between tonic-spiking and bursting activity.A canard-torus emerges from a periodic orbit with apair of multipliers e ± iφ ( φ ∼ µ ) on a unit circle. How-ever, the question about its stability is more challeng-ing compared to that of a limit cycle canard emergingfrom an equilibrium state with a pair of complex conju-gate exponents ± iω ( ω ∼ µ ), because the torus stabilitycannot be determined a priori by a local stability anal-ysis. En route to bursting the stable periodic orbit, as-sociated with tonic spiking activity M SPO , may loose itsstability via a period doubling bifurcation (single multi-plier e iπ = − ? ? .In the case of the torus bifurcation, the emerging sta-ble or unstable (repelling in R or saddle in R ) toruscan be resonant or ergodic, depending on the angle φ of its multipliers on an unit circle. For example, strongresonances such as 1:4 and 1:3 occur when φ = π/ φ = 2 π/
3, respectively. The torus is born ergodic when φ is not commensurable with π . The example of a reso-nant 3:1 torus with a rational Poincar´e winding number3/1 is shown in Fig. 1. In simple terms this means that FIG. 2. Torus breakdown unfolding including the resonancezone that originates from the torus bifurcation curve ( µ = 0)and bounded by the saddle-node (SN) bifurcation curves forfixed points (FP) on the invariant circle (IC). At larger µ ,ICs become non-smooth first, causing the torus breakdown,followed by a period-doubling bifurcation (PD) of the stablefixed point (FP), and formations of homoclinic tangles (HT)of the saddle FP, and of the saddle-node FP on the SN-bordersof the resonant zone. there is a pair of periodic orbits, stable and unstable,on the 2D torus that correspond to a pair of periodicorbits of period-3 on the invariant circle in some localcross-section transverse to the torus. The winding num-ber remains constant within a synchronization zone (seeFig. 1), also known as an Arnold tongue, and takes vary-ing irrational values outside of the synchronization zone.In this example of the resonant torus, the frequency lock-ing ratio 3/1 means that the “horizontal” oscillations arethree times faster than the “vertical” slow-componentsdue to the timescales of x and y -variables of Eqs. (1).Andronov and Vitt ? were first to study synchroniza-tion or phase-locking in the periodically forced van derPol equation ? . They showed that the wedge-like shapeof each resonant zone is bounded by the saddle-nodebifurcation (SNIC) curves in the parameter plane sothat having crossed one, the resonant torus becomes er-godic with everywhere dense covering (trajectory) on it.Later, Afraimovich and Shilnikov ? ? ? ? proposed a phe-nomenological scenario of torus breakdown. It is illus-trated in Fig. 2, which depicts pivotal evolution stages ofthe invariant circle of a 2D Poincar´e return map and fixedpoints on it inside and outside of the resonant zone in theparameter plane. Here µ = 0 corresponds to the torusbifurcation and µ controls the angle of the multipliers ofthe bifurcating periodic orbit within the range from zerothrough π ( e ± iφ ). At small µ values, outside the reso-nant zone, the invariant circle (IC) (read the 2D torus)is initially smooth and rounded. It possesses a saddle-node fixed point (FP) on SN-curves so that the closureof its 1D unstable set constitutes the IC. Inside the res-onant zone, there are two FPs, stable and saddle, on theIC. As µ is increased the IC becomes “non-smooth” sothat the stable FP later undergoes a period-doubling bi-furcation (PD curves) that breaks the torus down. Next,the unstable sets of the saddle-node or saddle FP touches (HT-curves) and crosses its stable sets and creates homo-clinic tangles giving rise to complex chaotic dynamics inthe system. The Reader is welcome to consult with thereview ? and the reference texts ? for more details aboutthe torus breakdown.In-depth understanding of the universal mechanisms oftransitions between oscillatory activity patterns in singleneuron models is a challenge for the theory of bifurca-tions. It requires the use of wide range of advanced appa-ratus of the bifurcation theory, ranging from the blue skycatastrophe to various homoclinic bifurcations of saddleorbits and equlibria ? ? ? ? ? ? ? ? ? ? ? ? ? ? .In what follows we will present several examples ofbiologically plausible and formal mathematical slow-fastmodels to demonstrate the universality and complexityof the the torus bifurcations near the distinct fold onthe tonic-spiking manifold. We use the bottom-up ap-proach, starting with highly detailed models and endingwith simple toy models, all featuring same generic prop-erties. These bifurcations can result in the onset of sta-ble, repelling and saddle tori as well as their coexistencein the phase space leading to bi-stability of tonic spikingand bursting activities. We discuss the mechanisms oftorus breakdowns through various resonances, and con-ditions for period-doubling bifurcations alternatively oc-curring on the same fold. We note the torus-canards havebecome a new trend in mathematical neuroscience witha focus on transition between tonic-spiking and burstingin various slow-fast models of individual neurons see andreferences therein ? ? ? ? ? ? ? ? . II. A MODEL FOR SPONTANEOUS VOLTAGEOSCILLATIONS IN HAIR CELL
Our first example is a detailed Hodgkin-Huxley typemodel, whose 12 ODEs describe spontaneous dynamicsof membrane potential of a hair cell. Hair cells are pe-ripheral receptors which transform mechanical stimuliinto electrical signals in the senses of hearing and bal-ance of vertebrates. These sensors rely on nonlinear ac-tive processes to achieve astounding sensitivity, selectiv-ity and dynamical range ? . In particular, it was proposedthat ciliary bundles of hair cells may self-tune to oper-ate on the verge of (degenerate) AH bifurcations caus-ing a sharpened selectivity to weak periodic mechanicalforcing ? ? .In amphibians some hair cells demonstrate variousspontaneous oscillatory activities. In particular, severalexperimental studies have documented various oscillationregimes of membrane potential in the frog sacculus ? ? .Although physiological implications of the membrane po-tential oscillations on sensory function of hair cells arenot clear, experimental study ? showed that spontaneousvoltage oscillations drive sensory neurons resulting in pe-riodic firing. Furthemore, the membrane potential of thehair cell affect the dynamics of its apical hair bundle com-partment, as documented experimentally ? and in mod-eling work ? ? . Modeling work ? ? ? predicted that oscil-latory dynamics of the membrane potential may lead toimproved sensory performance of the hair cell. FIG. 3. Dynamical regimes of the hair cell model. A ( b, g K1 ) bifurcation diagram of the hair cell model is superimposedwith a color-coded map indicating the spike number per burst. It includes Andronov-Hopf (red) curves of supercritical (solidline) and subcritical (dashed line) bifurcations including a codimension-2 Bautin point (filled circle, BP) that bounds theregion of oscillatory activity of the hair cell. Black dashed line corresponds to the torus bifurcation (TB), and solid green linecorresponds to a period-doubling (PD) bifurcation of the periodic orbit (PO). B Gallery of voltage traces for the indicatedincreasing g K1 -values at fixed b = 0 . A. Hair cell model
The model is based on several experimental studiesof basolateral ionic currents in saccular hair cells inbullfrog ? ? ? ? and is a modification of a model devel-oped in ? . It includes 6 ion currents plus a leak currentwith corresponding equations for kinetics of ion channelsand for dynamics of Ca concentration, resulting in aHodgkin-Huxley type system of 12 coupled nonlinear or-dinary differential equations. A detailed description ofthe model is available in an open access publication ? anda code for the right hand sides of corresponding ODEs isprovided in the Appendix .The membrane potential dynamics in the hair cellmodel is described by this equation: C m V (cid:48) = − I K1 − I h − I DRK − I Ca − I BKS − I BKT − I L . (2)These six ionic currents can be grouped according totheir activation patterns. The K + inward-rectifier ( I K1 )and the cation h-type ( I h ) currents are hyperpolariza-tion activated. The rest four currents are depolariza-tion activated: the voltage activated K + direct-rectifier( I DRK ) and Ca ( I Ca ) currents; the calcium-activatedK + steady ( I BKS ) and transient ( I BKT ) currents. Finally,the leak current ( I L ) stands for the mechano-electricaltransduction from the hair bundle compartment of thecell. Previous computational study ? showed that a con-venient choice of control parameters is the maximal con-ductance of K1-current, g K1 , and the strength of theCa -activated potassium currents, b .Experimental work showed that depolarization acti-vated currents are responsible for the so-called phe-nomenon of electrical resonance, whereby the haircell shows fast oscillatory responses when knocked byan external current pulse ? ? ? ? ? . Furthermore, aself-sustained voltage oscillations were observed ex-perimentally and in a model which contained the voltage-activated Ca and the Ca -activated potas-sium currents ? . On the other hand, hyperpolariza-tion activated currents give rise to a slow (3–4 Hz)large-amplitude oscillations of the membrane potential,owing to slow kinetics of I h current ? . Taken to-gether, the membrane potential of the frog hair cellsshow diverse patterns of oscillatory activity documentedexperimentally ? , and well captured by the model ? . Inparticular, the two distinct oscillatory mechanisms men-tioned above may lead to quasi-periodic oscillations, witha torus bifurcation and torus breakdown ? , which will bestudied here in details. B. Quasi-periodicity & period-doubling pathways
To unravel the dynamical mechanisms behind transi-tion to quasiperiodic oscillations we built a bifurcationdiagram of the model for two control parameters, g K1 and b . The diagram shown in Fig. 3 was constructed us-ing the combined techniques of parameter sweeping andcontinuation.In Fig. 3 A , the AH bifurcation (red curve) demarcatesthe region of the quiescence state from the oscillatory dy-namics. The solid and dashed segments of the AH linecorrespond to super- and subcritical AH bifurcation. Thecriticality of the bifurcation means that either a stable ora saddle periodic orbit emerges and collapses into the sta-ble equilibrium state, after the corresponding segment ofthe AH-curve is crossed. The point, labeled BP, is wherethe criticality of the AH-bifurcation changes. This bifur-cation was first studied by N. Bautin ? who found thatits unfolding includes another curve corresponding to adouble, saddle-node periodic orbit around the equilib-rium state. The orbit results from a coalescence of thestable and saddle ones, that have emerged from a stablefocus through super- and sub-critical AH bifurcations, re- g K1 (nS) A B g K1 (nS) V ( m V ) FIG. 4. A Bifurcation diagram representing the g K1 -parameter sweep of V min -values reveals the stability loss of the tonic-spiking periodic orbit through a torus bifurcation en route to bursting at level b = 0 .
01. Note that the ergodic torus becomesresonant before its breakdown into large-amplitude bursts. B Diagram showing the cascade of period-doubling bifurcations oftonic-spiking orbits transitioning to bursting at level b = 0 . spectively, upon crossing the solid and dashed branchesof the (red) curve, AH, in the diagram in Fig. 3 A .A synonym of the saddle-node bifurcation of POs is afold where two branches, foliated by stable and saddleperiodic orbits, merge. Such a fold occurs on a 2D tonic-spiking manifold M PO that is traced down by the periodicorbits of the full system using the parameter continuation(more about this method in ? ). The notion of the foldis imperative here, as it is where a round tonic-spikingperiodic orbit loses stability and tonic spiking transformsinto bursting activity though either a series of perioddoubling bifurcations or through a torus bifurcation.There is a novel dynamic feature that makes the haircell model stand out in the list of conventional modelsof bursters. Namely, it is due to the fold on the tonic -65.4 -65.2 -65.0-65.4-65.2-65.0 -66.0 -65.5 -65.0 -64.5-66.0-65.5-65.0-64.5-67 -66 -65 -64-67-66-65-64 -67 -66 -65 -64-67-66-65-64 V n (mV) V n (mV) V n + ( m V ) V n + ( m V ) V n + ( m V ) V n + ( m V ) A BC D
FIG. 5. Poincar´e return maps, V ( n )min → V ( n +1)min , depictingfour pivotal stages of the torus onset and breakdown on thepathway b = 0 . A Convergence to a spiraling (reso-nant) fixed point (FP) corresponding to a tonic-spiking orbitat g K1 = 28 .
335 nS. B Stable smooth IC corresponding toa stable ergodic torus at g K1 = 28 .
345 nS with a repellingFP inside. C A non-smooth (distorted) IC for the yet er-godic torus at g K1 = 28 .
355 nS. D A stable 7-period orbit at g K1 = 28 . spiking manifold, M PO , that results from the saddle-node bifurcation of periodic orbits, originating in turnfrom the codimension-2 Bautin bifurcation on the AHcurve in the diagram in Fig. 3. The fold per se is nodynamic feature in slow-fast systems, unless the stableperiodic orbit for tonic oscillations slides down closer tothe fold to lose the stability to the 2D torus ? . Thismechanism, first described in ? , happens to occur typi-cally in slow-fast systems with such folds ? , particularlyin neuronal models ? ? including the reduced model ofPurkinje cells ? below.A limit cycle born through the AH bifurcation loses itsstability either through a torus bifurcation (TB, dashedblack line in Fig. 3 A ) or through a period-doubling bi-furcation (PD, solid green line in Fig. 3 A ). A regionbounded by these lines is characterized by various burst-ing patterns and by period-adding bifurcations. The col-ormap in Fig. 3 A , showing the number of spikes perburst, indicates that the number of spikes increases to-wards the torus bifurcation. To the right of the PD-curve, i.e. for 0 . < b <
1, the model produces robustlytonic oscillations ? , which underlines stabilizing physi-ological role of the Ca -activated potassium currents( I BKS , BKT ) ? .Next we investigate the torus bifurcation and the torusbreakdown using the slow-fast manifolds and Poincar´ereturn maps on a route from tonic spiking to bursting.We also examine the transition to bursting on anotherroute featuring the period-doubling bifurcation, and con-trast the corresponding return maps in both cases. ThePoincar´e return map is well suited for exploration of com-plex oscillations, revealing the mechanisms of underlyingbifurcations and their precursors. Since the voltage isthe only observable variable in many experimental se-tups, we will employ the Poincar´e return map T definedon consecutive minima of the membrane voltage trace asfollows: T : V ( n )min → V ( n +1)min . In the case of a sparse map(insufficiency of various points V min ), e.g., for periodic orweak chaotic bursting, a small random perturbation canbe added to the model, resulting in variability sufficientto obtain densely populated iterates, revealing some keydynamical features. We point out that sensitivity of so-lutions to external perturbations is a common feature ofslow-fast systems like neuronal models ? .The described two scenarios of the transition from -67 -67 -65 -65 -66 -66 -64 -64 -63 -75 -85 -65 -55 -75 -85 -65 -55 -63 V n + ( m V ) V n (mV) V n (mV) A B
FIG. 6. Poincar´e return maps, V ( n )min → V ( n +1)min , for the con-secutive V min -values in voltage traces generated by the haircell model. A Evolution of stable invariant circles (IC) fromergodic to resonant with further non-smooth torus break-down as the g K1 parameter is increased from 29.185 through29.2073 nS. B Chaotic bursting after the torus breakdownat g K1 = 29 .
213 nS. The flat, stabilizing section of the mapcorresponds to the hyperpolarized quiescence, while multiplesharp folds reveals a ghost of the non-smooth IC in the depo-larized range. tonic-spiking to bursting (through the torus or period-doubling bifurcations) can be differentiated with the aidof the Poincar´e maps for consecutive minimal values ofthe membrane voltage traces shown in Fig. 4.
C. From tonic spiking to bursting throughquasi-periodicity
Let us first examine the torus bifurcation scenario, as g K1 increases with a fixed value of b = 0 .
01. Typicalevolution of the voltage traces is shown in Fig. 3 B . Asthe parameter g K1 increases a small-amplitude periodicoscillations become quasi-periodic (slow amplitude mod-ulation) through the supercritical (stable) torus bifurca-tion. When the torus brakes down, it gives rise to large-amplitude bursting oscillations which may be chaotic orregular in alternation. This is illustrated in the bifur-cation diagram of Fig. 4 A , which was built using thePoincar´e return map of consecutive voltage minima. Inthis diagram, a single point for a given parameter valuecorresponds to a stable periodic orbit. Several pointsper a single parameter value may correspond to differ-ent regimes. Dense points above −
70 mV correspond toergodic tori, whereas multiple extended branches corre-spond to stable periodic orbit(s) on the resonant torusexiting within its resonant zone. Points below −
70 mV,refer to bursting orbits whereby the cell becomes stronglyhyperpolarized. Bursting orbits are stable when theirseveral branches in the bifurcation diagram are continuedwithout interruption within a stability window, followedby regions of chaotic bursting.One can see from the bifurcation diagram in Fig. 4 A that the stable periodic orbit, earlier emerging throughthe supercritical AH bifurcation, loses its stability when g K1 is increased through 29 .
19 nS. This stability loss givesrise to the emergence of a 2D stable torus, first ergodic,as a widening parabola-shaped solid region in the dia-gram suggests. In what follows, the breakdown of thetorus comes along with the generic scenario lines pro-
FIG. 7. Poincar´e return maps, T : V ( n )min → V ( n +1)min , depict-ing strong-resonant ICs (tori) as g K1 -parameter increases atvarious fixed b -values. Insets A - C illustrate the evolution ofinitially ergodic IC as it becomes strongly resonant, 1:4, 1:3and 1:2 at b = 0 . , . . D Weak reso-nant case at b = 0 . A – B
3D extended evolution of the ICs with in-creasing g K1 ; grey dots represent V min -values extracted fromvoltage traces; black contour lines represent the ICs of thereturn maps at specific g K1 -values; blue dots are connectedto show the iteration order; green and orange branches revealthe PD bifurcations of resonant 3- and 4-period orbits. posed by Afraimovich and Shilnikov ? . The torus be-comes resonant as g K1 increases through a synchroniza-tion boundary (of an Arnold tongue) ? ? . This asser-tion is supported by the occurrence of several branchescorresponding to the local minima of a stable periodicorbit coexisting with a saddle orbit on the torus aftera saddle-node bifurcation at g K1 (cid:39) .
199 nS. Furtherincrease of the control parameter causes the torus break-down, which leads to bursting oscillations with modula-tory (quasi-periodic) tonic spiking episodes interruptedby slow hyperpolarized quiescent periods (Fig. 3 B at g K1 = 29 .
213 nS). Fig. 4 A also reveals a number ofspike adding bifurcation of bursting. The cascade of spikeadding ends on the upper branch of the green bifurcationcurve, that corresponds to a reverse period-doubling bi-furcation (Fig. 3 A ). To the right of the curve, the haircell model produces stably a large-amplitude voltage os-cillations (Fig. 3 B at g K1 = 40 nS). The black dashedcurve in Fig. 3 A corresponds to the torus formation andtherefore quasi-periodic tonic spiking oscillations in themodel.Figures 5 and 6 represent the Poincar´e return mapsdepicting the characteristic stages of the torus formationand breakdown. In Fig. 6 A , a family of invariant cir-cles (ICs) corresponds to the growing of the first smoothand later non-smooth ergodic stable IC emerging froma tonic-spiking periodic orbit, as the latter has lost “itsskin” to the 2D torus (Fig. 5 A and B ). As g K1 increases,the torus becomes resonant with a stable periodic orbiton it. The degree of the resonance can be evaluated bythe number (here 7) of periodic points of the stable or-bit, which equals the number of the ”sleeves” originatingfrom an unstable FP towards the non-smooth invariantcurve. This agrees well with the theory of torus break-down that causes the onset of complex chaotic dynamics.In essence, the torus breakdown begins with the emer-gence of two periodic orbits, stable and saddle, on thetorus when it becomes resonant through a saddle-node bi-furcation. Later the torus becomes non-smooth after thestable and unstable manifolds (sets) of the saddle orbitstart making wiggles transitioning to homoclinic tangleswith the stability loss of the counterpart periodic orbitthrough a period-doubling bifurcation. This bifurcationcan be identified through end-branching on the resonanttraces left by the periodic-point as the control parameter g K1 is increased. One can observe from Fig. 5 D that thestable 7-period orbit has a pair of leading complex multi-pliers, which indicates the proximity of period-doublingbifurcations (green branches indicated in Fig. 7 below).We note that bursting can already co-exist with thestable torus in the hair cell model, so the torus breakdowngives rise to bursting, first regular and then chaotic at g K1 = 29 .
213 nS. The corresponding return map is shownin Fig. 6 B . It has a distinct shape, so characteristic forthe square wave bursters ? ? . The map has three compo-nents: the transient “toroidal” section rivaling a ghost ofthe disappeared non-smooth IC around V min = − min = − g K1 is depicted inFig. 3 B at g K1 = 29 .
213 nS; it demonstrates slowly mod-ulated chaotic tonic-spiking phases alternating with hy-perpolarized voltage sags.The resonance of a torus is determined by its windingnumber, calculated as the ratio of 2 π and the angle φ of the Floquet multipliers e ± iφ on a unit circle ? . Forexample, when the angle is π/ A ), the windingnumber is 4, meaning that the resonance of the torus is1:4. The irrational and rational winding numbers corre-spond to ergodic (quasi-periodic) and resonant (periodic)torus, respectively. The cases 1:1, 1:2, 1:3, 1:4 are consid-ered as strong resonances and others as weak. Figure 7shows the emergence, evolution and breakdown of the 2Dtori through three strong resonances, namely 1:2, 1:3, 1:4,and one weak case. When the torus becomes resonant,the corresponding IC (grey color) of the Poincar´e returnmap possesses stable periodic orbits (shown as blue dots).The resonant ICs becomes non-smooth and then startsbreaking down when the leading multipliers of the reso-nant periodic orbits become complex conjugate towards -85 -80 -75 -70 -65 -60-85-80-75-70-65-60 -85 -80 -75 -70 -65 -60-85-80-75-70-65-60 A B V n + ( m V ) V n (mV) V n (mV) FIG. 8. Period-doubling shown in Poincar´e return maps, V ( n )min → V ( n +1)min , with multiple unimodal branches. A Chaoticbursting in the deterministic model at g K1 = 26 . b = 0 . B Weak ( σ = 5 × − ) noise induces burstingthat overshadows a period-2 obit (shown by two green dots)at g K1 = 25 .
972 and b = 0 . forthcoming PD bifurcations (shown as green and orangedots in the return maps) with the increase of the parame-ter g K1 at the indicated values of parameter b , so that theunstable sets of the saddle periodic orbits no longer con-verge monotonically but start spiraling onto the stableorbits on the torus (see Fig. 5 D ). FIG. 9. A
3D (m h , m DRK , V)-phase space projection depict-ing a few exemplary tonic-spiking, quasi-periodic and burst-ing orbits being superimposed on the slow-motion manifolds:1D S-shaped quiescent M eq with two knees/folds, and the 2Dtonic spiking M PO . The fold on M PO corresponds to a saddle-node bifurcation giving rise to periodic and quasi-periodictonic-spiking [torus as dark-blue orbits] occurring on the out-ward stable branch of M PO . The grey surface, m (cid:48) h = 0, is theslow nullcline above/below which the slow variable increas-es/decreases during tonic-spiking/quiescent phases of burst-ing. B Voltage traces with matching colors of the correspond-ing orbits in A . D. Period-doubling en route to bursting
Next we consider the pathway at b = 0 .
015 passingthrough the period-doubling curve (green line) in the 2Dparameter plane of Fig. 3 A . Having crossed this curve,the stable periodic orbit (representing tonic-spiking os-cillations) loses its stability to the one with the doubledperiod. This bifurcation and the new periodic orbit corre-spond to the branching in the diagram shown in Fig. 4 B .A cascade of forthcoming period-doubling bifurcationsdevelops quite rapidly as g K1 increases, and terminateswith the onset of regular bursting in the hair cell model.Spike adding causes the onset of chaotic bursting in thetransition between stable 5-spike and 4-spike bursting asthis bifurcation diagram reveals.The Poincar´e return map of the chaotic bursting shownin Fig. 8 A has a wider drop-shaped section, comparedto chaotic bursting emerging through the torus scenario(Fig. 6 B ). Quantitatively, the width of the drop indicateshow fast the model switches back and forth between tonicspiking and hyperpolarized quiescence in bursting. Thewider the drop in the map is, the less rigid the switch-ing is, or, alternatively, the less contrast between fastand slow dynamics is. Since the switches are quite loose,not constrained, there is a wide variability in transitionsthat give rise to multiple branches of the return map inFig. 8 A . The multiple branches look like as they repre-sent snapshots of some return maps at several parame-ter values in progression. Each branch reveals a genericunimodal (non-sharp) shape of the return map that isneeded for a cascade of period-doubling bifurcations tooccur. We note that other Hodgkin-Huxley type mod-els with the fold on the tonic spiking manifold can ex-hibit a period-doubling cascade leading to a chaotic tonic-spiking attractor in the phase space ? . The progressionof the branches also suggests that the stable periodicorbit emerges through a saddle-node bifurcation whenthe graph of the map touches the 45 ◦ -degree line. Re-call that the corresponding bifurcation comes out of thecodimension-two Bautin point in the parameter plane inFig. 3 A . Next, as the map graph lowers further and be-comes steeper, the stable fixed point becomes unstablethus giving rise to a period-two orbit, as such indicatedby two green solid circles in Fig. 8 B . The bottom fourbranches of the map correspond to bursting where eachhyperpolarized voltage drop is followed by chaotic tonic-spiking phases developed through all steps of the perioddoubling cascade.Precursors of bursting can be revealed by adding smallrandom perturbations to the model. Gaussian whitenoise ξ ( t ) is injected to the right hand side of Eq.(2)with the term σξ ( t ), where σ is the standard deviation ofrandom perturbation. Fig. 8 B demonstrates that weakrandom perturbations can effectively reveal precursors oftransition to bursting. For the value of g K1 , the modelpossesses a stable period-2 orbit and so the Poincar´e re-turn map contains two points only. Weak noise inducesbursting resulting in the map which has a similar shapeas the deterministic chaotic bursting.We stop here to resume that without 1D Poincar´ereturn maps generated through voltage oscillations, it would be impossible to characterize the exact forma-tion mechanisms underlying transitions from simple tonicspiking to bursting and back in this model as well as inother Hodgkin-Huxley type neuronal models ? ? . FIG. 10. Scaling down the slow m h -variable shifts the emer-gent torus closer to the fold and changes its stability in thehair cell model. A – B
3D (m h , m DRK , V)-phase space pro-jection show the stable torus (blue line) emerging from a PO(yellow ring, detected by MATCONT) on stable section ofM PO , for the original and 0.55 slowed rate of m h -gating vari-able, resp. C : Further decreasing the rate by the factor oftwo makes the unstable torus (of the saddle type) emergingat the fold on M PO . The saddle torus creates bistability asbounds the stable tonic orbit (green circle) away from burst-ing (see Fig. 11). Insets A – C depict the correspondingvoltage traces of the tonic-spiking attractors above. E. Parameter continuation & slow-fast dissection
Bifurcations underlying the transitions between differ-ent oscillatory patterns in the hair cell model can be ex-plained geometrically using the slow-fast dissection as il-lustrated in Fig. 9.The essence of the approach is based on dissecting theoriginal dynamics of the model into slow and fast compo-nents. In the model the activation of the h -current, m h ,is the slowest dynamical variable, while all other vari-ables are fast compared to it. In particular, we choosethe membrane potential, V , and the activation of theK + DRK-current, m DRK , as representatives of the fastsubsystem. The feature of a slow-fast system that oftenreduces the complexity of its dynamics is that the phasetrajectory stays close to the slow-motion manifolds of lowdimensions: a 2D tonic-spiking manifold, M PO , foliatedby small-amplitude periodic orbits, and a 1D S -shapedquiescent manifold, M eq , comprised of hyper- and depo-larized equilibrium states of the model. This approachimplies that bursting, being a multiple time-scale pro-cess, is interpreted as repetitive fast switching betweenthese manifolds at their terminal phases, followed by slowpassages of trajectories turning around M PO and thenalong M eq , corresponding to the episodes of tonic-spikingand hyperpolarized quiescence, respectively. Figure 9 de-picts the topology of these slow-motion manifolds for theparameters close to the torus (dark blue trajectory) bi-furcation, as well as shows several bursting orbits su-perimposed (green, light blue and brown trajectories) inthe phase-space projection. These spiking and quiescent FIG. 11. A
2D (m h , V)-phase space projection shows bistabil-ity of tonic spiking orbit (green circle, PO) and bursting orbit(brown) being overlaid on the quiescent manifold, M eq , andthe tonic-spiking manifold, M PO , when the m h -gating vari-able of the hair cell model is scaled down two times slowerthan its original rate. These two stable states are separatedby a saddle torus, denoted here by the transient part of thetonic spiking orbit (pink). Inset B shows the voltage tracescorresponding to the tonic spiking and bursting orbits in A . manifolds are obtained by the parameter continuation ofperiodic orbits and equilibria of the model with the aux-iliary parameter shifting the slow nullcline m (cid:48) h = 0 upand down.In what follows we outline how the manifolds are de-tected in the 12D phase space of the hair cell model.We use a new and practical way for localization of bothslow motion manifolds, 1D quiescent, M eq , and 2D tonic-spiking, M PO , composed of equilibria and periodic orbitsof the full system. The method capitalizes on the slow-fast dissection as well as on the parameter continuationtechnique ? . The methods supplements the conventionalslow-fast dissection in the singular limit with the slowvariable frozen as a parameter. The main essence is theparameter continuation technique applied to the entireset of the model equations. The advantage of our ap-proach is that it yields the sought manifolds themselves inthe phase space of the intact model rather the manifoldsof its fast subsystem without taking into account slowyet finite component of the overall dynamics. This ap-proach is especially valuable for complex models of higherdimensions where slowfast dissections could be problem-atic because of multiple time-scales various ionic currentsinvolved in complex dynamics especially at transitions.Furthermore, this approach can detect and explains bi-furcations, often non-local due to reciprocal interactionsof slow and fast dynamical variables that underlie thesetransitions in the model in question.Specifically, we introduce a faux control parameter,which is a voltage offset that affects only the slowest dy-namics of the model; this offset shifts the voltage thresh-old at which the slow gating variable/channel is half- FIG. 12. A
3D ([Ca ] i , m iNa , V)-phase space projection ofthe pyramidal cell model shows the stable torus (light blueline) on the fold of slow-motion manifolds, M PO , superim-posed on the 1D S-shaped quiescent manifold, M eq , and theslow nullcline, [Ca ] (cid:48) i = 0 (yellow surface). The Poincar´ecross-section transversal to the torus highlights the dark bluedots constituting the stable IC on it. Inset B shows the slowlymodulated voltage trace corresponding to the torus in A . open. Geometrically, its variations lift or lower the slownullcline, m (cid:48) h = 0, in the 3D phase-space projection de-picted in Fig. 9. As the voltage offset changes, the in-tersection point of the slow nullcline (sigmoidal surface)with M eq [the equilibrium state of the model] is shiftedalong M eq , thus tracing M eq down and explicitly reveal-ing its position and shape in the phase space. Same istrue for the other manifold M PO . More details on thisapproach can be found in ? .Like many neuron models, the quiescent manifold inthe hair cell system has a characteristic S -shape in theprojection onto the slow-fast (m h , V)-plane. The shapeof M eq implicitly endows the model with a hysteresis re-quired for dynamic switching between hyper- and depo-larized phases in large-amplitude bursting. Reaching thelow, hyperpolarized fold on M eq indicates the beginningof a burst. In addition, the two real bifurcation param-eters, g K1 and b , of the model can shift the position ofthe manifold M eq relative to that of the slow nullcline,m (cid:48) h = 0 in the phase space, so that the intersection pointcan move onto the stable upper or low section of M eq .In those cases, the hair cell rests on the depolarized orhyperpolarized steady state, respectively.The hair cell model oscillates tonically when there is astable periodic orbit on M PO . This orbit emerges fromthe depolarized equilibrium state through the supercrit-ical AH bifurcation (discussed above). Moving the slownullcline by varying the faux parameter makes the pe-riodic orbit slides along M PO , thus revealing its shape.The manifold ends up through a homoclinic bifurcationoccurring after the periodic orbit of the increased size0 FIG. 13. A
3D (h , M , V)-phase space projection of the Purk-inje cell model at I app = − .
487 depicting the saddle (red)torus-canard slowly oscillating between the stable (green) in-ward M
SPO and unstable (grey) outward M
UPO branches merg-ing at the fold of the 2D spiking manifold. The yellow surfaceis the slow nullcline M (cid:48) = 0. The 1D quiescent manifold M eq (black line) has stable (solid) and unstable (dashed) branches. B Magnified saddle torus-canard T U (red) encloses a stabletonic-spiking PO (green circle). The phase points at maxi-mum and minimum voltage values trace out two circles onthe space space. C Slowly modulated voltage trace in A . and touches the middle, saddle branch of M eq to becomethe homoclinic obit of the saddle. This is a generic con-figuration for square-wave bursters. III. PURKINJE CELL MODEL
Whenever neither the stable branch of M PO nor thestable branches of the quiescent manifold M eq are cutthrough by the slow nullcline, m (cid:48) h = 0, i.e., both man-ifolds are freely transient, the hair cell model exhibitsbursting. Bursting is represented by solutions repeat-edly switching between the low, hyperpolarized branchof M eq and the tonic spiking manifold M PO . The num-ber of complete revolutions of the solution around M PO is that of spikes per burst in the voltage traces. Thiswinding number is used to classify the bursting activity.The larger the number is, the longer the burst durationlasts. The relative position of the torus to the fold onM PO can be changed by the rate of the slow variable m h .The slower the rate, the closer the stable torus emergesto the fold as shown in Figs. 10 A – B , as well as the closerthe torus bifurcation (TB) is detected near the fold bythe bifurcation package MATCONT ? ? . More interest-ingly, when the rate of m h is about two times slower,the torus bifurcation becomes sub-critical and generatesa saddle torus at the fold instead (see Fig. 10 C ). Thisgives rise to bistability in the system, where both tonicspiking and bursting can be produced by the model de-pending on the initial conditions, as illustrated n Fig. 11.A small basin of attraction of the stable tonic-spiking or-bit is bounded by the 2D saddle torus in the phase space FIG. 14. Bistability of tonic-spiking and bursting in thePurkinje cell model. A Bifurcation diagram representedby the bi-directional forward/backward I app -sweeps of V min -values (green/blue dots, resp., and red box for saddle torus-canard). White horizontal bar indicates the range of the co-existence of stable tonic-spiking and bursting, red vertical lineat I app = − . B
2D (V , M)-phase projection depicting the co-existing stabletonic-spiking (green) and bursting (blue) orbits separated bythe saddle torus-canard (red) orbit superimposed on the man-ifold M PO (grey), the slow nullcline M (cid:48) = 0 (yellow), and themanifold M eq (black line) at I app = − . C and D show the voltage traces corresponding to the tonic-spikingand bursting orbits in B , resp. from that of the stable bursting orbit. We note thoughthat similar bistability is also found in the Purkinje cellmodel, parabolic burster model and FitzHugh-Nagumo-Rinzel (FNR) model (see the following sections).In a system with drastically distinct slow and fast timescales, the periodic orbit would bifurcate into a torusor through period-doubling right at the location of thefold ? . The type of the bifurcation depends on globalstructure of the flow at the fold. Specifically, dependingon whether the phase space volume contracts or not, thetonic spiking orbit undergoes period-doubling or torus bi-furcation. While the verification of this condition seemsdoable in 3D slow-fast systems (see the FNR section be-low), the high-order models may present a challenge asone has to identify the set of the variables (from a 3Dsubspace containing the central manifold in which thebifurcation takes place) that are involved directly intothe ongoing bifurcation.1 IV. TORUS IN A PYRAMIDAL CELL MODEL
Our next example is another highly detailed 14DHodgkin-Huxley type two-compartmental (for soma anddendrite) model describing interactions of pyramidal cellsand fast-spiking inhibitory interneurons ? . The readercan find details of all currents constituting the modelin Ref. ? , which describes how electrogenic properties ofthe N a + / K + ATPase control transitions between normaland pathological brain states, specifically during epilep-tic seizures. The equations describing the dendritic andsomatic voltage dynamics are: C m V (cid:48) d = − I intd − I leakd − αI pumpd − g c ( V d − V s ) s d ,g c ( V d − V s ) s d = − I ints − I leaks − αI pumps , (3)where C m is the capacitance of the dendritic compart-ment, V is the membrane potential, I s are various ioniccurrents including Na + /K + pump current, s stands for asurface of each compartment, g c is the coupling conduc-tance of the soma (s) and the dendrite (d), α is a scalingcoefficient.Unlike in the hair cell model, where the dynamics ofintracellular calcium concentration, [Ca ] i , is relativelyfast, here [Ca ] i is the slowest dynamical variable inthe model. As for the hair cell model, we will use the3D phase projection to show the topology of the slow-motion manifolds that determine the shape of burstingoscillations in the model. Figure 12 depicts a similarlyshaped 1D quiescent manifold M eq . Its lower, hyper-polarizing knee corresponds to the homoclinic saddle-node bifurcation giving rise to the stable (green) outersection of the 2D tonic-spiking manifold M PO . This man-ifold wraps back inwardly so that its unstable (pink)shrinking section terminates on the depolarizing branchof M eq through a sub-critical AH bifurcation. The sur-face [Ca ] (cid:48) i = 0 (yellow) is the slow nullcline: [Ca ] i increases/decreases above/below it in the phase space.When this nullcline is slightly below its position depictedin Fig. 12, there would be a stable periodic orbit on thestable branch of M PO that corresponds to tonic-spikingactivity. As the nullcline [Ca ] (cid:48) i = 0 is shifted up, thestable orbit comes close to the fold and looses its stabil-ity through a supercritical torus bifurcation. The stablenewborn torus oscillates back and forth between the sta-ble and unstable branches of M PO near the fold. Thetorus is stable because it lingers longer on the stablebranch then on the unstable one on average. Other-wise it would be of the saddle type, i.e. repelling in them Na and [Ca ] i coordinates and stable in the all otherdimensions. One can observe from Fig. 12 A that gapsbetween the successive points filling in the stable invari-ant curve are large enough to assume that dynamics of[Ca ] i is not that slow compared to the torus-canardobserved Purkinje cell model discussed next. After thetorus breaks down as the slow nullcline is further shiftedup, the pyramidal cell model demonstrates bursting ac-tivity, see ? .Purkinje cells are the main output of the cerebellumand are involved in motor learning/conditioning ? . To FIG. 15. Averaging technique for the Purkinje cell model. A – B (V , M)-projections of the stable M
SPO (green) and un-stable M
UPO (grey) sections of the 2D tonic-spiking manifold,the 1D S-shaped quiescent manifold M eq (black line), slownullcline M (cid:48) = 0 (yellow line) and the 1D curve (cid:104) V (cid:105) (brown)made up of the averaged phase coordinates of the POs consti-tuting M PO . A Basin of the stable tonic-spiking PO (darkgreen) is bounded by the saddle torus, shown in B (red), atI app = − . UPO bend-ing outward from the fold terminates on the middle, saddlebranch of M eq after the homoclinic bifurcation (HB). Insets A – B disclose the graph (black line) of the average func-tion (cid:104) M (cid:48) (cid:105) determining the dynamics of the slow M-variableon M PO ; a single zero (dark green dot) of (cid:104) M (cid:48) (cid:105) and its slopedetermine the position and the stability of the correspondingPO in the M-direction. Inset B shows the stable IC (madeup of red dots), corresponding to the 2D torus in B , oscil-lating around the fold on the graph of (cid:104) M (cid:48) (cid:105) . study these cells, several highly detailed models were con-structed based on experimental data ? ? . A reduced yetquantitatively-similar 5D Hodgkin-Huxley type modelwas proposed and studied in ? . Its generic representa-tion reads as follows: V (cid:48) = − I app − g K n ( V + 95) − g Na m h ( V − − g L ( V + 70) − g Ca c ( V − − g M M ( V + 95) , x (cid:48) = ( f ∞ [ V + ∆] − x ) /τ x [ V ] , (4)where the vector x represents the gating varibles: fast n , h , c and slow M , for the following four ionic currents:a delayed rectifier potassium current, a transient inacti-vating sodium current, a high-threshold non-inactivatingcalcium current and a muscarinic receptor suppressedpotassium current; here the external current I app is thebifurcation parameter of the model. The Reader is wel-come to consult with the original paper ? for more de-tailed description of the model and Appendix for its equa-tions and Matlab code.In this 5D model, the gating variable M is the slowestone. A faux parameter, ∆, is introduced in the right-hand side of the equation (see Eqs. (4)) governing M (cid:48) toperform the parameter continuation letting us sweep theslow-motion manifolds using the package MATCONT ? without affecting the fast subsystem of the model .Figure 13 A shows that the tonic-spiking manifold withthe distinct fold consists of two components: the inner2stable (green) cyllinder-shaped M SPO and the outer un-stable (red) part M
UPO . Note that M
UPO vanishes via thehomoclinic bifurcation (HB) when it touches the mid-dle, saddle branch of the quiescent manifold M eq ( S -shaped black line); see also 2D phase-space projectionsin Figs. 14 and 15).It was first reported in the original paper ? that dur-ing the transition from tonic spiking to bursting, thePurkinje cell model shows slow amplitude-modulation oftonic-spiking activity in the voltage traces (Fig. 13 C ).This modulation corresponds to a torus-canard thatemerges around the fold of the tonic spiking manifoldM PO in the phase space of the model (see Fig. 13 A and B ). This torus-canard exists in a rather narrow param-eter interval beyond which the model quickly transformsinto a square-wave buster. A. Saddle torus-canard and bistability
To examine the stability of the torus in question, wetracked the stable state of the model by varying the ap-plied current I app in small incremental steps ( ∼ × − )in both directions - increasing and decreasing. Simula-tions for each applied current value were long enoughfor the model to reach the steady state. We then repre-sented the stable sates by minima of the voltage, shownin Fig. 14 A .Interestingly, for I app ∈ [ − . , − . A ), the stable state in I app -increasing direction is tonic spiking (green dots inFig. 14 A ), while in I app -decreasing direction is burst-ing (blue dots). That is, the system is bistablewithin that small window of I app . What separates thetwo stable states is a saddle torus-canard detected atI app = − .
487 (red box), whose corresponding voltagetrace shows slow amplitude-modulation (Fig. 13 C ). Thephase-space trajectories corresponding to the two sta-ble states and the separating torus-canard are shown inFig. 14 B and the corresponding voltage traces are shownin Fig. 14 C and D . The vertical (red) line indicated theI app where MATCONT detects torus bifurcation. B. Averaging technique
As for most slow-fast systems, the location of the tonic-spiking periodic orbit on the manifold in the phase spacein the Purkinje cell model at various parameters can beaccurately predicted by the averaging technique, see ? for details. The approach is illustrated in Fig. 15 A . Itscore is the use of the parameter continuation to detectperiodic orbit(s) on M PO where the average derivative ofthe slow variable vanishes as well as to find the corre-sponding function in the right-hand side of the averageslow equation to predict forthcoming bifutcations likethe blue sky catastrophe ? . Specifically for the Purk-inje cell model, first, (cid:104) M (cid:48) (cid:105) and (cid:104) M (cid:105) are calculated foreach periodic orbit on the tonic-spiking manifold basedon the continuation data; then (cid:104) M (cid:48) (cid:105) , is plotted against (cid:104) M (cid:105) (which represents each periodic orbit on the tonic FIG. 16. 3D (Ca , n , V)-phase-space projection of theparabolic burster model. A Saddle torus (red line), enclosinga stable tonic-spiking (blue line) and a stable bursting orbit(black line) orbits are superimposed on the 1D quiescent man-ifold, M eq (dark blue line), and the 2D tonic-spiking manifoldwith stable M SPO (green surface) and unstable M
UPO (purplesurface) sections merging at the fold. Magnified inset depictsa local Poincar´e cross-section (grey plane at n = 0 .
22) with arepelling (red) IC circle, T U , enclosing the basin of the stableFP (light blue dot) corresponding to the tonic-spiking PO.Inset B shows the corresponding voltage traces of the orbitsin A ; parameters are given in Appendix. spiking manifold) to pin down where the stable tonicspiking occurs (Fig. 15 A ). The prediction from the av-eraging technique showing where (cid:104) M (cid:48) (cid:105) = 0 (dark greendot in Fig. 15 A ) is quite precise, as it matches verywell the result of the parameter continuation simulations(dark green PO in Fig. 15 A ). V. FITZHUGH-NAGUMO-RINZEL MODEL
It is obvious though that the averaging technique fora 1D slow equation cannot predict the occurrence of thetorus bifurcation. Note that the averaging approach il-lustrated in Fig. 15 A is done at the parameter valuesI app and ∆ when the saddle-torus separates two stableoscillatory states: spiking and bursting. Although theaveraging approach is only applicable to detect tonic-spiking state orbits, it can be yet helpful to locate thetorus (Fig. 15 B ). Actually, the slow averaging appliedto the torus orbit in the phase space does reveal the cor-responding invariant circle in the (cid:104) M (cid:48) (cid:105) − (cid:104) M (cid:105) projectionas shown in Fig. 15 B . VI. 4D PARABOLIC BURSTER MODEL
The same kind of bistability due to a saddle torus atthe fold, separating the stable tonic spiking and bursting,can also be observed in a simpler 4D model of parabolic3
FIG. 17. A
3D ( y, w, v )-phase space of the FNR model at δ = 0 .
08 showing 2D tonic-spiking manifold with its unstableinward (purple surface), M
UPO , merging with stable outward section (green cylinder-shape), M
SPO at the fold, as well as the 1Dquiescent manifold M eq with stable and unstable (solid/dashed) branches separated by the sub-critical AH bifurcation. Shownare trajectories of period doubling (magenta circle PD at c = − . c = − .
5) and tonicspiking (dark green circle PO at c = − . v -traces are presented in B . C
2D repelling torus (grey,T U ) containing the nested stable torus (blue, T S ) enclosing a repelling periodic orbit (red circle PO U ) at c = − . D
2D cross-section (orange plane) depicting the nested repelling, stable ICs and repelling FP (matching colors in C ). E Slowlymodulated voltage traces of the tori in C . burster. This Hodgkin-Huxley type model is given by C m V (cid:48) = − I L − I K − I Ca − I KCa − I CaS + I app ,n (cid:48) = φ ( n ∞ ( V ) − n ) /τ n ( V ) , [Ca ] (cid:48) = ε ( µI Ca − [Ca ]) ,s (cid:48) = ε ( s ∞ ( V ) − s ) /τ s , ε = 0 . , (5)(see its full equations in the Appendix below). Thismodel for an endogenously parabolic burster, i.e, withthe inter-spike frequency maximized in the middle of eachburst, was proposed and examined in ? . Its bifurcationfeature is the homoclinic saddle-node bifurcation (SNIC)in the fast 2D subsystem that should demarcate the en-try and exit points of the corresponding tonic-spikingand quiescent phases of bursting activity in the full sys-tem (5). Figure 16 shows the topology of slow-notionmanifolds in the phase-space projection of the model.Here the 2D cylinder-shaped manifold M PO originatesfrom the depolarized branch of the 1D quiescent mani-fold M eq through a subcritical Andronov-Hopf bifurca-tion. Next, the unstable (purple) section M UPO mergeswith the stable (green) section M
SPO through the saddle-node bifurcation at the first fold. The second fold on M PO takes place at the higher values of the [Ca ]-variable,where the torus bifurcation takes place in the full system.Here, the saddle torus (red) bounds the attraction basinof 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 repellingIC, representing the saddle-torus, that surrounds a lightblue fixed point corresponding to the stable tonic-spikingorbit at the right fold on the manifold M UPO .The FitzHugh-Nagumo-Rinzel (FNR) system is a for- mal mathematical neuron model of an elliptic buster ? : v (cid:48) = v − v / − w − y + I ext ,w (cid:48) = δ (0 . v − . w ) ,y (cid:48) = µ ( c − y − v ) , (6)where v and w represent the fast “voltage” and “gating”variables, while y is the slow variable due to the smallnessof µ = 0 . I ext = 0 . c in the slow equation to continue and reveal thetopology of the slow-motion manifolds of this system; theother bifurcation parameter δ is used to demonstrate howtime scales of the variables of the model can determinethe bifurcations (torus and period-doubling) of periodicorbits of this model. Figure 17 A shows the tonic spik-ing manifold of the FNR model consisting of an unstableinner layer M UPO (purple surface) emerging though a sub-critical Andronov-Hopf bifurcation from the 1D quiescentmanifold, M eq , and a stable outer layer M SPO (green sur-face), both merging though a fold.The ( c, δ )-bifurcation diagram, shown in Fig. 18 A , ofthe FNR model is obtained using the parameter contin-uation package MATCONT ? . It includes several bifur-cation curves: AH stands for the Andronov-Hopf [super-critical] bifurcation of the only equilibrium state of theFNR model: TB stands for the torus bifurcation; alongthis curve the multipliers, e ± iφ , of the bifurcating pe-riodic orbit runs through strong resonances detected atpoints labeled as R 1:4 with φ = π/
4, and R 1:3 with φ = 2 π/
3, through the final resonance R 1:2 where φ completes the arch of the unit circle and reaches π . Thelast point is included in the PD-curve corresponding tothe period doubling bifurcation through which the peri-odic orbit looses its stability. The diagram also includestwo bifurcation curves, SN, of saddle-node periodic or-bits. These curves follow the PD-curve closely. They4cannot be further continued in the c -parameter any fur-ther down the bifurcation diagram for the given slopeof the slow nullcline v = c − y (being a 2D plane inthe 3D phase space) because it becomes tangent to thefolds of the 2D spiking manifold M PO . We hypothesizethat the bifurcation curves, SN’s, and PD, will terminateat the corresponding resonances, 1:1 and 1:2 occurringnear the point, labelled as a red dot R (with δ (cid:39) c, δ )-parameter plane. Geometrically, this pointcorresponds to a cusp underlying the transition fromthe hysteresis to the monotone, increasing dependenceof the Vmax / min- and (cid:104) V (cid:105) -coordinates of the periodicorbits on the c -parameter in the bifurcation diagrams inFigs. 18( B ) and ( C ), respectively.Observe the ∩ -shape of the TB-curve that indicatesthat depending on the level of the δ -parameter, the pe-riodic orbit can undergo two supercritical torus bifurca-tions, forward and backward, when its looses stabilityto the torus and re-gains for δ = 0 . c isincreased, for δ = 0 .
08, (lower dashed grey line).The bifurcation diagrams in Figs. 18 B - C , representingthe maximum, v max , and minimum values, v min , of thevoltage plotted against the c -parameter, illustrate howthe tonic-spiking manifold is shaped and where the peri-odic orbits that constitute it, bifurcate (indicated by ver-tical red lines) as c -parameter is varied for fixed valuesof δ = 0 .
08 and δ = 0 .
3, respectively. One can clearlysee that at a rather large value δ = 0 . PO unlikethe period-doubling bifurcation that occurs between twosaddle-node bifurcations (indicated by the vertical lineslabeled as SNs in Fig. 18 B ) corresponding to the twofolds on M PO at δ = 0 . c -parameter sweeping simulations of the limiting solu-tions of the model when δ = 0 .
08 and when δ = 0 . D and E .At the level δ = 0 .
3, after crossing AH-line in the bifur-cation diagram Fig. 18 A , trajectories of the FNR-modelare attracted into a whirlpool stirred by a stable PO inthe phase space shown in Fig. 18 F . After crossing theleft branch of TB-curve, this stable PO loses stability (redcircle) to a stable [ergodic and round] torus in the phasespace depicted in Fig. 18 F . The sweeping bifurcation di-agram in Fig. 18 E shows that with a further increase ofc, the torus becomes resonant (Fig. 18 F ) and its shapebecomes more complex giving rise to the so-called mixed-model oscillations (MMO) of alternating small-amplitudesub-threshold and large spiking ones. Once crossing theright branch of TB line, the torus goes into a strong res-onance 1:3 before it breaks down at larger values of the c -parameter giving rise to the stable period-3 orbit in thephase space, see Fig. 18 F .At δ = 0 .
08, after crossing the left branch of TB-curvein the bifurcation diagram in Fig. 18 A , the PO loses sta-bility (red line in Fig. 17 C ) to a stable torus (blue line)in the phase space. Interestingly, at c = − . D ). Note that in a 3D sys-tem a repelling torus can be detected and traced down inthe backward time when it becomes attracting, which isimpossible for saddle tori occurring in high-dimensionalsystems.The same approach employed to detect bistability inthe Purkinje cell model due to a hysteresis effect, is usedto visualize bistability in the FNR-model as its solutionsare swept in the opposite directions of the c -parameter.The obtained bifurcation diagram, with v min -values (asgreen/blue dots) plotted against increasing/decreasing c -values, is represented in Fig. 19. Here a single dot for thegiven c -value corresponds to a stable periodic orbit of asingle v min -value, while several vertical dots correspondto co-existing busting orbits. The basins of these orbitsare bounded by the repelling torus enclosing the stabletonic spiking orbit in the 3D phase space of the FNR-model.Again we note that the tori emerge in the FNR modelnear the smooth cone-shaped end of M UPO instead ofaround its fold (Fig. 17 A inset). Fig. 17 E shows theslowly-modulated “voltage” traces corresponding to thetwo tori. As c -parameter further increases to c = − . A and B ). The period-doublingcascade starts at the left branch of the PD-curve andends at its right branch in the bifurcation diagram inFig. 18 A and Fig. 18 D (magenta line in Fig. 17 A and B , c = − . A and B , c = − . A. Origin of slow-fast torus bifurcation
Let us outline a possible origin and stability of thetorus bifurcation in the FNR-model. Its start with asaddle-node [fold] bifurcation of two limit cycles, sta-ble 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. So, thevalue y = 0 . S and unstable (red circle) LC U limitcycles emerge through a saddle-node [fold] bifurcationin the ( v, w )-subspace, see Fig. 20 A . While the stablelimit cycle LC S increases in its size to become the large-amplitude orbit of the relaxation oscillator, the unstablelimit cycle LC U keeps shrinking as y increases, see thecorresponding ( v, w )-phase portrait in Figs. 20 B and C .As the result, in the full model, the slow increase/de-crease of the y -variable makes the LC U drag back andforth to form the 2D repelling torus in the 3D phasespace of the model, see Fig. 20 D .5 FIG. 18. A Bifurcation diagram of the FNR model depicts the curves for the following bifurcations: Andronov-Hopf (AH,black line), period-doubling (PD, green line), saddle-node (SN, purple line) and torus bifurcation (TB, blue line) with indicatedstrong resonances 1:4, 1:3 and 1:2 (red dots). Dashed grey lines: c-parameter pathways at δ = 0 .
08 and δ = 0 .
3, presented in B – C and in D – F . B – C Parameter continuations of v max , v min and (cid:104) v (cid:105) plotted against the c-parameter at δ = 0 .
08 and δ = 0 . D – E The c-parameter sweeping diagrams for v max at δ = 0 . δ = 0 .
3, resp., reveal: D PD-cascade en route to elliptic busting, and E the torus emergence, evolution and breakdownthrough 1:3-resonance. F − F
2D (y , v)-phase projections depicting a stable PO (green) with a whirlpool born after AH at c = − . c = − . c = − . c = − .
47; light-blue highlighted are intersection points of the trajectories with a local transverse plane.FIG. 19. Sweeping diagram of the attractors in the FNRmodel, represented by their V min -values (green/blue dots)plotted against the c -parameter being increased/decreased.Within the bistability window (horizontal black bar), tonic-spiking and bursting attractors co-exist and are separated bythe repelling torus. Vertical (red) line at c = − . B. Divergence-sign test to predict TB or PD
Either the torus bifurcation (TB) or the period-doubling bifurcation (PD) can underly the stability loss of the periodic orbit existing on the attracting cylinder-shape section of the tonic-spiking manifold M PO in thephase space of the FNR-model when it reaches the fold,see Fig. 17. We have already noticed depending on thevalue of the δ -parameter, this large amplitude orbit lossesstability though the torus bifurcation followed by its res-onances and breakdown for δ ≥ . δ = 0 .
08. Thisis indicative that wherever the time scales of the v - and w -variables becomes of a similar order the torus bifurcationcase prevails over the period doubling one where both w, y -variables become much slower than the v -variable,and vice versa.To predict the type of the bifurcation, torus or perioddoubling, the stable orbit will undergo, without comput-ing its multipliers, we can evaluate the sign of the di-vergence given by ∇ = (cid:104) ∂v (cid:48) ( t ) ∂v + ∂w (cid:48) ( t ) ∂w + ∂y (cid:48) ( t ) ∂y (cid:105) alongthe periodic orbit on the fold of the slow-motion mani-fold in the 3D phase space of the FNR model. The ar-guments are similar to the well-known Bendixson-Dulaccriterion for the existence of a close orbit, like a limitcycle and a separatrix loop of a saddle in a system ona 2D plane; note the Dulac’s argument follows from the6 FIG. 20. A
2D ( v, w )-projection depicting the stable (green) M
SPO , and unstable (purple) M
UPO branches of the tonic spiking-spiking manifold of the FNR system at δ = 0 .
08 that are superimposed A with the saddle-node limit cycle (which then bifurcatesinto a stable (blue) limit cycle LC S and an unstable (red) limit cycle LC U ) of the fast subsystem at y = 0 . B withlimit-cycle canards following the stable and unstable branches of the cubic fast nullclines at y = 0 . C
3D phase spacedepicting M
SPO and M
UPO being superimposed with stable and unstable LC-canards from A – B . Inset D zooms-in the torus(T U , grey), the unstable limit cycle LC U (red ring) around y = 0 . ∇ = (cid:104) ∂v (cid:48) ∂v + ∂w (cid:48) ∂w + ∂y (cid:48) ∂y (cid:105) along the peri-odic orbit (with period normalized to 1) superimposed withthe blue level-curve for the average divergence (cid:104)∇(cid:105) (its valueshown below) for the PD-bifurcation scenario ( A ) or for theTB scenario ( B ) at δ = 0 . Green’s theorem ? . Loosely speaking, it can be directlyiterated as follows: no closed orbit exists in a domain ofthe 2D phase plane where the vector field divergence isof a constant sign. As far as the 2D torus in a 3D phasespace is concerned, it can not emerge in a 3D subspacewhere the divergence does not change its sign. For exam-ple, this guarantees no torus bifurcation in the classicalLorenz-63 model with the chaotic attractor, which is adissipative system with a constant negative divergenceto collapse phase volumes. On the other hand, the torusbifurcation occurs in the 3D Lorenz-84 model ? . One ofthe features of the Lorenz-84 model is that the diver-gence of its vector filed is phase-coordinate dependentand hence can be sign-alternating. This is why the torusbifurcation is a part of bifurcation unfolding of an equilib-rium state with the characteristic exponents (0 , ± iω ) (thesum of the exponents vanishes), which is also referred toas the codimension-two Gavrilov-Guckenheimer or sim-pler fold-Hopf bifurcation, which is the key feature of theLorenz-84 model. Such a bifurcation may also occur inthe FNR-model near the tip of the cone-shaped manifoldM PO crossed out by the slow nullcline y (cid:48) = 0 with µ (cid:28) (cid:104)∇(cid:105) , i.e., all ∇ -valuessummed up along the periodic orbit and averaged overits period, on the fold of the manifold M PO , is close to zero then the stable periodic orbit may undergo the torusbifurcation at the stability loss. Otherwise if (cid:104)∇(cid:105) < (cid:104)∇(cid:105) in the FNR model (6)at various points on the PD- and TB-curves supports ourhypothesis. Figure 21 shows how ∇ varies over the nor-malized period of the orbit loosing the stability thoughthe period-doubling and torus bifurcations. Its average (cid:104)∇(cid:105) -value is about ∼ .
56 and 0 in the period-doublingand torus cases, respectively, when δ = 0 .
08. This de-facto validates the divergence-based approach and thatit can be used to predict the type of the bifurcations willoccur at the loss of stability of periodic orbits on the dis-tinct fold of the slow-motion manifold in the phase spaceof 3D slow-fast systems.
VII. NORMAL FORMS FOR TORUS-CANARDS
The amplitude modulation phenomenon that happensin various complex and realistic neuron models inspires usto build a simple and adjustable mathematical model togenerate and analyze torus bifurcations using evident ge-ometric constrains. Since such tori emerge near the foldof the tonic-spiking manifold, the first step is to createthe desired number of such folds in the 2D normal-formfast subsystem given by x (cid:48) = ωy + xQ ( x, y, ε ) ,y (cid:48) = − ωx + yQ ( x, y, ε ) (7)with a polynomial function Q = ε + (cid:80) ∞ n =1 l n ( x + y ) n .When Q ≡
0, Eqs. (7) describe a linear harmonic oscil-lator that has a single equilibrium state – the center atthe origin in the ( x, y )-phase plane. When Q = ε < ε >
0. When Q ( x, y, ε ) is nega-tively/positively defined, the nonlinear system has yet asingle equilibrium state that globally attracts/repels allother spiral-shaped trajectories. A level curve given by Q = 0 correspond to a closed orbit – limit cycle of thesystem. In case of a single zero (other then the origin),7 FIG. 22. The ( z, y, x )-phase space of the the 2-layer model (8) depicting the tonic-spiking manifold with an outward stablepart M
SPO (green surface) merging with the inward unstable section M
UPO (red surface) at the fold, and the quiescent manifoldM eq (black line) with stable (solid) and unstable sections on the z -axis. the slow nullcline z (cid:48) = 0 is represented by the greycylinder as its radius at selected ∆-values: -0.01, -0.3803, -0.9. An outside trajectory (blue line) in A quickly spirals ontoM eq , then gets “sucked” into the cylinder (where z (cid:48) > SPO outside the cylinder where ( z (cid:48) < B transitions into an elliptic burster in C . Other parameters: ω = 2 , l = 1 , l = − .
8. Insets A – C show the corresponding “voltage” traces.FIG. 23. The ( z, y, x )-phase space of the 3-layer [relaxation torus] model (9) showing the tonic-spiking manifolds with twofolds connecting its inner and outer stable (green surfaces, M SPO ) and the middle unstable (red surface, M
UPO ) branches. Thequiescent manifold, the z -axis (black line) has two sections: stable ( z <
0, solid) and unstable ( z >
0, dashed). The slownullcline z (cid:48) = 0 is a pair of parallel (grey) planes between which z (cid:48) >
0, and z (cid:48) < A – C Rearrangements of M PO with l -variations at indicated values reshape modulations of voltage traces shown in A – C at ∆ = 1.5, 3.5 and 1.092; otherparameters: l = − . , l = − . , ω = 10 , µ = 0 . the limit cycle is double or of the saddle-node type. Incase of two distinct zeros, they correspond to two nestedlimit cycles surrounding the origin that determines theirstability. So, if the origin is stable for ε <
0, then theinner LC is repelling, whereas the outer one is stable, orvice versa if ε >
0. In the special case ε = 0, the sta-bility of the limit cycle is determined by the sign of thecorresponding Lyapunov coefficient l n .This basically describes the algorithm for devising thedesired systems. Let us first consider the model with asingle fold. Its equations are given below: x (cid:48) = ωy + x (cid:16) z + l ( x + y ) + l ( x + y ) (cid:17) ,y (cid:48) = − ωx + y (cid:16) z + l ( x + y ) + l ( x + y ) (cid:17) ,z (cid:48) = µ (cid:0) − x − y + ∆ (cid:1) , (8)with ( x, y ) being the fast variable and z being a slow (as µ is small), dynamic variable replacing fixed ε . When µ = 0, the first two equations present a normal formfor the Bautin bifurcation of the equilibrium state at theorigin in case when ε = 0 and first Lyapunov coefficient, l = 0, of the cubic terms (note that l can be scaled toequal −
1, for simplicity). The last constrains implies thatthe outer LC will be stable and a nested inner one will berepelling as l >
0. The equilibrium state at the origin isstable as long as z <
0, and becomes unstable when z ≥ z = 0. This ends up bi-stabilityin the fast subsystem system and the stable LC is theonly global attractor. The left boundary for bistability isgiven by z sn = − (cid:104) l ( x + y ) + l ( x + y ) (cid:105) that corre-spond 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 phase8 FIG. 24. Torus-canard transformations in the 3-layer [relaxation-torus] model (9) as ∆-parameter is decreased at indicatedvalues. 3D ( z, y, x )-phase space depicts the tonic-spiking manifold M PO (grey surface) with two folds connecting its stableinner, outer and middle unstable sections. The quiescent manifold is the z -axis (black line) with stable ( z <
0, solid) andunstable ( z >
0, dashed) sections. The slow nullcline z (cid:48) = 0 is a pair of parallel (yellow) planes between which z (cid:48) > A , – E , Stable PO (green circle) loses stability to a torus-canard (blue IC on a Poincar´e cross-section given by y = 0) near the leftfold that becomes wider and weekly resonant. F , – H , Unstable PO (red cycle/FP) regains stability with a sub-critical TBgenerating a small repelling torus (red IC, T U ) that grows in size, merges with the stable torus-canard T S and both annihilate,leaving the stable PO on the inner branch of M PO . Other parameters: l = − . , l = 1 . , l = − . , ω = 10 , µ = 0 . . space.The slow equation in (8) is designed here as follows:the slow nullcline, where z (cid:48) = 0 is a cylinder given by x + y = 1+∆ of radius √ a + ∆ in the 3D phase space ofthe system. Inside this cylinder z (cid:48) >
0, and z (cid:48) < PO inMATCONT.Now, let us put all components together in the fullsystem. Figure 22 shows the 3D phase space of the2-layer model (8). Here is the 1D quiescent manifoldM eq (the z -axis) when stable ( z <
0) and unstable( z >
0) branches. The subcritical AH bifurcation gives rise to the unstable (red) branch M
UPO that folds backand continues further (rightwards) as stable (green)branch M
SPO . The slow nullcline is represented as agrey cylinder parallel to the z -axis. Let us pick aninitial point to the left quite far away from the foldon M UPO in the 3D phase space. Then, the trajectoryfast spirals onto the stable branch (solid line) of M eq inside the slow nullcline/cylinder along which it will beslowly dragged to the right toward the AH bifurcationin the fast subsystem. Passing throughout it, the phasepoint keeps following the unstable branch of M eq (theso-called delayed loss of stability) before it spirals awayto converge to the stable, outer branch M SPO . Provided9that radius of M
SPO is larger than that of the slow null-cline, the phase point, turning around M
SPO , will slideleftward toward the fold on the 2D spiking manifold. If∆ controlling the radius of the slow nullcline is largeenough, then the phase point will converge to a stableand round periodic orbit on M
SPO . Varying ∆ makes thisperiod orbit slide along M PO thus revealing its shape,including the fold. This is the essence of the parametercontinuation approach to localize slow-motion manifoldsof the full system by varying the control [faux] parameterof the slow subsystem in the phase space. By decreasing∆, the stable PO slides down towards the fold belowwhich the orbit continues as repelling both unstablein z -direction and in the fast ( x, y )-coordinates. Aswe pointed earlier, the transverse intersection of M PO with the slow nullcline/cylinder guarantees that the POdoes not vanish but loses stability at the fold through asuper-critical torus bifurcation that result in the emer-gence of the stable torus-canard (provided µ is smallenough) switching back and forth between the stableand unstable branches of the tonic spiking manifoldnear its fold, see Figure 22 B . Further decreasing ∆makes the stable torus morph into the “elliptic” busterwith amplitude modulations and episodes of quiescentmeta-states.Our last 3-layer (torus-relaxation) model below x (cid:48) = ωy + x (cid:16) z + l ( x + y )+ l ( x + y ) + l ( x + y ) (cid:17) ,y (cid:48) = − ωx + y (cid:16) z + l ( x + y )+ l ( x + y ) + l ( x + y ) (cid:17) ,z (cid:48) = µ (1 − x + ∆) (9)is designed to have two folds between two stable, inner( l <
0) and outer ( l <
0) and middle unstable ( l > PO , see itsphase space in Fig. 23. The distinction of the slow equa-tion in Eqs. (9) is that its slow nullcline is represented bytwo parallel planes x = ±√ z -axis, functionally provide the sameslow dynamics and means of its control by varying ∆ inthe full system.Like in the 2-layer model, by varying ∆ one can setthe 3-layer model to demonstrate various tonic spikingoscillations and torus-canards. In addition, by calibrat-ing values of l and l one can shifts the relative po-sitions of the folds to shape and produce various typesof elliptic bursters generated by this system, as Fig. 23demonstrates. Let us describe the torus-canard transfor-mations as ∆ is decreased. The principal stages of thetransformation are documented in Fig. 24.At the initial stage the stable PO exists on the outersurface M SPO at some large enough ∆, see Fig.24 A . Withan initial decrease, this stable PO slides down off the leftfold and loses stability to a torus-canard, T S (Fig.24 B ).With further decreasing ∆ the repelling PO slides downalong the middle branch M UPO to the right fold, mean-while the stable relaxation torus-canard with a head grows in size expanding the whole space between the in-ner and the outer layer, its stages are shown in Figs. 24 C - E . Next, via a subcritical TB, a repelling torus, T U ,emerges at the right fold from the PO that regains thestability, see Figs. 24 F - G . At this point there are twotori: stable (blue) and unstable (red IC in the cross-section). Note that the repelling torus separates the at-traction basins of the stable torus from the stable PO,thereby creating bistability in the system. Finally, thestable and the unstable tori merge and annihilate througha saddle-node bifurcation so that the stable PO at theright fold remains the only attractor in the phase space,see Fig. 24 H . VIII. CONCLUSION
Following the bottom-up approach we have consid-ered and analyzed a collection of chosen exemplary slow-fast models, starting off with the biologically plausibleHodgkin-Huxley ones up to light mathematical toy sys-tems that all feature the torus bifurcations occurring onthe characteristic fold on the slow-motion tonic-spikingmanifold. Unlike the flat canards occurring in 2D slow-fast systems whose stability can be evaluated analyticallyas it is dictated by the analytical properties of the func-tion on the right-hand side of the fast equation in thesingular limit, the question concerning stability of emer-gent tori is yet to be fully understood. We have shown allkinds of tori in the model list, from stable and repelling(made stable in the backward time in 3D systems) tosaddle ones with 3D unstable manifolds and XD stablemanifold, where X is determined by the phase space di-mension of the model in question. The fact that the torusemerges locally next to the fold of a low-dimensional sur-face let one use one’s skills to choose a suitable Poincar´ecross-section and to find two initial conditions on it todemonstrate that one trajectory goes inside such a torusto converge to a nested stable periodic, while the otherconverges to another attractor in the phase space. Thisbistability is a de-facto proof of the existence of 2D sad-dle tori in the phase space, and can be further supportedby bi-directional parameter sweeps to reveal the hystere-sis due to overlapping interval of the co-existence of twoattractors, like tonic-spiking and bursting orbits in thePurkinje cell model and the FNR-model.While parameter continuation packages like MAT-CONT and AUTO are very handy to detect torus bifur-cations and some strong resonances on the correspond-ing bifurcation curves, one has to use additional compu-tational tools to analyze torus bifurcations and break-downs to demonstrate that they follow up with the ex-isting mathematical theory. One such reduction tool, es-pecially in the context of neuroscience and neurophysi-ological models, is to construct reduced Poincar´e mapsfrom time series. This can be done using the map: V min n → V min n +1 for minimal or maximal values found inof voltage traces. Using this approach we have shownthat the emergence, evolution of invariant circles corre-sponding to tori in the phase space and slow modulationsin voltage traces in our simulation match very well the0 FIG. 25. Poincar´e return maps, τ n → τ n +1 , for recurrent timeintervals between consecutive spikes in the voltage traces. A Unstable (red) IC [corresponding to the saddle torus-canardT U ] encloses a stable (blue) FP (corresponding to tonic-spiking orbit PO S ) in Purkinje cell model at I app = − . B Stable (blue) IC (T S ) bounding the unstable (red) FP(PO U ) in the normal-form (8). Inset C shows a quasi-periodicvoltage trace with recurrent time intervals τ n employed in thereturn map in A for Purkinje cell model. non-local theory proposed in the classical works by V.S.Afraimovich and L.P. Shilnikov ? ? ? ? ? . Namely, suchan invariant circle, born smooth and round, first becomesnon-smooth, resonant and next gets distorted by devel-oping homoclinic tangles of unstable sets of resonant sad-dle periodic orbits on the torus. The latter leads to thetorus breakdown and to the onset of emergent hyperbolicdynamics. This prerequisite for deterministic chaos com-petes with regular dynamics due to the presence of stableperiodic orbits overshadowing its host, i.e. the resonanttorus. Weakly resonant tori exist in narrow ranges in theparameter space of the system. The same holds true forergodic tori as they quickly become non-smooth betweenresonant zones densely populating the parameter spaceof such systems.The other option for the visualization and analysis oftori is to use the return map, τ n → τ n +1 , for recurrenttime intervals between consecutive spikes in the voltagetraces. This approach is illustrated in Fig. 25. The mapin Fig. 25 A shows the saddle IC (red, T U ) enclosing a sta- ble FP (blue, PO S ) corresponding to robust and periodictonic spiking in the Purkinje cell model. The return mapin Fig. 25 B depicts the stable IC (blue, T S ) encircling anrepelling tonic-spiking FP (red, PO U ) in the 3F normal-form model (8). This demonstrates that both amplitude-and frequency-modulation approaches suite well to studyquasi-periodicity using the observables like recordings ofvoltage traces.We note that transitions from tonic-spiking to burst-ing in the given family of slow-fast models can also beunderlined by a cascade of period-doubling bifurcations,instead of the torus bifurcation. We proposed and testedthe divergence-sign assessment to predict the type of bi-furcation. The test works well for 3D model, while itsapplicability is problematic for higher dimensional sys-tems, as one has to single out a 3D subspace in restric-tion to which the analogous divergence should be eval-uated. It has become more or less clear that the choicebetween these two bifurcations is implicitly determinedby multiple or just two time-scale reciprocal interactionsof dynamic variables of the slow-fast model, as well asby the smoothness of the fold in the averaged system.The smoother such a fold and the closer the system tothe singular limit, the more likely it undergoes the torusbifurcation yet unclear weather it will be sub or super-critical. IX. ACKNOWLEDGMENTS
This work was partially funded by NSF grant IOS-1455527, RSF grant 14-41-00044 at the Lobachevsky Uni-versity of Nizhny Novgorod, RFFI grant 11-01-00001,and MESRF project 14.740.11.0919. H.J. and A.S. thankGSU Brain and Behaviors Initiative for the fellowshipand pilot grant support. We are grateful to H. Meyerfor his patient MATCONT guidance, and to all mem-bers of Shilnikov’s NeurDS lab for the helpful discus-sions and suggestions. A.N. acknowledges support fromthe Quantitative Biology Institute and from NeuroscienceProgram at Ohio University.
X. APPENDIX
The Purkinje cell equations are as follows:1 V (cid:48) = − I app − . n ( V + 95 . − . (cid:0) e − . V +34 . (cid:1) h ( V − . − . V + 70 . − . c ( V − . − . M ( V + 95 . h (cid:48) = 1 (cid:16) e V +59 . . (cid:17) (cid:18) .
15 + . e V +33 . . (cid:19) (1 − h ) − (cid:18) −
11 + e V +59 . . (cid:19) .
15 + . e V +33 . . hn (cid:48) = 1 (cid:0) e − . V +29 . (cid:1) (cid:0) .
25 + 4 . e − . | V +10 . | (cid:1) (1 − n ) (10) − (cid:32) − (cid:0) e − . V +29 . (cid:1) (cid:33) .
25 + 4 . e − . | V +10 . | nc (cid:48) = 1 . . − c )1 + e − . V − . − . c ( V + 8 . − e . V +8 . M (cid:48) = 0 .
021 + e − . V +20 . (1 − M ) − . e − ( V +43 . . M . Matlab codes
1. 12D hair cell model is given bygL = 0 . 1 7 5 ;gk1 = 2 9 . 5 ;gh = 2 . 2 ;gCa = 1 . 2 ;BK= 0 . 0 1 ;Ek= − − − − − − − − − − − − − − − zeros ( 1 2 , 1 ) ; m k1 s i nf =1.0/(1.0+ exp ( ( y ( 1 ) +0.110) / 0 . 0 1 1 )) ;t k 1 f =0.7 e − ∗ exp ( ( y ( 1 ) +0.120) /( − − − ∗ exp ( ( y ( 1 ) +0.120) /( − − ∗ − ∗ (y ( 1 ) − Ek ) ∗ ( 0 . 7 ∗ y ( 2 ) +0.3 ∗ y( 3 ) ) ;Ih=gh ∗ − ∗ (y ( 1 ) − Eh) ∗ ( 3 ∗ y ( 4 ) ∗ y ( 4 ) ∗ (1.0 − y( 4 ) )+y ( 4 ) ∗ y ( 4 ) ∗ y ( 4 ) ) ;mhinf =1.0/(1+ exp ( ( y ( 1 ) +0.087+ I s c a n )/ 0 . 0 1 6 7 ) ) ;th =0.0637+0.1357 ∗ exp ( − (( y ( 1 ) +0.0914)/ 0 . 0 2 1 2 ) ∗ ( ( y ( 1 ) +0.0914) / 0 . 0 2 1 2 ) ) ;ex= exp ( − c o n s t 2 ∗ y ( 1 ) ) ;permiab=c o n s t 1 ∗ y ( 1 ) ∗ ( Ki − Ko ∗ ex ) /(1.0 − ex ) ;mdrkinf =1.0/ sqrt (1.0+ exp ( − (y ( 1 ) +0.0483)/ 0 . 0 0 4 1 9 ) ) ;a l f d r k = 1 . 0 / ( 0 . 0 0 3 2 ∗ exp ( − y ( 1 ) / 0 . 0 2 0 9 )+0.003) ;b e t d r k = 1 . 0 / ( 1 . 4 6 7 ∗ exp ( y ( 1 ) / 0 . 0 0 5 9 6 )+0.009) ;taudrk =1.0/( a l f d r k+b e t d r k ) ;I d r k=y ( 5 ) ∗ y ( 5 ) ∗ Pdrk ∗ permiab ;mCainf =1.0/(1.0+ exp ( − (y ( 1 ) +0.055)/ 0 . 0 1 2 2 ) ) ;tauCa =0.046 e − − ∗ exp ( − (( y ( 1 )+0.077) / 0 . 0 5 1 6 7 ) ∗ ( ( y ( 1 ) +0.077)/ 0 . 0 5 1 6 7 ) ) ;I c a=gCa ∗ − ∗ (y ( 1 ) − Eca ) ∗ y ( 6 ) ∗ y ( 6 ) ∗ y ( 6 ) ;ex1= exp ( d d e l 1 ∗ z ∗ c o n s t 2 ∗ y ( 1 ) ) ;k1=km1/ k01 ∗ ex1 ;k2=km2/ k02 ;k3=km3/ k03 ∗ ex1 ;a l p h c=a l p h c 0 ∗ exp ( y ( 1 ) /VA) ;yy7 =1.0 − (y ( 7 )+y ( 8 )+y ( 9 )+y ( 1 0 ) ) ;h b k t i n f =1.0/(1.0+ exp ( ( y ( 1 ) +61.6 e −
3) / 3 . 6 5e −
3) ) ;taubkt =2.1 e − − ∗ exp ( − (( y ( 1 ) +66.9 e −
3) / 1 7 . 7 e − ∗ ( y ( 1 ) +66.9 e −
3) / 1 7 . 7 e −
3) ;2IBKS=PBKS ∗ BK ∗ permiab ∗ ( y ( 9 )+y ( 1 0 ) ) ;IBKT=PBKT ∗ BK ∗ permiab ∗ ( y ( 9 )+y ( 1 0 ) ) ∗ y ( 1 2 ) ;IL=gL ∗ − ∗ y ( 1 ) ;dy ( 1 ) = − (I k 1+Ih+I d r k+I c a+IBKS+IBKT+IL ) /Cm;dy ( 2 ) =( mk1sinf − y ( 2 ) ) / t k 1 f ;dy ( 3 ) =( mk1sinf − y ( 3 ) ) / t k 1 s ;dy ( 4 ) =(mhinf − y ( 4 ) ) / th ;dy ( 5 ) =(mdrkinf − y ( 5 ) ) / taudrk ;dy ( 6 ) =(mCainf − y ( 6 ) ) / tauCa ;dy ( 7 )=k1 ∗ y ( 1 1 ) ∗ yy7+km2 ∗ y ( 8 ) − (km1+k2 ∗ y( 1 1 ) ) ∗ y ( 7 ) ;dy ( 8 )=k2 ∗ y ( 1 1 ) ∗ y ( 7 )+a l p h c ∗ y ( 9 ) − (km2+b e t c) ∗ y ( 8 ) ;dy ( 9 )=b e t c ∗ y ( 8 )+km3 ∗ y ( 1 0 ) − ( a l p h c+k3 ∗ y( 1 1 ) ) ∗ y ( 9 ) ;dy ( 1 0 )=k3 ∗ y ( 1 1 ) ∗ y ( 9 ) − km3 ∗ y ( 1 0 ) ;dy ( 1 1 )= − U ∗ I c a / ( z ∗ f a r a d ∗ Cvol ∗ e e p s ) − Ks ∗ y( 1 1 ) ;dy ( 1 2 ) =( h b k t i n f − y ( 1 2 ) ) / taubkt ;2. 5D Purkinje Cell model is given byV’= − Iapp − (95.0+V) ∗ ( 0 . 7 5 ∗ M) − ∗ cˆ2 ∗ ( − − ∗ (70.0+V) − ∗ (1.0/(1.0+ exp (( − V − ∗ h ∗ ( − − ∗ n ˆ 4 . 0 ∗ ( 9 5 . 0 +V) ;h ’ = 1 . 0 / ( 1 . 0 + exp ( (V+59.4) / 1 0 . 7 ) )/ ( 0 . 1 5 + 1 . 1 5 / ( 1 . 0 + exp ( (V+33.5) / 1 5 . 0 ) ) ) ∗ (1.0 − h ) − (1.0 − exp ( (V+59.4) / 1 0 . 7 ) ) ) / ( 0 . 1 5 + 1 . 1 5 / ( 1 . 0 + exp ( (V+33.5)/ 1 5 . 0 ) ) ) ∗ h ;n ’ = 1 . 0 / ( 1 . 0 + exp (( − V − ∗ exp ( − abs (V+10.0) / 1 0 . 0 ) ) ∗ (1.0 − n ) − (1.0 − exp (( − V − ∗ exp ( − abs (V+10.0)/ 1 0 . 0 ) ) ∗ n ;c ’ = ( 1 . 6 ∗ ( 1 . 0 − c ) ) /(1.0+ exp ( − ∗ ( − − (0.02 ∗ c ∗ (8.9+V) ) /( − exp ( ( 8 . 9 +V) / 5 . 0 ) ) ;M’ = 0 . 0 2 / ( 1 . 0 + exp (( − − s h i f t − V) / 5 . 0 ) ) ∗ (1.0 − M) − ∗ exp (( − − V) / 1 8 . 0 ) ∗ M;3. 4D parabolic burster model is given byV1’=( − gL ∗ (V1 − EL) − (gK ∗ n1+gKCa ∗ ( Ca1/(1+Ca1) ) ) ∗ (V1 − EK) − (gCa ∗ ( 0 . 5 ∗ ( 1 + tanh ( ( V1+1.2) / 1 8 ) ) )+gCaS ∗ x1 ) ∗ (V1 − ECa)+Iapp ) /Cm;n1 ’= p h i ∗ cosh ( ( V1 −
12) / 3 4 . 8 ) ∗ ( ( 0 . 5 ∗ ( 1 + tanh ( ( V1 −
12) / 1 7 . 4 ) ) ) − n1 ) ;x1 ’= eps ∗ ( ( 0 . 5 ∗ ( 1 + tanh ( ( V1 −
12) / 2 4 ) ) ) − x1 )/ 0 . 0 5 ;Ca1’= eps ∗ ( − mu ∗ ( ( gCa ∗ ( 0 . 5 ∗ ( 1 + tanh ( ( V1+s h i f t +1.2) / 1 8 ) ) )+gCaS ∗ x1 ) ∗ (V1+C a s h i f t − ECa) ) − Ca1 ) ;p a r a m e t e r s : gL=2, gK=8, gKCa=1; gCa=4,gCaS=1, EL= −
60, EK= −
84, ECa=120 , p h i=0.06666 , epseps