# A Synchronization Criterion for Two Hindmarsh-Rose Neurons with Linear and Nonlinear Coupling Functions Based on the Laplace Transform Method

Authors: Chunlin Su; Bin Zhen; Zigen Song
Journal: Neural Plasticity
 (2021)

---

## Abstract

In this paper, an analytical criterion is proposed to investigate the synchronization between two Hindmarsh-Rose neurons with linear and nonlinear coupling functions based on the Laplace transform method. Different from previous works, the synchronization error system is expressed in its integral form, which is more convenient to analyze. The synchronization problem of two HR coupled neurons is ultimately converted into the stability problem of roots to a nonlinear algebraic equation. Then, an analytical criterion for synchronization between the two HR neurons can be given by using the Routh-Hurwitz criterion. Numerical simulations show that the synchronization criterion derived in this paper is valid, regardless of the periodic spikes or burst-spike chaotic behavior of the two HR neurons. Furthermore, the analytical results have almost the same accuracy as the conditional Lyapunov method. In addition, the calculation quantities always are small no matter the linear and nonlinear coupling functions, which show that the approach presented in this paper is easy to be developed to study synchronization between a large number of HR neurons.


---

## Body

## 1. Introduction

The study of neurons plays a very important role in many applications in neural science, brain science, and so on. The famous Hodgkin-Huxley (HH) equation, proposed by Hodgkin and Huxley in 1952, [1, 2] usually was used to construct neural systems or exhibit the neural dynamic behavior. Since then, other neuronal models, such as the FizHugh-Nagumo (FHN) model [3], the Hindmarsh-Rose (HR) model [4], the Chay model [5], and the Morris-Lecar model [6], also have been published. The HR model was constructed from voltage clamp data to provide a simple description of the patterned activity observed in molluscan neurons. Although the HR model is not based on physiology wholly, it can exhibit some features observed in neuronal bursting. The HR model has been investigated from bifurcation analysis to the firing mechanism [7].Synchronization processes are ubiquitous in nature. Synchronization of coupled networks has been studied with increasing interest over the last few decades due to its numerous potential applications [8–11]. One way to gain a deeper understanding of synchronization in complex networks is to investigate the stability of the completely synchronized state of the population of identical oscillators [12, 13], which has been extensively investigated based on Lyapunov’s direct method [14, 15]. One of the most popular and widely used methods to determine stability of the synchronized state is using the Lyapunov exponents as average measurements of expansion or shrinkage of small displacements along synchronized trajectory, which is called the conditional Lyapunov exponent method [16–18].Pecora and Carroll [19] proposed the master stability function (MSF) method to discuss the local stability of the synchronization manifold. This method provides an objective criterion to characterize the stability of the global synchronized state, independent of particularities of oscillators. Relevant insights about the structure-dynamics relationship have been obtained using this technique. Chen [20] and Li et al. [21] obtained conditions for synchronization using the matrix measure approach. Their criteria do not depend on solutions of the synchronous state and have fewer limitations on network connections. He [22] presented a new approach for chaotic synchronization of HR neural networks with a special nonlinear coupling function. Synchronization can be achieved without the requirement to calculate the maximum Lyapunov exponents when the coupling strengths are taken as reference values. Lu and Chen [23] used a variation method near the projection on the synchronization manifold to analyze synchronization of linearly coupled ordinary differential equations. The uncoupled dynamical behavior at each node is general and can be chaotic or otherwise; the coupling configuration is also general, with the coupling matrix not assumed to be symmetric or irreducible.The above methods greatly further one’s understanding of synchronization in complex networks. However, the results obtained by these methods are limited. Conditions for synchronization obtained by using the Lyapunov function method are often too conservative. It is also difficult to find a proper Lyapunov function. For most cases, calculating the conditional Lyapunov exponents is a cumbersome procedure. It should also be noted that the negativity of the conditional Lyapunov exponents is only a necessary condition for the stability of the synchronized state. Thus, additional conditions should be added to ensure synchronization in a necessary and sufficient way [24]. The MSF method [19] proposed by Pecora and Carroll requires that the Laplacian matrices are diagonal or block diagonal. Nishikawa and Motter [25] developed an extension of the master stability framework to study local synchronization in any symmetric and asymmetric network with any output function. The method has some variations, but generally the stability analysis is embarrassed by the need to calculate the conditional Lyapunov exponents. Using the matrix measure approach requires calculating the Lyapunov exponents of the map to achieve a sufficient condition for synchronization. The new method proposed by He [22] is only valid for a special nonlinear coupling function and coupling strengths. The local stability condition presented in [23] for synchronization is related to the synchronization trajectory. Since the synchronization trajectory is unknown in advance, the criterion is hard to be used in applications. To sum up, in addition to Lyapunov’s direct method, nearly all synchronization criteria require the computation of Lyapunov exponents, eigenvalues of the coupling matrix, or the synchronization trajectory, which is inconvenient and ineffective. Therefore, more convenient and effective approaches to determine synchronization need to be developed.In this paper, we use the Laplace transform to investigate the synchronization conditions for two HR neurons with linear and nonlinear coupling functions. Different from other studies, the synchronization error system derived from the original coupled neural system is converted into sets of Volterra integral equations based on the convolution theorem in the Laplace transform. The synchronization problem of two HR coupled neurons is ultimately converted into the stability problem of roots to a nonlinear algebraic equation. Then, an analytical criterion for synchronization in the two HR neurons can be given by using the Routh-Hurwitz criterion. Numerical simulations show that the synchronization criterion derived in this paper is valid, regardless of the periodic spikes or burst-spike chaotic behavior of the two HR neurons. Additionally, the analytical results have almost the same accuracy as the conditional Lyapunov method. The rest of the paper is organized as follows. In Section2, the synchronization conditions for two HR neurons with a linear coupling function are discussed by using the Laplace transform. In Section 3, the case of a nonlinear coupling function is investigated. In Section 4, numerical simulations are carried out to verify the effectiveness of the analytical criterion derived in Sections 2 and 3. Conclusions are drawn in Section 5.


