Dynamic mechanism of multiple bursting patterns in a whole-cell multiscale model with calcium oscillations

The dynamic mechanism of a whole-cell model containing electrical signalling and two-compartment Ca signalling in gonadotrophs is investigated. The transition from spiking to bursting by Hopf bifurcation of the fast subsystem about the slow variable is detected via the suitable parameters. When the timescale of K gating variable is changed, the relaxation oscillation with locally small fluctuation, chaotic bursting and mixed-mode bursting (MMB) are revealed through chaos. In addition, the bifurcation of with regard to is analysed, showing periodic solutions, torus, period doubling solutions and chaos. Finally, hyperpolarizations and torus canard-like behaviours of the full system under a set of specific parameters are elucidated.


Introduction
Calcium signal is essential in cells, participating in and controlling many important cellular biological processes and mechanisms, such as the contraction of muscle cells, the release of neurotransmitters from neurons and astrocytes, and the activation of eggs, the healing and metabolism of liver and pancreatic trauma [1][2][3], as well as cell maturation, differentiation and death [4][5][6]. Since the Ca 2+ oscillation was first observed, the discussion of its molecular mechanism from the aspects of experiment and model has never stopped [7][8][9].
Intracellular interactions between a variety of elements regulate Ca 2+ signals: for example, the entry and discharge through plasma membrane (PM), the release and re-extraction from intracellular stores, and the buffering of endogenous Ca 2+ buffers. In detail, the link of hormones or other signalling molecules (glutamate, ATP, etc.) and the G protein-coupled receptor (GPCRs) on PM activates intracellular phospholipase C (PLC), hydrolyses phospholipid, especially phosphatidylinositol-4, 5-diphosphate (PIP2), located in the inner layer of PM. Hydrolysis of PIP2 produces two products: diacylglycerol (DAG) and inositol 1, 4, 5-triphosphate (IP 3 ). IP 3 and Ca 2+ are combined to IP 3 receptor (IP 3 R ), and IP 3 R in turn regulates the release of Ca 2+ from endoplasmic reticulum (ER). A Ca 2+ leakage (leak), from ER to the cytoplasm, may exist. At the same time, Ca 2+ in the cytoplasm can be pumped into ER through the Ca 2+ -ATP enzyme (SERCA) on the sarcoplasmic reticulum. Extracellular Ca 2+ ions enter the cells through the source-operated Ca 2+ entrance, and intracellular Ca 2+ is discharged into the extracellular matrix through the PM ATP enzyme (PMCA). Parts of the principle of intracellular Ca 2+ signal transduction pathway are shown in Figure 1 [10].
Based on the biological principles of Ca 2+ signals, the concept of Ca 2+ 'toolbox' [11] can be used to understand Ca 2+ oscillation in various cells. The common opinion is that the combination of IP 3 and Ca 2+ opening IP 3 R channels is the most important toolbox component, which is modulated to get kinds of Ca 2+ signals with distinguished frequencies and patterns, presented as various nonlinear terms in the mathematical models [12][13][14]. As a result, the dynamic system theory is a powerful tool to discuss the behaviour of Ca 2+ signals, especially the bifurcation theory can give the mechanism of different Ca 2+ signal patterns [12,15,16]. Meanwhile, the change of intracellular Ca 2+ inevitably causes the difference of membrane potential, forming a variety of spiking modes of neurons, in turn, the change of membrane potential also causes distinguished Ca 2+ oscillation patterns. The motivation gives multiple mathematical models [17,18] describing the interaction between the membrane potential and Ca 2+ oscillations. Because of the multiple timescales in the two parts, such models generally correspond to the slow-fast systems, with which the bifurcation analysis shows the mechanism to generate complex spiking patterns and Ca 2+ oscillations. As typical examples, [19] has shown a complete topological classification of bursting about codimension-1 bifurcation, including 120 different topological types. Additionally, it is obtained that the neuron has multiple kinds of bursting: 'fold/fold', 'Hopf/Hopf', 'Hopf/homoclinic' and 'fold cycle/homoclinic' burstings by bifurcation analysis [18,20,21].
Li et al. [10,22] proposed a whole model including PM and ER parts in gonadotrophs in order to show the bursting caused by Ca 2+ oscillations on the basis of the experimental description. The outstanding feature of the model is to consider Ca 2+ signals in both the cytosol and ER, which together with the membrane potential give a system with a fast and two slow time scales, and the slow scales according to Ca 2+ are not at the same time level. Such model extremely tends to declare chaotic bursting, discussed in the similar system in [17]. Our purpose is to investigate the model with two compartments containing only K + , Ca 2+ and calcium dependent potassium currents in order to find the interesting spiking, bursting, mixed-mode bursting (MMB), chaotic bursting, behaviour like torus canard [23,24] and so on combined with bifurcation theory.
The paper is organized as follows. Section 2 introduces the form of the model to be investigated. In Section 3, the spiking patterns according to [Ca 2+ ] i and time scale are explored. In Section 4, the spiking transition and MMB depending on the time scale are discussed. The dynamics of [Ca 2+ ] i is studied in Section 5, which shows periodic solutions, torus and chaos. In Sections 6-7, hyperpolarizations related to [Ca 2+ ] i and torus canardlike pattern are declared. The last part is the conclusion in Section 8.

