Study on cell apoptosis based on bifurcation and sensitivity analysis

1.   Introduction

Apoptosis is a natural physiological process of programmed cell death. It helps an organism to eliminate cells which are no longer needed or irreparably damaged. Apoptosis plays a critical role in embryonic development and adult tissue homeostasis. Deregulation of apoptosis is an important aspect of cancer pathogenesis ^[1], and has been widely recognized as a hallmark of most types of cancer ^[2]. In general, apoptosis is strictly regulated by the body to maintain the stability of tissues and organs throughout the life cycle. However, when the regulation of apoptosis is out of balance, excessive proliferation or excessive apoptosis may occur, thus leading to the occurrence of related diseases such as ovarian cancer and endometriosis ^[3,4]. Therefore, deeper understanding the mechanisms of apoptosis contributes to the study of the pathogenesis of related diseases and plays an important role in the prevention and treatment of them.

There are three major signal transduction pathways in cell apoptosis: the death receptor, the mitochondrial, and the endoplasmic reticulum pathways, of which the former belongs to the outer pathway and the latter two belong to the intrinsic pathways. Here we only study the mechanism of apoptosis related to the intrinsic pathway. Caspases are a class of cysteine-aspartyl proteases that are synthesized as inactive precursor enzymes or proenzymes. These proteases normally lie dormant in healthy cells. In response to diverse stimuli they become activated, by either proteolytic cleavage or by recruitment into large complexes. Once activated, caspases cleave their substrates typically after conserved aspartate residues and are responsible for most of biochemical and morphological features of apoptosis. The p53 tumor suppressor gene, which acts as a key conduit in cells between DNA damage detection and apoptosis initiation, is frequently inactivated by mutations in cancer. Inactivation of p53 enables cancer cells not only to bypass intrinsic apoptotic responses to their underlying genomic aberrations, but also to block apoptosis initiation in response to DNA damage they incur upon exposure to conventional therapies ^[5]. Therefore, agents that can restore apoptosis in cancer cells hold promise for therapy and have been the focus of many preclinical drug discovery studies. Developmental cues or severe cell stress, such as DNA or cytoskeletal damage, stimulate the intrinsic pathway by inducing transcription or post-translational activation of BH3-only proapoptotic BCL-2 family proteins, which, in turn, drive activation of other family members such as BAX and BAK. In healthy cells, antiapoptotic family members such as BCL-2, BCL-XL, or myeloid cell leukemia sequence 1 (MCL-1) bind to proapoptotic members such as BAX and BAK to prevent their activation. Upon stimulation, BH3-only proteins sequester the antiapoptotic proteins away, allowing BAX and BAK to act. BAX and BAK engage the mitochondria to release factors that promote caspase activation in the cytoplasm. One such factor is cytochrome c, which cooperates with another protein called Apaf-1 to activate caspase 9. This caspase, in turn, activates the effector caspases 3, 6, and 7, which carry out apoptosis ^[6].

A mathematical model of the intrinsic pathway that addresses specifically the dynamical basis of the three physiological features of apoptosis was presented ^[7]. And they subdivided the intrinsic pathway into three modules (initiator, amplifier, and executioner modules). Based on prior work, we identify important regulatory processes in apoptosis through single-parameter sensitivity analysis and bifurcation theory. Parameter sensitivity analysis has been used by many researchers ^[8] and has been supported by some experiments. The work of ^[9] is encouraging in this regard. The authors studied the mitochondrial outer membrane permeabilization network that controls apoptosis, and found that the oncogenic mutation of each protein functional domain is strongly correlated with the parameter sensitivity of the bifurcation point. In ^[10], the authors determined the location of the bifurcation point rather than the concentration of some specified proteins by identifying parameters that have a large influence on the bifurcation points. The bifurcation-based sensitivity is a better measure for determining the function of the bio-network and plays an important role in the field of genetic mutation ^[8,9,10]. Fortunately, we found that modules or feedback loops but not specific parameters determine the network functions although similar bifurcation-based sensitivity analysis is performed.

2.   Model and methods

