Chaotic Synchronization of memristive neurons: Lyapunov function versus Hamilton function

We study the dynamical behaviors of this improved memristive neuron model by changing external harmonic current and the magnetic gain parameters. The model shows rich dynamics including periodic and chaotic spiking and bursting, and remarkably, chaotic super-bursting, which has greater information encoding potentials than a standard bursting activity. Based on Krasovskii-Lyapunov stability theory, the sufficient conditions (on the synaptic strengths and magnetic gain parameters) for the chaotic synchronization of the improved model are obtained. Based on Helmholtz's theorem, the Hamilton function of the corresponding error dynamical system is also obtained. It is shown that the time variation of this Hamilton function along trajectories can play the role of the time variation of the Lyapunov function - in determining the asymptotic stability of the synchronization manifold. Numerical computations indicate that as the synaptic strengths and the magnetic gain parameters change, the time variation of the Hamilton function is always non-zero (i.e., a relatively large positive or negative value) only when the time variation of the Lyapunov function is positive, and zero (or vanishingly small) only when the time variation of the Lyapunov function is also zero. This clearly therefore paves an alternative way to determine the asymptotic stability of synchronization manifolds, and can be particularly useful for systems whose Lyapunov function is difficult to construct, but whose Hamilton function corresponding to the dynamic error system is easier to calculate.

biological neurons, the original 2D HR neuron model has undergone few modifications. The main objective of Hindmarsh and Rose in 1984 (after the formulation of their 2-equations model [3]) was to model the synchronization of firing of two snail neurons in a simple way, without the use of the HH equations [3,8]. Hence, with the goal to design a neuron model that exhibits triggered firing, some modifications were done on the 2-equations model (by adding an adaptation variable, representing the slowly varying current, that changed the applied current to an effective applied one) to obtain the 3-equations model [8,9]. This 3D model has been very popular in studying biological properties of spiking and bursting neurons, including their chaotic dynamics. Over the last decades, some detailed investigations and studies of bifurcations and the dynamics of the 3D HR model has be done [10,11,12], showing many behaviors observed of real biological neurons.
However, not only the 3D HR model fails to take into account the dynamics of calcium ions across the membrane, it can only capture a relatively small domain of the chaotic regimes of real biological neuron. A few years after the 3D HR model was proposed, Selverston et al. [13] studied a computational and electronic model of stomatogastric ganglion neurons. They used the standard 3D HR model, and discovered that in spite of the fact that the 3D model can produce several modes of spiking-bursting behaviors seen in biological neurons, its parameter space for chaotic activity is much more limited than the one observed in real neurons. For this reason, they modified the 3D HR model, by adding a fourth variable (a slower process) representing the calcium dynamics [13]. The complexity of the 3D model increased and it was then capable of mimicking the complex dynamical (spiking, bursting and chaotic) behavior of pyloric central pattern generator neurons of the lobster stomatogastric system [13,14]. Thus, 4D HR model could capture more complex behavior than the 3D model [13,15] by simply taking into account the calcium exchange between intracellular warehouse and the cytoplasm, and in particular to completely produce the chaotic behavior of the stomatogastric ganglion neurons that the 3D model cannot. Some detailed investigations and studies of bifurcations and the dynamics of the 4D HR model has be done [16,17,18,19].
Before the experimental confirmation of the nonnegligible effects of the magnetic flux (generated by the flow of ions across membrane) on the action potential of neurons, all previous studies on the dynamics of neuron models (including the HH, FHN, Izhikevich, 3D and 4D HR models), have been done without taking into account magnetic flux effects. Recently, M. Lv et al. [20] proposed a modified HR model which takes into account the effect of magnetic field by adding an additional variable for the magnetic flux into the 3D HR model. This modified model not only can generate a variety of modes in electric activities by changing the external forcing current and the magnetic flux parameters, but could also be useful in the investigation of the effect of electromagnetic fields on biological tissue [21]. All previous works on the memristive HR neuron model have considered only the 3D version of the HR model [20,22,23,24,25,26,27], leading to a 4D model including the memristive (magnetic flux) variable. Thus, for the first time, the memristive property of the 4D version of the HR model would be considered in this paper. Including the memristive dynamics into the 4D HR model leads to a 5D HR model on which the study in this paper will focused on. The aim of this paper is not to present a detailed bifurcation analysis of this improved 5D HR neuron model, but just to show the richer dynamical behavior the model can reproduce. On the other hand, because synchronization of memristive nonlinear systems has become an active area of research, and has widespread applications in secure communication and neuromorphic circuits [28,29], the paper mainly focuses on the synchronization dynamics of the memristive 5D HR model.
Using the memristive 5D HR neuron model, an alternative way to determine the asymptotic stability of the synchronization manifold of coupled chaotic systems is presented. We show that the time variation of the Hamilton function of the error dynamical system associated to a pair of coupled chaotic 5D HR neurons can be used to determine the asymptotic stability of the synchronization manifold, just as the time rate of change of Lyapunov function along trajectories of the system would do. It is shown that as the synaptic coupling strengths and memristive gain parameters of the model are varied, we always have a non-zero time variation of this Hamilton function only when the time variation of the Lyapunov function is positive (indicating an unstable synchronized state), and zero (or vanishingly small) only when the time variation of the Lyapunov function is also zero (indicating an asymptotically stable synchronized state). This indicates that the time variation of the Hamilton function associated to the error dynamical system of coupled oscillators can be used as an asymptotic stability function.
The rest of this paper is organized as follows: In Sect.2, we present the improved 5D HR neuron model and show its rich dynamical behavior including spiking, bursting, super-bursting an chaos in time series, phase portraits and bifurcation diagrams. In Sect.3, we investigate the chaotic synchronization dynamics of a pair of neurons in the chaotic super-bursting regime and coupled via both time-delayed electrical and chemical synapses. The first part of this section is devoted to the analysis of the asymptotic stability of the synchronization manifold of the coupled neurons via the Lyapunov stability theory. The second part is devoted to the energy of the synchronization dynamics via the Hamilton function. Sect.4 deals with numerical simulations, showing the sign correlation between time variations of the Lyaponuv and the Hamilton function of the system. In Sect.5, we summary and concluding remarks.