System model
The paper is devoted to a whole-cell model containing electrical signalling and twocompartment Ca 2+ signalling in gonadotrophs of the form: where V is the membrane potential, ω is the gating variable for the voltage-gated potassium. The functions I Ca , I K and I KCa represent Ca 2+ current, K + current and Ca 2+ dependent K + current, respectively. The flux balance equations (see Figure 1) for ER and PM are shown as follows: where J in PM = −αI Ca , [Ca 2+ ] i is the free cytosolic Ca 2+ between PM and ER. [Ca 2+ ] T is the 'total free Ca 2+ ' of the cell. J PMCA and J SERCA are pump fluxes, J LEAK is the unregulated leak conductance, and J RyR is the RyR channel flux. h is the fraction of channels not inactivated by Ca 2+ . Compared with the model in [22], the system (1)-(2) ignores the external and leakage currents, and only considers the influence of K + and Ca 2+ ions on the membrane potential. Without the evolution according to the space and emphasizing the different time scales, such a multi-scale model will produce complex dynamics. From a biological point of view, carrying signals in the form of Ca 2+ oscillations can not only achieve intracellular signal transmission but also avoid the harm of high concentrations of Ca 2+ to cells. Therefore, intracellular Ca 2+ oscillations and the interaction mechanism with membrane potential is an invaluable subject.

Transition from spiking to bursting
We begin with the relationship between V and [Ca 2+ ] i . Figure 2 shows PM and ER oscillations where [IP 3 ] = 0.9μM, φ = 12/s and other parameters can be found in Table A1.   instead, the decreasing of [Ca 2+ ] i gives the phenomenon of repetitive spiking. The results are consistent with the physiological mechanism.
The timescale curve of ω about the membrane potential V in Figure 2(d) tells us that ω is not a slow variable. Thus neglecting other system slow variables, Figure 2(e) is the projection of V and [Ca 2+ ] i phase plot (the green curves) on the bifurcation diagram, indicating the rapid and repetitive spiking of the PM oscillator during the low concentration oscillation of Ca 2+ . HB represents the supercritical Hopf bifurcation point, corresponding to [Ca 2+ ] i = 0.3139. With the increase of [Ca 2+ ] i , the equilibrium of the fast subsystem changes from the initial unstable node (the red curve) through the HB point to the stable (the black curve) and a stable limit cycle (SLC) appears on the left side of HB, forming a limit cycle of the whole system, as shown in Figure 2f, whose maximum and minimum amplitudes are denoted by the blue curve.

Relaxation oscillations caused by timescale changes
Next, we reduce the value of φ to 0.01 with the same other parameters in Figure 2. As a result, we get a system with one of fast variables changed into a slow one, as shown in Figure 3 (a). Figure 3(b) shows that when [Ca 2+ ] i increases slowly, the PM goes through depolarization, overshoot and the membrane potential V decreases slightly. With the sudden increase of [Ca 2+ ] i , V begins to oscillate and continues to fall down, until [Ca 2+ ] i reaches the maximum. After that, PM completes the repolarization and the hyperpolarization instantaneously. Finally, the oscillation decrease of [Ca 2+ ] i pushes V into a slight hyperpolarization vibration. Therefore, a relaxation oscillation according to the shape of 'a weak tip' of [Ca 2+ ] i becomes visible. The simulation in 3-d space (Figure 3c) shows that V has two parts of oscillation above and below the plane V = 0 with respect to [Ca 2+ ] i , corresponding to parts a and b in Figure 3d, while the two mutations across V = 0 correspond to c and d, respectively. When the time scale (Figure 3e) is reduced by 1/10, the relaxation oscillation phenomenon also appears, but 'the strong tip shape' of [Ca 2+ ] i , as shown in Figure 3f, makes the local oscillation of V stronger, enlarged in Figure 3 (g), and the concentrated oscillation of the two parts in the 3-d space is more obvious (Figure 3h). Choosing the plane ω = 0.4 as the Poincaré section, it can be seen from Figure 3(i) that the points on the section are stacked into two parts, where the small figure in the middle is the magnification of the upper part. There is a reason to consider that chaos occurs. In order to verify this conclusion, we calculate the maximum Lyapunov exponent of the system to find that its value is 4.28534.

