Parameter sensitivity analysis for a stochastic model of mitochondrial apoptosis pathway

Understanding how gene alterations induce oncogenesis plays an important role in cancer research and may be instructive for cancer prevention and treatment. We conducted a parameter sensitivity analysis to the mitochondrial apoptosis model. Both a nonlinear bifurcation analysis of the deterministic dynamics and energy barrier analysis of the corresponding stochastic models were performed. We found that the parameter sensitivity ranking according to the change of the bifurcation-point locations in deterministic models and the change of the barrier heights from a living to death state of the cell in stochastic models are highly correlated. For the model we considered, in combination with previous knowledge that the parameters significantly affecting the system’s bifurcation point are strongly associated with frequently mutated oncogenic genes, we conclude that the energy barrier height can be used as indicator of oncogenesis as well as bifurcation point. We provide a possible mechanism that may help elucidate the logic of cancer initiation from the view of stochastic dynamics and energy landscape. And we show the equivalence of energy barrier height and bifurcation-point location in determining the parameter sensitivity spectrum for the first time.


Introduction
It has been a consensus that cancer is a complex disease resulting from genomic alterations [1][2]. However, the mechanism of mutation-induced oncogenesis is not fully understood. In recent years, there are several interpretations have been proposed. Huang et al. [3] put forward the concept of 'cancer attractors' and thought that genetic aberrations would increase the region of cancer attractors or reduce the threshold of entering cancer basins. By constructing global potential landscape with simplified cellular networks, Li et al. [4] showed that gene mutations could trigger cell state transitions from normal to cancer states. Stites et al. [5] focused on the Ras functional module and demonstrated that gene mutations could alter the biochemical properties and further cause pathological changes. We presented a framework for mapping gene mutations to tumorigenesis [6][7][8] bifurcation behavior and investigated mutation enrichments in DNA damage induced apoptosis, Rb-related cell cycle and specific mitochondrial apoptosis pathways. We established that the location of the bifurcation point could be a threshold for cell cancerous changes [6][7], and subsequently verified the functional role of the bifurcation point with more precise relation between protein functional domain mutations and parameters [8].
The nonlinear dynamic bifurcation research in cancer-related networks is based solely on continuous and deterministic formalism. In reality, the biochemical reactions in cells are stochastic and discrete in nature. In particularly, when the number of reactant molecules is small, random fluctuations might be important, and cannot be negligible. What's more, the stochastic effect will become apparent when the system is close to the bifurcation point, where the quasipotential energy barrier becomes relatively low and the transitions become easier. This case happens even when the number of molecules is large. Thus, it is desirable to investigate the quantitative feature of the biological networks with stochastic dynamics and figure out what's the indicators of oncogenesis in the corresponding stochastic systems, which is the aim of this paper.
Here, we used mitochondrial apoptosis dynamics to perform our study. From the perspective of the energy landscape (Fig 1B), the living and death states of the cell correspond to two steady states in the landscape, and the process of apoptosis corresponds to a dynamic path in which the system at steady state climbs an energy barrier and transits to another steady state. The height of this barrier indicates the likelihood of apoptosis. The higher the barrier is, the more difficult it for the cell to undergo apoptosis. Thus, it is more likely for the cell to become cancerous. By constructing the energy landscape of apoptosis dynamics, we can quantitatively study the potential of cancerization.
In this paper, we first determine the nonlinear bifurcation and parameter sensitivity analysis results for the deterministic dynamics of the mitochondrial apoptosis model. Then, we construct the corresponding stochastic model for the mitochondrial apoptosis pathway. By taking advantage of the quasi-potential landscape concept in probability theory [9,10] (see Sec. 2.3 and Sec. 1.7 in S1 File for a brief introduction), we compute the energy barrier height required for a cell to transit from a normal to a death state using the Geometric Minimum Action The living and death states correspond to the two local minima of the potential function in terms of the apoptosis effector -Bax 4 Method (GMAM) [9,10] and direct Stochastic Simulation Algorithm (SSA) [11]. We then study the changes in the barrier height in response to parameter variations and compare these results with those of the nonlinear bifurcation analysis. In addition, we perform Sobol sensitivity analysis which is a variance-based measurement and compare the results with the outcomes of the local analysis.
Through comparison and analysis, we find that the parameter sensitivity ranking according to the changes in bifurcation-point locations in deterministic models and the change in the barrier heights from living to death states of the cell in stochastic models are highly correlated. In particular, when the system is close to the bifurcation point, the parameter sensitivity patterns are almost identical. For the model that we considered here, our group previously showed that the parameters significantly affecting the system's bifurcation point are strongly associated with frequently mutated oncogenes. Taken together, these results suggest that both the energy barrier height and bifurcation point position are indices of tumor formation.

