Next Article in Journal
Approximation of Generalized Ovals and Lemniscates towards Geometric Modeling
Next Article in Special Issue
Controlling Chaos in Van Der Pol Dynamics Using Signal-Encoded Deep Learning
Previous Article in Journal
A Mathematical Model for Early HBV and -HDV Kinetics during Anti-HDV Treatment
Previous Article in Special Issue
Quasi-Periodic Oscillations of Roll System in Corrugated Rolling Mill in Resonance
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computational Analysis and Bifurcation of Regular and Chaotic Ca2+ Oscillations

School of Mathematics and Physics, Guangxi University for Nationalities, Nanning 530006, China
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(24), 3324; https://doi.org/10.3390/math9243324
Submission received: 18 November 2021 / Revised: 14 December 2021 / Accepted: 17 December 2021 / Published: 20 December 2021
(This article belongs to the Special Issue Nonlinear Dynamics and Chaos Theory)

Abstract

:
This study investigated the stability and bifurcation of a nonlinear system model developed by Marhl et al. based on the total Ca2+ concentration among three different Ca2+ stores. In this study, qualitative theories of center manifold and bifurcation were used to analyze the stability of equilibria. The bifurcation parameter drove the system to undergo two supercritical bifurcations. It was hypothesized that the appearance and disappearance of Ca2+ oscillations are driven by them. At the same time, saddle-node bifurcation and torus bifurcation were also found in the process of exploring bifurcation. Finally, numerical simulation was carried out to determine the validity of the proposed approach by drawing bifurcation diagrams, time series, phase portraits, etc.

1. Introduction

Ca2+ is one of the vital ions for information processing in humans. It generally acts as a biological messenger in different cell types and is an indispensable ion for various physiological activities in the human body [1,2,3,4]. Ca2+ plays a particularly important role in the regulation and control of body functions. It is involved in the mechanism of heartbeat and thrombin formation, but also transmits neural signaling, modifies the enzyme activity, and regulates the mechanisms underlying the maturation and fertilization of germ cells. A phenomenon observed experimentally in most cell types is that intracellular Ca2+ concentration remains stable without stimulation. Some Ca2+ enter the cytosol from the Ca2+ store when the cell is stimulated, leading to an increment in Ca2+ concentration that triggers a series of physiological activities [5,6,7,8,9,10,11,12,13]. The cell adjusts itself and closes the Ca2+ channels when the stimulation increases to a certain extent, and then Ca2+ is taken up by Ca2+ stores. After a while, the intracellular Ca2+ concentration gains stability, and the cell returns to a resting state. The source of intracellular Ca2+ is mainly through the influx of extracellular Ca2+ and the release of Ca2+ stored in the endoplasmic reticulum or sarcoplasmic reticulum. The latter is mainly based on the mechanism of Ca2+ release from intracellular Ca2+ stores triggered by Ca2+ signaling, namely Ca2+-induced Ca2+ release [14,15,16].
Over the past few years, many investigations have been carried out on evoked Ca2+ dynamics. Oscillations of free Ca2+ concentrations are highly significant and ubiquitous control mechanisms in a large number of cells. Experimental evidence has demonstrated that Ca2+ oscillations have different characteristics such as regularity and chaos. Chaos refers to dynamic systems that are sensitive to initial values, resulting in unpredictable random motion. Chaotic motion has been associated with making long-term weather predictions and the discovery of Lorentz attractors. In physiology, various arrhythmias, atrioventricular block, and ventricular fibrillation may also be associated with chaos. A well-known manifestation of chaos is the butterfly effect. This refers to a small change that brings about widely varying consequences for the future as it continues to pass and has long been applied to the weather and the stock market. In psychology, for example, when a person receives a small psychological stimulus in childhood, as they grow older, this stimulus can have a large impact on them in adulthood. In addition, the approach to chaos mainly consists of period doubling bifurcations, quasi-periodic transitions, and intermittent chaos. The period doubling bifurcation can be seen in the bifurcation diagram.
After the discovery of oscillations by Wood et al. [17,18], Ca2+ oscillations have been extensively studied. Ca2+ oscillations can be classified into two types: the first type is simple Ca2+ oscillations, and the second is bursting. Simple Ca2+ oscillations have been explored extensively. In the experiments of Borghans et al. (1997), some more complex forms of Ca2+ oscillations have been found such as bursts and chaos. When the system enters chaos, the Ca2+ oscillations show an extremely sensitive dependence on the initial value. Borghans et al. and Houart et al. further proposed several possible mechanisms to explain the complex intracellular Ca2+ oscillations and performed mathematical analysis [19,20]. Chay also proposed a model for Ca2+ oscillations in excitable cells [21], and Shen and Larter (1995) studied complex Ca2+ oscillations and chaotic phenomena in non-excitable cells [22,23]. In these models, it is almost certain that the burst results from changes in IP3 yields because IP3 can stimulate Ca2+ channels and release a large number of Ca2+ into the cytosol or take up Ca2+ from the cytosol, thereby participating in Ca2+ regulation. This is one of the most straightforward explanations for the bursting of intracellular Ca2+. One mechanism proposed by Borghans et al. to explain complex Ca2+ oscillations is that bursting oscillations are associated with Ca2+-releasing channel activities in the endoplasmic reticulum. A further mechanism investigated in the same study is that there is not just one intracellular Ca2+ store. Both IP3-sensitive and non-sensitive Ca2+ stores have been taken into account. In all models, the endoplasmic reticulum is the primary Ca2+ store for intracellular Ca2+.
Marhl et al. suggested the significance of free Ca2+ concentrations in mitochondria for Ca2+ oscillations and proposed a likely explanation for Ca2+ oscillations in cells [24]. In this study, the model focuses on Ca2+ exchange in the middle of the cytosol, endoplasmic reticulum, and mitochondria. The Ca2+ considered in the model includes cytosol, endoplasmic reticulum, and mitochondrial Ca2+ [25,26] and binds Ca2+ in the cytosol. These models well simulate enriched Ca2+ oscillations in cells. Nevertheless, in past investigations, the changes induced by different values of Catot have not been discussed in detail. Moreover, most of the literature has focused on experiments, rarely analyzing the principles of Ca2+ oscillations theoretically. Therefore, this study takes Marhl et al.’s Ca2+ oscillation model as the research object, selects Catot as a bifurcation parameter, analyzes the features of equilibrium by the qualitative theory of differential equation, and discusses the bifurcation of equilibrium by the bifurcation and center manifold theory [27,28,29,30]. The qualitative theory of differential equations allows for analysis of the characteristics of the equilibrium point of the system including stability, coordinates, and types. Specifically, the analysis was conducted using eigenvalues and the Hurwitz criterion. In addition, when investigating the bifurcation of higher-order nonlinear dynamic systems, the center manifold theory can transform the study of various behaviors of n-dimensional dynamic systems approaching equilibrium into the study of equations on m-dimensional (m < n) center manifolds, simplifying the system. Finally, the model was simulated numerically using AUTO and MATLAB in this study. The advantage is that it confirms our conclusions [31,32].

