Dynamics Analysis of Neuron Bursting under the Modulation of Periodic Stimulation

A nonsmooth neuron model with periodic excitation which can reproduce spiking and bursting behavior of cortical neurons is investigated in this paper. Based on nonsmooth bifurcation analysis, the mechanism of the bursting behavior induced by slowchanging periodical stimulation as well as the associated evolution with the variation of the stimulation is explored.Themodulating character of the external excitation and the effect of the bifurcation occurring at the switching boundary of the vector field are presented.


Introduction
Recently, understanding of the brain and its behavior has been a research focus for its importance and various applications.As primary building blocks of brain, innumerable neurons compose the very complex information network.It can encode, transfer, and integrate information by firing activities [1].Therefore, deep exploration of single neuron's behavior is necessary to reveal how brain works.
The neuron is one of the most sophisticated nonlinear dynamical systems.A series of dynamical models which can describe the corresponding action potential shapes for different neuronal cell types has been developed [2][3][4].In general, these models can be assigned to two classes [5].One category are conductance-based models with high dimensional and high biological precision such as Hodgkin-Huxley and Morris-Lecar models.These models can imitate most experimental measurements to a high degree of accuracy, but they are mostly complicated and difficult to physically implement [6].The second are one-dimensional integrateand-fire (IF) type models, which are often advocated since they are easy to analyze.However, they do not quite capture the dynamics of a truly excitable system with gating variables.Thus it is natural to search for planar models possessing one voltage and one gating variable that can mimic the behavior of high dimensional conductance-based models [7].The FitzHugh-Nagumo model is a typical example.And piecewise linear (PWL) nullcline is introduced inspired by this to build a series of PWL models that can describe the dynamical behavior of many common cell types, such as rich spiking and bursting dynamics [8,9].Although many conclusions have been achieved, most work is confined to the autonomous case.Period sensory irritation is common in natural environment yet [10], and the periodical stimulation is confirmed to have an effect on information handling in the cortex on larger scale compared with aperiodic one [11].
To explore the influences of periodical stimulation on the dynamical behavior of the neuron, a two-dimensional PWL neuron model is discussed in this paper.We tackle the case with slowly varying periodic excitation and investigate the evolution of the bursting dynamics based on bifurcation analysis.

The PWL Model and Periodic Bursting
This model is developed from a PWL modification of the Izhikevich model [1], which approximates the quadratic part of the Izhikevich model with two crossed lines [12].This approximation can be formulated as where V represents the membrane potential of the neuron and u represents a membrane recovery variable, which accounts for the activation of K + ionic currents and inactivation of Na + ionic currents and it provides negative feedback to V.
Here, we focus on the case with a periodic stimulation; that is, the stimulation in model ( 1) is set as  =  0 cos( 0 ).In order to reveal the dynamical behaviors of this system under the influence of stimulation, the parameters are fixed at  = 1.8,  = 2.06,  1 = 2.8,  2 = 3.0, and  3 = 7.5, and the value of the angular frequency of the stimulation is set at  0 = 0.02.With the amplitude of the simulation being seen as a modulating parameter, the periodic bursting for  0 = 1.0 and  0 = 1.4 can be observed by numerical simulation (seen in Figure 1), the oscillating amplitude of which is noticed to be amplified with the increase of  0 and the duration time of spiking state is longer in the case that  0 = 1.4.
Remark 1.The transient states in these two cases are both eliminated.
Demonstrated as the red curve in Figure 1(a), timedependence of the stimulation is superimposed on the time series of V to investigate the internal connection.It could be seen obviously that the transition between a rest state and a spiking state of the model is closely related to the variation of the periodic stimulation.In other words, the variable V in system (1) oscillates with two frequencies and the value of the small one depends on the frequency of the periodic stimulation.In the following, mechanisms of the behavior presented above will be discussed quantitatively based on bifurcation analysis.

