Energy feedback and synchronous dynamics of Hindmarsh-Rose neuron model with memristor

K Usha and P A Subha

1 Department of physics,University of Calicut,Kerala 673635,India

2 Department of physics,Farook College,University of Calicut,Kerala 673632,India

Keywords:HR model,memristor,Hamilton energy,energy feedback,synchronization

1.Introduction

The transmission of nerve impulses in the brain occurs via the propagation of action potentials and a large fraction of the total energy consumed by brain is utilized to generate the firing patterns.The evaluation of required metabolic energy to maintain the signaling activity in neurons is an active research area.[1]The collective dynamics of neural networks have a considerable in fluence in the propagation of information from one region to another and their collective behavior can be greatly in fluenced by the energy demands.Neuron models are used to analyze the connection between energy demands and firing modes.[2]The Hamilton energy associated with the dynamical system can be derived using the generalized Hamiltonian approach.[3]Using this formalism,Sarasola et al.have derived the energy function associated with Lorenz,Rossler,and Chua systems.[4]Moujahid et al.have reported the energy consumption during the synchronization process in electrically coupled HR neurons.[5]The consumption of energy in transmitting nerve impulse using Hodgkin-Huxley model has also been reported.[6]

It is important to investigate the energy utilization of neurons subjected to different kinds of external force during information encoding.Biological experiments confirm that the electrical activities of neurons change by adjusting extracellular calcium or potassium concentrations.[7]The Hindmarsh-Rose(HR)neuron oscillates periodically for small and large external current,whereas for intermediate currents they become chaotic.[8]The periodic external forcing also in fluences the bursting modes and neuronal activities.[9]In addition,computational models have been developed to study the effect of noise.[10-12]The statistical features of Gaussian white noise seem to be appropriate to mimic the complex behavior shown by a neuron under the in fluence of other neurons and the environment.Zambrano et al.have analyzed the synchronization of uncoupled FitzHugh-Nagumo neurons with common noise both experimentally and numerically.[13]Noise-induced resonances in the HR model has also been studied.[14]Wang et al.have reported that,the forced HR neurons with white Gaussian noise leads to firing modes like multi-modal firing,intrinsic oscillation,and bi-modal firing.[15]Lindner et al.have studied the dynamics of mathematical models of excitable systems in the presence of white Gaussian noise.[16]

Recent studies confirm that a two terminal electric device called memristor can mimic the the key characteristics of synapses and neurons.A memristor models the effect of electromagnetic field created as a result of the exchange of ions across the nerve membrane.The electrical activities in the cardiac tissues exposed to electromagnetic radiation can be described using memristor models.This provides signi ficant clues about the mechanism of heart disorders induced by electromagnetic radiation.[17]The synchronization of coupled memristive neural network via pinning control has been proposed by Guan et al.[18]Recently Ma et al.have discussed the RCL-shunted junction circuit with memristor.[19]The control of dynamical systems with negative feedback in energy has recently been reported.[20]

In this paper,we analyze the energy aspects of single and coupled HR neuron model with a memristor.The HR neuron model with quadratic flux controlled memristor is presented in Section 2.The Hamilton energy of the system is derived in Section 3 and the energy change in the presence of different external stimuli has been discussed.In addition,the control of chaotic trajectories by applying a negative feedback is discussed in this section.The energy aspects during the synchronization process of electrically and chemically coupled HR neurons are studied in Section 4.Finally,Section 5 concludes the study.

2.Model

The continuous exchange ofcharged ions across the nerve membrane induce an electromagnetic field that controls the membrane potential of neurons.Recently Ma et al.proposed that memristors can be used to bridge the magnetic flux and membrane potential.This coupling in fluences signal transmission via the superposition of electric field.[9,21]Memristors are used to realize the coupling.The HR model with memristor is capable of producing biologically relevant dynamical states such as anti-phase oscillations,co-existence of resting and spiking state etc.[22]The dynamical equations of HR neuron model with quadratic flux controlled memristor have the form:[9]