Models
We focus on the mitochondrial apoptosis pathway, which is the most commonly deregulated form of cell death in cancer [12][13][14] and the cross-link of the intrinsic and extrinsic apoptosis pathways [15][16][17]. The mechanism of this pathway can be briefly summarized as follows: Cas-pase8 is an initiator caspase which can be activated in response to death stimulation, which originates from DNA replication stress, unfolded protein response or other causes. Activated Caspase8 further activates pro-apoptosis activators, such as Bid and BIM. The accumulation of activators will give rise to the oligomerization of the effectors Bax and Bak which results in pores in the mitochondrial membrane. The pores then trigger mitochondrial outer membrane permeabilization (MOMP). Subsequently intermembrane space proteins such as cytochrome c and Smac are released to the cytoplasm and eventually lead to downstream cascades and cell death [18][19][20]. We take advantage of the mitochondrial apoptosis model of Zhao et al [8] which conformed to the biological facts we have known. The entire regulation network is shown in Fig 1A. Proteins with redundant functions are compressed into one representative reactant. For example, Bax and Bak are the major effectors of apoptosis, either Bax or Bak alone is sufficient to form oligomer and then lead to MOMP, therefore, we used Bax as the representative. The ODEs in this network can be found in SI. We denote it by for short. x ! represent the reactants. The dynamic behavior of the system is governed by the function bð x ! Þ, which is derived according to the law of mass action.
In the progress of mitochondrial apoptosis, MOMP acts as a defining event that irreversibly commits cells to death. Consequently, for simplicity, cell MOMP is considered a death state without accounting for post-MOMP regulation. The concentration of Bax 4 MOM is the output variable, and the Caspase8 concentration is the input variable. In this continuous-state deterministic dynamics, the system has two stable fixed points and one saddle point in the parameter region that we consider, which gives rise to a bistability character.

Stochastic model and simulation of mitochondrial apoptosis
Starting from the deterministic model described above, we now discuss the stochastic setup of the system. We model the mitochondrial apoptosis dynamics as a Markov jump process. We denote the state of the system by a vector X = (X 1 ,X 2 ,. . .,X 20 ), whose i-th component represents the number of molecules for the i-th species in ODEs (1). Each reaction in Fig 1 can be described by a propensity function a j (X) and a state change vector ν j . The system volume size V plays an important role in describing the propensity function. There are three types of reactions in this system. The first type is when a certain molecule is added to the system from outside with a constant rate, such as 0! k 5 tBid. In this case, the propensity function is a(X) = Vk 5 . The second type is the involvement of only one species in the reaction, such as Bid ! Cas8 k 1 tBid. The propensity function of this reaction is a(X) = cas8k 1 X 1 . For the last type of reaction, two types of molecules (including two of the same molecule) are involved, such as The propensity function of this reaction is a(X) = V −1 k 12 X 13 (X 13 − 1).
State change vector characterizes the change of system state if a reaction fires. For example, in the Bid! Cas8 k 1 tBid reaction, the state change vector is ν = (−1,1,0,. . .,0), which means that the Bid number will be decreased by 1 step reaction, and the tBid number will be increased by 1 by one step reaction. Once fired, the state of the system would be updated from X to X + ν. All other reactions can be described similarly. There are 53 total reactions in this system. Deterministic Eq (1) is the limit of the stochastic process X/V as V goes to infinity [21].

