CRA toolbox: software package for conditional robustness analysis of cancer systems biology models in MATLAB.

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 Object-Oriented 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 user-friendly 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 Prostate-specific Pten-/- mouse model, the Pulse Generator Network and the EGFR-IGF1R pathway.
CONCLUSIONS: The CRA Toolbox for MATLAB is an open-source 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/CRA-Matlab .

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)

The main underlying theoretical principle of CRA is the definition of conditional robustness proposed by Kitano in [10]: it is the quantitative measure of the ability of a system S to maintain a specific property τ against some perturbations of the parameter vector p of S. The mathematical formulation is the following: where f(p) is the probability density function of p, P is the parameter space and is a function that quantifies and represents the property τ that is under investigation.

The CRA is a stochastic approach for performing conditional robustness analysis of mathematical models, such as ODE models representing biochemical interaction networks. Its purpose is to quantify the influence of each model parameter on the behavior of a specific output node. Let denote with S the following ODE system: where x∈R is the state space vector that contains the species included in the biological model under study; p∈R denotes the parameter vector; u∈R and y∈R are the input and output vectors respectively. The key features of this algorithm are:

simultaneous perturbation of model parameters;

definition and computation of the evaluation function;

estimation of the conditional probability density functions (pfds) for each model parameter.

The parameter vector p is perturbed through Latin Hypercube sampling (LHS) and the model S is integrated for each one of the N samples generated in R. This procedure allows the collection of N vectors of the observables y. Then, the CRA is based on the definition and the in silico computation of an evaluation function on a specific output node i.e., on a specific observable y. In more details, the evaluation function can be formalized as follows:

where the index i represents the selected output variable of the model. Thus, the evaluation function, that depends on the time behavior of y (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. The set of computed evaluation functions, for a specific output node y, has cardinality equal to N, i.e. the number of sampled parameter vectors. Let denote with the pdf of the set of evaluation functions previously defined.

The CRA algorithm aims at quantifying the influence of each model parameter on a specific output node through . In more detail, it is interested in the estimation of the distribution of the parameter vector p only when the lower and upper tail of are selected. To this purpose, the domain of can be partitioned into two regions by the definition of L(α) and U(α) as:

where α is the level of probability that represents the area under the lower and upper tail of and a is the corresponding threshold value in the domain of the pdf.

The partition defined in the domain of allows the estimation of two conditional pdfs for each parameter, and , respectively. These two pdfs are the distributions of the parameters of S when the values of the evaluation function belong to the lower and upper tail of respectively. Here the purpose of the estimation of and is to select a subset of the N samples of the parameters that give rise to the most divergent behaviors of the evaluation function. The two conditional densities defined above are employed for the calculus of the Moment Independent Robustness Indicator (MIRI) according to the following formula:

The MIRI is an index that measures the level of separation between and for each parameter of the model. An high value of the MIRI for a parameter p means that the perturbation of the parameter space along the p 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 [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.

The source code of the CRA Toolbox is written according to the Object-Oriented programming paradigm as it is shown in the UML class diagram in Fig. 1. For a detailed description of all the components of the tool see (Additional file 1). The architecture of the software is modular because we implement it using software engineering design patterns to model relationships and interactions between classes. As an example, we use the behavioral Strategy pattern [12] twice between the three main classes of the tool: once between the TimeBehavior and EvaluationFunction classes and the other between TimeBehavior and Tail classes. This makes the code easily to extend because if a user wants to add another type of evaluation function or wants to implement another method for the pdf tail calculus, he does not have to change any other part of the code.

Fig. 1

UML class diagram. Classes implemented for the development of the CRA Toolbox. The diagram shows the relations between the different classes, the signatures of the methods and the applied design patterns

Results

In this section we show how to use the CRA Toolbox and we report the results of the application of CRA to three different ODE models: the Prostate-specific Pten−/− mouse model, the Pulse Generator Network and the EGFR-IGF1R pathway. The second and the third examples are used in order to verify that the CRA Toolbox produces results in agreement with those in [5]. Figure 2, Additional files 2 and 3 contain a flowchart and two images to guide the user in the use of the tool.

Fig. 2

Flowchart of the CRA Toolbox.The flowchart sums up the different steps necessary to perform the CRA

Prostate-specific Pten−/− mouse model

In this section, we use the ODE model proposed in [13] to illustrate the functionalities of the CRA Toolbox. This model was developed to study the interactions between prostate cancer and immune microenvironment. In more detail, it is a prostate-specific Pten−/− mouse model for analyzing the effect of the combined therapies with vaccines and Androgen deprivation therapy (ADT) in prostate tumor. The original model consists of two compartments, prostate and lymphoid, 11 state variables and four types of therapeutic strategies, resulting in a system of 15 ODEs and 29 parameters. The pathway of the model is shown in Fig. 3. In this example, we analyze a simplified version of this model because we consider only the off-treatment condition (sham-castration). As a result, the mathematical model consists of 12 ODEs and 23 parameters. Equations of the model are:

Fig. 3

Pathway of the Prostate-specific Pten−/− mouse model [13]

Parameter values and initial conditions of state variables are shown, respectively, in Tables 1 and 2. Since initial condition of androgen is set to 1, i.e. A0=1, the concentration of androgen keeps unchanged. Moreover, parameters p, p and p represent probabilities fixed all to 0.5 in [13] and thus they are not perturbed in the CRA procedure.

Table 1

List of the kinetic parameters of the Prostate-specific 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
α_D_CD_R	0.1
α_D_RR	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

Table 2

List of the initial conditions of the Prostate-specific 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

We apply the CRA Toolbox to the ODE model by setting tuning parameters of the procedure as follows: number of samples N equal to 10000, lower boundary and upper boundary of the LHS equal to 0.1 and 10 respectively. We run 100 independent realizations to verify the reliability and stability of the procedure. We choose different output nodes and evaluation functions in order to show results of the CRA Toolbox in a complete and comprehensive way. Specifically, we select as output nodes both variables C2 and C1, which represent cytotoxic T lymphocyte (CTL) in prostate and lymphoid respectively. For C2, we measure all the three evaluation functions offered by the software, i.e. the area under the curve, the maximum and the time of maximum reached by the time behavior of CTL, as shown in Fig. 4. Conversely, for C1, we perform robustness analysis using as evaluation function only the area under the curve. In all cases, we set equal to 1000 the dimension of the upper and lower tail of the evaluation function pdf, in order to guarantee a stable estimation of the conditional parameter pdfs [5].

Fig. 4

Evaluation functions available in the CRA Toolbox. Nominal time behavior of C2 (blue line) and the three evaluation functions measured for robustness analysis

Results of CRA applied to variable C2 to measure the area are shown in Fig. 5. Parameters α and k have MIRI values above 1, thus they have a great impact on the chosen evaluation function. Other parameters influencing the selected output are μ, α and π, having MIRI values between 0.5 and 1.

Fig. 5

Results of the Prostate-specific Pten−/− mouse model: C2 area. Output of the CRA Toolbox when the area under the curve of C2 time behavior is chosen as evaluation function. a Pdf of the area under the curve of C2. b Boxplot of the 100 realizations of the MIRIs. c-d Conditional pdfs of the parameters with a MIRI value above 1

Similar results are obtained for the maximum value of C2, as shown in Fig. 6. Parameters with the highest MIRI values (∼1) are k and α as before, while parameters r, α and μ have all MIRIs around 0.5.

Fig. 6

Results of the Prostate-specific Pten−/− mouse model: C2 maximum value. Output of the CRA Toolbox when the maximum value of C2 time behavior is chosen as evaluation function. a Pdf of the maximum value of C2. b Boxplot of the 100 realizations of the MIRIs. c-d Conditional pdfs of the parameters with a MIRI value above 1

Figure 7 shows results obtained for the time of maximum of variable C2. In this case, MIRIs have values lower if compared to the previous examples. The most influential parameters are α, k and μ, with MIRI values around 0.6.

Fig. 7

Results of the Prostate-specific Pten−/− mouse model: C2 time of maximum. Output of the CRA Toolbox when the time of maximum of C2 time behavior is chosen as evaluation function. a Pdf of the time of maximum of C2. b Boxplot of the 100 realizations of the MIRIs. c-d-e Conditional pdfs of the parameters with a MIRI value above 0.5

As regards variable C1, results of the measured area are shown in Fig. 8. Two are the parameters with MIRIs above 1: r and α. All the remaining parameters have low values, except for α and with values around 0.5.

Fig. 8

Results of the Prostate-specific Pten−/− mouse model: C1 area. Output of the CRA Toolbox when the area under the curve of C1 time behavior is chosen as evaluation function. a Pdf of the area of C1. b Boxplot of the 100 realizations of the MIRIs. c-d Conditional pdfs of the parameters with a MIRI value above 1

Table 5 contains the time necessary to perform all the simulations.

Table 5

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
Prostate-specific Pten−/− mouse model		Time of Maximum	1706
	C1	Area	1708
EGFR-IGF1R pathway in lung cancer	ERK	Area	102

The times are expressed in seconds and they refer to a single realization of the CRA

Pulse generator network

We test the CRA Toolbox on a small toy system belonging to the field of Synthetic Biology because synthetic models are one of the best examples of the importance of theoretical modeling in the biological reality [14]. The pulse generator network consists of three nodes, representing three genes, aimed at producing a transient output response to a persistent input stimulus [14]. Figure 9 shows the interaction schema of the model. Node S1 is the input step signal that activates both R2 and Y. R2 is the so called repressor because it acts as a deactivator of the product Y. The corresponding ODE model has two state variables, eight kinetic parameters and one input signal. The following ODEs are written using a Hill function for the activation and repression function [5, 15]:

Fig. 9

Pulse generator network. Three nodes network of the pulse generator model: solid arrows stand for activation and the T-dashed line represents a deactivation function

Nominal values for the parameters in Eq. 6 are k1=5 nM/min, k12=20 nM/min, λ2=0.01 nM/min, λ=0.04 nM/min, K1=1 nM, K2=100 nM and n1= n2=3. We run the CRA Toolbox setting the tuning parameters as reported in [5]. More in detail, we set the number of realizations, the lower boundary, the upper boundary and the number of samples equal to 100, 0.1, 10 and 10000 respectively. Parameters n1 and n2 remain fixed to their nominal values since they are not included in robustness analysis. We select as reference node the observable Y and we progressively perform the CRA using all the three evaluation functions provided by the tool, i.e. the area under the curve, the maximum value and the time of maximum. Figure 10 shows the pdf of the area and the corresponding boxplot of MIRIs when choosing 1000 as the dimension of both the lower and higher tails. The other subfigures of Fig. 10 show the conditional pdfs of each parameter that are used for the calculus of MIRIs. Figure 11 shows the results of robustness analysis when the maximum of output variable Y is selected as evaluation function. Figure 12 reports the output generated by the CRA Toolbox about the time of maximum of variable Y. From Figs. 10b, 11b and 12b, it is clear how parameter λ has a great impact on the output behavior of Y, since it influences all the three evaluation functions with an high value of the corresponding MIRI. Moreover, parameter k12 is the one with highest value of the MIRI in Fig. 11b and thus it especially influences the maximum value reached by the Y node. Finally, parameter K2 is the most relevant when dealing with the time of maximum of variable Y because it has the highest MIRI value in Fig. 12b. Table 5 reports the time required to complete all the simulations described in this section.

Fig. 10

Pulse generator network results: Y area. Output of the CRA Toolbox selecting the area of Y as evaluation function. a Pdf of the area of Y. b Boxplot of the MIRIs for the 100 realizations. c-d-e-f-g-h Conditional pdfs of parameters k1, K1, K2, k12, λ2 and λ respectively

Fig. 11

Pulse generator network results: maximum of Y. Output of the CRA Toolbox selecting the maximum of Y as evaluation function. a Pdf of the maximum of Y. b Boxplot of the MIRIs for the 100 realizations. c-d-e-f-g-h Conditional pdfs of parameters k1, K1, K2, k12, λ2 and λ respectively

Fig. 12

Pulse generator network results: time of maximum of Y. Output of the CRA Toolbox selecting the time of maximum of Y as evaluation function. a Pdf of the time of maximum of Y. b Boxplot of the MIRIs for the 100 realizations. c-d-e-f-g-h Conditional pdfs of parameters k1, K1, K2, k12, λ2 and λ respectively

EGFR-IGF1R pathway in lung cancer

The last example of application of the CRA Toolbox is the EGFR-IGF1R pathway, which is one of the most relevant interaction network for the study of cancer pathogenesis and progression in Non-Small Cell Lung Cancer (NSCLC). Figure 13 depicts the pathway. For detailed information about the biological importance of this network and the role of the different nodes see [16]. The corresponding mathematical ODE model is published in [17] (id:MODEL1209230000). The dynamical ODE model is composed of ten equations (Eq. 7), containing two types of kinetic laws: the law of mass action and the Michaelis-Menten kinetics.

Fig. 13

Pathway of the EGFR-IGF1R model in lung cancer [16]

Table 3 explicits names and initial concentrations of the ten state variables, of the three input signals and of the eight total protein parameters. Table 4 lists all the kinetic model parameters and their nominal values. We run the CRA Toolbox according to the guidelines in [5]: the number of realizations, the lower boundary, the upper boundary and the number of samples are set equal to 100, 0.1, 10 and 10000 respectively. We select as reference node the active form of ERK protein, since it is one of the best indicator of the proliferation attitude of lung cancer. More in detail, we perform robustness analysis of the model choosing as evaluation function the area under the curve of the time behavior of ERK. Figure 14a shows the pdf of the evaluation function and Fig. 14b the boxplot for each one of the 39 MIRI parameters. This figure clearly shows that four parameters have a MIRI higher than 0.6 and thus they significantly influence the dynamical behavior of the output node. More in detail, these parameters are: k, KM, k and KM and Fig. 14c, d, e and f show the corresponding conditional pdfs.

Table 3

List of species included in the EGFR-IGF1R 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
x3T	DSOS	120000
x4T	Ras	120000
x5T	Raf	120000
x6T	MEK	600000
x7T	Erk	600000
x8T	p90Rsk	120000
x9T	PIK3	120000
x10T	Akt	120000

Table 4

List of the kinetic parameters of the EGFR-IGF1R 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 p90Rsk-Erk	0.0213697
p 5	K M p90Rsk-Erk	763523
p 6	k SOS-EGFR	694.731
p 7	KM SOS-EGFR	6086070
p 8	k Ras-SOS	32.344
p 9	KM Ras-SOS	35954.3
p 10	k Erk-MEK	9.85367
p 11	KM Erk-MEK	1007340
p 12	k DSOS-p90Rsk	161197
p 13	KM DSOS-p90Rsk	896896
p 14	k SOS-IGFR	500
p 15	KM SOS-IGFR	1000000
p 16	k PIK3-IGF1R	10.6737
p 17	KM PIK3-IGF1R	184912
p 18	k PIK3-EGFR	10.6737
p 19	KM PIK3-EGFR	184912
p 20	k Akt-PIK3	0.0566279
p 21	KM Akt-PIK3	653951
p 22	kd Akt	0.005
p 23	k Erk-PP2A	9.85367
p 24	KM Erk-PP2A	1007340
p 25	k PIK3-Ras	0.0771067
p 26	KM PIK3-Ras	272056
p 27	k Raf-Ras	0.884096
p 28	KM Raf-Ras	62464.6
p 29	k Raf-MEK	185.759
p 30	KM Raf-MEK	4768350
p 31	k Raf-Akt	15.1212
p 32	KM Raf-Akt	119355
p 33	k Ras-RasGab	1509.36
p 34	KM Ras-RasGab	1432410
p 35	k MEK-PP2A	2.83243
p 36	KM MEK-PP2A	518753
p 37	k Raf-RafPP	0.126329
p 38	KM Raf-RafPP	1061.71
p 39	kd p90Rsk	0.005

Fig. 14

EGFR-IGF1R pathway results: ERK area. Output of the CRA Toolbox selecting as reference node ERK. a Pdf of the area under the curve of ERK. b Boxplot of the 100 realizations of the MIRIs. c-d-e-f Conditional pdfs of the parameters with a MIRI value above 0.6

Table 5 reports the time required to complete the robustness analysis through the CRA algorithm.

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 variance-based 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 Object-Oriented 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 i7-4700HQ CPU, 2.40GHz 8, 16-GB 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 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 times. Moreover, the choice of N influences the estimation of the evaluation function pdf, whose cost increases with N. 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/CRA-Matlab

Operating system(s): Platform independent

Programming language: MATLAB

Other requirements: MATLAB

License: MATLAB

Any restrictions to use by non-academics: License needed for MATLAB

This .pdf file is a detailed description of all the classes and methods implemented in the tool and it is a useful guide for users to understand how to run and use the software package. (PDF 62 kb)

This .png file is a screenshot of the first part of the GUI of the CRA Toolbox. (PNG 47 kb)

This .png file is a screenshot of the second part of the GUI of the CRA Toolbox. (PNG 60 kb)

This .pdf file is the user guide of the CRA Toolbox. (PDF 327 kb)