Controlling switching between birhythmic states in a new conductance-based bursting neuronal model

A birhythmic conductance-based neuronal model with fast and slow variables is proposed to generate and control the coexistence of two different attracting modes in amplitudes and frequencies. However, periodic bursting, chaotic spiking and bursting have not been clearly observed there. The control of bistability is investigated in a three-dimensional birhythmic conductance-based neuronal model. We consider slow processes in neuron materialized by an adaptation variable coupled to system in the presence of an external sinusoidal current applied. By using the harmonic balance method, we obtain the frequency–response curve in which membrane potential resonance with his corresponding frequency is controlled by varying a specific parameter. At the resonance frequency, bifurcation and Lyapunov exponent diagrams versus a control parameter are obtained. They reveal a coexistence of two different complex attractors, namely periodic and chaotic spiking, periodic and chaotic bursting. By using the control parameter as the slow variable, the system can switch from bistable to monostable behavior. This is done by destroying subthreshold (small) oscillation of the neuron. The role of adaptation variable in neuron system is then to filter many existing electrical processes and permit to adapt the system by the multiple transitions states in the chosen electrical mode. A fairly good agreement is observed between analytical and numerical results.


Introduction
Multistability is a common dynamical feature of many natural systems [1]. It has been observed in modeling studies [2] such as coexistence of several attractors in the phase space of a system. Bistability is an interesting manifestation of multistability which occurs very frequently and deserves special attention due to their importance in physics and biology. It plays a significant role in living systems in maintaining different modes of oscillations to various biochemical processes in varying environment [3]. It is a fundamental attribute of the dynamic of neurons and neuronal networks [4] because neurons can be viewed as strongly nonlinear dynamical capable of generating complex patterns of activity. The coexistence of different attracting regimes of activity that reflect their intrinsic ionic channel behaviors within a specific parameter range, is not uncommon for a such system. The particular attractor exhibited by it, is determined by initial conditions. Multistability by bifurcations phenomena sets a remarkable framework for potential dynamical plasticity of neurons with implications for dynamical memory, information processing and motor control [5].
The coexistence of neuronal activity regimes has been observed and demonstrated by theoretical and experimental studies [6] and recently by Tagne Nkounga et al. [7]. One of the most studies coexistence states is the coexistence of tonic spiking and silence [8][9][10]. Compared with the Hodgkin Huxley (HH) [11] and the 2D Morris-Lecar (ML) [8] models, the 2D birhythmic conductance-based neuronal model is so simple. In fact, this mathematical model can generate various action voltage potentials [7]. It cannot trigger bursting firings or complex dynamical transition from periodic to chaotic motions [12] that appear generally in biological neuronal systems with fast and slow variables. Because the using of multistability can permit to develop therapeutic strategies of some diseases such as epilepsy, it is so important to identify an intrinsic property to control multistability which occurs in neuron. Many studies have showed a way to control multistability in living systems. Guttman et al. [13] have shown in experiment on the squid giant axon in saline with low calcium concentration that a switching between two regimes can be executed by a pulse of current. Ghosh et al. [14] reported an effective control mechanism of birhythmic behavior in a modified van der Pol system by using a variant of Pyragas technique of time delay control [15]. They have shown that depending upon the time delay, system can induce monorhythmic oscillations out of birhythmicity. Ghosh et al. [14] established that their technique can suppress the effective birhythmic zone.
HR have shown in their 3D neuronal system that the third linear equation coupled with their previous twodimensional model [10], controls the rhythmic states transitions of attractors. This exhibited periodic spiking, periodic and chaotic bursting [16]. In the Jerk systems, the control of multistability was based on the linear augmentation scheme [17]. It is also showed in a leech model that the adding of slow variable to fast subsystem can switch the system from the coexistence states [6]. Although the coexistence of complex dynamics presented in some biophysical models, these models are complex and the demonstration of these phenomena is generally done through numerical simulations. In the present work, we propose and analyze a 3D birhythmic conductance-based neuronal model. The model derived from the 3D Prescott model [18] introduced as an extension of the 2D ML model. We focus on the control of switching between coexisting attractors. The importance of the present work resides on the fact that it possesses multi-limit cycles. More importantly, we analytically define the system parameter's domain where these many attractors coexist. On the contrary of other models where the coexistence between two periodic regimes (birhythmicity), random oscillations (chaos) and coexistence of a stable limit-cycle oscillations are done numerically [19][20][21], our model offers a room for analytical analysis. Note that the understanding of the duality (birhythmicity) of neurons has motivated many researchers in proposing some models [17,[19][20][21]. The faithfulness of these models is to underline that biological neuron and it complexity are crucial in describing the information processing and transmission in neuronal system [5]. Thus, by considering slow processes as an adaptation variable in our model and by applying an external sinusoidal current, we analyze analytically and numerically the coexistence of periodic or/and chaotic regimes and the control of these various activities in the neuron.
The work is organized as follows: In Sect. 2, we use a 3D Prescott model to model a 3D birhythmic conductance-based neuronal model by adding a slower variable than the first two in the 2D model. In Sect. 3, we proceed to a deterministic dynamical analysis of the model in absence of a third variable. Amplitudes, frequencies and the stability of oscillatory states are studied using, respectively, the perturbation method and the Floquet theory. We then study the effect of the third added variable without and with a forcing term on a 3D birhythmic conductance-based neuronal model. These use firstly analytical methods as harmonic decomposition and harmonic balance method, respectively, and secondly numerical methods as the bifurcation and Lyapunov exponent. In Sect. 4, we conclude this work.