In previous studies, there were many models of apoptotic pathways. In this paper, we use the model described in Figures 1 and 2 to analyze. For more detailed information of the model, please refer to literature ^[7]. The model of the intrinsic pathway describes the dynamics of the three characteristic physiological features of apoptosis. In Figure 1, the intrinsic path is subdivided into three modules (initiators, amplifiers, and actuator modules). It also shows how the positive feedback in the initiator module is responsible for signal thresholds and time delays, and how the actuator module ensures commitment to cell death. The initiator module describes BAX activation by BH3 proteins. The amplifier module describes how BAXm causes the release of CytoC and SMAC. The executioner module describes how CytoC and SMAC cooperate to activate caspase. The schematic draft of the network is shown in Figure 2. Here, we mainly analyze three modules: initiator module, and executioner modules and the how p53 module governs cell fates after DNA damage.

Figure 1.  Wiring diagrams for a three-module decomposition of the intrinsic apoptotic pathway. Initiator module in the green box, amplifier module in the red box, and executioner module in the blue box.

DownLoad: Full-Size Img PowerPoint

Figure 2.  Wiring diagrams for processes regulating p53 and those regulated by p53.

DownLoad: Full-Size Img PowerPoint

We perform bifurcation and sensitivity analysis to investigate dynamics of the system. Because the BAXm concentration, i.e., [BAXm], is an important measure of the initiator module (startup the executioner module), we use stress as the parameter and [BAXm] as the output, then we can get bifurcation diagram, as shown in Figure 3(a). As the stress increases from low to high, [BAXm] stays at low levels until the stress exceeds a threshold SN1, at which a saddle-node bifurcation occurs and the [BAXm] increases suddenly. Such a bifurcation corresponds to the conversion of the initiator module to the executioner module. Similarly, the other saddle-node bifurcation occurs at SN2 if the stress is decreased from high to low.

Figure 3.  Bifurcation diagrams with stress as a control parameter. (a) Two solid lines indicate stable states and dashed line shows unstable states. SN1 is one transition point when the stress increases from low to high and SN2 is the other transition point as stress decreases instead. (b) The change of the locations of the bifurcation points SN1 with increasing the parameter kd3. Black, control; blue, 0.5-fold increase; green, 1-fold increase.

DownLoad: Full-Size Img PowerPoint

At the high stress levels, the system starts producing lots of BH3 protein. At first, most [BH3] molecules bind with [BCL2], so they cannot activate [BAX]. Eventually, BH3 titrates out the BCL2 pool, and excess BH3 molecules are now free to convert [BAX] to [BAXm]. BAXm displaces additional BH3 from the BH3/BCL2 dimer pool, thereby speeding up its own production from inactive [BAX]. This positive feedback causes the [BAXm] level to rise quickly, which is sufficient to activate the amplifier module ^[7]. When the stress decreases from high to low, there is another saddle-node bifurcation point SN2. When the stress lies between SN1 and SN2, the module is bistable. The bistable states represent two different functional states. The higher level of [BAXm] stands for the initiator module and the lower level indicates the amplifier module.

It is known that bifurcations exhibit qualitative changes in biological systems and are found in many biological functions ^[11,12]. In this system, the saddle-node bifurcation provides a mechanism of switching between modules, i.e., a switch mechanism that ultimately indicates life-to-death. Bifurcation point location has an important effect on apoptosis. However, some parameter changes may lead to large changes in positions of bifurcation points, resulting in abnormal system functions. Therefore, it is an important issue to determine which parameters have more significant influence on the positions of bifurcation points. In this paper, we determine the important parameters and therefore important regulatory processes in fate decisions through single parameter sensitivity analysis of the bifurcation point. For instance, the bifurcation point SN1 is clearly shifted to the right by increasing the parameter kd3, i.e., the BH3 degradation rate, by 0.5 fold, as shown in Figure 3(b). Increasing kd3 to 1 fold of the basal value, SN1 continues to move to the right, and high [BAXm] can be achieved at higher stress, which means that if BH3 accelerates self-degeneration due to mutations in certain genes, even if the stress is very high, apoptosis will be difficult to carry out. Sufficient stress has no way to induce apoptosis, providing the possibility of tumorigenesis. Similarly, we performed bifurcation analysis of the executor and the p53 modules. This conclusion is similar to the initiator module, as shown in Figures 4 and 5.