Mathematical model and dynamics behavior
The dynamical equation for the memristive 5D HR neuron model is described by where x ∈ R represents the membrane potential variable, y ∈ R is the recovery current variable associated to fast ions, z ∈ R is a slow adaption current associated to slow ions, w ∈ R represents an even slower process. I 0 is the amplitude of a harmonic stimulus with frequency Ω and phase ψ. The other constant parameters have standard values: a = 1.0, b = 3.0, p = 0.99, c = 1.01, d = 5.0128, σ = 0.0278, r = 0.00215, s = 3.966, x 0 = 1.605, µ = 0.0009, γ = 3.0, y 0 = 1.619, δ = 0.9573. The parameters µ < r 1 play a very important role in neuron activity; r represents the ratio of timescales between fast and slow fluxes across the neurons membrane and µ controls the speed of variation of the slower dynamical process w, in particular, the calcium exchange between intracellular warehouse and the cytoplasm [13,15]. The µ timescale induces richer dynamical behaviors including chaos in parameter regimes that the 3D HR shows only periodic dynamics [17].
The fifth variable φ ∈ R, describes the magnetic flux across the neuron's cell membrane. The term ρ(φ) = α + 3βφ 2 is the memory conductance of a magneticflux controlled memristor [30,31,32,33], where α and β are fixed parameters which we will fixed throughout the paper at α = 0.1, β = 0.02. k 1 and k 2 are feedback gain; k 1 bridges the coupling and modulation on membrane potential x from magnetic field φ, and k 2 describes the degree of polarization and magnetization by adjusting the saturation of magnetic flux [34]. Following the Faraday's law of electromagnetic induction and the basic properties of a memristor, the term k 1 ρ(φ)v could be regarded as additive induction current on the membrane potential. The dependence of electric charge q on magnetic flux φ is defined by the memory-conductance as follows [25,35] Moreover, we know that current i is defined by the time rate of charge q. Hence, the physical significance for the term ρ(φ)v could be described as follows where the variable V denotes the induced electromotive force, which holds a same physical unit, and parameter k 1 is the feedback gain. The ion currents of sodium and potassium contribute to the membrane potential and also the magnetic flux across the membrane; thus, a negative feedback term −k 2 φ is introduced in the fifth equation of Eq. (1). We select the amplitude of harmonic forcing current I 0 , its phase ψ, and the memristive gain parameters k 1 and k 2 as: I 0 = 1.6, ψ = 0.1, k 1 = 1.0, and k 2 = 0.5. In Fig.1, the time series for membrane potential x under different values of the frequency Ω of the harmonic forcing are shown. The electrical activity of the 5D HR neuron show a rich dynamical behavior. In Fig.1(a), the model displays a simple periodic spiking activity with Ω = 0.2. In Fig.1(b), Ω decreases to 0.02 and the simple periodic spiking activity changes to the standard bursting activity, with each burst consisting of eleven or twelve spikes. In Fig.1(c), when the frequency is further reduced to Ω = 0.003, the standard bursting activity changes to a super-bursting activity, where each super-burst consist of three standard bursts which in turn consist of different number of spikes. In Fig.1(d), the frequency is increased by a little, i.e., to Ω = 0.0036, the pattern of the super-burst changes, with each super-burst consisting this time of only two standard bursts which have different spiking patterns.
The timing precision of the information processing in neural systems is very important because the information in encoded at the spiking and bursting times [36]. This means that in our 5D neuron model in a super-bursting regime, one kind of information could be encoded at super-burst time; then a second kind of information encoded at the standard burst times of each of these super-bursts; and then a third kind of information encoded at the spiking time of each standard burst of a super-burst. Whereas in a simple spiking regime, only one kind of information could be encoded at a time. It is well known that bursts provide a more reliable mode of information transfer than spikes [37].
In order to investigate the dependence of the system's behavior on its parameters -the driving harmonic current (I 0 , Ω) and the memristive gain parameters (k 1 , k 2 ) -several bifurcation diagrams, each fully traced by its corresponding maximum Lyaponuv exponent, have been computed by using the peaks of the spikes in the super-bursting regime of Fig.1(c). The diagrams shown have been chosen to illustrate the general dynamical structure of the memristive 5D HR neuron model. The maximum Lyaponuv exponent, Λ max , of the model which is defined as where L(τ ) = δx 2 + δy 2 + δw 2 + δz 2 + δφ 2 1/2 is computed numerically by solving, simultaneously, the system in Eq. (1) and its corresponding variational system of equations. Starting from an initial condition, the system of Eq. (1) is numerically integrated with the fourth-order Runge-Kutta algorithm and the data recorded after some transient time. A strictly negative maximum Lyapunov exponent characterizes an asymptotically stability system and the more negative the exponent the greater the stability. When Λ max = 0, we have a marginally stable system with quasi-periodic trajectories. And when Λ max > 0, the trajectories are un-stable, meaning two nearby trajectories would diverge and the evolution of the system becomes sensitive to infinitesimal perturbations of initial conditions and hence chaotic [38].
show the bounded phase space of a chaotic attractor of the 5D memristive HR neuron in a the super-bursting regime, projected onto the 3D subspaces of the 5D phase space: xyφ-space, xzφ-space, and xwφ-space, respectively. In Fig.3(a) and (b), we respectively show a bifurcation diagram and its corresponding Lyapunov spectrum Λ max . Different bifurcation sequences occur as amplitude of the external harmonic current is varied. The system dynamics is mainly chaotic and it is intermingled with thin windows of periodic orbits. In Fig.3(c) and (d), different bifurcation sequences occur as of the frequency of the external harmonic current is varied. The system dynamics is mainly chaotic for lower frequency values.
In Fig.4(a) and (b), we respectively show a bifurcation diagram and its corresponding Lyapunov spectrum Λ max for the memristive parameters k 1 . The system dynamics is mainly chaotic for intermediate values of k 1 . In Fig.4(c) and (d), chaotic and periodic dynamics are intermingled as k 2 is varied. In the next section of the paper, the values of the memristive parameters shall be fixed in a chaotic regime.

Synchronization dynamics
This section deals with the main topic of this paper -the chaotic synchronization dynamics of the coupled 5D HR neuron. Here, we consider a pair HR neurons coupled via both time-delayed electrical and inhibitory chemical synapses. Electrical and chemical synapses are the two ways through which neurons connect to each other [39]. Electrical synapses connects the cytoplasm of nearby neurons directly and as a result the transmission of electrical impulses occur relatively quickly. The corresponding functional form of the bidirectional in-teraction mediated by the electrical synapses is defined as the difference between the membrane potentials of two adjacent neurons. On the other hand, with chemical synapse, the transmission of information takes place via the release of a neurotransmitter. The functional form of this synaptic interaction is considered as a nonlinear sigmoidal input-output function [40]. The system of coupled memristive 5D HR neurons is given as with i, j = 1, 2. The parameters g e and g c are the electrical and chemical coupling strengths, respectively. The chemical synaptic function is modeled by a sig- where the parameter λ = 10.0 determines the slope of the function and θ s = −0.25 denotes the synaptic firing threshold. V syn represents the synaptic reversal potential. For V syn < x i the chemical synaptic interaction has a depolarizing effect that makes the synapse inhibitory, and for V syn > x i , the synaptic interaction has a hyperpolarizing effect making the synapse excitatory. For the 5D HR neuron model, the membrane potentials are bounded as |x i (t)| ≤ 2.0 (i = 1, 2) for all times t. For the choice of fixed V syn = −2.5 (maintained throughout computations), the term (x i − V syn ) in Eq. (5) is always positive. So, the inhibitory and excitatory natures of chemical synapses will depend only on the sign in front of the synaptic coupling strengths g c . To make the chemical synapse inhibitory, we chose a negative sign. The connectivity matrices ξ i and η ij are such that:

Asymptotic stability of synchronous states and Lyapunov function
The complete synchronization of the coupled system in Eq. (5) occurs when the two neurons asymptotically exhibit identical behavior, that is as t → ∞, for initial conditions chosen in some neighborhood of the synchronization manifold M s on which The synchronization solution in Eq. (7) which satisfies the system of dynamical equations where x j − x i = 0, is also always a solution of coupled system in Eq. (5). However, this synchronization solution might be stable only under some conditions. By introducing coordinates transformation defined by directions transverse to the synchronization manifold M s as we obtain the dynamics of the transverse perturbations to the synchronization manifold as The following theorem can be obtained.
The time derivative of the Lyapunov function V along trajectories of the error dynamical system in Eq. (10) yields Since chaotic systems have bounded phase space, there exists a positive constant J, such that |x(t)| < J, thus for every points of the attractor and can be compactly written as where E T = e x e y e z e w e φ is a row vector with transpose E, and the matrix M is given by with D 1 = (3aJ − 2b)J + k 1 α + g c 1 + e −λ(J−θs) −2 . To ensure that the origin of the error dynamical system in Eq. (10) is asymptotically stable, three conditions have to be met in Eq. (13): a ≥ 0 (which is always the case since a = 1.0 in the model), g e ≥ 0, and most importantly, the matrix M must be positive definite. M is positive definite if it satisfies Sylvesters criterion given in Appendix A.2 -in which case dV /dt in Eq. (13) will be negative semi-definite. According to the Lyapunov stability theory [41] and Barbalats lemma [42], all the transverse perturbations decay to the synchronization manifold without any transient growth, i.e., one obtains as t → ∞. It follows that the coupled system in Eq. (5) synchronizes when the inequalities in Eq. (A.2) are satisfied. This completes the proof.

Synchronization energy and Hamilton function
Here, we determine the Hamilton energy function H(e x , e y , e z , e w , e φ ) associated to the error system in Eq. (10). Using this energy function, we analytically evaluate the energy variation of the coupled system. This energy variation is an important measure because it gives us the flow of energy in the process of synchronization and hence, it can be considered to be the amount of energy per unit time needed to maintain a particular degree of synchrony [43]. We focus on the effects of the synaptic couplings strengths (g e , g c )and magnetic flux parameters (k 1 , k 2 ) on this energy variation and compare it to the time variation of the Lyapunov function previously calculated.
Based on Helmholtz's theorem [44], we express the velocity vector field F (e x , e y , e z , e w , e φ ) of the error dynamical system in Eq. (10) as the sum of two vector fields: a divergence-free vector and a gradient vector field, i,e., F (e x , e y , e z , e w , e φ ) = f c (e x , e y , e z , e w , e φ ) + f d (e x , e y , e z , e w , e φ ), with the conservative part f c (e x , e y , e z , e w , e φ ), containing the full rotation and the dissipative part, f d (e x , e y , e z , e w , e φ ), containing the whole divergence. A divergence-free vector and a gradient vector field of the error dynamical system in Eq. (10) are given by where For the conservative vector field, the equation where ∇H denotes the transpose gradient of function H, defines a partial differential equation from which a Hamilton energy function (a generalized Hamiltonian) H(e x , e y , e z , e w , e φ ) can be evaluated (see also [45] for a detailed and general description of the method). Thus, Hamilton energy function of the error dynamical system in Eq. (10) satisfies the partial differential equation given by By the method of separation of variables, one solution of Eq. (19) is given by where Q is an arbitrary constant. Since any positive definite quadratic form can always be a solution for the energy partial differential equation compatible with the generalized Hamiltonian formalism, independently of the system at hand, the same trivial positive definite quadratic form can always be assigned to different chaotic systems. However, assigning the same form of energy function to all chaotic systems fails to reveal the individual features of its own dynamics. Hence, the approach used to obtain the Hamilton energy function H requires additional hypothesis in order to be able to assign to the error dynamical system a specific energy function. This additional hypothesis is set by introducing a connection between the the change in the volume of the phase space of the coupled system and time rate of change of energy, as one cannot occur without the other. Hence, since the energy function in Eq. (20) is not unique for the system, we compute the time rate of change of this energy function along trajectories of the error dynamical system. This time rate of change of the Hamilton energy function is related to the divergence of the vector field, f d (e x , e y , e z , e w , e φ ), responsible for contraction of the volume of the phase space. Thus, the time variation of the Hamilton energy function H(e x , e y , e z , e w , e φ ) associated to the error dynamical system in Eq. (10) now becomes uniquely related to the specific dynamics of the error dynamical system. This energy is dissipated via the dissipative component of the velocity vector field f d (e x , e y , e z , e w , e φ ) according to the equation [45] dH dt = ∇H T f d .
Irrespective of the value of the arbitrary constant Q (which we will fixed at Q = −1.0 throughout the rest of the paper), Eq. (22) gives us an energy function which is unique to the error dynamical system of our coupled neurons, and it gives us information about the energy dissipated during the synchronization dynamics. In the next section, we show (and provide a theoretical explanation) that the time variation of the Hamilton function of the error dynamical associated to the coupled neuron system paves an alternative way of determining the asymptotic stability of the synchronized state of the system just as the time variation of the Lyapunov function would do.

Numerical simulations and discussion
In this section, we compare the derivatives of the Lyaponuv function and the Hamilton energy function given in Eq. (12) and Eq. (22), respectively, as the synaptic strengths (g e and g c ) and the memristive gain parameters (k 1 and k 2 ) vary. To simulate these equations, we simultaneously integrate Eq. (8) and Eq. (10) for a very large time interval using the fourth-order Runge Kutta algorithm. After discarding the transient time, we compute the mean values of dV /dt and dH/dt. For a weak chemical synaptic strength g c , Fig.5(a) shows the variations of dV /dt and dH/dt with respect to electrical synaptic strength g e . It is observed that synchronization manifold M s is always unstable as indicated by dV /dt > 0 (blue curve). While for these same values of g e , we always have dH/dt < 0 (red curve).
In Fig.5(b), we show the variations of dV /dt and dH/dt with respect to the chemical synaptic strength g c for a weak electrical synaptic strength g e . Here, we see that dH/dt < 0 only when M s is unstable as indicated by dV /dt > 0; and dH/dt = 0 only when M s is stable as indicated by dV /dt = 0. In Fig.6(a) and (b), we show a color-coded global behavior of the dV /dt and dH/dt with respect to the synaptic coupling strength parameters g e and g c for fixed memristive gain parameters. Here we oberve that the sign correlation between dV /dt and dH/dt persist for all values of the parameters, and not just for particular values. Here, we have 0 < dH/dt < 0 only when dV /dt > 0, indicating an unstable synchronization manifold. And dH/dt = 0 (or vanishingly small, dH/dt ≈ 0) only when dV /dt = 0, indicating a stable synchronization manifold.
In Fig.6(a), the color-bar shows the variation of dV /dt for the given range of values of g e and g c and indicates that for a sufficiently strong chemical coupling strength, i.e., g c > 1.25, the synchronization manifold becomes and stays stable, (i.e., dV /dt = 0) irrespective of the electrical synaptic strength g e .
In Fig.6(b), where the color-bar shows the variation of dH/dt, we observed that dH/dt = 0, for g c > 1.25 irrespective of the value of g e just as with dV /dt in Fig.6(a). For weak chemical coupling, i.e., g c < 1.25, we have 0 < dH/dt < 0 (in Fig.6(b)) for the same parameter values of g e and g c for which dV /dt > 0 (in Fig.6(a)), indicating unstable synchronized dynamics. However, in this weak chemical coupling regime (g c < 1.25), we can have some values of g c (e.g., g c = 0.25 and g c = 1.0) for which we have stable synchronized dynamics, i.e., dV /dt = 0 with dH/dt = 0 or dH/dt ≈ 0. Fig. 6 Color-coded variations of dV /dt in (a)) and dH/dt in (b)) with respect to the synaptic coupling strengths for fixed memristive gain parameter values. The synchronization manifold M s is unstable when dV /dt > 0 which occurs only when 0 < dH/dt < 0; and stable when dV /dt = 0 which occurs only when dH/dt = 0 or dH/dt ≈ 0 (i.e., vanishingly small). Other parameters: k 1 = 1.0, k 2 = 0.5.
In Fig.7, we show the variations of dV /dt and dH/dt with respect to the memristive gain parameters k 1 and k 2 in a weak and strong synaptic coupling regimes. In Fig.7(a), we consider weak synaptic coupling strength and we always have dH/dt < 0 whenever dV /dt > 0 for all the values of k 1 , indicating unstable synchronized states. In Fig.7(b), we switch to a stronger synaptic coupling strengths and we get stable synchronized states as dV /dt = 0 and dH/dt = 0 for intermediate values of k 1 , and unstable synchronized states (dV /dt > 0 and dH/dt = 0) for small and large values of k 1 .
In Fig.7(c) and (d), we also consider k 2 in weak and strong synaptic coupling regimes, respectively. And again we always the same behavior, i.e., synchronized states for dV /dt = 0 and dH/dt = 0, and unstable synchronized state when dV /dt > 0 and dH/dt = 0. Especially, in Fig.7(d), where dV /dt and dH/dt vary a lot, we can see that the signs of dV /dt and dH/dt are always mirror images of each other about zero-symmetry line, i.e., dH/dt = 0 only when dV /dt = 0 (stable synchronized states), and dH/dt < 0 only when dV /dt > 0 (unstable synchronized states). To have a global view on the behavior of dV /dt and dH/dt with respect to the memristive parameters k 1 and k 2 , we computed these functions in a two-parameter space. Fig.8(a) and (b) show a color-coded dV /dt and dH/dt as a function of k 1 and k 2 , respectively, in a strong synaptic coupling regime. In Fig.8(a), we observed synchronized states (the black regions) where dV /dt = 0 corresponds to the black regions in Fig.8(b) where dH/dt = 0 or or dH/dt ≈ 0. And wherever dV /dt > 0 in Fig.8(a), we have dH/dt = 0 (i.e., either positive or negative), indicating unsynchronized states.
In Fig.8(c) and (d), we switch to a weak synaptic coupling regime. Here, panels (c) and (d) also show a color-coded dV /dt and dH/dt as a function of k 1 and k 2 , respectively. Here, we mostly have dV /dt > 0 (unsynchronized state) in regions where dH/dt = 0. However, many black spots (where dV /dt = dH/dt = 0, indicating synchronized states) can also be observed. For example, in panels (c) and (d), around the coordinates (k 1 , k 2 ) = (1.5, 0.5), a black patch of synchronized states (dV /dt = dH/dt = 0) can be observed. Fig. 8 Color-coded variations of dV /dt in panel (a) and the corresponding dH/dt in panels (b), with resp. to k 1 and k 2 in a strong coupling regime with g e = 4.0, g C = 2.75. We see that stable synchronized states emerge at values of k 1 and k 2 for which dV /dt > 0 and 0 < dH/dt < 0, while the unstable synchronized state at dV /dt = 0 and dH/dt = 0. In a weak coupling regime with g e = 1.5, g c = 0.75 panels (c) and (d) respectively show the color-coded variations of dV /dt and dH/dt. Here, we also have 0 < dH/dt < 0 only when dV /dt > 0 and dH/dt = 0 (or dH/dt ≈ 0) only when dV /dt = 0.
A theoretical explanation for this sign correlation between time rate of change of the Lyaponuv dV /dt and Hamilton dH/dt functions and hence, the ability of dH/dt to also indicate whether or not a synchronized state is asymptotically stable is the following: First, we notice that dH/dt can be either positive or negative (only when dV /dt > 0, indicating unstable synchronized states). When dV /dt > 0 , the chaotic system (i.e., the coupled neurons) initially located outside its synchronization manifold would gain (dH/dt > 0) or lose energy (dH/dt < 0) in its movement towards synchronization manifold, where dV /dt = 0 (indicating an asymptotically stable synchronous state) and dH/dt = 0 or dH/dt ≈ 0 (i.e., or vanishingly small). This is so because on the synchronization manifold , the trajectory will repeatedly return to arbitrarily close states in the bounded phase space of the chaotic attractor and as a result to arbitrarily close energy values. Hence, on synchronization manifold the average time rate of of Hamilton energy function will be zero, i.e., dH/dt = 0 or vanishingly small, i.e., dH/dt ≈ 0. This implies that, in general, all the different regimes of synchronization that the two neurons attain at different values of the system's parameter would occur at zero (or very low) net dissipation of energy (i.e., at dH/dt = 0 or dH/dt ≈ 0). But, there could be ranges of parameter values where the activity of the coupled system is more demanding energetically, that is when 0 < dH/dt < 0 which as we result from the numerical simulation occur only when dV /dt > 0, indicating unstable synchronized states. This means that when the neurons are out of the synchronization manifold, a non-zero energy dissipation (0 < dH/dt < 0) is thus necessary to drive the coupled neurons to the synchronization manifold where the net energy dissipated become zero or very low (dH/dt = 0 or dH/dt ≈ 0) and where the time rate of change of the Lyaponuv function is also zero (dV /dt = 0) -indicating an asymptotically stable synchronized state following Krasovskii-Lyapunov stability theory. Hence, the time rate of change of the Hamilton function of an error dynamical system associated to a coupled system can be used as a synchronization stability function.