where ρ(φ)= αφ2+βφ +γ.The membrane potential of the neuron is represented using the variable x.y denotes the recovery variable representing the rate of change of fast current of K+or Na+ions and z denotes the adaptation variable that capture the slower dynamics of other ion channels.The parameters a and b denote activation and inactivation of the fast ion channel.R and xedescribe activation and inactivation of the slow ion channel.The speed of variation of z is controlled by r.[23-25]The parameter I represents the external stimuli.[26,27]The fourth variable φ denotes magnetic flux across the nerve membrane.A memristor with memductance ρ(φ)=d q(φ)/dφ is used to realize the coupling between magnetic flux of the field and membrane potential.It is possible to model ρ(φ)using a quadratic term.The memductance after suitable scaling is taken as ρ(φ)= αφ2+βφ +γ,α,β,and γ are parameters.[28]The term k1ρ(φ)x denotes the induced current through electromagnetic induction,where k1represents the modulation intensity of electromagnetic field.Relation between induced current and flux change can be understood using Faraday’s law of electromagnetic induction.[29]The term k2x represents the change in magnetic flux induced by membrane potential of the cell and k3φ denotes the leakage of magnetic flux.The parameters used are a=3.0,b=5.0,R=4.0,r=0.006,I=3.1,and xe=-1.61.[30]

We have analyzed the dynamics of the system in Eq.(1)by varying the external current I.The fourth order Runge-Kutta algorithm is applied for numerical calculations.The memristor parameters are taken as β =1.0,γ=1.0,k2=0.9,and k3=0.5.The value of k1is fixed at 0.1 and by increasing the value of α the inter spike interval bifurcation(IS I)diagrams have been plotted.Figures 1(a),1(b),and 1(c)show the distribution of ISI for α=0.1,0.4,and 0.8 respectively.It is found that as the value of α is increased,the system eventually transforms to its normal response state.The study has been extended by fixing α=0.1 and by varying k1.Figures 1(d),1(e),and 1(f)represent the I versus IS I for k1=0.8,0.4,and 0.1 respectively.From the plots it is clear that,for constant α,a decrease in k1is needed for normal firing.The response of a neuron with electromagnetic induction described by quadratic flux controlled memristor to external signals can be enhanced by properly selecting the memristor parameters.[21]

Fig.1.The ISI bifurcation diagram of HR neuron with memristor.Figures 1(a),1(b),and 1(c)shows I versus ISI for k1=0.1 and α=0.1,0.4,and 0.8.Figures 1(d),1(e),and 1(f)are drawn for α =0.1,k1=0.8,0.4,and 0.1.The parameters β and γ are set as 1.0.

3.Energy aspects

In this section,we derive the Hamilton energy of HR model with quadratic flux controlled memristor.The differential equation of an autonomous dynamical system is of the form:˙x=f(x).According to Helmholtz’s theorem the velocity vector field f(x)can be written as:

where fccomponent of the vector field is conservative.This does not contribute to the energy change along any trajectory of the system and satis fies the following equation,

The function H(x)is the generalized Hamiltonian for the conservative system as long as it can be rewritten in the form˙x=J(x)?H,where J is a skew symmetric matrix that satis fies Jacobi’s closure condition.[31]The component fdis composed of velocity-dependent terms and contribute to the divergence.[9]The dissipation of energy due to the fdpart obeys the relation:

The conservative and dissipative part of HR model in Eq.(1)can be expressed in the form:

where

and

Then,according to Eq.(3),the Hamilton energy associated with the system will satisfy the following partial differential equation:

A general solution for Eq.(5)is of the form:

The rate of change of Hamilton energy function is:

simplifying and rearranging Eq.(9),the expression for˙H takes the form,and hence obeys the relation in Eq.(4).

The average energy of the system is evaluated using the expression,