Chaotic and mixed-mode burstings according to φ
This section focuses on different forms of cell spiking as a result of variable φ. Figure 4 gives several solutions with regard to the membrane potential V. Using the parameters in Table A1 with φ = 0.7, a bursting solution in Figure 4(a) consisting of one spiking and a subthreshold oscillation is generated, further, when φ = 0.8, the number of spikes per burst has changed, ranging from 1 to 4, as shown in Figure 4 (b). When the value of φ is increased to 0.9, the stable bursting (Figure 4c) with 3 spikes is obtained, and then Figure 4(d) indicates that φ = 1 leads to the increase and instability of spiking number. We also note that when the value of φ continues to increase, there is a stable periodic spiking solution. The process of spiking number increase is extremely complex, and one of attractive reasons for this transition is chaotic dynamics. Poincaré regression map can be used to analyse chaotic behaviour. Starting from the bursting orbit of φ = 1, we draw the (V n , V n+3 ) plane ( Figure 4e) and observe its intersection with the diagonal. It is seen that the points of period 3 exist, which indicates that chaos occurs.
Re-examining the situation of φ = 1, after some periodic spikes, the cell entered the bursting stage, including some small-amplitude slow oscillations (SAO) and a set of repetitive large-amplitude fast oscillations (LAO). As shown in Figure 5 (a), a bursting containing 9 spikes occurs after an SAO, and across another SAO, a 7-spiking bursting is found.  Although the amplitude of the subsequent SAO has changed, the pattern of bursting with 7 spikes lasted 5 times, replaced by 5 continuous spikes, that remained for 13 times. Keep the time moving so that we get a new bursting with 7 spikes, followed by 5 LAOs between 2 SAOs, repeating 12 times. Finally, the time series ends with a 7-spike bursting when t is almost 1500. Thus it can be seen that the system has MMB with different number of LAOs.
In general, φ has an important influence on the spiking transition. Figure 6 shows the bifurcations of the ISI series with respect to φ. It can be seen that in addition to the period-1 behaviour of the membrane potential, there also be multi-period behaviour.

Dynamics of intracellular Ca 2+
In addition to complex bursting, intracellular Ca 2+ oscillation can also cause oscillatory hyperpolarization. As early as 1984, Ince et al. [25] declared that L cells (a mouse fibroblast cell line) and macrophages showed relatively low membrane potentials and slow oscillatory hyperpolarization. By injecting Ca 2+ into P388D1 macrophages, Ca 2+ -dependent oscillatory hyperpolarization of L cells and macrophages was reported, particularly, no oscillatory and single hyperpolarization reaction in Ca 2+ -free medium, so it is deduced that Ca 2+ was the activator of intracellular hyperpolarization.
The last several sections of our paper mainly discuss the dynamic changes of [Ca 2+ ] i and the induced hyperpolarization reaction based on the simulation and the parameters in Table A2. The ([IP 3 ], A) bifurcation diagram is shown in Figure 7 (a), including bifurcation curves: the blue HB is for Hopf bifurcation, the green NS represents Neimark-Sacker bifurcation and the orange PD is period doubling bifurcation, where A is a variable to control the relative time scale between [Ca 2+ ] i and h. It is obtained that the HB and NS bifurcation curves almost overlap for small A, and as A ≤ 2 PD bifurcation exists with the curve like a circle in ([IP 3 ], A) plane.     the stable periodic orbit. When [IP 3 ] = 0.157, the torus vanishes and more complex trajectories in Figure 8(c) arise standing for MMOs. Finally, for the PD bifurcation, we give the period-2 trajectory shown in Figure 8 (d), however, the period-3 orbit indicating the chaotic behaviour is not found, so we discuss the sensitive dependence of [Ca 2+ ] i on the initial value, where the black and red orbits in Figure 7(e) demonstrate different initial value solutions with small disturbances. The distribution of the points on the Poincaré cross section h = 0.6 is also given in Figure 7 (f), forming a scatter plot of dispersion accumulation. These herald chaos.  whereas distinguished from the previous spiking or busting of PM oscillator, under the parameters in Table A2, the significantly low Ca 2+ concentration is true and the oscillatory hyperpolarizations are evoked by periodic oscillation of [Ca 2+ ] i , that is, the membrane potential V is still negative for any time, the time series shown by the black curve in Figure 9 Figure 9(e) by the red curve, which leads to the peak-shaped hyperpolarizations (the black curve). The amplification of peak stage is shown in Figure 9(f) including the other two transient oscillations with large and small amplitudes, seen in Figure 9(g) combined with the pulse oscillation.

Torus canard-like behaviour
Canard refers to the solution of a singular perturbed system which contains both the attractive and repulsive parts. The one-dimensional isolated canard is also known as limit cycle canard [23], appearing in systems with Hopf bifurcation. In a singular perturbed system, canard passes through the fold point (or near the fold point) along the attraction part of the slow manifold, and then continues to form the periodic solution along the repulsive part of the slow manifold. Since the trajectories corresponding to these periodic solutions are similar in shape to canard in phase space in two-dimensional systems, these solutions are called canard solutions. According to the shape of periodic trajectory in phase space, it can be further divided into two cases: canard with head and without head. Similarly, a torus canard implies a family of attracting and repelling limit cycles of the fast system meets in a saddle-node, which is a quasiperiodic orbit from a torus bifurcation in the full system.   Figure 10(b), which tends to be like the headless torus canard. Continue to increase [IP 3 ] to 0.157 and the peak-shaped solutions in Figure 10(c) appears analogous to a torus canard with head, whose transformation trend are similar to that of [Ca 2+ ] T , in other words, the peaks hold as a trough of [Ca 2+ ] T . At this time, the period of [Ca 2+ ] T is slightly long, gives some sparse oscillation circles near H point. Further Figure 10(d) is  Figure 11(a) reveals the complete bifurcation diagram with stable (the black curves) and unstable (the red curve) equilibria, two Hopf bifurcation points ([Ca 2+ ] T = 30.78, 67.17 ), stable (the green curves) and unstable (the blue curves) periodic solutions. In detail, the upper H point bifurcates a periodic solution transforming from unstable one with large amplitude to stable one with small amplitude by the fold bifurcation of limit cycles, which meets the other limit cycle coming from the lower H point so that we can get a green closed curve in the middle of Figure 11 (a). We choose two different initial values of the full system close to the upper and lower H points respectively and find that the orbit with the initial value near the upper point fills about a half of the green circle, however, the second orbit is full of the circle, as shown in Figure 11 (b), that coincides with the time series in Figure 11(c), that is to say wherever the initial point is located at, the orbits distribute in the upper half of the circle.

Conclusion
The article focuses on a whole-cell model consisting of both membrane potential and Ca 2+ oscillations. We show that the system has spiking and bursting behaviour as well as a transition from spiking to bursting. By adjusting the voltage gating of K + , the relaxation oscillation and MMB can be exhibited. Changing the parameters, the variable [Ca 2+ ] i presents the complex dynamics, generating the hyperpolarizations and torus canard-like pattern.
Specially, we first find that in the case of Table A1, whether or not IP 3 exists, the system always spikes, however, the status spiking to bursting happens when its concentration is not zero as Hopf bifurcation of V about [Ca 2+ ] i . The move of parameter φ results in the variety of solutions, i.e. bursting (φ = 12/s), relaxation oscillation (φ = 0.01, 0.001/s), chaotic bursting (φ = 0.7 to 1 with the step 0.1) and MMB (φ = 1).
Parameters in Table A2 give the conclusion that Hopf, NS and PD bifurcations can be revealed according to [IP 3 ], with a stable torus, period doubling solutions and chaos derived. In particular, taking specific values of [IP 3 ] near the bifurcation points, it is shown that hyperpolarization oscillations and torus canard-like behaviour holds.
Mathematical models make it possible to study biological systems theoretically. We can use a relatively simple framework to guide experiments or verify hypotheses, and are not limited by experimental feasibility, cost, ethics and other factors. In the biological experiments of gonadotrophs, the dynamic behaviour of our model including bursting, relaxation oscillation, chaotic bursting and MMB has been found. Furthermore, for cardiomyocytes, it was found that Ca 2+ oscillations would prolong the cells long action potential, leading to arrhythmia, and the dynamic mechanism of this phenomenon is canard [26]. Certainly, due to the complexity of the real biological system, the mathematical model cannot fully capture its characteristics or replace the biological experiment. The two should complement and confirm each other.
Parameters in the paper are as follows: