 Software
 Open Access
CRA toolbox: software package for conditional robustness analysis of cancer systems biology models in MATLAB
 Fortunato Bianconi^{1}Email author,
 Chiara Antonini†^{2},
 Lorenzo Tomassoni†^{2} and
 Paolo Valigi^{2}
 Received: 26 April 2019
 Accepted: 5 June 2019
 Published: 9 July 2019
Abstract
Background
In cancer research, robustness of a complex biochemical network is one of the most relevant properties to investigate for the development of novel targeted therapies. In cancer systems biology, biological networks are typically modeled through Ordinary Differential Equation (ODE) models. Hence, robustness analysis consists in quantifying how much the temporal behavior of a specific node is influenced by the perturbation of model parameters. The Conditional Robustness Algorithm (CRA) is a valuable methodology to perform robustness analysis on a selected output variable, representative of the proliferation activity of cancer disease.
Results
Here we introduce our new freely downloadable software, the CRA Toolbox. The CRA Toolbox is an ObjectOriented MATLAB package which implements the features of CRA for ODE models. It offers the users the ability to import a mathematical model in Systems Biology Markup Language (SBML), to perturb the model parameter space and to choose the reference node for the robustness analysis. The CRA Toolbox allows the users to visualize and save all the generated results through a userfriendly Graphical User Interface (GUI). The CRA Toolbox has a modular and flexible architecture since it is designed according to some engineering design patterns. This tool has been successfully applied in three nonlinear ODE models: the Prostatespecific Pten^{−/−} mouse model, the Pulse Generator Network and the EGFRIGF1R pathway.
Conclusions
The CRA Toolbox for MATLAB is an opensource tool implementing the CRA to perform conditional robustness analysis. With its unique set of functions, the CRA Toolbox is a remarkable software for the topological study of biological networks. The source and example code and the corresponding documentation are freely available at the web site: http://gitlab.ict4life.com/SysBiOThe/CRAMatlab.
Keywords
 Ordinary differential equation models
 Conditional robustness analysis
 MATLAB package
 Signaling networks
Background
In Systems Biology, mathematical modeling and computational software are important tools to unravel the complexity of biological systems and predict their behavior under different perturbations [1]. Typically, many models consist of a set of Ordinary Differential Equations (ODEs) which allow understanding and reproducing the dynamic behavior of molecular interactions through simulations and integration of the ODEs [2]. To support mathematical modeling of biological networks, the use of software tools has grown substantially in recent years. These software are designed to assist the users at different stages of the modeling process, from model generation to parameter estimation and model analysis.
In cancer research, the use of systems biology approaches is particularly useful to elucidate mechanisms of tumorigenesis and tumor resistance. Computational predictive models, integrated with patient data, help scientists in the validation of new and durable therapies [3]. In order to discover effective and targeted drugs, robustness is one of the most relevant properties of cancer signaling networks to investigate. Robustness is defined as the ability of a biological system to maintain its functionalities against internal and external perturbations [4]. In more detail, cancer robustness is a quantitative measure of the tumor proliferation attitude against extracellular inputs. Thus, understanding new ways to reduce robustness of the cell proliferation activity is a key issue in cancer systems biology. Since cell growth is driven by protein interaction networks, the proliferation activity can be quantified by looking at the activation of a protein involved in the proliferation process [5]. In mathematical modeling, this can be done by perturbing model parameters and analyze how the concentration of the protein of interest changes over time. An algorithm developed for this purpose is the Conditional Robustness Algorithm (CRA) proposed in [5]. This algorithm, through computational perturbations and simulations, identifies a small number of nodes in the cancer network which influences most the activity of the proliferation indicator. As a result, by conditioning these nodes with specific drugs, it may be possible to reduce the tumor robustness.
Robustness of mathematical model problem is well studied in literature except when the models are based on nonlinear ODE. A class of methods is aimed at analyzing the geometry and the volume of the feasible region, which is the region in the parameter space that allows the system to properly work [6]. Moreover, there are other algorithms that infer the robustness of a model looking at its topology, such as the number and the structure of positive and negative loops [7]. However, these techniques are typically applied to mathematical models whose parameters are not kinetic. Another and different category of methods is the Global Sensitivity Analysis (GSA) class of algorithms [8]. They are similar to CRA and are involved in the analysis of the uncertainties in kinetic model parameters through the sampling of the parameter space. Despite that, GSA and CRA have clear distinct goals. GSA is typically interested in the variation of a performance index with the respect to the model parameters. Since many times this is basically implemented through the derivatives, GSA tools are useful when an optimization of the system is required [9]. In the CRA algorithm, on the other hand, the purpose is not to maximize or minimize a certain optimization function, but the main interest is in finding the conditioning set, which is the subset of parameters that more likely is able to impose a specific behavior to one output of the model.
In order to facilitate the study of cancer robustness and the application of CRA, we develop the CRA Toolbox, a software package for MATLAB. It is an open source tool which allows non expert users to apply the CRA to any ODE model in a simple and quick way. The CRA Toolbox consists in a easy to use Graphical User Interface (GUI) and in a set of functions which can be easily extended by the users in order to achieve specific requirements.
In the following subsection, we briefly describe the CRA method implemented in the toolbox. A detailed description of the method can be found in [5].
Conditional Robustness Algorithm (CRA)

