Confidence from uncertainty - A multi-target drug screening method from robust control theory
© Luni et al; licensee BioMed Central Ltd. 2010
Received: 26 May 2010
Accepted: 24 November 2010
Published: 24 November 2010
Robustness is a recognized feature of biological systems that evolved as a defence to environmental variability. Complex diseases such as diabetes, cancer, bacterial and viral infections, exploit the same mechanisms that allow for robust behaviour in healthy conditions to ensure their own continuance. Single drug therapies, while generally potent regulators of their specific protein/gene targets, often fail to counter the robustness of the disease in question. Multi-drug therapies offer a powerful means to restore disrupted biological networks, by targeting the subsystem of interest while preventing the diseased network from reconciling through available, redundant mechanisms. Modelling techniques are needed to manage the high number of combinatorial possibilities arising in multi-drug therapeutic design, and identify synergistic targets that are robust to system uncertainty.
We present the application of a method from robust control theory, Structured Singular Value or μ- analysis, to identify highly effective multi-drug therapies by using robustness in the face of uncertainty as a new means of target discrimination. We illustrate the method by means of a case study of a negative feedback network motif subject to parametric uncertainty.
The paper contributes to the development of effective methods for drug screening in the context of network modelling affected by parametric uncertainty. The results have wide applicability for the analysis of different sources of uncertainty like noise experienced in the data, neglected dynamics, or intrinsic biological variability.
Biological systems are hierarchically organized, from genes to proteins up to the organism level. At the cellular level, complex interconnected networks include metabolic signalling, signal transduction, and transcriptional regulatory networks . Some general features of biological networks have been explored computationally, such as robustness , modularity , control coefficients , and connectivity properties . Robustness is defined as the ability to maintain functional performance in the presence of uncertainty [2, 6], and it is particularly relevant in therapy design as drug effectiveness should be independent from predictable sources of variability.
Complex diseases often exploit the same strategies present in healthy networks to gain a robust status . Diseases such as diabetes, cancer, bacterial and viral infections, represent multiple disruptions within the host network structure rather than single events, such as a DNA point mutation . Signalling redundancy, feedback, and other network strategies adopted by the disease, ensure that it will be robust to disturbances within its architecture. Hence, single-target therapies fail in many cases because network characteristics are not accounted for during target identification [8, 9]. On the other hand, multi-drug therapies (MDT) have been proven to be effective for many complex diseases [10, 11]. Network-based design of MDTs can improve current drug regimes [11–14] by identifying targets that both moderate the immediate characteristics of the disease while disarming its robustness strategies . Furthermore, synergy within MDTs may reduce the required drug load, hopefully minimizing side effects [15, 16].
Some MDTs are currently used to treat chronic diseases and to boost antibiotic potency. AIDS infections routinely require a drug regimen of reverse-transcriptase inhibitors and protease inhibitors . Oncological chemotherapeutic regimens often involve the combination of cyclophosphamide, hydroxydaunorubicin, oncovin, and prednisone, abbreviated as CHOP . Augmentin, an amoxicillin-based antibiotic, contains clavulanic acid to inhibit a known mechanism of amoxicillin degradation . In comparison to their single-perturbation counterparts, these MDTs often show an order of magnitude greater efficacy . Most MDTs to date have been identified in an ad hoc fashion, relying on observational studies of previously available drug lines. Many pharmaceutical companies are now embracing the idea of a priori design of MDTs using in silico modelling and analysis to rapidly identify candidate targets .
Optimizing drug combinations and concentrations produces an unmanageable number of possible therapies to explore, demanding efficient computational methods of screening . Furthermore, it is unreliable to extrapolate the therapeutic efficacy from the necessarily few conditions tested. For example, a potential concentration-dependent synergistic behaviour may occur at intermediate concentrations not considered during experimentation. This situation is not unlikely considering that strongly nonlinear behaviours have been recognized in biological systems, such as switching or bistabilities .
For drug screening to succeed, additional insight into the biological mechanisms of drug action at the cellular level is needed to increase the predictability of the therapy. High-throughput experimental techniques are providing the data required to understand the connections between the biochemical nodes in the cellular sub-networks underlying specific functions. The causal relationships between these components are being explored by dynamic modelling through a continuous process of expansion and refinement. The most widely-used representation of the biochemical reaction network is a dynamic and continuous description, based on a system of ordinary differential equations (ODEs) , where the variables represent the concentrations of the components, and their change over time is simulated. Many ODE models are currently under development to gain insight into complex diseases, such as diabetes [23–25], and will be invaluable for future drug discovery, as reviewed elsewhere [26–28]. More than 200 network models from the literature have been curated and included in publicly accessible databases, such as Biomodels, BioPax and the CellML Model Repository. Systems Biology Markup Language (SBML) was created to standardize the description of biochemical networks, enabling communication between people and software , and paving the way for a biochemically detailed artificial organism reconstruction . ODE models can be interrogated to test hypotheses of cellular response to drug combinations, considering whole sets dosage permutations and used to discover optimal points of manipulation within the network [13, 16]. These models have the potential for in silico testing MDTs at reduced cost and time . Despite improvements in the accuracy of biological models, their reliability is often limited by parameter uncertainty. Even at best, parameter values can be inferred by experimental data as a range of values, rather than a fixed one. While increasingly precise experimental measurement methods are being developed, cell-cell heterogeneity in tissues and stochastic noise, the consequence of the small copy number of some intracellular components, are intrinsic sources of uncertainty and require ad hoc methods of analysis.
Case study description
Healthy performance and potential therapies
In practice, the performance bounds are derived by the standard deviation of the experimental data. In this work, the system "noise" is artificially generated simulating the system with the Stochastic Simulation Algorithm, SSA . A smooth performance envelope is then defined to approximately contain the concentration profiles resulting from these simulations, as shown in Figure 3A and explained in the Methods section.
Multiple therapeutic approaches can be investigated that aim at restoring the normal output concentration in the presence of a diseased input condition. A drug effect on the system can be modelled as a parameter perturbation, i.e., modifying a component's rate of synthesis, degradation, or interaction with other elements in the network. We first inferred the set of potential therapies by fitting the healthy output curve in the presence of a diseased input, u tot,d , targeting up to 4 parameters at a time. Thus, each therapy model is in the form of the equation system presented in Figure 2B, with a diseased input, u tot,d , and a different parameter set obtained solving a least-square optimization problem that minimizes the deviation of its output from the healthy one. A total of 56 possible therapies, i.e. combinations of the n = 6 parameters, were obtained. A comparison between the outcome, y dt , from each therapy and the performance envelopes is shown in Figure 3B, where the simulations were performed starting from the diseased steady-state in absence of any source of uncertainty.
Selection of therapies for nominal performance
Uncertainty description and robust performance
where k ∈ [kmin, kmax] is a generic parameter of the model, k mean = (kmin + kmax)/2, r k = (kmax - kmin)/(kmax + kmin), and δ k ∈ ℝ and |δ k | ≤ 1. In this case study, we assume that all parameters have a relative fluctuation of 45% about their mean value (i.e., r k = 45%).
The conditions for nominal performance (without uncertainty) can be extended to the case of an uncertain model. Specifically, a therapy meets the criterion for robust performance if, for any set of parameters within the defined uncertainty range, no output trajectory crosses the healthy performance envelope boundaries, when using the healthy steady-state as initial condition. A direct comparison between each therapy's output trajectories and the healthy performance envelope, as in the previous section, is not feasible for all the values of the uncertain parameters. The advantage of employing SSV analysis becomes apparent in this situation.
Rearrangement of the model in M-Δ form
SSV analysis is a tool developed in control theory to study the performance of systems affected by uncertainty . We provide here an intuitive understanding of how it works, and refer to textbooks in the field for a more technical explanation . Before SSV application, some preliminary steps are needed to recast the model in a suitable form, including model Jacobian linearization, and Laplace transforms. They are well-known techniques in control theory and numerical algorithms to perform them are readily available in technical software such as Matlab.
General aspects on SSV
where I is the identity matrix, and k m a scalar factor. A result well-known in control theory, simplistically stated here, is that, when det(I - MΔ) = 0, then the loop in Figure 5 becomes unstable, i.e., in our case, the performance is not fulfilled. As Δ is a matrix whose elements are uncertain, the above minimization problem is solved over all possible Δ's that are normalized and have the structure that we mentioned in the previous section.
The value min(k m ) represents the smallest perturbation that destabilizes the system, and μ RP is its reciprocal. Thus, μ RP = 1 means that there exists a perturbation, within the uncertainty description, that is large enough to pull the output exactly to the limit of the performance envelope. The model meets the conditions for robust performance if and only if μ RP < 1. Details on how μ RP is computed are available in the literature , and algorithms are also included in technical software, like Matlab.
Selection of therapies for robust performance by SSV
Increase* or decrease# of parameter values in the robust therapies respect to their nominal ones
Figure 6B shows the output dynamics of therapy 30. First, the nominal values of the parameters were randomly changed by up to ± 45%. 100 parameter sets were generated, and used in the linearized model to simulate the diseased treated system, with diseased steady-state as initial condition. In all cases the output steady-state value falls within the healthy performance envelope, as a confirmation of SSV analysis results. Then, the nominal nonlinear model of therapy 30 was simulated by SSA. SSV analysis does not guarantee performance in this case, as noise in component concentrations was not included in the uncertainty description, which was only parametric. Nonetheless, the stochastic envelope generated falls into the performance envelope even in this case after a transient (Figure 6B), increasing our confidence on the efficacy of this therapy in the presence of unexpected uncertainty.
In this paper we have proposed a new method for MDT selection, taking advantage of SSV analysis, a tool already successfully applied in other fields such as aeronautics . We have evaluated the feasibility of using this tool for drug screening by a simple case study, essentially a network given by an enzymatic reaction negatively regulated by its own product. While therapies can be easily selected based on a criterion of nominal performance, the importance of SSV application is apparent in presence of parametric uncertainty, when, to the best of our knowledge, alternative methods are not available.
Through the case study, we demonstrated the relevance of considering the effect of structured uncertainty, i.e. parametric noise, as only 5 therapies, out of 41 showing nominal performance, were robust. From a network perspective, the results emphasize how MDTs offer greater potency in regulating specific targets. In fact, all the 5 therapies passing the screening involved multiple perturbations. Furthermore, these resulted in therapies that are also less susceptible to internal biological fluctuations, as demonstrated by the SSA simulation of therapy 30, whose results are shown in Figure 6B.
If a general unstructured multiplicative uncertainty (namely, a full Δ matrix) had been included in the model, an analysis of performance would have produced conservative results and some robust therapies might have been discarded. In fact, this definition of uncertainty is not directly connected to the physical phenomena occurring in the system and will generally include physically unfeasible perturbations. Defining ranges for parameter values includes a structured uncertainty in the model, preserving a closer physical interpretation, related to the stochastic noise and the experimental error inherent to biological networks.
Several extensions of SSV analysis exists which can be invaluable to drug discovery. As parameter fitting can be computationally expensive, by reversing the idea of robust performance and searching for target combinations which most easily destabilize the diseased state, the number of parameter fittings performed early in the analysis can be greatly reduced . Furthermore, in this paper we considered a single-input, single-output system with uncertainty being limited to the interactions (parameters) of the network. In multiple input and multiple output (MIMO) networks, more complex performance envelopes can be considered during robust performance analysis. The analysis can also be extended to include other clinically-interesting sources of uncertainty, such as dosages issues, blood clearance, etc.
The complexity and size of biological systems make observation-based approaches to combinatorial drug therapy discovery prohibitive due to the associated financial burden and time requirements. Many companies are now aware of the value of using in silico techniques to guide discovery, but these analyses may rely too heavily on model accuracy. Using tools such as SSV analysis, biological networks can be screened for MDTs that are robust to various uncertainties. These uncertainties may be noise experienced in data, neglected dynamics, or even intrinsic biological variability. Furthermore, the performance of the network can be user-defined to cover several drug-related concerns such as drug efficacy and known potential side effects. MDTs identified by SSV analysis are robust to all model hypotheses expressed in the uncertainty description, and are more likely to be effective during experimentation. In conclusion, SSV analysis can prioritize target combinations by quantifying treatment efficacy given uncertainty in a systematic way.
All the numerical simulations were performed in MATLAB 7.7.0 R2008b (MathWorks, Inc.). The healthy and diseased steady-states were calculated using Matlab function fsolve to solve the system of equations in Figure 2B after equating to zero the time derivative terms. The nonlinear ODE model in Figure 2B was solved using Matlab function ode15s.
Stochastic Simulation Algorithm (SSA) and performance envelopes
The nonlinear model was expressed in terms of number of molecules, instead of concentrations, and SSA was applied according to . The initial conditions for the results shown in Figure 3A were the healthy and diseased steady-states, respectively. 100 trajectories were generated and interpolated at 100 regular time points. At each time point mean and standard deviation values among the trajectories were calculated.
where y ss,h is the healthy output at steady-state, s mean is the time-averaged standard deviation from the stochastic simulations, and f is a weighting factor chosen to have the stochastic envelopes reasonably contained in the performance ones. A value of f equal to 1.6 was used. The diseased performance envelope shown in Figure 3A was calculated analogously.
Derivation of potential therapies
where t i are 20 regularly-spaced time points in the simulated time span. The initial conditions were given by the diseased steady-state. Up to 4 parameters were simultaneously allowed to change in the range ±100% of their nominal value.
Test for nominal performance
The model of each therapy was run using the healthy steady-state as initial condition. The absolute deviation of the therapy model output from the healthy one was calculated, and the therapies that met the condition in (2) at each time point, t i , defined above, were selected as respondent to the requirements of nominal performance.
where , , , and A, B, C, and D are constant matrices. As A, B, C, and D depend on the parameter values, they are different for each therapeutically-treated diseased model and for the healthy one. Furthermore, the therapy models have input , while the healthy one has null input in deviation variables.
Model rearrangement and SSV analysis
An uncertainty r k = 45% was applied to each parameter in the model, according to the definition in (3). Multiple deviation models were then defined as the difference between each therapy model and the healthy one in the form described by (9). In practice these deviation models were obtained numerically in Matlab. Their output was normalized by the performance weighting factor, w p , defined in (3). The models were numerically converted into an M-Δ form by using the Robust Control Toolbox in Matlab.
SSV analysis was applied to all the deviation models, after a nominal stability check (data not shown). μ RP was calculated using Matlab function mussv.
CL, JES, and FJD were supported by Pfizer Inc. through Contract No. DFP01, and the Institute for Collaborative Biotechnologies through Grant No. W911NF-09-D-0001 from the U.S. Army Research Office. KRS and LRP were funded by Pfizer Inc. and the Institute for Collaborative Biotechnologies under Grant No. DFR3A-8-447850-23002 from the Army Research Office. LRP was also supported by Grant R01EB007511 from the National Institute of Biomedical Imaging and Bioengineering and DOE Contract No. DE-FG02-04ER25621. KRS was also supported by a National Science Foundation Graduate Research Fellowship.
- Quackenbush J: Extracting biology from high-dimensional biological data. J Exp Biol. 2007, 210: 1507-1517. 10.1242/jeb.004432View ArticlePubMedGoogle Scholar
- Kitano H: Biological robustness. Nat Rev Genet. 2004, 5: 826-837. 10.1038/nrg1471View ArticlePubMedGoogle Scholar
- Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology. Nature. 1999, 402: C47-C52. 10.1038/35011540View ArticlePubMedGoogle Scholar
- Hornberg JJ, Bruggeman FJ, Binder B, Geest CR, de Vaate A, Lankelma J, Heinrich R, Westerhoff HV: Principles behind the multifarious control of signal transduction - ERK phosphorylation and kinase/phosphatase control. FEBS J. 2005, 272: 244-258. 10.1111/j.1432-1033.2004.04404.xView ArticlePubMedGoogle Scholar
- Barabasi AL, Oltvai ZN: Network biology: Understanding the cell's functional organization. Nat Rev Genet. 2004, 5: 101-U115. 10.1038/nrg1272View ArticlePubMedGoogle Scholar
- Stelling J, Sauer U, Szallasi Z, Doyle FJ, Doyle J: Robustness of cellular functions. Cell. 2004, 118: 675-685. 10.1016/j.cell.2004.09.008View ArticlePubMedGoogle Scholar
- Schadt EE, Friend SH, Shaywitz DA: A network view of disease and compound screening. Nat Rev Drug Discov. 2009, 8: 286-295. 10.1038/nrd2826View ArticlePubMedGoogle Scholar
- Kola I, Landis J: Can the pharmaceutical industry reduce attrition rates?. Nat Rev Drug Discov. 2004, 3: 711-715. 10.1038/nrd1470View ArticlePubMedGoogle Scholar
- Sams-Dodd F: Target-based drug discovery: is something wrong?. Drug Discov Today. 2005, 10: 139-147. 10.1016/S1359-6446(04)03316-1View ArticlePubMedGoogle Scholar
- Durrant JD, Amaro RE, Xie L, Urbaniak MD, Ferguson MAJ, Haapalainen A, Chen ZJ, Di Guilmi AM, Wunder F, Bourne PE, McCammon JA: A Multidimensional Strategy to Detect Polypharmacological Targets in the Absence of Structural and Sequence Homology. PLoS Comput Biol. 2010, 6:Google Scholar
- Zimmermann GR, Lehar J, Keith CT: Multi-target therapeutics: when the whole is greater than the sum of the parts. Drug Discov Today. 2007, 12: 34-42. 10.1016/j.drudis.2006.11.008View ArticlePubMedGoogle Scholar
- Csermely P, Agoston V, Pongor S: The efficiency of multi-target drugs: the network approach might help drug design. Trends Pharmacol Sci. 2005, 26: 178-182. 10.1016/j.tips.2005.02.007View ArticlePubMedGoogle Scholar
- Lehar J, Krueger A, Zimmermann G, Borisy A: High-order combination effects and biological robustness. Mol Syst Biol. 2008, 4: 215- 10.1038/msb.2008.51PubMed CentralView ArticlePubMedGoogle Scholar
- Loscalzo J, Kohane I, Barabasi AL: Human disease classification in the postgenomic era: A complex systems approach to human pathobiology. Mol Syst Biol. 2007, 3: 124- 10.1038/msb4100163PubMed CentralView ArticlePubMedGoogle Scholar
- Agoston V, Csermely P, Pongor S: Multiple weak hits confuse complex systems: A transcriptional regulatory network as an example. Phys Rev E. 2005, 71: 10.1103/PhysRevE.71.051909.Google Scholar
- Araujo RP, Liotta LA, Petricoin EF: Proteins, drug targets and the mechanisms they control: the simple truth about complex networks. Nat Rev Drug Discov. 2007, 6: 871-880. 10.1038/nrd2381View ArticlePubMedGoogle Scholar
- Bonhoeffer S, May RM, Shaw GM, Nowak MA: Virus dynamics and drug therapy. Proc Natl Acad Sci USA. 1997, 94: 6971-6976. 10.1073/pnas.94.13.6971PubMed CentralView ArticlePubMedGoogle Scholar
- Fisher RI, Gaynor ER, Dahlberg S, Oken MM, Grogan TM, Mize EM, Glick JH, Coltman CA, Miller TP: Comparison of a standard regimen (CHOP) with 3 intensive chemotherapy regimens for advanced non-hodgkins-lymphoma. N Engl J Med. 1993, 328: 1002-1006. 10.1056/NEJM199304083281404View ArticlePubMedGoogle Scholar
- Stein GE, Gurwith MJ: Amoxicillin potassium clavulanate, a beta-lactamase-resistant antibiotic combination. Clin Pharm. 1984, 3: 591-599.PubMedGoogle Scholar
- Schoeberl B, Eichler-Jonsson C, Gilles ED, Muller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002, 20: 370-375. 10.1038/nbt0402-370View ArticlePubMedGoogle Scholar
- Tyson JJ, Albert R, Goldbeter A, Ruoff P, Sible J: Biological switches and clocks. J R Soc Interface. 2008, 5: S1-S8. 10.1098/rsif.2008.0179.focusPubMed CentralView ArticlePubMedGoogle Scholar
- Cho CR, Labow M, Reinhardt M, van Oostrum J, Peitsch MC: The application of systems biology to drug discovery. Curr Opin Chem Biol. 2006, 10: 294-302. 10.1016/j.cbpa.2006.06.025View ArticlePubMedGoogle Scholar
- Fridlyand LE, Harbeck MC, Roe MW, Philipson LH: Regulation of cAMP dynamics by Ca2+ and G protein-coupled receptors in the pancreatic beta-cell: a computational approach. American Journal of Physiology-Cell Physiology. 2007, 293: C1924-C1933. 10.1152/ajpcell.00555.2006View ArticlePubMedGoogle Scholar
- Kim J, Saidel GM, Kalhan SC: A computational model of adipose tissue metabolism: Evidence for intracellular compartmentation and differential activation of lipases. J Theor Biol. 2008, 251: 523-540. 10.1016/j.jtbi.2007.12.005PubMed CentralView ArticlePubMedGoogle Scholar
- Sedaghat AR, Sherman A, Quon MJ: A mathematical model of metabolic insulin signaling pathways. Am J Physiol Endocrinol Metab. 2002, 283: E1084-E1101.View ArticlePubMedGoogle Scholar
- Kreeger PK, Lauffenburger DA: Cancer systems biology: a network modeling perspective. Carcinogenesis. 2009, 31: 2-8. 10.1093/carcin/bgp261PubMed CentralView ArticlePubMedGoogle Scholar
- Lalonde RL, Kowalski KG, Hutmacher MM, Ewy W, Nichols DJ, Milligan PA, Corrigan BW, Lockwood PA, Marshall SA, Benincosa LJ, et al.: Model-based drug development. Clin Pharmacol Ther. 2007, 82: 21-32. 10.1038/sj.clpt.6100235View ArticlePubMedGoogle Scholar
- Michelson S, Sehgal A, Friedrich C: In silico prediction of clinical efficacy. Curr Opin Biotechnol. 2006, 17: 666-670. 10.1016/j.copbio.2006.09.004View ArticlePubMedGoogle Scholar
- Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, Cornish-Bowden A, et al.: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19: 524-531. 10.1093/bioinformatics/btg015View ArticlePubMedGoogle Scholar
- Alterovitz G, Muso T, Ramoni MF: The challenges of informatics in synthetic biology: from biomolecular networks to artificial organisms. Briefings in Bioinformatics. 2010, 11: 80-95. 10.1093/bib/bbp054PubMed CentralView ArticlePubMedGoogle Scholar
- Doyle J: Analysis of feedback-systems with structured uncertainties. IEE Proc-D. 1982, 129: 242-250.View ArticleGoogle Scholar
- Jacobsen EW, Cedersund G: Structural robustness of biochemical network models-with application to the oscillatory metabolism of activated neutrophils. Iet Systems Biology. 2008, 2: 39-47. 10.1049/iet-syb:20070008View ArticlePubMedGoogle Scholar
- Ma L, Iglesias PA: Quantifying robustness of biochemical network models. BMC Bioinformatics. 2002, 3:Google Scholar
- Shoemaker JE, Doyle FJ: Identifying fragilities in biochemical networks: Robust performance analysis of Fas signaling-induced apoptosis. Biophys J. 2008, 95: 2610-2623. 10.1529/biophysj.107.123398PubMed CentralView ArticlePubMedGoogle Scholar
- Ghaemi R, Sun J, Iglesias PA, Del Vecchio D: A method for determining the robustness of bio-molecular oscillator models. Bmc Systems Biology. 2009, 3:Google Scholar
- Brandman O, Meyer T: Feedback loops shape cellular signals in space and time. Science. 2008, 322: 390-395. 10.1126/science.1160617PubMed CentralView ArticlePubMedGoogle Scholar
- Gillespie DT: General method for numerically simulating stochastic time evolution of coupled chemical-reactions. J Comput Phys. 1976, 22: 403-434. 10.1016/0021-9991(76)90041-3.View ArticleGoogle Scholar
- Skogestad S, Postlethwaite I: Multivariable feedback control: analysis and design. 2005, Chichester, West Sussex, England: John Wiley & Sons, 2Google Scholar
- Bates D, Postlethwaite I: Robust Multivariable Control of Aerospace Systems. 2002, Ios Pr IncGoogle Scholar
- Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S: Global Sensitivity Analysis. 2008, Chichester, West Sussex, England: John Wiley & SonsGoogle Scholar