Conditional robustness analysis for fragility discovery and target identification in biochemical networks and in cancer systems biology 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. Electronic supplementary material The online version of this article (doi:10.1186/s12918-015-0216-5) contains supplementary material, which is available to authorized users. 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. Figure 1 shows the context of interest of this paper. Signal transduction pathways are complex networks based on protein interactions that determine the propagation of extracellular inputs through the cytoplasm driving the timing of cellular survival, apoptosis and proliferation. Fig. 1 Cancer cell proliferation. Green cells are normal cells and red cells are tumor cells. The proliferation activity of normal and tumor cells can be measured looking at the activation of a proliferation protein, which is driven by a complex network based on protein interactions. In a population of cells the proliferation activity can be described by means of probability density for the proliferation protein; e.g., phosphorylated form of the extracellular signal-regulated kinase (ERK). The plots at the bottom show an example of probability density of a proliferation indicator in tumor (red line) and normal cells (green line), respectively 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. Figure 2 illustrates the fragility problem: we propose a methodology to shift the probability density function of a proliferation signal in cancer cells toward the density describing the normal cells. The solution proposed here stems from robustness analysis tools, extending our previous results in [34]. Several authors have introduced the visionary idea of applying robustness in drug development [9, 10, 23, 35], although the application of the general definition is still a key problem in Systems Biology. Fig. 2 Problem formulation. The fragility problem in the oncology context is related to the cancer cells signaling network, and the goal is to reduce the cell proliferation attitude acting in few target molecules. The problem is related to the conditional robustness problem, namely the problem of shifting the probability density function of a proliferation signal in cancer cells toward the density describing the normal cells 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 In this paper, we start from the general definition of robustness proposed by Kitano [36]: robustness is a property that allows a system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {S}$\end{document}S to maintain its property/capability τ against internal and external perturbations p. Kitano [36] proposes the following measure: (1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ R^{\mathcal{S}}_{\tau,\mathbb{P}} := \int_{\mathbb{P}} \psi(p) \zeta^{\mathcal{S}}_{\tau}(p) dp, $$ \end{document}Rτ,ℙS:=∫ℙψ(p)ζτS(p)dp, where ψ(p) is the probability of parameter vector p, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {P}$\end{document}ℙ is the parameter space and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\zeta ^{\mathcal {S}}_{\tau }(p)$\end{document}ζτS(p) is the evaluation function of capability τ for the system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {S}$\end{document}S. Kitano’s definition is very general and can be used in a wide number of applications and problems. The evaluation function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\zeta ^{\mathcal {S}}_{\tau }(\cdot)$\end{document}ζτS(·) for the systems \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {S}$\end{document}S can be considered as an input-output map relating the parameters p to the function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$z, z \in \mathbb {R}$\end{document}z,z∈ℝ, measuring the capability of interest τ: (2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ \zeta:\mathbb{P} \longrightarrow \mathbb{R}, \quad z=\zeta^{\mathcal{S}}_{\tau}{(}{p}{)}. $$ \end{document}ζ:ℙ→ℝ,z=ζτS(p). The function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\zeta ^{\mathcal {S}}_{\tau }(p)$\end{document}ζτS(p) is not necessarily known analytically, and in Systems Biology applications it is often computed in silico through a simulation of the mathematical model of system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {S}$\end{document}S for a given value p of the parameter vector. In this paper we will assume that the biological system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {S}$\end{document}S can be modeled by a system of ordinary differential equations of the form: (3) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ \mathcal{S} : \left\{\begin{array}{rcl} \dot{x} & = & f(x,u,p), \quad x(0)=x_{0}, x \in \mathbb{R}^{n}, \; p \in \mathbb{P} \\ z & = & h(x,p), \quad y \in \mathbb{R}^{m} \end{array}\right. $$ \end{document}S:x˙=f(x,u,p),x(0)=x0,x∈ℝn,p∈ℙz=h(x,p),y∈ℝm 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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {P}$\end{document}ℙ, a subset of the positive orthant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {R}_{>0}^{q}$\end{document}ℝ>0q of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {R}^{q}$\end{document}ℝ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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\zeta ^{\mathcal {S}}_{\tau }(\cdot)$\end{document}ζτS(·), 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(\mathbb {P},\mathcal {A})$\end{document}(ℙ,A). We denote by FP(p) the cumulative distribution function of P, which is a model of the “a priori” knowledge on P. Let fP(p) denote the corresponding probability density function. Similarly, the evaluation function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\zeta ^{\mathcal {S}}_{\tau }(\cdot)$\end{document}ζτS(·) can be interpreted as a transformation from the random variable P to the random variable Z, whose realizations are given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$z=\zeta ^{\mathcal {S}}_{\tau }(p)$\end{document}z=ζτS(p). Hence, let FZ(z) and fZ(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 FP(p) and fP(p) through the mapping ζ. Partitioning definition domain Z The system behavior, in addition to the output Z, can also be characterized by means of the above mentioned density function fZ(z) and by means of the associated probability of having values of the measure Z within a given interval of its definition domain. The red curve in Fig. 3 is an example of probability density function of the model output; i.e., the proliferation activity in a population of cancer cells. Fig. 3 General problem formulation of conditional robustness. The problem addressed 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) having values of the output function in a given lower set L i To further characterize system behavior, we consider some subset T(α) of the definition domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {D}_{Z}$\end{document}DZ of the output Z. For example, in the case of the normal distributions, examples of subsets T(α) are the lower quartile and the upper quartile. Here, we will use a subset T(α) defined according to one of the following conditions: (4a) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$\begin{array}{@{}rcl@{}} T =T(\alpha)&:=& \left\{ z \le a : {\int_{0}^{a}}f_{Z}(z)dz=\alpha\right\}, \end{array} $$ \end{document}T=T(α):=z≤a:∫0afZ(z)dz=α, (4b) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$\begin{array}{@{}rcl@{}} T =T(\alpha)&:=& \left\{ z \ge a : \int_{a}^{\infty}f_{Z}(z)dz=\alpha\right\}. \end{array} $$ \end{document} T = T ( α ) : = z ≥ a : ∫ a ∞ f Z ( z ) dz = α . 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 fZ(z) curve. More generally, we can define a set T(α,a,b) such that: (5) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$\begin{array}{*{20}l} T(\alpha, a, b) :&=\left\{{\vphantom{{\int_{a}^{b}}}} z\in \mathcal{D}_{Z} : Pr\{a0}^{6}$\end{document}ℝ>06. We fixed the number of sample to N=10000 and we generated 100 realizations of the subset PS,i,i=1,…,100, of the parameters space with L2HS. The kernel distribution estimates for the 100 realizations of the three evaluation functions are shown in Fig. 6b, 6c and 6d, 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{P_{i}|U}(p_{i})$\end{document}fPi|U(pi) (red line) and the lower distributions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{P_{i}|L}(p_{i})$\end{document}fPi|L(pi) (blue line) for the whole parameter set, i.e., k1, K1, k12, K2, λ2 and λ, for a single realization of PS. MIRI was evaluated for each parameter (Fig. 7a) and the most significant parameters turn out to be λ, k12 and K2, as can be inferred, for a single realization of PS, from Additional file 1: Figure S1F, Additional file 1: Figure S1C and Additional file 1: Figure S1D, respectively. MIRI for the six model parameters was evaluated for the 100 realizations and the result is presented in the box plot in Fig. 7b, confirming λ, k12 and K2 as the parameters most relevant for the pulse area control. Fig. 7 Moment Independent Robustness Indicator (MIRI) for the Y signal area. a Single realization. b Box plot for 100 realizations 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 λ, k12 and K2 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. Conditional robustness was performed setting each parameter to the values chosen according to the procedure above. Figure 8a shows the conditional probability density estimates setting k12, K2 and λ parameters as shown in Table 1. Also \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|k_{12},K_{2},\lambda }(z)$\end{document}fZ|k12,K2,λ(z) is compared with the unconditional probability density for the evaluation function. The parameter λ plays the major role in maximizing the pulse area and the use of the three parameters selected allows us to obtain the expected behavior. Figure 8b presents an example of simulation where the parameters k12, K2 and λ are fixed and the other parameters are randomly selected: the network generates a version of the Y signal with increased area. Fig. 8 Conditional robustness for Y area (N=10000). a Probability density function. b Simulation examples conditioning parameters as shown in first row of Table 1 Table 1 Parameters setting for the pulse generator and EGFR-IGF1R in silico experiments 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 L2HS 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: k12 and λ. We fixed k12 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|k_{12}}(z)$\end{document}fZ|k12(z): its first moment is higher than the unconditioned distribution fZ(z); the yellow curve is the conditional density function, fZ|λ(z) (see Table 1) and the distribution spread at higher values of the evaluation function; the red line is the fZ|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 k1, K2 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 k1, K2 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). Finally, we performed conditional robustness analysis with the goal of obtaining a fast amplified pulse generation; therefore, we minimized the time to maximum of Y and maximized both the area and the intensity of Y. For each of the above indexes we chose the parameter with the higher MIRI according to the previous analysis (see Fig. 7b, Additional file 3: Figure S3B and Additional file 5: Figure S5B). K2 minimizes the time to maximum of Y and k12 and λ maximizes the area and the intensity of the pulse Y, respectively. The modes are presented in the forth row of Table 1. Figure 9 shows the results of in silico experiment through conditional probability densities for the time to maximum of Y (Fig. 9a), the intensity of Y (Fig. 9b) and the area of Y (Fig. 9c). An example of time behavior of the model states with the selected values for the k12, K2 and λ parameters is shown (see Fig. 9d). Fig. 9 Multiobjective conditional robustness for area, time to maximum and amplification of Y (N=10000 and parameters values are shown in the forth row of Table 1). a Pdf for time to maximum of Y. b Pdf for area of Y. c Pdf for Intensity of Y. d Simulations example with fixed k 12, K 2 and λ EGF-IGF network in lung cancer The proposed algorithm has been applied to the EGFR-IGF1R network, which is our key motivating example [44]. This network is important in cancer pathogenesis and progression, mostly in the development of Non Small Cell Lung Cancer (NSCLC). Although this signaling pathway is a therapeutic target, recent clinical trials have exhibited limited effects [45]. These complex networks include interactions between epidermal growth factor receptor (EGFR), insulin-like growth factor 1 receptor (IGF1R) along with their downstream effectors of the Mitogen-activated protein kinases (MAPK) cascade and the phospholipids inositol kinases axis (PIK3). One of the main downstream effectors of this network is ERK, a kinase able to phosphorylate both cytosolic and nuclear substrates, including several transcription factors (see pathway in Fig. 10a and [20]). Fig. 10 The EGFR-IGF1R pathways in NSCLC. a Pathways graph as presented in [46]. b State variables representation of the EGFR-IGF1R pathways The biochemical network is sketched in Fig. 10b. The activation of the EGFR and IGF1R (x1 and x2, respectively) induced by the ligands binding triggers a cascade of complex proteins interactions that leads to the activation of the ERK protein (x7). 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]. Here, we rewrite the model in [20] considering the conservation law and assuming (as is common in this field) that the total amount of each protein and receptor is constant and equal to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${x_{i}^{T}}$\end{document}xiT, i=1,⋯,10, as presented in [48]. The complete system has 10 state variables, 3 input signals and a parameter vector comprising 39 entries. The output function z of interest, the evaluation function, is the area under the ERK∗ curve over the whole time span studied (30 min). This evaluation function is selected because ERK∗ can be measured with several experimental methodologies such a quantitative immunoblotting, immunofluorescence, in live-cell by means of fluorescence microscopy and mass spectrometry [49] and in patient samples with immunohistochemistry [50]. In [11] the measure of ERK concentration is used to infer the silencing of its activity in tumor cells. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ \begin{array}{rcl} \dot{x}_{1} \!& \!= \! &\! -p_{1}\,{x_{1}}\\ \dot{x}_{2} \!& \!= \! &\! -p_{2}\,{x_{2}}\\ \dot{x}_{3} \!& \!= \! &\! p_{6}{x_{1}}\frac{{{x^{T}_{3}}} - {x_{3}}}{p_{7}+{{x^{T}_{3}}} - {x_{3}}} + p_{14}{x_{2}}\frac{{{x^{T}_{3}}} - {x_{3}}}{p_{15}+{{x^{T}_{3}}} - {x_{3}}} \\ &&- p_{12}{x_{13}}\frac{{x_{3}}}{p_{13}+{x_{3}}} \\ \dot{x}_{4} \!& \!= \! &\! p_{8}{x_{3}}\frac{{{x^{T}_{4}}}-{x_{4}}}{p_{9}+{{x^{T}_{4}}}-{x_{4}}} - p_{33}{u_{3}}\frac{{x_{4}}}{p_{34}+{x _{4}}}\\ \dot{x}_{5} \!& \!= \! &\! p_{27}{x_{4}}\frac{{{x^{T}_{5}}}-{x_{5}}}{p_{28}+{{x^{T}_{5}}}-{x_{5}}} - p_{37}{u_{1}}\frac{{x_{7}}}{p_{38}+{x_{7}}} \\ && -p_{31}{x_{10}}\frac{{x_{5}}}{p_{32}+{x_{5}}}\\ \dot{x}_{6} \!& \!= \! &\! p_{29}{x_{5}}\frac{{{x^{T}_{6}}}-{x_{6}}}{p_{30}+{{x^{T}_{6}}}-{x_{6}}} - p_{35}{u_{2}}\frac{{x_{6}}}{p_{36}+{x_{6}}}\\ \dot{x}_{7} \!& \!= \! &\! p_{10}{x_{6}}\frac{{{x^{T}_{7}}}-{x_{7}}}{p_{11}+{{x^{T}_{7}}}-{x_{7}}} - p_{23}{u_{2}}\frac{{x_{7}}}{p_{24}+{x_{7}}}\\ \dot{x}_{8} \!& \!= \! &\! p_{4}{x_{7}}\frac{{{x^{T}_{8}}}-{x_{8}}}{p_{5}+{{x^{T}_{8}}}-{x_{8}}} - p_{39}{x_{8}}\\ \dot{x}_{9} \!& \!= \! &\! p_{25}{x_{4}}\frac{{{x^{T}_{9}}}-{x_{9}}}{p_{26}+{{x^{T}_{9}}}-{x_{9}}}+ p_{16}{x_{2}}\frac{{{x^{T}_{9}}}-{x_{9}}}{p_{17}+{{x^{T}_{9}}}-{x_{9}}}\\ &&+ p_{18}{x_{1}}\frac{{{x^{T}_{9}}}-{x_{9}}}{p_{19}+{{x^{T}_{9}}}-{x_{9}}}-p_{3}{x_{9}}\\ \dot{x}_{10} \!& \!= \! &\! +p_{20}{x_{9}}\frac{{x^{T}_{10}}-{x_{10}}}{p_{21}+{x^{T}_{10}}-{x_{10}}}-p_{22}{x_{10}}\\ \end{array} $$ \end{document}x˙1=−p1x1x˙2=−p2x2x˙3=p6x1x3T−x3p7+x3T−x3+p14x2x3T−x3p15+x3T−x3−p12x13x3p13+x3x˙4=p8x3x4T−x4p9+x4T−x4−p33u3x4p34+x4x˙5=p27x4x5T−x5p28+x5T−x5−p37u1x7p38+x7−p31x10x5p32+x5x˙6=p29x5x6T−x6p30+x6T−x6−p35u2x6p36+x6x˙7=p10x6x7T−x7p11+x7T−x7−p23u2x7p24+x7x˙8=p4x7x8T−x8p5+x8T−x8−p39x8x˙9=p25x4x9T−x9p26+x9T−x9+p16x2x9T−x9p17+x9T−x9+p18x1x9T−x9p19+x9T−x9−p3x9x˙10=+p20x9x10T−x10p21+x10T−x10−p22x10 Following an approach widely used in literature, we assume that each element pi of the parameter vector p, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$p \in \mathbb {R}^{39}$\end{document}p∈ℝ39, can assume values in a range defined by two multiplicative bounds w.r.t. a nominal (wild-type) value pwt,i: a lower bound bℓ,i=cℓ,ipwt,i and an upper bound bu,i=cu,ipwt,i, where the coefficients cℓ,i and cu,i are the multiplicative perturbations, cℓ,i<1, ∀i and cu,i<1, ∀i. Hence, the parameter space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {P}$\end{document}ℙ is given by the cartesian product \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {P} =\prod _{i=1}^{39}[c_{\ell,i} p_{wt,i}, \, c_{u,i} p_{wt,i}]$\end{document}ℙ=∏i=139[cℓ,ipwt,i,cu,ipwt,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). We apply conditional robustness analysis for the model in (9) generating in silico measures of the active ERK. The parameter space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {P}$\end{document}ℙ was sampled with the L2HS method, the number of samples and the number of realizations of the parameter space sampling were fixed to N=10000 and 100, respectively. Figure 11a shows a realization of MIRI for the EGF-IGF model. The highest MIRI values are achieved by p23, p24, p35 and p36: they represent the activation rates and the Michealis-Menten constants of protein phosphatase 2A (PP2A) with ERK and MAPK/ERK kinase (MEK), respectively. The activation rate of Ras and Raf interaction (p27) also has relevant MIRI. The MIRI box plot for 100 realizations confirmed the MIRI ordering for the 39 parameters (Fig. 11b). Fig. 11 Moment Independent Robustness Indicator (MIRI) for ERK activity in EGFR-IGF1R model. a Single realization. b Box plot for 100 realizations The estimated kernel density for the values of the area under ERK∗ signal is shown in Fig. 12a: it has a mean, a variance and a mode presented in Table 2. We evaluate the conditional robustness for the single value p23, p24,p35, p36 and p27 obtained fixing the lower L and upper U set with a parameter α=0.1, and setting them to the highest density values for the corresponding densities conditioned on set L, (see last line in Table 1). We evaluate conditional robustness fixing only the activation rates p23, p35 and p27 and also fixing all 5 parameters selected. Figure 12b shows the conditional densities for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{23}}(z)$\end{document}fZ|p23(z), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{35}}(z)$\end{document}fZ|p35(z), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{27}}(z)$\end{document}fZ|p27(z), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{24}}(z)$\end{document}fZ|p24(z), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{36}}(z)$\end{document}fZ|p36(z), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{23},p_{35},p_{27}}(z)$\end{document}fZ|p23,p35,p27(z) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{Z|p_{23},p_{24},p_{35},p_{36}p_{27}}(z)$\end{document}fZ|p23,p24,p35,p36p27(z) and Table 2 presents the statistical descriptors for the conditional densities. The lowest mean, variance are obtained conditioning all the selected parameters. Figure 12c shows a comparison between the wild type simulation of ERK activity (black line) and the conditioned simulation with parameters with fixed value as in Table 1 (red line). The ERK∗ activity is strongly reduced as showed in Table 2. Fig. 12 a Unconditional distribution of evaluation function E R K ∗. b Conditional robustness for ERK activity (N=10000). c E R K ∗ time simulation: wild type (black line) vs simulation at conditioned value as in Table 1 Table 2 Proliferation output pdf measures It is known from biochemical and clinical literature that the concentration of the receptors EGFR and IGF1R has a major role in the proliferation attitude (see, among others, [20, 44] and the references therein). To exploit this characteristic, we also investigated perturbation on the initial value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${x_{1}^{0}}=EGFR$\end{document}x10=EGFR and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${x_{2}^{0}}=IGF1R$\end{document}x20=IGF1R. We use a receptor space X1,2, obtained in a manner similar to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb {P}$\end{document}ℙ, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${X}_{1,2} := \left [c_{\ell,1}x^{0}_{wt,1}, \,c_{u,1}x_{wt,1}1^{0}\right ]\,\times \, \left [c_{\ell,2}x_{wt,2}^{0}, \,c_{u,2}x_{wt,2}^{0}\right ]$\end{document}X1,2:=cℓ,1xwt,10,cu,1xwt,110×cℓ,2xwt,20,cu,2xwt,20, where the c’s coefficients are the perturbations and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$x_{wt,i}^{0}$\end{document}xwt,i0, i=1,2, are the wild-type initial conditions. Figure 13a shows the unconditioned measure distribution with perturbed EGFR and IGF1R (red line) that shows higher mean and variance than is the case with fixed initial conditions (see Table 2). The MIRI box plot in Fig. 13b, generated with our algorithm by perturbing both the 39 parameters and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${x_{1}^{0}}$\end{document}x10 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${x_{2}^{0}}$\end{document}x20 initial conditions, confirmed the same conditioning set obtained in the previous analysis. Conditioning the p23, p24,p35, p36 and p27 parameters we are still able to reduce the mean and variance (see Table 2 and Fig. 13c)). Fig. 13 a Unconditional E R K ∗ activity at EGFR and I G F1R initial condition wild type and at perturbed EGFR and I G F1R. b MIRI of the parameters for 100 realizations. c Effect of conditioning of p 23, p 24,p 35, p 36 and p 27 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 (L2HS) 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- IGF1R 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. 13b). 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 12b can be interpreted as the best combination of target to be used and Fig. 13b 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 cu,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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {Z}$\end{document}Z of the evaluation function; the probability α defining the subsets of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {Z}$\end{document}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 fZ(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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{P_{i}|L}(p)$\end{document}fPi|L(p) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{P_{i}|U}(p)$\end{document}fPi|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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\frac {1.000}{\alpha }$\end{document}1.000α. 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 fZ(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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{P_{i}|L}(p)$\end{document}fPi|L(p) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$f_{P_{i}|U}(p)$\end{document}fPi|U(p) that turn out to be sufficiently different and to the overall effect on the conditional density fZ|K(z). The computational cost of our CRA has two main sources: the L2HS 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {Z}$\end{document}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.