simultaneous perturbation of model parameters;

definition and computation of the evaluation function;

estimation of the conditional probability density functions (pfds) for each model parameter.
where the index i represents the selected output variable of the model. Thus, the evaluation function, that depends on the time behavior of y_{i} (which, in turn, depends on the selected parameter vector p), can be considered as a user defined summary statistic that stands for a specific property of the chosen output node y_{i}. The set of computed evaluation functions, for a specific output node y_{i}, has cardinality equal to N_{S}, i.e. the number of sampled parameter vectors. Let denote with\(f_{Z_{i}}(z_{i})\) the pdf of the set of evaluation functions previously defined.
where α is the level of probability that represents the area under the lower and upper tail of \(f_{Z_{i}}(z_{i})\) and a is the corresponding threshold value in the domain of the pdf.
The MIRI is an index that measures the level of separation between \(f_{p_{i}U}\) and \(f_{p_{i}L}\) for each parameter of the model. An high value of the MIRI for a parameter p_{i} means that the perturbation of the parameter space along the p_{i} direction leads to high variation of the evaluation function. Thus, the higher the value of a MIRI, the higher is the influence of that parameter on the dynamical behavior of the selected output node.
Finally, the output of CRA is the vector μ that contains the value of the MIRI associated to each parameter of the model. For further details about CRA, see [5].
Implementation
The CRA Toolbox is an open source software developed as a MATLAB package. Its core is a set of functions that the user can run locally in a MATLAB environment by downloading the folder containing the toolbox. This software automates all the functionalities required by CRA to perform robustness analysis of an ODE model.
The CRA Toolbox includes a GUI that runs within MATLAB to encourage the employment of the software also for non expert users. In more detail, the tool firstly performs the import of a mathematical model written in Systems Biology Markup Language (SBML) and saved in .xml file format. Then, it allows the user to set the tuning parameters to regulate the parameter space perturbation and the model integration, such as the specific ODE solver to use. Before selecting the reference node from a scrollbar that lists all the outputs of the model, it is necessary to start the simulations by clicking on a specific button. Once the in silico measures are completed, the tool requires the selection of a specific evaluation function in a predefined list and the method for the computation of the lower and higher tail of the pdf of the evaluation function. Finally, it is possible to plot and save in a user defined directory all the in silico measures, the estimated pdfs and the boxplot of MIRIs. In order to guarantee the reliability of results, the toolbox supports the generation of multiple realizations of the entire procedure and of the resulting MIRIs and pdfs. In order to speed up the model simulation we use parallel processing through the Parallel Computing Toolbox ^{TM} [11].
Moreover, we also provide an alternative implementation of the CRA Toolbox that allows the user to run the algorithm in batch mode directly from the command line. The core functions and the architecture of the software remains unchanged, but for this version we removed the GUI and we also avoided the use of Simbiology to enhance the portability of the code. Indeed, in this version of the software, the mathematical model can be given in input as a Matlab function where all the ODEs are specified and it is not required to use the SBML language and the corresponding Simbiology Object.
Results
Prostatespecific Pten ^{−/−} mouse model
List of the kinetic parameters of the Prostatespecific Pten^{−/−} mouse model and their corresponding nominal values [13]
Parameter name  Value 

