Hysteresis bifurcation and application to delayed Fitzhugh-Nagumo neural systems
HHysteresis bifurcation and application to delayedFitzhugh-Nagumo neural systems (cid:73)
L. Chen, S. A. Campbell ∗ Department of Applied Mathematics, University of Waterloo, Waterloo, ON, N2L 3G1,Canada
Abstract
Hysteresis dynamics has been described in a vast number of biological experi-mental studies. Many such studies are phenomenological and a mathematicalappreciation has not attracted enough attention. In the paper, we explore thenature of hysteresis and study it from the dynamical system point of view byusing the bifurcation and perturbation theories. We firstly make a classificationof hysteresis according to the system behaviours transiting between differenttypes of attractors. Then, we focus on a mathematically amenable situationwhere hysteretic movements between the equilibrium point and the limit cycleare initiated by a subcritical Hopf bifurcation and a saddle-node bifurcation oflimit cycles. We present a analytical framework by using the method of multipletimescales to obtain the normal form up to the fifth order. Theoretical resultsare compared with time domain simulations and numerical continuation, show-ing good agreement. Although we consider the time-delayed FitzHugh-Nagumoneural system in the paper, the generalization should be clear to other systemsor parameters. The general framework we present in the paper can be naturallyextended to the notion of bursting activity in neuroscience where hysteresis isa dominant mechanism to generate bursting oscillations.
Keywords: hysteresis, bifurcation, Fitzhugh-Nagumo neuron, time delayed, (cid:73)
This work was supported by the Natural Sciences and Engineering Research Council ofCanada. ∗ Corresponding author
Email addresses:
[email protected] (L. Chen), [email protected] (S. A.Campbell)
Preprint submitted to Elsevier September 30, 2020 a r X i v : . [ n li n . AO ] S e p ursting
1. Introduction
Hysteresis widely exists in biology from microscopic cell biology [1] , ge-netics [2] and neuroscience [3] up to macroscopic bio-mechanical properties oforgans such as the eye [4] and muscle [5]. In particular, hysteresis is one es-sential mechanism to generate bursting oscillations which play important rolesin communication between neurons [6]. However, there had been few mathe-matical investigations of this biological process until the discovery of a numberof molecular mechanisms with bistable dynamical behavior by the early 1990s[7]. In addition, influenced by mathematical treatments to physical and en-gineering systems, most studies in biology concentrate on identification andmodelling by inserting hysteresis operators into mathematical equations, e.g.the Preisach model of ATP hysteresis [3] and the models for bacteria growthor prey-predator systems [8]. However, there do exist a variety of biologicalmodels without explicitly embedded hysteresis operators, but still distinctlydemonstrating hysteresis, e.g. systems introduced in [1, 9, 10, 11]. In addition,hysteresis is not new and has been widely observed in a variety of disciplines,such as material science, mechanics, electronics and economics. As a result ofyears of interdisciplinary work, the definitions of hysteresis are useful but dif-ferent in specific contexts. A stringently mathematical and universal definitionhas not yet appeared.Three essential components are usually used to characterize hysteresis: lag-ging, rate-independence and looping behaviour [12]. These can be easily under-stood from a simple input/output plot of hysteresis shown in Fig. 1. Laggingmeans that the output lags the input; rate-independence indicates that the out-put only depends on the values, not the rates of change, of the input; and thelooping behaviour implies that the output is affected by the previous values ofthe input, demonstrating a memory effect. Although all the three componentsare generally regarded as crucial features of hysteresis, contradictory examples2 igure 1: Binary hysteresis (also called relay) with output ∈ {− , } and width a + b . are not uncommon. Therefore, we need to understand the nature of hysteresis.Recently, a new definition was proposed from the dynamical system point ofview. Definition [12].
A hysteretic system is one which has (1) multiple stableequilibrium points and (2) dynamics that are considerably faster than the timescale at which inputs are varied.
The definition points out two main features of a dynamical system withthe property of hysteresis. One is multistability, another is dramatic changeswith respect to the slower input. Further, it implies that hysteresis by naturecan be understood by analysis of the multistability displayed in the bifurcationdiagram where dramatic transitions occur between multistable attractors byvarying the relatively constant bifurcation parameters. From this perspective,hysteresis dynamics has a strong link to the notion of bursting oscillations.Bursting oscillations, as an important neural activity, are usually studied viabifurcation theory and analysis of fast-slow systems, where the slow variablesare treated as parameters of the fast dynamics [6, 13]. In some examples, thefast subsystem exhibits multistability, which leads to a hysteretic loop visitingalternately one of two different attractors corresponding to resting and spikingstates, respectively.Moreover, in bifurcation theory terminology, hysteresis dynamics above has3n equivalent name, hysteresis bifurcation which is a type of reversible catas-trophe. Catastrophic bifurcation occurs when a microscopic variation of a pa-rameter triggers a macroscopic movement from one attractor to another. Ifthe system can be driven back to the initial attractor, the catastrophe is calledreversible. Fig. 2 depicts one kind of hysteresis bifurcation induced by twosaddle-node bifurcations at critical points p ∗ and p ∗ , respectively. Between p ∗ and p ∗ , the system is bistable with two stable equilibrium points. As the bi-furcation parameter p increases, the trajectory of the system slowly slides upalong the lower stable path (the solid curve from D ) until it reaches the rightknee, A . At this moment, it quickly jumps to point B , another attractor lead-ing to a higher stable branch. This jump is ”considerably faster than the timescale at which” the bifurcation parameter is varied. Likewise, the backwardsprocedure goes down along the upper stable branch. Upon reaching the leftknee, C , the system jumps to the lower branch and slide left. If the parameter p is varied back and forth, the trajectory of the system follow closely the loop A → B → C → D resulting in a reversible catastrophe. The loop is calledhysteretic loop, or briefly, hysteresis, which is analogous to Fig. 1 with the threecharacteristics: lagging, rate-independence and looping behaviour. Generally,attractors in the hysteretic loop could be a stable equilibrium point, a stable pe-riodic orbit, an attractive torus, even a strange attractor. Therefore, we suggestreplacing the phrase ”equilibrium points” in the definition with a more generalexpression, ”attractors”.Hysteresis is easy to be understood conceptually, but some of the attributesare quite difficult to study mathematically [13]. Our work aims to mathemat-ically investigate hysteresis from the dynamical system point of view. Bifur-cation and perturbation theories are used to analytically study qualitative ortopological changes of the trajectories of the nonlinear dynamics. We startwith classification of hysteresis initiated from all possible codimension-one bi-furcations of equilibrium points. Then, we perform a specific analysis on thetime-delayed FitzHugh-Nagumo neural system to show how the subcritical Hopfbifurcation and saddle-node bifurcation of limit cycles generate hysteresis. It4 igure 2: Hysteresis generated by two saddle-node bifurcations in the equation ˙ x = p + x − x .Modified from [7]. may be the simplest instance to form a hysteretic loop transiting between equi-librium points and limit cycles. The method of multiple timescales is used toderive a normal form up to the fifth order. While we focus on a specific system,the methods we use can be applied to any model involving ordinary or delaydifferential equations.The paper is organized as follows. In Sec. 2, we summarize the possible situ-ations where hysteresis bifurcation occurs and classify them. Section 3 presentsthe theoretical framework for hysteresis analysis of the time-delayed FitzHugh-Nagumo neuron. In Sec. 4, we validate our analytical results against solutionsobtained with the time domain simulation and numerical continuation. Finally,we conclude our findings in Sec. 5.
2. Classification of hysteresis bifurcation
A bifurcation indicates a transition from one qualitative type of dynamicsto another [14, 15]. Thus, we classify a hysteresis bifurcation by its generationmechanism, that is, what kinds of attractors are involved in this transition.
The neural system has two classic types of attractors: the resting state(quiescence) and periodic spiking. These two states correspond to a stable5quilibrium point and a limit cycle attractor, respectively. Switching betweentwo stable resting states has been observed in many experiments, e.g. [16].Hysteresis formed by transitions between equilibrium points can also be seen inthe Hodgkin-Huxley model for the squid axon where the transmembrane voltageis the bifurcation variable and the external potassium concentration acts as thebifurcation parameter [17].Three codimension-one bifurcations involve equilibrium points in a dynam-ical system: saddle-node bifurcation, transcritical bifurcation and pitchfork bi-furcation. Hysteretic loops generated by the three bifurcations are summarizedin Fig. 2, 3 and 4. Mathematically each example can be described by a one-dimensional nonlinear equation with the bifurcation parameter p .The bifurcation diagram in Fig. 2 is generated from a dynamical systemexpressed as ˙ x = p + x − x , where ˙ x = d x/ d t is the derivative of the variable x with respect to time t .From p + x − x = 0 (equilibrium condition) we find p = x − x . By using theextreme value theory, letting dp/dx = 0, one derives the mirrored critical valuesat p ∗ = − √ and p ∗ = √ . Then, through bifurcation analysis, we know thatthe system has two saddle-node bifurcations at p ∗ and p ∗ , respectively, with zeroeigenvalues at equilibria A ( x = − / √
3) and C ( x = 1 / √ χ = p ∗ − p ∗ = 2 √ − ( − √
39 ) = 4 √ . Consider the nonlinear equation˙ x = px + x − x . The bifurcation diagram in Fig. 3 shows multi-stability and hysteresis of thisdynamical system. By bifurcation analysis one can derive that two saddle-nodebifurcations occur at p ∗ = − / p ∗ = 0.Two symmetric hysteretic loops are generated with the range calculated as χ = p ∗ − p ∗ = 0 − ( − /
4) = 14 . igure 3: Two hysteresis bifurcations generated by two saddle-node bifurcations at p ∗ and asubcritical pitchfork bifurcation at p ∗ in the equation ˙ x = px + x − x . In addition, variables of biological systems mostly are positive, which canlead to different bifurcations. The hysteretic curve of Fig. 4 looks like the flippedcopy of the upper part of Fig. 3. However, the hysteresis generation mechanismis not the same. Let us consider the following nonlinear equation with x > x = − px + 4 x − x . By bifurcation analysis, we can see that a transcritical bifurcation at p ∗ = 0and a saddle-node bifurcation at p ∗ = 4 complete the hysteretic loop in Fig. 4with a width of χ = p ∗ − p ∗ = 4 − . A similar bifurcation diagram, except shifting to the right some units, can befound in the exploited population model [10].Besides codimension-one bifurcations, the cusp catastrophe, a codimension-two bifurcation, can give rise to hysteresis. Fig. 5 depicts a two-parameterbifurcation diagram of cusp from the equation˙ x = p + p x − x . Within the cusp-shaped grey region illustrated in the ( p , p ) parameter plane,there are three equilibrium points present. Outside of this region, there is onlyone equilibrium point. Compared with two macroscopic jumps in Fig. 2, 3 and7 igure 4: Hysteresis generated by a saddle-node bifurcation at p ∗ and a transcritical bifurca-tion at p ∗ in the equation ˙ x = − px + 4 x − x .Figure 5: Cusp bifurcation of the equation ˙ x = p + p x − x . Modified from [18].
8, the hysteretic loop here is formed by a smooth movement along the arrow C → D → A → B and a catastrophic transition from B to C . From the examples above we can see that the saddle-node bifurcation fre-quently appears in forming a hysteretic loop. Thus, it should not be surprisingthat the counterpart saddle-node bifurcation of limit cycles can also be involvedin hysteresis. Hysteresis involving movement between an equilibrium point anda limit cycle has been found experimentally in the squid axon and numericallyin the Hodgkin-Huxley model in response to the variation of the injected biascurrent [19]. This has been explained by the combination of a subcritical Hopfbifurcation and a saddle-node bifurcation of limit cycles which initiate hystereticdynamics in the model.A mathematical appreciation of such a hysteresis bifurcation can be achievedby reducing the system model to a fifth order normal form with the equation ofthe amplitude of periodic orbits,˙ r = α r r + β r r + c r r , (1)where α r , β r and c r are real values [14, 15]. The solutions of (1) are r = 0 , r , = (cid:118)(cid:117)(cid:117)(cid:116) − β r ± (cid:113)(cid:0) β r ) (cid:1) − α r c r c r , (2)where r = 0 corresponds to the equilibrium point, and the periodic orbit existswhen either r , or both have positive real values. The stability of the solutionsis evaluated by the sign of the Jacobian J = α r + 3 β r r + 5 c r r . (3)Fig. 6 illustrates a sketch of bifurcation diagram of (1). The system undergoesa subcritical Hopf bifurcation at the critical point p ∗ , where α r ( p ∗ ) = 0, anda saddle-node bifurcation of limit cycles at p ∗ , where the local extremum of α r ( p ) with respect to r reaches, that is, α r ( p ∗ ) = β r / (4 c r ). When p > p ∗ , the9ystem has only one stable equilibrium; when p < p ∗ , the system has an unstableequilibrium point and a stable limit cycle; a bistable region appears between p ∗ and p ∗ , where the system trajectory transits between a stable equilibrium pointand a stable limit cycle. In the next section, we will show how to derive thenormal form (1) and investigate such a hysteresis bifurcation by application toa time-delayed neural model. The relevant methods can be generalized to othersituations. Figure 6: Sketch of the bifurcation diagram in the equation ˙ r = − pr + r − r . Hysteresisinitiated by a saddle-node bifurcation of limit cycles at p ∗ and a subcritical Hopf bifurcationsat p ∗ . The arrows show one possible movement. Solid (dash) lines correspond to stable(unstable) solutions. Hysteresis may also occur due to a sequence of bifurcations that occurs in aparticular model. For example, [18] introduces a more complex hysteresis foundin a tritrophic food chain model. Here, catastrophic transitions between theequilibrium point and the prey-predator limit cycle are initiated by a transcrit-ical bifurcation and a homoclinic bifurcation.
We have seen that hysteresis may result from the coexistence of two sta-ble equilibrium points, it is also possible that two stable limit cycles coexist.The possible corresponding behaviours in a neural system are spikes fired with10ifferent periods. For example, it has been shown that bistability between in-phase and anti-phase oscillations can occur in models for systems of two coupledneurons [20, 21, 22, 23]. This has been linked to pitchfork bifurcations of limitcycles [20, 22] and subcritical Hopf bifurcations [23]. Similar phenomena havebeen found in many biological systems, including an ionic model of ventricularmembrane, where hysteretic transitions between periodic orbits with respective1:1 and 2:1 rhythms occurs at different driving frequencies [24]. Mathematically,one possibility to generate such a hysteresis bifurcation can be achieved by twosaddle-node bifurcations of limit cycles, similar to Fig. 2.
3. Hysteresis of the time-delayed FitzHugh-Nagumo neurons
In the section, we develop a theoretical analysis of hysteresis induced froma subcritical Hopf bifurcation and a saddle-node bifurcation of limit cycles byapplication to a time-delayed FitzHugh-Nagumo (FHN) neural system.The FHN model [25] is a two-dimensional simplification of the Hodgkin-Huxley equations describing spike generation. Although not clearly derivablefrom biology, the model has becomes a central model in mathematical neuro-science and is simple enough to allow analytical developments. Further, theinfluence of synaptic delays on system dynamics cannot be ignored or underes-timated in the field of neuroscience [26]. This has motivated many time-delayedneuron models, including the one proposed in [27].
The time-delayed FHN model introduced in [27] is modelled by a system ofdelay differential equations,˙ v = v ( t ) − v ( t ) − w ( t ) + µ (cid:0) v ( t − τ ) − v (cid:1) , ˙ w = ρ (cid:0) v ( t ) + a − bw ( t ) (cid:1) , (4)where ρ represents the timescale ratio between the membrane potential v andthe recovery variable w , the time delay τ > µ is the strength of thefeedback, positive for excitatory and negative for inhibitory feedback.11or µ = 0, the system (4) has an equilibrium point at ( v , w ) given by0 = v − v − b ( v + a ) ,w = ( v + a ) /b. (5)Moreover, under the following conditions0 < ρ < , < b < , − b/ < a < , (6)and 1 − bρ < v < bρ + 2 √ ρ, the equilibrium is unique and a stable focus [28]. Define x = v − v , y = w − w and the vector u = [ x, y ] T ( (cid:48) T (cid:48) means transpose), the equilibrium point istransformed to zero in the transformed model:˙ u = A u ( t ) + µB u ( t − τ ) + f (cid:0) u ( t ) (cid:1) , (7)where A = − v , − ρ, − ρb , B = , f (cid:0) u ( t ) (cid:1) = − v x ( t ) − x ( t )0 . In the section, we demonstrate the hysteresis bifurcation structure of (7).The method of multiple scales [29, 30] is used to obtain the normal form byexpanding the evolution of the dynamical system (7) around the Hopf location.Let us take µ as the bifurcation parameter and define µ = µ c + ε δ , (8)where µ c is the Hopf bifurcation point, 0 < ε (cid:28) µ c and δ takes the values ± µ = µ c in theform u ( t, ε ) = (cid:88) k =1 ε k U k ( T , T , T ) = (cid:88) k =1 ε k X k ( T , T , T ) Y k ( T , T , T ) . (9)12ere, T = t is the fast timescale, T = ε t and T = ε t are the first and secondslow timescales, respectively. The derivative with respect to t is transformedinto ddt = ∂∂T + ε ∂∂T + ε ∂∂T (10)Given by u = [ x, y ] T and (9), f (cid:0) u ( t ) (cid:1) in (7) is rewriten as f (cid:0) u ( t ) (cid:1) = ε − v X + ε − v X X − X + ε − v X − v X X − X X + ε − v X X − v X X − X X − X X ≡ (cid:88) k ≥ ε k f k ( U , U , U , U , U ) (11)In addition, the delay term u ( t − τ ) in (7) is expressed in terms of the scales T , T and T as u ( t − τ, ε ) = ε U τ + ε U τ + ε (cid:0) U τ − τ D U τ (cid:1) + ε ( U τ − τ D U τ ) + ε ( U τ − τ D U τ − τ D U τ ) (12)where U iτ = U i ( T − τ, T , T ) , i = 1 , ,
3. By substituting (8)-(12) into (7) andmatching these terms by their ε order, we obtain five differential equations asfollows: ∂ U ∂T − A U − µ c B U ( T − τ ) = 0 (13) ∂ U ∂T − A U − µ c B U ( T − τ ) = f (14) ∂ U ∂T − A U − µ c B U ( T − τ ) = − ∂ U ∂T − τ µ c B ∂ U ∂T ( T − τ )+ δ B U ( T − τ ) + f , (15)13 U ∂T − A U − µ c B U ( T − τ ) = − ∂ U ∂T − τ µ c B ∂ U ∂T ( T − τ )+ δ B U ( T − τ ) + f (16) ∂ U ∂T − A U − µ c B U ( T − τ ) = − ∂ U ∂T − τ µ c B ∂ U ∂T ( T − τ ) − ∂ U ∂T − τ µ c B ∂ U ∂T ( T − τ ) + δ B U ( T − τ ) − δ Bτ ∂ U ∂T ( T − τ ) + f (17) Next, we will solve the foregoing three equations (13)-(15) one by one toderive the third-order normal form of the Hopf bifurcations.Solving (13) is a typical nonlinear eigenvalue problem that has a generalsolution, U = W ( T , T ) q e iw c T + c.c. (18)Here, c.c. stands for the complex conjugate of the preceding terms and has theform of W ( T , T ) q e − iw c T ; W ( T , T ) is the complex amplitude depending onthe slow timescales, T and T and will be determined in later steps; s = iw c isthe eigenvalue at the Hopf point µ = µ c where the system is marginally stableand can be obtained by solving the characteristic equation det( M s ) = 0 with M s ≡ sI − A − µ c Be − sτ . (19) q is the corresponding eigenvector, which is not unique. Here, we use a generalnotation to the 2-D eigenvector, q = X W Y W . (20) q can be also taken to be q = ρρb + iw c , (21)for the specific FHN system (7). 14ubstituting (18) into (14) yields ∂ U ∂T − A U − µ c B U ( T − τ ) = | W | F | W | + (cid:0) W F W e iw c T + c.c. (cid:1) , (22)where F | W | = − v | X W | = − v | X W | (23a) F W = − v ( X W ) = − v ( X W ) . (23b)When using (21), then gives F | W | = − v , F W = − v . (24)Here, the superscripts of F show the dependence on the amplitude W and sub-scripts indicate the corresponding equations (13)-(17). The notation is borrowedfrom [31] and will be employed throughout the paper.Assume U has the same form as the forcing term in (22) U = | W | U | W | + (cid:0) W U W e iw c T + c.c. (cid:1) . (25)Substituting it into (22) and balancing similar terms yield specific expressions(see details in the appendix).Next, substituting (18) and (25) into the differential equation (15), we have ∂ U ∂T − A U − µ c B U ( T − τ ) = (cid:18) − ∂W∂T H q + W F W + | W | W F | W | W (cid:17) e iw c T + W F W e iw c T + c.c. (26)where H = I + τ µ c Be − iw c τ . (27)To guarantee (26) has solutions, a solvability condition has to be satisfied. Thecondition is that the sum of the resonant forcing terms, that is, the terms with15 iw c T on the right-hand side of (26), should be orthogonal to every solution ofthe adjoint homogeneous problem [32]. In this case, the adjoint problem is M † iw c q † = 0 , where M † iw c is the Hermitian of the matrix M iw c and has the form M † iw c ≡ − iw c I − A T − µ c B T e iw c τ . (28)Taking the inner product of the resonant forcing terms of (26) with q † yieldsthe solvability condition, (cid:68) q † , − ∂W∂T H q + W F W + | W | W F | W | W (cid:69) = 0 , (29)which is rewritten as ∂W∂T = α W + β | W | W. (30)Here, the complex values α and β are calculated as α = (cid:104) q † , F W (cid:105)(cid:104) q † , H q (cid:105) , β = (cid:104) q † , F | W | W (cid:105)(cid:104) q † , H q (cid:105) . (31)For easy of calculation, we can choose a unique q † by imposing the followingcondition, (cid:104) q † , q (cid:105) = (cid:104) q , q † (cid:105) = ¯ q T q † = 1 . In addition, the solution of (26) has the form U = W U W e iw c T + | W | W U | W | W e iw c T + W U W e iw c T + c.c. (32)The appendix gives the specific expressions.Eq. (30) is a third-order normal form usually used to understand the Hopfbifurcation [14, 15]. Substituting W ( t ) = r ( t ) e iθ ( t ) into (30) and taking realand imaginary parts of the resulting equation, we get the expressions in polarcoordinates as ˙ r = Re( α ) r + Re( β ) r , (33a)˙ θ = Im( α ) + Im( β ) r , (33b)16here α and β are called the Landau coefficients. The amplitude equation(33a) has solutions r = 0 , r = (cid:115) − Re( α )Re( β ) . (34)The stability of the solutions is determined by the sign of the eigenvalue λ evaluated at the solutions. These are λ ( r ) = Re( α ) , λ ( r ) = − α ) (35)Based on bifurcation theory, a subcritical Hopf bifurcation occurs whenRe( β ) > Figure 7: Sketch of the bifurcation diagram obtained from ˙ r = r − pr . Substituting (18), (25) and (32) into the differential equation (16), we have ∂ U ∂T − A U − µ c B U τ = | W | F | W | + | W | F | W | + (cid:16) W F W e iw c T + | W | W F | W | W e iw c T + W F W e iw c T + c.c. (cid:17) , (36)17ne can see that there are no resonant terms on the right-hand side of theequation. Then, used the ansatz U = | W | U | W | + | W | U | W | + (cid:16) W U W e iw c T + | W | W U | W | W e iw c T + W U W e iw c T + c.c. (cid:17) , (37)a set of expressions is readily derived in the appendix.Similarly, by substituting the previous solutions of U to U , the differentialequation (17) can be rewritten as ∂ U ∂T − A U − µ c B U ( T − τ ) = (cid:18) − ∂W∂T H q + W F W + | W | W F | W | W + | W | W F | W | W (cid:17) e iw c T + (cid:16) | W | W F | W | W + W F W (cid:17) e iw c T + W F W e iw c T + c.c. (38)One can see that there are resonant terms with e iw c T on the right-hand side ofthe equation. Therefore, by applying the solvability condition on the resonantterms as before, we obtain the fifth-order normal form at the timescale T , ∂W∂T = α W + β | W | W + c | W | W, (39)where the Landau coefficients are α = (cid:104) q † , F W (cid:105)(cid:104) q † , H q (cid:105) , β = (cid:104) q † , F | W | W (cid:105)(cid:104) q † , H q (cid:105) , c = (cid:104) q † , F | W | W (cid:105)(cid:104) q † , H q (cid:105) . (40)Eventually, the final fifth-order normal form is derived by combining (30),(39) and using the scaling T = ε T [31] asd W d T = ∂W∂T + ∂W∂T ∂T ∂T = (cid:0) α + ε α (cid:1) W + (cid:0) β + ε β (cid:1) | W | W + ε c | W | W (41)Substituting the polar representation W = re iθ , we have the normal form withthe amplitude and phase of limit cycle solutions,d r d T = Re( α ) r + Re( β ) r + Re( c ) r , (42a)d θ d T = Im( α ) + Im( β ) r + Im( c ) r . (42b)18ere, α = α + ε α , β = β + ε β and c = ε c . One can see that (42a) hasthe same expression as (1), where hysteretic transitions between a equilibriumpoint and a limit cycle is generated by a subcritical Hopf bifurcation and asaddle-node bifurcation of limit cycles, as shown in Fig. 6.The generation mechanism of hysteresis is simpler when only involving equi-librium points. Complexity increases with the involvement of limit cycles. Here,we have focused on a relatively analytically tractable case and presented a an-alytical framework by applying the method of multiple scales to the delayedFHN neuron model. Our approach can be easily extended to other systems, orto other bifurcation parameters, such as τ , to investigate the impact of timedelays on dynamics.
4. Numerical analysis
Using the method of multiple scales to derive the normal form, even a low-order one, may be a lengthy and tedious process. However, such a procedure isstandard and can be automatised with symbolic solvers. See for example [33].In this section, we show some numerical results to confirm the analyticalexpressions. The parameter values are chosen as a = 0 . b = 0 .
8, commonlyused in the literature. As for ρ = 0 .
08 and τ = 60, originally used in [27], wehave found by DDE-BIFTOOL [34] that the system undergoes more complicatedbifurcations, including not only the subcritical Hopf bifurcation and saddle-nodebifurcation of limit cycles, but also the period-doubling bifurcation of limitcycles and torus bifurcation. Therefore, we set ρ = 0 . τ = 15 to obtainthe relatively simpler and illustrative bifurcation structure to show hystereticdynamics.The system with large time delay generally has several Hopf bifurcationsover short parameter intervals [35]. It also occurs in our system. For the nu-merical analysis, we choose the Hopf point µ c = − . µ increases and passes through µ c . The correspondingeigenvalues at the point are ± w c i = ± . i and the first Lyapunov coeffi-19ient is L = 0 . µ f = − . Figure 8: Branches of periodic solutions and equilibrium points from numerical continuationwith respect to µ . A hysteresis bifurcation is induced by a subcritical Hopf bifurcation (SH) at µ c = − . µ f = − . a = 0 . b = 0 . ρ = 0 . τ = 15. The symbol ’s’ denotes stable and ’u’ denotesunstable. On the other hand, we found the Landau coefficients of the 3rd-order normalform from (31): when µ < µ c , α = 0 . − . i and β = 0 . − . i ;when µ > µ c , α change the sign and becomes − . . i , whereas β keep the same. Based on (34) and (35), we can see that in the vicinity of thecritical point µ c = − . µ < µ c : α = − . . i , β = − . − . i and c = − . e − − . e − i . When µ > µ c , α and c remain the same, whereas β changes sign across the Hopf point. We choose ε = | µ c − µ f | / . (cid:28) µ < µ c : J ( r ) = 0 . > , J ( r ) = − . < µ > µ c : J ( r ) = − . < , J ( r ) = − . < , J ( r ) = 0 . > . This shows that before the Hopf point µ c , there are one unstable equilibriumpoint and one stable periodic orbit; after µ c , the equilibrium point becomesstable and there exist two periodic solutions, one is stable and the other isunstable. The analytical results are consistent with the numerical continuationin Fig. 8, where a hysteretic loop is formed between µ c = − . µ f = − . Figure 9: Numerical simulations showing stable behaviours for (a) µ = − .
81 and (b) µ = − .
47. The initial conditions are v ( t ) = 0 and w ( t ) = w . Other parameter values are as inFig. 8. In addition, Fig. 9 shows the numerical solutions for µ = − .
81 and for µ = − .
47, respectively to show stable behaviours beyond the hysteresis region.21 igure 10: Numerical simulations showing bistability between the equilibrium point solutionand periodic orbit solution. Parameter values are as in Fig. 8 and µ = − .
6. Switchingsbetween the attractors are achieved by applying two perturbations to the parameter a asfollows: ∆ a = 0 .
01, 501 ≤ t ≤ a = − .
2, 912 ≤ t ≤ Fig. 10 gives numerical simulations to demonstrate bistability for µ = − .
5. Conclusions
Over the years, there have been a substantial number of purely experimen-tal work with phenomenological descriptions of the remarkable dynamical be-haviour: hysteresis. A mathematical appreciation of such dynamics must dealwith the analysis from the dynamical system point of view by using bifurcationand perturbation theories. In this paper, we have summarized some types of hys-teresis bifurcations and shown biological examples to illustrate these phenom-ena. We have classified hysteresis in terms of catastrophic transitions betweendifferent types of attractors. Hysteretic dynamics can be easily appreciatedwhen only involving equilibrium points. Situations become complicated wheninvolving cycles, multiple attractors and/or complex, even global bifurcations.Correspondingly, the theoretical analysis becomes more difficult.22e have theoretically investigated the instance where hysteretic movementsbetween the equilibrium point and the limit cycle are initiated by a subcrit-ical Hopf bifurcation and a saddle-node bifurcation of limit cycles. We haveapplied the method of multiple scales in the time-delayed FitzHugh-Nagumoneural system close to the Hopf point and reduced the governing equations toa fifth-order normal form without delays. From the normal form, we can pre-dict the amplitude and frequency of stable and unstable limit cycles, and theregion of hysteresis with bistability. Before the expansion, we need informationabout the value of the bifurcation parameter at the Hopf point, the marginallystable eigenvalues and the corresponding direct and adjoint eigenvectors. Thelater process of analytical expansion may be lengthy and tedious, but the proce-dure is standard and can be automatically realized by symbolic solvers, such asMaple [33]. Our theoretical results have shown good agreement with numericalsimulations and continuation.In addition, we should point out that the normal form derived from the pa-rameter expansion is strictly valid only for the vicinity of the Hopf point, where ε (cid:28) cknowledgments This work benefited from the support of the Natural Science and EngineeringResearch Council of Canada.
Appendix A. Expressions when solving (14)-(17)
1. Expressions of (25) U | W | = M − F | W | ≡ X | W | Y | W | , (A.1a) U W = M − iw c F W ≡ X W Y W , (A.1b)where the non-singular matrices M and M iw c are from (19) with s = 0and s = 2 iw c , respectively.2. Expressions of (32) U W = M +1 iw c (cid:0) F W − α H q (cid:1) , (A.2a) U | W | W = M +1 iw c (cid:0) F | W | W − β H q (cid:1) , (A.2b) U W = M − iw c F W , (A.2c)where M +1 means pseudo inverse of the matrix M because the matrix M iw c is singular. F W = δ X W e − iw c τ , (A.3a) F | W | W = − v (cid:16) X W X | W | + X W X W (cid:17) − X W | X W | , (A.3b) F W = − v X W X W − (cid:0) X W (cid:1) . (A.3c)Here, X means complex conjugate.24. Expressions of (37) F | W | = − I + τ µ c B )Re( β ) U | W | + v (cid:16) X | W | (cid:17) + 2 v (cid:12)(cid:12)(cid:12) X W (cid:12)(cid:12)(cid:12) − v X W X | W | W − (cid:0) X W (cid:1) X W + 2 (cid:12)(cid:12) X W (cid:12)(cid:12) X | W | , (A.4a) F | W | = − I + τ µ c B )Re( α ) U | W | + δ X | W | − v X W X W , (A.4b) F W = − I + τ µ c Be − iw c τ ) α U W − v X W X W , (A.4c) F | W | W = − I + τ µ c Be − iw c τ ) β U W − v X | W | X W − v X W X | W | W + 2 v X W X W − X | W | (cid:0) X W (cid:1) + 2 (cid:12)(cid:12) X W (cid:12)(cid:12) X W , (A.4d) F W = − v (cid:16) X W (cid:17) − v X W X W − (cid:0) X W (cid:1) X W (A.4e) U | W | = M − F | W | , (A.5a) U | W | = M − F | W | , (A.5b) U W = M − iw c F W , (A.5c) U | W | W = M − iw c F | W | W , (A.5d) U W = M − iw c F W . (A.5e)25. Expressions of (38) F W = α (cid:0) I − τ µ c Be − iw c τ (cid:1) U W + δ Be − iw c τ ( U W − α τ U W ) , (A.6a) F | W | W = (cid:0) I − τ µ c Be − iw c τ (cid:1) (cid:16) β U W + (cid:0) α + 2Re( α ) (cid:1) U | W | W (cid:17) + δ BU | W | W e − iw c τ − β τ δ BU W e − iw c τ − v X | W | X W + X W X W + X | W | X W + X W X W − | X W | X W + ( X W ) X W , (A.6b) F | W | W = (cid:0) I − τ µ c Be − iw c τ (cid:1) (cid:0) β + 2Re( β ) (cid:1) U | W | W − v X | W | X | W | W + X W X | W | W + X W X W − v X | W | X W + X W X | W | W − (cid:16) X W (cid:17) X W + X W (cid:16) X | W | (cid:17) + 2 X | W | X W X W (A.6c) References [1] D. Angeli, J. J. Ferrell, E. D. Sontag, Detection of multistability, bifurca-tions, and hysteresis in a large class of biological positive-feedback systems,Proc Nat Acad Sci U S A 101 (2004) 1822–1827.[2] B. T. Andrews, D. T. Capraro, J. I. Sulkowska, J. N. Onuchic, P. A. Jen-nings, Hysteresis as a marker for complex, overlapping landscapes in pro-teins, J Phys Chem Lett 4 (2013) 180–188.[3] H. R. Noori, Substantial changes in synaptic firing frequencies induced byglial atp hysteresis, Biosystems 105 (2011) 238–242.264] D. W. DelMonte, T. Kim, Anatomy and physiology of the cornea, JCataract Refract Surg 37 (2011) 588–598.[5] J. Ramos, S. Lynch, D. Jones, H. Degens, Hysteresis in muscle, Interna-tional Journal of Bifurcation and Chaos 27 (1) (2017) 1730003–(1–16).[6] E. M. Izhikevich, Neural excitability, spiking and bursting, InternationalJournal of Bifurcation and Chaos 10 (6) (2000) 1171–1266.[7] H. R. Noori, Hysteresis phenomena in biology, 1st Edition, SpringerBriefsin Mathematical Methods, Springer-Verlag Berlin Heidelberg, 2014.[8] Kopfov´a., Hysteresis in biological models, Journal of Physics: ConferenceSeries 55 (2006) 130–134.[9] T. S. Gardner, C. R. Cantor, J. J. Collins, Contrustion of a genetic toggleswitch in escherichia coli, Nature 403 (2000) 339–342.[10] I. Noy-Meir, Stability of grazing systems: An application of predator-preygraphs, J. Ecol. 63 (1975) 459–483.[11] Y. A. Kuznetsov, S. Muratori, S. Rinaldi, Homoclinic bifurcations in slow-fast second-order systems, Nonlinear Analysis 25 (7) (1995) 747–762.[12] K. Morris, What is hysteresis, Applied Mechanics Reviews 64 (5) (2012)050801–(1–14).[13] E. M. Izhikevich, Dynamical systems in neuroscience : the geometry ofexcitability and bursting, Computational neuroscience, MIT Press, Cam-bridge, Mass, 2007.[14] J. Guckenheimer, p. Holmes, Nonlinear Oscillations, Dynamical Systemsand Bifurcations of Vector Fields, Springer-Verlag, New York, 1983.[15] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer, NewYork, 1998. 2716] I. Tasaki, Demonstration of two stable states of the nerve membrane inpotassiumrich media, Journal of Physiology 148 (1959) 306–331.[17] K. Aihara, G. Matsumoto, Two stable steady states in the hodgkin-huxleyaxons, Biophysical Journal 41 (1) (1983) 87–89.[18] F. Dercole, S. Rinaldi, Dynamical systems and their bifurcations, in: Ad-vanced Methods of Biomedical Signal Processing, John Wiley & Sons, Ltd,2011, pp. 291–325.[19] R. Guttman, S. Lewis, J. Rinzel, Control of repetitive firing in squid axonmembrane as a model for a neuroneoscillator, The Journal of Physiology305 (1980) 377–395.[20] S. A. Campbell, I. Kobelevskiy, Phase models and oscillators with time de-layed coupling, Discrete & Continuous Dynamical Systems-A 32 (8) (2012)2653–2673.[21] S. Kunec, A. Bose, Role of synaptic delay in organizing the behavior ofnetworks of self-inhibiting neurons, Phys. Rev. E. 63 (2) (2001) 021908.[22] Y. Park, B. Ermentrout, Weakly coupled oscillators in a slowly varyingworld, Journal of computational neuroscience 40 (3) (2016) 269–281.[23] H. Ryu, S. A. Campbell, Effect of time delay on the synchronization ofexcitatory-inhibitory neural networks, bioRxiv doi:10.1101/2020.08.30.274662 .[24] A. R. Yehia, D. Jeandupeux, F. Alonso, M. R. Guevara, Hysteresis andbistability in the direct transition from 1:1 to 2:1 rhythm in periodicallydriven single ventricular cells, Chaos 9 (4) (1999) 916–931.[25] R. Fitzhugh, Impulses and physiological states in theoretical models ofnerve membrane, Biophysical Journal 1 (6) (1961) 445–466.2826] A. Ghosh, Y. Rho, A. McIntosh, R. K¨otter, V. Jirsa, Cortical networkdynamics with time delays reveals functional connectivity in the restingbrain, Cognitive neurodynamics 2 (2) (2008) 115.[27] R. E. Plant, A fitzhugh differential-difference equation modeling recurrentneural feedback, SIAM Journal of Applied Mathematics 40 (1) (1981) 150–162.[28] R. Fitzhugh, Mathematical models of excitation and propagation in nerve,in: H. P. Schwan (Ed.), Biological Engineering, McGraw Hill, New York,1969, pp. 1–86.[29] S. L. Das, A. Chatterjee, Multiple scales without center manifold reductionsfor delay differential equations near hopf bifurcations, Nonlinear Dynamics30 (2002) 323–335.[30] R. S. Johnson, Singular Perturbation Theory Mathematical and AnalyticalTechniques with Applications to Engineering, 1st Edition, Mathematicaland Analytical Techniques with Applications to Engineering, Springer US,New York, 2005.[31] A. Orchini, G. Rigas, M. P. Juniper, Weakly nonlinear analysis of ther-moacoustic bifurcations in the rijke tube, Journal of Fluid Mechanics 805(2016) 523–550.[32] J. T. Oden, L. F. Demkowicz, Applied functional analysis, 2nd Edition,CRC Press, 2010.[33] N. E. Sanchez, The method of multiple scales: asymptotic solutions andnormal forms for nonlinear oscillatory problems, Journal of Symbolic Com-putation 21 (1996) 245–252.[34] J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, D. Roose, Dde-biftoolv.3.1.1 manual bifurcation analysis of delay differential equations.URL http://arxiv.org/abs/1406.7144http://arxiv.org/abs/1406.7144