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

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.


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][9][10][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][17][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 neu-rons 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 Section 2, 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.

Two HR Neurons with a Linear Coupling Function
A HR neuron is described by the following equation of motion: 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, s 0 , r, and x 0 are parameters, and I is the external current input. Consider two coupled HR neurons with a linear coupling function: where β is the coupling strength.
, synchronization between two HR neurons is achieved. Assume that ðx 01 , y 01 , z 01 Þ is the equilibrium point of system (1). For simplicity, substituting Neural Plasticity into equation (2), one has where ϕ 1 = −3ax 2 01 + 2bx 01 , ϕ 2 = −3ax 01 + b, and ϕ 3 = −2d x 01 . By letting system (4) can be written as where If e 1,2,3 → 0 for t → ∞, the two HR neurons will achieve synchronization. Consider the Laplace transform defined as follows: Taking the Laplace transforms of both sides of equation (6) and arranging them yield where e i0 , i = 1, ⋯, 6, are given initial values of system (6), andF Since only e 1,2,3 need to be considered, solving the first three equations in system (9) leads tô where From the convolution theorem, taking the inverse Laplace transform on both hands of equations in equation (11), one has 3 Neural Plasticity where Φ i ðtÞ, i = 1, 2, ⋯, 9, denote the following inverse Laplace transforms, respectively, The kernels in Φ i ðtÞ, i = 1, 2, ⋯, 9, are true fractions which can be partitioned into the partial fraction expansions; therefore, where ξ i are undetermined constants and s = s i , i = 1, 2, 3, are the roots to the equation g 1 ðsÞ = 0 (g 1 ðsÞ is defined in equation (11)). If all the real parts of roots of g 1 ðsÞ = 0 are negative, Φ 1,⋯,9 ðtÞ → 0 when the time approaches to infinity. In this case, equation (13) becomes From the successive approximation method introduced in references [26,27], e 1,2,3 → 0 in equation (16) with time approaching to infinity if all the roots of g 1 ðsÞ = 0 have negative real parts. g 1 ðsÞ has the form of According to the Routh-Hurwitz criterion, 3 Hurwitz matrices using the coefficients p i , i = 1, 2, 3, of g 1 ðsÞ are given by Since det H 1 = p 1 > 0, the necessary and sufficient condition that all the roots of g 1 ðsÞ = 0 have negative real parts can be written as 4 Neural Plasticity where ϕ 1,3 are defined in equation (4). Equation (21) just is the synchronization condition for the two coupled HR neurons in system (2).
The numerical results are presented in Figure 1, which demonstrates the effectiveness of the synchronization condition (21). When r = 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).

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.

Data Availability
No data were used to support this study.

Disclosure
The funding sponsors had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, and in the decision to publish the results.

Conflicts of Interest
The authors declare no conflict of interest.