# Hopf Bifurcation for a FitzHugh–Nagumo Model with Time Delay in a Network

Authors: Suxia Wang
Journal: Complexity
 (2021)

---

## Abstract

A reaction diffusion system is used to study the interaction between species in a population dynamic system. It is not only used in a population dynamic system with the diffusion phenomenon but also used in physical chemistry, medicine, and animal and plant protection. It has been studied by more and more scholars in recent years. The FitzHugh–Nagumo model is one of the most famous reaction-diffusion models. This article takes a deeper look at a FitzHugh–Nagumo model in a network with time delay. Firstly, we studied the linear stability of the equilibrium, then the existence of Hopf bifurcation is given, and finally, the stability of the Hopf bifurcation is introduced.


---

## Body

## 1. Introduction

Reaction-diffusion equations are widely used in ecology and biology fields [1–3]. The Hopf bifurcation is very useful in practice; in [4–7], the turing instability is analyzed. Turing [8] first introduced the coupled reaction-diffusion equations; he pointed out that the stable process could evolve into an instability with diffusive effects. As a time- and space-breaking symmetry, Turing–Hopf bifurcation is a well-known mechanism for generating spatiotemporal patterns [9, 10].The FitzHugh–Nagumo model is one of the most famous reaction-diffusion models; Hodgkin and Huxle first introduced the FitzHugh–Nagumo model. They used the model to study electrical impulses on nerve fibers [11–13]. In the Fitzhugh–Nagumo model, the variables of equations for the voltage r and s have many forms, and the following form is commonly used:(1)∂r∂t=∂2r∂x2−r+I+fr,∂s∂t=ar−bs.The functionfr is positive feedback, and it is a third-order polynomial; however, the function r is negative feedback. In [14–18], the pattern formation in the FitzHugh–Nagumo system is studied frequently. In 2001, Medvedev and Kopell [19] introduced the FitzHugh–Nagumo oscillators with coupled electrically, and the synchronization and transience are obtained; in 2004, Kostova et al. [20] investigated the revised FitzHugh–Nagumo equation. The FitzHugh–Nagumo model with time delays is considered by C. J. in 2013 [21]. Based on the abovementioned research, we discussed the FitzHugh–Nagumo model with delay as follows:(2)∂r∂t=−r3t,x+α+1r2t,x−αrt,x−st−τ,x+DrΔr,∂s∂t=aβrt,x−ast,x+bDsΔs.In the abovementioned system,rt,x is the membrane potential, st,x is the recovery variable, Dr and Ds are the diffusive coefficients of rt,x and st,x, respectively,α is the excitatory threshold, b is the diffusion coefficient, a is the excitability and β,a are parameters, they can change the rest-state and dynamics, △ denotes the Laplacian operator which describes the diffusion, and τ represents the time delay. For detail, one can see [22]. This model is based on a spatially continuous domain, that is, the entire landscape is spatially continuous; in this paper, we introduce a network that owns system (2).(3)r˙it=−ri3t+α+1ri2t−αrit−sit−τ+Dr∑j=1nLijrjt,s˙it=aβrit−asit+Ds∑j=1nLijsjt,where G=gijn×n is an adjacency matrix, whose entry gij equals 1 if there exists an edge from node j to node i and 0 otherwise. It is obvious that gii=0 for i=1,2,…,n. We assume that G is strongly connected, i.e., for any pair of distinct nodes, there exists a path from one node to the other. It is known that the adjacency matrix G is irreducible [23]. The definition of the Laplacian matrix L is as follows:(4)L=−∑k=1ng1kg12…g1ng21−∑k=1ng2k…g2n…………gn1gn2…−∑k=1ngnk.The rest of this paper is designed as follows. In Section2, we give the linear stability of the equilibrium and the existence of Hopf bifurcation. In Section 3, the stability of the Hopf bifurcation is introduced.


## 2. Stability of the Equilibrium and Hopf Bifurcation