The conductance-based neuronal model
Many mathematical models were used to describe neural activity. A two-dimensional Prescott neuron model is an improvement of ML neuron model (describes the oscillations in the barnacle giant muscle fiber) [18]. It consists of a fast activation variable and a slow recovery variable or activation of potassium current [22]. A three-dimensional Prescott model introduces an adaptive ion current due to sodium channel inactivation since it operates on the same timescale as delayed rectified channel activation. This is in order to make the model more realistic (a real description to emission of membrane potential about the long timescale between action potential). Since adaptability is a prominent characteristic of neural dynamics, the dynamical equations are as follows: where System (1) consists of voltage-gated N + a current, delayed rectified K + current, adaptive currents and the leak current. Here v, w, z represent the membrane potential of the neuron, the activation variable of K + ion channels (the recovery variable of slow ionic currents) and the gating variable of the adaptive currents, respectively. The parameters g Na , g K , g a and g L indicate the maximum conductance functions of ions N + a , K + , adaptive currents and leak currents, respectively. E Na , E K , E a and E L are reversal potentials of the different tonic current functions. C means the membrane capacitance and it is considered as unity. φ w and φ z represent the temperature scaling factor for the K + channel opening and ions adaptive currents, respectively. β m , γ m , β w , γ w , β z , γ z are other parameters of model, and I is the external current applied to the model. In the presence of adaptive current, neuron can exhibit periodic and chaotic bursting, coexisting multistable firing and the transition to each other as shown in the 3D ML model [23]. In the present work, the three trigonometric hyperbolic functions present in a system are expanding using the Taylor series. Considering the higher-order nonlinear terms and taking into account some selected terms, we can write eq.(1) in the form: where the parameters a, b, B, α, θ 1 , θ 2 , c, r , s, γ 1 , γ 2 , γ 3 , γ 4 and their range values, are given in appendix.
It is also noted as in [7] that θ 1 , θ 2 , α are the function of the maximum conductance of ions channels. By setting a = 2.0, c = B = γ 3 = 1.0, b = 1/3, γ 4 = 0.0, s = 4.0 and γ 2 = 5.0, we obtain the new neuronal mathematical model below which will be used in all this analysis: The neuronal model described by eq.(4) is obtained by an approximation of the biophysical improved Morris-Lecar model, namely Prescott's neuronal model. Contrary to a simple three-dimensional mathematical model as HR where the function f is taken as a cubic polynomial function, we establish this simple model in which f is a septic polynomial function and g is a quadratic polynomial function. In the simple tridimensional model of HR [16] as in the three-dimensional Prescott model [22], the third variable is shown as a variable of bifurcation control of a multi-stable dynamic. Thus, a brief current given by the variable z called adaptation variable whose phase and amplitude given in an appropriate range, can provide hysteresis and switching of the system between coexisting stable states [24]. It also may suppress bistability and maintain the desired dynamical behaviors in some cases.