Figure 4.  Bifurcation diagrams with [CytoC] as a control parameter. (a) Two solid lines indicate stable states and dashed line shows unstable states. SN1 is one transition point when [CytoC] increases from low to high and SN2 is the other transition point as [CytoC] decreases instead. (b) The change of the locations of the bifurcation points SN1 with increasing the parameter kdac9. Black, control; blue, 0.5-fold increase; green, 1-fold increase.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Bifurcation diagrams with [p53killer] as a control parameter. (a) Two solid lines indicate stable states and dashed line shows unstable states. SN1 is one transition point when [p53killer] increases from low to high and SN2 is the other transition point as [p53killer] decreases instead. (b) The change of the locations of the bifurcation points SN1 with increasing the parameter k2s3dr. Black, control; blue, 0.5-fold increase; green, 1-fold increase.

DownLoad: Full-Size Img PowerPoint

3.   Simulation calculation

As far as we know, the influence of parameter variations on dynamics is uneven in many biological systems ^[13,14]. And the effects of parameter variations on the positions of the bifurcation points have not been well investigated. There are many parameters in the model that can affect the positions of bifurcation points. In order to determine which parameters have more significant influences on the locations of bifurcation points in the model, we perform bifurcation analysis and record the percentage change in bifurcation points by increasing or decreasing each parameter by 10% off its base value. Thus, single parameter sensitivity of bifurcation points are derived and further the sensitivity spectrum of all parameters can be achieved.

3.1.   Initiator module

BCL2 is an anti-apoptotic gene that prevents apoptosis caused by various stimulating factors and promotes cell survival. When a stress signal is applied to the initiator module, the BH3 level increases. At first, these BH3 molecules are taken out of service by binding to the inhibitory protein, BCL2. If the stress signal is sustained, enough BH3 will accumulate to saturate the pool of anti-apoptotic [BCL2] molecules. Free and active BH3 induces conformational changes in BAX which facilitates its insertion into the mitochondrial outer membrane. The membrane-localized form of BAX, i.e., BAXm, binds strongly to BCL2, releasing additional free BH3 to promote the conversion of BAX to BAXm ^[7,15,16,17].

The initiator module can be described by four ordinary differential equations (ODEs), once we assume constant total concentrations of BAX and BCL2 proteins. All numerical simulations and bifurcation analysis are performed based on these equations ^[7]. The parameter sensitivity analysis is shown in Figure 6. All the saddle-node bifurcation points are distributed under the standard parameter set and the perturbed parameter set with each parameter increased or decreased by 10% off its base value, as shown in Figure 6(a). Notably, one sees that the all SN2 are always positive, which means the switch is robustly reversible under the given perturbations.

Figure 6.  Parameter sensitivity analysis of the bifurcation points. (a) The level of [BAXm] at the two saddle-node bifurcation points SN1-2 when each of 10 parameters in the feedback loops is increased or decreased by 10% with respect to the standard parameter set. (b) Percentage of changes in the levels of [BAXm] at SN1-2 compared to that under the standard parameter set.

DownLoad: Full-Size Img PowerPoint

We aim to identify specific proteins through sensitive analysis, and then find out important feedback loops or modules in apoptosis. Due to very few parameters in the initiator module, we analyze the role of individual proteins. We then further calculate the percentage of changes in the two thresholds, as shown in Figure 6(b). Since the relative influence of the parameters on the bifurcation is basically unchanged, we distinguish the strong and weak influences according to the overall change of the parameters. Here we choose 10% as the standard. Most of the parameters in initiator module (approximately 70%) are insensitive, i.e., perturbations to these parameters result in small changes in the bifurcation. But some (approximately 30%) are sensitive and perturbation to them has a strong influence on the bifurcation, e.g., kd3 (decomposition rate of BH3), ks32 (stress stimulates the speed of BH3), and kb (the conversion rate from BAX to BAXm). Increasing parameter kd3 causes the SN1 point to move to the right, which requires a higher stress level to achieve a high level of BAXm than the basal rate. In other words, more stress is needed to start the amplifier module.

