A Test of Highly Optimized Tolerance Reveals Fragile Cell-Cycle Mechanisms Are Molecular Targets in Clinical Cancer Trials

Robustness, a long-recognized property of living systems, allows function in the face of uncertainty while fragility, i.e., extreme sensitivity, can potentially lead to catastrophic failure following seemingly innocuous perturbations. Carlson and Doyle hypothesized that highly-evolved networks, e.g., those involved in cell-cycle regulation, can be resistant to some perturbations while highly sensitive to others. The “robust yet fragile” duality of networks has been termed Highly Optimized Tolerance (HOT) and has been the basis of new lines of inquiry in computational and experimental biology. In this study, we tested the working hypothesis that cell-cycle control architectures obey the HOT paradigm. Three cell-cycle models were analyzed using monte-carlo sensitivity analysis. Overall state sensitivity coefficients, which quantify the robustness or fragility of a given mechanism, were calculated using a monte-carlo strategy with three different numerical techniques along with multiple parameter perturbation strategies to control for possible numerical and sampling artifacts. Approximately 65% of the mechanisms in the G1/S restriction point were responsible for 95% of the sensitivity, conversely, the G2-DNA damage checkpoint showed a much stronger dependence on a few mechanisms; ∼32% or 13 of 40 mechanisms accounted for 95% of the sensitivity. Our analysis predicted that CDC25 and cyclin E mechanisms were strongly implicated in G1/S malfunctions, while fragility in the G2/M checkpoint was predicted to be associated with the regulation of the cyclin B-CDK1 complex. Analysis of a third model containing both G1/S and G2/M checkpoint logic, predicted in addition to mechanisms already mentioned, that translation and programmed proteolysis were also key fragile subsystems. Comparison of the predicted fragile mechanisms with literature and current preclinical and clinical trials suggested a strong correlation between efficacy and fragility. Thus, when taken together, these results support the working hypothesis that cell-cycle control architectures are HOT networks and establish the mathematical estimation and subsequent therapeutic exploitation of fragile mechanisms as a novel strategy for anti-cancer lead generation.


Introduction
The capability to gather protein-protein and protein-DNA interaction data, for example using the Yeast Two-Hybrid (Y2H) system [1,2], Fluorescence Resonance Energy Transfer (FRET) techniques [3], quantitative Mass Spectrometry (MS) proteomic or Chromatin Immunoprecipitation (ChIP)-DNA micro-array techniques [4,5] has far outstripped our ability to understand it. Transforming large-scale interaction data into a better understanding of the biomolecular networks underlying disease progression and eventually to new therapies requires integrative tools and strategies. Perhaps one strategy to leverage our knowledge of interaction networks into efficacious therapies would be to identify and exploit weak or fragile mechanisms while avoiding the manipulation of robust network interactions.
Robustness, a long-recognized property of living systems and networks, allows function in the face of uncertainty while fragility, i.e., extreme sensitivity, can potentially lead to catastrophic failure following seemingly innocuous perturbations [6][7][8][9][10]. Different factors can influence why elements of a network are robust or fragile. Venkatasubramanian and co-workers demonstrated that the structure of complex networks can result from a trade-off between efficiency and robustness [11] while You and Yin explored how the environment has shaped the robust properties of bacteriophage T7 [12]. Leibler computationally predicted and later experimentally verified robust features of chemotaxis control networks [13] and Stelling et al., reviewed several examples of robust biological networks [9]. Perhaps no better example of robustness can be found than cell division. The cell-cycle is one of the most fundamental and highly controlled processes in biology. The decision to divide is tightly regulated integrating extracellular signals, such as growth factors and hormones, with intracellular cues that coordinate events leading to division. However, despite extensive control and surveillance subsystems guiding the progression of cells through the division cycle, malfunctions do occur as evidenced by the uncontrolled proliferation underlying many cancers [14]. Thus, while evolutionary pressure may have programmed cells to be robust to shifting nutritional environments or varying growth factor availability, perhaps rare challenges could result in unforeseen consequences. For example, exposure to radiation, exotic chemicals (carcinogens) or even Single Nucleotide Polymorphisms (SNPs) may cause seemingly innocuous changes which manifest themselves in the breakdown of cell-cycle logic. Carlson and Doyle have hypothesized that highly-evolved networks can be resistant to some perturbations while extremely sensitive to others. The ''robust yet fragile'' duality of networks and systems has been termed Highly Optimized Tolerance (HOT) and has been the basis of new lines of inquiry in computational and experimental biology [10].
Sensitivity analysis is an enabling tool for the investigation of robustness and fragility in networks relevant to human health and more generally for model-based knowledge discovery. Cho et al., used sensitivity analysis to study TNF-a-mediated NF-kb signaling where parametric uncertainty was addressed using a monte-carlo parameter sampling protocol; a family of random parameter sets, generated from the best parameter guess, was used to calculate the sensitivity profile in a region of parameter space [15]. Bullinger and coworkers explored the robustness of models of programmed cell death or apoptosis [16] while Stelling et al., computationally identified points of robustness and fragility, using monte-carlo sensitivity analysis and Overall State Sensitivity Coefficients (OSSCs), in models of circadian rhythm [17]. Mahdavi et al., employed sensitivity analysis to better understand stem cell differentiation [18], while Luan et al., used an uncertain mechanistic model of the coagulation cascade in combination with monte-carlo sensitivity analysis, to show that computationally derived sensitive mechanisms were consistent with anticoagulation therapeutic strategies [19]. Sensitivity analysis has also been used to integrate model identification and discrimination with optimal experimental design. Several optimal experimental design and model identification studies are resident in the literature [20][21][22][23][24] along with many techniques to estimate sensitivity coefficients for models composed of ordinary differential equations, differential algebraic and stochastic equations [25][26][27][28].
In this study, we employ mathematical modeling and montecarlo sensitivity analysis to explore the working hypothesis that cell-cycle control architectures are HOT networks. If our working hypothesis is true, then fragile cell-cycle mechanisms (reaction steps) should be overrepresented among experimentally observed malfunctions underlying solid and hematological cancers. Moreover, the manipulation of fragile mechanisms in a therapeutic context, which has been suggested by Kitano [29] to be more likely to elicit an efficacious response from a network or system, should also be prevalent in the treatment literature. We test our working hypothesis by computationally screening three overlapping qualitative models of cell-cycle control architectures; we employ monte-carlo sensitivity analysis and k-means clustering to rank-order mechanisms in cell-cycle and then contrast the predicted fragile and robust mechanisms with literature. If cellcycle control architectures obey the HOT paradigm, then computational identification of fragile mechanisms using proteinprotein or protein-DNA network models could be a novel strategy for anti-cancer lead generation or more broadly as a strategy to identify and exploit weakness in arbitrary networks relevant to human health.

Results
The whole-cycle model of Novak and Tyson (Fig. 1), the G1-S model of Qu et al., (Fig. 2A) and the G2/M-DNA damage model of Aguda (Fig. 2B) were implemented from literature and screened for fragile mechanisms using monte-carlo sensitivity analysis [30][31][32]. The Novak and Tyson model, which employed a complex description of the G1/S and G2/M checkpoints, programmed protein expression and degradation, was composed of 18 dynamic species, 4 species constraints and 74 parameters. The mass-action G1/S and G2/M-DNA damage models described only the molecular logic in their respective checkpoints; the G1/S model was composed of 16 dynamic protein balances, 2 species constraints and 44 parameters while the G2/M-DNA damage model consisted of 15 dynamic protein balances,1 constraint and 40 parameters. Parameter values for each model were taken from literature. Unreported initial conditions were adjusted so that simulated model trajectories were qualitatively consistent with published values (Supplementary Material Figure S1). The published parameter sets, with fixed initial conditions, were used to generate random parameter sets (N = 500, unless otherwise noted) where each nominal parameter was perturbed by up to 650%, 61-order, or 62-orders of magnitude. Overall State Sensitivity Coefficients (OSSCs) were calculated over the random parameter families for each cell-cycle model using three different numerical algorithms. For each model, the mean OSSC values were ranked-ordered and plotted. The Area Under the Curve (AUC) was used to measure the cumulative sensitivity contribution of each parameter. A cumulative cutoff of 95% of the overall sensitivity was used to establish the list of mechanisms (Supplementary Material Figure S2) which were clustered into three groups (high, medium and low sensitivity) using a k-means algorithm.
Approximately 65% of the G1/S mechanisms (reaction steps) were responsible for 95% of the sensitivity, conversely, the G2-DNA damage network showed a stronger dependence on a few interactions. Of the 44 G1/S reactions steps, 29 were responsible for 95% of the sensitivity (Supplementary Material Figure S2). The distribution of fragility was not specific to any single class of interaction ( Table 1). The dephosphorylation of CDC25, the expression of cyclin E, the degradation of the cyclin E-CDK2 complex, and the concentration of the transcription factor E2F were classified as the most fragile reaction steps in the G1/S checkpoint (Table 1, cluster I). A previous model of G1/S by Aguda et al., [33] found that although pRB and cyclin E-CDK2 formed a positive feedback loop, they did not form a sharp robust switch at the restriction point, i.e., the increase in active cyclin E-CDK2 concentration was gradual and sensitive to model parameters. However, addition of CDC25 to the positive cyclin E-CDK2-pRB feedback loop, made the restriction point robust to model parameter variation, thus supporting our findings of the importance of CDC25 interactions. The synthesis, activation and degradation of CKIs, the expression and degradation of CDC25, pRB concentration, the expression of cyclin D and cyclin E-CDK2 mechanisms dominated the second-tier of G1/S fragility (Table 1, cluster II). Tier-three of G1/S fragility involved several cyclin D mechanisms, cyclin E-CDK2 activity and E2F mediated cyclin E expression ( Table 1, cluster III). When taken together, the most heavily implicated G1/S protein was cyclin E, with 11 of 29 mechanisms, followed by CKIs with six, CDC25 and cyclin D were each involved in five fragile mechanisms and E2F and pRB were each listed twice. Moreover, 16 of the 29 fragile parameters were functionally associated with cyclin E and cyclin E-CDK2 activity. As expected, the expression and degradation of the G1/Sphase cyclins and their associated CKIs were predicted to be important. However, the expression and degradation of cyclin E and other it's interactions were ranked higher than the corresponding cyclin D mechanisms with the exception of the dissociation of the cyclin E-CDK2-CKI complex. The G2-DNA damage network showed a stronger dependence on a few mechanisms when compared with G1/S; ,32% or 13 of 40 mechanisms accounted for 95% of the sensitivity (Supplementary Material Figure S2). Consistent with G1/S, no single class of mechanism dominated the fragility list. The most sensitive mechanisms were related to the generation and degradation of the cyclin B-CDK1 complex otherwise known as the Maturation Promoting Factor (MPF) ( Table 2). The top five mechanisms were either directly or closely associated with the formation and activity of MPF while mechanisms leading the deactivation of MPF, e.g., the expression, degradation and activity of p21, 14-3-3s and Wee1 phosphorylation dominated the remaining eight mechanisms (Table 2, cluster III). Activation of inactive MPF complex, whose expression is negatively regulated by p53, was the most sensitive G2 mechanism (Table 2, cluster I), followed by preMPF generation, activation and transport of CDC25 into the nucleus ( Table 2, cluster II). The finding that all CDC25 related mechanisms were more fragile than Wee1, is consistent with earlier work by Aguda [34] which showed that even though both Wee1 and CDC25 form a phosphorylation-dephosphorylation (PD) loop with MPF, only CDC25 coupling gave rise to qualitatively different behavior. Interestingly, while the generation of p53 itself was not predicted to be sensitive, interactions involving p53 were prevalent, e.g., the expression of inactive MPF and p21, both of which are regulated by p53, were predicted to be sensitive. Approximately 77% of the Novak and Tyson parameters (57 of 74) were responsible for 95% of the sensitivity (Supplementary Material Figure S2). Both global and local components of the model were predicted to be fragile. The most sensitive global mechanism was the translational efficiency while local mechanisms such as activation of IE (hypothetical protein which activates the E3-ligase CDC20), expression of cyclin B and CDH1 degradation were also predicted to be fragile (Table 3, cluster I). The secondtier mechanisms were associated with deregulation of programmed proteolysis (Table 3, cluster II). Interestingly, while the percentage of mechanisms responsible for 95% of the sensitivity of the Novak and Tyson model was the largest of the three models, several mechanisms in cluster III had small OSSC values, including most of the G1/S checkpoint logic. Thus, sampling the complex Novak and Tyson model produced less information than the mechanistic mass-action based G1/S and G2-DNA damage models.
The qualitative conclusions drawn from sampling the cell-cycle models were robust to the choice of solution method and the size of the parameter perturbation but sensitive to the number of parameter sets sampled. Three different numerical techniques were used to solve the sensitivity equations to control for possible numerical artifacts. The ODE15s routine of Matlab (The Mathworks, Natick MA), a third-order backward-difference implicit method (BDF3; see Supplementary Material S1) and forward finite difference (FD), generated qualitatively similar sensitivity results (Fig. 3). The lowest Spearman rank between any two methods (ODE15s versus FD for the G1/S model) was 0.91 indicating a worse case correlation of approximately 91%. Interestingly, while the Spearman rank indicated good agreement between the solution methods, there were statistically significant shifts in OSSC values indicating the solution methods systematically shifted mechanisms, i.e., different OSSC values were calculated but the order or ranking of mechanisms was maintained (see Supplemental Material Table S1). Two additional sampling controls were conducted to verify the robustness of the qualitative conclusions drawn from our analysis. First, the perturbation size used to generate the random parameter families was varied to test if different conclusions would have been drawn with different perturbation sizes; OSSC values computed over random parameter families generated using 650%, 61-order and 62-orders of magnitude showed no qualitative difference as quantified by the Spearman rank correlation for the G1/S model (Fig. 4). The worst case correlation of 0.90 was observed between the 650% and 62orders of magnitude cases indicating on average 90% of the conclusions drawn between the two cases were consistent (Fig. 4C). Such a strong correlation in Spearman ranks across 2-orders of magnitude in the parameter values might suggest that network structure (connectivity) is more important than parameter values. Comparison of exactly similar mechanisms across the three models supported the hypothesis of connectivity dominance where mechanisms classified as either fragile or robust in the G1/S and G2-DNA damage models were also predicted to be important in the Novak and Tyson model, albeit with different ranks (Table 4). There were 11 mechanisms which appeared exactly in each model, 10 mechanisms were classified similarly while one was ranked inconsistently. Second, the cumulative Spearman rank correlation between sensitivity results generated using the ODE15s, BDF3 and FD methods for each model was calculated as a function of the number of parameter sets sampled. While the cumulative Spearman rank converged to the population mean as the number of parameter sets increased, a population size dependence was observed (Fig. 5). For each model, the results reported were obtained in the region of convergence; hence, no new information would have been gained if additional random parameter sets were sampled.

Discussion
Literature evidence supports the hypothesis that computationally identified fragile cell-cycle interactions represent efficacious targets. Consider the fragility of CDC25 mechanisms. Boutros et al., recently reviewed the role of CDC25 phosphatases and CDC25 inhibitors in human cancer progression and treatment [35]. While the inhibition of CDC25 as a cancer treatment strategy is still in the laboratory stage, several CDC25 inhibitors in development have shown promising results. The CDC25 inhibitor PM20 inhibited growth in human hepatoma-derived Hep3B celllines at a inhibitory concentration (IC) .700 nM, PM-20 also inhibited the growth of several other cell-lines, albeit at higher ICs [36]. BN82685, which inhibited CDC 25A, B and C in-vitro and in- vivo and repressed the growth of HeLa and human pancreatic tumor Mia PaCa-2 xenografts in athymic nude mice, also inhibited the growth of human cell lines resistant to cytotoxic drugs e.g., the human myeloblastic leukemia cell-line HL-60 [37]. The CDC25 antagonist, CPD-5, inhibited the growth of the rat hepatoma cell-line JM-1 in-vitro and the mouse cancer cell-line tsFT210 through selective inhibition of CDC25 [38]. Thus, inhibition of CDC25 represents a viable treatment option which could be pursued further in the clinic. Inhibition and degradation of the active cyclin E-CDK2 complex, the second ranked mechanism in the G1/S network, has also been exploited as a treatment strategy. Bristol-Myers Squibb (BMS) developed BMS-387032, a cyclin E-CDK2 inhibitor, with an IC50 of 95 nM [39]. Preclinical and phase I ovarian cancer studies demonstrated that BMS-387032 possessed better efficacy than Flavopiridol, a promiscuous CDK inhibitor [40]. Flavopiridol, the first cyclin dependent kinase inhibitor in clinical trials, alone or in combination with other drugs is currently being investigated in 52 active phase I or II trials [41]. Flavopiridol has been proposed for the treatment of recurrent, locally advanced, or metastatic soft tissue sarcoma [42], lymphoma and multiple myeloma [43], metastatic breast cancer (with Trastumuzumab) [44] or in combination with other drugs (Cisplatin and Carboplatin) for the treatment of advanced solid tumors [45]. Cyclin E expression, the fourth ranked mechanism in the G1/S model, has also been explored therapeutically for the treatment of pancreatic and lung cancers [46,47]. The correlation between fragility and treatment strategy was also found to hold for the G2/M-DNA damage network. The activation of preMPF (cyclin B-CDK1 complex), catalyzed by CDC25, was predicted to be the most sensitive mechanism in the G2/M-DNA damage model while three of the four tier-two G2/M-DNA mechanisms were associated with CDC25 activity. Bryostatin-1, a protein kinase C (PKC) inhibitor and antagonist of the cyclin B-CDK1 complex, has been explored in the clinic for the treatment of multiple myeloma [48], relapsed non-Hodgkin's lymphoma and chronic lymphocytic leukemia [49]. In preclinical models, Bryostatin-1 has demonstrated singleagent activity against B16 melanoma, M5076 reticulum sarcoma and L10A B-cell lymphoma [50] and has been shown to disrupt cyclin B-CDK1 complex formation and activity by several different mechanisms [51,52]. When taken together, the top fragile mechanisms for both the G1/S and G2/M phases of the cell-cycle, estimated by monte-carlo sensitivity analysis, were found to be consistent with on-going preclinical and clinical trials for the treatment of a broad spectrum of human cancers. Modulation of translational efficiency and the manipulation of programmed proteolysis, prominently featured among the group of fragile mechanisms across all the models, are also active areas of therapeutic development. Initiation of translation in eukaryotes is thought to be rate limiting [53] and overexpression of initiation components, for example the initiation factor elF4E, occurs frequently in human cancers [54]. Arnqvist and coworkers explored translation inhibition in MCF-7 breast cancer cells following cycloheximide, puromycin or emetine exposure in the presence and absence of Insulin-like Growth Factor1 (IGF-1) [55]. Addition of puromycin, cycloheximide and emetine in the absence of IGF-1 resulted in increased apoptosis at 48 hr relative to the control, however, when IGF-1 was present, a concentration dependent reduction in apoptosis was observed. Bjornsti and Houghton recently reviewed another small molecule translation inhibitor, Ramapycin [56], which inhibits the Target of Ramapycin (TOR) protein, a serine/threonine kinase involved in translation and other functions. While Ramapycin has FDA approval as an immunosuppressant, development of anticancer therapies has been slow despite anti-tumor activity against established solid-tumor models [57,58]. Ramapycin analogs have been evaluated in clinical trials for the treatment of different indications including pediatric patients with relapsed or refractory acute leukemia and renal-cell carcinoma [56,59]. Peptide inhibitors have also been used to downregulate translation e.g., BL22, an immunotoxin developed for the treatment of Chronic Lymphocytic Leukemia (CLL) [60], consists of the variable FV   [61], has been the target of several different therapeutic developments [62]. The ubiquination of target proteins involves the coordinated activity of the ubiquitin activating enzyme family (E1), the ubiquitin-conjugating enzyme family (E2) and the ubiquitin ligase family (E3) [63]. While E1 malfunctions have not been observed in cancer, deregulation of E3 and to a lesser extent E2 activity has been directly linked to cancer progression [63]. The Novak and Tyson model has only a skeleton representation of UPS, however, it does explicitly represent Cell Division Cycle protein 20 (CDC20), CDH1 and Anaphase Promoting Complex/Cyclosome (APC/C), all of which are E3 components. APC/C is the core subunit to which the adapter proteins CDC20 and CDH1 bind [64][65][66]. Inhibition of specific E3 ligases remains a technical challenge [67], however, cisimidazoline analogs called Nutlins have been developed which inhibit MDM2, an E3-ligase responsible for the recognition of p53. Activity of Nutlins-3 against a human osteosarcoma xenograft model in nude mice showed 90% inhibition of tumor growth relative to control [68]. While multiple lines of experimental evidence support the assertion that malfunctions in fragile mechanisms are implicated in solid and hematological cancers, some evidence is contradictory. CDC25 activity, cyclin E expression and activity of cyclin E-CDK2 were the largest group of fragile G1/S mechanisms. Traditionally, cyclin E expression and cyclin E-CDK2 activity were thought to be critical for cell-cycle progression [69]. Ohtsubo et al., have shown that cyclin E-CDK2 activity was maximum during the G1/S phase and overexpression of cyclin E accelerated cell-cycle progression [70]. Lucas et al., showed that abnormal cyclin E but not Cyclin D1 expression was able to override G1 arrest by the INK4a family of CKIs [71]. Keyomarsi et al., found that cyclin E expression plays a strong role in human breast cancer tumors and the cyclin E-CDK2 complex remains active throughout the cell-cycle suggesting the now established hypothesis that truncated (deregulated) cyclin E variants were responsible for the constitutive function of cyclin E-CDK2 in breast cancer tumors [72,73]. Recent studies, however, have challenged the traditional role of cyclin E. Deletion of both cyclin E genes was lethal in-utero but deletion of cyclin E1 or cyclin E2 was tolerated with no obvious abnormalities [74]. Interestingly, double cyclin E knockout mice were born alive if cyclin E was restored in the embryonic component of the placenta [74] and CDK2 null mice were born viable and healthy [75]. Thus, while the cyclin E and CDK2 knockout studies seem to contradict the essential role of cyclin E, clinical evidence suggests further studies are required. Evidence supporting the involvement of other fragile components, such as the concentration of E2F and pRB (constraints in the G1/ S and Novak and Tyson models), is also prevalent in the literature [76,77]. However, contradictory evidence suggests that the role of cyclin D mechanisms maybe complex. Sensitivity analysis suggested that cyclin D-CDK4/6 and cyclin D-CDK4/6-CKIs trimer mechanisms were robust or only moderately sensitive while cyclin D expression was fragile in the G1/S checkpoint. While Keenan et al., demonstrated in IIC9 Chinese hamster embryonic fibroblasts that cyclin E expression renders cyclin D-CDK4 dispensable [78], overexpression of cyclin D variants, particularly cyclin D1, has been observed in several human cancers [79,80]. Moreover, cyclin D1, D2 or D3 2/2 mice displayed tissue specific phenotypes including defective proliferation [81][82][83]. However, while mice lacking all the cyclin D genes died by day E17.5 of gestation, most tissue and organs were formed by day E13.5 indicating that cyclin D was not required for embryogenies [84]. When taken together, the retrospective cyclin E studies in breast cancer patients and the CDC25 studies support the hypothesis that malfunctions in fragile mechanisms are strongly implicated in cancer progression. However, the cyclin E and CDK2 knockout studies as well the confusing role of cyclin D suggests a more nuanced perspective in which redundant proteins or subsystems might be able to compensate for malfunctions in fragile mechanisms.
Consistent with the conjecture of Kitano, the anecdotal comparison between predicted fragile mechanisms and literature suggested that cell-cycle control architectures are HOT networks [29]. However, while different controls were conducted to ensure the fidelity of the monte-carlo sampling protocol, the mathematical models being explored were coarse-grained and not structurally complete. While quantifying the impact of structural uncertainty remains a critical challenge, the general correlation  between efficacy and fragility appears to be model independent as other studies have yielded similar results [19]. Moreover, initial results presented here suggest that while the quantitative values of sensitivity coefficients calculated using different models with overlapping biology will change between models, the qualitative conclusions drawn may be invariant. However, this conclusion is likely true only for a subset of mechanisms. One possible strategy to explore structural uncertainty would be to construct detailed subsystem models of the coarse-grained components which were determined by our analysis to be fragile, e.g., translation or UPS. While this top-down strategy does not specifically address structural uncertainty, it does allow us to determine the molecular interactions which are perhaps mediating fragility in the coarsegrained model. A second critical issue is the choice of sensitivity metric. OSSCs quantify the overall impact that a parameter has; however, other measures of sensitivity might be better suited for analysis of the cell-cycle. Doyle and colleagues have established tools for the analysis of mammalian circadian rhythm that could prove useful in understanding how fragility influences phenotypic properties such as division frequency [17,85,86]. A third critical issue not addressed in this study was safety. Highly efficacious strategies have resulted in unwanted and possible harmful side effects, e.g., the association of rofecoxib with adverse cardiovascular events [87]. While there may not be an obvious linkage between fragility and safety for single agents, initial retrospective studies by Luan et al., using combinations of coagulation inhibitors, have suggested that shifts in mechanism rank could be used to understand molecular drug-drug synergies [19]. were calculated for each cell-cycle model over a family of random parameters sets (N = 500 unless otherwise noted) generated by randomly perturbing the published set by 61-order of magnitude. Three different numerical methods were used to solve the sensitivity equations to control for numerical artifacts. A-C: Sensitivity results for the Novak and Tyson model [32]. D-F: Sensitivity results for the G1/S checkpoint model of Qu et al., [31]. G-I: Sensitivity results for the G2/M-DNA damage model of Aguda [30]. The different numerical techniques used to solve the sensitivity equations yield qualitatively similar results as quantified by the Spearman rank correlation between any two methods (lower right-hand corner of each plot). doi:10.1371/journal.pone.0002016.g003

Model formulation and sensitivity analysis
The cell-cycle models used in this study [30][31][32] were represented as systems of Differential Algebraic Equations (DAEs): where xMR m denotes the concentration vector, f(x, p)MR m denotes the mass balance equation vector describing the kinetics and connectivity of the cell cycle network and pMR p denotes the parameter vector. The diagonal m6m matrix H contains 1's for dynamic elements of the concentration vector, 0 otherwise (constraints).
The fragile elements of the cell-cycle networks were determined by computing Overall State Sensitivity Coefficients (OSSC) [17]. OSSC values were calculated by first calculating the first-order sensitivity coefficients (at time t k ): which are solutions of the equation: subject to the initial condition s j (t 0 ) = 0. The quantity j denotes the parameter index, P denotes the number of parameters and s j denotes the m61 vector of first-order sensitivity coefficients with respect to parameter j. The Jacobian matrix (A) and the matrix of first derivatives of the mass balances w.r.t the parameter values (B) (whose columns are denoted by b j ) are given by: where x denotes a point along the nominal or unperturbed system solution. We solved the sensitivity equations for each parameter using three different numerical methods to control for possible artifacts; a 3-order Backward Difference (BDF3) method was compared with forward Finite Difference (FD), and the fifth-order variable step-size ODE15s routine of Matlab (The Mathworks, Natick MA). The matrices A and B were estimated numerically at each time step using a generalized gradient algorithm [88]. Overall State Sensitivity Coefficients (OSSC), first used by Stelling et al., to characterize mechanisms in circadian rhythm as fragile or robust [18], were calculated for each parameter j: The quantity N T denotes the number of time points used in the simulation while N s denotes the number of proteins/protein complexes in the model. To account for parametric uncertainty, the OSSC values (S oj ) were calculated over a family of random   parameter sets; we randomly perturbed each nominal parameter by up to 61-order of magnitude then solved the sensitivity balances for each family member. To control for perturbation effects, two other random parameter families were also tested (650% and 62-orders of magnitude, N = 500).

Statistical and clustering analysis of OSSC values
Three different tests were performed to identify large statistically significant shifts in the OSSC values. The OSSC values calculated over the family of parameter sets were assumed to be normally distributed. The statistical significance of shifts in OSSC values for each algorithm relative to ODE15s (control) were determined by performing a Welch t-test with the null hypothesis that the means of the OSSC values were equal at a 1% significance level [89]. The list of significant OSSC values was further restricted to only those shifts with a magnitude larger than a specified z-score (1.0) away from the squared mean displacement over the significant OSSC values. We defined the displacement of an OSSC value relative to the control as: where S S C oj denotes the mean OSSC value over the family of parameter sets for parameter j in the control while S S q oj denotes the same quantity for algorithm q. A significant shift in OSSC value was accepted if: where z denotes a desired z-score, s dq denotes the standard deviation of the total displacement over all significant OSSC values for the q th numerical algorithm and m dq denotes the mean of the significant displacements for algorithm q. Large statistically significant shifts in OSSC values, while perhaps indicative of the shifting importance of mechanisms, do not guarantee that mechanisms are qualitatively different between the algorithms considered (see Supplementary Material Table S1). The Spearman rank correlation denoted by r and defined as: was used to measure the difference in qualitative ranking of mechanisms between algorithms considered. The quantity d i denotes the difference in the ordinal rank of mechanisms between algorithms or perturbation size, N denotes the number of pairs of values and P denotes the number of parameters considered. The Spearman rank is bounded by 21$r$1; a Spearman rank of one indicates that two ranked lists are identical, a Spearman rank of negative one indicates a perfect negative correlation, while a Spearman rank of zero indicates that two ranked lists are uncorrelated.
The distributions of OSSC values obtained from monte-carlo sampling were clustered using a k-means algorithm [90]. The mean and standard deviation obtained from the monte-carlo sensitivity analysis was used to estimate the underlying OSSC distribution (N = 500 points) where the OSSC values were assumed to be normally distributed. One-hundred different clustering attempts were run for each model to control for clustering artifacts. The most probable configuration was reported.

Supporting Information
Material S1  [30]. E-F: Concentration profiles for the Cdh1 protein and the Cdk1:CycB complex versus time for the reimplementation (E) and published whole-cycle model of Novak and Tyson [32]. In all cases the reimplemented models were qualitatively consistent with published results. Table S1 Statistically significant shifts of Overall State Sensitivity Coefficients (OSSCs) between solution methods computed using the Welch t-test. The mean and one standard deviation of the OSSC score computed over the family of random parameter sets is reported. Only shifts recorded with a p-value of 0.01 and zscore of 1 are shown. Found at: doi:10.1371/journal.pone.0002016.s004 (0.04 MB DOC)