2. Description of the Model

The model established by Marhl et al. investigated the oscillation of free Ca2+ concentrations in the cytosol, endoplasmic reticulum, and mitochondria. A mathematical model was established to simulate the oscillation of intracellular Ca2+ concentration, and the specific system is as follows [24]:
d C a c y t d t = J c h + J l e a k + J o u t J p u m p J i n + k C a P r k + C a c y t P r d C a e r d t = β e r ρ e r ( J p u m p J c h J l e a k ) d C a m d t = β m ρ m ( J i n J o u t ) ,
where Cacyt is the cytosolic Ca2+ concentration; Caer is the Ca2+ concentration in the endoplasmic reticulum; Cam is the Ca2+ concentration in mitochondria; and Catot is the conservation of total extracellular Ca2+. Details of the specific meaning and value of each parameter can be referred to in Table 1 and [24]. These parameters have the following relationships:
J p u m p = k p u m p C a c y t ,   P r t o t = P r + C a P r ,   J o u t = ( k o u t   C a c y t 2 K 3 2 + C a c y t 2 + k m   ) C a m ,   J i n = k i n   C a c y t 8 K 2 8 + C a c y t 8 , J c h = k c h   C a c y t 2 K 1 2 + C a c y t 2   C a e r C a c y t ,   C a P r = C a t o t C a c y t ρ e r β e r C a e r ρ m β m C a m , J l e a k = k l e a k C a e r C a c y t .

3. Results

3.1. Analysis of Stability

