1. Introduction
Ca
2+ 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]. Ca
2+ 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 Ca
2+ concentration remains stable without stimulation. Some Ca
2+ enter the cytosol from the Ca
2+ store when the cell is stimulated, leading to an increment in Ca
2+ concentration that triggers a series of physiological activities [
5,
6,
7,
8,
9,
10,
11,
12,
13]. The cell adjusts itself and closes the Ca
2+ channels when the stimulation increases to a certain extent, and then Ca
2+ is taken up by Ca
2+ stores. After a while, the intracellular Ca
2+ concentration gains stability, and the cell returns to a resting state. The source of intracellular Ca
2+ is mainly through the influx of extracellular Ca
2+ and the release of Ca
2+ stored in the endoplasmic reticulum or sarcoplasmic reticulum. The latter is mainly based on the mechanism of Ca
2+ release from intracellular Ca
2+ stores triggered by Ca
2+ signaling, namely Ca
2+-induced Ca
2+ 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], Ca
2+ oscillations have been extensively studied. Ca
2+ oscillations can be classified into two types: the first type is simple Ca
2+ oscillations, and the second is bursting. Simple Ca
2+ oscillations have been explored extensively. In the experiments of Borghans et al. (1997), some more complex forms of Ca
2+ oscillations have been found such as bursts and chaos. When the system enters chaos, the Ca
2+ 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 Ca
2+ oscillations and performed mathematical analysis [
19,
20]. Chay also proposed a model for Ca
2+ oscillations in excitable cells [
21], and Shen and Larter (1995) studied complex Ca
2+ oscillations and chaotic phenomena in non-excitable cells [
22,
23]. In these models, it is almost certain that the burst results from changes in IP
3 yields because IP
3 can stimulate Ca
2+ channels and release a large number of Ca
2+ into the cytosol or take up Ca
2+ from the cytosol, thereby participating in Ca
2+ regulation. This is one of the most straightforward explanations for the bursting of intracellular Ca
2+. One mechanism proposed by Borghans et al. to explain complex Ca
2+ oscillations is that bursting oscillations are associated with Ca
2+-releasing channel activities in the endoplasmic reticulum. A further mechanism investigated in the same study is that there is not just one intracellular Ca
2+ store. Both IP
3-sensitive and non-sensitive Ca
2+ stores have been taken into account. In all models, the endoplasmic reticulum is the primary Ca
2+ store for intracellular Ca
2+.
Marhl et al. suggested the significance of free Ca
2+ concentrations in mitochondria for Ca
2+ oscillations and proposed a likely explanation for Ca
2+ oscillations in cells [
24]. In this study, the model focuses on Ca
2+ exchange in the middle of the cytosol, endoplasmic reticulum, and mitochondria. The Ca
2+ considered in the model includes cytosol, endoplasmic reticulum, and mitochondrial Ca
2+ [
25,
26] and binds Ca
2+ in the cytosol. These models well simulate enriched Ca
2+ 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 Ca
2+ oscillations theoretically. Therefore, this study takes Marhl et al.’s Ca
2+ 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].
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:
One can immediately see that the equilibrium point of system (2) meets the following equation:
Then,
can be calculated to obtain
Consequently, substituting the expression for
y and
z obtains
Therefore, we have the following equation:
where
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 =
x −
x0,
y1 =
y −
y0,
z1 =
z −
z0, then we can obtain the following representations [
30]:
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:
The linearized system corresponding to the equilibrium (0, 0, 0) of system (5) is
where
There is also a simple way to obtain the following eigen equation of the Jacobian matrix:
where
The following formula is obtained by calculating:
where
The Hurwitz matrix is obtained as follows.
If the determinants of the Hurwitz matrix are all positive, then it can be confirmed that the system is stable.
The specific expression of the Hurwitz criteria is
According to this inequality and MATLAB, the values of
r can be obtained:
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 d
r1/d
t = 0 into system (5), we assume that when
r1 =
r −
r0 when
r =
r0, system (5) can be sorted out as follows [
30]:
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.4880
i,
η3 = −0.4880
i,
η4 = 0. The eigenvector is
System (6) can vary in shape such as the following:
and
where
After the calculation, we can obtain the following equations:
Due to the center manifold theory, it can be conclusively demonstrated that system (6) has a central manifold, which can be expressed as follows:
Let
h1(
v,
w,
s) =
a1v2 +
b1w2 +
c1s2 +
d1vw +
e1vs +
f1ws +
…, the center manifold of system (6) can be expressed as follows:
Therefore, the high-order partial derivatives can be applied to obtain the values of
a1 to
f1. The equation is given below:
From the center manifold theory, one can see that
After downscaling system (6), it will be confined to a two-dimensional system as follows:
where
Hence, it is easy to verify that
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.1872
i,
η3 = −4.1872
i,
η4 = 0, and the eigenvector is
System (6) can vary in shape such as the following:
and
where
After calculation, we can obtain the following equations:
Let
h2 (
v,
w,
s) =
a2v2 +
b2w2 +
c2s2 +
d2vw +
e2vs +
f2ws + …, the center manifold of system (6) can be expressed as follows:
Therefore, the high-order partial derivatives can be applied to obtain the values of
a2 to
f2. The equation is given below:
From the center manifold theory, one can know that
After downscaling system (5), it will be confined to a two-dimensional system as follows:
where
Hence, it is easy to verify that
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 HB
1 and HB
2, 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 Ca
2+ 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
=
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
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.