BACKGROUND
apoptosis, the major form of programmed cell death, is a conserved cell suicide process critical for the health and survival of multicellular organisms  <cit> .
apoptosis plays a fundamental role in animal development, by sculpting tissues and structures, as well as in tissue homeostasis, by regulating and maintaining balanced cell number  <cit> .
dysregulation of apoptosis is associated with various human diseases, ranging from developmental disorders, neurodegeneration to cancer  <cit> .
apoptosis is regulated by two interrelated signaling pathways: the extrinsic or death-receptor pathway, and the intrinsic or mitochondrial pathway  <cit> .
they converge on the execution pathway, mediated intracellularly by a cascade of cysteine proteases, termed caspases  <cit> .
caspases are specialized cysteine proteases found in animal cells as inactive procaspases .
through proteolytic cleavage, procaspases are activated to carry out its apoptotic mission.
the intrinsic pathway begins with the release of cytochrome c  from mitochondria through membrane permeabilization  <cit> , triggered by intracellular stress such as dna damage and hypoxia  <cit>  .
once cc translocates to the cytosol, it binds to apoptotic protease activating factor  <dig>  to form a multimeric protein complex called the apoptosome.
this apoptosome complex then activates the initiator procaspase, called procaspase- <dig>  the activated caspase  cleaves the executioner procaspase  to form active executioner/effector caspase , whereby the apoptotic response is irreversibly triggered  <cit> .
experiments have shown that the activation of effector caspases occurs in an all-or-none fashion, emphasizing the functional role of the apoptosis system as a molecular switch.
in the past several years, advances in experimental skills have allowed the measurement of apoptosis dynamics in individual cells  <cit> , confirming the switch-like dynamics, while revealing another feature of prominent stochasticity in the apoptotic responses at single-cell level.
since the key constituents and molecular interactions of apoptosis pathways have been experimentally identified, the approach of mathematical modeling and computer simulations have been employed extensively to help elucidate the complicated regulatory network and dynamic responsiveness related to apoptosis at average cellular population level  <cit> .
nevertheless, recent experiments at single-cell resolution have discovered noisy phenotypic diversity of apoptosis activity in that significant cell-to-cell heterogeneity of the dynamic apoptosis responses exist across a genetically-identical cell population  <cit> .
toward understanding such single-cell variability in apoptosis response, some theoretical efforts have been taken recently to model the stochastic response of receptor-mediated apoptotic pathway.
the stochastic behavior of intrinsic apoptosis pathway, on the other hand, has been the subject of relatively little mathematical modeling to date.
in this work, we will focus on addressing the intrinsic apoptosis pathway under stochastic perturbations by developing theoretical and computational models at single-cell level.
the models will be exploited to investigate the heterogeneous behavior of intrinsic apoptosis network among individual cells.
deterministic model based on ordinary differential equations  is the most widely used mathematical approach to describe the molecular kinetics during cell death signaling.
fussenegger et al.
developed a well-accepted ode model that integrates components of the extrinsic as well as the intrinsic apoptosis pathways  <cit> .
qualitatively the fussenegger model compares reasonably well with published experimental kinetics of caspase activation at average cell population level.
nevertheless, there is lack of understanding of the nonlinear stability and systems properties of this model, which hinders deeper understanding of the system behavior.
for instance, studies have suggested that bistability is a key system feature for apoptotic signaling networks  <cit> , which could achieve the all-or-none responses and in addition confer robustness to the apoptosis system  <cit> .
it is unclear whether the fussenegger model presents the property of bistability.
since then, there have been considerable theoretical efforts on modeling and systems analysis using ode models of death-receptor mediated apoptosis  <cit> , mitochondria-mediated apoptosis  <cit> , or integrated extrinsic and intrinsic apoptosis pathways  <cit> .
the past few years have seen increasing efforts in stochastic modeling to address the heterogeneous apoptosis responses at single-cell level.
specifically, these efforts incorporate cellular noise perturbations into the apoptosis framework.
cellular noise is defined as stochastic fluctuations of biomolecular processes within and between cells.
it can be divided into intrinsic noise and extrinsic noise  <cit> .
intrinsic noise in genetically identical cells refers to random deviation of the molecular processes from their average deterministic kinetics within a cell, mostly due to probabilistic biochemical reactions associated with low copy number of molecular quantities  <cit> .
extrinsic noise arises from global perturbation factors such as cellular environment and organelle distribution, which results in cell-to-cell variation in rate constants of biochemical reactions, expression levels of genes and proteins, and other parameters of biochemical processes  <cit> .
towards the analysis of cellular noise, several statistical measures of noise have been proposed to quantify the level of stochastic fluctuations of biomolecular processes  <cit> .
two measures of noise are commonly used to characterize the stationary averages and variances of random cellular components.
in particular, noise strength can be quantified by fano factor, which is defined as the steady-state variance over average and has a value of  <dig> for poisson process.
the fano factor of an arbitrary stochastic system reveals deviations from poissonian behavior  <cit> .
a more standard and frequently used measure of cellular noise is the dimensionless coefficient of variation, which is defined as standard deviation divided by mean.
it measures the inverse signal-to-noise ratio and has been widely employed to characterize intrinsic and extrinsic noises of gene and protein expression and their determining factors with respect to cellular network organization  <cit> .
for the latter measure, the coefficient of variation squared may be alternatively used  <cit> .
in this work, we use the coefficient of variation to quantify the noise of the random distribution of molecular components and stochastic cellular response time.
in the aspect of mathematical modeling of noise in apoptosis pathway, several previous stochastic apoptosis models have taken into account of the intrinsic noise by either applying gillespie’s stochastic simulation algorithm to the ode models or constructing monte carlo models from first principles  <cit> .
with regard to modeling the impact of extrinsic noise on apoptosis pathway, there have been a few studies notably only on receptor-mediated apoptosis pathway  <cit> .
in this study, we attempt to develop mathematical and computational models of the intrinsic apoptosis pathway at single-cell level, and to identify plausible sources of non-genetic heterogeneity of apoptosis dynamics observed in a cell population using stochastic simulations.
we start with a deterministic ode model of intrinsic apoptosis pathway adapted from the fussenegger model and find that bistability is missing.
by adding positive feedback regulations that are supported by previous experimental evidences, we develop a model of intrinsic apoptosis pathway that functions as a bistable switch.
we are particularly interested in understanding the stochastic behavior of this apoptosis switch under perturbation of intrinsic noise and/or extrinsic noise.
stochastic modeling and simulations of intrinsic apoptosis pathway indicate that noise could enhance robustness of the bistable switch.
in addition, we show that intrinsic noise is not sufficient to induce the observed level of cell-to-cell variability of apoptosis response at biologically relevant level of molecular numbers, while the extrinsic noise of protein variations is plausibly the main source giving rise to the degree of heterogeneous responses of intrinsic apoptosis pathway between single cells.
RESULTS
deterministic model of intrinsic apoptosis pathway
to build a single-cell level deterministic model of intrinsic apoptosis pathway, we adopt the ode modeling framework of intrinsic apoptosis proposed by fussenegger et al.
<cit>  as it has the “minimal” model complexity while preserving all the critical interactions of intrinsic apoptosis pathway.
the model scheme is shown in figure 1b and summarized here: the intrinsic apoptosis pathway is initiated by the release of cytochrome c  from mitochondria to cytosol, and binding to apaf- <dig> to form the apoptosome complex .
the apoptosome catalyzes procaspase- <dig> , the precursor of initiator caspase, to its active form, caspase  <dig> .
the executioner procaspase- <dig>  is then activated by c9a to form the active executioner caspase  <dig> , whose response represents the onset of the irreversible apoptosis fate.
finally, cea is subject to the direct inhibition by iap protein  <cit>  .
simulations of this ode model demonstrate that the time trajectories of cea starting from an initially low concentration first undergo some time delay and then switch to a high steady state with a relatively sharp transition .
nevertheless, the time trajectories of cea from different initial conditions all converge to the same steady-state level, indicating the existence of single equilibrium point.
indeed, bifurcation analysis of this ode model with respect to the parameter of cc concentration shows that the steady-state of cea is monostable with a sigmoidal input-output relationship .
this means that regardless of different input strengths and different initial conditions, the response would asymptotically settle at only one steady state.
further analytical analysis of the steady-state response of the fussenegger model proves that the mapping between the input of cc signal and the output of cea is either one-to-one or one-to-two and thus there exist at most two solutions of the output signal at equilibrium .
as a result, the fussenegger model of intrinsic apoptosis pathway cannot be bistable as bistability of a system requires three distinct steady-state solutions in some range of input.
the fussenegger model of intrinsic apoptosis pathway presents a system property of so-called ultrasensitivity  rather than bistability  <cit> .
this, however, is inconsistent with the current consensus of bistability feature of apoptosis system, and thus extra model refining step is needed to resolve the discrepancy.
to this end, we modify the existing model by incorporating a positive feedback from the executioner caspase  to the initiator caspase , an auto-catalysis loop to the activation kinetic of c9a, and a mild cooperativity in cea activation induced by c3p, which are supported by previous experimental and computational findings   <cit> .
the existing mathematical models accounting for intrinsic apoptosis pathway differ from our proposed model in choice of apoptosis pathway.
specifically, we follow the fussenegger model to assume that the inhibition of caspase  <dig> by iap is not essential and ignored in our model.
this choice is the main distinction from the model proposed by legewie et al <cit> , where iap inhibits both caspase  <dig> and caspase  <dig> thus giving rise to an iap-mediated positive feedback and bistability.
the harrington et al model  <cit> , the zhang et al model <cit> , and the kutumova et al model  <cit>  all inherit the schematic of the intrinsic apoptosis model proposed by legewie et al, and thus also include the double inhibitions of caspase  <dig> and caspase  <dig> by iap and the resulting implicit positive feedback.
the governing odes of the modified model of intrinsic apoptosis pathway are listed in the methods section.
simulations of the modified ode model show that when given a step input of low concentration of cc, the time courses of cea gradually settle at a near-zero steady state starting from different initial cea concentrations, while given a relatively high concentration of cc, cea eventually settles at a high steady state .
such behavior with two stable output steady states indicates that bistability is achieved by the modified ode model.
in addition, the time trajectories agree with experimental results in that the cea response is not elicited until after a few hours of delay time .the switching-on kinetic of cea activation is sigmoidal shape and completed within ∼1hr, presenting all-or-none switch-like behavior  <cit> .
indeed, one-parameter bifurcation analysis of the modified ode model confirms that the steady-state response of cea is bistable with respect to the input signal of cc, where two stable steady states coexist between the input concentrations of  <dig>  μm and  <dig>  μm .
the bifurcation curve has two saddle-node bifurcation points sn <dig> and sn <dig>  giving rise to a middle unstable branch and two stable branches, where the upper and lower branches correspond to the apoptosis and survival fates, respectively.
when the concentration of cc approaches sn <dig> and sn <dig>  hysteretic behavior occurs: the system remains at near-zero cea activity at low amount of cc, until an on threshold  is reached, whereby cea activity switches to the apoptosis state abruptly; inversely, the cea activity switches from the apoptosis state to the survival state only if the cc concentration falls below the off threshold .
the system properties of bistability and hysteresis could confer robustness to the apoptotic responsiveness by allowing cells that are not committed to apoptosis to remain at survival state, even in the event of mildly noisy input.
in addition, two-dimensional bifurcation analysis with respect to four selected parameters, namely the hill constant that regulates the cea-mediated positive feedback loop , the cooperativity of activation of cea , the degradation rate of c9a and the inhibition rate of cea by iap , show that the bistability property of the modified model exists in extended parameter space around the nominal parameter set  and is hence a robust behavior.
stochastic model of intrinsic apoptosis pathway under intrinsic noise perturbation
the deterministic model of the intrinsic apoptosis pathway has allowed us to analyze nonlinear properties of the system and quantify region of robust behavior.
nevertheless, the ode modeling approach accounts for average cellular dynamics while ignoring the inevitable unpredictability embedded in biomolecular reactions and in intra- and extra-cellular environments.
it has been observed by several different experimental techniques that the apoptosis response at single-cell resolution is subject to inherent stochastic perturbations, giving rise to pronounced cell-to-cell variability even in genetically identical cell population  <cit> .
therefore, it is necessary for us to develop a stochastic modeling framework of intrinsic apoptosis pathway, which can be used to explore the plausible origin of the stochasticity underlying the phenotypic heterogeneity.
in this study, we investigate the impact of both intrinsic-noise and extrinsic-noise perturbations on the intrinsic apoptosis dynamics.
first, to model the stochasticity under the intrinsic-noise perturbation due to low copy number of biomolecules, we assume that the deterministic model represents the nominal single-cell behavior, and apply the standard gillespie ssa   <cit> .
the detail of applying ssa to the deterministic model can be found in the methods section.
to reflect the plausible low molecular copy numbers arising in the intrinsic apoptosis pathway within single cells, we perform simulations under different combinations of number of molecules that fall into the experimentally observed ranges  <cit> .
furthermore, in choosing number of molecules used for stochastic simulation of intrinsic noise, we limit the number to be within  in general, as previous experimental and computational work have suggested that the coefficient of variation of cellular signal due to intrinsic noise reaches a minimal level when the mean abundance of molecules increases to  <dig> <cit> .
using this limit and previous experimental observation that procaspase- <dig>  and procaspase- <dig>  could have relatively low number of molecules per cell, ranging respectively from 5× <dig> to  <dig> × <dig> and from undetectable to  <dig> × <dig> <cit> , we choose to run ssa of the model at molecule numbers of c9p = {5× <dig> } and c3p = { <dig> , <dig> }, which correspond to a total of eight combined conditions.
note that each realization of stochastic simulation corresponds to the apoptosis response of one single cell.
figure 4a shows  <dig> stochastic time trajectories of output cea response with  =, representing the behavior of  <dig> cells, under varying amount of cytochrome c input.
each time course of cea activation exhibits sigmoidal switch-like shape, converging to an elevated steady state.
to illustrate the cellular variability in the stochastic response due to intrinsic noise, histograms of the steady-state cea copy numbers in the same  <dig> cells is plotted in figure 4b against different input cc levels.
the histogram shows that if the cc molecule is above  <dig> copy number the distribution of steady-state cea is bimodal, with a low mean steady-state value of ∼ <dig> number of cea molecules and a high mean steady-state value of ∼ <dig> cea molecules.
such bimodal distribution of cea response indicates that the stochastic response of intrinsic apoptosis pathway subject to intrinsic noise under this particular c9p/c3p condition is bistable  <cit> .
the bistability behavior persists even when the copy number of cc increases to ∼ <dig> molecules, showing that the corresponding fold change of bistability domain under intrinsic noise is above five times that of the deterministic model, where the bistability region of cc is   as shown in figure 2d.
such phenomenon of enhanced robustness induced by intrinsic noise is in agreement with previous computational work which suggests that stochastic signaling networks may perform more robustly than their deterministic counterpart  <cit> .
it is noteworthy that the existence and range of bistable cea response are dependent on the abundances of c9p and c3p in that the bistability boundary shrinks as the copy numbers of procaspases increase, and the bistable range becomes almost undetectable at significantly high amount of procaspases  .
this result implies that the bistability property of apoptosis might only be observable at proper condition of molecular abundance.
since figure 4b has demonstrated that the apoptosis state of cea varies from cell to cell with a normal distribution, next we evaluate the degree of stochasticity of the cea response at equilibrium caused by intrinsic noise.
we use the coefficient of variation  to quantify the level of stochasticity, following previous studies of isogenic cell-to-cell variability of apoptosis response  <cit> .
as shown in figure 5a and 5b, under all the eight conditions of molecule numbers as described above, the mean value of the stochastic cea response increases while its standard deviation decreases as the input signal cc increases, both with exponential pattern.
note that the condition of higher c9p and c3p abundances give rise to elevated mean response of cea at equilibrium as well as sharper stimulus-response curve, agreeing with the pro-apoptotic roles played by these two components.
the consequent stochasticity of the output cea value, quantified by cv, monotonically decreases with increasing amount of input and reaches a minimum level as the number of cc molecules is above  <dig> for all the eight conditions .
this indicates that the intrinsic noise in the steady-state behavior of intrinsic apoptosis response at single cell level can be suppressed by enhancing the strength of input signal.
note also that lower abundances of c9p and c3p induce more prominent stochasticity in output when the input amount of cc is sufficiently low.
we then focus on an important quantity that has been experimentally measured to evaluate the cell-to-cell variability in apoptosis response.
as demonstrated in figure 6a, each stochastic time trajectory of cea presents a time delay between time zero, the cc release time, and the cea activation time, which is conventionally defined as the half-maximal caspase  <dig> cleavage time.
such delay time  is defined analogously to the quantity previously used in experiments to characterize the response time of extrinsic apoptosis pathway in single cells  <cit> .
and experimental quantifications suggest that td varies from cell to cell, even in genetically identical clones  <cit> .
our single-cell stochastic simulations with intrinsic noise, under the eight conditions of c9p and c3p abundances, show that in general the mean value of td monotonically decreases with increasing cc level .
for instance, when c9p and c3p are respectively  <dig> and  <dig> molecules per cell, the mean value of td decreases from around  <dig> min at  <dig> number of cc molecules to ∼ <dig> min at  <dig> number of cc molecules.
understandably, lower abundances of pro-apoptotic proteins c9p and c3p lead to longer response delay time.
we then quantify the degree of cell-to-cell stochasticity by cv of td.
figure 6c shows that the cv of td due to intrinsic noise is, under all the conditions, a monotonically decreasing function with regard to the number of cc molecules.
such findings of attenuated cellular variability of response under the circumstance of elevated input signal are consistent with previous experimental and computational results on the stochastic apoptosis pathway in genetically identical cells as well as other cellular processes  <cit> .
the above stochastic simulations show that the cell-to-cell variability of intrinsic apoptosis due to intrinsic noise is especially evident when the number of c3p is relatively small , where the cv of td persist at the level of non-genetic noise of apoptosis observed by experiments   <cit>  .
moreover, the abundances of c9p and c3p inversely regulate the cellular variability among the eight conditions.
specifically, the lowest c9p/c3p copy-number condition leads to a cv slightly above  <dig>  while the highest c9p/c3p copy-number condition leads to a cv slightly below  <dig>  , suggesting that the wide copy-number range of procaspases, especially procaspase- <dig>  under the sub- <dig> condition is likely a source of cell-to-cell variability under the perturbation of intrinsic noise.
previously, there has been study of apoptosis pathway implying that the molecular numbers of participating biochemical species seem to reside in a regime much higher than  <dig>   <cit> .
therefore, our stochastic model with intrinsic noise is further explored under the various copy numbers of c9p and c3p with sufficiently strong cc input.
as shown in figure 6d, under all the eight conditions the cv of td drops dramatically when the number of cc molecules increases above  <dig>  in particular, when the copy numbers of c9p, c3p and cc molecules are all raised to  <dig> and above, the cv of td tends to an almost negligible level of ∼ <dig> .
such low degree of cv of td at physiologically more plausible condition indicates that the perturbation by intrinsic noise alone seems insufficient to induce the observed degree of cell-to-cell stochasticity of apoptosis dynamics , and other sources of uncertainty needs to be taken into account.
stochastic model of intrinsic apoptosis pathway under extrinsic noise perturbation
recently, natural protein variations across cell population has been identified as a major source of extrinsic fluctuations for apoptosis pathway  <cit> .
experiments have suggested that the concentration of a protein naturally varies among different cells following a log-normal distribution with a typical cv value of  <dig>  to  <dig>   <cit> .
to address the impact of this kind of extrinsic noise on intrinsic apoptosis pathway, we assume that in the deterministic single-cell model of intrinsic apoptosis pathway described above certain protein of interest has log-normally distributed concentration with cv equal to  <dig> .
first, we simulate the stochasticity due to variation in the concentration of the protein cc, the critical input signal of intrinsic apoptosis pathway, using the single-cell deterministic model with a randomly generated cc concentration.
again, as a measure of cell-to-cell variability, the cv of td among  <dig> cells is calculated at varying mean values of cc protein .
the variability peaks  at mean concentration of cc =  <dig>  μm, and then declines as the mean abundance of cc protein increases, settling at a value of cv = ∼ <dig>  eventually.
it indicates that the natural variation in the cc protein can induce some degree of cell-to-cell stochasticity to the intrinsic apoptosis response, especially at relatively low abundance of cc protein.
second, we evaluate the impact of extrinsic noise driven by the natural variation of iap protein, which is one of the most pivotal anti-apoptotic proteins tightly regulating apoptosis via antagonizing the activity of cea  <cit> .
the single-cell deterministic model plus a log-normally distributed amount of iap with cv equal to  <dig>  is simulated for  <dig> cells at different concentrations of input cc protein.
here, we plot the resultant cv of td evaluated among  <dig> cells in figure 8a as a 2-dimensional  heat map, versus the concentration of cc protein and the mean concentration of iap protein.
it shows that given the same cc input, higher mean concentration of iap, albeit with the same dispersion of random distribution, leads to relatively larger cv of td and hence higher degree of cell-to-cell variability in the timing of cea activation.
most likely, large amount of iap significantly represses the activity of cea, while it is generally known that the lower amount of a biochemical component the higher degree of its stochasticity.
additionally, figure 8a indicates that at the same mean abundance of stochastic iap protein, higher input level of cc protein yields lower cell-to-cell variability in the apoptosis response.
the above findings together imply that the extrinsic noise in cc protein has opposite influence on the induction of cell-to-cell variability to intrinsic apoptosis response than the extrinsic noise in iap protein.
such contrary effect is in line with their opposite regulatory roles in apoptosis pathway, in that cc is pro-apoptotic while iap is anti-apoptotic.
notably, figure 8a also shows that if the mean concentration of iap is moderately high , the cv of td can achieve a value of  <dig>  or even higher, which is the level of cell-to-cell stochasticity measured by experiments.
hence the extrinsic noise in moderate amount of iap protein may serve as a significant source of cell-to-cell variability in intrinsic apoptosis response.
but if iap is so high that paramount of inhibition is produced, the apoptosis response of cea is then constantly switched off, denoted by the upper-left non-responsive region in the 2d heat map .
it is interesting to find out that when the 2d heat map is superimposed with the boundary of the two-parameter bistability diagram of the stochastic model, which is estimated using the histograms of cea response , the cv of td with relatively high values  is located between the low-threshold curve  and the high-threshold curve .
such behavior is confirmed under systematic perturbations of the parameter kc  and the other three parameters selected for the sensitivity analysis in figure  <dig> .
compared to the threshold lines for the bistability domain of the corresponding deterministic model , the extrinsic noise in cc seems not significantly affect the area of bistability domain, albeit the thresholds of cc shift toward smaller values.
this trend also holds under different values of kc .
as a further exploration of the impact of extrinsic noise, we allow the concentrations of both the cc and iap proteins to be log-normally distributed random numbers with cv equal to  <dig> .
the resulting 2d heat map of the cv of td in figure 8b shows that the level of stochasticity across cells is in general higher than the case under individual extrinsic perturbation of cc protein or iap protein.
in most of the area in the 2-parameter region of figure 8b, the cv of td achieves a value above  <dig> .
similar to figure 8a, the cv of td with double extrinsic noises attains highest values  between the low-threshold and high-threshold curves estimated for the stochastic bistability diagram, which is almost twice the cv value in the scenario under the extrinsic noise in iap protein only.
moreover, the activation of cea response can now be elicited even in the non-responsive region of figure 8a due to the additional degree of fluctuation in the input signal.
collectively, the 2d heat-map in figure 8b suggests that the extrinsic noises in cc and iap proteins are sufficient to yield the experimentally observed degree of cell-to-cell variability in the apoptosis response.
the bistability domains under the nominal value of kc in figure 8b and those under systematically perturbed values of kc in additional file 4: figure s <dig> demonstrate the same trend of behavior as those in figure 8a and additional file 3: figure s <dig>  close comparison between figure 8a/additional file 3: figure s <dig> and figures 8b/additional file 4: figure s <dig> indicate that the area of bistability domain on average is slightly larger in the latter group.
therefore, additional degree of extrinsic noise may induce extra robustness for bistability of apoptosis.
stochastic model of intrinsic apoptosis pathway under combined intrinsic and extrinsic noise perturbations
we have so far analyzed the model of intrinsic apoptosis pathway subject to either intrinsic noise or extrinsic noise independently.
additionally, we would like to find out if both types of noises are present in the intrinsic apoptosis pathway, how the cell-to-cell variability of the delay time of cea response is influenced.
to simulate a model under the perturbation of combined intrinsic and extrinsic noises, we implement the stochastic simulation algorithm of the apoptosis model as described above to mimic the intrinsic noise, and simultaneously allow certain protein concentrations to be log-normally distributed random variables to represent the extrinsic noise.
in order to make comparison with the results of the intrinsic-noise only case, we again use the eight conditions of molecule numbers at c9p = {5× <dig> } in combination with c3p = {  <dig> , <dig> }.
first, the impact of extrinsic fluctuation in the amount of cc protein on top of the intrinsic noise is simulated.
the resulting cv of td in response to the abundance of cc less than  <dig> is shown in the inset of figure 7b.
note that the unit of cc protein is now copy number of molecules rather than μm as a requirement by the implementation of ssa.
comparing it to figure 6c, we find that the behavior of cell-to-cell variability due to the combined types of noises is almost the same as that under intrinsic noise alone.
that is, under all the eight molecule-number conditions, the cv of td monotonically decreases in a near exponentially-decaying fashion as the number of cc molecules increases, and it falls into similar range of value .
and similar to figure 6c, the variability value is inversely correlated with the abundances of c9p/c3p, achieving cv of td above  <dig>  if c9p/c3p have sub- <dig> copy numbers.
this result indicates that at the limit of low number of molecules , the intrinsic noise seems to make the dominating contribution to the cell-to-cell stochasticity of intrinsic apoptosis response.
we subsequently run the same stochastic model under the perturbation of intrinsic noise plus extrinsic noise in cc, while increasing the mean abundance of cc molecules up to  <dig> , a more plausible reacting scale for apoptosis pathway.
we find that the stochasticity curves under different c9p/c3p abundance conditions relatively converge at cc equal to  <dig>  copy number , due to diminished intrinsic noise.
the comparison between figure 7b and figure 6d enables us to see the relative contributions of intrinsic and extrinsic noises: except for the case of almost undetectable abundance of c3p , which induce severe intrinsic noise, the cv of td due to combined noises is significantly higher than that due to intrinsic-only noise.
specifically, the contrast between the case of combined noises and the case of intrinsic-only noise is highest  when the copy numbers of c9p, c3p and cc molecules are all above  <dig>  therefore, at higher molecular numbers, the cellular stochasticity seems to be majorly contributed by the source of extrinsic noise.
with respect to the possible impact of the extrinsic noise in iap, we simulate the stochastic model with additional perturbation in the amount of iap protein.
that is, both the abundances of iap and cc proteins are assumed to be log-normally distributed around their respective mean value with a cv of  <dig> , alongside with intrinsic noise .
at the limit of low cc abundance the additional extrinsic noise of iap does not introduce extra cell-to-cell stochasticity to apoptosis response .
when the mean number of cc molecules is in the physiologically plausible range of  <dig>  to  <dig> , simulations show that td presents a high cv, mostly above the value of  <dig>  and even going up to a value of  <dig> in the 2d heat map .
we compare the cell-to-cell variability under combined noises  with that under extrinsic-only noise , and find that the cv of td increases only slightly in the former case .
this is consistent with the minor role of intrinsic noise played at elevated abundances of molecules.
as a consequence, the boosted degree of stochasticity in figure  <dig> has to owe predominantly to the extrinsic fluctuations in the amount of cc and iap proteins rather than the intrinsic noise.
simulations under the eight conditions of c9p and c3p abundances demonstrate that the overall cell-to-cell variability is, the same as what we see earlier, inversely correlated with the copy numbers of molecules, where the cv of td on average declines through row  <dig> to row  <dig> as the abundance of c3p increases from  <dig> to  <dig>  while the cv of td on average slightly declines from column  <dig> to column  <dig> as the abundance of c9p increases from 5× <dig> to  <dig>  this is consistent with the aforementioned observation of the evident variability induced by the procaspase- <dig> under sub- <dig> copy-number condition through the channel of intrinsic noise.
it is noteworthy that the cv of td in figure  <dig> is generally larger at higher mean value of iap molecular number, agreeing with the inference from figure 8a that the abundance of iap protein is capable of promoting the stochasticity of apoptosis pathway.
it highlights the interesting role of iap protein in controlling the cell-to-cell variability of intrinsic apoptosis response, and suggests that in treating diseases exploiting the apoptosis mechanism, such as cancer, the iap protein offers a potential therapeutic target not only for effective modulation of apoptosis  <cit>  but also for eliminating the undesired cell-to-cell heterogeneity, a major obstacle to effective cancer treatment and personalized medicine  <cit> .
this finding is in line with the current view of the therapeutic function of iap based on the study of apoptosis pathway at cell population level  <cit> .
CONCLUSIONS
the recently observed heterogeneous apoptosis phenotypes at single cell level have drawn increasing attention from researchers.
mathematical modeling and computer simulation provide an efficient approach to gain deep insight into the dynamic behavior of apoptosis network.
this paper develops a theoretical and computational framework for single-cell stochastic modeling of the intrinsic apoptosis pathway.
using this modeling framework we explore the stochastic behavior of the intrinsic apoptosis response at single-cell level and seek to understand the plausible sources underlying the experimentally observed cell-to-cell variability of apoptosis response.
we show that in the presence of noise, the bistable response of intrinsic apoptosis pathway can be more robust than its deterministic behavior.
the coefficient of variation  of the delayed timing of the activity of executioner caspase is utilized to quantify the stochasticity in the apoptosis dynamics.
we find that the intrinsic noise can introduce significant cell-to-cell variability if the abundances of reacting biomolecules are relatively low.
the level of cellular stochasticity solely due to intrinsic noise decreases dramatically to a negligible level of cv equal to ∼ <dig>  when the copy number of cytochrome c is raised to  <dig> , which is the amount suggested by a previous study.
in addition, the extrinsic noise caused by the natural variations in protein concentrations of two key components in the intrinsic apoptosis pathway, cytochrome c and the inhibitor of apoptosis  proteins, is also accounted for without or with intrinsic noise.
series of simulations indicate that the extrinsic noise is plausibly the major source of the cell-to-cell variability of intrinsic apoptosis response at high number of biomolecules.
furthermore, we find that the mean abundance of fluctuating iap is positively correlated with the degree of cell-to-cell variability, thus making iap a potential target for therapeutically suppressing the stochasticity of intrinsic apoptosis response across cell population in treating diseases such as cancer.
in summary, this study based on our theoretical and computational models characterizes the behavior of the intrinsic apoptosis pathway under complex stochastic perturbations, and suggests that certain deterministic features, such as system bistability and iap as potential therapy target, still remain in the presence of noise.
altogether, the work can enable us to gain deeper understanding toward the experimentally observed uncertainty in cellular decision making.

 apoptosis, the major form of programmed cell death, is a conserved cell suicide process critical for the health and survival of multicellular organisms. apoptosis plays a fundamental role in animal development, by sculpting tissues and structures, as well as in tissue homeostasis, by regulating and maintaining balanced cell number. dysregulation of apoptosis is associated with various human diseases, ranging from developmental disorders, neurodegeneration to cancer. apoptosis is regulated by two interrelated signaling pathways: the extrinsic or death-receptor pathway, and the intrinsic or mitochondrial pathway. they converge on the execution pathway, mediated intracellularly by a cascade of cysteine proteases, termed caspases. caspases are specialized cysteine proteases found in animal cells as inactive procaspases. through proteolytic cleavage, procaspases are activated to carry out its apoptotic mission. the intrinsic pathway begins with the release of cytochrome c from mitochondria through membrane permeabilization, triggered by intracellular stress such as dna damage and hypoxia. once cc translocates to the cytosol, it binds to apoptotic protease activating factor to form a multimeric protein complex called the apoptosome. this apoptosome complex then activates the initiator procaspase, called procaspase- the activated caspase cleaves the executioner procaspase to form active executioner/effector caspase, whereby the apoptotic response is irreversibly triggered. experiments have shown that the activation of effector caspases occurs in an all-or-none fashion, emphasizing the functional role of the apoptosis system as a molecular switch. in the past several years, advances in experimental skills have allowed the measurement of apoptosis dynamics in individual cells, confirming the switch-like dynamics, while revealing another feature of prominent stochasticity in the apoptotic responses at single-cell level. since the key constituents and molecular interactions of apoptosis pathways have been experimentally identified, the approach of mathematical modeling and computer simulations have been employed extensively to help elucidate the complicated regulatory network and dynamic responsiveness related to apoptosis at average cellular population level. nevertheless, recent experiments at single-cell resolution have discovered noisy phenotypic diversity of apoptosis activity in that significant cell-to-cell heterogeneity of the dynamic apoptosis responses exist across a genetically-identical cell population. toward understanding such single-cell variability in apoptosis response, some theoretical efforts have been taken recently to model the stochastic response of receptor-mediated apoptotic pathway. the stochastic behavior of intrinsic apoptosis pathway, on the other hand, has been the subject of relatively little mathematical modeling to date. in this work, we will focus on addressing the intrinsic apoptosis pathway under stochastic perturbations by developing theoretical and computational models at single-cell level. the models will be exploited to investigate the heterogeneous behavior of intrinsic apoptosis network among individual cells. deterministic model based on ordinary differential equations is the most widely used mathematical approach to describe the molecular kinetics during cell death signaling. fussenegger et al. developed a well-accepted ode model that integrates components of the extrinsic as well as the intrinsic apoptosis pathways. qualitatively the fussenegger model compares reasonably well with published experimental kinetics of caspase activation at average cell population level. nevertheless, there is lack of understanding of the nonlinear stability and systems properties of this model, which hinders deeper understanding of the system behavior. for instance, studies have suggested that bistability is a key system feature for apoptotic signaling networks, which could achieve the all-or-none responses and in addition confer robustness to the apoptosis system. it is unclear whether the fussenegger model presents the property of bistability. since then, there have been considerable theoretical efforts on modeling and systems analysis using ode models of death-receptor mediated apoptosis, mitochondria-mediated apoptosis, or integrated extrinsic and intrinsic apoptosis pathways. the past few years have seen increasing efforts in stochastic modeling to address the heterogeneous apoptosis responses at single-cell level. specifically, these efforts incorporate cellular noise perturbations into the apoptosis framework. cellular noise is defined as stochastic fluctuations of biomolecular processes within and between cells. it can be divided into intrinsic noise and extrinsic noise. intrinsic noise in genetically identical cells refers to random deviation of the molecular processes from their average deterministic kinetics within a cell, mostly due to probabilistic biochemical reactions associated with low copy number of molecular quantities. extrinsic noise arises from global perturbation factors such as cellular environment and organelle distribution, which results in cell-to-cell variation in rate constants of biochemical reactions, expression levels of genes and proteins, and other parameters of biochemical processes. towards the analysis of cellular noise, several statistical measures of noise have been proposed to quantify the level of stochastic fluctuations of biomolecular processes. two measures of noise are commonly used to characterize the stationary averages and variances of random cellular components. in particular, noise strength can be quantified by fano factor, which is defined as the steady-state variance over average and has a value of for poisson process. the fano factor of an arbitrary stochastic system reveals deviations from poissonian behavior. a more standard and frequently used measure of cellular noise is the dimensionless coefficient of variation, which is defined as standard deviation divided by mean. it measures the inverse signal-to-noise ratio and has been widely employed to characterize intrinsic and extrinsic noises of gene and protein expression and their determining factors with respect to cellular network organization. for the latter measure, the coefficient of variation squared may be alternatively used. in this work, we use the coefficient of variation to quantify the noise of the random distribution of molecular components and stochastic cellular response time. in the aspect of mathematical modeling of noise in apoptosis pathway, several previous stochastic apoptosis models have taken into account of the intrinsic noise by either applying gillespie’s stochastic simulation algorithm to the ode models or constructing monte carlo models from first principles. with regard to modeling the impact of extrinsic noise on apoptosis pathway, there have been a few studies notably only on receptor-mediated apoptosis pathway. in this study, we attempt to develop mathematical and computational models of the intrinsic apoptosis pathway at single-cell level, and to identify plausible sources of non-genetic heterogeneity of apoptosis dynamics observed in a cell population using stochastic simulations. we start with a deterministic ode model of intrinsic apoptosis pathway adapted from the fussenegger model and find that bistability is missing. by adding positive feedback regulations that are supported by previous experimental evidences, we develop a model of intrinsic apoptosis pathway that functions as a bistable switch. we are particularly interested in understanding the stochastic behavior of this apoptosis switch under perturbation of intrinsic noise and/or extrinsic noise. stochastic modeling and simulations of intrinsic apoptosis pathway indicate that noise could enhance robustness of the bistable switch. in addition, we show that intrinsic noise is not sufficient to induce the observed level of cell-to-cell variability of apoptosis response at biologically relevant level of molecular numbers, while the extrinsic noise of protein variations is plausibly the main source giving rise to the degree of heterogeneous responses of intrinsic apoptosis pathway between single cells.  deterministic model of intrinsic apoptosis pathway to build a single-cell level deterministic model of intrinsic apoptosis pathway, we adopt the ode modeling framework of intrinsic apoptosis proposed by fussenegger et al. as it has the “minimal” model complexity while preserving all the critical interactions of intrinsic apoptosis pathway. the model scheme is shown in figure 1b and summarized here: the intrinsic apoptosis pathway is initiated by the release of cytochrome c from mitochondria to cytosol, and binding to apaf- to form the apoptosome complex. the apoptosome catalyzes procaspase-, the precursor of initiator caspase, to its active form, caspase. the executioner procaspase- is then activated by c9a to form the active executioner caspase, whose response represents the onset of the irreversible apoptosis fate. finally, cea is subject to the direct inhibition by iap protein. simulations of this ode model demonstrate that the time trajectories of cea starting from an initially low concentration first undergo some time delay and then switch to a high steady state with a relatively sharp transition. nevertheless, the time trajectories of cea from different initial conditions all converge to the same steady-state level, indicating the existence of single equilibrium point. indeed, bifurcation analysis of this ode model with respect to the parameter of cc concentration shows that the steady-state of cea is monostable with a sigmoidal input-output relationship. this means that regardless of different input strengths and different initial conditions, the response would asymptotically settle at only one steady state. further analytical analysis of the steady-state response of the fussenegger model proves that the mapping between the input of cc signal and the output of cea is either one-to-one or one-to-two and thus there exist at most two solutions of the output signal at equilibrium. as a result, the fussenegger model of intrinsic apoptosis pathway cannot be bistable as bistability of a system requires three distinct steady-state solutions in some range of input. the fussenegger model of intrinsic apoptosis pathway presents a system property of so-called ultrasensitivity rather than bistability. this, however, is inconsistent with the current consensus of bistability feature of apoptosis system, and thus extra model refining step is needed to resolve the discrepancy. to this end, we modify the existing model by incorporating a positive feedback from the executioner caspase to the initiator caspase, an auto-catalysis loop to the activation kinetic of c9a, and a mild cooperativity in cea activation induced by c3p, which are supported by previous experimental and computational findings. the existing mathematical models accounting for intrinsic apoptosis pathway differ from our proposed model in choice of apoptosis pathway. specifically, we follow the fussenegger model to assume that the inhibition of caspase by iap is not essential and ignored in our model. this choice is the main distinction from the model proposed by legewie et al, where iap inhibits both caspase and caspase thus giving rise to an iap-mediated positive feedback and bistability. the harrington et al model, the zhang et al model, and the kutumova et al model all inherit the schematic of the intrinsic apoptosis model proposed by legewie et al, and thus also include the double inhibitions of caspase and caspase by iap and the resulting implicit positive feedback. the governing odes of the modified model of intrinsic apoptosis pathway are listed in the methods section. simulations of the modified ode model show that when given a step input of low concentration of cc, the time courses of cea gradually settle at a near-zero steady state starting from different initial cea concentrations, while given a relatively high concentration of cc, cea eventually settles at a high steady state. such behavior with two stable output steady states indicates that bistability is achieved by the modified ode model. in addition, the time trajectories agree with experimental results in that the cea response is not elicited until after a few hours of delay time.the switching-on kinetic of cea activation is sigmoidal shape and completed within ∼1hr, presenting all-or-none switch-like behavior. indeed, one-parameter bifurcation analysis of the modified ode model confirms that the steady-state response of cea is bistable with respect to the input signal of cc, where two stable steady states coexist between the input concentrations of μm and μm. the bifurcation curve has two saddle-node bifurcation points sn and sn giving rise to a middle unstable branch and two stable branches, where the upper and lower branches correspond to the apoptosis and survival fates, respectively. when the concentration of cc approaches sn and sn hysteretic behavior occurs: the system remains at near-zero cea activity at low amount of cc, until an on threshold is reached, whereby cea activity switches to the apoptosis state abruptly; inversely, the cea activity switches from the apoptosis state to the survival state only if the cc concentration falls below the off threshold. the system properties of bistability and hysteresis could confer robustness to the apoptotic responsiveness by allowing cells that are not committed to apoptosis to remain at survival state, even in the event of mildly noisy input. in addition, two-dimensional bifurcation analysis with respect to four selected parameters, namely the hill constant that regulates the cea-mediated positive feedback loop, the cooperativity of activation of cea, the degradation rate of c9a and the inhibition rate of cea by iap, show that the bistability property of the modified model exists in extended parameter space around the nominal parameter set and is hence a robust behavior. stochastic model of intrinsic apoptosis pathway under intrinsic noise perturbation the deterministic model of the intrinsic apoptosis pathway has allowed us to analyze nonlinear properties of the system and quantify region of robust behavior. nevertheless, the ode modeling approach accounts for average cellular dynamics while ignoring the inevitable unpredictability embedded in biomolecular reactions and in intra- and extra-cellular environments. it has been observed by several different experimental techniques that the apoptosis response at single-cell resolution is subject to inherent stochastic perturbations, giving rise to pronounced cell-to-cell variability even in genetically identical cell population. therefore, it is necessary for us to develop a stochastic modeling framework of intrinsic apoptosis pathway, which can be used to explore the plausible origin of the stochasticity underlying the phenotypic heterogeneity. in this study, we investigate the impact of both intrinsic-noise and extrinsic-noise perturbations on the intrinsic apoptosis dynamics. first, to model the stochasticity under the intrinsic-noise perturbation due to low copy number of biomolecules, we assume that the deterministic model represents the nominal single-cell behavior, and apply the standard gillespie ssa. the detail of applying ssa to the deterministic model can be found in the methods section. to reflect the plausible low molecular copy numbers arising in the intrinsic apoptosis pathway within single cells, we perform simulations under different combinations of number of molecules that fall into the experimentally observed ranges. furthermore, in choosing number of molecules used for stochastic simulation of intrinsic noise, we limit the number to be within in general, as previous experimental and computational work have suggested that the coefficient of variation of cellular signal due to intrinsic noise reaches a minimal level when the mean abundance of molecules increases to. using this limit and previous experimental observation that procaspase- and procaspase- could have relatively low number of molecules per cell, ranging respectively from 5× to × and from undetectable to ×, we choose to run ssa of the model at molecule numbers of c9p = {5× } and c3p = {, }, which correspond to a total of eight combined conditions. note that each realization of stochastic simulation corresponds to the apoptosis response of one single cell. figure 4a shows stochastic time trajectories of output cea response with =, representing the behavior of cells, under varying amount of cytochrome c input. each time course of cea activation exhibits sigmoidal switch-like shape, converging to an elevated steady state. to illustrate the cellular variability in the stochastic response due to intrinsic noise, histograms of the steady-state cea copy numbers in the same cells is plotted in figure 4b against different input cc levels. the histogram shows that if the cc molecule is above copy number the distribution of steady-state cea is bimodal, with a low mean steady-state value of ∼ number of cea molecules and a high mean steady-state value of ∼ cea molecules. such bimodal distribution of cea response indicates that the stochastic response of intrinsic apoptosis pathway subject to intrinsic noise under this particular c9p/c3p condition is bistable. the bistability behavior persists even when the copy number of cc increases to ∼ molecules, showing that the corresponding fold change of bistability domain under intrinsic noise is above five times that of the deterministic model, where the bistability region of cc is as shown in figure 2d. such phenomenon of enhanced robustness induced by intrinsic noise is in agreement with previous computational work which suggests that stochastic signaling networks may perform more robustly than their deterministic counterpart. it is noteworthy that the existence and range of bistable cea response are dependent on the abundances of c9p and c3p in that the bistability boundary shrinks as the copy numbers of procaspases increase, and the bistable range becomes almost undetectable at significantly high amount of procaspases. this result implies that the bistability property of apoptosis might only be observable at proper condition of molecular abundance. since figure 4b has demonstrated that the apoptosis state of cea varies from cell to cell with a normal distribution, next we evaluate the degree of stochasticity of the cea response at equilibrium caused by intrinsic noise. we use the coefficient of variation to quantify the level of stochasticity, following previous studies of isogenic cell-to-cell variability of apoptosis response. as shown in figure 5a and 5b, under all the eight conditions of molecule numbers as described above, the mean value of the stochastic cea response increases while its standard deviation decreases as the input signal cc increases, both with exponential pattern. note that the condition of higher c9p and c3p abundances give rise to elevated mean response of cea at equilibrium as well as sharper stimulus-response curve, agreeing with the pro-apoptotic roles played by these two components. the consequent stochasticity of the output cea value, quantified by cv, monotonically decreases with increasing amount of input and reaches a minimum level as the number of cc molecules is above for all the eight conditions. this indicates that the intrinsic noise in the steady-state behavior of intrinsic apoptosis response at single cell level can be suppressed by enhancing the strength of input signal. note also that lower abundances of c9p and c3p induce more prominent stochasticity in output when the input amount of cc is sufficiently low. we then focus on an important quantity that has been experimentally measured to evaluate the cell-to-cell variability in apoptosis response. as demonstrated in figure 6a, each stochastic time trajectory of cea presents a time delay between time zero, the cc release time, and the cea activation time, which is conventionally defined as the half-maximal caspase cleavage time. such delay time is defined analogously to the quantity previously used in experiments to characterize the response time of extrinsic apoptosis pathway in single cells. and experimental quantifications suggest that td varies from cell to cell, even in genetically identical clones. our single-cell stochastic simulations with intrinsic noise, under the eight conditions of c9p and c3p abundances, show that in general the mean value of td monotonically decreases with increasing cc level. for instance, when c9p and c3p are respectively and molecules per cell, the mean value of td decreases from around min at number of cc molecules to ∼ min at number of cc molecules. understandably, lower abundances of pro-apoptotic proteins c9p and c3p lead to longer response delay time. we then quantify the degree of cell-to-cell stochasticity by cv of td. figure 6c shows that the cv of td due to intrinsic noise is, under all the conditions, a monotonically decreasing function with regard to the number of cc molecules. such findings of attenuated cellular variability of response under the circumstance of elevated input signal are consistent with previous experimental and computational results on the stochastic apoptosis pathway in genetically identical cells as well as other cellular processes. the above stochastic simulations show that the cell-to-cell variability of intrinsic apoptosis due to intrinsic noise is especially evident when the number of c3p is relatively small, where the cv of td persist at the level of non-genetic noise of apoptosis observed by experiments. moreover, the abundances of c9p and c3p inversely regulate the cellular variability among the eight conditions. specifically, the lowest c9p/c3p copy-number condition leads to a cv slightly above while the highest c9p/c3p copy-number condition leads to a cv slightly below, suggesting that the wide copy-number range of procaspases, especially procaspase- under the sub- condition is likely a source of cell-to-cell variability under the perturbation of intrinsic noise. previously, there has been study of apoptosis pathway implying that the molecular numbers of participating biochemical species seem to reside in a regime much higher than. therefore, our stochastic model with intrinsic noise is further explored under the various copy numbers of c9p and c3p with sufficiently strong cc input. as shown in figure 6d, under all the eight conditions the cv of td drops dramatically when the number of cc molecules increases above in particular, when the copy numbers of c9p, c3p and cc molecules are all raised to and above, the cv of td tends to an almost negligible level of ∼. such low degree of cv of td at physiologically more plausible condition indicates that the perturbation by intrinsic noise alone seems insufficient to induce the observed degree of cell-to-cell stochasticity of apoptosis dynamics, and other sources of uncertainty needs to be taken into account. stochastic model of intrinsic apoptosis pathway under extrinsic noise perturbation recently, natural protein variations across cell population has been identified as a major source of extrinsic fluctuations for apoptosis pathway. experiments have suggested that the concentration of a protein naturally varies among different cells following a log-normal distribution with a typical cv value of to. to address the impact of this kind of extrinsic noise on intrinsic apoptosis pathway, we assume that in the deterministic single-cell model of intrinsic apoptosis pathway described above certain protein of interest has log-normally distributed concentration with cv equal to. first, we simulate the stochasticity due to variation in the concentration of the protein cc, the critical input signal of intrinsic apoptosis pathway, using the single-cell deterministic model with a randomly generated cc concentration. again, as a measure of cell-to-cell variability, the cv of td among cells is calculated at varying mean values of cc protein. the variability peaks at mean concentration of cc = μm, and then declines as the mean abundance of cc protein increases, settling at a value of cv = ∼ eventually. it indicates that the natural variation in the cc protein can induce some degree of cell-to-cell stochasticity to the intrinsic apoptosis response, especially at relatively low abundance of cc protein. second, we evaluate the impact of extrinsic noise driven by the natural variation of iap protein, which is one of the most pivotal anti-apoptotic proteins tightly regulating apoptosis via antagonizing the activity of cea. the single-cell deterministic model plus a log-normally distributed amount of iap with cv equal to is simulated for cells at different concentrations of input cc protein. here, we plot the resultant cv of td evaluated among cells in figure 8a as a 2-dimensional heat map, versus the concentration of cc protein and the mean concentration of iap protein. it shows that given the same cc input, higher mean concentration of iap, albeit with the same dispersion of random distribution, leads to relatively larger cv of td and hence higher degree of cell-to-cell variability in the timing of cea activation. most likely, large amount of iap significantly represses the activity of cea, while it is generally known that the lower amount of a biochemical component the higher degree of its stochasticity. additionally, figure 8a indicates that at the same mean abundance of stochastic iap protein, higher input level of cc protein yields lower cell-to-cell variability in the apoptosis response. the above findings together imply that the extrinsic noise in cc protein has opposite influence on the induction of cell-to-cell variability to intrinsic apoptosis response than the extrinsic noise in iap protein. such contrary effect is in line with their opposite regulatory roles in apoptosis pathway, in that cc is pro-apoptotic while iap is anti-apoptotic. notably, figure 8a also shows that if the mean concentration of iap is moderately high, the cv of td can achieve a value of or even higher, which is the level of cell-to-cell stochasticity measured by experiments. hence the extrinsic noise in moderate amount of iap protein may serve as a significant source of cell-to-cell variability in intrinsic apoptosis response. but if iap is so high that paramount of inhibition is produced, the apoptosis response of cea is then constantly switched off, denoted by the upper-left non-responsive region in the 2d heat map. it is interesting to find out that when the 2d heat map is superimposed with the boundary of the two-parameter bistability diagram of the stochastic model, which is estimated using the histograms of cea response, the cv of td with relatively high values is located between the low-threshold curve and the high-threshold curve. such behavior is confirmed under systematic perturbations of the parameter kc and the other three parameters selected for the sensitivity analysis in figure . compared to the threshold lines for the bistability domain of the corresponding deterministic model, the extrinsic noise in cc seems not significantly affect the area of bistability domain, albeit the thresholds of cc shift toward smaller values. this trend also holds under different values of kc. as a further exploration of the impact of extrinsic noise, we allow the concentrations of both the cc and iap proteins to be log-normally distributed random numbers with cv equal to. the resulting 2d heat map of the cv of td in figure 8b shows that the level of stochasticity across cells is in general higher than the case under individual extrinsic perturbation of cc protein or iap protein. in most of the area in the 2-parameter region of figure 8b, the cv of td achieves a value above. similar to figure 8a, the cv of td with double extrinsic noises attains highest values between the low-threshold and high-threshold curves estimated for the stochastic bistability diagram, which is almost twice the cv value in the scenario under the extrinsic noise in iap protein only. moreover, the activation of cea response can now be elicited even in the non-responsive region of figure 8a due to the additional degree of fluctuation in the input signal. collectively, the 2d heat-map in figure 8b suggests that the extrinsic noises in cc and iap proteins are sufficient to yield the experimentally observed degree of cell-to-cell variability in the apoptosis response. the bistability domains under the nominal value of kc in figure 8b and those under systematically perturbed values of kc in additional file 4: figure s demonstrate the same trend of behavior as those in figure 8a and additional file 3: figure s close comparison between figure 8a/additional file 3: figure s and figures 8b/additional file 4: figure s indicate that the area of bistability domain on average is slightly larger in the latter group. therefore, additional degree of extrinsic noise may induce extra robustness for bistability of apoptosis. stochastic model of intrinsic apoptosis pathway under combined intrinsic and extrinsic noise perturbations we have so far analyzed the model of intrinsic apoptosis pathway subject to either intrinsic noise or extrinsic noise independently. additionally, we would like to find out if both types of noises are present in the intrinsic apoptosis pathway, how the cell-to-cell variability of the delay time of cea response is influenced. to simulate a model under the perturbation of combined intrinsic and extrinsic noises, we implement the stochastic simulation algorithm of the apoptosis model as described above to mimic the intrinsic noise, and simultaneously allow certain protein concentrations to be log-normally distributed random variables to represent the extrinsic noise. in order to make comparison with the results of the intrinsic-noise only case, we again use the eight conditions of molecule numbers at c9p = {5× } in combination with c3p = {, }. first, the impact of extrinsic fluctuation in the amount of cc protein on top of the intrinsic noise is simulated. the resulting cv of td in response to the abundance of cc less than is shown in the inset of figure 7b. note that the unit of cc protein is now copy number of molecules rather than μm as a requirement by the implementation of ssa. comparing it to figure 6c, we find that the behavior of cell-to-cell variability due to the combined types of noises is almost the same as that under intrinsic noise alone. that is, under all the eight molecule-number conditions, the cv of td monotonically decreases in a near exponentially-decaying fashion as the number of cc molecules increases, and it falls into similar range of value. and similar to figure 6c, the variability value is inversely correlated with the abundances of c9p/c3p, achieving cv of td above if c9p/c3p have sub- copy numbers. this result indicates that at the limit of low number of molecules, the intrinsic noise seems to make the dominating contribution to the cell-to-cell stochasticity of intrinsic apoptosis response. we subsequently run the same stochastic model under the perturbation of intrinsic noise plus extrinsic noise in cc, while increasing the mean abundance of cc molecules up to, a more plausible reacting scale for apoptosis pathway. we find that the stochasticity curves under different c9p/c3p abundance conditions relatively converge at cc equal to copy number, due to diminished intrinsic noise. the comparison between figure 7b and figure 6d enables us to see the relative contributions of intrinsic and extrinsic noises: except for the case of almost undetectable abundance of c3p, which induce severe intrinsic noise, the cv of td due to combined noises is significantly higher than that due to intrinsic-only noise. specifically, the contrast between the case of combined noises and the case of intrinsic-only noise is highest when the copy numbers of c9p, c3p and cc molecules are all above therefore, at higher molecular numbers, the cellular stochasticity seems to be majorly contributed by the source of extrinsic noise. with respect to the possible impact of the extrinsic noise in iap, we simulate the stochastic model with additional perturbation in the amount of iap protein. that is, both the abundances of iap and cc proteins are assumed to be log-normally distributed around their respective mean value with a cv of, alongside with intrinsic noise. at the limit of low cc abundance the additional extrinsic noise of iap does not introduce extra cell-to-cell stochasticity to apoptosis response. when the mean number of cc molecules is in the physiologically plausible range of to, simulations show that td presents a high cv, mostly above the value of and even going up to a value of in the 2d heat map. we compare the cell-to-cell variability under combined noises with that under extrinsic-only noise, and find that the cv of td increases only slightly in the former case. this is consistent with the minor role of intrinsic noise played at elevated abundances of molecules. as a consequence, the boosted degree of stochasticity in figure  has to owe predominantly to the extrinsic fluctuations in the amount of cc and iap proteins rather than the intrinsic noise. simulations under the eight conditions of c9p and c3p abundances demonstrate that the overall cell-to-cell variability is, the same as what we see earlier, inversely correlated with the copy numbers of molecules, where the cv of td on average declines through row to row as the abundance of c3p increases from to while the cv of td on average slightly declines from column to column as the abundance of c9p increases from 5× to this is consistent with the aforementioned observation of the evident variability induced by the procaspase- under sub- copy-number condition through the channel of intrinsic noise. it is noteworthy that the cv of td in figure  is generally larger at higher mean value of iap molecular number, agreeing with the inference from figure 8a that the abundance of iap protein is capable of promoting the stochasticity of apoptosis pathway. it highlights the interesting role of iap protein in controlling the cell-to-cell variability of intrinsic apoptosis response, and suggests that in treating diseases exploiting the apoptosis mechanism, such as cancer, the iap protein offers a potential therapeutic target not only for effective modulation of apoptosis but also for eliminating the undesired cell-to-cell heterogeneity, a major obstacle to effective cancer treatment and personalized medicine. this finding is in line with the current view of the therapeutic function of iap based on the study of apoptosis pathway at cell population level.  the recently observed heterogeneous apoptosis phenotypes at single cell level have drawn increasing attention from researchers. mathematical modeling and computer simulation provide an efficient approach to gain deep insight into the dynamic behavior of apoptosis network. this paper develops a theoretical and computational framework for single-cell stochastic modeling of the intrinsic apoptosis pathway. using this modeling framework we explore the stochastic behavior of the intrinsic apoptosis response at single-cell level and seek to understand the plausible sources underlying the experimentally observed cell-to-cell variability of apoptosis response. we show that in the presence of noise, the bistable response of intrinsic apoptosis pathway can be more robust than its deterministic behavior. the coefficient of variation of the delayed timing of the activity of executioner caspase is utilized to quantify the stochasticity in the apoptosis dynamics. we find that the intrinsic noise can introduce significant cell-to-cell variability if the abundances of reacting biomolecules are relatively low. the level of cellular stochasticity solely due to intrinsic noise decreases dramatically to a negligible level of cv equal to ∼ when the copy number of cytochrome c is raised to, which is the amount suggested by a previous study. in addition, the extrinsic noise caused by the natural variations in protein concentrations of two key components in the intrinsic apoptosis pathway, cytochrome c and the inhibitor of apoptosis proteins, is also accounted for without or with intrinsic noise. series of simulations indicate that the extrinsic noise is plausibly the major source of the cell-to-cell variability of intrinsic apoptosis response at high number of biomolecules. furthermore, we find that the mean abundance of fluctuating iap is positively correlated with the degree of cell-to-cell variability, thus making iap a potential target for therapeutically suppressing the stochasticity of intrinsic apoptosis response across cell population in treating diseases such as cancer. in summary, this study based on our theoretical and computational models characterizes the behavior of the intrinsic apoptosis pathway under complex stochastic perturbations, and suggests that certain deterministic features, such as system bistability and iap as potential therapy target, still remain in the presence of noise. altogether, the work can enable us to gain deeper understanding toward the experimentally observed uncertainty in cellular decision making. 
our deterministic model, adapted from the well-accepted fussenegger model, shows that an additional positive feedback between the executioner caspase and the initiator caspase plays a fundamental role in yielding the desired property of bistability. histograms of the steady-state output at varying input levels show that the intrinsic noise could elicit a wider region of bistability over that of the deterministic model. cellular networks, however, are inherently stochastic and significant cell-to-cell variability in apoptosis response has been observed at single cell level. to address the inevitable randomness in the intrinsic apoptosis mechanism, we develop a theoretical and computational modeling framework of intrinsic apoptosis pathway at single-cell level, accounting for both deterministic and stochastic behavior. furthermore, the extrinsic noise represented by random variations of two key apoptotic proteins, namely cytochrome c and inhibitor of apoptosis proteins, is modeled separately or in combination with intrinsic noise. however, the system stochasticity due to intrinsic fluctuations, such as the noise of steady-state response and the randomness of response delay, shows that the intrinsic noise in general is insufficient to produce significant cell-to-cell variations at physiologically relevant level of molecular numbers. we then examine the impact of intrinsic fluctuations of biochemical reactions, viewed as intrinsic noise, and natural variation of protein concentrations, viewed as extrinsic noise, on behavior of the intrinsic apoptosis network. dependent on the type of stimulus, apoptosis can be propagated by extrinsic pathway or intrinsic pathway. the comprehensive understanding of the molecular mechanism of apoptotic signaling allows for development of mathematical models, aiming to elucidate dynamical and systems properties of apoptotic signaling networks. finally, simulations illustrate that the mean abundance of fluctuating iap protein is positively correlated with the degree of cellular stochasticity of the intrinsic apoptosis pathway. in addition, it predicts that the iap protein could serve as a potential therapeutic target for suppression of the cell-to-cell variation in the intrinsic apoptosis responsiveness. the resultant stochasticity in the timing of intrinsic apoptosis response shows that the fluctuating protein variations can induce cell-to-cell stochastic variability at a quantitative level agreeing with experiments. our theoretical and computational study shows that the pronounced non-genetic heterogeneity in intrinsic apoptosis responses among individual cells plausibly arises from extrinsic rather than intrinsic origin of fluctuations. intrinsic apoptosis pathwaystochastic modelintrinsic noiseextrinsic noise apoptosis is a cell suicide mechanism that enables multicellular organisms to maintain homeostasis and to eliminate individual cells that threaten the organism’s survival.
 apoptosis is a cell suicide mechanism that enables multicellular organisms to maintain homeostasis and to eliminate individual cells that threaten the organism’s survival. dependent on the type of stimulus, apoptosis can be propagated by extrinsic pathway or intrinsic pathway. the comprehensive understanding of the molecular mechanism of apoptotic signaling allows for development of mathematical models, aiming to elucidate dynamical and systems properties of apoptotic signaling networks. there have been extensive efforts in modeling deterministic apoptosis network accounting for average behavior of a population of cells. cellular networks, however, are inherently stochastic and significant cell-to-cell variability in apoptosis response has been observed at single cell level.  to address the inevitable randomness in the intrinsic apoptosis mechanism, we develop a theoretical and computational modeling framework of intrinsic apoptosis pathway at single-cell level, accounting for both deterministic and stochastic behavior. our deterministic model, adapted from the well-accepted fussenegger model, shows that an additional positive feedback between the executioner caspase and the initiator caspase plays a fundamental role in yielding the desired property of bistability. we then examine the impact of intrinsic fluctuations of biochemical reactions, viewed as intrinsic noise, and natural variation of protein concentrations, viewed as extrinsic noise, on behavior of the intrinsic apoptosis network. histograms of the steady-state output at varying input levels show that the intrinsic noise could elicit a wider region of bistability over that of the deterministic model. however, the system stochasticity due to intrinsic fluctuations, such as the noise of steady-state response and the randomness of response delay, shows that the intrinsic noise in general is insufficient to produce significant cell-to-cell variations at physiologically relevant level of molecular numbers. furthermore, the extrinsic noise represented by random variations of two key apoptotic proteins, namely cytochrome c and inhibitor of apoptosis proteins, is modeled separately or in combination with intrinsic noise. the resultant stochasticity in the timing of intrinsic apoptosis response shows that the fluctuating protein variations can induce cell-to-cell stochastic variability at a quantitative level agreeing with experiments. finally, simulations illustrate that the mean abundance of fluctuating iap protein is positively correlated with the degree of cellular stochasticity of the intrinsic apoptosis pathway.  our theoretical and computational study shows that the pronounced non-genetic heterogeneity in intrinsic apoptosis responses among individual cells plausibly arises from extrinsic rather than intrinsic origin of fluctuations. in addition, it predicts that the iap protein could serve as a potential therapeutic target for suppression of the cell-to-cell variation in the intrinsic apoptosis responsiveness. intrinsic apoptosis pathwaystochastic modelintrinsic noiseextrinsic noise 