Catot corresponds to the total intracellular Ca2+ concentration including free Ca2+ in the cytosol, endoplasmic reticulum, and mitochondria. The model ignores the exchange between extracellular and cytosolic Ca2+ and focuses on the regulation of Ca2+ concentration in the cytosol and organelles. In cells, Ca2+ generally exists in the organelles as free or protein-bound. Ca2+ can be induced to increase or decrease either by binding to buffer proteins or by the dissociation of binding Ca2+.
There are two principal ways in which Ca2+ is regulated: absorption and liberation of Ca2+ by the endoplasmic reticulum and mitochondria. Therefore, during the process, a series of experimental manipulations can be carried out to change the Ca2+ concentration in organelles. Then, Catot will be changed in value. For example, intracellular injection of IP3 leads to the liberation of Ca2+ in the endoplasmic reticulum, or Ca2+ pump inhibitors are used to inhibit Ca2+ uptake from the endoplasmic reticulum, thereby increasing free cytosolic Ca2+. Based on the above reasons, Catot is a reasonable selection as the parameter that can be changed in the experiment. This is a practical guide for the experiment. Therefore, Catot was chosen as the bifurcation parameter to discuss the existence, quantity, type, and bifurcation of the equilibrium.
To make the calculations easier, we can make x = Cacyt, y = Caer, z = Cam, r = Catot. Therefore, model (1) can be transformed into the following expression:
d x d t = 0.01 r + 0.01 y 20.06 x 0.04 z 0.1 x ( x r + 4 y + 4 z + 120 ) 1500 x 2 ( x y ) x 2 + 25 300 x 8 x 8 + 0.16777216 + z ( 125 x 2 x 2 + 25 + 0.00625 ) d y d t = 5.0125 x 0.0125 y + 375 x 2 ( x y ) x 2 + 25 d z d t = 0.25 z ( 125 x 2 x 2 + 25 + 0.00625 ) + 75 x 8 x 8 + 0.16777216 .
One can immediately see that the equilibrium point of system (2) meets the following equation:
0.01 r 20.06 x + 0.01 y 0.04 z 0.1 x ( x r + 4 y + 4 z + 120 ) 1500 x 2 ( x y ) x 2 + 25 300 x 8 x 8 + 0.16777216 + z ( 125 x 2 x 2 + 25 + 0.00625 ) = 0 5.0125 x 0.0125 y + 375 x 2 ( x y ) x 2 + 25 = 0 0.25 z ( 125 x 2 x 2 + 25 + 0.00625 ) + 75 x 8 x 8 + 0.16777216 = 0  
First, by calculating
0.25 z ( 125 x 2 x 2 + 25 + 0.00625 ) + 75 x 8 x 8 + 0.16777216 = 0 ,
we can obtain
z = 75 x 8 ( 31.25 x 2 x 2 + 25 + 0.0015625 ) ( x 8 + 0.167772 ) .
Then,
5.0125 x 0.0125 y + 375 x 2 ( x y ) x 2 + 25 = 0
can be calculated to obtain
y = 30401 x 3 + 10025 x 30001 x 2 + 25 .
Consequently, substituting the expression for y and z obtains
0.01 r 20.06 x + 0.01 y 0.04 z 0.1 x ( x r + 4 y + 4 z + 120 ) 1500 x 2 ( x y ) x 2 + 25 300 x 8 x 8 + 0.16777216 + z ( 125 x 2 x 2 + 25 + 0.00625 ) = 0 .
Therefore, we have the following equation:
f ( x , r ) = 0.01 r 20.06 x 0.1 x ( x r + 4 σ 2 σ 3 + 300 x 8 σ 1 + 120 ) 3 x 8 σ 1 + 0.01 σ 2 σ 3 300 x 8 x 8 + 0.16777216 1500 x 2 ( x σ 2 σ 3 ) x 2 + 25 + 75 x 8 ( 125 x 2 x 2 + 25 + 0.00625 ) σ 1 = 0 y = 30401 x 3 + 10025 x 30001 x 2 + 25 z = 75 x 8 σ 1 ,  
where
σ 1 = ( 31.25 x 2 x 2 + 25 + 0.0015625 ) ( x 8 + 0.167772 ) ,   σ 2 = 30401 x 3 + 10025 x ,   σ 3 = 30001 x 2 + 25 .
Depending on the practical implications of x, y, z, and r, whether Equation (2) has an equilibrium point was considered to meet our special needs when r [0, 150]. Suppose (x0, y0, z0) is equilibrium. Let x1 = xx0, y1 = yy0, z1 = zz0, then we can obtain the following representations [30]:
d x 1 d t = 0.01 r 20.06 x 1 + x 0 + 0.01 y 1 + y 0 0.04 z 1 + z 0 0.1 x 1 + x 0 ( 4 ( y 1 + y 0 ) + 4 z 1 + z 0 + 120 + x 1 + x 0 r ) 1500 ( x 1 + x 0 ) 2 x 1 + x 0 y 1 y 0 ( x 1 + x 0 ) 2 + 25 300 ( x 1 + x 0 ) 8 ( x 1 + x 0 ) 8 + 0.16777216 + z 125 ( x 1 + x 0 ) 2 + 0.00625 ( x 1 + x 0 ) 2 + 25 ( x 1 + x 0 ) 2 + 25 d y 1 d t = 0.0125 y 1 + y 0 + 5.0125 x 1 + x 0 + 375 ( x 1 + x 0 ) 2 x 1 + x 0 y 1 y 0 ( x 1 + x 0 ) 2 + 25 d z 1 d t = 75 ( x 1 + x 0 ) 8 ( x 1 + x 0 ) 8 + 0.16777216 0.25 z 1 + z 0 125 ( x 1 + x 0 ) 2 ( x 1 + x 0 ) 2 + 25 + 0.00625 .
Obviously, (0, 0, 0) is the equilibrium of system (5), which, according to bifurcation theory, has the same features as that of the equilibrium of system (2) in terms of type, stability, and bifurcation type. The Jacobian matrix of system (5) can be easily calculated as follows:
J = ( b i j ) 3 × 3 = b 11 b 12 b 13 b 21 b 22 b 23 b 31 b 32 b 33 .
The linearized system corresponding to the equilibrium (0, 0, 0) of system (5) is
d x 1 d t = b 11 x 1 + b 12 y 1 + b 13 z 1 d y 1 d t = b 21 x 1 + b 22 y 1 + b 23 z 1 d z 1 d t = b 31 x 1 + b 32 y 1 + b 33 z 1 ,
where
b 11 = 0.1 r 0.2 x 0 0.4 y 0 0.4 z 0 32.06 2400 x 0 7 x 0 8 + 0.16777216 + 2400 x 0 15 ( x 0 8 + 0.16777216 ) 2 + z 0 ( 250 x 0 x 0 2 + 25 250 x 0 3 ( x 0 2 + 25 ) 2 ) 1500 x 0 2 x 0 2 + 25 3000 x 0 ( x 0 y 0 ) x 0 2 + 25 + 3000 x 0 3 ( x 0 y 0 ) ( x 0 2 + 25 ) 2 , b 12 = 0.4 x 0 + 0.01 + 1500 x 0 2 x 0 2 + 25 ,   b 13 = 0.4 x 0 0.03375 + 125 x 0 2 x 0 2 + 25 , b 21 = 5.0125 + 375 x 0 2 x 0 2 + 25 + 750 x 0 ( x 0 y 0 ) x 0 2 + 25 750 x 0 3 ( x 0 y 0 ) ( x 0 2 + 25 ) 2 ,   b 22 = 0.0125 375 x 0 2 x 0 2 + 25 , b 23 = 0 ,   b 31 = 600 x 0 7 x 0 8 + 0.16777216 600 x 0 15 ( x 0 8 + 0.16777216 ) 2 0.25 z 0 ( 250 x 0 x 0 2 + 25 250 x 0 3 ( x 0 2 + 25 ) 2 ) , b 32 = 0 ,   b 33 = 0.0015625 31.25 x 0 2 x 0 2 + 25 .
There is also a simple way to obtain the following eigen equation of the Jacobian matrix:
P ( λ ) = λ 3 + C 1 λ 2 + C 2 λ + C 3 = 0 ,
where
C 1 = b 11 b 22 b 33 , C 2 = b 11 b 22 + b 11 b 33 + b 22 b 33 b 13 b 31 b 12 b 21 b 32 b 23 , C 3 = b 31 b 13 b 22 + b 12 b 21 b 33 + b 32 b 23 b 11 b 11 b 22 b 33 b 12 b 23 b 31 b 13 b 21 b 32 .
The following formula is obtained by calculating:
C 1 = 0.1 r + 0.2 x 0 + 0.4 y 0 + 0.4 z 0 + 2400 x 0 7 x 0 8 + 0.16777216 2400 x 0 15 σ 4 z 0 σ 5 + 1906.25 x 0 2 x 0 2 + 25 + 32.0740625 + 3000 x 0 ( x 0 y 0 ) x 0 2 + 25 3000 x 0 3 ( x 0 y 0 ) σ 6 , C 2 = ( σ 3 + 0.0125 ) σ 7 + σ 1 σ 7 + ( σ 3 + 0.0125 ) σ 1 ( σ 2 0.4 x 0 + 0.01 ) ( σ 3 + 750 x 0 ( x 0 y 0 ) x 0 2 + 25 750 x 0 3 ( x 0 y 0 ) σ 6 + 5.0125 ) ( 0.4 x 0 125 x 0 2 x 0 2 + 25 + 0.03375 ) ( 600 x 0 15 σ 4 600 x 0 7 x 0 8 + 0.16777216 + 0.25 z 0 σ 5 ) , C 3 = σ 1 σ 3 + 0.0125 ( 0.2 x 0 0.1 r + 0.4 y 0 + 0.4 z 0 + 2400 x 0 7 x 0 8 + 0.16777216 2400 x 0 15 σ 4 z 0 σ 5 + σ 2 + 3000 x 0 x 0 y 0 x 0 2 + 25 3000 x 0 3 x 0 y 0 σ 6 + 32.06 ) σ 3 + 0.0125 ( 0.4 x 0 125 x 0 2 x 0 2 + 25 + 0.03375 ) 600 x 0 15 σ 4 600 x 0 7 x 0 8 + 0.16777216 + 0.25 z 0 σ 5 σ 1 σ 2 0.4 x 0 + 0.01 ( σ 3 + 750 x 0 x 0 y 0 x 0 2 + 25 750 x 0 3 x 0 y 0 σ 6 + 5.0125 ) ,
where
σ 1 = 31.25 x 0 2 x 0 2 + 25 + 0.0015625 ,   σ 2 = 1500 x 0 2 x 0 2 + 25 ,   σ 3 = 375 x 0 2 x 0 2 + 25 , σ 4 = ( x 0 8 + 0.16777216 ) 2 ,   σ 5 = 250 x 0 x 0 2 + 25 250 x 0 3 σ 6 ,   σ 6 = ( x 0 2 + 25 ) 2 , σ 7 = 0.2 x 0 0.1 r + 0.4 y 0 + 0.4 z 0 + 2400 x 0 7 x 0 8 + 0.16777216 2400 x 0 15 σ 4 z 0 σ 5 + σ 2 + 3000 x 0 ( x 0 y 0 ) x 0 2 + 25 3000 x 0 3 ( x 0 y 0 ) σ 6 + 32.06 .
The Hurwitz matrix is obtained as follows.
ϕ 1 = C 1 ,   ϕ 2 = C 1 1 C 3 C 2 ,   ϕ 3 = C 1 1 0 C 3 C 2 1 0 0 C 3 .
If the determinants of the Hurwitz matrix are all positive, then it can be confirmed that the system is stable.
det ( ϕ i ) > 0 , i = 1 , 2 , 3 .
The specific expression of the Hurwitz criteria is
C 1 > 0 ,   C 3 > 0 ,   C 1 C 2 > C 3 .
According to this inequality and MATLAB, the values of r can be obtained:
r 1 = 58.61 , r 2 = 101.83 .
In general, the system has a stable node when the eigenvalues are all positive. When the eigenvalues have positive and negative values, the system has saddle. When the eigenvalues are a pair of pure virtual roots, the equilibrium point is a non-hyperbolic equilibrium point. The type of equilibrium points can be determined by deriving the eigenvalues with different r.
The investigation related to differential equations provides a theoretical basis for the conclusions drawn below:
(1)
r < 58.61, system (2) is characterized by a unique equilibrium, and the equilibrium is stable (stable node);
(2)
r = 58.61, system (2) has a unique equilibrium M1 = (0.0454, 5.2741, 0), and is a non-hyperbolic equilibrium;
(3)
58.61 < r < 101.83, system (2) is characterized by a unique equilibrium, and the equilibrium is unstable (saddle);
(4)
r = 101.83, system (2) has a unique equilibrium M2 = (0.3533, 1.2951, 0.6915), and is a non-hyperbolic equilibrium; and
(5)
r > 101.83, system (2) is characterized by a unique equilibrium, and the equilibrium is stable (stable node).