As the free BAX increases, the oligomers form and open channels in the outer membrane of the granules. CytoC and SMAC leave the mitochondria and enter the cytoplasm through these channels, where they cooperate to activate the execution daughter caspase 3 and then eventually degenerate. Given the parameter settings in the model, the amplifier module only transfers the signal from the mitochondrial BAX to the caspase in the cytoplasm ^[7,18,19]. Therefore, we do not conduct in-depth research on the amplifier module.

3.2.   Executioner module

The executioner module describes caspase activation by [CytoC] and [SMAC]. In order to further study the executive module, we use parameter sensitivity analysis again. The executioner module uses caspase 3 activity as an indicator of cell death ^[7]. Therefore we use CytoC as the input variable and caspase 3 activity as the output. The results of parameter sensitivity are shown in Figure 7. To map cancer-associated mutations to network dynamics through protein-protein interaction dynamics, we have narrowed our research to the parameters of protein interactions. We separated the parameters into the following three clusters according to their impact on the bifurcation point locations: most sensitive (red and blue), sensitive (green), and insensitive (black). In other words, the thresholds for the switch are most sensitive to the parameters in both the aC9 autocrine expression and C9/XIAP feedback loop, whereas the thresholds for the switch are sensitive to the parameters in the C3/XIAP feedback loop. Caspase 9 plays a more important role in the intrinsic pathway of apoptosis. A large number of experimental articles have also shown that apoptosis can be promoted by inducing caspase 9. For example, BMP4 inhibits caspase 9 via the PI3K/AKT signaling pathway, thereby inhibiting apoptosis of GCs ^[20]. Caspase 9 and AIF play important roles in the elongation-induced apoptosis of myoblasts, and AIF participates in this process in a manner independent of caspase-9 ^[21].

Figure 7.  Parameter sensitivity analysis. 10% decrease or increase of each parameter in the executioner module causes a percentage change in the positions of the bifurcation points SN1 and SN2. Most sensitive (red and blue): Red, aC9 autocrine expression, blue, C9/XIAP feedback loop; sensitive (green), C3/XIAP feedback loop; and insensitive (black), the rest proteins.

DownLoad: Full-Size Img PowerPoint

3.3.   The p53 module governs dynamical cell fates after DNA damage

After DNA damage, the p53 module dominates the cell fates. We want to explore which proteins play a decisive role in the process. From the study of the initiator module, we know that the level of [BAXm] is an important measure of the initiator module related to the executioner module, we use DNA damage as the input variable and [BAXm] as the output, then we can get the bifurcations diagram. Similar to previous analysis, we perform the bifurcation analysis and record the percentage change in bifurcation locations by increasing each parameter in the model by 10%. The sensitivity spectrum of all parameters is achieved in Figure 8. Similar results can be obtained for SN2 and the case of decreasing each parameter 10% off its base value.

Figure 8.  Parameter sensitivity analysis. 10% increase of each parameter in the p53 module causes a percentage change in the position of the bifurcation points SN1. Sensitive (red, yellow and green): Red, BH3AC/BH3DR feedback loop, yellow, MDM2/ARF/ASPP feedback loop, green, p53helper/p53killer feedback loop; and insensitive (blue), the rest proteins.

DownLoad: Full-Size Img PowerPoint

In this module we choose 4% as the standard. Of course, other criteria can be selected, but the relative influence of the parameters on the bifurcation point is substantially unchanged. From Figure 8, most of the parameters (approximately 82%) are insensitive: perturbations in these parameters result in small changes in the bifurcation. But some (approximately 18%) are sensitive: perturbing these parameters individually has a strong influence on the bifurcation. For sensitive parameters, we observed that the proteins represented by the sensitive parameters mainly focus on three feedback loops: Red, BH3AC/BH3DR feedback loop; yellow, MDM2/ARF/ASPP feedback loop; green, p53helper/p53killer feedback loop. In other words, these three feedback loops play a key role in the cell fate of p53-regulated after DNA damage. It has been experimentally proved that selection of a suitable combination of BH3 mimetic compounds can induce apoptosis of basophils and mast cells in mice ^[22]. UCHL1 activates the MDM2-ARF-p53 signaling pathway, which has two opposite enzymatic activities in the ubiquitin pathway, further leading to tumor cell apoptosis ^[23]. The results of the above experiments also indicate that BH3 plays an important role in inducing apoptosis, consistent with the conclusion of BH3AC/BH3DR feedback loop in this paper.

4.   Discussion

To investigate the connection between core functions of biological networks and oncogenesis, we conducted sensitivity analysis of all the parameters in the apoptosis model. We found the sensitive parameters that have greater influence on the bifurcation locations. The positions of bifurcation point represent threshold of the initiator module triggering the amplifier module. Only when the stress signal exceeds a specific value, the system will induce the transition from the initiator module to the amplifier module. If BH3 accelerates self-degeneration due to mutations in certain genes, even if the stress is very high, apoptosis will be difficult to carry out. Sufficient stress has no way to induce apoptosis, providing the possibility of tumorigenesis.

In this paper, we determine the important regulatory processes through single parameter sensitivity analysis of the bifurcation points. In executioner module, we found that aC9 autocrine expression and C9/XIAP feedback loop is more sensitive than the C3/XIAP feedback loop. In addition, we found that BH3AC/BH3DR feedback loop, MDM2/ARF/ASPP feedback loop, and p53helper/p53killer feedback play more important roles in cell fates after DNA damage. These results mean that modules or feedback loops but not detailed parameter values determine the network functions. Such results are partially supported by experimental findings. Our research provides a rational basis for the design of clinical trials of anticancer drugs.

Acknowledgments

The authors thank the anonymous reviewer for his valuable comments and suggestions which are helpful to improve this paper. This work was supported by the National Science Foundation of Shanghai (Grant No. 17ZR1410800).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

Initiator module

We use $x_{1}$, $x_{2}$, $x_{3}$, and $x_{4}$ to denote [BAXm], [BAXm/BCL2], [BH3], and [BH3/BCL2], so we can obtain the dynamical equations as follows

where $xf = xt-x_{1}$, $mf = x_{1}-x_{2}$, $bh3 = x_{3}-x_{4}$, and $bcl2 = rt-x_{2}-x_{4}$.

The parameter values are chosen as

$stress = 0.1$, $rt = 80$, $xt = 100$, $kf1 = 1$

$kf2 = 3$, $kb = 2$, $kam2 = 90$, $kdm2 = 0.05$

$ka32 = 10$, $kd32 = 0.01$, $ks31 = 0.1$, $ks32 = 0.6$, and $kd3 = 0.01$.

The initial values are

$x_{1} = 33.4$, $x_{2} = 33.4$, $x_{3} = 16$, and $x_{4} = 16$.

Executioner module

We use $x_{5}$, $x_{6}$, $x_{7}$, $x_{8}$, $x_{9}$, and $x_{10}$ to denote [C9], [aC9], [C3], [XIAP/C9], [XIAP/C3], and [SMAC/XIAP], respectively. Then we get the dynamics as follows

where $cf = smaccytot-x_{10}$, $pf = pt-x_{10}-x_{9}-x_{8}$.

Parameter values are

$cytoc = 0.1$, $smaccytot = 2$, $k1ac9 = 0.001$, $k2ac9 = 0.5$

$ksc9 = 0.001$, $c9t = 1$, $kdc9 = 0.002$, $kdac9 = 0.003$

$k1ac3 = 0.001$, $k2ac3 = 0.02$, $k3ac3 = 0.5$, $c3t = 1$

$kdc3 = 0.002$, $kas3x = 0.2$, $kds3x = 0.5$, $kd3x = 0.1$

$kas9x = 0.1$, $kds9x = 0.6$, $kd9x = 0.2$, $kassx = 2$

$kdssx = 0.01$, $kdsx = 0.007$, $pt = 6$, and $n = 2$.

The initial values are

$x_{5} = 0$, $x_{6} = 0$, $x_{7} = 0$, $x_{8} = 0$

$x_{9} = 0$, and $x_{10} = 0$.

P53 module