By simple calculation, we know that equation (3) has three equilibriums 0,0;r^1,s^1; and r^2,s^2, where(5)r^1=α+1+α+12−4β2,s^1=βr^1,r^2=α+1−α+12−4β2,s^2=βr^2.The interior equilibrium exists if and only if(6)α+12>4β.Now, we focus on the stability properties of the coexisting equilibriumr^1,s^1,r^2,s^2, we define the small perturbations δri=ri−r^;,δsi=si−s^, and system (3) can be expressed as follows:(7)δr˙iδs˙i=J1δriδsi+J2δri,τδsi,τ+D∑j=1nLijδrj∑j=1nLijδsj,where δri,τ=δrit−τ and(8)J1=a11−1bβ−b,J2=0−100,D=Dr00Ds,where a11=−3r^2+2α+1r^−α.As in [24], we use 0=Λ1>Λ2>⋯>Λn to express the eigenvalues of the Laplacian matrix L and Φγ:γ=1,2,…,n to express each eigenvector relevant to a topological eigenvalue Λγ, so we can get ∑j=1nLijΦjγ=ΛγΦiγ. Let T be the transpose and δrit,δsitT express any small disturbance in the equilibrium state; we can have the following basis decomposition:(9)δritδsit=∑α=1nc1γc2γeλγtΦγα,i=1,2,…,n.Inserting (9) into (7), using the orthogonality of the eigenvectors, we get(10)λγc1γc2γeλγt=J1c1γc2γeλγt+J2c1γc2γeλγt−τ+ΛγDc1γc2γeλγt.So, the characteristic equation can be written as(11)Δλγ,τ=detλγI−J1+e−λγτJ2+ΛγD.Substituting (8) into (11), the characteristic equation of (3) can be written as(12)λγ2−b+ΛγDs+a11+ΛγDrλγ+bβ+bβb−λγτ=0.Lemma 1.
If(13)3r^2+α>b,2bβ>A12,and (6) holds, then we can get the following:(1)
Whenτ=0, all characteristic roots of (12) have negative real parts(2)
Whenτ>0, the characteristic equation (12) has a pair of purely imaginary roots ±iω∗ at τ=τj, where(14)ω∗2=2bβ−A12,A1=b+ΛγDs+a11+ΛγDr,τj=1ω∗2jπ+arccosω∗2−bβbβ,j=0,1,2….Proof.
(1) Whenτ=0, the characteristic equation (12) can be reduced to(15)λγ2−λγe+ΛγDs−3r^2+2α+1r^−a+ΛγDr+2bβ=0.
Considering (13), we can obtain that b+ΛγDs−3r^2+2α+1r^−α+ΛγDr<0 and 2bβ>0. Therefore, the characteristic equation (12) has a root with a negative real part (3). Suppose ±iω are a pair of pure imaginary roots of (12). Substituting ±iω into (12), we can get(16)−ω2+iA1ω+bβ+bβcosωτisinωτ=0,and through separating real and imaginary parts of the root, we can get(17)bβcosωτ+bβ−ω2=0,A1ω±bβsinωτ=0.
It is noticed thatcos2ωτ+sin2ωτ=1, and we can get(18)ω2ω2+A12−2bβ=0.
Equation (13) has a unique positive real root ω∗2 if and only if (13) holds.Theorem 1.
(1)
Whenτ∈0,τ0, system (3) is locally stable at the equilibrium u^,v^ if (6) and (13) hold(2)
Whenτ∈τ0,∞, system (3) is unstable at the equilibrium u^,v^ if (6) and (13) hold(3)
We can find the constantσ>0 and a smooth curve λτ:τ0−σ,τ0+σ⟶C when λτ0=iω∗ and Δiω∗,τ0=0 for all τ∈τ0−σ,τ0+σProof.
We refer to the Theorem 3.6 in the work of Ruan [25]. In consideration of that τ0 is the smallest value when the root of equation (12) has zero real part, when t∈0,τ0, the real parts of the characteristic roots are all negative. Through using G. J. Buther’s Lemma [26], we know that u^,v^ is locally stable when τ∈0,τ0.
Now, we check the transversality condition(19)ddτReλτ|τ=τ0>0.
Supposeλγ=μ+iω is the eigenvalue of (12), and by substituting λγ into (12), we can obtain(20)μ2−ω2+A1μ+bβ+bβb−μτcosωτ=0,2μω+A1ω−bβb−μτsinωτ=0.
Differentiating (20) with respect to τ, we have(21)2μ+A1−bβτb−μτcosωτdμdτ−2ω+bβτb−μτsinωτdωdτ=bβb−μτωsinωτ+μcosωτ,2ω+bβτb−μτsinωτdμdτ+2ω+A1−A2τb−μτcosωτdωdτ=bβb−μτωcosωτ+μsinωτ.
Eliminating the termsdω/dτ from the abovementioned equation, we get(22)2μ+A1−bβτe−μτcosωτ2+2ω+bβτe−μτsinωτ2dμdτ=bβb−μτωsinωτ+μcosωτ2μ+A1−bβτb−μτcosωτ+bβb−μτμsinωτ+ωcosωτ2ω+bβτb−μτsinωτ.
To verifydμ/dτ>0, we just need to verify that(23)ωsinωτ+μcosωτ2μ+A1−bβτb−μτcosωτ+μsinωτ+ωcosωτ2ω+bβτb−μτsinωτ>0.
Whenτ=τ0,μ=0,ω=ω∗, (20) becomes(24)bβcosω∗τ0=ω∗2−bβ,bβsinω∗τ0=A1ω∗.
By some simple calculations, we can get(25)ω∗4+A12−2bβω∗2=0.
So, (23) becomes(26)ω∗sinω∗τ0A1∓bβτ0cosω∗τ0+ω∗cosω∗τ02ω∗±bβτ0sinω∗τ0>0.
We must verify(27)A1ω∗sinω∗τ0+2ω∗cosω∗τ0>0.
Substituting (25) into (27) yields(28)1ebA12ω∗2−1eb2ω∗2A2−ω∗2>0.
Sincebβ>0, we need only to verify that(29)2ω∗4+A12−2bβω∗2=2ω∗4−ω∗4=ω∗4>0.
The abovementioned inequality holds, so system (3) undergoes a Hopf bifurcation at τ=τ0.
By applying the implicit function toΔλ,τ, we can find the constant σ>0 and a smooth curve λτ:τ0−σ,τ0+σ→C when λτ0=iω∗ and Δiω∗,τ0=0 for all τ∈τ0−σ,τ0+σ.
Now, by using proposition 6.5 in [28], we obtain thatu^,v^ is locally unstable for τ∈τ0,∞.