3.2. Bifurcation of Equilibria

According to the above, when r is divided into 58.61 and 101.83, the equilibria of system (2) are both non-hyperbolic. Therefore, in the following, we will analyze the dynamics near these equilibrium points.
When a system with parameters uses the center manifold theory to reduce dimensions, the parameters are required as new variables of the system. After shifting the equilibriums, it is obvious that the equilibriums of system (2) at r1 = 0 is M (x1, y1, z1) = (0, 0, 0). Taking r as another dynamic variable and adding dr1/dt = 0 into system (5), we assume that when r1 = rr0 when r = r0, system (5) can be sorted out as follows [30]:
d x 1 d t = 0.01 r 1 + r 0 20.06 x 1 + x 0 + 0.01 y 1 + y 0 0.04 z 1 + z 0 0.1 x 1 + x 0 ( 4 ( y 1 + y 0 ) + 4 z 1 + z 0 + 120 + x 1 + x 0 r 1 r 0 ) 1500 ( x 1 + x 0 ) 2 x 1 + x 0 y 1 y 0 ( x 1 + x 0 ) 2 + 25 300 ( x 1 + x 0 ) 8 ( x 1 + x 0 ) 8 + 0.16777216 + z 125 ( x 1 + x 0 ) 2 + 0.00625 ( x 1 + x 0 ) 2 + 25 ( x 1 + x 0 ) 2 + 25 . d y 1 d t = 0.0125 y 1 + y 0 + 5.0125 x 1 + x 0 + 375 ( x 1 + x 0 ) 2 x 1 + x 0 y 1 y 0 ( x 1 + x 0 ) 2 + 25 d z 1 d t = 75 ( x 1 + x 0 ) 8 ( x 1 + x 0 ) 8 + 0.16777216 0.25 z 1 + z 0 125 ( x 1 + x 0 ) 2 ( x 1 + x 0 ) 2 + 25 + 0.00625 d r 1 d t = 0
System (6) has entirely the same features as the equilibriums of system (2). Next, the specific r0 value was analyzed. For r0 = 58.61, (0, 0, 0, 0) is the corresponding equilibrium of system (6). According to the calculation, we can easily obtain that the eigenvalue of the equilibrium of system (6) is η1 = −0.0041, η2 = 0.4880i, η3 = −0.4880i, η4 = 0. The eigenvector is
0.0065 0.0193 0.2287 i 0.0193 + 0.2287 i 0.0026 0.3411 0.9733 0.9733 0.1258 0.9399 0 0 0 0 0 0 0.9921 .
Suppose
x 1 y 1 z 1 r 1 = T u v w s ,
where
T = 0.0065 0.0193 0.2287 0.0026 0.3411 0.9733 0 0.1258 0.9399 0 0 0 0 0 0 0.9921 .
System (6) can vary in shape such as the following:
u ˙ v ˙ w ˙ s ˙ = 0.0041 0 0 0 0 0 0.4880 0 0 0.4880 0 0 0 0 0 0 u v w s + g 1 g 2 g 3 g 4
and
x ˙ 1 y ˙ 1 z ˙ 1 r ˙ 1 = T u ˙ v ˙ w ˙ s ˙ u ˙ v ˙ w ˙ s ˙ = T 1 x ˙ 1 y ˙ 1 z ˙ 1 r ˙ 1 = T 1 f 1 f 2 f 3   f 4 ,
where
f 1 = 0.01 q 14 20.06 q 11 + 0.01 q 12 0.04 q 13 0.1 q 11 ( q 11 q 14 + 4 q 12 + 4 q 13 + 120 ) 1500 q 11 2 ( q 11 q 12 ) q 11 2 + 25 300 q 11 8 q 11 8 + 0.16777216 + q 13 ( 125 q 11 2 q 11 2 + 25 + 0.00625 ) , f 2 = 5.0125 q 11 0.0125 q 12 + 375 q 11 2 ( q 11 q 12 ) q 11 2 + 25 , f 3 = 0.25 q 13 ( 125 q 11 2 q 11 2 + 25 + 0.00625 ) + 75 q 11 8 q 11 8 + 0.16777216 ,   f 4 = 0 , q 11 = 0.0026 s 0.0065 u 0.0193 v + 0.2287 w + 0.0454 , q 12 = 0.3411 u 0.1258 s + 0.9733 v + 5.2741 , q 13 = 0.9399 u ,   q 14 = 0.9921 s .
Furthermore,
g 1 g 2 g 3 g 4 = T 1 f 1 f 2 f 3 f 4 0.0041 0 0 0 0 0 0.4880 0 0 0.4880 0 0 0 0 0 0 u v w s ,
where
T 1 = 0 0 1.063942973 0 0 1.027432446 0.3728664831 0.1302802154 4.372540446 0.08670505559 0.001227344998 0.0004647809355 0 0 0 1.007962907 .
After the calculation, we can obtain the following equations:
g 1 = 1.063942973 f 3 + 0.0041 u , g 2 = 1.027432446 f 2 0.3728664831 f 3 + 0.1302802154 f 4 + 0.488 w , g 3 = 4.372540446 f 1 + 0.08670505559 f 2 0.001227344998 f 3 0.0004647809355 f 4 0.488 v , g 4 = 0 .
Due to the center manifold theory, it can be conclusively demonstrated that system (6) has a central manifold, which can be expressed as follows:
W l o c c ( M ) = u , v , w , s R 4 | u = h v , w , s ,   h 0 , 0 , 0 = 0 ,   D h 0 , 0 , 0 = 0 .
Let h1(v, w, s) = a1v2 + b1w2 + c1s2 + d1vw + e1vs + f1ws + , the center manifold of system (6) can be expressed as follows:
N ( h 1 ) = D h 1 v ˙ w ˙ s ˙ + 0.0041 h 1 g 1 0 .
Therefore, the high-order partial derivatives can be applied to obtain the values of a1 to f1. The equation is given below:
0.0006 0 0 0.9764 0 0 0 0.0030 0 0.9760 0 0 0 0 0.0083 0 0.0001 0 0.9757 0.9759 0 0.0006 0 0 0.0001 0 0 0 0.0019 0.4880 0 0 0 0 0.4880 0.0028 a 1 b 1 c 1 d 1 e 1 f 1 = 0 .
From the center manifold theory, one can see that
a 1 = 0.005121190007 ,   b 1 = 0.005120438342 ,   c 1 = 0.0000001949388882 , d 1 = 0.000003482559566 ,   e 1 = 0.000000272409307 ,   f 1 = 0.000001352303558 .
After downscaling system (6), it will be confined to a two-dimensional system as follows:
v ˙ w ˙ = 0 0.4880 0.4880 0 v w + B 1 ( v , w ) B 2 ( v , w ) ,
where
B 1 ( v , w ) = 0.01500565088 s 0.111890991 v + 1.665806175 w + , B 2 ( v , w ) = 1.237979119 v 0.1889085743 s 19.9606049 w + .
Hence, it is easy to verify that
a = 1 16 [ B v v v 1 + B v w w 1 + B v v w 2 + B w w w 2 ] | ( 0 , 0 , 0 ) + 1 16 × 0.4880 [ B v w 1 ( B v v 1 + B w w 1 ) B v w 2 ( B v v 2 + B w w 2 ) B v v 1 B v v 2 + B w w 1 B w w 2 ] | ( v = 0 , s = 0 , w = 0 ) = 25.8314 < 0 , d = d Re ( η ( s ) ) d s | ( 0 , 0 , 0 ) = 0.0500 > 0 .
As a consequence of the above discussion, the following conclusion is summarized.
Conclusion 1: When r0 = 58.61, a supercritical Hopf bifurcation occurs at the equilibrium M1 = (0.0454, 5.2741, 0). With the increase in r, when r > r0, the equilibrium changes from stable to unstable and loses its stability, whereas a stable periodic solution occurs near the equilibrium point. System (2) begins to oscillate.
For r0 = 101.83, we can easily obtain that the eigenvalue of the equilibrium of system (6) is η1 = −0.0762, η2 = 4.1872i, η3 = −4.1872i, η4 = 0, and the eigenvector is
0.0435 0.7846 0.7846 0.0099 0.0727 0.2101 0.4697 i 0.2101 + 0.4697 i 0.0159 0.9964 0.0127 + 0.3457 i 0.0127 0.3457 i 0.1169 0 0 0 0.9929 .
Suppose
x 1 y 1 z 1 r 1 = T u v w s ,
where
T = 0.0435 0.7846 0 0.0099 0.0727 0.2101 0.4697 0.0159 0.9964 0.0127 0.3457 0.1169 0 0 0 0.9929 .
System (6) can vary in shape such as the following:
u ˙ v ˙ w ˙ s ˙ = 0.0762 0 0 0 0 0 4.1872 0 0 4.1872 0 0 0 0 0 0 u v w s + g 1 g 2 g 3 g 4
and
x ˙ 1 y ˙ 1 z ˙ 1 r ˙ 1 = T u ˙ v ˙ w ˙ s ˙ u ˙ v ˙ w ˙ s ˙ = T 1 x ˙ 1 y ˙ 1 z ˙ 1 r ˙ 1 = T 1 f 1 f 2 f 3   f 4 ,
where
f 1 = 0.01 q 14 20.06 q 11 + 0.01 q 12 0.04 q 13 0.1 q 11 ( q 11 q 14 + 4 q 12 + 4 q 13 + 120 ) 1500 q 11 2 ( q 11 q 12 ) q 11 2 + 25 300 q 11 8 q 11 8 + 0.16777216 + q 13 ( 125 q 11 2 q 11 2 + 25 + 0.00625 ) , f 2 = 5.0125 q 11 0.0125 q 12 + 375 q 11 2 ( q 11 q 12 ) q 11 2 + 25 , f 3 = 0.25 q 13 ( 125 q 11 2 q 11 2 + 25 + 0.00625 ) + 75 q 11 8 q 11 8 + 0.16777216 , f 4 = 0 , q 11 = 0.0099 s 0.0435 u 0.7846 v + 0.3533 , q 12 = 0.0727 u + 0.2101 v + 0.4697 w 0.1258 s + 1.2951 , q 13 = 0.9964 u 0.0127 v 0.3457 w + 0.1169 s + 0.6915 , q 14 = 0.9929 s .
Furthermore,
g 1 g 2 g 3 g 4 = T 1 f 1 f 2 f 3 f 4 0.0762 0 0 0 0 0 4.1872 0 0 4.1872 0 0 0 0 0 0 u v w s ,
where
T 1 = 0.1902682299 0.7741178608 1.051788138 0.1133338854 1.263985893 0.04291884648 0.05831351516 0.006424634991 0.594839124 2.22963832 0.1367113649 0.01367789647 0 0 0 1.00715077 .
After calculation, we can obtain the following equations:
g 1 = 0.1902682299 f 1 + 0.7741178608 f 2 + 1.051788138 f 3 0.1133338854 f 4 + 0.0762 u , g 2 = 1.263985893 f 1 + 0.04291884648 f 2 + 0.05831351516 f 3 + 0.006424634991 f 4 + 4.1872 w , g 3 = 0.594839124 f 1 + 2.22963832 f 2 + 0.1367113649 f 3 + 0.01367789647 f 4 4.1872 v , g 4 = 0 .
Let h2 (v, w, s) = a2v2 + b2w2 + c2s2 + d2vw + e2vs + f2ws + …, the center manifold of system (6) can be expressed as follows:
N ( h 2 ) = D h 2 v ˙ w ˙ s ˙ + 0.0762 h 2 g 1 0 .
Therefore, the high-order partial derivatives can be applied to obtain the values of a2 to f2. The equation is given below:
0.1678 0 0 25.1159 0 0 0 0.1401 0 25.1209 0 0 0 0 0.1523 0 0.0003 0.0003 25.1209 25.1159 0 0.0769 0 0 0.0003 0 0 0.0001 0.0800 12.5579 0 0.0003 0 0.0001 12.5604 0.0731 a 2 b 2 c 2 d 2 e 2 f 2 = 0 .
From the center manifold theory, one can know that
a 2 = 27.10743425 , b 2 = 27.11200115 , c 2 = 0.006549229838 , d 2 = 0.1512035103 , e 2 = 0.0006281394405 , f 2 = 0.006719709047 .
After downscaling system (5), it will be confined to a two-dimensional system as follows:
v ˙ w ˙ = 0 4.1872 4.1872 0 v w + B 1 ( v , w ) B 2 ( v , w ) ,
where
B 1 ( v , w ) = 0.2467195946 s 20.06617206 v 4.210867326 w + , B 2 ( v , w ) = 4.776376224 v 0.004015215297 s 0.002071369214 w + .
Hence, it is easy to verify that
a = 1 16 [ B v v v 1 + B v w w 1 + B v v w 2 + B w w w 2 ] | ( 0 , 0 , 0 ) + 1 16 × 4.1872 [ B v w 1 ( B v v 1 + B w w 1 ) B v w 2 ( B v v 2 + B w w 2 ) B v v 1 B v v 2 + B w w 1 B w w 2 ] | ( v = 0 , s = 0 , w = 0 ) = 66.0608 < 0 , d = d Re ( η ( s ) ) d s | ( 0 , 0 , 0 ) = 0.0495 > 0 .