## 2. Two HR Neurons with a Linear Coupling Function

A HR neuron is described by the following equation of motion:(1)ẋ=y−ax3+bx2−z+I,ẏ=c−dx2−y,ż=rs0x−x0−rz,where x represents the membrane potential, y is a recovery variable associated with fast current, z is a slowly changing adaptation current. a, b, c, d, s0, r, and x0 are parameters, and I is the external current input. Consider two coupled HR neurons with a linear coupling function:
(2)ẋ1=y1−ax13+bx12−z1+I+βx1−x2,ẏ1=c−dx12−y1,ż1=rs0x1−x0−rz1,ẋ2=y2−ax23+bx22−z2+I+βx2−x1,ẏ2=c−dx22−y2,ż2=rs0x2−x0−rz2,where β is the coupling strength. If ∣x1−x2∣→0, ∣y1−y2∣→0, and ∣z1−z2∣→0 for t→∞, synchronization between two HR neurons is achieved. Assume that x01,y01,z01 is the equilibrium point of system (1). For simplicity, substituting
(3)u1=x1−x01,u2=y1−y01,u3=z1−z01,v1=x2−x01,v2=y2−y01,v3=z2−z01,into equation (2), one has
(4)u̇1=ϕ1u1+ϕ2u12−au13+u2−u3+βu1−v1,u̇2=ϕ3u1−du12−u2,u̇3=rs0u1−u3,v̇1=ϕ1v1+ϕ2v12−av13+v2−v3+βv1−u1,v̇2=ϕ3v1−dv12−v2,v̇3=rs0v1−v3,where ϕ1=−3ax012+2bx01, ϕ2=−3ax01+b, and ϕ3=−2dx01. By letting
(5)e1=u1−v12,e2=u2−v22,e3=u3−v32,e4=u1+v12,e5=u2+v22,e6=u3+v32,system (4) can be written as
(6)ė1=2β+ϕ1e1+e2−e3+F1,ė2=ϕ3e1−e2+F2,ė3=rs0e1−e3,ė4=ϕ1e4+e5−e6+F3,ė5=ϕ3e4−e5+F4,ė6=rs0e4−e6,where
(7)F1=−a3e1e42+e13+2ϕ2e1e4,F2=−2de1e4,F3=−a3e12e4+e43+ϕ2e12+e42,F4=−de12+e42.Ife1,2,3→0 for t→∞, the two HR neurons will achieve synchronization. Consider the Laplace transform defined as follows:
(8)e^i=Lei=∫0+∞eie−stdt,ei=L−1e^i=12πi∫σ−i∞σ+i∞e^isestds,i=1,⋯,6.Taking the Laplace transforms of both sides of equation (6) and arranging them yield
(9)s−2β−ϕ1e^1−e^2+e^3=e10+F^1,−ϕ3e^1+s+1e^2=e20+F^2,−rs0e^1+s+re^3=e30,s+ϕ1e^4−e^5+e^6=e40+F^3,−ϕ3e4+s+1e^5=e50+F^4,−rs0e^4+s+re^6=e60,where ei0, i=1,⋯,6, are given initial values of system (6), and
(10)F^j=∫0+∞Fjte−stdt,j=1,2,3,4.Since onlye1,2,3 need to be considered, solving the first three equations in system (9) leads to
(11)e^1=g2g1+s2+r+1s+rF^1g1+s+rF^2g1,e^2=g3g1+r+sϕ3F^1g1+s2+r−2β−ϕ1s+rs0−2β−ϕ1F^2g1,e^3=g4g1+rs0s+1F^1g1+rs0F^2g1,where
(12)g1=s3−2β−r+ϕ1−1s2−2r+1β+rϕ1−1−s0+ϕ1+ϕ3s−r2β−s0+ϕ1+ϕ3,g2=e10s2+e10r+1+e20−e30s+re10+e20−e30,g3=e20s2+ϕ3e10+r−2β−ϕ1e20s+e10rϕ3+rs0−2β−ϕ1e20−ϕ3e30,g4=e30s2+rs0e10+1−2β−ϕ1e30s+rs0e10+e20−2β+ϕ1+ϕ3e30.From the convolution theorem, taking the inverse Laplace transform on both hands of equations in equation (11), one has
(13)e1=Φ1+∫0tΦ2t−τF1τ+Φ3t−τF2τdτ,e2=Φ4+∫0tΦ5t−τF1τ+Φ6t−τF2τdτ,e3=Φ7+∫0tΦ8t−τF1τ+Φ9t−τF2τdτ,where Φit, i=1,2,⋯,9, denote the following inverse Laplace transforms, respectively,
(14)Φ1t=L−1g2g1,Φ2t=L−1s2+r+1s+rg1,Φ3t=L−1s+rg1,Φ4t=L−1g3g1,Φ5t=L−1r+sϕ3g1,Φ6t=L−1s2+r−2β−ϕ1s+rs0−2β−ϕ1g1,Φ7t=L−1g4g1,Φ8t=L−1rs0s+1g1,Φ9t=L−1rs0g1.The kernels inΦit, i=1,2,⋯,9, are true fractions which can be partitioned into the partial fraction expansions; therefore,
(15)Φ1,⋯,9=L−1ξ1s−s1+ξ2s−s2+ξ3s−s3=∑i=13ξiesit,s1≠s2≠s3,L−1ξ1s−s12+ξ2s−s1+ξ3s−s3=ξ1t+ξ2es1t+ξ3es3t,s1=s2≠s3,L−1ξ1s−s13+ξ2s−s12+ξ3s−s1=ξ12t2+ξ2t+ξ3es1t,s1=s2=s3,where ξi are undetermined constants and s=si, i=1,2,3, are the roots to the equation g1s=0 (g1s is defined in equation (11)). If all the real parts of roots of g1s=0 are negative, Φ1,⋯,9t→0 when the time approaches to infinity. In this case, equation (13) becomes
(16)e1=∫0tΦ2t−τF1τ+Φ3t−τF2τdτ=∫0tΦ2t−τ−a3e4τ2+e1τ2+2ϕ2e4τ−2dΦ3t−τe4τe1τdτ,e2=∫0tΦ5t−τF1τ+Φ6t−τF2τdτ=∫0tΦ5t−τ−a3e4τ2+e1τ2+2ϕ2e4τ−2dΦ6t−τe4τe1τdτ,e3=∫0tΦ8t−τF1τ+Φ9t−τF2τdτ=∫0tΦ8t−τ−a3e4τ2+e1τ2+2ϕ2e4τ−2dΦ9t−τe4τe1τdτ.From the successive approximation method introduced in references [26, 27], e1,2,3→0 in equation (16) with time approaching to infinity if all the roots of g1s=0 have negative real parts. g1s has the form of
(17)g1s=s3+p1s2+p2s+p3.According to the Routh-Hurwitz criterion, 3 Hurwitz matrices using the coefficientspi, i=1,2,3, of g1s are given by
(18)H1=p1,H2=p110p2,H3=p110p3p2p100p3.All of the roots of the polynomialg1s are negative or have negative real parts if and only if the determinants of all Hurwitz matrices are positive:
(19)detHj>0,j=1,2,3.Since(20)detH1=p1>0,detH2=p1p2>0,detH3=p1p2−p3p3>0,the necessary and sufficient condition that all the roots of g1s=0 have negative real parts can be written as
(21)2β−r+ϕ1−1<0,2r+1β+rϕ1−1−s0+ϕ1+ϕ3<0,r2β−s0+ϕ1+ϕ3<0,2β−r+ϕ1−12r+1β+rϕ1−1−s0+ϕ1+ϕ3+r2β−s0+ϕ1+ϕ3>0,where ϕ1,3 are defined in equation (4). Equation (21) just is the synchronization condition for the two coupled HR neurons in system (2).


