Systems Analysis of Cancer Cell Heterogeneity in Caspase-dependent Apoptosis Subsequent to Mitochondrial Outer Membrane Permeabilization

Background: Computational systems models may allow understanding of progression and impairment of apoptosis based on quantitative protein profiles. Results: We found apoptosis to be impaired in some cell lines after MOMP, which can be understood by our systems model “APOPTO-CELL”. Conclusion: APOPTO-CELL helps to understand heterogeneity in apoptosis execution after MOMP. Significance: APOPTO-CELL may help to assess cancer-specific efficacy of chemotherapeutics that induce apoptosis.

To study the mitochondrial pathway of apoptosis subsequent to MOMP, we previously developed the computational systems model APOPTO-CELL (17,19).As input, APOPTO-CELL requires absolute levels of five proteins: the pro-apoptotic proteins APAF-1, Smac, procaspase-3, and procaspase-9 and the anti-apoptotic protein XIAP.Assuming that MOMP has occurred, APOPTO-CELL calculates the time profile of concentrations for all proteins.It further quantifies apoptosis execution kinetics (speed) by calculating the amount of cellular substrate that gets cleaved over time.
In this study, we found that apoptosis can be impaired in some cancer cell lines subsequent to MOMP and that this impairment can be explained using APOPTO-CELL.We further present APOPTO-CELL as a tool to study cell-specific molecular mechanisms during apoptosis execution that arise from cell line individual protein profiles.We finally use APOPTO-CELL to test the robustness of apoptosis execution against fluctuations in enzyme activities and provide a workflow for model optimization.
Transfection and Single Cell Microscopy-Caspase-3-dependent FRET disruption in DLD-1 WT, DLD-1 XIAP 0/Ϫ , HCT-116 WT, HCT-116 XIAP 0/Ϫ , MCF-7/C3, and HeLa single cells was obtained from previous publications (17,20,23).De novo measurements of apoptosis kinetics for LoVo, HT-29, and HCT-116 Smac Ϫ/Ϫ were performed using the method described previously (17).Briefly, cells were equilibrated with 30 nM TMRM in Hepes-buffered medium (15 mM, pH 7.4), covered with mineral oil, and placed in a heated (37 °C) incubation chamber that was mounted on the microscope stage.The membrane-permeable, cationic probe TMRM was used to measure the mitochondrial membrane potential.The onset of TMRM loss that indicates depolarization of the mitochondrial membrane potential was used as a marker for cytochrome c release and MOMP.
Apoptosis was induced with 3 M staurosporine (STS).To report caspase-3 activation, cells were transfected with plasmid DNA (pmyc-CFP-DEVD-YFP (24)).Cleavage kinetics of this substrate were detected at the single-cell level.CFP/FRET emission ratio traces were obtained by dividing the average fluorescence intensity values of single cells after background subtraction.A Zeiss LSM 710 META inverted microscope (Carl Zeiss) attached to a confocal laser-scanning unit equipped with a 405-nm diode laser, 488-nm argon laser, and 543-nm helium/neon laser was used.CFP, YFP, and TMRM fluorescence and FRET were monitored with a ϫ63, numerical aperture 1.4 oil immersion objective with the optical slice thickness set to 1.5 m (full width at half-maximum) and detected using optimized filter and mirror sets.Subsequently, TMRM average intensity and CFP/FRET emission ratio traces were obtained for fluorescent cells.Experiments were terminated 16 -24 h after STS administration.Onset of TMRM decrease and DEVD-FRET substrate cleavage were determined visually against the base-line signal.Image processing and analysis were performed with MetaMorph software (Molecular Devices, Wokingham, UK).
Experimental Apoptosis Execution Kinetics from Literature-FRET traces for substrate cleavage for HeLa and XIAP-overexpressing HeLa cells were obtained from Ref. 17. Traces for DLD-1, DLD-1 XIAP 0/Ϫ , HCT-116, and HCT-116 XIAP 0/Ϫ cells were remodeled from data specifying onset and duration in (20).For MCF-7/C3 cells, data were taken from Refs. 17 and 22. FRET traces were represented by a sigmoid Boltzmann function, Here dt is the width, A min and A max are minimum and maximum, and t 0 is the point of inflection of the FRET signal.
Statistics-Statistical analysis was performed using MATLAB (The Mathworks Inc., Natick, MA).Protein concentrations are given as means Ϯ S.E.Single cell data were analyzed by analysis of variance and Student's t test when normally distributed and given as mean Ϯ S.E.When normality could not be assumed, Kruskal-Wallis and Mann-Whitney test were used, and data were given as median with interquartile range.p values of Ͻ0.05 were considered significant.
Computational Model of Mitochondrial Apoptosis-The original model was described previously (see main text and supplemental material of Ref. 17).Briefly, MOMP was assumed to be invoked, eventually culminating in activation of caspase-3 and cleavage of cellular substrate (model output).Quantified protein levels of APAF-1, XIAP, procaspase-3, procaspase-9, and Smac were used as model input.Each molecular interaction was modeled by mass action kinetics and transformed into a set of ordinary differential equations.Ordinary differential equations were solved by MATLAB (The MathWorks Inc.).Mechanistic details can be found in supplemental Table 1.The model code of APOPTO-CELL can be obtained from the authors by request.
It was assumed that the S.E.value of quantified proteins from three Western blot experiments partly reflects information on cell-to-cell heterogeneity in a particular cell line (each experiment is considered as a slightly biased sample) and is partly a sign of variations in the experimental handling.As a consequence, the S.E. was assumed to contain biological information as well as experimental uncertainties, both representing valuable information for the analysis of apoptosis in cell populations.Hence, three scenarios were considered (Table 1).The first scenario ("average scenario") used the average concentrations of three experiments as input.For the population-based analysis of Fig. 2, an additional second and third scenario not only took into account the average protein concentrations for each of the above proteins and for each cell but also the S.E.over three experiments.For the second scenario (defined as the "slow scenario"), the concentration of the pro-apoptotic proteins procaspase-3, procaspase-9, Smac, and Apaf-1 was set to the aver-age concentration minus one S.E., and the concentration of the anti-apoptotic protein XIAP was set to average concentration plus one S.E.For the third scenario ("fast scenario"), the concentration of the pro-apoptotic proteins was set to the average concentration plus one S.E., and the concentration for XIAP was set to the average concentration minus one S.E.For the optimization workflow (Figs. 5 and 6), the S.D. was used instead of the S.E. for the slow and fast scenario, in order not to overly restrict the criteria for parameter screening.
Methods of Parameter Assessment-Those parameters were selected whose changes would most likely have the highest impact on the model output.Two partially complementary methods were employed to select the most influential parameter sets and are described below.In both methods, the values of the kinetic model parameters ͗k j 0 ͘ were changed from the reference set k 0 ϭ ͗k j 0 ͘ published in Ref. 17 to a set k ϭ ͗k j ͘ with changed parameters.It was then assessed how these changes of parameter values affected the kinetics of the predicted substrate cleavage.The assessment was done either for individual parameters or for parameter combinations, depending on the method used.
Individual Parameter Assessment (Parameter Set "I")-In the first assessment method (individual parameter assessment), each parameter was multiplied separately by a series of factors 8, 4, 2, 1 ⁄ 2, 1 ⁄ 4, and 1 ⁄ 8. Subsequently, changes of the time point at which the model predicted 80% cleavage of the cell substrate were calculated (denoted as t 80 ).Whenever a parameter change resulted in less than 80% substrate cleavage after 300 min, the time was censored at 300 min.Parameters were ranked, giving the highest rank to that parameter where the time point t 80 changed most decisively between an 8-fold up-and down-regulation.The procedure was performed for HeLa, HCT-116, and DLD-1 cells separately assuming the "average scenario" (Fig. 4, A-C).
For studying the effects of parameter variations in all cells and scenarios ("slow scenario," "average scenario," and "fast scenario"; Table 1) simultaneously, the shifts in the predicted t 80 values between the reference APOPTO-CELL and the model assuming parameter variations were calculated.Values were subsequently averaged over all cells and all scenarios and ranked according to their change between 8-fold up-and down-regulation with the highest rank given to the parameter with highest influence.The 10 parameters with the highest rank were selected as parameter set I ("individual"; Table 2).
Assessment of Parameter Cooperation (Parameter Set "C")-In a second method, the effect on model output of simultaneous but local changes was assessed using a procedure described

Parameter combination for fast, average, and slow apoptosis execution scenarios based on protein variability
Scenarios were defined to take into account variability in experimentally detected protein levels within a population (Fig. 2, A-E).A "slow" scenario assumes that the concentration of the anti-apoptotic protein XIAP is one S.E.value higher than the average of three experiments, whereas all levels of pro-apoptotic proteins are one S.D. value lower.An "average scenario" assumes the concentrations to be the average concentrations that were measured for the respective protein.A "fast scenario" assumes XIAP concentration to be one S.E.value lower than the average concentration detected for XIAP and all other proteins one S.E.value higher than the measured concentrations.To extend the range of acceptable parameter sets, the S.D. was used instead of the S.E. for the parameter screening in the optimization workflow (Figs. 5 and 6).
previously (25).Briefly, pairs of two parameters were varied, and changes in model predictions (amount of substrate cleavage over time) were assessed.Results of pairwise variations were assigned as matrix elements, and principal components of this matrix were calculated.The principal components associated with the highest eigenvalues were regarded as parameter combinations with the highest influence on the model output.Only the first principal component was evaluated due to its dominant influence on substrate cleavage (85% explained variance where all cells and scenarios were pooled).Parameters were ranked according to the highest influence on this first principal component.The 10 parameters with highest rank were defined as parameter set C ("cooperative"; Table 2).
Workflow to Assess the Robustness of Apoptosis against Parameter Variations-A workflow for testing the robustness of apoptosis execution kinetics against variations of kinetic parameters was developed (Fig. 5A).The workflow employed parameter sets I and C, which were identified as those parameters whose variations had, either alone or in combination, the most decisive influence on the model output (caspase-dependent DEVD-FRET substrate cleavage over time) for HCT-116, DLD-1, and HeLa cells in combination.Values of these parameters were assumed to vary within a biologically realistic range, obtained from concurring values described in literature (supplemental Table 2).Several sets of parameter values were assigned to each parameter set.Values were taken from the above realistic range using an equidistant set (brute force screening) or randomly drawn from a uniform distribution (Monte Carlo screening).
A set of parameter values was assumed to be realistic when the APOPTO-CELL prediction of substrate cleavage over time fulfilled two experimentally motivated filter criteria.A first filter criterion (criterion "PX" for "protein expression"; Fig. 5B) was assumed to reflect the experimental fact that HeLa cells with about 4.5-fold expression of XIAP levels (HeLa XIAP Adv cells) did not show full cleavage of their substrate upon apoptosis, whereas HeLa wild type cells did (17).Therefore, a series of calculations for HeLa cells was performed with varied initial levels of XIAP.It was required that the model predict at least 95% of the cell substrate after 300 min for cells with less than 1.6-fold overexpression (0.1 M) of XIAP.Likewise, it was required that more than 50% of cell substrate remain present when XIAP was expressed 7.9 times or higher (0.5 M).Parameter value sets where the predicted substrate cleavage did not fulfill this criterion were immediately discarded.
A second filter criterion (criterion "S80," for "time to 80% substrate cleavage"; Fig. 5C) was employed to check whether the predicted apoptosis execution kinetics for DLD-1, HCT-116, and HeLa cells were in a range reasonably close to the experimentally detected values.Uncertainties of predictions imposed by experimental variations in input proteins were further taken into account.Fulfillment of the criterion required that the time t 80 fast (time point of 80% substrate cleavage for the "fast scenario") was predicted to be earlier and the time t 80 slow for the "slow scenario" to be later than the time where 80% FRET substrate was observed to be cleaved experimentally (Table 1).The compliance to this criterion was first checked for all investigated cells (DLD-1, HCT-116, and HeLa) individually.Subsequently, it was checked what parameter value sets fulfilled this criterion for all combinations of two or for all three cell types.
Brute force screening was performed on the five most influential parameters of parameter set C. For each of those five parameters, 10 equidistant values were taken over the entire parameter range and permuted, resulting in 100,000 parameter value sets and, thus, as many test calculations.
Monte Carlo screening was performed for the five and ten most influential parameters of parameter set I and parameter set C. Therefore, stochastically selected parameter values were used and taken from a uniform distribution of the five-dimensional or 10-dimensional parameter space (Table 2).For Monte Carlo simulations with five parameters, calculations with either 1050 or 100,000 parameter value combinations were performed to confirm that Monte Carlo screening with less density of even a 1050-parameter value set would lead to acceptable results.For Monte Carlo screening of 10 parameters, 1,050,000 parameter value sets were used.
Assessment of Screening Results to Quantify Alternative APOPTO-CELL Parameterizations-For each screening attempt (i.e.brute force screening with five parameters and Monte Carlo screening with five and 10 parameters), sets of parameter values were obtained such that they fulfilled both filter criteria defined above.These resulting sets of parameter values were assessed according to how well the predicted substrate cleavage (here denoted as calculated FRET trace, FRET calc ) agreed with experimental FRET traces (FRET exp ).For this assessment, the following goodness of fit function was used, which was adjusted for the degrees of freedom,

TABLE 2 List of selected kinetic parameters of the model
Parameters identified as most decisive regulators of model output (predicted FRET substrate cleavage over time) by individual assessment (parameter set I) and pairwise assessment (parameter set C). ID, parameter ID as defined in supplemental Table I.Description, short description of biochemical role of parameter.Set, parameter set in which this parameter is used.Further details of the mechanistic role of these parameters (as defined by model equations) can be found in supplemental Table I.In Equation 2, n denotes the number (n ϭ 101) of time points (t i ), and m is the number of screening parameters (m ϭ 5 or 10, respectively).All R adj 2 values belonging to a particular screening and to a particular cell line (HCT-116, DLD-1, and HeLa cells) were sorted in descending order.Then all R adj 2 values were scaled to the range from 1 (value set with best match) to 0 (value set with least match).For each screening, R adj 2 values for HeLa, HCT-116, and DLD-1 were added.The 10 sets of parameter values with the highest total values (R adj 2 over these three cells) were selected.With all different screening methods, different parameter sets, and different subsets of influential parameters, 60 sets of parameter values were used (see "Results" for the chosen combinations).

Execution of Caspase-dependent Apoptosis after MOMP Is
Impaired in HCT-116 Smac Ϫ/Ϫ and LoVo Cells-MOMP is often regarded as a point of no return in programmed cell death.However, evidence has been found that proteins involved in apoptotic processes subsequent to MOMP are also predictive markers of outcome of chemotherapy treatment (26 -29).These findings suggest that deregulation of these proteins may impact the efficiency of chemotherapeutically induced apoptosis.We therefore hypothesized that in some cell systems, apoptosis signaling subsequent to MOMP can be impaired.
Previous single cell microscopy experiments of our group demonstrated that HeLa cervical carcinoma cells, MCF-7 breast cancer cells expressing caspase-3, and the colorectal cancer cell lines DLD-1 and HCT-116 as well as their XIAP-deficient clonal lines DLD-1 XIAP 0/Ϫ and HCT-116 XIAP 0/Ϫ showed rapid cleavage of cellular substrate after MOMP (17,20,22,23).In contrast, HeLa cells that were infected to overexpress the anti-apoptotic protein XIAP (HeLa XIAP Adv ) showed slow effector caspase activation that led to submaximal cleavage of cellular substrate (17).We therefore sought further evidence for whether two types of responses in cancer cell lines after MOMP exist that can be associated with either rapid and complete or delayed and/or submaximal (i.e.impaired) cleavage of cellular substrate.
To investigate this, we first exposed HT-29 human colon adenocarcinoma cells to a 3 M concentration of the broad spectrum kinase inhibitor STS, a stimulus that induces cell death through the mitochondrial apoptosis pathway.From the cells analyzed, MOMP was observed in 35 of 42 cells (83%; n ϭ 2 independent experiments), as indicated by mitochondrial fluorescence reduction of the cationic dye TMRM.Moreover, HT-29 cells that underwent MOMP also showed rapid and pronounced caspase activation, as indicated by the cleavage of the caspase-sensitive DEVD-FRET substrate (see Fig. 1A for two typical cells).In particular, from the 35 cells showing MOMP, 31 showed onset of FRET disruption, indicating the beginning of caspase-induced substrate cleavage.Furthermore, 29 of those 35 cells showing MOMP also reached a plateau in the FRET response, indicating that this DEVD-FRET substrate cleavage was complete.Finally, onset of FRET disruption was already observed 5 Ϯ 3 min (mean Ϯ S.E.) after MOMP (Fig. 1, E and F), and complete cleavage was reached 26 min later (median, with interquartile range 14 -30 min).Due to the rapid kinetics of FRET substrate cleavage and the large amount of cells that underwent robust FRET disruption after MOMP, we denoted HT-29 and other cell lines showing these features as "post-MOMP responder cell lines." We next wondered whether the absence of a pro-apoptotic protein, such as Smac, which acts as enhancer of caspase activity, may impair caspase-dependent apoptosis.Using the same treatment protocol as in this study, we previously demonstrated that parental HCT-116 colorectal cancer cells underwent fairly rapid apoptosis (20).We therefore performed the same experiments in the Smac-deficient HCT-116 clone (HCT-116 Smac Ϫ/Ϫ ).As result, 34 of 41 analyzed cells (n ϭ 2 experiments) still underwent MOMP, yet fewer cells than in the HT-29 cell line that underwent MOMP also showed onset of FRET cleavage (76%; Fig. 1, B (for a typical cell) and D for all analyzed cells).Even fewer cells that underwent MOMP reached a plateau in the FRET signal (65%, Fig. 1D), which would have indicated complete cleavage of DEVD-FRET substrate during the time frame investigated here.Most strikingly, on average, onset (165 Ϯ 151 min) and duration (79 min; interquartile range 44 -122 min) of DEVD-FRET substrate cleavage was significantly delayed compared with HT-29 cells (Fig. 1F).As an aside, we note that the markedly slower apoptosis kinetics in HCT-116 Smac Ϫ/Ϫ as opposed to parental HCT-116 calls cannot be explained by differences in expressions of BCL-2 family proteins because both clones had similar expression levels of these proteins compared with the post-MOMP responder cell line HeLa (data not shown).To this end, BCL-2 proteins have been shown to be dispensable for processes after MOMP (1).
We finally sought a cell line that shows impairment of DEVD-FRET substrate cleavage but is not known to have any deficiencies or deregulations in proteins responsible for regulating apoptosis subsequent to MOMP.We speculated that LoVo adenocarcinoma cells may be such a cell line, supported by the fact that in-house flow cytometry proved them to be remarkably resistant to genotoxic stress. 4Indeed, repeating the above single cell microscopy experiments in LoVo cells demonstrated that from those cells that underwent MOMP (54 of 67 analyzed cells), less than half (25 of 54 cells) showed a detectable onset of DEVD-FRET substrate cleavage (Fig. 1, C (for typical cells with delayed cleavage and for such cells showing no DEVD-FRET cleavage) and D (for cell populations)), and even fewer cells (20 cells, 37%) reached a plateau, indicating that in the remaining 63% of cells, only little or no FRET substrate was cleaved.This was also corroborated by acceptor photobleaching, which demonstrated incomplete substrate cleavage in LoVo cells subsequent to MOMP (data not shown).Similar to HCT-116 Smac Ϫ/Ϫ cells, the onset of caspase activation was late (onset 209 Ϯ 196 min after MOMP; Fig. 1F), and substrate was cleaved slowly (80 min, interquartile range 26 -94 min; Fig. 1F).We concluded that, similar to the previously characterized XIAP-overexpressing HeLa XIAP Adv cells (17), HCT-116 Smac Ϫ/Ϫ and LoVo cells also showed a significantly delayed caspase activation, a high cell-to-cell variability in the kinetics of FRET cleavage, and often no FRET cleavage at all after MOMP.We therefore defined such cell lines with delayed and/or lack of FRET substrate cleavage as "post-MOMP-impaired responder cell lines." Systems Modeling Explains Impairment of Apoptosis Execution after MOMP-Delayed or blocked caspase activation can transform cells into a potentially malignant state or make them resistant to chemotherapeutically induced apoptosis (1, 2, 30).Indeed, once caspase activation is inhibited, cells have been shown to be able to recover from MOMP (15,(31)(32)(33).A correct classification as to whether or not a cell line is a potential post-MOMP-impaired responder cell line may therefore be a valuable parameter in assessing the success of chemotherapeutic treatments that induce the mitochondrial apoptosis pathway.We therefore investigated whether our model APOPTO-CELL was able to correctly distinguish between responder and impaired responder cell lines.Because APOPTO-CELL required absolute concentrations of the proteins APAF-1, Smac, procaspase-3/-9, and XIAP as input, we determined these concentrations by quantitative Western blots using previously determined values in HeLa cells as a standard (Fig. 2,  A-E, n ϭ 3 experiments, mean Ϯ S.E.).To allow us to model cell-to-cell variability within populations of each cell line, we further included the S.E. of the detected protein levels into our analysis.We constructed, for each cell line, three different sets of input protein data that are associated with different speeds and efficiency of apoptosis execution (Table 1; "slow scenario," "average scenario," and "fast scenario"; see "Experimental Pro- ).F-L, post-MOMP responder cell lines were classified by the model to have rapid kinetics of substrate cleavage in a typical cell of the population (dashed line), a rather synchronized response in cleavage over the population (smaller gray area), and a complete substrate cleavage in almost all cells (Ͼ80% after 300 min).Predictions aligned well with the fast and robust apoptosis execution kinetics for HT-29 cells in Fig. 1A and with previous measurements in the other cells where apoptosis execution was determined to require less than 80 min (17,20,23)).M-O, experimentally classified post-MOMP-impaired responder cell lines were correctly predicted to have slow kinetics of substrate cleavage in an average cell of a population (dashed line) and a low amount of substrate cleavage after 300 min, consistent with the reduced amount of cells showing robust DEVD-FRET substrate cleavage after MOMP.Kinetics of substrate cleavage were predicted to be highly heterogeneous over the population consistent with the findings in Fig. 1, E and F, and Ref. 17.
cedures").Using APOPTO-CELL, we calculated the amount of substrate cleavage that is cleaved over time.We assumed that any realistic prediction of substrate cleavage within a population of a particular cell line must lie in the range between the slow and fast scenario (gray areas in Fig. 2

, F-O, average scenario indicated by a dashed line).
By comparison with experimental data, we found that APOPTO-CELL correctly classified HT-29 cells as a post-MOMP responder cell line.This was indicated by the fact that the calculated kinetics of substrate cleavage, based on the average protein concentration for this cell line, were found to be fast (Fig. 2F, dashed line), that population-based variability was modest (small gray area), and that after 300 min, 80 -100% of the substrate was cleaved.By similar means, APOPTO-CELL correctly classified HCT-116 WT, HeLa WT, DLD-1 WT, DLD-1XIAP 0/Ϫ , HCT-116 XIAP 0/Ϫ , and MCF-7/C3 cells as responder cell lines (Fig. 2, H-L).This finding was consistent with previously published single cell data demonstrating that caspase activation was fast and robust in these cells (17,20,23).Further, these calculations confirmed experimental findings that lack of XIAP accelerated apoptosis execution (20).
In contrast, LoVo, HCT-116 Smac Ϫ/Ϫ , and HeLa XIAP Adv cells were correctly classified as post-MOMP-impaired responder cells.In particular, APOPTO-CELL predicted a high cell-to-cell variability in the apoptotic kinetics, as indicated by the size of the gray area in all of these cell lines (Fig. 2, M-O).This mirrored the experimental cell-to-cell variability in onset and duration of DEVD-FRET substrate cleavage (Fig. 1, E and F) (17).In further agreement with experimental data presented here and in Ref. 17, the DEVD-FRET substrate cleavage based on the average protein concentrations was correctly predicted to be slow and to be incomplete after 300 min (dashed lines in Fig. 2, M-O, substrate cleavage Ͻ10% for LoVo, ϳ50% for HCT-116 Smac Ϫ/Ϫ , and ϳ25% for HeLa XIAP Adv ).
After having successfully employed APOPTO-CELL to separate post-MOMP responder and non-responder cell lines, we finally wanted to investigate how well APOPTO-CELL resembled the onset and duration of FRET cleavage in the investigated cell lines.As can be seen in Fig. 3A, agreement between predicted and experimental onset was highly statistically signif-icant (p ϭ 0.002, ϭ 0.852 Pearson correlation).Predicted and experimentally determined duration of FRET cleavage (5-95%; only cell lines were shown where at least 95% substrate cleavage was predicted) matched well for HeLa and HT-29 cells as well as for DLD-1 cells, HCT-116 WT cells, and their XIAP-deficient clones.For MCF-7/C3 cells, we observed deviations (Fig. 3B), suggesting that additional molecular mechanisms may be relevant in this cell line that may have been underestimated by APOPTO-CELL.Nevertheless, the good agreement of apoptosis kinetics in the other cell lines advocates for further employing APOPTO-CELL for a more detailed investigation of molecular mechanisms of apoptosis as performed below.
Apoptosis Kinetics Are Influenced by Different Protein Interactions in HeLa, HCT-116, and DLD-1 Cells-Above, we demonstrated that APOPTO-CELL can distinguish between two classes of responder cells subsequent to MOMP and compared experimentally determined and predicted cleavage of FRET substrate.Encouraged by the agreement between model predictions and experiments, we focused in the remainder of this study on how APOPTO-CELL may be employed to study molecular mechanisms of apoptosis execution and heterogeneity between cells within a population of one cell line and between cell lines.We focused on HeLa, HCT-116, and DLD-1 cells because APOPTO-CELL not only correctly classified them as post-MOMP responder cells but also fairly accurately predicted the experimental kinetics of FRET onset and duration.
We first employed APOPTO-CELL to demonstrate how cell line-specific protein profiles (Fig. 2, A-E) give rise to a differential importance of the molecular interactions in the apoptosis pathway subsequent to MOMP in each cell line.In the model, the strength of each molecular (protein) interaction is defined by kinetic parameters.Changing their values and assessing the influences on apoptosis kinetics may therefore assess how robust apoptosis would be when these interactions are modulated by mutations or pharmacological intervention.
We used individual parameter assessment as described under "Experimental Procedures" to vary all 69 parameters of APOPTO-CELL between 8-fold down-and up-regulation.We then assessed how the time point where 80% cell substrate was FIGURE 3. APOPTO-CELL resembles experimental onset and kinetics of apoptosis execution subsequent to MOMP.A, comparison of predicted and experimental onset of DEVD-FRET substrate cleavage after MOMP for 10 cell lines.Onset in APOPTO-CELL was defined as when 5% of substrate is predicted to be cleaved, whereas experimental onset was determined visually from the traces.Results aligned along the identity line, and correlation was significant (Pearson correlation ϭ 0.860, p ϭ 0.001).B, comparison of predicted and experimental duration of DEVD-FRET substrate cleavage after MOMP aligned for most cell lines with sufficient accuracy.Only cell lines that were predicted to reach 95% were included.Duration in the model was specified as the time between 5 and 95% predicted substrate cleavage.
predicted to be cleaved (t 80 ) changes upon parameter variation.The most influential kinetic parameters were assumed to be those that most altered apoptosis execution kinetics.Results for the above cell lines are shown in Fig. 4, A-C, (using the "average scenario").Parameters were indexed by number, ranked according to the difference of t 80 between maximum and minimum values after down-and up-regulation, and associated with protein interactions according to supplemental Table 1.In all three cell lines, parameter 3, which described caspase-9 activity (p35/p12 fraction cleaving procaspase-3 (33,34)), was most important.In HCT-116 and DLD-1 cells, the speed of XIAP production and the degradation of XIAP⅐caspase-9 complex (parameters 2 and 54 (35,36)) were further important regulators of apoptosis execution kinetics.In turn, the substrate cleavage by active caspase-3 and the speed of apoptosome formation (parameters 68 and 69 (37,38)) were more important in HeLa cells.Together, these systems modeling data suggested that in different cell lines, different protein interactions may be influential for apoptosis execution.
We alternatively employed a method to investigate simultaneous variation of kinetic parameters.Such a method was proposed recently (25).It assumes pairwise variations of model parameters while studying the influence on apoptosis kinetics.Strikingly, results of this method agreed with the individual parameter assessment above in putting emphasis on caspase-9 activity, the speed of XIAP production, and the degradation of XIAP⅐caspase-9 complex with respect to apoptosis regulation.Similar to the individual parameter assessment, it draws attention to the cell-specific importance of different protein interactions (Table 2).
A Workflow to Assess the Robustness of Apoptosis Execution against Changes of Kinetic Parameters-Kinetic parameters depend on cell-specific properties, such as cellular pH or protein motility (39).We therefore first wanted to employ systems modeling to investigate how robust apoptosis execution kinetics were when kinetic parameters changed over a biologically realistic range.Second, we asked whether or not we could find sets of parameter values that improved the overall predictability of APOPTO-CELL for all cell lines.
To address both questions, we constructed a workflow (Fig. 5A).Because it would have required overly high computational costs to assess variations of all 69 kinetic parameters at the same time, we first sought a way to assess what kinetic parameters have the most decisive influence on apoptosis kinetics.We therefore employed both parameter assessment methods described above and assessed the cumulative changes in apoptosis kinetics for all three cell lines under investigation.Using the individual parameter assessment method, we obtained those parameters that exercised the highest influence on apoptosis execution kinetics when they were varied independently over a larger, 8-fold range (set I in Table 2).In contrast, by using the parameter cooperation assessment method, we identified a parameter set to assess the robustness of apoptosis kinetics against variations in those parameters that in combination with each other have the most decisive influence on apoptosis execution kinetics (set C in Table 2).
After identifying two partly overlapping sets of parameters that most decisively influenced apoptosis kinetics (Table 2), we looked for an appropriate range of associated values by scanning reports where those parameters were experimentally determined and which varied over a certain range between publications.The obtained values were sampled in sets (several sets of values belonging to each parameter set I or C; i.e. sets of sets), and each set was included in APOPTO-CELL.The resulting predictions of apoptosis execution kinetics were verified using two filter criteria for each cell.As a first criterion (criterion PX, for protein expression criterion), we required that for XIAPoverexpressing HeLa cells (HeLa XIAP Adv ), the predicted substrate cleavage for cells after 300 min (0.5 M XIAP) was less than 50% (Fig. 5B).This criterion was motivated by previous findings in our laboratory that FRET substrate in HeLa XIAP Adv cells was not entirely cleaved (17).Second, for each of the three cell lines, DLD-1 WT, HCT-116 WT, and HeLa WT, we required that the experimental time point where 80% substrate was cleaved occurred later than the predicted time point for 80% substrate cleavage for the fast scenario and earlier than the predicted time point for 80% substrate cleavage for the slow scenario (criterion S80, for 80% substrate cleavage; Fig. 5C).In conclusion, we devised a workflow that assessed the robustness of model predictions against changes of the most important kinetic parameters under constraints obtained from experimental single cell data.

Prediction of Apoptosis Execution Kinetics Is Robust against Changes in Kinetic Parameters in HeLa and HCT-116 Cells and
More Sensitive in DLD-1 Cells-Having established a workflow, we validated the APOPTO-CELL prediction for several sets of values that belonged to each parameter set and that were chosen from a biologically realistically range screened from the literature (supplemental Table 2).At first, we took sample val-  1. Caspase-9 (p35/p12) activity (cleaving procaspase-3) was the most decisive regulator of apoptosis execution kinetics in all cells (parameter 3).Cleavage of cell substrate by caspase-3 and speed of the apoptosome formation (parameters 68 and 69) were important in HeLa cells (A), whereas XIAP production and the derepression of apoptosome blocking by XIAP (parameters 2 and 54) were important in HCT-116 (B) and DLD-1 cells (C).
ues for each parameter obtained from dividing the entire range of values into uniformly distributed steps (brute force screening).However, assuming 10 parameters and 10 values/parameter led to 10 10 combinations, which would have led to an estimated computation effort of 200,000 days (assuming 1.7 s per calculation) on a single core CPU.We therefore performed brute force screening for only the five most influential parameters of set C (100,000 combinations of parameters, 48-h calculation time).This results in 84% of all sets of parameter values fulfilling criterion PX.We then checked whether other cells or cell combinations also fulfilled (the HCT-116-, DLD-1-, and HeLa-specific) criterion S80 in addition to criterion PX.Interestingly, 81% of all sets passed both criteria for HCT-116 and HeLa cells (Venn diagram in Fig. 6A).In contrast, only 24% of all samples fulfilled both criteria for DLD-1 cells, the same 24% for DLD-1 and HCT-116 cells in combination, and only 22% were found to be valid for HeLa and DLD-1 simultaneously.As a consequence, a low rate of only 22% fulfilled both criteria for all cell lines.
We next evaluated whether brute force screening could be replaced by using screenings of stochastic samples of parameter value sets ("Monte Carlo screening").In particular, we investigated whether or not stochastic sampling allows reducing the sampling size per parameter and therefore allows more parameters to be screened.We first applied Monte Carlo screening to samples of parameter values for the same five parameters of set C as for brute force screening.Using the same numbers of sample sets (100,000), a similar percentage of sets fulfilled both filter criteria for all cells (24%).As before, a high number of sets passed both criteria for HCT-116 and HeLa, whereas only a few sets passed both criteria when DLD-1 cells were taken into account in any combination with the other cell lines (Fig. 6B).Due to the high numbers of screenings and the low numbers of parameters, the good agreement between brute force screening and Monte Carlo screening was not unexpected.However, the agreement between both screening types was still good when we reduced the number of Monte Carlo samples to 1050 (Fig. 6C).This is interesting because only ϳ4 values/parameter range (4 5 ϭ 1024) were tested.This agreement therefore cannot be explained only by stochasticity and suggests that these five parameters (and generally all parameters in this model) do not represent separate degrees of freedom.
The fact that only 4 values/parameter were enough to achieve an acceptable Monte Carlo screening allowed us to extend the screening to 10 parameters, which in brute force screening would have exceeded any reasonable computation time.Applying Monte Carlo screening to the 10 parameters of parameter set C and using 1,050,000 (i.e.ϳ4 values/parameter, 4 10 ϭ 1,048,576), we ascertained that 15% of all parameter value sets fulfilled the filter criteria for all cell lines (Fig. 6D).Again, only comparably few parameter value sets fulfilled them for DLD-1 cells, irrespective of whether DLD-1 cells were tested alone (23%) or in combination with HCT-116 (23%) or HeLa cells (15%).
We then asked whether this dissimilarity of DLD-1 cells compared with HCT-116 and HeLa cells was due to the particular choice of the parameter set C. We therefore repeated Monte Carlo screening using parameter set I with five and 10 parameters (1050 and 1,050,000 screened value sets, respectively).Also here, the filter criteria were only fulfilled for 20% of sets (24%) of parameter values when criterion S80 was applied only to DLD-1 cells.In the screening of five parameters, the number of valid sets was not detectably decreased when this criterion was imposed on all cells (Fig. 6E) but was more pro-FIGURE 5.A workflow to assess the robustness of apoptosis execution and to calibrate systems models.A, a workflow was devised that allowed validation of the robustness of apoptosis execution kinetics against changes in values of kinetic parameters over a biologically justified range (supplemental Table 3).The workflow was confined to study changes in values for the parameter sets I and C according to Table 2. Sets of values that correspond to both parameter sets were obtained using either brute force or Monte Carlo screening.A set of five (or 10) parameter values was considered to explain the experimental apoptosis execution kinetics when two filter criteria were met.These criteria aim to filter those parameter value sets for which the alternative APOPTO-CELL gives reasonable predictions for experimental apoptosis execution kinetics in HCT-116, DLD-1, and HeLa cells and in XIAP-overexpressing HeLa cells under uncertainties in input protein levels.B and C, detailed filter criteria are given.B, a first filter criterion (criterion PX) required that HeLa cells with XIAP less than 0.1 M show full substrate cleavage, whereas HeLa cells with more than 0.5 M show substrate cleavage less than 50%, consistent with our experimental findings for HeLa XIAP Adv cells (17).C, a second filter criterion (criterion S80) required that the time point where 80% of cell substrate was experimentally determined to be cleaved be located between the predicted time points for 80% substrate cleavage of the fast scenario and for 80% substrate cleavage of the slow scenario of apoptosis (see Table 1).Criterion S80 was specifically tested for HeLa, HCT-116, and DLD-1 cells.
nouncedly decreased when 10 parameters were screened (9% of all value sets) (Fig. 6F).Thus, similar to set C, the DLD-1 cell line posed the most severe restriction to fulfill filter criterion S80, whereas HCT-116 and HeLa cells showed more tolerance against parameter variation.These results suggest that the latter two cell lines are more robust in apoptosis execution against variations in enzyme kinetics caused by fluctuations in pH or temperature.
Using different sets of parameter values to parameterize APOPTO-CELL as obtained from the different screening methods, we finally investigated whether any of these alternative parameterizations would better predict experimentally observed apoptosis execution kinetics.We therefore calculated a residual function R 2 and ranked the parameter sets according to how well they describe the experimental FRET traces for the three cell lines under investigation (see "Experimental Proce-  2), choosing kinetic parameters either from a uniform discrete grid (brute force; A) or stochastically drawn from a unique distribution (Monte Carlo; B-F) of the value range of supplemental Table 2. Screening was performed for different parameter sets and different numbers of sets of values.Results of screenings are given by Venn diagrams.A, brute force screening of 100,000 variations of the five most influential kinetic parameters from set C. Of the screened sets, 22% result in acceptable model predictions of apoptosis execution kinetics in all three cell lines.B, Monte Carlo screening of 100,000 variations of the five most influential kinetic parameters from parameter set C leads to a slightly higher amount (24%) of value combinations than brute force screening of the same parameter value range (see A). C, reducing the screening density to only 1050 variations of the five most influential kinetic parameters from parameter set C leads to similar Venn sets as with 100,000 variations.D, Monte Carlo screening of 1,050,000 variations of the 10 most influential kinetic parameters from parameter set C. E, Monte Carlo screening of 1050 variations of the five most influential kinetic parameters from parameter set I. F, Monte Carlo screening of 1,050,000 variations of the 10 most influential kinetic parameters from parameter set I. Raw data and results for the average of all screening can be found in supplemental Table 3. G and H, APOPTO-CELL using alternative parameter values did not improve the agreement between predicted and experimental onset (G) and duration (H) of DEVD-FRET substrate cleavage.A set of values that passed both filter criteria in the screening of Fig. 5 and that was ranked second best with respect to the goodness of fit (R 2 ) in the screening of D was used (see "Experimental Procedures").
dures").The 10 best parameters for each of the six screenings were further analyzed.
However, when the predicted single cell traces that emerged from those 60 sets were visually compared against the experimental FRET traces, no obvious improvement was observed (results not shown).Nevertheless, as exemplified by the value set that was ranked second best for the screening in Fig. 6D, the requirement to better describe the duration of experimental DEVD-FRET substrate cleavage as imposed by filter criterion S80 led to slightly better prediction of duration, which came at the cost of a deviation in the prediction of the onset (Fig. 6, G and H).Interestingly, this improvement in the prediction of the duration was also found for several of those cells that were not included in the filter criteria and which therefore served as an independent test set.

DISCUSSION
Apoptosis impairment is a major factor in cancer cells' resistance to chemotherapy and in clinical relapse (3).The heterogeneity of apoptosis execution between different cancer types and different patients and within tumor tissue presents a major obstacle to the provision of efficient treatments (40).We therefore employed systems modeling to understand the sources of heterogeneity and apoptosis impairment in different cancer cell lines that can arise subsequent to MOMP.
Apoptosis as a Two-step Process Prior and Subsequent to MOMP-Control points of apoptosis lie prior to and subsequent to MOMP (3-5, 26 -29).Prior to MOMP, control of apoptosis is governed by the interplay of antagonizing BCL-2 proteins (41).To this end, Letai and co-workers (4, 5) demonstrated that in situ cancer cells may accumulate tonic pro-apoptotic signals without committing MOMP.Due to having accumulated these stimuli, cancer cells are therefore thought to become more susceptible than non-transformed cells to undergoing MOMP in response to chemotherapeutics, a notion termed mitochondrial priming.It is important to note that within the context of their work, a higher susceptibility manifests in the requirement of a lower chemotherapeutic dose to induce MOMP.In contrast, considering that MOMP is an allor-none-process (42), any prior-to-MOMP dose dependence and any information on mitochondrial priming is unlikely to influence signaling after MOMP.Therefore, in the context of our work, the decision as to whether or not apoptosis is executed subsequent to MOMP is assumed to be dose-independent.This decision is therefore considered to depend solely on the expression levels of the pro-and anti-apoptotic proteins during apoptosis execution during this phase (17,43).Due to these specificities, independent studies of processes prior and subsequent to MOMP can bring more context-dependent and potentially complementary insights into mechanisms of impairment of cancer cell death and its therapeutic implications.
Apoptosis Can Be Impaired Subsequent to MOMP and Can Be Explained by a Systems Analysis Using Absolute Protein Quantifications-APOPTO-CELL assumes MOMP to be present and analyzes the subsequent protein interactions that govern cell death execution.Our study therefore focused on the sources of heterogeneity in apoptosis execution subsequent to MOMP that are introduced in different cancer cell lines and caused by differential protein expression of the proteins that are relevant in this downstream pathway.Our experimental and theoretical results indicate that certain cancer cells can evade apoptosis by tipping the balance toward higher expressions of anti-apoptotic proteins and that this evasion can be understood by a systems analysis.Indeed, loss of APAF-1, overexpression of XIAP, or down-regulation of Smac, which are likely to attenuate the apoptotic pathway subsequent to MOMP, have been linked to carcinogenesis and to mediating chemo-resistance (26 -29).To this end, we recently demonstrated that APOPTO-CELL can predict responsiveness to MOMP-inducing, chemotherapeutic stimuli using patient-individual protein profiles as input (17,43).
Heterogeneity in Apoptosis Signaling Arising from Variations in Protein Expression within One Cell Line-We have studied the effect of protein variations within cell populations and employed APOPTO-CELL to assess the impact of such variations on the kinetics of apoptosis within different cell lines.A representative variation of protein levels in a population was estimated from quantitative Western blotting data.Although this variability was similar in all cell lines, the predicted apoptosis execution kinetics from these data differed.As such, apoptosis in HeLa cells was predicted to be more robust against fluctuations in input proteins.In turn, HCT-116 WT cells showed a higher variability, and LoVo and HCT-116 Smac Ϫ/Ϫ cells exercised a huge variation in apoptosis kinetics, consistent with experimental findings.Because heterogeneity in apoptosis can obstruct the action of chemotherapeutics, APOPTO-CELL may help to assess how entire populations of different cancer cell types respond to such stimuli.
Heterogeneity in Apoptosis Signaling Arising from Cell Linespecific Protein Profiles-We further demonstrated that, due to different protein expression in different cancer cell lines, different protein interactions within the apoptosis pathway may also be important.This fact becomes important when designing and applying chemotherapeutic drugs that aim to enhance or restore apoptosis by targeting a particular protein or a particular protein interaction.As such, Smac mimetics (44) enhance apoptosis by blocking the apoptosis inhibitor XIAP, whereas some compounds directly activate caspases (45)(46)(47), and others target protein degradation (48 -50).Because these drugs target specific protein interactions and because the relative importance of protein interactions may differ between cancer cells, the impact of these drugs on cancer cell apoptosis may vary between different cell lines, different patients, and different cancers.
Heterogeneity and Robustness of Apoptosis Caused by Fluctuations in Enzymatic Activities-Kinetic parameters are factors that relate the speed of a chemical reaction to the concentrations of the reaction partners and therefore are a means to quantify reaction speed of chemical reactions (8).Kinetic parameters are assumed in APOPTO-CELL to describe enzymatic cleavage or enzyme inhibition (17).Although most system analyses assumed these parameters to be constant, enzyme activity is dependent on the cellular environment and influenced by post-translational modifications, temperature, pH, enzyme motility, and cellular viscosity (39,51).Thus, kinetic parameters may either differ between cells (even of the same population) or differ within one cell depending on the physiological context.We therefore provided a workflow to assess the robustness of apoptosis execution kinetics against such changes in kinetic parameters.Using two independent methods, brute force and Monte Carlo screening, to screen the effect of kinetic parameter changes over a range of biologically justified values from literature, we found that HeLa and HCT-116 cells are more robust against parameter variations than DLD-1 cells.This lack of robustness of DLD-1 cells against parameter variations may indicate that certain protein interactions that may have been neglected when constructing APOPTO-CELL may be of particular importance for DLD-1 cells to describe the exact apoptosis execution kinetics.Indeed, during model construction, we disregarded some players and interactions in mitochondrial apoptosis for simplicity.For example, we did not consider the pro-apoptotic protein Omi/HtrA2 (52), which, like Smac, gets released from the mitochondria and binds to and antagonizes XIAP.
Optimization of APOPTO-CELL for Uses as a Clinical Tool to Study Cancer Heterogeneity and Treatment Prognosis-The APOPTO-CELL prediction of the onset of apoptosis in the investigated cell lines correlated well with experimental data, whereas deviations were found when predicting the duration.We therefore wondered whether a better prediction could be achieved by optimizing the set of kinetic parameters used in APOPTO-CELL employing a numerical calibration method that assumed realistic values for these parameters as input.Despite the fact that we investigated simultaneous improvement only in DLD-1, HCT-116, and HeLa cells and did not restrict the search filter too severely, parameter screening and subsequent optimization did not improve the APOPTO-CELL prediction of apoptosis execution kinetics for the studied cell lines.This lack of improvement may be due to the nonexistence of a common set of parameter values for these cells.However, it may also draw attention to necessary refinements of the model, such as the inclusion of further protein-protein interactions.In addition, this study required time-consuming single cell microscopy measurements to detect the rapid apoptosis kinetics at sufficiently high temporal resolution.Therefore, the study focused on a reasonable amount of six parental cell lines and four mutated cell lines to remodel rather than prospectively predict apoptosis sensitivity.However, to better assess the predictive power of APOPTO-CELL, additional, prospective validations to other cell lines would be necessary.Nevertheless, this work along with previous studies from our laboratory suggests that APOPTO-CELL explains cancer cell death or apoptosis resistance subsequent to MOMP in several cancer cells with sufficient accuracy (17,43).In conclusion, our study is the first systems analysis of heterogeneity and impairment of apoptosis subsequent to MOMP, emphasizing the importance of this pathway component in determining the efficacy of apoptosisinducing chemotherapy.

FIGURE 1 .
FIGURE 1. Two classes of response to apoptosis execution subsequent to MOMP.A-C, representative single cell traces of the post-MOMP responder cell line HT-29 (A) and the post-MOMP-impaired responder cell lines HCT-116 Smac Ϫ/Ϫ (B) and LoVo (C) after exposure to 3 M STS.For each cell, time series of DEVD-FRET substrate (CFP/FRET ratio) and mitochondrial polarization measured by TMRM are given.The left panels show time series of single cell fluorescence images, whereas the right panels show quantifications of these values for the representative cells.Onset of loss of TMRM was regarded as the time point of MOMP and indicated by a black arrow.A, HT-29 cells showed rapid DEVD-FRET substrate cleavage after MOMP.B, typical HCT-116 Smac Ϫ/Ϫ cells showed a significant delay.C, a high amount of LoVo cells showed a delay and often no DEVD-FRET substrate cleavage at all.Cells showing no DEVD-FRET cleavage after MOMP were alive at termination of the experiment or underwent necrotic cell death (time point 929).D, from those HT-29, HCT-116 Smac Ϫ/Ϫ , and LoVo cells that underwent MOMP, the fraction of cells that showed onset (left column) and the fraction that showed robust DEVD-FRET substrate cleavage (indicated by a plateau; left column) are depicted.Heat map coloring from dark green to red indicates decrease in the amounts of cells.E and F, time of onset (E) and duration (F) of DEVD-FRET substrate cleavage in HT-29, HCT-116 Smac Ϫ/Ϫ , and LoVo cells indicates a significantly longer duration and a higher cell to cell variability in the latter two cell lines (*, p Ͻ 0.001, analysis of variance followed by pairwise t test for onset; p Ͻ 0.001, Kruskal-Wallis followed by pairwise Mann-Whitney for duration).a.u., arbitrary units.Bars showing meanϩS.E.(E) and median with upper quartile (F).

FIGURE 2 .
FIGURE 2. APOPTO-CELL explains impairment of apoptosis execution subsequent to MOMP.A-E, procaspase-3 (A), APAF-1 (B), XIAP (C), Smac (D), and procaspase-9 (E) concentrations were quantified in HCT-116, DLD-1, LoVo, MCF-7/C3, and HT-29 cells.HeLa cells whose absolute concentrations were quantified previously(17) were used as reference.Error bars, S.E. of three experiments.F-K, APOPTO-CELL predictions for kinetics of substrate cleavage in HT-29-(F), HCT-116-(G), HeLa-(H), DLD-1-(I), DLD-1 XIAP 0/Ϫ -(J), HCT-116 XIAP 0/Ϫ -(K), MCF-7/C3-(L), LoVo-(M), HCT-116 Smac Ϫ/Ϫ -(N), and XIAP-overexpressing (O) HeLa XIAP Adv cells using above protein concentrations as input.Dashed lines indicate the predicted kinetics of substrate cleavage over time after MOMP, calculated from above average protein concentrations.The gray area indicates plausible kinetics for cells of the entire population of the respective cell line (see "Experimental Procedures").F-L, post-MOMP responder cell lines were classified by the model to have rapid kinetics of substrate cleavage in a typical cell of the population (dashed line), a rather synchronized response in cleavage over the population (smaller gray area), and a complete substrate cleavage in almost all cells (Ͼ80% after 300 min).Predictions aligned well with the fast and robust apoptosis execution kinetics for HT-29 cells in Fig.1Aand with previous measurements in the other cells where apoptosis execution was determined to require less than 80 min(17,20,23)).M-O, experimentally classified post-MOMP-impaired responder cell lines were correctly predicted to have slow kinetics of substrate cleavage in an average cell of a population (dashed line) and a low amount of substrate cleavage after 300 min, consistent with the reduced amount of cells showing robust DEVD-FRET substrate cleavage after MOMP.Kinetics of substrate cleavage were predicted to be highly heterogeneous over the population consistent with the findings in Fig.1, E and F, and Ref. 17.

FIGURE 4 .
FIGURE 4. Individual parameter assessment in DLD-1, HCT-116, and HeLa cell lines.A-C, bars show the predicted time from MOMP to 80% caspase-3 substrate cleavage based on average protein concentrations and up to 8-fold down-and up-regulation of individual kinetic parameters.The 10 kinetic parameters that resulted in the strongest change in model predictions between 8-fold down-and up-regulation are depicted.Parameters are associated with protein interactions and indicated by their number according to supplemental Table1.Caspase-9 (p35/p12) activity (cleaving procaspase-3) was the most decisive regulator of apoptosis execution kinetics in all cells (parameter 3).Cleavage of cell substrate by caspase-3 and speed of the apoptosome formation (parameters 68 and 69) were important in HeLa cells (A), whereas XIAP production and the derepression of apoptosome blocking by XIAP (parameters 2 and 54) were important in HCT-116 (B) and DLD-1 cells (C).

FIGURE 6 .
FIGURE 6. Brute force and Monte Carlo screening determined apoptosis execution kinetics in HCT-116 and HeLa cells as robust and in DLD-1 cells as more sensitive to changes of kinetic parameters.A-F, the workflow of Fig.5was applied for parameter sets I and C (Table2), choosing kinetic parameters either from a uniform discrete grid (brute force; A) or stochastically drawn from a unique distribution (Monte Carlo; B-F) of the value range of supplemental Table2.Screening was performed for different parameter sets and different numbers of sets of values.Results of screenings are given by Venn diagrams.A, brute force screening of 100,000 variations of the five most influential kinetic parameters from set C. Of the screened sets, 22% result in acceptable model predictions of apoptosis execution kinetics in all three cell lines.B, Monte Carlo screening of 100,000 variations of the five most influential kinetic parameters from parameter set C leads to a slightly higher amount (24%) of value combinations than brute force screening of the same parameter value range (see A). C, reducing the screening density to only 1050 variations of the five most influential kinetic parameters from parameter set C leads to similar Venn sets as with 100,000 variations.D, Monte Carlo screening of 1,050,000 variations of the 10 most influential kinetic parameters from parameter set C. E, Monte Carlo screening of 1050 variations of the five most influential kinetic parameters from parameter set I. F, Monte Carlo screening of 1,050,000 variations of the 10 most influential kinetic parameters from parameter set I. Raw data and results for the average of all screening can be found in supplemental Table3.G and H, APOPTO-CELL using alternative parameter values did not improve the agreement between predicted and experimental onset (G) and duration (H) of DEVD-FRET substrate cleavage.A set of values that passed both filter criteria in the screening of Fig.5and that was ranked second best with respect to the goodness of fit (R 2 ) in the screening of D was used (see "Experimental Procedures").