Energy landscape and large deviation theory
To quantitatively understand mitochondrial apoptosis, we construct the energy landscape for apoptosis dynamics. For a general deterministic dynamical system, there is typically no natural potential landscape. In recent studies, Li et al. [9][10] showed that the quasi-potential in large deviation theory [22][23] is considered as one mathematical realization of the energy landscape. The quasi-potential landscape is defined as follows. For any stable steady state x 0 of Eq (1), the local quasi-potential S(x;x 0 ) with respect to x 0 is defined as the minimization of an action function where φ(t) is a continuous path connecting x 0 and x [24]. The function L(x,y), which can be considered as the Lagrangian, is the Legendre transform of a Hamiltonian H(x,p). Namely, For the Markov jump process described in Sec.2.2, the Hamiltonian has the form while L(x,y) does not have a closed form, where the summation is taken over all reactions. Based on the local quasi-potentials constructed from different stable steady states, the global quasi-potential S(x) is a proper combination of all local versions. One may refer to [25] for technical details.
The local quasi-potential S(x;x 0 ) characterizes the difficulty of transition from x 0 to x. Indeed, if we denote the FPT from x 0 to x by τ, the large deviation theory tells us [21] that when ! 0, where Eτ is the mean first passage time. Thus, we can use the barrier height of the quasi-potential landscape to study the oncogenic potential quantitatively. The higher this barrier is, the more difficult it is for the cell to undergo apoptosis and the easier it is for a normal cell to change into a cancer cell. Let us denote the live state as the lower Bax 4 M OM level by x low , and the death state as the higher Bax 4 M OM level by x high . The minimum action S for the transition from x low to x high can be computed efficiently using the GMAM algorithm [11,24]. We then utilize the obtained barrier heights to study the sensitivity of the oncogenic potential with respect to different parameter choices. The code for computing the barrier height with GMAM algorithm can be downloaded from the website (http://dsec.pku.edu.cn/~tieli/code/ Apoptosis-Code.zip).
Notably, if we select x 0 = x low and x = x high , then the mean time Eτ is the mean transition time from x low to x high . This time may be slightly different from the MFPT defined in Sec. 2.2. If we select a threshold less than x high , i.e. x high 2 fBax MOM 4 > thresholdg, then Eτ always overestimates MFPT in SSA simulation. Fortunately, when the volume size V is sufficiently large, the transition from x low to x high can be rare. Most of the transition time is spent climbing up the energy barrier. Once the system climbs over the barrier, it rapidly relaxes to x high . In such a case, both MFPT in the SSA simulation and Eτ are approximately equal to the time required to climb the energy barrier. Therefore, in practice, MFPT % Eτ. With this approximation, MFPT also quantifies the energy barrier height, depicted as Arrhenius' law. Thus, SSA can also be used to estimate the energy barrier height.

The bifurcation and parameter sensitivity analysis of the deterministic model
We take the concentration of Bax 4 M OM as the output and the Caspase8 concentration as the input to develop a bifurcation diagram (as shown in Fig 2) of the deterministic model. Under the biological meaningful conditions, the system has bistable characteristics. The low stable steady state branch indicates the cell living state and the high stable steady state represents the death state. The location of the saddle-node (SN) bifurcation point is Cas8 = 23.59nM which is in the concentration range of Caspase8 to activate downstream effector caspases [26]. Here, we use two different methods, local single-parameter analysis and variance-based Sobol analysis [27][28], to obtain a parameter sensitivity spectrum. Local sensitivity analysis is a classical method that is used to examine the impact of small perturbations of parameters on model outputs, and Sobol sensitivity analysis is a global sensitivity analysis method that is used to study how large variations of parameters affect the outputs [29]. The local sensitivity analysis is typically efficient in computer time but may be inadequate for nonlinear models compared with the global sensitivity analysis methods [30]. Fig 3A presents the parameter sensitivity spectrum calculated using single-parameter sensitivity analysis with a 2% change of each parameter. Here, we mainly focus on the parameter changes that can increase cell oncogenic potential. A right shift of the bifurcation point suggests an increasing oncogenic potential [12], that is, moving the bifurcation point from left to right will make it difficult for cells' death and be advantageous to the carcinogenic process. Therefore, the direction of parameter perturbation (increasing or decreasing) that leads to a right shift of the critical point is chosen. We also alter the perturbation size of the parameters (Fig A, B and C in S1 File) and find similar patterns of sensitivity as in Fig 3A. The variance-based Sobol sensitivity analysis makes no assumptions about the relations between models' inputs and outputs; it is a way to obtain global information on parameter sensitivity. In this calculation, we randomly generate parameters in a large parameter space using Sobol sequences and ten thousand sets of parameters that can render bistable behavior. Then, we use standard Sobol calculation method to get parameter sensitivity spectrum, as shown in Fig 3B. Here, we focus mainly on the total effect index, which quantifies the overall effects of a parameter. The principle theory of the Sobol sensitivity analysis are shown in SI. Spearman's correlation of two different sensitivity analysis methods is 0.981, with a significance level less than 0.001. and 90% of mean value of high state that characterize the change of low to high steady states are also used. There is almost no difference in sensitivity spectra using different thresholds (Fig D and E in S1 File). Fig 4A shows   Considering a volume size of V = 100 μm 3 Cas8 = 23 nM, the other parameters are the same as those in the deterministic model. Other volume sizes can also be chosen, the spectrum didn't change much (Fig F  and G in S1 File). We see that apoptosis also be activated without triggering the bifurcation with high levels of Cas8, as an important property obtained from our stochastic model. Ten thousand dynamic change trajectories with the same condition are obtained for statistical purposes. We calculate the mean first passage-time (MFPT) in response to a 2% change of each parameter. From this perspective, a longer MFPT indicates an enlarged oncogenic potential. Therefore, the direction of each parameter perturbation (increase or decrease) is chosen to make the MFPT longer. The parameter sensitivity ranked according to the influence that each parameter has on MFPT is shown in Fig 4B. The red and blue histograms illustrate the changes of MFPT in response to a 2% increase and decrease of the parameters, respectively.

Parameter sensitivity of energy barrier height
Based on the GMAM method described in Sec. 2.3, we can compute the energy barrier of the quasi-potential energy landscape using the same parameter setup. Fig 5A shows the results of the quasi-potential energy landscape computation in terms of log Bax 4 MOM when Cas8 = 4.9nM. The living and death states correspond to two local minima of the potential function. The result that the local landscape of death state is much deeper than that of living state is consistent with biological facts. The quasi-potential energy landscape in terms of Bax 4 MOM when Cas8 = 22 nM is also shown (Fig H in S1 File). In Fig 5B, the energy barrier changes with Cas-pase8 input. Indeed, when Caspase8 increases to the bifurcation point, the energy needed to exit the normal state decreases until reaching 0. When Cas8 is far from bifurcation point, the transition from normal state to death state is rare. Direct SSA simulation cannot generate a reliable statistic of MFPT in a reasonable time. Luckily, the energy barrier can still be calculated. As shown in Sec. 3.1 and Sec. 3.2, we change each parameter by 2% and observe the changes in the barrier height. Fig 5C and 5D show the results. The rank of sensitivity when Cas8 is near the bifurcation point ( Fig 5C) coincides quite with the SSA simulation and ODEs analysis. The Spearman's correlation between sensitivity of barrier height and MFPT in SSA is 0.973, with a significance level less than 0.001. Spearman's correlation between the sensitivity of the barrier height and bifurcation point in ODEs is 0.981, with a significance level less than 0.001. Fig 5D shows the sensitivity of the barrier height change in response to the change in each parameter when Cas8 = 4.9 nM. One may find that the barrier height is more than 50 times larger than that of case Cas8 = 22 nM (shown in Fig  5B). The rank of sensitivity when Cas8 is relatively far from the bifurcation point shows slightly different from the results when Cas8 is close to the bifurcation point. This difference is mainly between parameters k5 (production rate of Bid) and kf12 (dissociation rate of membranebinding Bax). This distinction suggests that the relative importance of the two parameters is different in low Cas8 and high Cas8 cells. The barrier height of our model can reflect the potential of cancer initiation. Thus, in cells with high Cas8 concentrations, the oncogenic potential of Bax MOM dimer mutation is lower than that for Bid. Inversely, the oncogenic potential of Bax MOM dimer mutation is higher than for Bid in low Cas8 cells. Spearman's correlation between these two results shown in Fig 5C and 5D is 0.826, suggesting that the rank bifurcation-point location in response to the corresponding changes in the parameters. Details of the correspondence between parameters and representative interactions are shown in Supporting Information. Asterisk, association; deg, degradation rate; onMOM/offMOM, membrane translocation and membrane separation; pro, production rate; minus, dissociation. (b) Total effect index of each parameter calculated by using Sobol sensitivity methods.
https://doi.org/10.1371/journal.pone.0198579.g003  (Table B in S1 File). By comparing the rank correlations, we can see that the parameter sensitivity spectrum according to bifurcation point position is similar to the spectrum obtained from the energy barrier analysis near bifurcation point.

Conclusion
For a bistable system, people may intuitively think that when the location of bifurcation point shifts, the energy landscape changes. However, the specific quantitative changes have scarcely been studied, particularly in complex biological systems. Using mitochondrial apoptosis as an example, we specifically compared the changes of bifurcation point and energy barrier height. We successfully calculated MFPT and energy barrier height using the SSA and GMAM algorithms. The resulting stochastic model provides the transition from different basins of attraction of metastable states, which cannot be described from the deterministic model. We observe activation of apoptosis without high levels of Cas8 triggering bifurcation from the stochastic model. Additionally, when a parameter change causes a shift of the bifurcation point, there is a corresponding change of the energy barrier height and vice versa. The higher the barrier height is, the more difficult it is for a cell to evolve to death, indicating increasing oncogenesis potential. A right shift of the bifurcation point also implies higher oncogenesis potential. These results confirmed the consistency between them.
We rank the parameters according to their influence on MFPT and energy barrier. Compared with the results of parameter ranking obtained from deterministic dynamic studies, we find a high consistence of parameter sort orders. The 'equivalence' of bifurcation point and energy barrier height in determining parameter sensitivity spectrum is also shown in Schlögl model which is a simple tri-molecular reaction model that can give bistability. The corresponding analysis process can be found in the Supporting Information. From sensitivity analysis shown above, we find that the production and degradation rates of Bcl-2 and Bax (bcl2_pro, bcl2_deg, bax_pro, and bax_deg) are the most sensitive, thus illustrating the crucial role of these two proteins. The increase of bcl2_pro, bax_deg and the decrease of bcl2_deg, bax_pro correspond to higher oncogenesis potential. This result is consistent with those of previous studies [31,32]. For example, chromosomal translocation-induced Bcl-2 overexpression and other amplification cases have been observed in many cancer types [33][34], so it is not surprising that the parameter describing Bcl-2 production stands out as the most sensitive one. Similarly, the frameshift mutations in a (G)8 mononucleotide tract of the Bax gene, which lead to synthesis of a truncated protein and reduced expression, commonly occur in colon and gastric cancers [35], suggesting that Bax inactivation during tumorigenesis may facilitate tumor progression by enhancing escape from apoptosis. Parameter bcl2_pro in the legend stands for production rate of Bcl2 and bcl2_deg represents degradation rate of Bcl2 that are the two most sensitive parameters. Blue, wild type parameters; green, bcl2_pro decreased by 2%; red, bcl2_pro increased by 2%; cyan, bcl2_deg decreased by 2%; and purple, bcl2_deg increased by 2%. (c) Parameter sensitivity spectrum of the quasi-potential barrier height when Cas8 = 22 nM, which is close to the bifurcation -point location of the corresponding deterministic dynamics. A 2% decrease or increase of each parameter induces the percentage change of energy barrier height according to the local sensitivity analysis method. (d) Parameter sensitivity spectrum of the quasi-potential barrier height when Cas8 = 4.9 nM, which is relatively far from the bifurcation point. A 2% decrease or increase of each parameter induced the percentage change of energy barrier height according to the local sensitivity analysis method. https://doi.org/10.1371/journal.pone.0198579.g005 For the model we considered here, we have already shown that the parameters significantly affecting the system's bifurcation point are closely correlated with high frequency mutations [8]. Thus, the bifurcation-point location can be considered as indicator of the potential of cell to become cancerous. Taken together, we conclude that energy barrier height can also be regarded as an indicator of oncogenesis.
In the mitochondrial apoptosis pathway, right shifts of bifurcation point render higher oncogenic potential. Consequently, a far-more-right position of the bifurcation point may render a normal cell cancerous. Pre-treatment of mitochondria that moves a cell near the threshold of apoptosis could be useful as an anti-cancer strategy.