where T is the energy calculation period(1000 time units)and t0is the starting time to calculate the average energy.Figure 2 shows the variations in the average energy with the external forcing current(I).The parameters used are taken as:a=3.0,b=5.0,R=4.0,r=0.006,xe=-1.61,k1=0.1,k2=0.9,and k3=0.5,α =0.4,β =0.02,and γ=0.1.[9]From the plot it is clear that as I is increased,the average energy decreases.

Fig.2.The variations in average energy of HR neuron with quadratic flux controlled memristor for different values of external current.

We have further analyzed the rate of change of Hamilton energy in the resting and bursting state of membrane potential by applying various external stimuli.

Case 1:Constant external stimulus

A constant external current(I)has been applied to the system.The value of I is changed from 2.0 to 3.0 at t=1000 and then switched to 1.0 at t=1500,as depicted in Fig.3(a).

Fig.3.(a)Variation of external current(b)rate of change of energy function and(c)time evolution of membrane potential of neuron.The external current I is changed from 2.0 to 3.0 at t=1000 time units and switched to 1.0 at t=1500 time units.

The corresponding variations in energy utilization and membrane potential are shown in Figs.3(b)and 3(c)respectively.From the plots it is clear that as the external current is varied the bursting mode of neuron change and for generating each action potential energy is consumed.The energy demand is a maximum during the repolarization period of the spike and at a minimum during the refractory period between two spikes.The energy utilization approaches zero when the membrane potential is close to the quiescent state.

Case 2:Periodic external stimulus

The effect of periodic stimulus in neural activity has been analyzed by applying an inputofthe form I=I1+A sin(0.05t),where I1=3.1 and A represents the amplitude.The membrane potential and energy utilization during the electrical activities are plotted.The value of A is changed from 0.1 to 1.0 at t=1000 and then switched to 2.0 at t=1500,as shown in Fig.4(a).In the case of periodic input,A acts as a control parameter for generating different types of electrical responses.As A is increased,the number of spikes per burst in membrane potential is also increased;as shown in Fig.4(c).The transition in bursting mode induce some transition in Hamilton energy.In the case of state with less number of spikes per burst,the Hamilton energy consumed also become smaller as depicted in Fig.4(b).

Fig.4.(a)Variation of periodic input(b)rate of change of energy function and(c)time evolution of membrane potential of neuron.The external periodic forcing has the form I=I1+A sin(0.05t).I1 is fixed at 3.1 and A is changed from 0.1 to 1.0 at t=1000 and then switched to 2.0 at t=1500.

Case 3:White Gaussian noise

The in fluence of noise in the firing pattern is discussed in Fig.5.

Fig.5.External noise variations(b)rate of change of energy function and(c)time evolution of membrane potential of neuron.The noise applied is I=I1+ζξ(t).I1=3.1,ζ is changed from 0.01 to 0.1 at t=1000 and then switched to 1.0 at t=1500.

White Gaussian noise is added through the electrical potential of the membrane;i.e.,the effective current imposed to the neuron contains a random term.[32]The noise is of the form I=I1+ζξ(t),with parameters 〈ξ(t)〉=0,〈ξ(t)ξT(t+τ)〉=δ(τ)and ζ defines the noise intensity.The value of I1is fixed as 3.1 and ζ is changed from 0.01 to 0.1 at t=1000 and then switched to 1.0 at t=1500 as shown in Fig.5(a).The corresponding variations in membrane potential are shown in Fig.5(c).From the plotitisclearthatthe bursting mode changes with the increase in noise intensity.The energy utilization in the presence of noise is depicted in Fig.5(b).It is confirmed that noisy external stimulus can trigger complex discharge in energy utilization and the effect of noise is more evident in the refractory period of the action potential.

The bifurcation diagrams in(A-Xmax)and(A-ISI)planes for periodic input are shown in Figs.6(a)and 6(b)respectively.The amplitude of periodic forcing is varied as 0≤A≤4.The plots show that,an alternate sequence bursting states occur with the increase in A.The bifurcations in the presence of the noisy external forcing are shown in Figs.6(c)and 6(d).The noise intensity is varied in the range 0≤ζ≤1.It is found that an increase in ζ produces the complex rhythm and the electrical discharge of the nerve cell becomes more complex leading to chaos.