## 3. The Stability of the Hopf Bifurcation

As a matter of convenience, letτ=τ0+ν. Then, ν=0 is the Hopf bifurcation value for system (3). Let ri−r^=r¯i,si−s^=s¯i. For convenience, we leave off the bars of r¯iand s¯i. System (3) can be rewritten as(30)r˙it=a11rit−si−sit−τ+2α+1−6r^ri2t+Dr∑j=1nLijrjt,s˙it=bβrit−bsit+Ds∑j=1nLijsjt.We rewrite the abovementioned equation as follows through using the orthogonality of the eigenvectors.(31)r˙it=a11rit−si−sit−τ+2α+1−6r^ri2t+DrΛαrit,s˙it=bβrit−bsit+DsΛαsit.We useXt=rit,sitT to denote the solution of (31) and let Xtθ=Xt+θ,θ∈−τ,0, where T denotes the transpose, Xt∈C−τ,0,R2≜C. System (31) can become the functional differential equation as follows:(32)X˙t=LνXt+RνXt,where Lν:C→R2,Rν:C→R2; they have the following forms:(33)Lνϕ=M0ϕ10ϕ20+M1ϕ1−τϕ2−τ,M0=a11+DrΛγ−1bβ−b+DsΛγ,M1=0−100,Rνϕ=−6r^+2α+1ϕ1200.We can let(34)ηθ,ν=M0δθ+M1δθ+τ,where δ⋅ is the Dirac function, and we can get(35)Lνϕ=∫−τ0ϕθdηθ,ν,forϕ∈C,where(36)dηθ,ν=M0δθdθ+M1δθ+τdθ.Then, two operatorsAν and Fν on C1−τ,0,R2≜C1 can be defined by(37)Aνϕθ=dϕdθ,θ∈−τ,0,∫−τ0ϕξdηξ,ν,θ=0,Fνϕθ=0,θ∈−τ,0,Rνϕ,θ=0.Then, (32) is transformed into(38)X˙t=AνXt+FνXt,(39)Aν∗ψθ∗=−dψdθ∗,θ∗∈0,τ,∫−τ0ψ−ξdηTξ,ν,θ∗=0,where ψ∈C10,τ,R2. We can define the bilinear form by taking the solution space of (38) as the complex space C2.(40)ψθ∗,ϕθ=ψ¯T0ϕ0−∫θ=−τ0∫ξ=0θψ¯Tξ−θdηθϕξdξ,for ϕ∈C−τ,0,C2,ψ∈C10,τ,C2, where ηθ=ηθ,0.From the proof in Lemma1, we can get that ±iω∗ are eigenvalues of (12); thus, ±iω∗ are the eigenvalues of A0 and A0∗, respectively. We set(41)A0qθ=iω∗qθ,andA0∗q∗θ∗=−iω∗q∗θ∗,and we know that q∗,q=1 and q∗,q¯=0. Let(42)qθ=q0biw∗θ=q1q2biw∗θ,andq∗θ∗=q∗0biw∗θ∗=q1∗q2∗biw∗θ∗,for θ∈−τ,0,θ∗∈0,τ.So, we have(43)qθ=1bβiω∗+bbiw∗θ,andq∗θ∗=ρbβ−iω∗−a110biw∗θ∗.In order to ensureq∗s,qθ=1, we should determine the value of ρ, and from (40), we get(44)q∗s,qθ=q¯∗0q0−∫θ=−τ0∫ξ=0θq¯∗ξ−θdηθqξdξ,=ρ¯q¯1∗q1,q2T−∫−τ0∫0θρ¯q¯1∗,q¯2∗b−iω∗ξ−θdηθq1,q2Tbiω∗ξdξ,=ρ¯q¯1∗+q2+τ0b−iω∗τ0a13q2,=ρ¯q¯1∗+q2−τ0b−iω∗τ0q2.Therefore, we can let(45)ρ¯=1q¯1∗+q2−τ0b−iω∗τ0q2.Then, we construct a coordinate for the center manifoldC0 at v=0 by using q and q∗; the method is similar to that in [27]. Letxt=xtθ be a solution of (38), and we define(46)zt≜q∗,xt,(47)Ut,θ≜xtθ−2Reztqθ.We further get(48)Ut,θ=Uzt,z¯t,θ≜U20θz22+U11θzz¯+U02θz¯22+U30θz36+⋯,where z and z¯ are the local coordinates of the center manifold C0 in directions q and q∗, respectively.We only take real solutions into account. For a solutionXt∈C0 of (26) with ν=0, from (47), we have(49)z˙t=q∗,x˙t,=iω∗z+q∗¯0F0Uz,z¯,θ+2Reztqθ,≜iω∗z+q∗¯0F0z,z¯.Thus, (49) can become(50)z˙t=iω∗z+gz,z¯,(51)gz,z¯=q∗¯0F0z,z¯,≜g202z2+g11zz¯+g022z¯2+g306z3+g212z2z¯+g122zz¯2+g036z¯3+⋯.Lemma 2.
If we set(52)k1=−6u^+2α+1,then gij can be explicitly expressed by (55).Proof.
Noticing (47) and (48), we obtain(53)xtθ=Ut,θ+2Reztqθ,=U20θz22+U11θzz¯+U02θz¯22+zq+z¯q¯+⋯,=U20θz22+U11θzz¯+U02θz¯22+1,q2Tbiω∗θ+1,q¯2Te−iω∗θz¯+⋯.
In view of (52), substituting (53) into (51),(54)gz,z¯=q∗¯F0z,z¯=q∗¯0F0xtθ,=ρ¯q¯1∗,1K1φ1200,(55)g20=2ρ¯q¯∗k1,g11=2ρ¯q¯∗k1,g02=2ρ¯q¯∗k1,g30=6ρ¯q¯∗k1U2010,g21=2ρ¯q¯∗k1U2010+2U1100,g12=2ρ¯q¯∗k1U0210+2U1110,g03=6ρ¯q¯∗k1U0210,where U20k,U11k, and U02k are the k-th components of U20,U11, and U02, respectively, for −τ≤θ≤0. The results in the theorem are valid.Next, we will get the values ofU20θ and U11θ. From (38) and (47), we obtain(56)U˙=x˙t−z˙q−z¯˙q¯,=A0U−2Req¯∗0F0qθ,−τ≤θ<0,A0U−2Req¯∗0F0qθ+F0,θ=0,where(57)Kz,z¯,θ=K20z22+K11zz¯+K02z¯22+⋯.From (48), we have(58)U˙=U˙zz˙t+U˙z¯z¯˙t=U20θz+U11θz¯+⋯iω∗zt+gz,z¯+U11θz+U02θz¯+⋯−iω∗z¯t+gz,z¯.Substituting (48) into (56), we can get that(59)U˙=A0U20θz22+U11θzz¯+U02θz¯22+⋯+K20θz22+K11θzz¯+K02θz¯22+⋯,=A0U20θ+K20θz22+A0U11θ+K11θzz¯+A0U02θ+K02θz¯22+⋯.By comparing the coefficients ofz2 and zz¯ from (58) and (59), we get(60)A0−2iω∗IU20θ=−K20θ,A0U11θ=−K11θ.Forθ∈0,τ, it follows from (47), (48), (57), and (58) that(61)Kz,z¯,θ=−2Req¯∗0F0qθ,=−q¯∗F0qθ−q∗0F¯0q¯θ,=−gz,z¯qθ−g¯z,z¯q¯θ,=−g20z22+g11zz¯+g02z¯22+⋯qθ−g¯20z¯22+g¯11zz¯+g¯02z22+⋯q¯θ.By comparing the coefficients ofz2 and zz¯ from (57) and (61), we get(62)K20θ=−g20qθ−g¯02q¯θ,K11θ=−g11qθ−g¯11q¯θ.Thus, we can get(63)U˙20θ=2iω∗U20θ+g20qθ+g¯02q¯θ.Sinceqθ=q1,q2Teiω∗θ, we obtain(64)U20θ=ig20ω∗q0eiω∗θ+ig¯02ω∗q¯0e−iω∗θ+G1e2iω∗θ,where G1=G11,G12T is a constant vector. We also can have(65)U11θ=−ig11ω∗q0eiω∗θ+ig¯11ω∗q¯0e−iω∗θ+G2,where G2=G21,G22T is a constant vector.From the definition ofA0 and (60), we have(66)∫−τ0dηθU20θ=2iω∗U200−K200,and(67)∫−τ0dηθU11θ=−K110,where ηθ=ηθ,0. From (56), when θ=0, we get(68)Kz,z¯,0=−2Req¯∗0F0q0+F0,=−gz,z¯q0−g¯z,z¯q¯0+F0.Then, we obtain(69)Kz,z¯,0=−2Req¯∗0F0q0+F0,=−gz,z¯q0−g¯z,z¯q¯0+F0.Therefore, we obtain(70)K20z22+K11zz¯+K02z¯22+⋯=−g20z22+g11zz¯+g02z¯22+⋯q0−g¯20z¯22+g¯11zz¯+g¯02z¯22+⋯q¯0+F0.Comparing both sides of (70), we have(71)K20=−g20q0−g¯02q¯0+2k10,(72)K11=−g11q0−g¯11q¯0+22k10.Sinceiω∗ is the eigenvalue of A0 and q0 is the corresponding eigenvector, we get(73)iω∗I−∫−τ0eiω∗θdηθq0=0,−iω∗I−∫−τ0e−iω∗θdηθq¯0=0.Then, substituting (65) and (67) into (72), we obtain(74)2iω∗I−∫−τ0e2iω∗θdηθG1=2k10,that is,(75)K∗G1=2k10,where(76)K∗=2iω∗−a11−DrΛαe−2iω∗τ02iω∗−DsΛα.Thus,(77)G11=2k1e−2iω∗τ02iω∗−DsΛαDetK∗,G12=2iω∗−a11−DsΛα2k100DetK∗.Therefore, we can get the values ofU20θ and U11θ from (64) and (65).



---

*Source: 9931662-2021-07-09.xml*