λ _{ A}  0.6931 
r _{ p1}  0.323432 
r _{ a1}  0.1 
r _{ m}  6.702665 
r _{ p2}  1.426189 
r _{ a2}  1.061691 
α _{ XD}  2.096368 
α _{ DC}  3.077551 
α _{ DR}  0.739121 
α _{ XR}  0.144144 
α _{ CI}  0.408862 
α _{ IR}  0.354526 
\(\alpha _{D_{C}D_{R}}\)  0.1 
\(\alpha _{D_{R}R}\)  0.4722 
k _{ CX}  0.1 
k _{ RC}  1.459617 
π _{ D}  5.051123 
π _{ C}  4.111938 
π _{ R}  0.1 
μ _{ C}  0.978832 
μ _{ R}  0.564691 
μ _{ I}  0.213819 
μ _{ D}  0.650308 
List of the initial conditions of the Prostatespecific Pten^{−/−} mouse model state variables and their corresponding values [13]
State variable  Species  Value 

A  Androgen  1 
X _{1}  CSPC  1 
X _{2}  CRPC  0 
D _{ m}  Mature DC  1 
C _{2}  CTL in prostate  1 
R _{2}  Treg in prostate  1 
I _{2}  IL2 in prostate  1 
D _{ C}  Functional DC  1 
D _{ R}  Regulatory DC  1 
C _{1}  CTL in lymphoid  1 
R _{1}  Treg in lymphoid  1 
I _{1}  IL2 in prostate  1 
Table 5 contains the time necessary to perform all the simulations.
Pulse generator network
EGFRIGF1R pathway in lung cancer
List of species included in the EGFRIGF1R model and the corresponding initial concentrations and total protein amount
State Variable  Species  Value 

x _{1}  EGFR ^{∗}  8000 
x _{2}  IGF1R ^{∗}  650 
x _{3}  SOS  0 
x _{4}  Ras ^{∗}  0 
x _{5}  Raf ^{∗}  0 
x _{6}  MEK ^{∗}  0 
x _{7}  Erk ^{∗}  0 
x _{8}  p90Rsk ^{∗}  0 
x _{9}  PIK3 ^{∗}  0 
x _{10}  Akt ^{∗}  0 
u _{1}  RafPP  120000 
u _{2}  PP2A  120000 
u _{3}  RasGab  120000 
\(x_{3}^{T}\)  DSOS  120000 
\(x_{4}^{T}\)  Ras  120000 
\(x_{5}^{T}\)  Raf  120000 
\(x_{6}^{T}\)  MEK  600000 
\(x_{7}^{T}\)  Erk  600000 
\(x_{8}^{T}\)  p90Rsk  120000 
\(x_{9}^{T}\)  PIK3  120000 
\(x_{10}^{T}\)  Akt  120000 
List of the kinetic parameters of the EGFRIGF1R model and their corresponding nominal values
Parameter  Name  Value 

