# The Study for Synchronization between Two Coupled FitzHugh-Nagumo Neurons Based on the Laplace Transform and the Adomian Decomposition Method

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

---

## Abstract

The synchronization between two coupled FitzHugh-Nagumo (FHN) neurons with or without external current is studied by using the Laplace transform and the Adomian decomposition method. Different from other researches, the synchronization error system is expressed as sets of Volterra integral equations based on the convolution theorem in the Laplace transform. Then, it is easy to analytically obtain the conditions that synchronization errors disappear based on the successive approximation method in integral equation theorem, the correctness of which is verified by numerical simulations. Furthermore, the synchronous dynamics of the two coupled FHN neurons also can be written in the form of Volterra integral equations, which is more convenient to analytically solve by using the Adomian decomposition method. It is found that the occurrence of synchronization between the two FHN neurons only depends on the coupling strength and is irrelevant to the external current. Only synchronous rest state in the two FHN neurons without external current can be achieved, while synchronous spikes appear if the external current is not zero.


---

## Body

## 1. Introduction

The neurons are considered as basic units to form neural systems, which have received much interest because of the great significance and applications in neural science, brain science, medical technology, and so on. To simplify matters, the Hodgkin–Huxley (HH) equation [1–3] is usually used to investigate the dynamics of neural networks. Although the HH equation and its modified versions are simplification of dynamics of a realistic neuron, it is still complex in some problems. A simplified model of the HH equation called the FitzHugh-Nagumo (FHN) model [4] is convenient in the investigation of neural behavior by using nonlinear dynamical theory [5, 6]. The FHN model has cubic nonlinearity, which can show excitability of a neuron. The qualitative nature of nerve impulse propagation and neural activity can exhibit in several coupled FHN models, such as separatix loops and bifurcation of equilibria and limit cycles [7, 8].Synchronization in neural systems is a crucial phenomenon. In many regions of the brain, synchronization is considered to be related to cognition as well as the correlate of behavior [3, 9]. Phase synchronization [10] may play a very important role in the brain’s ability to store, process, and communicate information. Synchronization caused by parameter mismatch [11] is a more general case in the real world. Synchronization and corresponding dynamics of coupled neurons have been studied in [12, 13] to understand information processing in the brain. The effects of Gaussian white noise on synchronization were researched in [14].In the past decade, many criteria have been given to judge the occurrence of synchronization in dynamical systems. The Lyapunov function technology is one of the earliest methods to establish the synchronization conditions in coupled oscillators. The famous master stability function (MSF) method is very effective to analyze the local stability of the synchronization manifold [15]. Chen [16] proposed the matrix measure approach (MMA) to analytically derive sufficient conditions for synchronization in time-varying networks. Lu and Chen [17] established the synchronization conditions for linearly coupled ordinary differential systems by considering the distance between trajectories and the synchronization manifold. The key to the Lyapunov function method is to find an appropriate Lyapunov function. Usually, the synchronization conditions by using such technology are sufficient and highly conserved. The study in [18] showed that the method may be invalid for systems with multiple components of different types. In the MSF method and the MMA, the conditional Lyapunov exponents have to be calculated numerically. For the new approach given by Lu and Chen, some information of the synchronous dynamics needs to be known previously. However, the trajectories in the synchronization manifold are unknown before solving the coupled systems. According to the statements above, the analytical criterion for synchronization occurrence is still a question worthy to further discuss.In this paper, we use the Laplace transform to investigate the synchronization conditions for two coupled FHN neurons. 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 Laplace transform. It should be noted that the technique that converts an ordinary differential equation into sets of Volterra integral equations has been successfully applied to other nonlinear systems. From References [19, 20], it was found that the form of the Volterra integral equations is clearer to show the influence of nonlinear parts in the governing equation of a beam under a moving load. According to the successive approximation method [21] in the integral equation theory, it is convenient to obtain the conditions that synchronization errors disappear. Furthermore, the governing equation describing the synchronous dynamics of the two FHN neurons also can be written in the form of the Volterra integral equation, which is easy to analytically solve by using the Adomian decomposition method [22]. Numerical simulations are carried out to demonstrate the correctness of the analytical approach proposed in this paper. It is analytically verified that two coupled FHN neurons without external current can only exhibit synchronization rest regardless of the coupling strength. When the external current is not zero, synchronization spikes between the two FHN neurons occur, and the synchronization conditions for coupling strength remain the same as that for the case of no external current.The rest of the paper is organized as follows. In Section2, the synchronization conditions for two coupled FHN neurons are discussed by using the Laplace transform. In Section 3, the synchronous dynamics of the coupled FHN neurons is analytically solved by using the Adomian decomposition method. In Section 4, numerical simulations are carried out to verify the effectiveness of the analytical results. Conclusions are given in Section 5.