4. Numerical Simulations

Variations of Catot and Cacyt were analyzed from the point of view of bifurcation dynamics. The system presents an equilibrium state and an oscillation state. The equilibrium bifurcation diagram of system (2) is schematically illustrated in Figure 1. In the figure, stable equilibrium points are denoted as the solid line, and unstable equilibrium points are denoted by the dashed line. Additionally, red hollow circles represent unstable limit cycles, and red-filled circles indicate stable limit cycles. Two Hopf bifurcations, namely HB1 and HB2, can also be found with corresponding values of r1 = 58.61 and r2 = 101.83. When r passes through r1 = 58.61 and r2 = 101.83, two supercritical bifurcations occur in the system.
First, the stability of equilibria is discussed. The stability undergoes a series of variations. To be specific, from r = 0 to 58.61 and from r = 101.83 to 150, there were stable equilibria. Between r = 58.61 and 101.83, there were unstable equilibria. This is followed by a discussion of the stability of limit cycles. A stable limit cycle formed at r = 58.61. When r increased to 58.62, saddle-node bifurcation (LP) of limit cycles occurred, and stable and unstable limit cycles had a meeting at this point. Thus, stable limit cycles became unstable. Similarly, LP of the limit cycles was at r = 65.09 and 89.96. The stability of limit cycles was constantly changing. Limit cycles go from unstable to stable. Then, unstable limit cycles were obtained through r = 89.96. At different values of r, for example, r = 101.8, torus bifurcation (TR) was in the system and period doubling bifurcation (PD) could also be found when r = 101.3. Period doubling bifurcation is a typical approach to chaos, which can be considered as a way to enter chaos from period windows.
When the limit cycles pass through TR, the stability of the system limit cycles transforms further, and unstable limit cycles of the system turn into stable limit cycles. Torus bifurcation (TR) is a way to cause chaos in the system and is shown in the following images. Figure 2 and Figure 3 represent the bifurcation diagram in the (r, y) and (r, z) planes, respectively. One can easily see that two Hopf bifurcation points appeared due to the variation in parameter r.
Some dynamic behaviors of system (2) are schematically presented in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10. Figure 4a, Figure 5a, Figure 6a, Figure 7a, Figure 8a, Figure 9a and Figure 10a represent the temporal evolution for different values of parameter r. Figure 4b, Figure 5b, Figure 6b, Figure 7b, Figure 8b, Figure 9b and Figure 10b correspond to different state trajectories in 3D phase space by varying different values of r. Among these, Figure 4, Figure 5, Figure 6, Figure 8 and Figure 9 represent the regular Ca2+ oscillations. Illustrations of chaos can be observed in Figure 7 and Figure 10. These phenomena are abundant and need further study.
Figure 11 is the return map for when r = 101.8. In the figure, the pioneer points x were compared to the successor points x (n + 1). When the track repeatedly traverses the same section, the points on the left of the section are very scattered. The image shows a scatterplot. Therefore, there is chaos in the system at this point. Furthermore, when the images drawn in Figure 11 were observed together with torus bifurcation (TR) and period doubling bifurcation (PD) in Figure 1 and the sequence diagram of time and phase diagram in Figure 10, we can see that the characteristics of system (2) displayed in these images were consistent and chaos was found in all of these diagrams, further validating the availability of the method.
Figure 12 investigates the Lyapunov exponential spectrum of system (2). Lyapunov exponent is a description of the stability of a dynamic system. When the Lyapunov exponent is greater than zero, it means chaos occurs in the system. For n-dimensional continuous dynamical systems x ˙ = g (x), the system forms an n-dimensional sphere with x0 as the center and δ x(x0,0) as the radius when t = 0. With the evolution of time, the sphere deforms into an n-dimensional ellipsoid at time t. Assume that the length of the half-axis in the direction of the ith-axis of the ellipsoid is δ xi (x0, t) , then the ith-Lyapunov exponent of the system is
λ i = lim t 1 t ln δ x i ( x 0 , t ) δ x ( x 0 , 0 ) .
Solving the formula will give the Lyapunov exponent. In Figure 12, the yellow curves represent LE1, the green curves represent LE2, and the blue curves represent LE3. The Lyapunov exponent was more than 0 within a certain parameter range, indicating the presence of chaos in system (2). For example, when r = 101.8, the corresponding Lyapunov exponents were LE1 = 0.1378, LE2 = −0.4361, and LE3 = −27.6405, respectively. LE1 > 0, so the system was in chaos at this point, which is consistent with Figure 1, Figure 10, and Figure 11.