p _{1}  γ _{ EGFR}  0.02 
p _{2}  γ _{ IGF1R}  0.02 
p _{3}  k d _{ PIK3}  0.005 
p _{4}  k _{ p90RskErk}  0.0213697 
p _{5}  K M _{ p90RskErk}  763523 
p _{6}  k _{ SOSEGFR}  694.731 
p _{7}  KM _{ SOSEGFR}  6086070 
p _{8}  k _{ RasSOS}  32.344 
p _{9}  KM _{ RasSOS}  35954.3 
p _{10}  k _{ ErkMEK}  9.85367 
p _{11}  KM _{ ErkMEK}  1007340 
p _{12}  k _{ DSOSp90Rsk}  161197 
p _{13}  KM _{ DSOSp90Rsk}  896896 
p _{14}  k _{ SOSIGFR}  500 
p _{15}  KM _{ SOSIGFR}  1000000 
p _{16}  k _{ PIK3IGF1R}  10.6737 
p _{17}  KM _{ PIK3IGF1R}  184912 
p _{18}  k _{ PIK3EGFR}  10.6737 
p _{19}  KM _{ PIK3EGFR}  184912 
p _{20}  k _{ AktPIK3}  0.0566279 
p _{21}  KM _{ AktPIK3}  653951 
p _{22}  kd _{ Akt}  0.005 
p _{23}  k _{ ErkPP2A}  9.85367 
p _{24}  KM _{ ErkPP2A}  1007340 
p _{25}  k _{ PIK3Ras}  0.0771067 
p _{26}  KM _{ PIK3Ras}  272056 
p _{27}  k _{ RafRas}  0.884096 
p _{28}  KM _{ RafRas}  62464.6 
p _{29}  k _{ RafMEK}  185.759 
p _{30}  KM _{ RafMEK}  4768350 
p _{31}  k _{ RafAkt}  15.1212 
p _{32}  KM _{ RafAkt}  119355 
p _{33}  k _{ RasRasGab}  1509.36 
p _{34}  KM _{ RasRasGab}  1432410 
p _{35}  k _{ MEKPP2A}  2.83243 
p _{36}  KM _{ MEKPP2A}  518753 
p _{37}  k _{ RafRafPP}  0.126329 
p _{38}  KM _{ RafRafPP}  1061.71 
p _{39}  kd _{ p90Rsk}  0.005 
Computational cost to run the CRA algorithm in all the examples provided in the Results section
Model  Output node  Evaluation Function  Time (sec.) 
Area  45  
Pulse Generator Network  Y  Maximum  45 
Time of Maximum  42  
Area  1702  
C2  Maximum  1703  
Prostatespecific Pten^{−/−} mouse model  Time of Maximum  1706  
C1  Area  1708  
EGFRIGF1R pathway in lung cancer  ERK  Area  102 
Discussion
The CRA is an algorithm to study the robustness of complex biological networks and it allows the identification of few parameters that have a major impact on the behavior of a selected output variable. One of its main innovations is the introduction of a sensitivity measure, the MIRI, that takes advantage of all the conditional distributions of the parameters, without reference to a specific moment. In [18] there are the mathematical details of this class of indicators and the comparison with the variancebased uncertainty importance measures.
Here we present the CRA Toolbox, a software package for MATLAB aimed at performing the robustness analysis based on the paradigm of the CRA. The tool consists of a set of MATLAB functions that automate all the necessary mathematical steps from the integration of an ODE model until the calculus of the MIRIs. We decide to use the ObjectOriented programming paradigm because it allows us the development of an extendable and flexible architecture through the implementation of different engineering design patterns. In this way, other users can add blocks of software to define novel evaluation functions and other methods for the identification of the pdf tails without modifying the structure of the software and the source code of the existing classes.
Here we tested the CRA on three ODE models that contain different kinetic laws and different number of state variables. All the simulations were performed on a Intel Core i74700HQ CPU, 2.40GHz 8, 16GB memory, Ubuntu 16.04 LTS (64bit). From Table 5 it is clear how the time to complete a realization of the CRA highly depends both on the dimension of the model and on the number of samples N_{S} of the LHS. Indeed, from a computational point of view, the most intensive part of the algorithm is the integration of the mathematical model N_{S} times. Moreover, the choice of N_{S} influences the estimation of the evaluation function pdf, whose cost increases with N_{S}. On the other hand, the computation of the MIRIs is pretty fast, independently of the type of evaluation function selected. For more details about the choice of the tuning parameters and the computational costs of the CRA for different settings of tuning parameters see Additional file 4.
One of the key ideas in [5] is that the theoretical statements of the CRA do not depend on the specific modeling technique used to represent a biological phenomenon. For this reason, one of the possible future development of the tool is to augment the number and types of formats of mathematical models taken in input. Finally, the code of the CRA Toolbox can be easily adapted to other open source programming languages such as Octave and Python.
Conclusions
The CRA Toolbox is unique in the category of the robustness analysis tools because it is the specific implementation of the CRA with all its features. It has the scope to enlarge and facilitate the usage of this algorithm and to disclosure it also to non expert users. This can significantly help the oncological research of physicians in discover novel targeted therapies. Moreover, in [19], the Conditional Robustness Calibration (CRC) algorithm is presented, which is an upgrade of the CRA that allows the generation of MIRIs for a given model taking into account multiple nodes simultaneously. This novel algorithm suggests that the CRA Toolbox can be modified and extended in the future in order to implement also the functionalities of CRC.
Availability data and requirements
Project name: CRA Toolbox
Project home page: http://gitlab.ict4life.com/SysBiOThe/CRAMatlab
Operating system(s): Platform independent
Programming language: MATLAB
Other requirements: MATLAB
License: MATLAB
Any restrictions to use by nonacademics: License needed for MATLAB
Notes
Declarations
Acknowledgements
The authors would like to thank ICT4life srl group for the helpful discussions about implementation, potential extensions/uses of modules within the algorithm and make freely available a public project with GitLab ICT4life repository.
We thank Italian Association of Cancer Research (AIRC) for supporting this work.
Funding
CA and TL are funded by Italian Association for Cancer Research (AIRC) grant 15713/2014. The funding body had no role in the design of the study and collection, analysis, interpretation of data and in writing the manuscript.
Authors’ contributions
FB, CA, and LT conceived and designed the toolbox. FB, CA and LT wrote the code. CA and LT wrote the user’s documentation and developed the examples. FB, CA, LT and PV wrote the paper. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable
Competing interests
The authors declare that they have no competing interests.
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
 Kitano H. Systems biology: a brief overview. Science. 2002; 295(5560):1662–4.PubMedView ArticleGoogle Scholar
 Ghosh S, Matsuoka Y, Asai Y, Hsin KY, Kitano H. Software for systems biology: from tools to integrated platforms. Nat Rev Genet. 2011; 12(12):821.PubMedView ArticleGoogle Scholar
 Werner HM, Mills GB, Ram PT. Cancer systems biology: a peek into the future of patient care?. Nat Rev Clin Oncol. 2014; 11(3):167.PubMedPubMed CentralView ArticleGoogle Scholar
 Kitano H. Biological robustness. Nat Rev Genet. 2004; 5(11):826.PubMedView ArticleGoogle Scholar
 Bianconi F, Baldelli E, Luovini V, Petricoin EF, Crinò L, Valigi P. Conditional robustness analysis for fragility discovery and target identification in biochemical networks and in cancer systems biology. BMC Syst Biol. 2015; 9(1):70.PubMedPubMed CentralView ArticleGoogle 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.View ArticleGoogle Scholar
 Kwon YK, Cho KH. Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics. Bioinformatics. 2008; 24(7):987–94.PubMedView ArticleGoogle Scholar
 Iooss B, Lemaître P. A review on global sensitivity analysis methods. In: Uncertainty management in simulationoptimization of complex systems: 2015. p. 101–22.View ArticleGoogle 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.View ArticleGoogle Scholar
 Kitano H. Towards a theory of biological robustness. Mol Syst Biol. 2007; 3(1):137.PubMedPubMed CentralGoogle Scholar
 Perform Parallel Computations on Multicore Computers, GPUs, and Computer Clusters. https://it.mathworks.com/products/parallelcomputing.html.
 Helm R, Johnson RE, Gamma E, Vlissides J. Design Patterns: Elements of Reusable Objectoriented Software: Addison Wesley; 2000.Google Scholar
 Peng H, Zhao W, Tan H, Ji Z, Li J, Li K, Zhou X. Prediction of treatment efficacy for prostate cancer using a mathematical model. Sci Rep. 2016; 6:21599.PubMedPubMed CentralView ArticleGoogle Scholar
 Szallasi Z, Stelling J, Periwal V. System Modeling in Cellular Biology: the MIT Press; 2006.Google Scholar
 Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007; 8(6):450.PubMedView ArticleGoogle Scholar
 Bianconi F, Baldelli E, Ludovini V, Crinò L, Flacco A, 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.PubMedView ArticleGoogle Scholar
 Li C, Donizelli M, Rodriguez N, Dharuri H, Endler L, Chelliah V, Li L, He E, Henry A, Stefan MI, et al.Biomodels database: An enhanced, curated and annotated resource for published quantitative kinetic models. BMC Syst Biol. 2010; 4(1):92.PubMedPubMed CentralView ArticleGoogle Scholar
 Borgonovo E. A new uncertainty importance measure. Reliab Eng Syst Saf. 2007; 92(6):771–84.View ArticleGoogle Scholar
 Bianconi F, Antonini C, Tomassoni L, Valigi P. Robust calibration of high dimension nonlinear dynamical models for omics data: An application in cancer systems biology. IEEE Trans Control Syst Technol. 2018; 99:1–12.View ArticleGoogle Scholar