THE MIXED-MODE OSCILLATIONS IN AV-RON-PARNAS-SEGEL MODEL

. Mixed-mode oscillations (MMOs) as complex ﬁring patterns with both relaxation oscillations and sub-threshold oscillations have been found in many neural models such as the stellate neuron model, HH model, and so on. Based on the work, we discuss mixed-mode oscillations in the Av-Ron-Parnas- Segel model which can govern the behavior of the neuron in the lobster cardiac ganglion. By using the geometric singular perturbation theory we ﬁrst explain why the MMOs exist in the reduced Av-Ron-Parnas-Segel model. Then the mixed-mode oscillatory phenomenon and aperiodic mixed-mode behaviors in the model have been analyzed numerically. Finally, we illustrate the inﬂuence of certain parameters on the model.


1.
Introduction. Lobster cardiac ganglion consists of nine neurons, five large motor neurons at the anterior end and four smaller neurons at the posterior end [5]. The number of neurons in the nervous system of the heart of the lobster is fewer and axon diameter is larger and their enormous neurons in function and chemical constitution are quite easy to identify. Owing to these advantages, the nervous system has aroused people's attention, and they are often used to study on synaptic mechanism models concerned with learning and memory [25,17,8,31].
Some lobster cardiac ganglion models which are based on various ion channels which are dependent on the voltage and the concentration of the ions in the cell membrane are modeled by Av-Ron's team. Different types of firing patterns have been gotten by changing parameters in models. They took the mathematical model with minimal ion channels to simulate the neural activity.
Models of Av-Ron's team are built on the Rinzel model [20]. The initial model based on conductance is the Hodgkin-Huxley model which was obtained by studying the firing behavior of the squid giant axons [13]. It contains three gating variables, the sodium ion activation gating variable m, inactivation gating variables h, and the potassium ion activation gating variable n. Rinzel simplified the Hodgkin-Huxley model by combining h and n to recovery variable w. Av-Ron et al further simplified the model and adjusted the parameters so that the model could better measure the lobster cardiac ganglion neurons. It is called the minimal biophysical model [1]. Furthermore, the bursting pattern can be produced from the model by adding a calcium ion channel, and the new model is called the basic bursting model, and it is also called the Av-Ron-Parnas-Segel model,or the ARPS model for short [2]. On the basis of the two models, the lobster cardiac nerve network model can be constructed [3,4]. Considering the influence of a transient potassium current on both models, the rich dynamic properties of the models mentioned above have been researched, but it seems that the more complex firing patterns such as mixed-mode oscillations have not been presented much in the literature.
Mixed-mode oscillations (MMOs) consist of small amplitude oscillations and relaxation oscillations. The patterns were first found in the Belousov-Zhabotinsky reaction [12], and they had been observed in experiments and models of chemical and biological rhythms [7], and their relation to canards was first proved by Milik et al [16]. A geometric analysis of a basic mechanism producing MMOs in a simple model of a chemical oscillator is presented by the geometric singular perturbation theory [10,14] and canard solutions [30]. A more detailed analysis of MMOs and generation of the canard phenomenon was done by Wechselberger and Brøns et al [28,6,27,9]. Using the geometric singular perturbation theory, especially the blow-up technique [26], some conditions associated with the existence of canards and relaxation oscillations in singularly perturbed systems with one fast and two slow variables were presented, and the geometric explanation of relaxation oscillations and canards was studied [30,29]. In addition, Rubin et al used the same theory to study the MMOs patterns in the Hodgkin-Huxley model [22,23].
In this paper, based on the same method in [6,22], our goal is to understand why the MMO phenomenon of the Av-Ron-Parnas-Segel model appears and explore the influence of different ionic channels and reversal potentials. In Section 2, the specific model and analytical methods are presented. In Section 3, we investigate the mechanism of MMOs in the reduced ARPS model, exhibit the simulation results of the ARPS model and study the bifurcation phenomenon in the ARPS model. Finally, we summarize our conclusions.

2.
Model and methodology. The model that we discuss in this paper includes five types of currents: an inward sodium current I N a , an outward potassium current I K , a calcium-dependent potassium current I K(Ca) , an inward calcium current I Ca and a leak current I L . The Av-Ron-Parnas-Segel model has four dynamical variables: V, w, x and c. V is the membrane potential. w is called a recovery variable, and it is a linear combination of two variables h and n in the Hodgkin-Huxley model.
x is the calcium activation variable which induces the calcium channels to remain open for a long time between firing periods. c is the intracellular concentration of calcium ions. The model is governed by the following equations: where C m is the membrane capacitance and I is the applied current. K p denotes the conversion rate constant from current to the concentration, and is taken as 0.00025 µM · µA −1 . R is the removal rate constant, and is taken as 0.0045. An inward sodium current I N a and an outward potassium current I K are from the Rinzel model [20]. That is, . A calcium-dependent potassium current I K(Ca) is described as the following equation [19]: The model of an inward calcium current I Ca is where V N a , V K , V Ca and V L are the reversal potentials for N a + , K + , Ca 2+ and leakage ions, respectively. g N a , g K , g Ca , g K(Ca) and g L are the maximal conductances for the voltage-sensitive N a + channel, the voltage-sensitive K + channel, the Ca 2+ -sensitive K + channel, voltage-sensitive Ca 2+ channel and the leakage channel, respectively. c K d + c denotes a saturating function of intracellular calcium and constant K d = 0.5. w is a linear combination of gating variables n and h in the Hodgkin-Huxley model Here h 0 and n 0 represent the quiescent values of h and n. s is taken as 1. m ∞ , w ∞ and x ∞ denote the sodium activation steady-state, the recovery steady-state and the calcium activation steady-state, respectively. τ w and τ x are relaxation time.
The explicit expressions are given by the following equations: where V (·) 1 2 denotes the voltage for the half-maximal value. The parameter a (·) represents the slope of the curve at this midpoint. λ is a temperature factor. Parameters of the kinetics equations are listed in Table 1 and 2.
The model was simulated by Matlab 2013a and Matcont, and all charts were plotted by using Matlab 2013a and Mathematica 9. The fourth-order Runge-Kutta method is selected to solve the ARPS models and the step size is taken as 0.005. Value Main results. Two types of firing patterns in the ARPS model are discovered by adjusting parameters of the K + channel and τ x . And on this basis, we further change other parameters to give more firing patterns in the model (see Fig. 1). Table 2 gives corresponding parameter values. In Table 2, the units of various parameters are as follows: the unit of maximal conductance is mS · cm −2 ; the unit of reversal potentials is mV; the unit of τ x is ms; and λ is a constant. From   The symbol L s stands for this pattern [21]. As seen above, MMOs may be produced in the ARPS model. Then we further analyze their generating mechanism and evolutive law by parameter variation in the reduced model. The causes of MMOs need to consider two aspects: the generation of large amplitude of oscillation and the generation of small amplitude oscillations. For large amplitude oscillations, they often require systems with a global return mechanism, and the folded node, the singular Hopf bifurcation and different time scales are likely to make the system to produce small amplitude oscillations [9]. Brøns et al give the existence theorem of MMOs in a three dimensional singular perturbation system [6]. They stress that if the system has a cubic-shaped critical manifold and a singular funnel caused by the folded node, and the global return point to fall again within the singular funnel, the system can produce MMOs. By using the method and numerical calculation, the existence of MMOs of a reduced three dimensional ARPS model will be proved first.
Because the ARPS model is a 4D model, it has to be reduced. Observing the model we can see that, V is a fast variable, where w and x are slow variables. However, c is a slower variable than others between 0.07 − 0.09 µM. In equation (4), K p = 0.00025 1, R = 0.0045 1, so let K p → 0 and R → 0, and then equation (4) becomes dc dt = 0.
So the calcium ion concentration in the ARPS model c is considered as a constant.
Combined with the range of calcium ion concentration, we take c = 0.086 µM.
where g = 100 ms · cm −2 , k v = 100 mV, k t = 1 ms,ḡ i = g i /g,V i = V i /k v , and i is taken as N a, K, K(Ca), Ca, respectively (whereV K(Ca) =V K ).Ī = I/(k v g). Because k t = 1 ms, for simplicity, τ w /k t and τ x /k t still denote by τ w and τ x , but they are dimensionless here. Set then system (1)-(3)can be rewritten as: So a three dimensional singular perturbation system (5)-(7) can be found, and ε is the perturbation parameter. This system is called the reduced Av-Ron-Parnas-Segel model, i.e. the reduced ARPS model for short.
In system (5)- (7), by setting ε → 0 we get the slow subsystem (reduced problem): Using the time scale transformation t 1 = τ · ε, system (5)-(7) is rewritten as: By setting ε → 0, the fast subsystem (layer problem) is Solution of fast subsystem are called fast fibers. Obviously, they are lines which parallel to the (w, x) plane. The critical manifold is composed of equilibrium points for the fast subsystem , and it's denoted by S. That is, We can prove that f (v, w, x) and exists. Therefore, considering continuity of f (v, w, x), S is a folded surface as the graph of the critical manifold S as shown in Fig.2 a). So we have: Proposition 1. The critical manifold of system (5)-(7) S 0 ⊂ S is (locally) a folded cubic-shaped surface.  (5)-(7) at τ x = 2.6. b) Folded curves on the critical manifold.
As shown in Fig.2, the critical manifold is locally a folded cubic-shaped surface. That is, are two attracting branches of S 0 . In addition, are folded curves of S 0 , which are shown in Fig.2 b) explicitly. In Fig.2 b), S 0 denotes the critical manifold and intersects Definition 3.1. Suppose that a point p on folded curves satisfies the transversality condition (normal switching condition) Then, the point p is called the jump point.
Therefore, there is the global return mechanism in system (5)-(7) because of the existence of folded cubic-shaped critical manifold: trajectories starting from the attracting branch S + a reach a folded curve L + along the critical manifold in finite time. Then they leave the attracting branch S + a along fast fibers at a jump point p in the folded curve L + and reach the opposite attracting branch S − a , and then arrive at the other folded curve L − along the critical manifold, and return to S + a along fast fibers. By repeating the above-mentioned process, large amplitude oscillations can be reached.
Folded singularities are those points which do not satisfy the normal switching condition. The existence is an important mechanism to produce small amplitude oscillations. Because slow subsystem is degenerated at folded singularities, the system (8)-(10) should be desingularized. Taking the derivative of equation (8) with respect to t 1 , we can obtain Using transformation dt 1 = −f v · dt 2 , we rewrite system (8)- (10) aṡ x where" · " denotes defferentiation with respect to t 2 . Compared with the slow subsystem(8)-(10), system (14)-(16) has the same phase portrait, but the orientation of trajectories is reversed on the repelling branch. Thus, we can study system (14)-(16) on the phase portrait. From equation (8) we have Then system (14)- (16) can be discussed by projecting this system on the (v, w) plane, or on the (v, x) plane by the using implicit function theorem. Here we takē V K = 0.7102, τ x = 2.6, and other parameters take values as in Table 1 and the first column in Table 2 to compute whether there is any folded singularity in system (14)- (16). The following results are obtained by a direct calculation.   Fig.3 shows that the folded curves of the critical manifold S 0 are projected on the (x, v) plane. Fig.3 b) shows magnification of a) . In Fig.3, the red curves L ± denote two folded curves. The black point p 0 is a folded node in system (8)- (10). The dark green curve SC is the strong canards [30] at the folded node, and the singular funnel bounded by SC and L + (shadow area in Fig.3 b)). In the singular funnel, small amplitude oscillations can arise. The blue curve Γ 1 and the purple curve Γ 2 denote trajectories of system (5)-(7) at τ x = 2.5 and τ x = 2.7, respectively. Finally, we are looking for the singular periodic orbit to connect both of two types of oscillations. According to [22], let P (L ± ) be the projection along the fast fibers of L ± onto the opposite attracting branch. If there is a subset N ± ⊂ P (L ∓ ) such that all trajectories of the slow subsystem with initial conditions in N ± arrive at the fold curve L ± , we can define the associated maps Π ± : N ± ⊂ P (L ∓ ) → L ± . Furthermore, we can also define the return map Π P •Π + •P •Π − : N − → P (L + ), and Π(N − ) ⊂ N − . The existence of a singular periodic orbit can be proved by the contraction mapping principle or Brouwer's fixed point theorem. If the image Π(q) of q ∈ N − falls into the singular funnel by the return map, MMOs can be generated. Otherwise, relaxation oscillation [11] will be produced.
From Fig.3, we can see that, when τ x = 2.5, the trajectory (blue curve) does not fall into the singular funnel by the return map, so the relaxation oscillation is generated in the system. When τ x = 2.7, the trajectory (purple curve) drops into the singular funnel by the return map, so MMOs can be observed in the system. The trajectory (purple curve) can not reach the folded curve in finite time at τ x = 2.9, so the system keeps resting. Moreover, we can obtain the following conclusion by numerical calculations.

Proposition 3.
When τ x ∈ (2.52633, 2.81142), there is a singular periodic orbit Γ = Γ s Γ g for system (5)-(7), where Γ s lies in the singular funnel, a Γ − f , and the ending point of Γ g drops into the singular funnel or on the strong canards, where Γ ± a is a trajectory which connects P (L ± ) and L ∓ , and Γ ± f is a fast fiber which connects L ± and P (L ± ).
According to [6], we can derive the following result immediately.

Theorem 3.2.
Suppose that all assumptions in Propositions 1-3 are satisfied, then for a sufficiently small ε, there exists a stable MMO with signature 1 s in system (5)-(7) for some s > 0.
We focus on the existence of the reduced ARPS model above. Next, the effect of parameters on MMOs in the reduced ARPS model is explored. When we solve system (5)-(7), we find that different types of MMOs can be generated in the reduced ARPS model. Thus, InterSpike Intervals (ISIs) [15] of small amplitude oscillations of the membrane potential and the firing number can be used to analyze the transition of MMOs in system (5)-(7). The period adding bifurcation of ISIs is shown in Fig. 4. When 0 < τ x < 2.52633, relaxation oscillations arise in the system, i.e. the 1 0 MMO pattern. There is no any small amplitude oscillation in the model. When τ x > 2.52633, small amplitude oscillations begin to appear and the 1 1 MMO pattern is generated in the model. With the further increase of τ x , the number of small amplitude oscillations is increasing. The 1 2 pattern, the 1 3 pattern, ... can be observed. When τ x approaches 2.8, the number of small amplitude oscillations increase rapidly. When τ x > 2.81142, oscillations of the system terminate and the system enters the resting state. For the L s MMO patterns, we can define their firing number as [18]: Therefore, we summarize the transition of MMOs in the model by the change of firing number with the parameter τ x . As shown in Fig.5 a), the firing number of MMOs monotonically increases and remains stable in some intervals with the increase of τ x . We also observe that a point whose firing value approximates to 1 appears around τ x = 2.62. By calculating the firing number further we can find when τ x ∈ [2.62054, 2.62066], MMOs with signature 1 2 and 1 3 alternate, and a non-periodic region occurs. One can see that the chaotic region appears when τ x ∈ [2.62054, 2.62066] by calculating the largest Lyapunov exponent (LLE).

3.2.
MMOs in the ARPS model. There is no mechanism of the generation of MMOs in the four dimensional dynamical systems, so we analyze MMOs in the ARPS model by the numerical method. Intracellular calcium ion concentration is changing in the ARPS model, and influences the membrane potential. At the same time, the change of membrane potential influences the intracellular calcium ion concentration, so MMOs in the four dimensional ARPS model are also more complex than the three dimensional system. Not only the 1 s pattern MMO (s > 1), but the L s patten MMO (L > 1, s > 1) can also arise in the ARPS model. MMOs appeared during experiments involving the lobster cardiac ganglion [8]. In our simulations, we also observe these patterns in the ARPS model. Both the time constant and maximum ion conductance are the key parameters to produce the MMOs. In this part, we first investigate the effect of time constant on MMO patterns. We adjust the values of parameters to generate the firing pattern 1 as shown in Fig.1 a). In this case, fixing V K = −71.02 mV (the chaotic region) and changing the values of time constant τ x induces the variation of MMOs.
From Fig.6, one can observe that MMOs can be produced by the basic bursting model. The time constant τ x has great influence on the patterns of the MMOs, and it affects the number of both L and s. At τ x = 2.7 the 1 1 MMOs are obtained. With the decreasing of τ x , the MMOs appear in the forms 2 1 , 4 2 , 5 3 and 11 7 . At τ x = 1.0003 the MMOs disappear and the neuron begins to produce an all spiking firing pattern. The mixed-mode wave forms can be summarized in the following diagram.
As shown in Fig.7, the firing number is horizontal like a piecewise function with τ x , and gives rise to a step on this staircase. It is shown that the MMOs will keep a fixed pattern for a particular duration. Between these steps are gaps, corresponding to the intermediate ranges of τ x where nonmixed-mode behaviors [24] are found.
In Fig.7, there are some gaps between the steps which can be seen when τ x increases. Specifically, the time courses of the membrane potential are calculated carefully between these gaps. For example, the values of τ x is taken between 2.078 and 2.079. At τ x = 2.079, the model produces the 2 1 MMO pattern firing. When the value of τ x is taken as 2.0787908, a bursting pattern which has the 3 1 3 2 form arises. While τ x is decreasing, more bursting patterns appear. At τ x = 2.0787, a  bursting pattern with a 4 2 form appears. At τ x = 2.078, the time courses of the membrane potential are non-periodic as shown by the following diagram. Because the bursts appear aperiodically and rarely appear over long time intervals, and the frequency of appearance of bursts increases, this gap may be an intermittently chaotic region.
Maximum ion conductance plays a key role in producing MMOs in the ARPS model. There are similar results with time constants which modulate their values. For example, when we change the values of g K and g K(Ca) , MMOs in the model also appear. Some cases are given in the following figures. The reversal potential of potassium parameter can also play the key role to transform the number of spikes of small amplitude (sub-threshold) oscillations in MMOs. The number of spikes of sub-threshold oscillations increases when the values of parameter V K decrease. For the interspike, where the spike means sub-threshold oscillations, interval sequences can also induce to the period doubling bifurcation phenomenon in the ARPS model (see Fig. 9).
3.3. The influence of parameters to the ARPS model. In order to understand the behavior of the individual neuron in the lobster giant axon, the property of the neural coding is studied in this subsection. Evyatar Av-Ron et al considered whether the model could account for the experiments of Dalton, concerning changes in external concentrations of K + [1]. So the reversal potential of potassium parameter is considered first. The effect on the membrane potential during changes in the reversal potential of potassium is discussed. Specifically, we investigate the abundant dynamical properties of bursting in the model through the transformation of the ISIs, and it is observed that changing the reversal potential of potassium has great influence on the response of the lobster cardiac ganglion. When the parameters values are taken to generate the first firing pattern (Fig. 1 a)), and changing the values of V K , the bifurcation figure of the ISIs can be gotten (Fig.10) . When the reversal potential of K + decreases, the firing patterns constantly change and the interspike interval increases in the model. After a chaotic bursting region arises, tonic spiking appears in the model. In the above diagram, the bifurcation pattern in neuronal action potential ISIs can be observed. The sequences of the action potential ISIs show a clear period doubling bifurcation phenomenon with the perturbation of the reversal potential of K + . With the decreasing of V K from −70.94 to −71.02 mV, the ISIs sequences of this neuron model appear as the period doubling bifurcation, and then transform into a complex region. Then by taking the parameters values to generate the last firing pattern (Fig.  1h)), another bifurcation pattern in ISIs sequences is obtained. In this case, the bursting sustains presence, and the number of the spikes in one burst increases along with a decrease in the reversal potential of K + .
In Fig.11, it presents the bifurcation pattern without chaos in which neuronal action potential ISIs can also be observed. The sequences of the action potential Figure 11. The period adding bifurcation diagram of V K vs ISIs at the firing pattern 8 as shown in Fig.1 h).
ISIs show a clear period adding bifurcation phenomenon with the perturbation of the reversal potential of the potassium ion. With the decreasing of V K from −60 to −85 mV, the time courses of the membrane potential of the period adding firing sequence vary from period-1 to 7. The numerical results show that the perturbation variables and reversal potential of the potassium ion have great influence on neuronal activities (action potentials). These rich bursting patterns and the switching between different firing patterns tell us that persistent potassium currents are very important in the generation of the rich bursting. In the description of the ARPS model, the change in the conductance of each ion is determined by three variables: w, x, and c, where the change of w affects sodium, potassium and calcium ion current, and the change of x and c affects the calcium and potassium-calcium mixed ion current. The speed with which these variables are altered determines the ion channel current's amplitude with regards to time, and then influences the firing pattern of the action potential, so the perturbation of the ion channel variable is critical to the activities of neurons. Now, we discuss the ion conductance in simple terms. As shown in Fig. 12, from spiking to bursting, three chaotic regions which start at 0.305, 0.322, and 0.357 ms · cm −2 , respectively, arise in the model. Furthermore, the clear period doubling bifurcation phenomenon of the sequences of action potential ISIs is observed with the conductance of the potassium-calcium mixed ion. With the increasing of g K(Ca) from 0.292 to 0.305 mS · cm −2 , the ISIs sequence of the model produces the period doubling bifurcations, and until it turns into a complex region.
4. Discussion. We have specifically explored MMOs and dynamical properties of the ARPS model in this paper through the geometric singular perturbation theory and simulation results, and found a rich dynamics of the model.
There are various firing patterns in the model. In this study we can find eight kinds of time courses in the model, when the electrophysiology parameters vary in a certain region, not only periodic spiking but also various bursting patterns. Then the MMO firing patterns in the ARPS model are studied in this paper. The existence of MMO patterns in the reduced ARPS model was presented. Then the MMOs can be observed by using the numerical simulation in this paper. Particularly, we calculate the non MMO patterns in the model. The non MMO patterns could potentially be a new chaotic phenomenon in the model.
Other changes of parameters can also cause the MMOs in the model, and the process of the MMOs change varies dependent on the values of the τ x . From Fig.13, one can see that the maximum ion conductance also affects the generation and variety of MMOs in the ARPS model. Finally, more intuitive results of how these firing patterns change from one pattern to another one can be indicated by the ISI bifurcation diagrams. The period doubling phenomena and the period-adding phenomena contained in the ISI bifurcation figures tell us about the intrinsic properties of corresponding parameters on the spiking transition modes and firing patterns in the lobster cardiac ganglion. In fact, many other parameters in this model can also affect the firing patterns transition, such as the potassium ion conductance, sodium ion conductance, and external stimulus et al.
When we study the existence of MMOs in the reduced ARPS model, the mechanism of generation of MMOs is given by the application of the folded node theory. But there is a fly in the ointment: the number of small amplitude oscillations s is not estimated. Because the MMOs are complex dynamic behavior, the relevant theories of MMOs are still not complete, and it is difficult to deal with some nonlinear problems. The ratio of eigenvalues at the folded node is µ ≈ 0.43, and the perturbation parameter ε = 0.01, thus µ = O( √ ε). It shows that the mechanism in which small amplitude oscillations arise in the system can be explained by the folded node theory or the singular Hopf bifurcation theory. From Fig. 2 we can see the folded curve is almost parallel to the x axis, thus it can be projected on the (v, x) or (w, x) plane . But because of the complexity of the critical manifold, it is impossible to give the explicit expression of w or v. Although the critical manifold can be written as the explicit function of x, if the critical manifold is projected on the (v, w) plane, a folded focus can only be obtained. So the folded node is calculated by the numerical method.
The ARPS model can better describe the lobster cardiac ganglia physiological activities, but because it is a foundational work, it seems that the model is neglected for a long time. Nevertheless, rich dynamic phenomena unearthed from the model, especially subthreshold membrane potential oscillations are ignored previously. Therefore, it is interesting and makes sense to study more on the physiological activities of the cardiac ganglia lobster in the futrue.