## 3. Two HR Neurons with Nonlinear Coupling Functions

Consider two coupled HR neurons with a nonlinear coupling function(22)ẋ1=y1−ax13+bx12−z1+I+αHx1−Hx2,ẏ1=c−dx12−y1+αdx12−dx22,ż1=rs0x1−x0−rz1,ẋ2=y2−ax23+bx22−z2+I+αHx2−Hx1,ẏ2=c−dx22−y2+αdx22−dx12,ż2=rs0x2−x0−rz2,where α is the coupling strength and Hx=ax3−bx2−x. Substituting equation (3) into equation (22) produces
(23)u̇1=ϕ1u1+ϕ2u12−au13+u2−u3+αau13−v13−ϕ2u12−v12−ϕ1+1u1−v1,u̇2=ϕ3u1−du12−u2+αdu12−v12−ϕ3u1−v1,u̇3=rs0u1−u3,v̇1=ϕ1v1+ϕ2v12−av13+v2−v3−αau13−v13−ϕ2u12−v12−ϕ1+1u1−v1,v̇2=ϕ3v1−dv12−v2−αdu12−v12−ϕ3u1−v1,v̇3=rs0v1−v3,where ϕ1,2,3 are defined in equation (4). Substituting equation (5) into equation (23) yields
(24)ė1=−2αϕ1+1+ϕ1e1+e2−e3+F5,ė2=1−2αϕ3e1−e2+F6,ė3=rs0e1−e3,ė4=ϕ1e4+e5−e6+F3,ė5=ϕ3e4−e5+F4,ė6=rs0e4−e6,where
(25)F3=−a3e12e4+e43+ϕ2e12+e42,F4=−de12+e42,F5=2α−1ae13+3ae1e42−2ϕ2e1e4,F6=2d2α−1e1e4.Ife1,2,3→0, synchronization appears in equation (23). Taking the Laplace transforms on both sides of the first three equations in equation (24), one has
(26)s+2αϕ1+1−ϕ1e^1−e^2+e^3=e10+F^5,2α−1ϕ3e1+s+1e^2=e20+F^6,−rs0e^1+s+re^3=e30,where ei0, i=1,2,3, are given initial values of system (24), and
(27)F^j=∫0+∞Fjte−stdt,j=5,6.Solving equation (26) leads to
(28)e^1=h2h1+s2+r+1s+rF^5h1+s+rF^6h1,e^2=h3h1+r+s1−2αϕ3F^5h1−h4F^6h1,e^3=h5h1+rs0s+1F^5h1+rs0F^6h1.where
(29)h1=s3+2αϕ1+1−ϕ1+r+1s2+2α−1ϕ1r+1+2αr+1+ϕ3+rs0+1−ϕ3s+r2α−1ϕ1+ϕ3+s0+2α,h2=e10s2+e10r+1+e20−e30s+re10+e20−e30,h3=−e202s2+e102α−1ϕ3−r+2αϕ1+1−ϕ1e20s−e20r2αϕ1+1−ϕ1+e10r−e302α−1ϕ3−rs0e20,h4=s2+2αϕ1+1−ϕ1+rs+r2αϕ1+1−ϕ1+s0,h5=e30s2+rs0e10+1+2αϕ1+1−ϕ1e30s+rs0e10+e20+2α−1ϕ1+ϕ3+2αe30.Taking the inverse Laplace transform on both sides of the three equations in equation (28), one yields
(30)e1=Ψ1t+∫0tΨ2t−τF5τ+Ψ3t−τF6τdτ,e2=Ψ4t+∫0tΨ5t−τF5τ−Ψ6t−τF6τdτ,e3=Ψ7t+∫0tΨ8t−τF5τ+Ψ9t−τF6τdτ,where Ψ1,⋯,9t represent the following inverse Laplace transforms:
(31)Ψ1t=L−1h2h1,Ψ2t=L−1s2+r+1s+rh1,Ψ3t=L−1s+rh1,Ψ4t=L−1h3h1,Ψ5t=L−1r+s1−2αϕ3h1,Ψ6t=L−1h4h1,Ψ7t=L−1h5h1,Ψ8t=L−1rs0s+1h1,Ψ9t=L−1rs0h1.Furthermore,Ψ1,⋯,9 can be expressed as the following forms:
(32)Ψ1,⋯,9=L−1η1s−s1+η2s−s2+η3s−s3=∑i=13ηiesit,s1≠s2≠s3,L−1η1s−s12+η2s−s1+η3s−s3=η1t+η2es1t+η3es3t,s1=s2≠s3,L−1η1s−s13+η2s−s12+η3s−s1=η12t2+η2t+η3es1t,s1=s2=s3,where ηi are undetermined constants and s=si, i=1,2,3, are the roots to the equation h1s=0 (h1s is defined in equation (28)). From the analysis for the case of a linear coupling function in the previous section, if Ψ1,⋯,9→0 when the time approaches to infinity, e1,2,3→0 in equation (24) will be achieved, which means that all roots of h1s=0 should have negative real parts. Based on the Routh-Hurwitz criterion, the synchronization condition in system (23) can be written as
(33)2αϕ1+1−ϕ1+r+1>0,2α−1ϕ1r+1+2αr+1+ϕ3+rs0+1−ϕ3>0,r2α−1ϕ1+ϕ3+s0+2α>0,2αϕ1+1−ϕ1+r+12α−1ϕ1r+1+2αr+1+ϕ3+rs0+1−ϕ3−r2α−1ϕ1+ϕ3+s0+2α>0.


## 4. Numerical Simulations

In this section, some numerical simulations are performed to illustrate the correctness of the synchronization criteria derived in previous sections. Here,a=1.0, b=3.0, c=1.0, d=5.0, s0=4, and x0=−1.60. First, the case of the linear function, which is system (2), is considered. If r=0.02 and I=3.6, the two HR neurons in system (2) with β=0 exhibit periodic spikes. From equation (21), the synchronization condition for system (2) is β<−0.55. To demonstrate the validity of the coupling strength, β=−0.6 and β=−0.3 are chosen to carry out the numerical simulations of system (2), respectively. The initial conditions are taken as x10=0.2, y10=0.1, z10=0.2, x20=0.1, y20=0.1, and z20=0.1. The numerical results are presented in Figure 1, which demonstrates the effectiveness of the synchronization condition (21).Figure 1
The dynamics of system (2) with r=0.02 and I=3.6: (a) β=−0.3 and (b) β=−0.6. Synchronization errors: (c) β=−0.3 and (d) β=−0.6. The initial conditions are x10=0.2, y10=0.1, z10=0.2, x20=0.1, y20=0.1, and z20=0.1.
(a)
dynamics forβ=−0.3(b)
dynamics forβ=−0.6(c)
errors forβ=−0.3(d)
errors forβ=−0.6Whenr=0.013 and I=3, the two HR neurons in system (2) with β=0 exhibit burst-spike chaotic behavior. According to equation (21), the synchronization condition for system (2) is β<−0.56. To demonstrate the validity of the coupling strength, β=−0.6 and β=−0.3 are taken to perform the numerical simulations of system (2), respectively. The initial conditions were kept the same as those in Figure 1. The results are presented in Figure 2, which illustrates the correctness of the synchronization condition (21).Figure 2
The dynamics of system (2) with r=0.013 and I=3.0: (a) β=−0.3 and (b) β=−0.6. Synchronization errors: (c) β=−0.3 and (d) β=−0.6. The initial conditions are x10=0.2, y10=0.1, z10=0.2, x20=0.1, y20=0.1, and z20=0.1.
(a)
dynamics forβ=−0.3(b)
dynamics forβ=−0.6(c)
errors forβ=−0.3(d)
errors forβ=−0.6For the case of the nonlinear coupling function, letx0=−1.56, r=0.006, I=3.0, and other parameters keep unchanged. According to [22], based on the conditional Lyapunov exponent method, the synchronization condition is approximately given as α∈0.265,0.61. From criterion (33), the synchronization condition is 0.275,0.7. Therefore, the synchronization criterion proposed in this paper is valid.


## 5. Conclusion

In this paper, the synchronization between two HR neurons with linear and nonlinear coupling functions is investigated by using the Laplace transform and the convolution theorem. Different from other researchers, the synchronization error system can be expressed in its integral form, which is more convenient to analyze. The synchronization problem is ultimately converted into the stability problem of roots to a nonlinear algebraic equation. Then, an analytical criterion for synchronization in the two HR neurons can be given by using the Routh-Hurwitz criterion. Numerical simulations show that the synchronization criterion derived in this paper is effective, regardless of the periodic spikes or burst-spike chaotic behavior of the two HR neurons. Furthermore, our analytical results have almost the same accuracy as the conditional Lyapunov method. Since the calculation quantities are very small, the approach presented in this paper is easy to be developed to study synchronization between a large number of HR neurons.



---

*Source: 6692132-2021-02-02.xml*