Conditional robustness analysis for fragility discovery and target identification in biochemical networks and in cancer systems biology
- Fortunato Bianconi^{1}Email author,
- Elisa Baldelli^{2},
- Vienna Luovini^{3},
- Emanuel F. Petricoin^{2},
- Lucio Crinò^{3} and
- Paolo Valigi^{4}
https://doi.org/10.1186/s12918-015-0216-5
© Bianconi et al. 2015
Received: 4 February 2015
Accepted: 16 July 2015
Published: 19 October 2015
The Erratum to this article has been published in BMC Systems Biology 2016 10:24
Abstract
Background
The study of cancer therapy is a key issue in the field of oncology research and the development of target therapies is one of the main problems currently under investigation. This is particularly relevant in different types of tumor where traditional chemotherapy approaches often fail, such as lung cancer.
Results
We started from the general definition of robustness introduced by Kitano and applied it to the analysis of dynamical biochemical networks, proposing a new algorithm based on moment independent analysis of input/output uncertainty. The framework utilizes novel computational methods which enable evaluating the model fragility with respect to quantitative performance measures and parameters such as reaction rate constants and initial conditions. The algorithm generates a small subset of parameters that can be used to act on complex networks and to obtain the desired behaviors. We have applied the proposed framework to the EGFR-IGF1R signal transduction network, a crucial pathway in lung cancer, as an example of Cancer Systems Biology application in drug discovery. Furthermore, we have tested our framework on a pulse generator network as an example of Synthetic Biology application, thus proving the suitability of our methodology to the characterization of the input/output synthetic circuits.
Conclusions
The achieved results are of immediate practical application in computational biology, and while we demonstrate their use in two specific examples, they can in fact be used to study a wider class of biological systems.
Keywords
Background
Most diseases, including cancer, involve a large number and variety of elements that interact via complex networks and, consequently, display highly nonlinear behavior. Traditional approaches to the study of biological phenomena consider single events, such as single mutations, single gene or protein alterations, and their corresponding effects on the biological phenotype. These reductionist approaches are intrinsically unable to fully capture the overall complexity of molecular interaction networks. The key main of systems biology is to take into account all the interactions within a given network, as emphasized by the system notion [1, 2].
The identification of specific molecular targets that play a central role in cancer cell proliferation and survival has led to the development of a targeted therapy approach for the treatment of cancer patients in the clinical setting. Nevertheless, knocking out one target molecule in a biochemical pathway may not be enough to treat a disease such as cancer, because tumor cells often find alternative molecular routes to escape the drug-induced blockage. This is one reason why current drug design strategies fail in some cases [3].
In the last few years, Systems Biology has received increasing attention as a promising approach toward personalized medicine and to assist the oncologist community [4]. Through the Systems Biology approach, for example, it could be possible to improve our understanding of the complex signalling networks involved in cancer. This methodology would allow for the development of smarter therapeutic strategies; e.g., disrupting simultaneously two or three key intersections crucial for tumor growth and progression in a biochemical network. This approach could lead to significant advances in the treatment of cancer and help in transforming traditional reductionism-based methods into systems-level ones for drug discovery [5–8].
The analysis of complete (or, at least, a large portion of) networks could be a way to understand the robustness property of cancer [9, 10]. In this study, we focus on cancer robustness as a quantitative measurement indicating the cells ability to maintain their functions against internal and external perturbations. In the framework of cancer research, it is relevant to discover how to reduce robustness of cell proliferation attitude; i.e., it is important to study the fragility of cancer cells and to discover ways to increase this fragility.
The proliferation activity of normal and tumor cells can be evaluated by looking at the activation of downstream proteins; e.g., phosphorylated form of extracellular signal-regulated kinase (ERK) [11]. To achieve quantitative measure of the proliferation activity, a suitable proliferation indicator can be used. For example, the intensity of phosphorylated ERK is one such indicator. If we describe the proliferation activity of a cell population by means of the probability density for the chosen proliferation indicator, the expected mean value and variance for cancer cells will be higher than normal cells. The plots at the bottom of Fig. 1 show an example of probability density of a proliferation indicator in tumor (red line) and normal cells (green line), respectively. Additional variability in the signaling networks, such as different topologies and settings, arises depending on cell type (e.g., lung, breast, colon,..) and their characteristics (normal vs cancer cell). The variability on drugs response can be seen as example [12].
To prove the suitability of our approach in a more general context, we investigated both an oncological network and a synthetic biology example. Synthetic biology aims to provide a way to synthesize living systems according to certain design specifications by means of strategies similar to the ones electrical engineers employ to construct electronic circuits[13]. Among the proposed applications of synthetic biology, there are medical applications, such as synthesis of biofuels, biosensing, selective destruction of cancerous cellular tissue and low-cost production of drugs [14]. At the heart of this new field are synthetic gene networks. Such circuits are frequently constructed in Escherichia coli by cutting and combining natural or engineered coding DNA regions and promoters. Widely studied examples of such circuits include a toggle switch [14], an auto-repressed circuit [15], an oscillating “repressilator” [16], as well as other electronics-inspired genetic devices, including pulse generators [17], digital logic gates, filters and communication modules [14].
Biochemical signal transduction networks, both natural and synthetic, can be modeled with a large spectrum of mathematical tools; here we focus on Ordinary Differential Equations (ODE) that allow us to describe the time evolution of a set of proteins of interest. Once the pathway structure is drawn, the corresponding equations are relatively easy to write down using widely accepted kinetic laws, such as the law of mass action, the Michaelis-Menten law or the Hill functions. The generated model will depend on several parameters, and the corresponding identification and model selection problems are relevant issues (see [18, 19] and the references therein).
The alterations giving rise to tumor development can be described by perturbations on the parameters characterizing the biochemical reactions; i.e., the parameters characterizing the whole set of ODE [20]. The lack of parameter identifiability in large-scale network models hamper translation of the results of modelling studies into the process of anti-cancer drug development. Hence, assuming large perturbations on those parameters and studying the corresponding Global Sensitivity Analysis (GSA) is a way to analyze the uncertainty of the model parameters and to generate valid predictions on parametric sensitivities [21, 22]. In [23] the authors applied GSA to an ErbB signalling network model with the goal of exploring how multi-parametric network perturbations affect signal propagation through cancer-related networks.
Other approaches to robustness analysis are those searching for the shape and volume of the region in the parameter space in which the system is functioning properly [24], and for the topology and geometry of such a region, which turns out to have important consequences for the robustness [25–27]. The topological properties of the networks and their relationship with robustness and fragility in large-scale bio-molecular networks have been studied [28], showing that networks with a larger number of positive feedback loops and a smaller number of negative feedback cycles are likely to be more robust against disturbances.
A related problem is model validation. In this case, a robustness index can be seen as a measure of plausibility in models of biochemical networks [29]. Along this line, [30] introduces the concept of “glocal” robustness analysis as a combination of global and local tools used for model validation and for identifying key causes of high or low robustness. A possibly related problem, studied in [31–33], is the core prediction.
The purpose of the present contribution arises from personalized therapy and is focused on computational analysis. We use the in silico model to generate samples of the proliferation indicator and some computational elaborations to select a small number of nodes in the cancer cells signaling network that reduces their proliferation robustness; i.e., that improves the fragility of cancer cell. It follows that these nodes are candidate drug targets and our algorithm permits us to discover how these nodes can be conditioned to obtain the desired behavior.
The rest of the paper is organized as follows. In the “Methods” section we introduce the complete theory associated with our procedure, namely the robustness problem, the conditional robustness concept, its use in the solution of the fragility problem, and the proposed computational algorithm. In the “Results” we illustrate the procedure with two examples drawn from molecular biology. The first example is a circuit from synthetic biology, the pulse generator circuit and the second one is a signal transduction network from Cancer Systems Biology. Finally, in the “Discussion and Conclusions” we summarize the new procedure, we give some additional remarks, and we point out how these findings should be of immediate interest to researchers in computational biology applied to cancer drug development.
Methods
Problem formulation
where ψ(p) is the probability of parameter vector p, \(\mathbb {P}\) is the parameter space and \(\zeta ^{\mathcal {S}}_{\tau }(p)\) is the evaluation function of capability τ for the system \(\mathcal {S}\). Kitano’s definition is very general and can be used in a wide number of applications and problems.
where the state space vector x denotes the concentration of some substances relevant to the biological phenomena under consideration; p denotes the vector of system parameters taking values in the parameter space \(\mathbb {P}\), a subset of the positive orthant \(\mathbb {R}_{>0}^{q}\) of \(\mathbb {R}^{q}\); u denotes the input vector; i.e., the external stimuli acting on the system; and z denotes the output response of the system; e.g., the evaluation function of interest.
Conditional robustness
Conditional robustness: short motivation
In several applications, and in Systems Biology as a special case, it is of interest to search for a (possibly small) subset of the parameter vector having strong influence on the evaluation function \(\zeta ^{\mathcal {S}}_{\tau }(\cdot)\), and allowing the system to exhibit extreme values for it.
For example, in the case of translational oncology and targeted therapy [6–8], the idea is that one has to carefully choose a few nodes along the pathway where to act pharmacologically, with the aim of achieving the best possible benefit [9, 10, 37]. Also in Synthetic Biology, the fine-tuning through a few parameters of a synthetic network is searched for to produce the input/output behaviors characteristic that can be used as the data sheet of biological circuits [14].
To this purpose, we introduce the notion of conditional robustness as follows.
Density probabilities on P and Z
We can consider a given vector p in the parameter space as a realization of a vector random variable P on a measurable space \((\mathbb {P},\mathcal {A})\). We denote by F _{ P }(p) the cumulative distribution function of P, which is a model of the “a priori” knowledge on P. Let f _{ P }(p) denote the corresponding probability density function.
Similarly, the evaluation function \(\zeta ^{\mathcal {S}}_{\tau }(\cdot)\) can be interpreted as a transformation from the random variable P to the random variable Z, whose realizations are given by \(z=\zeta ^{\mathcal {S}}_{\tau }(p)\). Hence, let F _{ Z }(z) and f _{ Z }(z) denote the cumulative distribution function and the probability density function of the model output Z, respectively, that is, the transformation of the corresponding measures F _{ P }(p) and f _{ P }(p) through the mapping ζ.
Partitioning definition domain Z
Notice that the above definitions assume that the output function spans the non negative real axes. We will call the class of subsets of type (4a) the “lower set”, and the ones given by (4b) the “upper set”, and we will also use the notation L=L(α) and U=U(α), respectively, to make it clearer. In Fig. 3 the lower set for the proliferation indicator is marked with dark red color under the f _{ Z }(z) curve.
Of course, the three values a, b and α in (5) are not independent.
A different approach to the partitioning of the output space \(\mathcal {D}_{Z}\) is based on population quartile; i.e., by collecting a fixed number of samples, say the 10 % of samples, with the lower values of the evaluation function. Numerical experiments with such an approach, not reported here for brevity, yield results pretty much equal to the ones achieved with the probability approach proposed in (4).
Conditioning on P _{ i }
Suppose now that we choose to fix the i-th component of the parameter vector at a given value \(\hat {p_{i}}\). We introduce the notion of conditional robustness of system \(\mathcal {S}\) given that the scalar variable P _{ i } is fixed at \(\hat {p_{i}}\). It will be studied by means of the conditional probability on the output Z; i.e., the probability density function \(f_{Z|P_{i}}(z)\) of Z upon selection of P _{ i } to a given value \(\hat {p_{i}}\). An example of conditional density distribution is shown in Fig. 3 with a green line.
The problem: qualitative formulation
With the above notation, the problem studied in this paper is that of selecting a suitable conditioning set K (a proper subset of parameter vector and the corresponding values), achieving a conditional probability density function f _{ Z|K }(z) properly placed on the support interval shared with the density f _{ Z }(z) (Fig. 3). With “properly placed” we mean, as an example, that the conditional probability of having values of the output function in a given (upper or lower) set T has to be much larger than the probability of the same set in the unconditioned case. This allows for a moment independent formulation of our problem (again, see Fig. 3).
The problem name: fragility problem
Taking into account that our main interest originates in the oncological context, and more precisely in reducing the cell proliferation attitude, we will call such a problem the fragility problem.
Solution to the fragility problem
The proposed solution: conditional density \(f_{P_{i}|T}(p_{i})\)
The proposed approach to the solution of the fragility problem relies on the availability of the probabilities \(f_{P_{i}|T}(p_{i})\) of each parameter P _{ i }, interpreted as a random variable, conditioned by the values of the output function Z belonging to a given subset T of the output space. Assuming the definitions (4) for such a subset, we use the density \(f_{P_{i}|T}(p_{i})\) to characterize the portion of the parameter space along the direction P _{ i } giving rise to the values of the output function belonging to such set T.
The proposed solution: selection of set K
In order to select a conditioning set K, the problem is now to search for the directions in the parameter space giving a sufficiently sharp separation of the conditional probabilities \(f_{P_{i}|L}(p_{i})\) and \(f_{P_{i}|U}(p_{i})\); i.e., the set K which has a stronger influence on the system output Z; i.e., on the evaluation function of interest. The basic idea is to compute for each of the components of the parameter vector the two conditional densities, \(f_{P_{i}|U}(p_{i})\) and \(f_{P_{i}|L}(p_{i})\). The separation properties of each parameter on Z are evaluated through an index measuring the intersection of the corresponding conditional densities.
The components of the parameter vector yielding the smaller intersection are candidates to participate in the conditioning set K.
The proposed solution: key steps
Computational algorithm
- 1.
Sampling the parameter space Generate a sufficiently large set of samples P _{ S } of the parameter vector; i.e., a sufficiently large number of points in the parameter space \(\mathbb {P}\).
- 2.
In silico analysis For each point p∈P _{ S }, generate a sample \( z=\zeta ^{\mathcal {S}}_{\tau }{(}{p}{)}\) of the output function by means of numerical simulation of model \(\mathcal {S} \), thus building up the set of samples \(\mathcal {Z} := \left \{z : z = \zeta ^{\mathcal {S}}_{\tau }{(}{p}{)}, \, \forall p \in P_{S}\right \}\).
- 3.
Partitioning output sample set \(\boldsymbol{\mathcal {Z}}\) Compute the upper and lower sets U=U(α) and L=L(α), respectively, of set \(\mathcal {Z}\), for a given bound α.
- 4.
Estimating conditional densities For each component p _{ i } of parameter ∈p, estimate the two conditional densities \(f_{P_{i}|U}(p_{i})\) and \(f_{P_{i}|L}(p_{i})\).
- 5.
Parameter selection For each parameter component p _{ i }∈p, compute the moment independent robustness index (defined in the following equation in (7) between the densities \(f_{P_{i}|U}(p_{i})\) and \(f_{P_{i}|L}(p_{i})\), and select the three parameters having the lowest index values.
- 6.
Conditioning set selection Set K equal to the three parameters selected above and choose, for each of them, the value corresponding to the maximum of the associated conditional density \(f_{P_{i}|L}(p_{i})\) (if one seeks to reduce the evaluation function).
The proposed solution: step by step discussion
The basic steps outlined above involve a number of ingredients and implications, described in the following.
Sampling the parameter space
In this paper, we assume Latin Hypercube Sampling with Linearly spaced samples (L ^{2}HS). The use of Latin Hypercube Sampling is common in global approaches to the analysis of complex biochemical networks such as the ones based on Global Sensitivity Analysis. For a comparison of different sampling approaches, see [21] and the references therein.
The issue of linear vs logarithmic sampling is not covered here for the sake of brevity (see [21, 34] for the reasons behind logarithmic sampling).
As for the number of samples, here it has been chosen by looking at the variation of the density f _{ Z }(z) while increasing the sample number, and stopping once a steady density has been reached.
In silico analysis
The computation of the evaluation function, of the model output z, is carried out by means of numerical integration of the ODE model \(\mathcal {S}\), for each point of the sample set P _{ S }, using tools such as Matlab, Octave or similar ones.
Partitioning the output sample set \(\mathcal {Z}\) and the parameter space P _{ S }
Estimating conditional densities
To estimate the conditional densities \(f_{P_{i}|U}(p_{i})\) and \(f_{P_{i}|L}(p_{i})\), we use a kernel density approach. Each density is estimated as a superposition of a suitable number of gaussian densities. The non-parametric nature of the approach is of special interest here, as well as its multivariate nature [38–40].
For each component p _{ i } of the parameter vector, we pick up conditioned values from the sets P _{ U } and P _{ L } (i.e., we pick up the values giving rise to output values in the two sets U and L) and use them to estimate the corresponding conditional densities.
Criteria to evaluate the accuracy of the estimated density turns out to be important whenever the raw data do not show a sufficiently large set of values.
Conditioning set selection
The MIRI indicator resembles the moment independent sensitivity indicator δ proposed by [41]. While the δ indicator is aimed at studying the shift between an unconditional case and a “sensitivity” one, here we compare two “conditioned” cases. Also, here the conditioning events, the lower and upper sets, have fixed probabilities α, while the δ indicator comprises an expectation over the conditioning variables.
The MIRI is then used to sort the parameter components: the selected ones are those with larger shifts between the two conditional densities, those with higher MIRI values. The conditioning values associated with those selected parameters are the values corresponding to the mode (the maximum of the associate conditional density).
Results
To verify the performance of the proposed conditional robustness approach, we applied it in two biologically relevant cases. The first example is taken from the field of Synthetic Biology: the incoherent feedforward network called “pulse generator network” [17]. We used this simple example to prove the general application of our approach and also to work with a tractable dimension of the parameter space.
The second example is an application from Cancer Systems Biology [3]: we investigate the robustness of the EGFR IGF1R pathway in lung cancer [20]. This model is directly related to the inspiration of the methodology we proposed in this paper.
Pulse generator network
where k _{1}=5 nM/min, k _{12}=20 nM/min, λ _{2}=0.01 nM/min, λ=0.04 nM/min, K _{1}=1 nM, K _{2}=100 nM, and n _{1}=n _{2}=3 (see [43]).
We applied our approach taking three different evaluation functions for the output function Y: the area under the curve, the maximum value and the time of occurrence of the maximum.
We set up in silico simulations for the model (8) to generate the measure of the area under the curve of signal Y as described in the “Methods” section: we selected the parameter space \(\mathbb {P}\) by using as the lower bounds \(\frac {1}{10}\) and as the upper ones 10 times the nominal value in \(\mathbb {R}_{>0}^{6}\). We fixed the number of sample to N=10000 and we generated 100 realizations of the subset P _{ S,i },i=1,…,100, of the parameters space with L ^{2} H S. The kernel distribution estimates for the 100 realizations of the three evaluation functions are shown in Fig. 6 b, 6 c and 6 d, respectively.
For the area of the Y signal we applied the proposed algorithm with the goal of selecting the parameters that are most significant to maximize such an indicator. According to the procedure we fixed the α parameter for the sets L and U, as in (4), to α=0.1. Then, the kernel distributions were estimated. Additional file 1: Figure S1 shows the upper set distributions \(f_{P_{i}|U}(p_{i})\) (red line) and the lower distributions \(f_{P_{i}|L}(p_{i})\) (blue line) for the whole parameter set, i.e., k _{1}, K _{1}, k _{12}, K _{2}, λ _{2} and λ, for a single realization of P _{ S }.
Since we are interested in maximizing the pulse area, we used the parameter densities conditioned on the upper sets, and chose as conditioning values the corresponding modes. Additional file 2: Figure S2 shows the λ, k _{12} and K _{2} lower and upper sets probability density for all the realizations and it confirms the low variability of the probability density shapes over the parameter space sampling.
Parameters setting for the pulse generator and EGFR-IGF1R in silico experiments
Experiment | Parameters setting |
---|---|
Y area | k _{12}=174.1101, K _{2}=789.4512, λ=0.0163 |
Y maximum | k _{12}=177.8501, λ=0.0204 |
Y time to maximum | k _{1}=43.4827, K _{2}=40.0512,λ=0.3360 |
Multiobjective of Y | k _{12}=177.8501, K _{2}=40.0512, λ=0.0163 |
E R K ^{∗} activity | p _{23}=73.9765, p _{24}=3.9577×10^{6}, p _{35}=23.5795, p _{36}=6.0595×10^{5}, p _{27}=1.0195 |
To further test our approach, we performed the robustness analysis using as evaluation function the maximum value reached by the signal Y. As shown in Additional file 3: Figure S3, MIRI values for a single realization of the parameters L ^{2} H S sampling (Additional file 3: Figure S3A) and MIRI box plot for 100 realization (Additional file 3: Figure S3B) suggested that two parameters are relevant for this evaluation function: k _{12} and λ.
We fixed k _{12} as in Table 1; i.e., the value giving rise to the maximum of the densities conditioned over the upper set and then we performed conditional robustness analysis. The simulation results are shown in Additional file 4: Figure S4A: the blue curve is the conditional density \(f_{Z|k_{12}}(z)\): its first moment is higher than the unconditioned distribution f _{ Z }(z); the yellow curve is the conditional density function, f _{ Z|λ }(z) (see Table 1) and the distribution spread at higher values of the evaluation function; the red line is the f _{ Z|k12,λ }(z) and the most probable output value is the one corresponding to the density maximum. An example of the conditioned behavior is showed in Additional file 4: Figure S4B.
For the third evaluation function, the time to maximum of the signal Y, in silico simulations were performed and MIRI for single realization of the parameter space and 100 realizations are shown in Additional file 5: Figure S5. The key parameters are k _{1}, K _{2} and λ.
The parameter values to be used for the conditioning operation are in Table 1 from the lower set parameters distributions. In this case, we are interested in minimizing the time to maximum of the signal Y; the results for the analysis are presented in Additional file 6: Figure S6A. The simulation example with the selected values of k _{1}, K _{2} and λ confirmed the ability of the proposed procedure to give a ready response of the signal Y, that upon conditioning reaches the maximum in 1.5 min (Additional file 6: Figure S6B).
EGF-IGF network in lung cancer
The biochemical network is sketched in Fig. 10 b. The activation of the EGFR and I G F1R (x _{1} and x _{2}, respectively) induced by the ligands binding triggers a cascade of complex proteins interactions that leads to the activation of the ERK protein (x _{7}). There is biological evidence that the time behavior of ERK concentration in its active form is strictly related to the proliferation attitude of the cell [11]. The complete model can be found in [20] and [46]. For this system the identification problem has been studied in [47].
Following an approach widely used in literature, we assume that each element p _{ i } of the parameter vector p, \(p \in \mathbb {R}^{39}\), can assume values in a range defined by two multiplicative bounds w.r.t. a nominal (wild-type) value p _{ w t,i }: a lower bound b _{ ℓ,i }=c _{ ℓ,i } p _{ w t,i } and an upper bound b _{ u,i }=c _{ u,i } p _{ w t,i }, where the coefficients c _{ ℓ,i } and c _{ u,i } are the multiplicative perturbations, c _{ ℓ,i }<1, ∀i and c _{ u,i }<1, ∀i. Hence, the parameter space \(\mathbb {P}\) is given by the cartesian product \(\mathbb {P} =\prod _{i=1}^{39}[c_{\ell,i} p_{wt,i}, \, c_{u,i} p_{wt,i}]\).
The key parameters characterizing the ODE model of a biochemical network are activation rates and Michaelis-Menten constants. In Table Additional file 7: Table S1, the Michaelis-Menten constants are all those whose names start with KM. While sampling the parameter space, we will consider the whole parameter vector with 39 entries, since the relationship between activation rates and Michaelis-Menten constants allows us to check the consistency of the analysis. We also assume that the total mean values of the number of molecules in the cell does not change significantly in NSCLC when compared to normal cells (Additional file 8: Table S2 Table).
Proliferation output pdf measures
Mean | Variance | Mode | |
---|---|---|---|
f _{ Z }(z) | 2.6×10^{8} | 1.5805×10^{17} | 2.1×10^{7} |
\(f_{Z|p_{23}}(z)\) | 1.5824×10^{8} | 7.6840×10^{16} | 1.0×10^{7} |
\(f_{Z|p_{35}}(z)\) | 1.6075×10^{8} | 7.0478×10^{16} | 9.6×10^{6} |
\(f_{Z|p_{27}}(z)\) | 1.6380×10^{8} | 9.5956×10^{16} | 7.2×10^{6} |
\(f_{Z|p_{24}}(z)\) | 1.3436×10^{8} | 6.6212×10^{16} | 7.085×10^{6} |
\(f_{Z|p_{36}}(z)\) | 1.4505×10^{8} | 7.0479×10^{16} | 7.41×10^{6} |
\(f_{Z|p_{23},p_{35},p_{27}}(z)\) | 4.2358×10^{7} | 8.0128×10^{15} | 2.528×10^{6} |
\(f_{Z|p_{23},p_{24},p_{35},p_{36},p_{27}}(z)\) | 8.2938×10^{6} | 5.3817×10^{14} | 5.131×10^{5} |
f _{ Z }(z) @ \({x_{1}^{0}}\) and \({x_{2}^{0}}\) perturbed | 3.2996×10^{8} | 2.0223×10^{17} | 2.451×10^{7} |
\(f_{Z|p_{23},p_{24},p_{35},p_{36}p_{27}}(z)\) @ \({x_{1}^{0}}\) and \({x_{2}^{0}}\) perturbed | 1.3099×10^{7} | 1.4868e×10^{15} | 2.528×10^{6} |
Discussion and Conclusions
We have presented a novel approach to study the robustness of complex biological networks and discover their fragility. We have used this approach as a base to design a new algorithm aimed at selecting a few of the system parameters that would allow us to achieve desired conditional robustness properties. The approach is based on the combination of robustness analysis, moment independent robustness indicator, and estimated kernel conditional densities. The approach takes its motivation from Cancer Systems Biology and the associated problem of drug therapy strategies, namely the selection of crucial proteins that need to be inhibited with the administration of therapeutic compounds with the goal to block cancer proliferation and progression [3]. Nevertheless, the proposed algorithm is general in nature and its suitability for a larger class of problems has been shown by means of a case study from Synthetic Biology.
Our algorithm is similar to Global Sensitivity Analysis (GSA) approaches as far as the sampling of the parameter space is concerned, while it is different from the point of view of the goals of analysis. Both GSA and our Conditional Robustness Algorithm (CRA) have to deal with the sampling of a large parameter space. Key approaches to this issue are Sobol Low Discrepancy Sequence, Monte Carlo techniques, Latin Hypercube and variants of all of the above. Here, we use Latin Hypercube sampling; additional results (data not shown) show that Sobol techniques yield similar results. Sampling techniques are discussed in [21, 23, 34] and the references therein. Uniform/linear vs logarithmic distribution is another relevant issue when sampling of large spaces has to be performed. The “priori” knowledge on parameter variability is related to to this issue. See [21, 34] for a more detailed discussion. In this paper, we focus on Latin Hypercube Sampling with Linearly spaced samples (L ^{2} H S) since our focus is on the key properties of the proposed algorithm.
As for the main goal of our approach, we are seeking a strategy to identify regions of the parameter space yielding desired system behavior. On the contrary, most of the GSA results are finalized to investigate the sensitivity of key performance indices with respect to variation on system parameters, usually in the classical form of derivative of the index with respect to parameters. Hence, GSA methods are quite useful as analysis tools and they provide “results that could potentially help in the optimization of the system” [27]. Similarly, they do not yield specific design criteria to enforce given system capabilities, as with this problem we are investigating on the problem we are investigating. Here our main interest is on the selection of some parameters and their values to achieve a desired behavior, and not on the deduction of the parameters having the largest influence on a given performance function.
We have shown how the CRA tool can be used both in Cancer Systems Biology and in Synthetic Biology through two examples, the EGFR- I G F1R network in non small cell lung cancer and a reduced model of pulse generator network.
The development of computational tools to study cancer disease is an emergent field in Cancer Systems Biology. The general formulation of robustness proposed in [36] and its extensions to drug development in cancer research [9] were explored in our paper by introducing the idea of conditional robustness.
Our approach assumes the knowledge of a mathematical model of the biological system and the associated parameter space. In addition, a proper evaluation function whose behavior over the whole parameter space can be described by means of an associated density function must be defined. The problem we are interested in, namely the Conditional Robustness Problem, is that of shifting such a density to a desired region in order to achieve the suitable global behavior.
The mathematical model used here is based on ordinary differential equations; nevertheless, several others types of mathematical tools, e.g. stochastic models, boolean models and many other, can be used as well with minor modification. The dependence of the Conditional Robustness Algorithm on the specific modeling approach is confined to the two blocks in dark blu in Fig. 4. More specifically, the key condition is the availability of a computable relationship between system parameters and evaluation function. Notice that both numerical and analytical maps can be used. From a mathematical point of view, any modeling approach allowing to compute the input-output map in (2) can be directly cast in our general framework. As an example, the Boolean model in [51] has been implemented by the authors through a proper simulation model, which depends on a number of rules, initial states and other parameters. Each execution of the block “in silico generation of ζ” in our CRA, for a given parameter configuration, turns out to be a corresponding simulation of the above boolean model.
In silico experiments are used to generate the density function describing the property of interest as well as to estimate the effect on such density of the parameter conditioning, measuring the amount of overlapping between the parameters conditioned density functions. From a technical point of view, the distance between the parameter densities is a separation measurement in the sense of [52], and it is evaluated trough the Moment Independent Robustness Index (MIRI) which is inspired to the measure introduced in [41]. Borgonovo et al. illustrated that variance is not always sufficient enough to describe uncertainty and pointed out that a sensitivity measure should refer to the entire output distribution instead of a particular moment. Therefore we introduce in our algorithm the MIRI. A higher value of the MIRI indicates a greater separation between the densities and a greater capacity of the parameter to discriminate between different behaviors of the evaluation function. The parameters with largest MIRI are selected as the conditioning ones, and their values are fixed to the modes of the conditional densities, allowing to achieve the desired system behavior.
The analysis of the EGFR-IGF1R model in non small cell lung cancer is an example of nontrivial biochemical network in cancer application. The output of the EGFR-IGF1R network is the active form of the ERK protein, which is related to the cancer proliferation. We considered the area under the active ERK curve as a measure of proliferation, and our goal is to reduce this network output almost to zero, namely to silence the proliferation activity.
In [11] the authors concluded that cell proliferation is effectively silenced only when active ERK protein level falls below a lower threshold. Their findings provide rationale for combined inhibition of multiple nodes in the ERK pathway, since acting on a single node turns out not to be enough to constrain ERK signal below the threshold required for a proper proliferation level. A multiple node action is the objective of the computational algorithm we are proposing, whose goal is finding to how to silence ERK protein with a conditioning action on multiple model parameters.
The CRA results suggest three points in the network that could be potentials targets for the inhibition of cancer cell proliferation: the activation rate of RAF/Ras in the MAPK cascade, and the activation rate and Michaelis Menten constants of PP2A both with MEK and ERK. In silico experiments show that ERK silencing is more effective whenever all the three nodes are conditioned (by fixing the corresponding five parameters, see Table 1), i.e., they are candidate that could be targeted.
The therapeutic strategies recently proposed, primarily pertaining to co-targeting EGFR and IGF1R receptors, exhibited notable advantages in overcoming resistance to standard chemotherapy. Furthermore, these techniques might offer benefits beyond the limited effects of single-agent targeting previously reported. The strategies of blocking these pathways in combination with the result obtained from our conditional robustness, suggest a novel potential approach to develop future and more effective therapies for cancer treatments [53].
It is well known that the over expression of both EGFR and IGR1R plays a crucial role in the proliferation of cancer cells [44]. The increase in the mean value and variance of proliferation index associated with an increase on the number of receptors showed in this study confirm their role in cancer (see Table 2 and Fig. 13 b).
The results obtained can be applied to other problems concerned with drug development. For example, if we have a set of drug candidates that target proteins in a pathway, our methods can be used to asses which is the best compound to use as a single agent or in combination. Figure 12 b can be interpreted as the best combination of target to be used and Fig. 13 b predicts the ability of the target to silence ERK activity.
We considered the pulse generator network from Synthetic Biology to study the validity of CRA in different frameworks. This analysis allows us to conclude that the algorithm is able to select the proper parameters in order to maximize the pulse transferring (and intensity of the pulse) and minimize the time to reach the pulse maximum. The example also demonstrates the solution for a multi objective evaluation function. For this type of synthetic networks, the knowledge of the most significant parameters and the corresponding values allowing achievement of specific input/output behaviors is of crucial relevance to characterize the biological circuits and is a key piece of information contained in their data sheet.
The algorithm proposed in this paper has four classes of configuration data: the multiplicative perturbation constants c _{ ℓ,i } and c _{ u,i } defining the parameter space; the number N of samples to be taken in the parameter space and used to run in silico experiments to compute the corresponding set of samples \(\mathcal {Z}\) of the evaluation function; the probability α defining the subsets of \(\mathcal {Z}\) used to carry out conditional robustness; and the number of parameters to be fixed by the conditioning operation. The selection of those inputs is a user choice and depends on the problem at hand. Below we give some guidelines to select proper values.
In the framework of mathematical modeling of biological systems, and in particular Computational Biology, the nominal system parameters can be either estimated from experimental data with a given range of uncertainty, (and quite often the experiments are carried out in different settings and conditions, or derived using biochemical/biological a priori knowledge of the system, and therefore only the order of magnitude can be set). In this setting, a strategy widely used in the literature is to assume that the parameters space spans for two orders of magnitude below and above the nominal value [21, 27, 34].
The choice of the number N of samples in the parameters space and the probability α defining the extremal subsets U and L are somehow interdependent. The number N of samples has to be set in order to guarantee that the estimated output density function f _{ Z }(z), as a function of increasing N, has reached a steady state shape [34]. The subsets L and U build upon this distribution generate data used to estimate the conditional densities \(f_{P_{i}|L}(p)\) and \(f_{P_{i}|U}(p)\) of each parameter. To generate a satisfactory estimation of these densities at least 1.000 samples are required. Hence, on the average, a good estimate for a lower bound on the number N of overall samples is equal to \(\frac {1.000}{\alpha }\). As an example, in the EGFR-IGFR model we fixed α=0.1 and we generated N=10.000 samples in the parameter space. Notice that, in this case, such a value is quite higher than the number of samples required to achieve a satisfactory (i.e., stationary with respect to N) density f _{ Z }(z); hence we are over-sampling the parameter space.
The fourth configuration constant of the CRA algorithm is the number of parameters to be fixed, i.e., to be conditioned. Such a number can be 6set by trial and error, looking both to the number of densities \(f_{P_{i}|L}(p)\) and \(f_{P_{i}|U}(p)\) that turn out to be sufficiently different and to the overall effect on the conditional density f _{ Z|K }(z).
The computational cost of our CRA has two main sources: the L ^{2} H S sampling of the parameter space and the ODE numerical integration of the whole set of samples. As for the former, in [54] the computational cost of Matlab function lhsdesign, is compared with other optimization techniques with respect to the number of parameters and sampled points. As for the latter, the computational cost associated with the numerical integration of a single ODE depends heavily on model dimensions and model properties, such as its stiffness. Parallel implementations, such as CUDA models, can be used to deal with this issue.
We implemented a parallel versions of our algorithm: referring to the diagram in Fig. 4, the most significant computation cost is the in silico data generation associated with the block “For each sample”, which integrates the ODE model. Each ODE integration is independent of the others and therefore they can be executed in parallel on a cluster of work-stations: on a 7 core processor the computation time to generate the set of samples \(\mathcal {Z}\) can be reduced by a factor of 7.
In conclusion, our Conditional Robustness Algorithm is a new method to study the role of key parameters to discover the system robustness/fragility or in conditioning the systems to the wished behaviour. The CRA significantly contributes to formal methods in computational systems biology and introduces a new framework to identify target identification and to support drug discovery in oncology.
Notes
Declarations
Acknowledgements
This work was supported by an AIRC grant (15713/2014) to LC, LV and BF. We thank Mustafa Khammash and Gabriele Lillacci for helpful comments and discussion.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
- Karamouzis MV, Papavassiliou AG. Tackling the cancer signal transduction “labyrinth”: A combinatorial use of biochemical tools with mathematical models will enhance the identification of optimal targets for each molecular defect. Cancer. 2014; 120(3):316–22. doi:http://dx.doi.org/10.1002/cncr.28424.View ArticlePubMedGoogle Scholar
- Tabchy A, Ma CX, Bose R, Ellis MJ. Incorporating genomics into breast cancer clinical trials and care. Clin Cancer Res. 2013; 19(23):6371–379. doi:http://dx.doi.org/10.1158/1078-0432.CCR-13-0837.View ArticlePubMedGoogle Scholar
- Werner HMJ, Mills GB, Ram PT. Cancer systems biology: a peek into the future of patient care?Nat Rev Clin Oncol. 2014; 11:67–176.View ArticleGoogle Scholar
- Gonzalez-Angulo AM, Hennessy BTJ, Mills GB. Future of personalized medicine in oncology: A systems biology approach. J Clin Oncol Off J Am Soc Clin Oncol. 2010; 28(16):2777–783. doi:http://dx.doi.org/10.1200/JCO.2009.27.0777.View ArticleGoogle Scholar
- Wang E. Cancer systems biology. CRC Press. 2010.Google Scholar
- Janku F, Wheler JJ, Westin SN, Moulder SL, Naing A, Tsimberidou AM, et al. Pi3k/akt/mtor inhibitors in patients with breast and gynecologic malignancies harboring pik3ca mutations. J Clin Oncol Off J Am Soc Clin Oncol. 2012; 30(8):777–82. doi:http://dx.doi.org/10.1200/JCO.2011.36.1196.View ArticleGoogle Scholar
- Cancer Genome Atlas Research Network. Comprehensive molecular characterization of human colon and rectal cancer. Nature. 2012; 487(7407):330–7. doi:http://dx.doi.org/10.1038/nature11252.View ArticleGoogle Scholar
- Cancer Genome Atlas Research Network. Comprehensive molecular characterization of clear cell renal cell carcinoma. Nature. 2013; 499(7456):43–9. doi:http://dx.doi.org/10.1038/nature12222.View ArticleGoogle Scholar
- Kitano H. A robustness-based approach to systems-oriented drug design. Nat Rev Drug Discov. 2007; 6:202–10.View ArticlePubMedGoogle Scholar
- Sikkema AH, den Dunnen WFA, Diks SH, Peppelenbosch MP, de Bont ESJM. Optimizing targeted cancer therapy: Towards clinical application of systems biology approaches. Crit Rev Oncol Hematol. 2012; 82(2):171–86. doi:http://dx.doi.org/10.1016/j.critrevonc.2011.05.002.View ArticlePubMedGoogle Scholar
- Albeck JG, Mills GB, Brugge JS. Frequency-modulated pulses of erk activity transmit quantitative proliferation signals. Mol Cell. 2013; 49(2):249–61. doi:http://dx.doi.org/10.1016/j.molcel.2012.11.002.View ArticlePubMedPubMed CentralGoogle Scholar
- Basu A, Bodycombe NE, Cheah JH, Price EV, Liu K, Schaefer GI, et al. An interactive resource to identify cancer genetic and lineage dependencies targeted by small molecules. Cell. 2013; 154(5):1151–1161. doi:http://dx.doi.org/10.1016/j.cell.2013.08.003.View ArticlePubMedPubMed CentralGoogle Scholar
- Andrianantoandro E, Basu S, Karig DK, Weiss R. Synthetic biology: new engineering rules for an emerging discipline. Mol Syst Biol. 2006; 2:2006–0028.View ArticlePubMedPubMed CentralGoogle Scholar
- Khalil AS, Collins JJ. Synthetic biology: applications come of age. Nat Rev Genet. 2010; 11(5):367–79. doi:http://dx.doi.org/10.1038/nrg2775.View ArticlePubMedPubMed CentralGoogle Scholar
- Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in escherichia coli. Nature. 2000; 403(6767):339–42. doi:http://dx.doi.org/10.1038/35002131.View ArticlePubMedGoogle Scholar
- Elowitz M, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000; 403:335–8.View ArticlePubMedGoogle Scholar
- Basu S, Mehreja R, Thiberge S, Chen MT, Weiss R. Spatiotemporal control of gene expression with pulse-generating networks. Proc Natl Acad Sci USA. 2004; 101(17):6355–0. doi:http://dx.doi.org/10.1073/pnas.0307571101.View ArticlePubMedPubMed CentralGoogle Scholar
- Lillacci G, Khammash M. Parameter estimation and model selection in computational biology. PLoS Comput Biol. 2010; 6(3):1000696. doi:http://dx.doi.org/10.1371/journal.pcbi.1000696.View ArticleGoogle Scholar
- Oguz C, Laomettachit T, Chen K, Watson L, Baumann W, Tyson J. Optimization and model reduction in the high dimensional parameter space of a budding yeast cell cycle model. BMC Syst Biol. 2013; 7(1):53. doi:http://dx.doi.org/10.1186/1752-0509-7-53.View ArticlePubMedPubMed CentralGoogle Scholar
- Bianconi F, Baldelli E, Ludovini V, Crinò L, Valigi P. Computational model of EGFR and IGF1R pathways in lung cancer: A systems biology approach for translational oncology. Biotechnol Adv. 2012; 30(1):142–53. doi:http://dx.doi.org/10.1016/j.biotechadv.2011.05.010.View ArticlePubMedGoogle Scholar
- Rodriguez-Fernandez M, Banga JR, Doyle III FJ. Novel global sensitivity analysis methodology accounting for the crucial role of the distribution of input parameters: application to systems biology models. Int J Robust Nonlinear Control. 2012; 22(10):1082–1102. doi:http://dx.doi.org/10.1002/rnc.2797.View ArticleGoogle Scholar
- Tarantola S, Becker W, Zeitz D. A comparison of two sampling methods for global sensitivity analysis. Comput Phys Commun. 2012; 183(5):1061–1072. doi:http://dx.doi.org/10.1016/j.cpc.2011.12.015.View ArticleGoogle Scholar
- Lebedeva G, Sorokin A, Faratian D, Mullen P, Goltsov A, Langdon SP, et al. Model-based global sensitivity analysis as applied to identification of anti-cancer drug targets and biomarkers of drug resistance in the erbb2/3 network. Eur J Pharm Sci. 2012; 46(4):244–58. doi:http://dx.doi.org/10.1016/j.ejps.2011.10.026.View ArticlePubMedPubMed CentralGoogle Scholar
- Von Dassow G, Meir E, Munro EM, Odell GM. The segment polarity network is a robust developmental module. Nature. 2000; 406(6792):188–92.View ArticlePubMedGoogle Scholar
- Dayarian A, Chaves M, Sontag ED, Sengupta AM. Shape, size, and robustness: Feasible regions in the parameter space of biochemical networks. PLoS Comput Biol. 2009; 5(1):1000256. doi:http://dx.doi.org/10.1371/journal.pcbi.1000256.t002.View ArticleGoogle Scholar
- Chaves M, Sengupta A, Sontag E. Geometry and topology of parameter space: investigating measures of robustness in regulatory networks. J Math Biol. 2009; 59(3):315–58. doi:http://dx.doi.org/10.1007/s00285-008-0230-y.Geometry.View ArticlePubMedPubMed CentralGoogle Scholar
- Rizk A, Batt G, Fages F, Soliman S. A general computational method for robustness analysis with applications to synthetic gene networks. Bioinformatics. 2009; 25(12):169–78. doi:http://dx.doi.org/10.1093/bioinformatics/btp200.View ArticleGoogle Scholar
- Kwon YK, Cho KH. Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics. Bioinformatics (Oxford, England). 2008; 24(7):987–4. doi:http://dx.doi.org/10.1093/bioinformatics/btn060.View ArticleGoogle Scholar
- Morohashi M, Winn AE, Borisuk MT, Bolouri H, Doyle J, Kitano H. Robustness as a measure of plausibility in models of biochemical networks. J Theor Biol. 2002; 216(1):19–30. doi:http://dx.doi.org/10.1006/jtbi.2002.2537.View ArticlePubMedGoogle Scholar
- Hafner M, Koeppl H, Hasler M, Wagner A. ‘glocal’ robustness analysis and model discrimination for circadian oscillators. PLoS Comput Biol. 2009; 5(10):1000534. doi:http://dx.doi.org/10.1371/journal.pcbi.1000534.View ArticleGoogle Scholar
- Ceska M, Safránek D, Dražan S, Brim L. Robustness analysis of stochastic biochemical systems. PloS one. 2014; 9(4):94553. doi:http://dx.doi.org/10.1371/journal.pone.0094553.View ArticleGoogle Scholar
- Salerno L, Cosentino C, Merola A, Bates DG, Amato F. Validation of a model of the GAL regulatory system via robustness analysis of its bistability characteristics. BMC Syst Biol. 2013; 7:39. doi:http://dx.doi.org/10.1186/1752-0509-7-39.View ArticlePubMedPubMed CentralGoogle Scholar
- Cedersund G. Conclusions via unique predictions obtained despite unidentifiability: new definitions and a general method. FEBS J. 2012; 279(18):3513–527. doi:http://dx.doi.org/10.1111/j.1742-4658.2012.08725.x.View ArticlePubMedGoogle Scholar
- Bianconi F, Baldelli E, Ludovini V, Crinò L, Perruccio K, Valigi P. Robustness of complex feedback systems: application to oncological biochemical network. Int J Confl Manag. 2013; 86(7):1304–1321.Google Scholar
- Mathew S, Bartels J, Banerjee I, Vodovotz Y. Global sensitivity analysis of a mathematical model of acute inflammation identifies nonlinear dependence of cumulative tissue damage on host interleukin-6 responses. J Theor Biol. 2014; 358:132–48. doi:http://dx.doi.org/10.1016/j.jtbi.2014.05.036.View ArticlePubMedPubMed CentralGoogle Scholar
- Kitano H. Towards a theory of biological robustness. Mol Syst Biol. 2007; 3(137):137. doi:http://dx.doi.org/10.1038/msb4100179.PubMedPubMed CentralGoogle Scholar
- Huang CY, Huang CH, Chang PM-H, Wu MY, Ng KL. In silico identification of potential targets and drugs for non-small cell lung cancer. IET Syst Biol. 2014; 8(2):56–66.View ArticlePubMedGoogle Scholar
- Sheather SJ. Density estimation. Stat Sci. 2004; 19(4):588–97.View ArticleGoogle Scholar
- Borgonovo E, Tarantola S, Plischke E, Morris MD. Transformations and invariance in the sensitivity analysis of computer experiments. J R Stat Soc Ser B Stat Methodol. 2013; 76:925–47.View ArticleGoogle Scholar
- Hwang JN, Lay SR, Lippman A. Nonparametric multivariate density estimation: a comparative study. IEEE Trans Signal Process. 1994; 42(10):2795–810.View ArticleGoogle Scholar
- Borgonovo E. A new uncertainty importance measure. Reliab Eng Syst Saf. 2007; 92(6):771–84.View ArticleGoogle Scholar
- Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007; 8(6):450–61. doi:http://dx.doi.org/10.1038/nrg2102.View ArticlePubMedGoogle Scholar
- Szallasi Z, Stelling J, Periwal V. System modeling in cellular biology. from concepts to nuts and bolts. 2006.Google Scholar
- Ludovini V, Bellezza G, Pistola L, Bianconi F, Carlo LD, Sidoni A, et al. High coexpression of both insulin-like growth factor receptor-1 (igfr-1) and epidermal growth factor receptor (egfr) is associated with shorter disease-free survival in resected non-small-cell lung cancer patients. Ann Oncol. 2009 May. 2009; 20(5):842–9.Google Scholar
- Falconi A, Lopes G, Parker JL. Biomarkers and receptor targeted therapies reduce clinical trial risk in non–small-cell lung cancer. J Thorac Oncol. 2014; 9(2):163–9.View ArticlePubMedGoogle Scholar
- Bianconi F, Chelliah V. BIOMD0000000427 - Bianconi2012 - EGFR and IGF1R Pathway in Lung Cancer. http://www.ebi.ac.uk/biomodels-main/BIOMD0000000427.
- Bianconi F, Lillacci G, Valigi P. Dynamic modeling and parameter identification for biological networks: application to the dna damage and repair processes. 2010:478–510.Google Scholar
- Bianconi F, Baldelli E, Valigi P. An approach to the conditional robustness problem for biochemical networks. In: American Control Conference (ACC), 2014: 2014. p. 3417–424. doi:http://dx.doi.org/10.1109/ACC.2014.6859085.
- Ahmed S, Grant KG, Edwards LE, Rahman A, Cirit M, Goshe MB, et al. Data-driven modeling reconciles kinetics of erk phosphorylation, localization, and activity states.Mol Syst Biol. 2014; 10:718. doi:http://dx.doi.org/10.1002/msb.134708.View ArticlePubMedPubMed CentralGoogle Scholar
- Wu DW, Wu TC, Wu JY, Cheng YW, Chen YC, Lee MC, et al. Phosphorylation of paxillin confers cisplatin resistance in non-small cell lung cancer via activating ERK-mediated Bcl-2 expression. Oncogene. 2014; 33(35):4385–95. doi:http://dx.doi.org/10.1038/onc.2013.389.View ArticlePubMedGoogle Scholar
- Sahin O, Fröhlich H, Löbke C, Korf U, Burmester S, Majety M, et al. Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance. BMC Syst Biol. 2009; 3:1. doi:http://dx.doi.org/10.1186/1752-0509-3-1.View ArticlePubMedPubMed CentralGoogle Scholar
- Glick N. Measurements of separation among probability densities or random variables. Can J Stat. 1975; 3(2):267–76.View ArticleGoogle Scholar
- Vasan N, Boyer J, Herbst RS. A ras renaissance: Emerging targeted therapies for kras-mutated non-small cell lung cancer.Clin Cancer Res. 2014; 20(15):3921–930. doi:http://dx.doi.org/10.1158/1078-0432.CCR-13-1762.View ArticlePubMedGoogle Scholar
- Viana FA, Venter G, Balabanov V. An algorithm for fast optimal latin hypercube design of experiments. Int J Numer Methods Eng. 2010; 82(2):135–56.Google Scholar