Dynamic Behaviors in Coupled Neuron System with the Excitatory and Inhibitory Autapse under Electromagnetic Induction

The induced current produced by electromagnetic induction can adjust the membrane potential of neuron through the feedback of a magnetic flux-controlled memristor. We adopt the numerical simulation method with the aim of investigating the synchronous behavior in the neuronal system that is coupled by chemical and electrical synapses under electromagnetic induction. Within the improved model, the effects of electromagnetic induction on neurons are described with additive memristive current on the membrane variable, and the memristive current is dependent on the variation of magnetic flow. The simulation results show that the two coupling modes play an important role in the synchronization of the system. By increasing the chemical synaptic feedback gain, we observe a transition from mixed oscillatory to periodic state at a critical value. In addition, two Hopf bifurcation points are found with the change of the external stimuli, and the state of neuron discharge is influenced by initial values. Furthermore, there is a domain of coupling strength and feedback gain values, in which the two coupled neuron system is synchronized and longer time lag is not conducive to the system synchronization.


Introduction
A neural system, which is made up of a large number of neurons, is a complex information network.Different types of discharge patterns can be switched under the control of external stimulation or bifurcation parameter.In order to understand the regulating function of the nervous system, many models of neuronal electrical activity have been proposed.Commonly used models include the FitzHugh-Naguma model [1], Morris-Lecar neuron model [2], Hindmarsh-Rose model [3,4], Nagumo-Sato neuron model [5], and Wilson-Cowan neuron model [6].These models that describe neuron dynamics with a set of differential equations are almost derived from the Hodgkin-Huxley [7] model or the simplified version.Some results from biological experiments [8][9][10] can be explained by theoretical neuron models, such as the Morris-Lecar neuron model.In this model, the membrane potential exhibits quiescent, spiking, or bursting state by changing the external forced current [11].Neurons do not work in isolation, but they interact to affect the processing of information.There are two forms of synaptic coupling found in the real nervous system, namely, electrical synapse and chemical synapse.The synchronization phenomena are a typical manifestation of the rhythms of group movement; that is, all neurons in the system have a certain connection at the same time or rhythm [12][13][14].Bazhenov et al. [15] designed a coupled linear chain of Hindmarsh-Rose model neurons with reciprocal inhibition between neighboring neurons that exhibited synchronous oscillations.Zhang et al. [16] proposed a class of synchronization problems of nonlinear time-delay dynamic networks with a nonuniform impulse effect.Burić et al. [17] studied the synchronization of Hindmarsh-Rose neurons with a time-delayed fast threshold modulation synapse.Xu et al. [18] analyzed the synchronization behavior and mode selection in neural networks under the coupling of chemical or electrical synapses.Yao et al. [19] investigated the influence of coupling strength, time delay, and network topology on synchronization behavior in delay-coupled networks of chaotic pendulums.Gokul and Kapitaniak [20] studied the synchrony of coupling multistable systems which have hidden attractors with each other.In coupled oscillators or coupled neurons, synchronization may occur because of the appropriate coupling effect [21][22][23][24].Interestingly, the stochastic and coherence resonance [25][26][27] of the nervous system is induced by appropriate noise intensity and external periodic stimulus.The synchronization of the coupling system is an interesting research filed.It is challenging to analyze the dynamic mechanism caused by the variation of the coupling parameters and modes of the system.The synchronization phenomena in Hindmarsh-Rose (HR) neurons that are connected by electrical coupling and chemical coupling, moreover, complete synchronization, phase synchrony, and antisynchrony of neurons are realized [28,29].The neural electrical activity has also been widely studied and verified in the circuit [30][31][32][33][34][35].For example, Vaidyanathan et al. [36][37][38] designed electronic circuits to study the feasibility of the 3D novel jerk chaotic system with hyperbolic sinusoidal nonlinearity.Conti and Turchetti [39] performed a circuit to realize approximate identity neural network for the analog synthesis nonlinear dynamical system.Pham et al. [40] proved the existence of chaotic behavior in a threedimensional autonomous chaotic system with a circular equilibrium by using OrCAD PSpice software and experimental.
It is necessary to study the effects of electromagnetic induction on neuronal cells [41][42][43][44].The changes of membrane potential can induce electromagnetic induction between neurons.As reviewed in [45,46], the effects of electromagnetic radiation in Homo sapiens include electrical activity of neurons, energy metabolism, genomic responses, neurotransmitter balance, blood-brain barrier permeability, cognitive function, sleep, brain tumors, and other encephalopathy.Lu et al. [47][48][49] investigated the effects of highand low-frequency signal stimulus on neural activity under electromagnetic radiation.According to Faraday's law of induction, the magnetic field is a result of fluctuations in the action potential.That is, the distribution of electromagnetic field both inside and outside neurons can be changed by the fluctuation of the membrane potential.Therefore, a new three-variable ML neuron model is established by introducing an additional variable as magnetic flux which adjusts the membrane potential via a memristor [50,51].
The following study is based on the proposed Morris-Lecar neuron model with consideration of magnetic flux, in which the dynamic characteristics of the neurons are studied by using bifurcation diagrams and time series of the discharge.A preliminary synchronization analysis was conducted in the excitatory and inhibitory neural system.The study revealed that excitatory and inhibitory neurons can be synchronized under the appropriate coupling strength.The synchronization behavior of the system is also affected by the time lag when the coupling strength and the feedback gain are maintained.

Model and Scheme
The Morris-Lecar (ML) equations were originally developed as a mathematical model of muscle fiber.For the neuron, the effect of electromagnetic induction should be considered during the discharge process of the membrane potential.The electric activity will change because of the fluctuation of electromagnetic induction and ion concentration in the process of ion exchange.We modify the basic ML model, including the impact of the electromagnetic radiation.The improved ML neuronal model [44] contains three variables, and the dynamic properties are described as follows: where V and ω denote the variables for the membrane potential (mV) and gate channel, respectively.Parameter c is the capacitance of the membrane (μF/cm 2 ).The g Ca , g K , and g L denote the maximum conductance (mS) of calcium ion, potassium ion, and leak ion, respectively.V Ca , V K , and V L are the reversal potential (mV) corresponding to these channels.m ∞ V and ω ∞ V define the value of the opening probability for the calcium ion channel and the potassium ion channel in the steady state, where V 1 , V 2 , V 3 , and V 4 are the parameters of the steady system, and λ ω V defines the rate constant for the opening of potassium ion channel.The parameter ϕ is marked as the variation between the fast and the slow scales of neurons.As described in [52,53], the variations of the intercellular and extracellular ion concentration can induce electromagnetic induction, which can be expressed by magnetic flux according to Faraday's law of electromagnetic induction.The induced current produced by electromagnetic induction can adjust the membrane potential by the feedback of the memristor.The memristor in model (1) can be divided into two ways: the charge controlled and the magnetic controlled.For the potassium ion-channel memristor, the second term in the right of (1) can be rewritten as i and G K is the potassium memductance function.The fourth term in the right of (1) can be rewritten as i φ = G φ φ v φ with v φ ⇔V and G φ ⇔kρ φ , which defines another first-order memristor, and the conductance value of the memristor depends on the input current.The expression of ρ φ = α + 3βφ 2 denotes the memory conductance of a magnetic flux-controlled memristor [54], it is used to calculate the effect of feedback regulation on the 2 Complexity membrane potential when the magnetic flux is changed, and α and β are fixed parameters.Therefore, as in [54], the induced current and electromagnetic induction can be described by The variable i′ represents induction current.The term − kρ φ V represents the inhibitory modulation of membrane potential, and it describes the induced current induced by electromagnetic induction.The parameter k is the induction coefficient, and its value depends mainly on the medium itself.I ext is the external forcing current.The terms k 1 x and k 2 φ in the (1) mean the influence of membrane potential on magnetic flux and leakage of magnetic flux, respectively.
For the analysis of the possibility and stability of the synchronized dynamics between two neurons under bidirectional coupling, the dynamic equations are given by where the subscripts α and β are a pair of coupled ML neurons under electromagnetic radiation.C denotes coupling intensity between adjacent neurons.In order to simulate the chemical synapse feedback of neurons, we shall use the so-called fast threshold modulation scheme proposed by Somers and Kopell [55] and often used by others, for example, [56,57].This chemically feedback form, which clearly combines the time lag of the synapse, is provided by the following functions: The variable parameter H syn is the feedback gain at time t with itself connected at time t − τ.The symbol τ indicates the time lag (ms) of the signal propagation.C is coupling strength between two neurons.V syn represents the synaptic reversal potential (mV), which depends on the presynaptic neurons and receiver.The chemical coupling is characterized by the difference between the synaptic reversal potential and the synaptic potential.A positive or negative sign of the difference corresponds to an excitatory or inhibitory effect of the synapse.If the synapsis is excitatory, V syn = 15 mV, and if the synapsis is inhibitory, V syn = −10 mV.The parameter θ is a synaptic threshold.Considering that the neuron membrane potential value of the improved ML model is between −17 mV and 15 mV, θ = 4 mV is selected to ensure that the spike of the V is over the threshold, and the quiescent state of the V is less than the threshold.That is, the membrane potential of the presynaptic neuron is more than θ, and it can play a role in the postsynaptic neuron [58,59].σ is the ratio constant to the start of excitement or inhibition.In this paper, we focus on the collective behavior of the two coupled neuron system driven by the excitatory and inhibitory autapse, and the schematic diagram is shown in Figure 1.Parameters of the improved ML neuronal model are given as c = 20 μF, V Ca = 120 mV, To characterize the synchronization in the system of coupled spiking neurons, a method of calculating the error function is introduced in the following [24]: Equation ( 6) indicates that the lower the value of e corresponds with the better synchronization in the system.

Results and Discussion
In this section, first of all, the bifurcation is theoretically analyzed to reveal the dynamic mechanism of the discharge mode in the improved model (1).Then the fourth-order Runge-Kutta method is used to calculate the improved ML neuronal model, the step of time h is selected as 0.01, and the transient period for calculating is 5000 time units. Let If the parameter k = 0.1, the system (8) can be written in the form , 3 Complexity and thus, Obviously, the expression of the equilibrium point (V s , ω s , φ s ) can be obtained by (10).Then we analyze the stability of the equilibrium point with I ext as the bifurcation parameter.It is noted that the stability of the equilibrium point is determined by the eigenvalue of its Jacob matrix; that is, when all the eigenvalues are a negative real part, then the equilibrium point is stable; otherwise, it may be marginally stable or unstable [60].The appearance of a pair of pure imaginary eigenvalues signifies Hopf bifurcation [61].
At the equilibrium point (V s , ω s , φ s ), the linearization Jacobi matrix of the improved model is that Therefore, the characteristic determinant of this system at the equilibrium point is where 13 If the Jacobi matrix has one eigenvalue of the negative real part and two zero real parts at a critical value of the bifurcation parameter, then we usually say that a Hopf bifurcation occurs [61].According to the relationship between λ and I ext , we calculate the branch of the equilibrium point which undergoes two Hopf bifurcations at parameter I ext1 = −27 99 and I ext2 = 23 79.Obviously, they correspond to the bifurcation points HB in Figure 2(b).Similarly, the corresponding bifurcation points can be obtained by changing the parameter k.
The results of the bifurcation analysis for the system with the electromagnetic effect described by (1) are numerically simulated.In the bifurcation analysis, the influence of the external forcing current I ext on the discharge behavior of the neuron under different electromagnetic effects known as induction coefficient k is investigated.The system will have the process of "resting-exciting-silent" by increasing the external forcing current.In Figure 2, several typical bifurcation diagrams for the different external forcing currents, without (k = 0) and with (k ≠ 0) the electromagnetic effects, are plotted, respectively.
Figure 2 shows the bifurcation with I ext as the bifurcation parameter for six different values of the induction coefficient.The equilibrium point of the system has undergone five changes without electromagnetic radiation, as shown in Figure 2(a).In the beginning, the system has only a stable equilibrium point, and the action potential of the neuron will As the induction coefficient k increases, a similar phenomenon is observed and the corresponding bifurcations are summarized in Table 1.
It can be noted from Table 1 that parameters describe the interaction between membrane potential and magnetic flux which is further increased to k = 0.6, the threshold of excitability is increased to a higher value corresponding to I ext = −11.39μA (see Figure 1 and Table 1).However, the position of the second bifurcation point was reduced from 23.55 mV to 9.0 mV.That is to say, with the increase of the induction coefficient k, the region of the limit cycle in the system is gradually compressed.It should be pointed out that the system has a Hopf bifurcation point HOPF1 with a negative value of parameter I ext (see Table 1), which will not be considered in the following text since it loses the biophysical meaning.At the same time, in the numerical simulation, we find that when I ext > HOPF2, the selection of the initial value of the system is as important as the external current.Based on the analysis of the bifurcation diagram of V max and V min for the improved ML neuronal model, the external forcing currents I ext = 40 μA, I ext = 50 μA, I ext = 60 μA, I ext = 70 μA, and I ext = 80 μA are chosen to impose on the neuron.The type of membrane potential discharge at different initial values is calculated in Figure 3.
The results in Figure 3 show that the initial values of V 0 and ω 0 have a great influence on the discharge of neurons in the improved ML model.The oscillating areas of the membrane potential are reduced by the increase of the external stimulation current; that is, the red area in the picture becomes smaller.The oscillating region of the membrane potential is found to be banded, which indicates the changes in the initial value of φ 0 with little effects on the discharge of neurons as shown in Figure 3, a2-e2 and a3-e3.
The results in Figure 4 show that the initial value region that produces the oscillations of the membrane potential is 5 Complexity reduced by the increase of the external stimulation current.Most interesting, the study found that the selection of the initial value will also be influenced by the induction coefficient k under the same external stimulus current.This discovery is in accordance with the bifurcation diagram of Figure 2. Therefore, it is very important to choose the appropriate initial value with the different induction coefficients k.According to the initial value area shown in Figures 3 and 4, sampled time series for membrane potential and phase portraits are plotted in Figure 5.
The improved ML neural model can exhibit several kinds of oscillations.When the initial value is determined, the electric activity depends on the external forcing current.The change of electrical activity between the quiescent state and the spiking state can be observed by selecting different external forcing currents.The phase portraits of the external forcing currents (μA) are chosen as 0, 10, 20, 30, 40, and 80 in Figure 5 The numerical results in Figure 5(a) show that the regions of limit cycles become larger with increasing the external forcing current.There exist some thresholds in the system, which determines the conversion of the discharge mode of the neuron membrane potential.In the absence of external excitation, the neuron can still be discharged under the electromagnetic effect, as shown by the blue curve in Figure 5(b).With increasing the external forcing current, the amplitude of the neuron membrane potential will increase gradually; in other words, the amplitude of the periodic oscillation is related to the area size of the limit cycle.When the external forcing current is increased to a certain threshold, the system will change from the oscillating state to the resting state, as depicted in Figure 5(c).
Meanwhile, we studied the collective behavior of the two neurons driven by the excitatory autapse and inhibitory autapse in the case of electric coupling.The initial values are V 0 = 100 mV, ω 0 = −1.5 μA, and φ 0 = 0.1, and the external forcing current I ext = 40 μA is chosen for its simplicity; the induction coefficient is k = 0 1.At first, the inter-spike interval of the β-neuron membrane potential (as in ISI representation) at different feedback gain H syn is calculated, and the results are plotted in Figures 6 and 7.
Presented results show clearly that feedback gain H syn as well as coupling strengths plays an important role in the modes of electrical activities.When the value of H syn is smaller, the rich discharge modes are observed in Figure 6.The value of the bifurcation point is also affected by the coupling intensity.It is found that the mode of β-neuron discharge does not exhibit periodicity, unless larger feedback gain could be applied.In essence, there may be two types of autapse in the system, and the increase of H syn makes the difference in the discharge of the two neurons increasing; that is, (a1)  6 Complexity the inhibitory neurons tend to be quiescent and hence cannot affect the excitatory neurons.However, it can be observed that the modes of electrical activities depend significantly on the values of the synaptic feedback gain.
Within a certain coupling strength, if the error e tends to zero with the time increased, the coupled neurons are fully synchronized.According to the analysis of the results of Figure 6, the sampled time series for membrane potential  7 Complexity and phase portraits of the coupling strength C = 0.5 are calculated with different feedback gains H syn ; the results are plotted in Figure 8.
The phase portraits and time series of the membrane potential of two coupled neurons are illustrated in Figure 8.The limit cycle is shown in Figure 8(a1), indicating that the β-neuron is periodic discharge.The region of the two limit cycles is different in Figure 8, (a2), and the value of the V is within the range of −20 mV to 20 mV.The phase portrait of (V β , V α ) is located near the corner line of the first quadrant, which means the occurrence of approximate synchronization.The error e (blue line) is found to exhibit periodic oscillations in Figure 8, (a3).Interestingly, when the coupling strength C is further increased, the phase portrait of (V β , V α ) tends to have a straight line, and the results are shown in Figure 9.
The results in Figure 9 confirmed that the phase portrait of (V β , V α ) coincides with a straight line, which is located on the angle bisector of the first quadrant.With appropriate time lag and feedback gain, the two coupled neuron system will synchronize with the increase of the coupling strength between neurons.For the improved model in this paper, the numerical results show that phase synchrony can be achieved by selecting the appropriate coupling intensity of   8 Complexity the two neurons with electromagnetic radiation.The effect of time lag in autapse should also be considered; for example, the time lag is increased to 100 ms, and some results are found in Figures 7 and 10. Results presented in Figure 7 reveal that the oscillatory pattern is largely influenced by the synaptic delay.For suitably long values of τ, a complex oscillatory pattern can be observed.Interestingly, however, if the H syn is sufficiently high, we can observe the emergence of a periodic firing, which implies that there is a transition from chaos to periodic discharge in the system.Theoretically, under sufficiently long synaptic delay condition, neurons have enough time to fire more than once during a whole periodic cycle, before the synaptic currents caused by the first synchronous spiking within the same periodic cycle start to affect their firing [62].
It is found in Figure 10 that spiking and bursting discharge behaviors of β-neuron reappear depending on the gain and delay of the autapse.The error e is observed to increase obviously through the time lag from 50 ms to 100 ms.But what is more interesting is that the increment of τ makes the amplitude of the membrane potential of αneuron (excitatory) decrease in Figure 10, (a3).This trend may be that the role of the autapse is suppressed in the appropriate feedback gain and time lag.Therefore, the time lag plays a crucial role in the dynamics of the coupled system.Numerical studies on the synchronization of the two coupled   The trajectories when H syn = 0.03 indicate that the system of the two coupled neuron system is in the oscillation in Figure 11.When the time lag is 0.1, the membrane potential error e of the neurons is smaller.Although the electrical coupling plays a dominant role in the synchronization of the system, the feedback gain from the synapse and time lag are equally important.The increment of time lag is not conducive to the synchronization of the coupling neuron system.This conclusion is consistent with the result of Figure 10.

Conclusions
In this paper, the dynamics of the improved Morris-Lecar neuron model under electromagnetic induction were investigated using bifurcation diagrams and time series of discharge; phase portraits of the neuron under different conditions are investigated in a numerical manner as well.By analyzing the simple numerical simulation of the improved model, the basic dynamic behaviors are obtained by introducing an external forcing current.In the case of the electromagnetic induction, the mechanism of neuron firing has been changed.That is, two Hopf bifurcation points Comparing these results with a previous work [33,58], the bifurcation diagram has an obvious difference due to the consideration of electromagnetic induction based on the ML neural model.In fact, the fluctuation in membrane potential and signal propagation in the neuronal system can generate an induced electrical field and additive current in the media due to electromagnetic induction.As a result, the membrane potential of a neuron can be adjusted slightly by induction field and induced current associated with the variation of magnetic flux.By analyzing the interspike interval series of neural firing, we find that the improved model can generate electrical activity with multiple modes.These results are consistent with the observation observed in the experiment [8].Meanwhile, the preliminary synchronization analysis of a system of excitatory and inhibitory neurons was conducted.
In this aspect, it was unveiled that the neurons in the system can be synchronized by selecting an appropriate coupling strength.A longer time lag is not conducive to the system synchronization, and the higher the feedback gain H syn and the longer the time lag τ are, the more obvious the electrical mode changes in the two coupled neuron system; this conclusion is in accordance with previous experiments [63].Synchronization phenomena are associated with either brain functions [64] or pathological brain states in the neural system.For example, Stam and Bruin used synchronization likelihood to characterize statistical interdependencies between EEG and MEG (magneto encephalography) signals in early and mild Alzheimer's disease [65,66].Rubchinsky et al. [67] presented extensive experimental documentation of the relevance of synchronized oscillations to motor behavior in Parkinson's disease, and they confirmed that the real pathological state is not completely synchronous but showed a complex weak synchronization and highly intermittent dynamics.These results could provide potential theoretical supports for the treatment of neurological diseases.

12 V syn = 15 Neuron 10 Figure 1 :
Figure 1: A schematic diagram of the coupled neuron system.
(a).And Figures 5(b) and 5(c) show the time series of the membrane potential for I ext = 0, 40, and 80.

Figure 9 :
Figure 9: Phase portraits and sampled time series of membrane potential are calculated at coupling strength C = 20, feedback gain H syn = 0 03, and time lag τ = 50 ms.

Table 1 :
Summary of the type and position of bifurcation point.