5. Summary

In this study, the stability and bifurcation phenomena of a 3D Ca2+ oscillation model were investigated theoretically and simulated numerically. Oscillations of Ca2+ in the cytosol, endoplasmic reticulum, and mitochondria were studied using the total concentration of Ca2+ in the cytosol, mitochondria, and endoplasmic reticulum as the bifurcation parameter. Two Hopf bifurcation points were obtained after varying the value of the parameter for investigation. There was a stable periodic solution near the supercritical Hopf bifurcation point. The Hopf bifurcation point was the source of oscillation.
After the theoretical analysis of the system, the correctness of the theoretical results was verified by the numerical method. During the numerical simulation, there were some interesting oscillations. In addition to the usual simple oscillations, there were also some Ca2+ oscillations of bursting. The peak values and periods of oscillations of different bursting were also different, which may be related to bifurcation. Such findings provided ideas for further research. Moreover, by exploring the bifurcation diagram, the Hopf bifurcation case calculated in this study well matched that of the image. Moreover, there were also limit cycle LP and TR points. Where the system occurred, there was chaos at TR. This provided further insights into the dynamic behavior between the two points where the Hopf bifurcation occurs. In general, the total Ca2+ concentration greatly influences the formation and characteristics of Ca2+ oscillations in cells.
Some of the special phenomena found above need to be explored in more detail. Future works should select different models, explore different models, and conduct in-depth research. In addition, in many intracellular dynamic models, time delay has a certain effect on the dynamic behavior of cells. However, the model established in this study did not consider the effect of time delay. Therefore, in future works, the study of cell models with time delay is also worth considering.