Amplitudes, frequencies of oscillatory states and their stability
In optic to find concordance of analytical and numerical results, we proceed to reduce equations of system (4) in the absence of adaptation variable z (α = 0) and with a constant external current, in a nonlinear equation that we can solve through the nonlinear analytical analysis methods known by all. That is why, by eliminating w in system (4), we obtain after setting t = ω 0 τ in the found equation, the following equation similar to what found in [25] describing the multi-limit cycles lineartype equation: where This equation is of the Leanard form and can be solved analytically using the Lindsted's perturbation method [26]. In order to find amplitudes and frequencies of the system, we posed τ = ωt , where ω is a unknown frequency. The periodic solution of (5) can be expressed by: where x n (τ ) (n = 0, 1, 2, ...) are periodic functions of τ of period 2π . Therefore, the frequency ω can be represented by the following expression: where ω n (n = 0, 1, 2, ...) are unknown constants at this step. By substituting Eqs. (6) and (7) in (5) and finding the secular terms in different order of μ, a solution of Eq. (5) can take this following approximation: x(t) =A cos ωt + μ( 2 sin ωt + 3 sin 3ωt where the amplitude A satisfies the following equation: and the frequency ω is given by with Fig. 1 Region of existence to one (unhatched part) and three limit cycles (hatched part) in the plane (θ 1 , θ 2 ) obtained by solving algebraic Eq. (9) All data refer to the case γ 1 = 100.0 The parameter used is: Amplitude Eq. (9) can be solved using an algorithm of Newton-Raphson [24]. The results appearing in Fig. 1 in the (θ 1 , θ 2 ) plane describe the area of one (unhatched) and three amplitudes (hatched), respectively. As shown experimental by Lechner and al. [6], neuronal system (4) can emerge in two different stable attractors separate each to the other by an intermediate unstable attractor. This phenomenon called bistability has been mentioned by many authors [7,27]. The particular attractor (electrical mode) exhibited by the neuron is determined by initial conditions. Therefore the initial conditions of the model represent membrane potential and the states of the ion conductances [27]. Thus the attractor with small amplitude and high frequency characterizes a subthreshold potential, while the attractor with large amplitude and small frequency characterizes an action potential as shown in [28].
A verification of the convergence of results is evaluated in Table 1 where the comparison between analytical and numerical results has been done, as well as the variation of the amplitude versus of θ 1 and θ 2 (such as in Fig. 2) for γ 1 = 100.0. The neuronal activity The coexistence between two regimes of self-sustained oscillations corresponding to the smaller (amplitude A 1 , and frequency ω 1 ) and larger (amplitude A 3 , and frequency ω 3 ) limit cycles. The set of parameters is the same as for the attractors shown in Table 1. The parameters used are: γ 1 = 35.0 and (d) has been given in rest, spiking or bursting [16,29,30] and in different forms of coexistence states [7]. In our system, the flexibility in the behavior of the nerve cell is shown by the coexistence of two different oscillatory modes in amplitude and frequency. Figure 3 shows the coexistence between two regimes of self-sustained oscillations, as numerically shown by Lechner and al [6]. Bistability characterizing a existence to bifurcation by a coexistence of different modes (such as here electrical activities), it gives to neurons a rich forms in his information processing to generating pulse. We have in Fig. 4 the representation of the unstable attractor in dotted line and hatched part obtained using the Floquet theory [7].

Statement of the problem
The problem of bidimensional neuron's models is that it does not give a real description to emission of membrane potential about the long timescale between action potential. It is shown that bursting behavior cannot appear in two variables [24]; then, the slow modulation of spiking during a burst (limitation to the time of emission to nerve impulse by a rest state or subthreshold oscillations) requires additional physical mechanism and dynamic variables. It is the reason for what HR [16] following Izhikevich [23] have added a third variable to stopped spiking modes. That is why in the goal to have a real system, we add a slow processes ( channel affected by slow ion accumulation) to our idealized two-variable model and we obtain three-dimensional birhythmic conductance-based neuronal model eq.(4).

Effects of adaptation variable with a constant external current
Let us recall that (α = 0), the evolution of the dynamic of our neuronal system is now represented by the following three nonlinear differential system: where z(t) is an adaptation variable characterizing the slow flux ions current in the membrane. This parameter can be viewed as a low-pass filter. α is the self-feedback strength and r the cutoff frequency. Other parameters have the same significance as upper. It is clear that if (α = 0) the slow z variable becomes a parameter independent to the fast subsystem v(t) and w(t), we therefore have a two-dimensional system without a third variable. By considering system (4) and taking same considerations as the study of the precedent two dimensional equation, we obtain the following equation: where An appropriate analytical tool to understand the effect of r and α on the dynamic of the system, is the harmonic decomposition method [25]. Let us assume the approximate solution of (13) as: and where A, ω, φ 1 are, respectively, the natural amplitude, frequency and phase of the oscillator. B, ω, φ 2 are, respectively, the amplitude, frequency and phase of the filter. Therefore, considering a phase difference φ = φ 2 − φ 1 and letting φ 1 = 0 for convenience, we can write solution of eq. (13) as: and Substituting of (16) and (17) into (13), not taking into account the higher harmonics terms considered as forcing terms [31], we obtain: and Let us assume as [31] that the amplitude B of the filtered signal, z(t), remains almost the same as the original of the signal x(t), but different phase φ. Using this last consideration, (19) into (18) leads us, respectively, to the following frequency and amplitude equation we observe that for α = 0 or r = 1, ω 2 1 and Eq. (21) is the same to amplitude Eq. (9) obtained in birhythmic oscillator (5) in the absence of feedback. For other case, Eqs. (20) and (21) prove that the behavior of the system is influenced by α and r; therefore, the system has an effective natural frequency and amplitude. The frequency in the harmonic limit reduces to ω = 1, and the three limit cycles (two stable, one unstable) are results of the three roots of the amplitude equation. To get an idea of the amplitude of the limit cycles, we analyze the stability of the system using the energy balance method [31]. In order to analyze our system (5), we rewrite Eq. (13) as the following equation: The solution of Eq. (22) for μ eff = 0 is given by: and The change energy E in one period 0 ≤ t ≤ T is found by considering term (25) as the external forcing term.
Therefore the change energy is given by: So, Eqs. (23) and (24) into Eq. (26) give: The presence of limit cycle occurs for E = 0, that's why the presence of limit cycles and their stability can be determined by solving Eq. (27). Then to do this, we plot h(A 2 ) versus A. A stability of the limit cycle is determined by the slope of the curve h(A 2 ) versus A at the zero crossing points. The negative slope determines the stable limit cycle; thus, we can write as the condition of a stable limit cycle. It is so evident as shown by the Floquet theory that for a birhythmic dynamic, the second amplitude is unstable. Figure 5 shows the curves of h(A 2 ) versus amplitude A for different values of α, with r = 0.0001 (see Fig.5a) and for different values of r, with α = 7.0 (see Fig.5b). One finds that for fixed range values of r , the graph is displaced upward with the increases of r . However, an opposite effect is observed when increasing the coupling parameter α. It is important to note that the slow variable called here feedback variable plays an important role on the existence and stability of the three amplitudes in our new model. The map in Fig.  6a resumes the interaction between r and α and shows the couple (α ,r) in which we can have one solution (monorhythmic mode) or three solutions (birhythmic Fig. 12 Times series of a coexistence to period-11 subthreshold oscillation and spiking mode α = 0.28 (a), coexistence of chaotic subthreshold oscillation and period-2 bursting α = 1.03 (b), coexistence to chaotic subthreshold oscillation and period-3 bursting α = 0.56 (c), coexistence to chaotic subthreshold oscillation and period-4 bursting α = 0.44 (d), coexistence to chaotic subthreshold oscillation and chaotic bursting α = 0.46 (e), existence to a monostable chaotic bursting α = 1.71 (f); for the two initial conditions (2,2,5) and (5,5,5). The parameters used are:  Fig. 6a show, respectively, the area of existence of three solutions (not effect of filtered feedback then permeability to potassium K + or sodium N a + ions) and one solution (effect of filtered then strong permeability only to sodium N a + ions). It is also shown in Fig. 6a-c that when α > α c , it is possible to have effect of filtered feedback. But when r > r c , there is no effect of filtered feedback for all values of α. Therefore, effect of z(t) on the dynamic of the system strongly depends on value of r . Large values of r have no filtering effect, that is the neuron can stay to one subthreshold oscillation or to emission of action potential depending on a state of ion conductance. For small values of r , there is a strong possibility of filtering depending to the self-strength feedback α. That is why H.R [16], Rinzel and Ermentrout [24], Andrey and Cymbalyak [27] said in their dynamic's studies that only r taken in range 0 < r 1 can permit to make z(t) a slow variable than the fast variable v and w to have an effect on the dynamic.
In order to verify our analytical results, we firstly show in Fig. 6d the bifurcation diagram obtained numerically which confirms the coexistence of only regular states and the switching at the critical value of α. Secondly, we show in Fig. 7, the evolution of membrane potential for two couple of (α,r) taken in map of Fig. 6 with initial condition (v 0 ,w 0 ,z 0 ). When (α, r ) = (1.0, 0.0001) and with the two following Fig. 13 Times series of a coexistence to period-11 subthreshold oscillation and spiking mode α = 0.28 (a), coexistence of chaotic subthreshold oscillation and period-2 bursting α = 1.03 (b), coexistence to chaotic subthreshold oscillation and period-3 bursting α = 0.56 (c), coexistence to chaotic subthreshold oscillation and period-4 bursting α = 0.44 (d), coexistence to chaotic subthreshold oscillation and chaotic bursting α = 0.46 (e), existence to a monostable chaotic bursting α = 1.71 (f); for the two initial conditions (2,2,5) and (5,5,5). The parameters used are:  Fig. 3 consequence, no effect of filtered feedback. When (α, r ) = (7.0, 0.0001), the oscillation with small amplitude A 1 called (subthreshold potential) is suppressed and the dynamic of system is only in the attractor of large amplitude A 3 (called nerve impulse). This is only in spiking mode for any initial condition as shown by Andrey and Cymbalyak [27] and in a Jerk systems [17].

Simultaneous effect of filtered feedback and external variable current
In the goal to approach the near of possible a threedimensional autonomous Morris-Lecar model [23] about characteristics as many complex forms of coexisting and control of multistable states, our system Eq. (12) needs an additional modulation. We consider here that the intrinsic nonlinear terms avoided in our simple model (3) will be taken into account here as a forcing terms which modulate the constant external current With a sinusoidal force, the amplitude of the forced harmonic oscillatory states can be found using the harmonic balance method [31]. Assuming that the fundamental component of the solutions has the period of the external excitation, we express the solution x as: By substituting Eq. (29) into (28), we obtain a particular solution of a variable z as: where 1 = r r 2 + 2 , 2 = r 2 + 2 . By substituting Eq. (30) into system (28), it becomes: where Inserting Eq. (29) into (31) and equating the coefficient of the cosine and sine terms separately, the amplitude of the oscillatory states satisfies the following algebraic equation: where In the absence of adaptation variable (α = 0), we plot in Fig. 8 the analytical and numerical results of the frequency responses curves versus the frequency . This curve at = 1 reveals the resonant properties to the neuronal system due to the presence of subthreshold oscillations. It is appeared that although a forced terms considered, the model can exhibit one or three limit cycle as the autonomous system. In the presence of an adaptation variable (α = 0), we also present the response curve of the system versus frequency for different values of the self-feedback strength α, Fig. 9a, when E is fixed. It appears that the parameter α controls the resonance frequency and amplitude's peak of the system. In fact, when α increases, the resonance's peak decreases until damped completely, while his corresponding frequency increases. Figure 9b gives the response curve versus the self-feedback strength α for some values of frequency of external stimulation . It appears there that the system has a hysterical behavior when the external frequency is not taken in resonance ( = 1). Figure 9c gives the response curve versus α for some values of intensity of external stimulation E. The phenomenon of hysteresis and the switching to the system by suppressing to one stable attractor when α reaches a critical value α c (for example, α c ∼ 1.5 for = 1 and E = 2.25), is observed. Figure 9d shows in the map of r versus α, a coexistence of two different electrical modes (bistability in hatched part) or a existence of only one electrical mode (monostability in nonhatched part). In order to confirm the behavior obtained by the analytical study, a numerical study using system (28) is analyzed by the means of bifurcation and his corresponding Lyapunov exponents diagrams represented in Fig. 10a, b. These figures show clearly the coexistence of two different forms of attractors (periodic or chaotic) got, respectively, with the initial conditions (2, 2, 5) (red) and (5, 5, 5) (black). When α varies, the neuron presents a coexistence with various complex and captivating behavior such as period, perioddoubling cascades and chaos which prove a variability to frequency. At the critical value of α c ≈ 1.6 the system turns from the bistability to monostability for any initial condition, which proves the switching predicts in analytical analysis.
These results were founded by Cymbalyuk and Shilnikov [27] and in the three-dimensional autonomous ML model [12], but coexisting attractors were in the same amplitude. They were also founded in a simple and chaotic Jerk system [17], in the Chua's system [19] and in a Yu-Wang four-wing chaotic system [20]. By taking some value of the control parameter α, we show in Fig. 11 with some phase portrait, the dynamic of the neuron as: a coexistence to period-11 subthreshold oscillations and tonic spiking, Fig. 11a; a coexistence to tonic chaotic subthreshold oscillations and periodic bursting with two, three, four spike per burst, Fig. 11bd; a coexistence to tonic chaotic subthreshold oscillations and chaotic bursting, Fig. 11e. Figure 11f shows a existence to only one chaotic attractor specially with large amplitude. This large amplitude is therefore desirable because it can become more chaotic. The corresponding times series are done in Figs. 12 and 13. Figure 14 shows other profiles obtained in the others phase plane for some values of a control parameter α and initial conditions. Variable z(t) characterizing slow flux ions in membrane, plays an important role in neuron because it can be considered as a filtered. It permits to the neuronal system through his bifurcation parameter α, to moved in attractor who can change significatively and where his efficiency is not limited in his process of transmission of information. This control of the neuronal dynamic may throw a new light on the diagnosis and treatment of neurological diseases.

Conclusion
In summary, we have proposed firstly in this paper, a three-dimensional conductance-based neuronal model inspired by an improvement of M-L system [12]. We have shown that his autonomous dynamic without adaptation variable is the same as a modified van der Pol system, and his dynamic depends strongly to nonlinear coefficients θ 1 and θ 2 function of conductance parameters. By using the Lindsted's perturbation method [26], with an appropriate parameters θ 1 and θ 2 , we have found a coexistence of two different attracting modes in frequency and amplitude (bistability) function to initial condition given in Figs. 1, 2. This behavior characterizes the neuronal system as a system possessing rich forms of information process. This is principally here in subthreshold oscillation (small amplitude A 1 ) or in emission of nerve impulse (large amplitude A 3 ) separated by an unstable attractor A 2 called threshold Fig.  3.
Secondly, in order to make the system more realistic, we inspired to the concept of slow-fast variable in coupling our first bidimensional mathematical model with an appropriate variable called adaptation variable. This is with the coupling parameter α due to ions accumulation in the membrane cellular. We have then quantified the effect of this third variable with a forcing term on our autonomous dynamic by the harmonic balance method [31]. We have shown in Figs. 8, 9 that our system admits a membrane potential resonance with a resonance frequency which can be controlled by a control parameter α. This existence of membrane poten-tial resonance is a feature of a neuron-type resonator. At the resonance, the bifurcation and Lyapunov exponent diagrams versus α given at Fig. 10, have shown a coexistence of two complex attractors, namely periodic, period doubling and chaotic.
These period-doubling cascades are an evidence to the presence of biological behaviors as tonic spiking, periodic bursting, chaotic spiking and bursting, respectively. They are given in Figs. 11, 12, 13. For α ≥ 1.6, neuron judges necessary to switch his activity from a state of two coexisting attractors to monostability also a attractor with large amplitude (action potential), given in Figs. 11f, 12f, 13f. This is by a suppressing of the attractor with small amplitude and high frequency ( subthreshold oscillation ); or a maintains his bistable activity for α < 1.6 when need is it. This last result is so important to neuron because it is permit to be efficient in his function to transmit information by a filtering of many process within before transmission to nervous message.

Data availability
The data from simulations that support the findings of this study are available on request from the corresponding author RY.

Conflict of interest
The authors declare that they have no conflict of interest.