Summary and concluding remarks
In this work, we have investigated the dynamics of an improved version of the standard 3D Hindmarsh-Rose neuron model by taking into account not only the exchange of calcium ions across the cell membrane, but by also considering the dynamics of the magnetic flux induced across the membrane as a result of the motion of ions. The electric activity of the improved 5D Hindmarsh-Rose model shows upon variations of the external input current and the magnetic flux (memristive gain) parameters, a rich dynamical behavior including periodic/chaotic spiking and bursting, chaotic super-bursting -dynamical behaviors observed in real biological neurons upon variations of the corresponding parameters.
To investigate the synchronization dynamics of two coupled 5D Hindmarsh-Rose neurons, we considered a pair neurons coupled via both a instantaneous electrical and inhibitory chemical synapses. Using the Krasovskii-Lyapunov stability theory, we prove that the synchronization manifold of the coupled neurons can be asymptotically stable for suitable values of the electrical and chemical coupling strengths g e and g c , and the memristive gain parameters k 1 and k 2 . Moreover, we used Helmholtzs theorem to calculate the Hamilton energy function associated to the error dynamical of system of the coupled neurons. Numerical computations indicated that we always have a non-zero (0 < dH/dt < 0) time variation of the Hamilton function only when the time variation of the Lyapunov function is positive (dV /dt > 0), and zero (dH/dt = 0) (or vanishingly small, dH/dt ≈ 0) only when the time variation of the Lyapunov function is also zero (dV /dt = 0). Thus, the time variation of the Hamilton energy function of the error dynamical system associated to a coupled system, can be used as an asymptotic stability function for synchronization. This result which might be also useful for general engineering purposes, paves an alternative way of determining the asymptotic stability of the synchronized states in coupled systems from an energy perspective without necessarily having to construct a Lyapunov function which might be a difficult task for systems modeled with more complicated mathematical equations.
The work presented in this paper could be extended in two directions. First, by considering synaptic and/or channel noise (e.g., Gaussian white noise or non-Gaussian colored noises), due to their presence and relevance in real neural dynamics. And secondly, by investigating the variation of the synchronization energy between the layers of a multiplex neural network -a relevant network structure ubiquitous in the brain.
The five ordered-main sub-determinants of matrix M in Eq. (14) are given by It is easy to see that, for the given values of a ,b, d, r ,s ,p, µ, α, δ, γ, σ, and for suitable values of g c , k 1 , and k 2 , we can always have