Author Contributions

Conceptualization, Q.J.; methodology, X.Q.; software, X.Q.; validation, Q.J.; formal analysis, Q.J.; investigation, Q.J. and X.Q.; resources, Q.J. and X.Q.; data curation, X.Q.; writing—original draft preparation, X.Q.; writing—review and editing, Q.J. and X.Q.; visualization, X.Q.; supervision, Q.J.; project administration, Q.J.; funding acquisition, Q.J. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Natural Science Foundation of China under Grant Nos. 12062004, 11872084 and the Guangxi Natural Science Foundation under Grant No. GXNSFAA297240.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All the data utilized in this article have been included, and the sources from where they were adopted were cited accordingly.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Jaiswal, J.K. Calcium—How and why? J. Biosci. 2001, 26, 357–363. [Google Scholar] [CrossRef]
  2. Clapham, D.E. Calcium signaling. Cell 2007, 131, 1047–1058. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Leake, I. Microglial calcium activity in ischaemic brain injury. Nat. Rev. Neurol. 2021, 17, 2. [Google Scholar] [CrossRef]
  4. Bootman, M.D.; Collins, T.J.; Peppiatt, C.M.; Prothero, L.S.; Mackenzie, L.; Smet, P.D.; Travers, M.; Tovey, S.C.; Seo, J.T.; Berridge, M.J.; et al. Calcium signalling—An overview. In Seminars in Cell & Developmental Biology; Academic Press: Cambridge, MA, USA, 2001; Volume 12, pp. 3–10. [Google Scholar]
  5. Toyoshima, C.; Nomura, H. Structural changes in the calcium pump accompanying the dissociation of calcium. Nature 2002, 418, 605–611. [Google Scholar] [CrossRef] [PubMed]
  6. McAinsh, M.R.; Pittman, J.K. Shaping the calcium signature. New Phytol. 2009, 181, 275–294. [Google Scholar] [CrossRef]
  7. Blum, I.D.; Keleş, M.F.; Baz, E.-S.; Han, E.; Park, K.; Luu, S.; Issa, H.; Brown, M.; Ho, M.C.W.; Tabuchi, M.; et al. Astroglial calcium signaling encodes sleep need in Drosophila. Curr. Biol. 2021, 31, 150–162.e7. [Google Scholar] [CrossRef] [PubMed]
  8. Novikova, I.N.; Manole, A.; Zherebtsov, E.A.; Stavtsev, D.D.; Vukolova, M.N.; Dunaev, A.V.; Angelova, P.R.; Abramov, A.Y. Adrenaline induces calcium signal in astrocytes and vasoconstriction via activation of monoamine oxidase. Free Radic. Biol. Med. 2020, 159, 15–22. [Google Scholar] [CrossRef] [PubMed]
  9. Lewis, K.J.; Frikha-Benayed, D.; Louie, J.; Stephen, S.; Spray, D.C.; Thi, M.M.; Seref-Ferlengez, Z.; Majeska, R.J.; Weinbaum, S.; Schaffler, M.B. Osteocyte calcium signals encode strain magnitude and loading frequency in vivo. Proc. Natl. Acad. Sci. USA 2017, 114, 11775–11780. [Google Scholar] [CrossRef] [Green Version]
  10. Siddiqui, M.H.; Al-Whaibi, M.H.; Basalah, M.O. Interactive effect of calcium and gibberellin on nickel tolerance in relation to antioxidant systems in Triticum aestivum L. Protoplasma 2011, 248, 503–511. [Google Scholar] [CrossRef] [PubMed]
  11. Lautner, S.; Fromm, J. Calcium-dependent physiological processes in trees. Plant Biol. 2010, 12, 268–274. [Google Scholar] [CrossRef]
  12. Ahmad, P.; Sarwat, M.; Bhat, N.A.; Wani, M.R.; Kazi, A.G.; Tran, L.-S.P. Alleviation of cadmium toxicity in Brassica juncea L. (Czern. & Coss.) by calcium application involves various physiological and biochemical strategies. PLoS ONE 2015, 10, e0114571. [Google Scholar]
  13. Shi, J.Q.; Wu, Z.X.; Song, L.R. Physiological and molecular responses to calcium supplementation in Microcystis aeruginosa (Cyanobacteria). N. Z. J. Mar. Freshw. Res. 2013, 47, 51–61. [Google Scholar] [CrossRef] [Green Version]
  14. Endo, M. Calcium-induced calcium release in skeletal muscle. Physiol. Rev. 2009, 89, 1153–1176. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Carter, A.G.; Vogt, K.E.; Foster, K.A.; Regehr, W.G. Assessing the role of calcium-induced calcium release in short-term presynaptic plasticity at excitatory central synapses. J. Neurosci. 2002, 22, 21–28. [Google Scholar] [CrossRef] [Green Version]
  16. Stern, M.D.; Cheng, H. Putting out the fire: What terminates calcium-induced calcium release in cardiac muscle. Cell Calcium 2004, 35, 591–601. [Google Scholar] [CrossRef]
  17. Wood, C.D.; Darszon, A.; Whitaker, M. Speract induces calcium oscillations in the sperm tail. J. Cell Biol. 2003, 161, 89–101. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Wood, A.; Wing, M.G.; Benham, C.D. Specific induction of intracellular calcium oscillations by complement membrane attack on oligodendroglia. J. Neurosci. 1993, 13, 3319–3332. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Borghans, J.A.M.; Dupont, G.; Goldbeter, A. Complex intracellular calcium oscillations A theoretical exploration of possible mechanisms. Biophys. Chem. 1997, 66, 25–41. [Google Scholar] [CrossRef] [Green Version]
  20. Dupont, G.; Houart, G.; De Koninck, P. Sensitivity of CaM kinase II to the frequency of Ca2+ oscillations: A simple model. Cell Calcium 2003, 34, 485–497. [Google Scholar] [CrossRef]
  21. Chay, T.R. Electrical bursting and luminal calcium oscillation in excitable cell models. Biol. Cybern. 1996, 75, 419–431. [Google Scholar] [CrossRef]
  22. Perc, M.; Marhl, M. Different types of bursting calcium oscillations in non-excitable cells. Chaos Solitons Fractals 2003, 18, 759–773. [Google Scholar] [CrossRef] [Green Version]
  23. Shen, P.; Larter, R. Chaos in intracellular Ca2+ oscillations in a new model for non-excitable cells. Cell Calcium 1995, 17, 225–232. [Google Scholar] [CrossRef]
  24. Marhl, M.; Haberichter, T.; Brumen, M.; Heinrich, R. Complex calcium oscillations and the role of mitochondria and cytosolic proteins. BioSystems 2000, 57, 75–86. [Google Scholar] [CrossRef]
  25. Gellerich, F.N.; Gizatullina, Z.; Trumbeckaite, S.; Nguyen, H.P.; Pallaas, T.; Arandarcikaite, O.; Vielhaber, S.; Seppet, E.; Striggow, F. The regulation of OXPHOS by extramitochondrial calcium. Biochim. Biophys. Acta (BBA)-Bioenerg. 2010, 1797, 1018–1027. [Google Scholar] [CrossRef] [Green Version]
  26. Grubelnik, V.; Larsen, A.Z.; Kummer, U.; Olsen, L.F.; Marhl, M. Mitochondria regulate the amplitude of simple and complex calcium oscillations. Biophys. Chem. 2001, 94, 59–74. [Google Scholar] [CrossRef]
  27. Carr, J. Applications of Centre Manifold Theory; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  28. Erhardt, A.H. Bifurcation analysis of a certain Hodgkin-Huxley model depending on multiple bifurcation parameters. Mathematics 2018, 6, 103. [Google Scholar] [CrossRef] [Green Version]
  29. Sun, M.; Li, Y.; Yao, W. A Dynamic Model of Cytosolic Calcium Concentration Oscillations in Mast Cells. Mathematics 2021, 9, 2322. [Google Scholar] [CrossRef]
  30. Zuo, H.; Ye, M. Bifurcation and Numerical Simulations of Ca2+ Oscillatory Behavior in Astrocytes. Front. Phys. 2020, 8, 258. [Google Scholar] [CrossRef]
  31. Etter, D.M.; Kuncicky, D.C.; Hull, D.W. Introduction to MATLAB; Prentice Hall: Hoboken, NJ, USA, 2002. [Google Scholar]
  32. Venkataraman, P. Applied Optimization with MATLAB Programming; John Wiley & Sons: Hoboken, NJ, USA, 2009. [Google Scholar]