Nonsmooth Bifurcation Mechanism
System (1) is a nonautonomous system since the value of  changes periodically with .Based on dynamical system theory, any -dimensional nonautonomous system, that is, ẋ  =   (,  1 , . . .,   ),  = 1, . . ., , can be transformed into a ( + 1)-dimensional autonomous system which is described as system (2) by introducing a variable quantity  = () = : (2) Thus, system (1) could be written in an autonomous form, namely (3)-( 5), by introducing the transformation  0  = , where 0 <  0 ≪ 1 for the slowly varying periodic excitation: Since the value of  0 is small, ( 5) represents dynamics of a relatively slowly changing process while (3) and ( 4) describe the relatively fast ones.The full system can be divided into slow and fast subsystems.The slow subsystem is given by (5), while the fast subsystem is given by ( 3) and (4).Thus  and even  =  0 cos() could be treated as a slow-varying parameter of the fast subsystem.That is, the fast subsystem is an autonomous system at a certain moment which could be described as follows: where For this PWL autonomous system, one switching boundary exists, denoted by ∑ : V = − 2 .Thus, two cases of equilibria can be obtained for a given value of  0 , which can be expressed as Further analysis indicates that the number of equilibrium points is not constant under the above parameters due to piecewise linearity of (V).Specifically speaking, there is a critical value of  0 ; that is, two equilibria could be observed for  0 < 1.32, only one equilibrium could be found for  0 = 1.32, and no equilibrium could be obtained for  0 > 1.32.For example, Figure 2 shows the distribution of the equilibria with a few different values of  0 .The stabilities of these equilibria of the associated autonomous systems are determined by the eigenvalues of the related linearization matrix.Here, we take the case of  0 = 1.0 as an example.The equilibria are calculated as  1 {−3.0658,−6.3156} and  2 {−2.5676, −5.2892}.The stability of  1 and  2 is determined by the eigenvalues of  1 and  2 , respectively, which are expressed as follows with the above parameters: with eigenvalues   1 = 1.7578 and   2 = −0.7578.Obviously, the equilibrium point  1 is a stable focus, while  2 is an unstable node.Now we focus on the full system (3)-( 5) which is equivalent to system (1) describing the nonautonomous case.From the analysis above, though the initial system (1) is a nonautonomous system, it can be regarded as an autonomous one at any certain moment and the associated solution for which at the moment still could be described as formula (7).What is special is that the value of  0 in ( 7) is variable periodically and the solution will vary with it.That is, the full system, as well as system (1), could be regarded as a "generalized autonomous system" and the related solution could be regarded as a "generalized equilibrium point" as well [13].Next, we investigate the evolution mechanism of this periodic bursting via bifurcation of these "generalized equilibria." Due to the linear structure of the system, though nonsmooth, the generalized solution for system (1) can be assumed as the following form and the related coefficient can be calculated by substituting this "solution" into system (1): The details can be described as follows: (1) For V < − 2 The "solution" expressed by ( 10) is just the "generalized equilibrium point" of the "generalized autonomous system" mentioned above, the respective stability of which shows a coincidence with that of the genuine equilibrium of the associated autonomous system because the expression of the Jacobian matrixes of the solution is irrelevant to the value of .It means that the "solution" is a stable focus for V < − 2 , while it is an unstable node for V > − 2 .The case of  0 = 1.0 is still taken as an example.The distribution of the "solutions" for system (1) is presented in Figure 3, where the red segment in region (I) denotes the stable "solutions" while the pink segment in region (II) represents the unstable "solutions."The green line   is the nullcline of .Above (below)   , we have u  < 0 ( u  > 0).
Remark 2. The set of the "generalized equilibrium points" is not straight lines but piecewise curves (seen in the magnification of the curves).
It should be pointed out that the classical continuous Jacobian matrix cannot be obtained due to the nonsmoothness of the vector filed of this system.According to differential inclusions, we can use the generalized differential of Clarke to set up a "generalized Jacobian matrix" to explore the bifurcation of the equilibrium at the switching boundary and the "generalized Jacobian matrix" of system is expressed by [14].The eigenvalues of   (), denoted by  1,2 , are set-valued and form a path in the complex plane with  as path parameter.
Figure 4 described the path of eigenvalues of the generalized Jacobian at the switching boundary ∑ : V = −3.0 for  0 = 1.0, from which it may be found that the eigenvalues of   () are purely imaginary for  = 0.8214.The path of the eigenvalues of   () shows that discontinuous Hopf bifurcation occurring at the switching boundary [14].And the oscillating frequency generated from this bifurcation is indicated at   = 0.665.Now we explain the mechanism of the periodic bursting based on nonsmooth bifurcation theory.For this purpose, the phase portrait shown in Figure 1 is projected onto the distribution diagram of the "generalized equilibrium" of the nonautonomous system demonstrated in Figure 3 (shown in Figure 5).
In Figure 5, the stable "generalized solution" located in region (I) is observed to match the phase portraits exactly while the "solution" in region (II) does not agree with the portraits.The "generalized equilibrium points" are stable ones in region (I), and they lose their stability at the switching boundary ∑ with the variation of the value of stimulation .Thus the trajectory diffuses and oscillates with the frequency generated from the discontinuous Hopf bifurcation (demonstrated in Figure 4), which is much higher than that of stimulation.In other words, this bifurcation leads the system transfer from rest state to spiking state.Based on the qualitative analysis, the oscillation of the system in one period is described in detail quantitatively as follows.
Point  located on the portrait in region (I) is set as the starting point, which is above the nullcline of .At this point, the values of  will decrease since u  < 0. It means that the trajectory runs downward at point .The anticlockwise direction of the trajectory can be demonstrated by the time series plotted in Figure 1.
The orbit may move down exactly along the stable "solution" with the variation of the stimulation, that is, the angular frequency  0 = 0.02, until it turns around as the stimulation reaches its extremum value.Then the trajectory continues to move up along the "stable solution" with the variation of .When the orbit crosses the switching boundary ∑, the stable focus loses its stability through the discontinuous Hopf bifurcation occurring at the boundary.Thus the trajectory no longer follows the unstable "generalized solutions" in region (II) but oscillates with the frequency generated from the bifurcation; namely,   = 0.665, which approaches the value   = 2/11 ≈ 0.571 calculated from the numerical simulation (shown in Figure 1(a)).Note that this passage is much faster relative to that described previously.It can be called spiking state while the passage along the "stable solutions" is the slow rest state.At the same time, the value of stimulation still varies periodically with its own frequency during this oscillation process.Thus, the orbit converges with the slowly varying stimulation until it meets the stable "solution" at the switching boundary ∑.Then the trajectory moves along these "solutions" slowly again and another period starts.
We now account for the evolution of the phase portrait with the increases of the stimulation amplitude.
For the purpose of comparison, Figure 6 presents the partial enlarged structures of the phase portraits for  0 = 1.0 and  0 = 1.4.For description convenience, a pair of most remote intersections of the trajectory and the switching boundary, that is,   and   ( = 1, 2), is introduced here.As the analysis presented above, the stability and dynamical properties of the "generalized equilibrium points" are irrelevant to the external stimulation.That means that the eigenvalue of the associated equilibrium points for  0 = 1.4 remains the same as that for  0 = 1.0.It is worth pointing out that, with the increase of the stimulation amplitude, the "equilibrium points" may disappear when the value of the stimulation exceeds the threshold 1.32.It leads to the divergency of the orbit to some extent.Due to the periodicity of the stimulation, the trajectory still could converge to the "stable solution" again with the stimulation varying.Thus the oscillation amplitude for  0 = 1.4 is enlarged.Thus, the distance from  2 to  2 is larger than that from  1 to  1 with the increase of  0 .On the other hand, the imaginary parts as well as the real parts of the complex eigenvalues keep constant for different  0 .It means that the oscillating frequency of the trajectory as well as the related convergence rate is constant.It is exactly the reason

Figure 2 :
Figure 2: Distribution of the equilibria of system (6) with different values of  0 .