Fig.6.The bifurcation diagrams.Figures 6(a)and 6(b)show X max and IS I by varying the amplitude of periodic forcing.Figures 6(c)and 6(d)show X max and ISI by varying the intensity of external noise.

4.Energy feedback

The Hamilton energy of the system depends on all system parameters,and hence the changes made to energy function causes significant change in the phase space of the dynamic system.[20]The change in energy is realized by giving a negative feedback as follows:

where k4is the feedback gain.The phase space dynamics in the(X-H)plane by varying k4has been plotted.

Case 1:Constant external stimulus

The formation of attractors in the presence of constant current(I=3.1)is illustrated in Fig.7.The phase space for k4=0.0,k4=1.0,k4=5.0,k4=10.0 are shown in Figs.7(a),7(b),7(c),and 7(d),respectively.As the feedback gain is increased,the number of dense orbits in the attractor is reduced and the chaotic trajectories are controlled.The results are further confirmed by plotting the largest Lyapunov exponent(LLE).The variation of LLE with k4is shown in Fig.8.It is found that the LLE is decreased below zero with increase in k4and ensure the stabilization of chaotic trajectories.

Fig.7.The phase space of dynamics for I=3.1.The energy feedback obeys Eq.(12).(a)k4=0,(b)k4=1.0,(c)k4=5.0,and(d)k4=10.0.

Fig.8.Transition of LLE for different feedback gains in energy function in the presence of constant external current.The inserted figure is the enlarged version.

Case 2:Periodic external stimulus

The phase space dynamics in the presence periodic input I=I1+2sin(0.05t)for k4=0,k4=1.0,k4=5.0,and k4=10.0 are shown in Figs.9(a),9(b),9(c),and 9(d)respectively.

Fig.9. The phase space of dynamics for periodic external forcing I=I1+A sin(0.05t),where I1=3.1 and A=2.0.The energy feedback is according to Eq.(12).(a)k4=0,(b)k4=1.0,(c)k4=5.0,and(d)k4=10.0.

The plots show that,as the feedback gain k4is increased the Hamilton energy function which is composed of the variables and bifurcation parameters controls the evolution of the system and the chaotic trajectories are stabilized.The variation of LLE is shown in Fig.10.The maximum value of LLE is nearly equal to 0.008,greater than compared to the one obtained for constant external current.

Fig.10.Transition of LLE with k4.Periodic variations are applied in the external current.The inserted figure is the enlarged version.

Case 3.White Gaussian noise

The phase portrait of the system with energy feedback in the presence of the external noise I=I1+ζξ(t)has also been studied with I1=3.1 and ζ=1.As the feedback gain k4is increased,the noisy trajectories are controlled.Figures 11(a),11(b),11(c),and 11(d)represent the mechanism of the chaos control for k4=0,k4=1.0,k4=5.0,and k4=10.0 respectively.The LLE for noisy external forcing has the largest value in comparison with the previous cases as shown in Fig.12.The increased Lyapunov exponent indicates sensitivity to initial conditions.

Fig.11.Phase space in the presence of white Gaussian noise,(I=2+ξ(t)).The energy feedback obeys Eq.(12).(a)k4=0,(b)k4=1.0,(c)k4=5.0,and(d)k4=10.0.

Fig.12.Variations of LLE with k4.White Gaussian noise is added as the external stimuli.The inserted figure is the enlarged version.

5.Synchronous dynamics

5.1.Electrical coupling

The dynamic equations for electrically coupled HR neuron model with quadratic flux controlled memristor has the form,

where

The parameter geis the coupling strength of synaptic junction and D describes the field coupling strength.The parameters used are a=3.0,b=5.0,R=4.0,r=0.006,xe=-1.61,k1=0.1,k2=0.9,and k3=0.5,α=0.4,β=0.02,and γ=0.1.[9]We have evaluated the effect of modulation intensity of electromagnetic field in regulating the average energy of neurons in the synchronized state.

Figure 13(a)shows the average energy variations for k1=0.1.The plot verifies that with increase in coupling strength the〈H〉changes in a waving pattern and it suddenly stabilizes at the point of synchronization.The membrane potential of both neurons are equal in the synchronized state and this leads to the vanishing of coupling term.As a result the energy in the synchronized state returns to its initial uncoupled value.Our results with quadratic flux controlled memristor are in accordance with the results obtained for HR model without memristor in Ref.[1].The variation in〈H〉for an increased value of k1(k1=0.5)is shown in Fig.13(b).From the plot it is clear that as the value of k1is increased,the onset of synchronization occurs at a low value of ge.The results imply that,an autonomous chaotic system with linear feedback coupling will move to its natural oscillatory regime in the synchronized state by gaining or dissipating energy.If the system continues in the same state,then the change in total average energy will be zero due to the repeating nature of trajectories with arbitrarily close energy values.

Fig.13.Average energy of electrically coupled HR model by varying the coupling strength:(a)k1=0.1 and(b)k1=0.5.

The transition to the synchronized state is further confirmed by plotting the transverse Lyapunov exponents(TLE).[26].Figure 14(a)shows the variations of two of the largest TLEs with gefor k1=0.1.With increase in ge,the largest TLE(λ⊥1)also starts to increase,reaches a maximum,and then starts to decrease.λ⊥1crosses zero at ge=0.41 indicating a transition from desynchronized state to synchrony at this point.

Fig.14.TLEs of electrically coupled HR model by varying the coupling strength:(a)k1=0.1 and(b)k1=0.5.

Figure 14(b)shows the variations in TLEs for k1=0.5.Here λ⊥1crosses zero at a low value of ge(ge=0.37).These results are consistent with the average energy changes shown in Fig.13.

5.2.Chemical coupling

The equationsgoverning the dynamicsofchemically coupled HR neuron model with quadratic flux controlled memristor has the form,

where

The parameter gcis the coupling strength of synaptic junction and D describes the field coupling strength.Vs,the reversal potential is always greater than x for all neuron at all times.For each neuron to reach the threshold,we choose θ=-0.25,Vs=2,and λ=7.5.[33]

The change in average energy with gcfor k1=0.1 is shown in Fig.15(a).From the plot it is clear that the fluctuating nature of average energy disappears at the point of synchronization.This occurs at gc=1.55.After that,the system shows an interesting behavior;i.e.,〈H〉linearly increases with increase in gc.

Fig.15.Average energy of chemically coupled HR model by varying g c for(a)k1=0.1.The green,red,black,and blue lines in the inset plot shows the AD state obtained for g c=1.55,1.7,1.8,and 2,respectively.(b)Average energy variations for k1=0.2.The inset plot represents the time series of membrane potential corresponding to periodically oscillating(red)and AD(black)states.

To unravel the reason for increase in〈H〉after gc=1.55,we have examined the time series of the system in the linearly increasing regime of 〈H〉.The results are shown in the inset plot of Fig.15(a).The ‘X’axis represents time and ‘Y’axis denotes the membrane potential.The plot verifies that,even though both neurons are synchronized(gc≥1.55),their membrane potential is in the amplitude death state(AD).The green,red,black and blue lines show the AD state which has been obtained for gc=1.55,1.7,1.8,and 2 respectively.As the value of gcis increased the value of membrane potential at which AD takes place also increases.Thus,in chemical synaptic coupling where the interaction terms do not go to zero permits a change in average energy in the synchronized state.We have further analyzed the effect of k1in regulating〈H〉.Figure 15(b)shows the 〈H〉variations for k1=0.2.For 0＜ gc≤ 0.55,〈H〉shows some fluctuations and this region corresponds to the desynchronized state.After that,〈H〉remainsconstantfor0.55＜gc≤1.484 and linearly increasesfor 1.484＜gc≤2.0.The time evolution of membrane potential corresponds to these two cases are shown in the inset plot of Fig.15(b).The periodically oscillating(red)time series is for the constant〈H〉regime,i.e.,in the synchronized state there occurs a transition to stable orbit and the total average energy change is zero due to the repeating nature of trajectories.The black line in the inset plot is the time evolution of membrane potential corresponding to the linearly increasing regime of〈H〉.Here also AD occurs at different points for different values of gcleading to a net average energy change.Thus it can be concluded that,an extra amount of energy is needed when a coupled system is forced to oscillate in different regions of the phase space where the average energy change is not zero.The extra demand of energy required for the collective dynamics is provided by the coupling mechanism.[34]In the central nervous system,some specialized structures are located at the postsynaptic sites for producing ATP molecules to balance the energy demands.[35]

Fig.16.TLEs of chemically coupled HR model by varying the coupling strength:(a)k1=0.1 and(b)k1=0.2.

The transition to the synchronized state isverified by plotting the TLEs.Figure 16(a)depicts the variations of two of the largest TLEs with chemical coupling strength for k1=0.1.The λ⊥1crosses zero at gc=1.55 corresponds to the synchronized state.For k1=0.2,this transition occurs at a low value of gc(gc=0.55),as shown in Fig.16(b).

6.Results and conclusions

We have analyzed the energy aspects of single and coupled HR neuron models with a quadratic flux controlled memristor.The bifurcation analysis of single system by increasing the value of I suggest that the response of the system in the presence of external stimuli can be improved by properly modulating the electromagnetic induction.Based on Helmholtz theorem the rate of change of Hamilton energy of HR model with memristor has been derived.It is found that the average Hamilton energy decreases with increase in I.The time evolution of membrane potential and the rate of change of energy function for different external stimuli have been analyzed.In the case of constant external current,the electrical mode of neuron changes with the external forcing and energy is consumed for generating each action potential.In the presence of periodic stimuli,the firing mode changes with the change in amplitude of external forcing.The energy consumption of bursting state with less number of spikes per burst is found to be low compared to burst with more number of spikes.The analysis of the system by applying Gaussian white noise reveals that,bursting mode changes with increase in noise intensity.The bifurcation analysis corresponding to periodic input shows the presence of intermittently occurring states with the variations in A.The bifurcation diagram in the presence of noise reveals that as the noise intensity is increased,the system shows complex chaotic dynamics.The dependence of Hamilton energy function on system parameters is used to control and stabilize chaotic trajectories by giving a negative energy feedback.The suppression of chaotic trajectories and the stabilization of phase space of the system for periodic input and noise are discussed.As the feedback gain in the energy function is increased,the initially positive LLE become negative which in turn ensure the stabilization of chaotic trajectories.

In the case of electrically coupled neurons,as geis increased the〈H〉changes in a fluttering pattern and it stabilizes at the point of synchronization.The energy in the synchronized state returns to its initial uncoupled value due to the vanishing of coupling term in the synchronized state.This study has been repeated by increasing k1and found that for an increased value of k1the onset of synchronization occurs at a low value of ge.The average energy variations exhibit three important regions when the neurons are coupled via chemical synapse.The fluctuating region indicates desynchrony.In the region where the〈H〉remains constant,the system shows synchronization with periodically oscillating dynamicsand the total average energy change is zero due to the repeating nature of trajectories.In the linearly increasing regime,the dynamics are AD.As the value of gcis increased,the value of membrane potential at which the system stabilizes is also increased and leads to a net average energy change.We conclude that if two neurons are coupled and forced to oscillate,then their phase space may contain different oscillatory regimes.As a result,the change in average energy of the system will not be zero and an additional amount of energy is used to sustain the synchronized state.The proposed method will be useful to study the energy aspects of other coupled chaotic and hyperchaotic systems.Possible extensions to neural networks can providefiner insight to the energy modulation mechanism of various biological systems.

Acknowledgments

UK would like to acknowledge University Grants Commission,India for providing financial assistance through JRF scheme for doing the research work.PAS would like to acknowledge DST,India for their financial assistance through the FIST program.