Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Methodology article

Confidence from uncertainty - A multi-target drug screening method from robust control theory

Camilla Luni1, Jason E Shoemaker1, Kevin R Sanft2, Linda R Petzold2 and Francis J Doyle1*

Author Affiliations

1 Department of Chemical Engineering, University of California, Santa Barbara, CA 93106-5080, USA

2 Department of Computer Science, University of California, Santa Barbara, CA 93106-5070, USA

For all author emails, please log on.

BMC Systems Biology 2010, 4:161  doi:10.1186/1752-0509-4-161


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/4/161


Received:26 May 2010
Accepted:24 November 2010
Published:24 November 2010

© 2010 Luni et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

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.

Results

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.

Conclusions

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.

Background

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 [1]. Some general features of biological networks have been explored computationally, such as robustness [2], modularity [3], control coefficients [4], and connectivity properties [5]. 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 [2]. 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 [7]. 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 [7]. 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 [17]. Oncological chemotherapeutic regimens often involve the combination of cyclophosphamide, hydroxydaunorubicin, oncovin, and prednisone, abbreviated as CHOP [18]. Augmentin, an amoxicillin-based antibiotic, contains clavulanic acid to inhibit a known mechanism of amoxicillin degradation [19]. In comparison to their single-perturbation counterparts, these MDTs often show an order of magnitude greater efficacy [17]. 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 [20].

Optimizing drug combinations and concentrations produces an unmanageable number of possible therapies to explore, demanding efficient computational methods of screening [11]. 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 [21].

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) [22], 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 [29], and paving the way for a biochemically detailed artificial organism reconstruction [30]. 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 [11]. 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.

We propose the use of Structured Singular Value (SSV) analysis as a powerful tool for drug target discrimination in biological models also accounting for uncertainty. This technique was developed in the control theory field [31], but has already been successively applied in the analysis of biological systems [32-35]. In the proposed methodology, SSV allows the discrimination of highly effective MDTs from a large pool of potential candidates, according to the robust response of the diseased network in the face of known or inferred sources of uncertainty (Figure 1). For illustration purposes, we explain the methodology through a case study, given by a negative feedback network motif. Moreover, we discuss strengths, limitations, and extensions (Figure 1) of the proposed method and its application, with respect to other existing ones.

thumbnailFigure 1. Range of applicability of SSV analysis for robust therapy design. Once a model has been established that satisfactorily explains the dynamics of the diseased state, SSV analysis can be used to identify potent and robust multi-drug therapy candidates. SSV analysis first identifies which therapies can best manipulate the protein(s) of interest. Then, the candidate list is further filtered to therapies which are robust to known or perceived uncertainty affecting the treatment. The uncertainty may include parameter uncertainty and uncertainty generated during model development, but also disturbances occurring during the actual treatment, such as failure to properly adhere to a drug regimen schedule.

Results

Case study description

A schematic description of the case study used to illustrate the proposed methodology is shown in Figure 2A. The component X is converted to Y through an enzymatic reaction, catalyzed by U, that includes the intermediate production of the complex UX. The production of X and the degradation of X and Y are also considered. All reaction rates are modelled by mass action. The product Y regulates its own production via autoinhibition. This negative feedback mechanism is modelled as a multiplicative factor dependent on the concentration of Y. Negative feedback is a widespread strategy in biological networks that strongly contributes to their spatial and temporal complexity [36]. The equation system is shown in Figure 2B, and the arbitrary set of nominal parameter values are provided in the figure caption.

thumbnailFigure 2. Case study model. (A) Topology of the network. U is the enzyme catalyzing the conversion of X to product Y, through the intermediate UX. The open arrows indicate the chemical reactions, and the oval arrow a negative regulation. k's represent the parameters involved in each step. Ø is the null component to indicate production and degradation. Input to the system is given by the total enzyme, Utot, concentration, constant over time and given by utot = u+ux (lower-case component names indicate the corresponding concentrations). Output of the system is y. Inputs and outputs are highlighted in red. (B) Nonlinear model equations. The reaction rates are given by mass-action, negative feedback is described by the multiplicative term containing k5. The nominal values of the parameters are: k1 = 1, k2 = 2, k3 = 10, k4 = 0.5, k5 = 0.5, k-1 = 3, k-3 = 1.

The first requirement for any drug investigation is to identify the appropriate inputs and outputs of the system. These choices depend on which components (such as cytokine concentration, mRNA level, marker expression, etc.) are significant for defining the healthy and diseased states, and are measurable by available biological assays. This example analyzes a single-input single-output system. We assume the input is the total concentration of enzyme, utot, constant over time. A disease state emerges when the input has a high concentration, utot,d, compared to the basal healthy state, utot,h. As a consequence, the output concentration, y, is up-regulated in this condition (Figure 3A). The goal is to re-parameterize the diseased model to obtain a therapeutically treated model that, with a diseased input utot,d, allows recovery from the diseased output to the healthy one, even in presence of uncertainty.

thumbnailFigure 3. Performance envelopes and therapeutic fitting. (A) Temporal simulation of the nonlinear model in Figure 2B with nominal parameters, under healthy conditions (utot,h = 0.2, and initial conditions given by the healthy steady-state), and diseased conditions (utot,d = 2, and initial conditions given by the diseased steady-state). Performance envelopes are generated to contain the stochastic envelopes, resulted by the Stochastic Simulation Algorithm (mean ± standard deviation). Nominal results are also shown in dashed lines. (B) Comparison between the performance envelopes and the results obtained from the nonlinear model, starting from the diseased steady-state, with parameter values modified according to the 56 therapies (blue curves).

Healthy performance and potential therapies

Due to biological variability, the healthy performance is given by an envelope that defines upper, yub, and lower, ylb, bounds on the output, rather than an idealized, nominal single trajectory (Figure 3A). Thus, the system meets the requirements for healthy performance when:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M1">View MathML</a>

(1)

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 [37]. 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, utot,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, utot,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. <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M2">View MathML</a> combinations of the n = 6 parameters, were obtained. A comparison between the outcome, ydt, 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

According to the definition in (1), the criterion for nominal performance requires that the output of a therapeutic model does not cross the boundaries of the healthy performance envelope, when using the healthy steady-state as initial condition. It is formally convenient to re-formulate the problem defining an upper bound for the absolute deviation of the therapy output from the healthy one:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M3">View MathML</a>

(2)

We applied this preliminary screening method, based on direct trajectories comparison, to our case study. A total of 41 therapies, out of 56 potential, were selected as giving a healthy nominal performance in the presence of a diseased input, utot,d, i.e., their output without any parametric uncertainty falls within the healthy performance envelope (Figure 4).

thumbnailFigure 4. Nominal performance analysis results. Trajectories obtained with the nonlinear deviation model for each of the 56 therapies without parametric uncertainty. Green and red lines denote therapies that pass and do not pass the nominal performance selection criterion, respectively. The performance envelope is shaded in light green.

Uncertainty description and robust performance

A mathematical approximation of a complex biophysical system must account for multiple sources of uncertainty, due to stochastic noise, experimental error, or other possible fluctuations induced by the interaction of the system with its surrounding. A confidence interval can be assigned to each parameter, during the procedure of experimental data fitting, as a lumped measure of these multiple sources of uncertainty. Thus, each parameter in the model is represented in the following form:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M4">View MathML</a>

(3)

where k ∈ [kmin, kmax] is a generic parameter of the model, kmean = (kmin + kmax)/2, rk = (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., rk = 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 [31]. We provide here an intuitive understanding of how it works, and refer to textbooks in the field for a more technical explanation [38]. 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.

We defined a deviation model as the difference between a therapy model and the healthy one, and we normalized the output by a weighting factor, wp, representing the performance specification. The criterion for robust performance can now be expressed in terms of ratio between the normalized deviation model output, (ydt - yh)/wp, and its input,(utot,d - utot,h), i.e.:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M5">View MathML</a>

(4)

The performance weighting factor, wp, is related to the performance envelope bounds by the following relationship:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M6">View MathML</a>

(5)

Now the model includes parametric uncertainty. It is always possible, through a linear fractional transformation (LFT), to pull out the uncertain elements from the nominal model, and to recast it in a M-Δ form (Figure 5), where M is the matrix describing the nominal model and Δ the matrix containing all the uncertainty, namely the δk's defined above, and the performance specification. The Δ-matrix has a particular structure, due to the presence of a number of zeros in some positions, dependent on the specific starting model. As |δk| ≤ 1 and the output is normalized by wp, Δ is also normalized. From a control theory standpoint, putting the system in this form converts the performance analysis problem to the study of the stability of the loop in Figure 5.

thumbnailFigure 5. Block diagram representation of the deviation models. The deviation model is shown in M-Δ form. The vector of input and output between the two blocks are also indicated. uΔ and yΔ represent the uncertain components of the input and output, respectively, for the system M.

General aspects on SSV

SSV is a worst-case analysis that excludes a therapy if, even for a single parametric combination within the defined uncertainty ranges, it fails to meet the performance specifications defined by the envelope. It is based on the calculation of the structured singular value, μRP (where RP stands for robust performance), by solving the following minimization problem:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M7">View MathML</a>

(6)

where I is the identity matrix, and km a scalar factor. A result well-known in control theory, simplistically stated here, is that, when det(I - ) = 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(km) 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 [38], and algorithms are also included in technical software, like Matlab.

Selection of therapies for robust performance by SSV

Values of μRP for the 41 therapies are shown in Figure 6A. Only 5 therapies have μRP < 1 and passed this screening test. Table 1 summarizes the parameters involved in each one. Interestingly, no single-parameter therapies met the robust performance specification. This is a confirmation of the importance of a MDT approach, as opposed to a drug strategy having only one target of intervention. In fact, because of the interconnected structure of the network, a robust therapy was obtained by drug (i.e. parameter) combinations affecting at the same time processes with the same effect but dislocated in different points of the network. For example, therapy 11 increases X concentration by reducing its consumption in the forward reaction involving parameter k1, and by reducing its degradation dependent on parameter k-3. Therapy 12 instead decreases Y concentration by reducing its production from UX (parameter k2) and increasing its degradation (parameter k4).

thumbnailFigure 6. Robust performance analysis and results of one robust therapy in presence of parametric uncertainty. (A) Results of the SSV analysis applied for robust performance, in presence of parametric uncertainty. Green and red dots illustrate therapies that pass and do not pass the robust performance selection criterion, respectively. (B) Comparison between the performance envelopes described in Figure 3A and the results obtained from therapy no. 30. Blue curves are the simulation results by therapy 30 linearized model with 100 different parameter sets, sampled within ± 45% of the nominal values. Gray curves are the upper and lower bounds of the stochastic envelope generated by the Stochastic Simulation Algorithm of therapy 30 nonlinear model (mean ± standard deviation of 100 trajectories).

Table 1. 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.

Discussion

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 [39]. 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.

Sensitivity analysis is a possible alternative approach to identify parameters for therapy design. For comparison, the local sensitivity analysis (LSA) of the case study is reported in Figure 7. Sensitivity coefficients with respect to only one parameter were calculated, even if global sensitivity analysis methods exist [40]. Parameter k3 was not included, as it is not present in the linearized version of the model used in the SSV analysis. LSA results show the disease state is particularly sensitive to parameters k2 and k4, those involved in the single target therapies 2 and 3. Therapy 2 resulted in non-robust performance (Figure 6A), and therapy 3 failed even the nominal performance test, when parametric uncertainty is not accounted for (results not shown). The discrepancy in the results is due to the local character of LSA. Even when the model is linear, LSA is unable to account for system behaviour in the presence of large perturbations. Moreover, as it does not include a definition of performance, once sensitive parameters are identified, the information on the amount of perturbation needed to restore a healthy performance is not available.

thumbnailFigure 7. Local sensitivity analysis results. Sensitivity coefficients, Si, of the diseased output at steady-state, yss,d, respect to the parameter indicated on the x-axis. The coefficients are normalized with the nominal value of each parameter and with the steady-state output concentration yss,d.

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 [34]. 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.

Conclusions

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.

Methods

Nominal model

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 [37]. 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.

The upper, yub, and lower, ylb, bounds of the healthy performance envelope shown in Figure 3A were then calculated by:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M8">View MathML</a>

(7)

where yss,h is the healthy output at steady-state, smean 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

Each therapy is obtained by fitting the output of the nonlinear model with diseased input, utot,d, to the healthy output response. The fitting was performed using Matlab function fmincon by minimizing the following cost function, C:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M9">View MathML</a>

(8)

where ti 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, ti, defined above, were selected as respondent to the requirements of nominal performance.

Model linearization

The model of each therapy was transformed by analytical Jacobian linearization around the healthy steady-state. The following deviation variables, indicated by the over bar, were defined: <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M10">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M11">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M12">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M13">View MathML</a>. Then, the model was rearranged in state-space form:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M14">View MathML</a>

(9)

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M15">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M16">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M17">View MathML</a>, 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 <a onClick="popup('http://www.biomedcentral.com/1752-0509/4/161/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/4/161/mathml/M18">View MathML</a>, while the healthy one has null input in deviation variables.

Model rearrangement and SSV analysis

An uncertainty rk = 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, wp, 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.

Authors' contributions

CL participated in the design of the study, carried out the numerical simulations and drafted the manuscript. JES participated in the design of the study and helped to draft the manuscript. KRS helped with preliminary stochastic analyses in the design of the manuscript. LRP participated in the design of the study and helped to draft the manuscript. FJD participated in the design of the study and helped to draft the manuscript. All authors read and approved the final manuscript.

Acknowledgements

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.

References

  1. Quackenbush J: Extracting biology from high-dimensional biological data.

    J Exp Biol 2007, 210:1507-1517. PubMed Abstract | Publisher Full Text OpenURL

  2. Kitano H: Biological robustness.

    Nat Rev Genet 2004, 5:826-837. PubMed Abstract | Publisher Full Text OpenURL

  3. Hartwell LH, Hopfield JJ, Leibler S, Murray AW: From molecular to modular cell biology.

    Nature 1999, 402:C47-C52. PubMed Abstract | Publisher Full Text OpenURL

  4. 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. PubMed Abstract | Publisher Full Text OpenURL

  5. Barabasi AL, Oltvai ZN: Network biology: Understanding the cell's functional organization.

    Nat Rev Genet 2004, 5:101-U115. PubMed Abstract | Publisher Full Text OpenURL

  6. Stelling J, Sauer U, Szallasi Z, Doyle FJ III, Doyle J: Robustness of cellular functions.

    Cell 2004, 118:675-685. PubMed Abstract | Publisher Full Text OpenURL

  7. Schadt EE, Friend SH, Shaywitz DA: A network view of disease and compound screening.

    Nat Rev Drug Discov 2009, 8:286-295. PubMed Abstract | Publisher Full Text OpenURL

  8. Kola I, Landis J: Can the pharmaceutical industry reduce attrition rates?

    Nat Rev Drug Discov 2004, 3:711-715. PubMed Abstract | Publisher Full Text OpenURL

  9. Sams-Dodd F: Target-based drug discovery: is something wrong?

    Drug Discov Today 2005, 10:139-147. PubMed Abstract | Publisher Full Text OpenURL

  10. 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 PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. 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. PubMed Abstract | Publisher Full Text OpenURL

  12. 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. PubMed Abstract | Publisher Full Text OpenURL

  13. Lehar J, Krueger A, Zimmermann G, Borisy A: High-order combination effects and biological robustness.

    Mol Syst Biol 2008, 4:215. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  15. Agoston V, Csermely P, Pongor S: Multiple weak hits confuse complex systems: A transcriptional regulatory network as an example.

    Phys Rev E 2005., 71 Publisher Full Text OpenURL

  16. 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. PubMed Abstract | Publisher Full Text OpenURL

  17. Bonhoeffer S, May RM, Shaw GM, Nowak MA: Virus dynamics and drug therapy.

    Proc Natl Acad Sci USA 1997, 94:6971-6976. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. 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. PubMed Abstract | Publisher Full Text OpenURL

  19. Stein GE, Gurwith MJ: Amoxicillin potassium clavulanate, a beta-lactamase-resistant antibiotic combination.

    Clin Pharm 1984, 3:591-599. PubMed Abstract OpenURL

  20. 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. PubMed Abstract | Publisher Full Text OpenURL

  21. Tyson JJ, Albert R, Goldbeter A, Ruoff P, Sible J: Biological switches and clocks.

    J R Soc Interface 2008, 5:S1-S8. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. 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. PubMed Abstract | Publisher Full Text OpenURL

  23. 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. PubMed Abstract | Publisher Full Text OpenURL

  24. 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. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Sedaghat AR, Sherman A, Quon MJ: A mathematical model of metabolic insulin signaling pathways.

    Am J Physiol Endocrinol Metab 2002, 283:E1084-E1101. PubMed Abstract | Publisher Full Text OpenURL

  26. Kreeger PK, Lauffenburger DA: Cancer systems biology: a network modeling perspective.

    Carcinogenesis 2009, 31:2-8. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. 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. PubMed Abstract | Publisher Full Text OpenURL

  28. Michelson S, Sehgal A, Friedrich C: In silico prediction of clinical efficacy.

    Curr Opin Biotechnol 2006, 17:666-670. PubMed Abstract | Publisher Full Text OpenURL

  29. 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. PubMed Abstract | Publisher Full Text OpenURL

  30. 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. PubMed Abstract | Publisher Full Text OpenURL

  31. Doyle J: Analysis of feedback-systems with structured uncertainties.

    IEE Proc-D 1982, 129:242-250. OpenURL

  32. 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. PubMed Abstract | Publisher Full Text OpenURL

  33. Ma L, Iglesias PA: Quantifying robustness of biochemical network models.

    BMC Bioinformatics 2002., 3 PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  34. Shoemaker JE, Doyle FJ III: Identifying fragilities in biochemical networks: Robust performance analysis of Fas signaling-induced apoptosis.

    Biophys J 2008, 95:2610-2623. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  35. 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 PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  36. Brandman O, Meyer T: Feedback loops shape cellular signals in space and time.

    Science 2008, 322:390-395. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  37. Gillespie DT: General method for numerically simulating stochastic time evolution of coupled chemical-reactions.

    J Comput Phys 1976, 22:403-434. Publisher Full Text OpenURL

  38. Skogestad S, Postlethwaite I: Multivariable feedback control: analysis and design. 2nd edition. Chichester, West Sussex, England: John Wiley & Sons; 2005.

  39. Bates D, Postlethwaite I: Robust Multivariable Control of Aerospace Systems. Ios Pr Inc; 2002.

  40. Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S: Global Sensitivity Analysis. Chichester, West Sussex, England: John Wiley & Sons; 2008.