## 2. Synchronization Conditions for Two Coupled FHN Neurons

A single FHN neuron is given in the following form:(1)u̇1=−u1u1−1u1−a−u2+Isinωt,(2)u̇2=bu1−γu2,where u1 is the potential difference and u2 is a recovery current. I and ω are the magnitude and frequency of the external current, respectively. 0<a<0.5, b,γ>0, 1−a+a2/3>bγ, and 1−a+a2/3<1/γ. 0,0 is the unique equilibrium point needed to be considered in system (1). Two FHN neurons are coupled as follows:
(3)u̇1=−u1u1−1u1−a−u2+cu1−u3+Isinωt,(4)u̇2=bu1−γu2,(5)u̇3=−u3u3−1u3−a−u4+cu3−u1+Isinωt,(6)u̇4=bu3−γu4,where c is the coupling strength. The two FHN neurons are said to be synchronous if limt→+∞u1=u3 and limt→+∞u3=u4. By letting
(7)e1=u1−u32,e2=u2−u42,e3=u1+u32,e4=u2+u42,system (3) can be written as
(8)e·1=−e13+−3e32+2a+1e3−a+2ce1−e2,(9)ė2=be1−γe2,(10)ė3=−e33+a+1e32−3e12+ae3+a+1e12−e4+Isinωt,(11)ė4=be3−γe4.Then, the synchronization conditions of system (3) become limt→+∞e1,2=0. In the following, we use the Laplace transform method to obtain the condition that ei→0,i=1,2. The Laplace transform used in this paper is defined as follows:
(12)e^is=Lei=∫0+∞eite−stdt,eit=L−1e^i=12πi∫σ−i∞σ+i∞e^isestds,i=1,2,3,4.Taking the Laplace transforms of both sides of Equation (9) and arranging them yield
(13)s+a−2ce^1+e^2=e10+F^1,(14)s+bγe^2+e20=+be^1,(15)s+ae^3+e^4=e30+ωIs2+ω2+F^2,(16)s+bγe^4=e40+be^3,where ei0,i=1,2,3,4, are given initial values of system (9), and
(17)F^1=∫0+∞−e13−3e1e32+2a+1e1e3e−stdt≡∫0+∞F1te−stdt,F^2=∫0+∞−e33+a+1e32−3e12e3+a+1e12e−stdt≡∫0+∞F2te−stdt.Substituting the second and fourth equations into the first and third equations in system (14), respectively, produces
(18)e^1=s+α1e10−e20s2+a2s+a3+s+α1F^1s2+a2s+a3,(19)e^3=s+α1e30−e40s2+a4s+a5+s+α1ωI/s2+ω2s2+a4s+a5+s+a1F^2s2+a4s+a5,where α1=bγ, α2=a+bγ−2c, α3=ba−2cγ+1, α4=a+bγ, and α5=baγ+1. Since there exist the simply linear relations between e^2,4 and e^1,3, respectively, only the dynamical behaviors of e^1,3 need to be considered. Let Φit,i=1,2,3,4, denote the following inverse Laplace transforms, respectively,
(20)Φ1t=L−1s+a1s2+a2s+a3,Φ2t=L−11s2+a2s+a3,Φ3t=L−1s+a1s2+a4s+a5,Φ4t=L−11s2+a4s+a5.Then, Φit,i=1,2,3,4, can be given by
(21)Φ1t=e−α2t/2coshβ1t+β22β1sinhβ1t,β1≠0,e−α2t/21+2β2t,β1=0,Φ2t=e−α2t/2β3sinhβ1t,β1≠0,te−α2t/2,β1=0,(22)Φ3t=e−α4t/2coshβ3t+β42β3sinhβ3t,β3≠0,e−α4t/21+2β4t,β3=0,Φ4t=e−α4t/2β3sinhβ3t,β3≠0,te−α4t/2,β3=0,where β1=α22−4a3/2,β2=2a1−a2,β3=α42−4a5/2,β4=2a1−a4. From the convolution theorem, taking the inverse Laplace transform on both sides of equations in Equation (18) yields
(23)e1=e10Φ1−e20Φ2+∫0tΦ1t−τF1τdτ,(24)e3=e30Φ3−e40Φ4+I∫0tΦ3t−τsinωτdτ+∫0tΦ3t−τF2τdτ,whereF1,2 is defined in Equation (14). From Equations (21) and (22), if −a2/2±Reβ1<0, Reβ1 represents the real part of β1,Φ1,2 decay to zero quickly. In this case, the first equation in Equation (23) is simplified as
(25)e1t=∫otΦ1t−τF1τdτ=∫0tΦ1t−τ−e1τ2−3e3τ2+2a+1e3τe1τdτ.From the successive approximation method introduced in [21, 23], the recurrence relation
(26)e1n+1t=∫0tΦ1t−τ−e1nτ2−3e3τ2+2a+1e3τe1nτdτ,n≥0can be used to get the solution of Equation (25). Starting with an initial guess e10=0, it is clear that e1=0 is the solution of Equation (25) for any e3.If −a2/2±Reβ1≥0,Φ1,2 cannot approach to zero, then it means that e1 does not decay to zero. Thus, −a2/2±Reβ1<0 is a necessary and sufficient condition for e1→0 in system (9). It is easy to check that
(27)sign−α22+Reβ1=+,a2<0,sign−a3,a2>0&a22−4a3≥0,−,a2>0&a22−4a3<0,(28)sign−α22−Reβ1=signa3,a2<0&a22−4a3≥0,+,a2<0&a22−4a3<0,−,a2>0.Clearly, α2=a+bγ−2c>0 must be satisfied to guarantee −a2/2±Reβ1<0. Under such condition, α3=ba−2cγ+1>b1−bγ2>0 holds for the parameter ranges given in Equation (1). Therefore, α2>0c<a+bγ/2 is the necessary and sufficient condition for e1→0 in system (9). Furthermore, if α2>0 and a22−4a3<0,e1 approaches to 0 with oscillatory when t→+∞. If α2>0 and a22−4a3≥0,e1 will lose its oscillatory behavior when itconverges to zero.Lemma 1.
For two coupled FHN neurons in system (1) with or without external current, synchronization occurs when and only when c<a+bγ/2.Remark 2.
Only whenΦ1,2 defined in Equation (21) decay to 0, e1→0 for t→+∞. To guarantee Φ1,2 approach to 0, it is necessary and sufficient that the denominator in the first equation in Equation (18)
(29)s2+α2s+α3has two roots with negative real parts. According to the Routh-Hurwitz criterion, α2,3>0 are necessary and sufficient to guarantee that all roots of Equation (29) have negative real parts. Moreover, if polynomial (29) has two negative real roots, e1t is nonoscillatory.


## 3. Synchronous Dynamics of the Coupled FHN Neurons

In this section, we discuss the synchronous dynamics of system (3) with I=0 and I≠0, respectively.


### 3.1.I=0 in System (3)

In this case, from the second equation in Equation (23), the governing equation of synchronous dynamics is given by
(30)e3=e30Φ3−e40Φ4+∫0tΦ3t−τF2τdτ=e30Φ3−e40Φ4+∫0tΦ3t−τ−e3τ+a+1e3τ2dτ.Since α4=a+bγ>0 and −a4/2±Reβ3<0 hold according to the parameter ranges given in system (1), Φ3,4 defined in Equation (22) decay to zero very quickly when t→+∞. Using the successive approximation method, similar as that for Equation (25), the solution of Equation (30) is e3=0. Thus, if there is no external current in the coupled system, two FHN neurons only can achieve the synchronization of the resting state.


### 3.2.I≠0 in System (3)

Under such condition, by getting rid of the terms includinge−α4t in the second equation in Equation (23), the governing equation of synchronous dynamics can be simplified as
(31)e3=ϕ0eiωt+ϕ¯0e−iωt+∫0tΦ3t−τF2τdτ=ϕ0eiωt+ϕ¯0e−iωt+∫0tΦ3t−τ−e3τ+a+1e3τ2dτ,where
(32)β3≠0,ϕ0=I−2ω4ω2+a42+2a4β4+4β32+i−a43+4β32−ω2a4+−a42+4ω2+4β32β44ω2+a4−2β324ω2+a4+2β32,β3=0,ϕ0=I−2ω4ω2+8a4β4+a42+i44ω2−a42β4−a44ω2−a424ω+a422.ϕ¯0 is the complex conjugate of ϕ0, and Φ3 is defined in Equation (22).According to the Adomian decomposition method [22], we let
(33)e3=∑n=0∞yn,(34)F2=−e33+a+1e32=∑n=0∞An,where
(35)An=1n!dndλnF2∑n=0∞λnynλ=0,n=0,1,2,⋯,are Adomian polynomials of ordern, in terms of y0,y1,⋯,yn,⋯. The first several Adomian polynomials are given as
(36)A0=a+1y02−y03,A1=2a+1y0y1−3y02y1,A2=y02a+1y2−3y12−3y02y2+a+1y12,⋯,⋯,Substituting Equations (33) and (34) into Equation (31) yields
(37)∑n=0∞yn=ϕ0eiωt+ϕ¯0e−iωt+∫0tΦ3t−τ∑n=0∞Andτ.Then, the solution structure of Equation (37) can be written as
(38)y0=ϕ0eiωt+ϕ¯0e−iωt,y1=∫0tΦ3t−τa+1y02−y03dτ,y2=∫0tΦ3t−τ2a+1y0y1−3y02y1dτ,⋯,⋯,⋯.In the calculations ofy1,2, the terms including e−α4t can be discarded to simplify the expressions. For β3≠0,y1,2 is easily given as
(39)y1=w0+∑n=13wneinω+wn¯e−inω,wn=f1hn,nω,n=0,1,2,3,y2=f1q0,0+∑n=15f1qn,nωeinω+f1qn,nω¯e−−inω,where
(40)f1x,y=2xi2y+β4+a4i2β3−a4+2yi2β3+a4−2y,h0=2a+1ϕ02,h1=−3ϕ02ϕ0,h2=a+1ϕ02,h3=−ϕ03,q0=2a+1ϕ¯0w1+w1¯ϕ0−3ϕ02¯w2−6w0ϕ02−3w2¯ϕ02,q1=2a+1ϕ¯0w2+ϕ0w0−6w1ϕ02−3ϕ02¯w3−3w1¯ϕ02,q2=2a+1ϕ¯0w3+ϕ0w1−6w2ϕ02−3ϕ02w0¯,q3=2a+1ϕ0w2−6w3ϕ02−3ϕ02w1¯,q4=2a+1ϕ0w3−3ϕ02w2,q5=−3ϕ02w3.Then, the approximate solution of Equation (31) can be expressed as
(41)e3=y0+y1+y2.The approximate solution for the case ofβ3=0 can be derived in a similar way, which is not presented here. The validity of the above analysis results will be demonstrated in Section 4.


## 3.1.I=0 in System (3)

In this case, from the second equation in Equation (23), the governing equation of synchronous dynamics is given by
(30)e3=e30Φ3−e40Φ4+∫0tΦ3t−τF2τdτ=e30Φ3−e40Φ4+∫0tΦ3t−τ−e3τ+a+1e3τ2dτ.Since α4=a+bγ>0 and −a4/2±Reβ3<0 hold according to the parameter ranges given in system (1), Φ3,4 defined in Equation (22) decay to zero very quickly when t→+∞. Using the successive approximation method, similar as that for Equation (25), the solution of Equation (30) is e3=0. Thus, if there is no external current in the coupled system, two FHN neurons only can achieve the synchronization of the resting state.


## 3.2.I≠0 in System (3)

Under such condition, by getting rid of the terms includinge−α4t in the second equation in Equation (23), the governing equation of synchronous dynamics can be simplified as
(31)e3=ϕ0eiωt+ϕ¯0e−iωt+∫0tΦ3t−τF2τdτ=ϕ0eiωt+ϕ¯0e−iωt+∫0tΦ3t−τ−e3τ+a+1e3τ2dτ,where
(32)β3≠0,ϕ0=I−2ω4ω2+a42+2a4β4+4β32+i−a43+4β32−ω2a4+−a42+4ω2+4β32β44ω2+a4−2β324ω2+a4+2β32,β3=0,ϕ0=I−2ω4ω2+8a4β4+a42+i44ω2−a42β4−a44ω2−a424ω+a422.ϕ¯0 is the complex conjugate of ϕ0, and Φ3 is defined in Equation (22).According to the Adomian decomposition method [22], we let
(33)e3=∑n=0∞yn,(34)F2=−e33+a+1e32=∑n=0∞An,where
(35)An=1n!dndλnF2∑n=0∞λnynλ=0,n=0,1,2,⋯,are Adomian polynomials of ordern, in terms of y0,y1,⋯,yn,⋯. The first several Adomian polynomials are given as
(36)A0=a+1y02−y03,A1=2a+1y0y1−3y02y1,A2=y02a+1y2−3y12−3y02y2+a+1y12,⋯,⋯,Substituting Equations (33) and (34) into Equation (31) yields
(37)∑n=0∞yn=ϕ0eiωt+ϕ¯0e−iωt+∫0tΦ3t−τ∑n=0∞Andτ.Then, the solution structure of Equation (37) can be written as
(38)y0=ϕ0eiωt+ϕ¯0e−iωt,y1=∫0tΦ3t−τa+1y02−y03dτ,y2=∫0tΦ3t−τ2a+1y0y1−3y02y1dτ,⋯,⋯,⋯.In the calculations ofy1,2, the terms including e−α4t can be discarded to simplify the expressions. For β3≠0,y1,2 is easily given as
(39)y1=w0+∑n=13wneinω+wn¯e−inω,wn=f1hn,nω,n=0,1,2,3,y2=f1q0,0+∑n=15f1qn,nωeinω+f1qn,nω¯e−−inω,where
(40)f1x,y=2xi2y+β4+a4i2β3−a4+2yi2β3+a4−2y,h0=2a+1ϕ02,h1=−3ϕ02ϕ0,h2=a+1ϕ02,h3=−ϕ03,q0=2a+1ϕ¯0w1+w1¯ϕ0−3ϕ02¯w2−6w0ϕ02−3w2¯ϕ02,q1=2a+1ϕ¯0w2+ϕ0w0−6w1ϕ02−3ϕ02¯w3−3w1¯ϕ02,q2=2a+1ϕ¯0w3+ϕ0w1−6w2ϕ02−3ϕ02w0¯,q3=2a+1ϕ0w2−6w3ϕ02−3ϕ02w1¯,q4=2a+1ϕ0w3−3ϕ02w2,q5=−3ϕ02w3.Then, the approximate solution of Equation (31) can be expressed as
(41)e3=y0+y1+y2.The approximate solution for the case ofβ3=0 can be derived in a similar way, which is not presented here. The validity of the above analysis results will be demonstrated in Section 4.


## 4. Numerical Simulations

In this section, some numerical simulations for Equation (3) with a=0.1, b=0.08, and γ=3 are carried out to demonstrate the effectiveness of the analysis in the last section. For the case of I=0, the synchronization condition is c<a+bγ/2=0.1+3×0.08/2=0.17 according to the analysis results in Equation (28). To illustrate the validity of the critical value of the coupling strength, c=0.16 and c=0.18 are chosen to carry out the numerical simulations for Equation (3) with I=0. The initial conditions are taken as u10=0.01, u20=−0.01, u30=0.02, and u40=0.01. The results are presented in Figure 1, which show that synchronous spikes do not occur for any given coupling strength in two coupled FHN neurons without external current.Figure 1
The dynamics ofu1,2,3,4 in Equation (3) with I=0. (a) The dynamics of u1,3 with c=0.16. (b) The dynamics of u2,4 with c=0.16. (c) The dynamics of u1,3 with c=0.18. (d) The dynamics of u2,4 with c=0.18. The initial conditions are taken as u10=0.01, u20=−0.01, u30=0.02, and u40=0.01.
(a)
Dynamics ofu1,3 for c=0.16(b)
Dynamics ofu2,4 for c=0.16(c)
Dynamics ofu1,3 for c=0.18(d)
Dynamics ofu2,4 for c=0.18For the case ofI≠0, the current and frequency are chosen as I=0.01 and ω=0.1, respectively. The values of a, b, and γ remain the same. The initial conditions are taken as u10=0.01, u20=u30=u40=0. The dynamics of u1,2,3,4 and synchronization errors u1−u3, u2−u4 are depicted in Figure 2 for c=0.16 and c=0.18, which demonstrate the effectiveness of the synchronization criterion c<a+bγ/2.Figure 2
The dynamics ofu1,2,3,4 and synchronization errors in Equation (3) with I=0.01, ω=0.1. (a) The dynamics of u1,3 and the error between them for c=0.16. (b) The dynamics of u2,4 and the error between them for c=0.16. (c) The dynamics of u1,3 and the error between them for c=0.18. (d) The dynamics of u2,4 and the error between them for c=0.18. The initial conditions are taken as u10=0.01, u20=u30=u40=0.
(a)
u1,3 and their error for c=0.16(b)
u2,4 and their error for c=0.16(c)
u1,3 and their error for c=0.18(d)
u2,4 and their error for c=0.18Next, we verify the effectiveness of analytical results of Equations (38)–(41) describing the synchronous dynamics of the two FHN neurons derived by using the Adomian decomposition method. The parameters in Equation (3) are given by a=0.1, b=0.08, γ=3, I=0.01, ω=0.1, and c=0.1, according to Equations (38)–(41); the approximation of e3 is expressed as
(42)e3=−0.0000891cos0.1t+0.0003954cos0.274t+−0.0001033cos0.1t+0.00003sin0.1t+0.003753sin0.274te−0.17t+0.000853+0.001248cos0.1t−0.001233cos0.2t+0.00005934cos0.3t+0.02617sin0.1t−0.00003028sin0.2t−0.0001005sin0.3t.The comparison between the numerical results of Equation (3) using the Runge-Kutta method for different initial conditions and the analytical results from Equation (42) is given in Figure 3. Clearly, the analytical results coincide exactly with the numerical ones after enough time, which show the correctness of the analytical results expressed in Equation (42).Figure 3
The comparison between numerical and analytical results of synchronous dynamicse3=u1+u3/2 in Equation (3) with I=0.01, ω=0.1 for different initial conditions. (a) The synchronization dynamics of e3 and u1+u3/2. (b) The errors between e3 and u1+u3/2. The initial conditions are u20=u30=u40=0, while u10 varies.
(a)
Dynamics ofe3 and u1+u3/2(b)
Errors betweene3 and u1+u3/2


## 5. Conclusions

In this paper, the synchronization of two FHN neurons with or without external current is investigated by using the Laplace transform and the Adomian decomposition method. A linear transformation is carried out for the original system to obtain the synchronization error system. Different from other methods, the error system is converted into sets of Volterra integral equations by using the convolution theorem. Then, the successive approximation method in the integral equation theorem is used to judge whether the synchronization error approaches to zero. Thus, a criterion is obtained to determine the occurrence of synchronization in the two FHN neurons.It is found that the two coupled FHN neurons without external current can only maintain neural rest when synchronization is achieved. While the external current is not zero, synchronous spikes occur. The numerical simulations demonstrate that the criterion presented in this paper is very effective. The analytical result describing the synchronous dynamics of the two FHN neurons derived by using the Adomian decomposition method is accurate enough. The analytical expression helps us to deeply understand the synchronous dynamics of coupled FHN neurons. The method adopted in this paper is valid to discuss the synchronization between FHN neurons. Moreover, the calculation quantities of this method are small, which can be developed to judge synchronization betweenN coupled FHN neurons.



---

*Source: 6657835-2021-04-24.xml*