We use $x_{11}$, $x_{12}$, $x_{13}$, $x_{14}$, $x_{15}$, $x_{16}$, $x_{17}$, $x_{18}$, $x_{19}$, $x_{20}$, $x_{21}$, $x_{22}$, $x_{23}$, $x_{24}$, $x_{25}$, $x_{26}$, $x_{27}$, $x_{28}$, $x_{29}$, $x_{30}$, and $x_{31}$ to denote [DYRK2], [p53helper], [p53killer], [ASPP], [MDM2], [ARF], [ARF/MDM2], [CycE], [E2F1], [RBp], [p21], [p21/CycE], [BAX], [BAXm], [BCL2], [BH3AC], [BH3DR], [BAXm/BCL2], [BH3AC/BCL2], [BH3DR/BCL2], and [p73]. The dynamics of the P53 model can be expressed by the following equations

where

$p53t = x_{12}+x_{13}$, $rbe2f1 = e2f1t-x_{19}$, $rb = rbt-x_{20}-rbe2f1$

$ks3ac = k1s3ac+k2s3ac\times (1+x_{14})\times x_{13}+k3s3ac\times x_{21}+k4s3ac\times (1+x_{14})\times x_{12}$

$ks3dr = k1s3dr+k2s3dr\times (1+x_{14})\times x_{13}+k3s3dr\times x_{21}+k4s3dr\times (1+x_{14})\times x_{12}$.

Parameter values are

$dna = 0.1$, $kin = 0.00015$, $kout = 0.0003$, $dyrt = 2$

$ksp53 = 0.5$, $k1dp53 = 0.1$, $kpp53 = 1$, $jpp53 = 0.1$

$kdpp53 = 0.5$, $jdpp53 = 0.1$, $pp = 1$, $k1saspp = 0.01$

$k2saspp = 0.05$, $kdaspp = 0.1$, $k1smdm2 = 0.02$,

$k2smdm2 = 0.3$, $k1dmdm2 = 0.1$, $k2dmdm2 = 0.1$, $k1sarf = 0.01$

$k2sarf = 0.3$, $kdarf = 0.1$, $kasam = 10$, $kdsam = 2$

$k1se = 0.01$, $k2se = 0.5$, $kde = 0.12$

$kasre = 5$, $kdsre = 1$, $kprb = 1$, $jprb = 0.1$

$kdprb = 0.5$, $jdprb = 0.1$, $rbt = 2$, $e2f1t = 1$

$k1sp21 = 0.03$, $k2sp21 = 0.3$, $k3sp21 = 0.01$, $kdp21 = 0.2$

$kasp21e = 10$, $kdsp21e = 1$, $ksbax = 0.97$, $kdbax = 0.01$

$kf1 = 0.1$, $kf2 = 3.4$, $kb1 = 4.2$, $kb2 = 4.2$

$kdbaxm = 0.01$, $kasml = 87$, $kdsml = 0.05$, $ksbcl2 = 0.57$,

$kdbcl2 = 0.01$, $kas3acl = 11.14$, $kds3acl = 0.01$, $kas3drl = 11.14$

$kds3drl = 0.01$, $kd3ac = 0.01$, $kd3dr = 0.01$, $kdml = 0.01$

$kd3acl = 0.01$, $kd3drl = 0.01$, $k1sp73 = 0.01$, $k2sp73 = 0.01$

$kdp73 = 0.5$, $k1s3ac = 0.01$, $k2s3ac = 0.05$, $k3s3ac = 0.1$

$k4s3ac = 0.02$, $k1s3dr = 0.01$, $k2s3dr = 0.1$, $k3s3dr = 0.1$, and $k4s3dr = 0.02$.

The initial values are

$x_{11} = 0$, $x_{12} = 0$, $x_{13} = 0.78$, $x_{14} = 0.55$

$x_{15} = 2.53$, $x_{16} = 2.79$, $x_{17} = 1.99$, $x_{18} = 3.83$

$x_{19} = 0.9$, $x_{20} = 1.88$, $x_{21} = 1.32$, $x_{22} = 1.25$

$x_{23} = 95$, $x_{24} = 0$, $x_{25} = 47.49$, $x_{26} = 0$

$x_{27} = 0$, $x_{28} = 2.26$, $x_{29} = 3.78$, $x_{30} = 3.78$, and $x_{31} = 0.04$.