Figure 1. Curve of dynamical bifurcation of system (2) with r in the (r, x) plane. HB1 and HB2 refer to the Hopf bifurcations. LP refers to the saddle-node bifurcation. TR refers to the torus bifurcation. PD refers to the period doubling bifurcation.
Figure 1. Curve of dynamical bifurcation of system (2) with r in the (r, x) plane. HB1 and HB2 refer to the Hopf bifurcations. LP refers to the saddle-node bifurcation. TR refers to the torus bifurcation. PD refers to the period doubling bifurcation.
Mathematics 09 03324 g001
Figure 2. Curve of the dynamical bifurcation of system (2) with r in the (r, y) plane. HB1 and HB2 refer to the Hopf bifurcations.
Figure 2. Curve of the dynamical bifurcation of system (2) with r in the (r, y) plane. HB1 and HB2 refer to the Hopf bifurcations.
Mathematics 09 03324 g002
Figure 3. Curve of dynamical bifurcation of system (2) with r in the (r, z) plane. HB1 and HB2 refer to the Hopf bifurcations.
Figure 3. Curve of dynamical bifurcation of system (2) with r in the (r, z) plane. HB1 and HB2 refer to the Hopf bifurcations.
Mathematics 09 03324 g003
Figure 4. (a) Temporal evolution for r = 60. (b) State trajectories in 3D phase space for r = 60.
Figure 4. (a) Temporal evolution for r = 60. (b) State trajectories in 3D phase space for r = 60.
Mathematics 09 03324 g004
Figure 5. (a) Temporal evolution for r = 70. (b) State trajectories in 3D phase space for r = 70.
Figure 5. (a) Temporal evolution for r = 70. (b) State trajectories in 3D phase space for r = 70.
Mathematics 09 03324 g005
Figure 6. (a) Temporal evolution for r = 86. (b) State trajectories in 3D phase space for r = 86.
Figure 6. (a) Temporal evolution for r = 86. (b) State trajectories in 3D phase space for r = 86.
Mathematics 09 03324 g006
Figure 7. (a) Temporal evolution for r = 96. (b) State trajectories in 3D phase space for r = 96.
Figure 7. (a) Temporal evolution for r = 96. (b) State trajectories in 3D phase space for r = 96.
Mathematics 09 03324 g007
Figure 8. (a) Temporal evolution for r = 100.8. (b) State trajectories in 3D phase space for r = 100.8.
Figure 8. (a) Temporal evolution for r = 100.8. (b) State trajectories in 3D phase space for r = 100.8.
Mathematics 09 03324 g008
Figure 9. (a) Temporal evolution for r = 101.64. (b) State trajectories in 3D phase space for r = 101.64.
Figure 9. (a) Temporal evolution for r = 101.64. (b) State trajectories in 3D phase space for r = 101.64.
Mathematics 09 03324 g009
Figure 10. (a) Temporal evolution for r = 101.8. (b) State trajectories in 3D phase space for r = 101.8.
Figure 10. (a) Temporal evolution for r = 101.8. (b) State trajectories in 3D phase space for r = 101.8.
Mathematics 09 03324 g010
Figure 11. Return map of Cacyt with r = 101.8.
Figure 11. Return map of Cacyt with r = 101.8.
Mathematics 09 03324 g011
Figure 12. Lyapunov index of system (2) with parameter r.
Figure 12. Lyapunov index of system (2) with parameter r.
Mathematics 09 03324 g012
Table 1. Unless otherwise specified, the parameter values are used as given in this model.
Table 1. Unless otherwise specified, the parameter values are used as given in this model.
Prtot120 μMkpump20 s−1
ρer0.01kout125 s−1
ρm0.01km0.00625 s−1
βer0.0025k0.01 s−1
βm0.0025K15 μM
kch1500 s−1K20.8 μM
kleak0.05 s−1K35 μM
kin300 μM−1 s−1k+0.1 μM−1 s−1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Qie, X.; Ji, Q. Computational Analysis and Bifurcation of Regular and Chaotic Ca2+ Oscillations. Mathematics 2021, 9, 3324. https://doi.org/10.3390/math9243324

AMA Style

Qie X, Ji Q. Computational Analysis and Bifurcation of Regular and Chaotic Ca2+ Oscillations. Mathematics. 2021; 9(24):3324. https://doi.org/10.3390/math9243324

Chicago/Turabian Style

Qie, Xinxin, and Quanbao Ji. 2021. "Computational Analysis and Bifurcation of Regular and Chaotic Ca2+ Oscillations" Mathematics 9, no. 24: 3324. https://doi.org/10.3390/math9243324

APA Style

Qie, X., & Ji, Q. (2021). Computational Analysis and Bifurcation of Regular and Chaotic Ca2+ Oscillations. Mathematics, 9(24), 3324. https://doi.org/10.3390/math9243324

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop