A versatile mathematical workflow to explore how Cancer Stem Cell fate influences tumor progression
 Chiara Fornari^{1, 2},
 Gianfranco Balbo^{1, 3},
 Sami M Halawani^{3},
 Omar BaRukab^{3},
 Ab Rahman Ahmad^{3},
 Raffaele A Calogero^{2},
 Francesca Cordero†^{1} and
 Marco Beccuti†^{1}Email author
https://doi.org/10.1186/175205099S3S1
© Fornari et al.; licensee BioMed Central Ltd. 2015
Published: 1 June 2015
Abstract
Background
Nowadays multidisciplinary approaches combining mathematical models with experimental assays are becoming relevant for the study of biological systems. Indeed, in cancer research multidisciplinary approaches are successfully used to understand the crucial aspects implicated in tumor growth. In particular, the Cancer Stem Cell (CSC) biology represents an area particularly suited to be studied through multidisciplinary approaches, and modeling has significantly contributed to pinpoint the crucial aspects implicated in this theory.
More generally, to acquire new insights on a biological system it is necessary to have an accurate description of the phenomenon, such that making accurate predictions on its future behaviors becomes more likely. In this context, the identification of the parameters influencing model dynamics can be advantageous to increase model accuracy and to provide hints in designing wet experiments. Different techniques, ranging from statistical methods to analytical studies, have been developed. Their applications depend on casespecific aspects, such as the availability and quality of experimental data, and the dimension of the parameter space.
Results
The study of a new model on the CSCbased tumor progression has been the motivation to design a new workflow that helps to characterize possible system dynamics and to identify those parameters influencing such behaviors. In detail, we extended our recent model on CSCdynamics creating a new system capable of describing tumor growth during the different stages of cancer progression. Indeed, tumor cells appear to progress through lineage stages like those of normal tissues, being their division autoregulated by internal feedback mechanisms. These new features have introduced some nonlinearities in the model, making it more difficult to be studied by solely analytical techniques. Our new workflow, based on statistical methods, was used to identify the parameters which influence the tumor growth. The effectiveness of the presented workflow was firstly verified on two well known models and then applied to investigate our extended CSC model.
Conclusions
We propose a new workflow to study in a practical and informative way complex systems, allowing an easy identification, interpretation, and visualization of the key model parameters. Our methodology is useful to investigate possible model behaviors and to establish factors driving model dynamics.
Analyzing our new CSC model guided by the proposed workflow, we found that the deregulation of CSC asymmetric proliferation contributes to cancer initiation, in accordance with several experimental evidences. Specifically, model results indicated that the probability of CSC symmetric proliferation is responsible of a switchinglike behavior which discriminates between tumorigenesis and unsustainable tumor growth.
Keywords
Background
The use of mathematical models to investigate biological systems is becoming progressively crucial to better understand their complex behaviors [1]. One of the remarkable results obtained through mathematical models is the suggestion of when and how cell fate determines the outcome of the phenomenon under study.
For instance, in cancers it is crucial to characterize when and how cells provide a balance between stem cells and daughtercell lineages, since cell heterogeneity significantly contributes to tumor progression and maintenance. In this context, the acquisition of the tumor microenvironment concept has contributed to define a new generation of hallmarks of cancer [2], which extended the original ones [3] incrementing the complexity of the tumor biology. Specifically, different cell types make up tumor microenvironment which, on turn, preserves cell heterogeneity and provides regulations of cell individuality and cell collective functions. This continuous and finely tuned interplay can promote cancer outbreak, sustain tumor development and invasion, and provide niches for Cancer Stem Cells (CSCs) [4]. CSCs are defined as cells that possess the capacity to both selfrenew and to generate the heterogeneous lineage of cancer cells comprising the tumor [5]. Moreover, they are considered cancerpromoters thanks to their ability of developing new tumors upon inoculating them into host mice [6]. Many evidences point out that CSCs drive tumor growth and evolution of several human cancers, such as lung [7], brain [8], colon cancers [9], etc.
In this paper we focus on CSCbased tumors, which have been largely investigated through approaches combining wetlab experiments and mathematical techniques, as demonstrated by several papers [10–13]. CSCtumors are hierarchically structured and characterized by different subpopulations of cells: CSCs, Progenitor Cells (PCs), and Terminal Cells (TCs). This heterogeneity influences both cancer progression and response to treatments, making fundamental the full understanding of the mechanisms underlying the CSC hierarchy [14]. In particular, the alternation of symmetrical vs. asymmetrical CSC division and the way in which feedback mechanisms  induced by microenvironmental changes  affect the tumor growth have been investigated in several papers [15–17]. However, it is not clear how the balance between the CSC asymmetric and symmetric division rates is maintained in order to preserve a constant level of CSCs in tumors and, at the same time, generate more differentiated cells.
To investigate this issue, we expanded our linear Ordinary Differential Equation (ODE) model on the initial phase of CSCcancer growth [18] describing each stage of tumor progression. Our extended model is composed of CSCs, PCs, TCs and Dead Cells (DCs) subpopulations and it accounts for the tumor microenvironment effects, which can be modeled as mechanisms of auto growth limitation expressed by feedback controls on cell division. In our model tumor cells are assumed to progress through lineage stages like those of normal evolution, and a bounded cell division is introduced. These new features have introduced some nonlinearities in the model, making it more difficult to be studied by solely analytical techniques. Thus, we have proposed a new analysis workflow that helps to characterize system behaviors and to identify those parameters influencing such behaviors. The techniques that we have used range from statistical methods (as sensitivity analysis) to analytical studies (as bifurcation analysis). Their applications depend on casespecific aspects, such as the availability and quality of experimental data, and the dimension of the parameter space as well.
Therefore, the aim of this paper is twofold: (i) to extend our previous model on CSCtumor growth  including more complex dynamics in the ODE system  and (ii) to provide a workflow for the analysis of this class of models, using state of the art methodologies for sensitivity analysis [19].
Although the proposed methodology has been designed to study our tumorgrowth model, it is not cancer specific. It could be applied to study dynamical systems in general, and to help in assessing parameter influences on ODE model dynamics. Moreover, the identification of key model parameters could also provide hints to design new experiments that might enhance the knowledge of the phenomenon under study.
Materials and methods
Population model
In Fornari et al. [18] we presented a system of ODEs focused on the description of the initial phase of CSCtumor progression. The proposed model was linear and, consequently, its predictions relative to cell homeostasis were sensitive to small differences in parameter values. Indeed, to achieve homeostasis in any physical or biological system, feedback mechanisms are necessary to maintain stability in the face of infinitesimal parameter changes.
Here, we present an extension of our model, which accounts for stable cell homeostasis and considers cell subpopulation dynamics during the cancer growth. This has required to introduce: (i) a new cell subpopulation called Dead Cells (DCs) and (ii) a feedback control of cell number.
where N_{ CSC } , N_{ PC } , N_{ TC } and N_{ DC } are the total number of CSCs, PCs, TCs and DCs, respectively. Notice that we have modeled the PC subpopulation as made by two different levels: PC_{1} and PC_{2}, i.e. ${N}_{PC}={N}_{P{C}_{1}}+{N}_{P{C}_{2}}$. Possible cell behaviors are regulated by specific rates expressed through the following phenotypic parameters: P_{ sy } for the probability of symmetric CSC division; ω_{ CSC } , ω_{ PC } for CSC and PC_{1} proliferation; η_{1}, η_{2} and η_{3} for CSC, PC_{1} and PC_{2} differentiation; γ_{ PC } for PC_{1} dedifferentiation; δ_{1}, δ_{2} and δ_{3} for CSC, PC, and TC death, while δ_{4} is used to determine the DC lysis.
where h_{ CSC } and h_{ PC } correspond to the feedback intensities. Namely, CSC and PC_{1} proliferation rates now depend on the TC number, which is in agreement with the knowledge that the growth and progression of cancer cells depend not only on their intrinsic malignant potential, but also on a mutual and continuous dialog among them and the tumor microenvironment [4]. Indeed, the growth conditions of CSCs are influenced by blood supplies, growth factor and, in particular, local cell types. Considering several experimental observations [16, 15] we parameterize this knowledge by specifying with equations (2) a bounded cell division regulated by the density of TC subpopulation.
Finally, note that the introduction of the auto growth mechanism makes the ODE system (1) non linear.
Methodology
In this section we describe the methodology developed to study our complex model showing how it is sufficiently general and powerful to be effectively applied when the presence of uncertainties in experimental data makes the analysis of biological models extremely complicated.
 1input/output characterization

1.1 sampling of model inputs through the Latin Hypercube Sampling (LHS) technique;

1.2 investigation of possible model behaviors deriving from different parameter samples, i.e. creation of model outputs;

 2key parameter identification

2.1 evaluation of Partial Rank Correlation Coefficients (PRCCs) between model parameters (inputs) and model behaviors (outputs), at different time points;

2.2 identification of model parameters and time points which need to be carefully studied;

 3key parameter analysis

3.1 colored visualization of model outputs (i.e. variable values over time) with respect to the values of "key" parameters identified during the previous step;

3.2 creation of scatter plots representing model outputs (at designated time points) versus parameter values, colored in accordance with key parameter values;

3.3 accurate characterization of the role of the selected parameters through analytical techniques (as bifurcation analysis).

The application of the LHS method and PRCC analysis may be performed starting from Matlab functions described in [19], enhanced to extend their analysis capabilities of specific parameters by plotting graphs colored according to key parameter values. Bifurcation analysis may be performed using the graphical Matlab package MatCont.
Input/output characterization
Starting from these input values, each row of the LHS matrix is used as an input to numerically integrate the system over the time interval T = [t_{ ini }, t_{ end }], and thus producing n timedependent model solutions. The time points are selected to homogeneously cover the time interval of interest. If specific time subintervals are known to be relevant, they can be investigated increasing the number of time points in these subintervals of T . Model outputs (matrix Y) are then collected for each experiment considering the different time points of T , such that model temporal behaviors (i.e. model traces) can be derived. Specifically, a plot providing a graphical representation of these timedependent traces is produced for each parameter combination, as sketched in Figure 2, panel A.
Key parameter identification
To identify critical inputs and to quantify how they impact model outcomes, Partial Rank Correlation Coefficients (PRCCs) between model parameters (LHS matrix) and model outputs (Y matrix) are evaluated on the interval T . Indeed, PRCC measures monotonic relationships between outputs and inputs and it provides a measure of monotonicity after the removal of the linear effects of all but one variable [26]. Considering multiple time points, PRCCs between model variables and model parameters are evaluated and plotted, as shown in Figure 2, panel B. PRCC analysis is done in an exploratory way to identify any significant relationships throughout the entire time course, and to point out whether correlations occur either over the entire time interval, or at specific time points. Since even small correlations may be significant, statistical tests assessing if the PRCCs are different from zero are also performed [19]. Details on the indexes that can be used are provided in [27].
PRCC analysis and corresponding significance tests are then used to identify key model parameters and to select time points which need an indepth study. Specifically, among all model parameters, only those with high and significant PRCCs for all variables are further investigated, and the analysis starts from those parameters having the highest PRCCs values. More precisely, PRCC values close to 1 (1) identify positive (negative) monotone relationships between inputs and outputs, by definition. In addition, significance tests allow to discover those correlations that are important, despite having relatively small PRCC values, namely those correlations that have significant pvalues. In this phase of the study it is crucial the role of the analyst who, considering the PRCCs and pvalues, decides on which parameters focus the study, and whereas focus it on the whole time interval or on specific time points. A general criterion to guide this process is difficult to be defined, since it strongly depends on the peculiarities of the model under investigation.
Key parameter analysis
After the identification of the relevant model parameters, a detailed analysis of their effects on model outputs is performed as concluding step of our methodology. For the sake of simplicity we will describe this process considering only one relevant parameter, that we will call p_{1}. The proposed approach can however be easily extended to a set of relevant parameters, {p_{1}, ..., p_{ k } }, reproducing such analysis k times and then crossanalyzing the results. An example of crossanalysis is provided in our second case study (i.e. the apoptosis model), where two relevant parameters have been identified and their crossanalysis have been necessary to explain model bistability.
The variation range of parameter p_{1} is divided into r subintervals, such that the p_{1}'s possible values are partitioned into r levels. The choice of r depends on several factors, such as the variation range of the parameter and its expected qualitative behavior. Similarly to the adaptive numerical methods, an adaptive partition could be defined with a varying r, such that the most a parameter is critical in a region, the smallest r is set to cover that region. Then, as shown in Figure 2  panel C1, a specific color is assigned to each of these levels and model traces are colored accordingly to this classification. Specifically, each model trace is colored in accordance with the p_{1} value used for its computations, i.e. each trace has the same color of the subinterval whom its p_{1} value belongs to; see Figure 2, panel C2. This new representation turns out to be very effective in showing how particular model behaviors, such as switches, bistabilities, etc., are related to the p_{1}'s variation. This visualization of model behaviors can thus be considered as a preliminary, but effective, analysis of the role that p_{1} plays in the global model dynamics.
The subsequent step consists in expanding the analysis at selected interesting time points, i.e. at those points which are known to be crucial for the problem under investigation, or which have been previously detected as relevant in PRCC analysis. Scatter plots of model outputs versus p_{1} values (or other relevant parameters, or combinations of them) are evaluated and colored in accordance with the parameter variation range. As depicted in Figure 2, panel C3, colored scatter plots enable a simple and direct visual detection of correlations between model inputs and outputs, emphasizing the role of p_{1}.
Let us note that when few parameters are identified as relevant for model dynamics, the previous analysis is performed for each of them trying to discover multiple dependencies.
Finally, after having identified in a qualitative manner the key model parameters, an analytic study is performed to obtain an accurate characterization of the actual role of these parameters using more sophisticated mathematical techniques, such as bifurcation analysis.
Results
Before using our methodology to explore the dynamics of the extended CSCtumor model (1), we verified the effectiveness of this approach on two wellknown and experimentally validated models by Tyson et al [28]. These models  that describe the oestrogen signalling network in breast epithelial cells  have been widely studied by Tyson and coworkers, exhibit bistable switches, and are among the reference models in mathematical biology.
Cell cycle model
In [28], cell cycle is modeled describing the interactions among a set of key proteins which control the G1toS phase transition in mammalian cells, namely RetinoBlastoma protein (RB), Cyclin D (CycD), the E2F family of transcription factors (E2F), and Cyclin E (CycE). The cell choice between quiescience and proliferation is controlled by a bistable switch characterized by an OFF state (quiescient cell arrested in G1 phase) and an ON state (proliferating cell progressing through S, G2 and M phase). Analyzing this model, the authors characterized the role of each protein in maintaining alternation between the two stable steady states.
Input/output characterization
Key parameter identification
Key parameter analysis
Color codes.
Model  parameter  I_{1} (black)  I_{2} (blue)  I_{3} (red)  I_{4} (green) 

cell cycle  serum  [0.01, 1.50]  (1.50, 3.00]  (3.00, 4.49]  (4.49, 5.99] 
apoptosis  stress  [0.01, 0.25]  (0.25, 0.50]  (0.50, 0.75]  (0.75, 0.99] 
apoptosis  BCL2T  [56.00, 68.00]  (68.00, 80.00]  (80.00, 91.99]  (91.99, 103.99] 
population  P _{ sy }  [0.00, 0.25]  (0.25, 0.50]  (0.50, 0.75]  (0.75, 1.00] 
Despite the use of completely different approaches, the results obtained from this first case study reproduced exactly those presented by Tyson et al. in [28]. This agreement is a preliminary validation of our workflow and thus supports the usage of our methodology.
Apoptosis model
From the same paper [28], we selected a second model to validate our methodology. This second case study concerns apoptosis in mammalian cells and it is modeled by Tyson et al. through a bistable system. The irrevocable commitment to apotosis reaches a oneway decision point in which prodeath and prosurvival signals are processed determining cell fate. In [28], the authors focused on interactions among proteins BAXm, BCL2, and BH3, which are responsible of the bistable and irreversible switch governing apoptosis. Specifically, the switch is OFF or ON depending on the balance among BAX, BCL2 (brake) and BH3 (accelerator ): in the OFF state BAX is inactivated by binding to BCL2, while in the ON state BAX is active since BCL2 is displaced due to BH3 accumulation.
Input/output characterization
Key parameter identification & analysis
Contrary to the cell cycle model where we found only one relevant parameter, PRCC analysis on apoptosis model revealed two putative key parameters, namely stress and BCL2T (total BCL2 concentration). As reported in Additional File 9, they both have high PRCC values (1 and 0.8, respectively) in the whole time interval [0, 500]. Significance tests resulted in pvalues <0.01, so that we expected them both to be related with the bistable apoptosis switch. We defined a colorrange made of r = 4 variations for both stress and BCL2T (see Table 2 for values), and model traces were then colored and analyzed twice (i.e. one for each parameter). Specifically, as shown in Additional File 8, traces became clustered by colors and color sequences followed orders describing the increase/decrease of the two parameters. By means of this approach we were able to qualitatively detect the monotonic relationships between stress and BH3, stress and BCL2, BCL2T and BH3. We also pointed out how model bistability depends on variations in parameter values, being high/low states of BAXmT and BH3 characterized by specific parameter concentrations (i.e. colors). To further characterize these dependencies it was necessary to crossanalyze the correlations among variables (BAXmT, BH3) and parameters (stress, BCL2T), since they were not sufficient to explain model bistability when considered separately. In detail, we analyzed scatter plots of variables versus stress variations at the end time of our numerical experiments (i.e. t = 500), and we colored the resulting plots twice: one applying stress colorcode, the other one focusing on BCL2T variations. As shown by Figure 7 (panels A and B), BAXmT can be found in two concentrations (high and low), depending on the amount of both stress and BCL2T. Indeed, when cells are in extremely high or low stress conditions, BAXmT concentration is mainly in one configuration; see the points in the areas at the left or right margins of the plots which are mainly concentrated below or over 60. Specifically, as reported in Figure 7, panel A, very low stress values induce low concentrations of BAXmT (black points below 60), while highly stressed cells have high BAXmT concentrations (green points above 60). Instead, for intermediate values of stress (blue and red points), BAXmT has two stable steady states: a low and a high one (below and over 60, respectively). As shown in Figure 7, panel B, these configurations are mainly characterized by the concentration of BCL2T. When BCL2T is high, BAXmT is in the low configuration no matter how much stress is affecting the system, as resulted from the the green points below 60 in panel B.
Similarly, BH3 may evolve into two states: (i) a zero one in which the protein is almost absent, and (ii) a positive one, in which the protein concentration increases as stress does (see Figure 7, panels C and D). Let us note that, while the positive state is sensitive to stress variation, the zero one is independent of the same fluctuation. Indeed, as it is evident from Figure 7 panel C, BH3 may end in the zeroconfiguration for each stress concentration. In fact, what really defines this zerostate is the high concentration of BCL2T, as evident from the green points in Figure 7, panel D.
Summarizing, following our methodology we were able to detect bistabilities in the apoptosis model and to characterize different model evolutions in terms of parameter changes, reaching conclusions similar to those discussed in [28].
Application to our case study
The aim of the previous part of the study was to test (validate) the proposed methodology on wellknown case studies. Having gained confidence in the approach, we applied the method to our population model (1) to explore the phenotypic characteristics inherent to CSCs on tumor growth. Indeed, system (1) is governed by a number of cellular behaviors which are difficult to measure with biological experiments due to their high variability. The proposed approach helps to understand better how the system functions and which phenotype parameters mostly influence it.
Model input/output characterization
Following the indications of the first phase of our methodology, parameter values (model inputs) were generated with the LHS technique and timedependent subpopulation behaviors (model outputs) were evaluated for each sampled parameter combination.
As described in Materials and Methods, the LHS technique needs a set of baseline parameter values to start from. These rates and the initial cell concentrations were tuned starting from values found in the literature [18, 29, 30] and provided by experimental evidences. Moreover, initial conditions were fixed to reproduce specific subpopulation proportions: TCs were set as the largest subpopulation  as they represent the main part of the tumor mass  and CSCs as the smallest one [18]. Parameters were then retrieved by tuning system (1) to reproduce the growth trend of tumor mass observed in BALB/c mice after a subcutaneous injection of 10^{5} cancer cells [31]. Baseline parameter values and model initial conditions are reported in Supplementary Materials as Additional File 10.
Key parameter identification
Correlations between model inputs and model outputs were evaluated over the whole time interval [0, 200] to assess which phenotype parameters mostly influence subpopulation dynamics in terms of monotonic relationships.
We proceeded with the study firstly investigating the role of P_{ sy } and, then, exploring those of ω_{ CSC } and η_{1}.
Key parameter analysis
We firstly focused our study on P_{ sy } , fixing r = 4; see Table 2 for the partition of P_{ sy } values. Model traces were colored in accordance with their P_{ sy } values to detect how P_{ sy } variation affects model dynamics. Output data were also expressed using a logarithmic scale to reduce their wide range to a more manageable size, thus allowing to pinpoint some interesting behaviors that have been subsequently studied in detail.
Summarizing, these results suggested a critical CSC symmetric proliferation value at which tumor stabilization occurs. Moreover, we identified an intermediate region in which other phenotype parameters involved in CSC variation cause tumor growth and maintenance.
Conclusions
In this paper we have proposed a new workflow that helps to characterize ODE system behaviors and to identify those parameters which mostly influence such behaviors. Our workflow uses state of the art methodologies for sensitivity analysis, readapted to allow an easy identification, interpretation, and visualization of key model parameters.
When a limited amount of data is available, proper numerical estimates of model parameters are difficult to obtain. Estimated values must therefore be regarded as preliminary information, and alternative strategies must be identified to assess the quality of the results produced by the models. In these cases, our methodology can be followed to investigate the nature of the relationships between input parameters and output values. In particular, when parameters are difficult to be experimentally measured and may thus be affected by large variations, our methodology allows to study model behavior as parameters are varied with statistical techniques. By varying experimental settings, such as the distribution used, the parameter variation ranges, and their partition, our methodology allows to deeply investigate all system behaviors and to identify which parameters are relevant for explaining interesting output dynamics with a limited computational effort.
Having identified a restricted number of parameters which are relevant for the model, analytical approaches can then be applied within this simplified context. Our workflow is a practical and informative tool to approach the study of complex systems such as the biological models where metabolic pathways are described or where detailed kinetics need to be accounted for.
The effectiveness of this methodology has been verified on two wellknown models, whose results published in the literature have been accurately reproduced using our approach. After this preliminary validation, our methodology has been applied to a CSCtumor model which extends a previous representation of the same phenomenon that we have published in [18] and that aims to figure out which are the phenotypic parameters that drive cancer growth.
Our CSCtumor based model describes all phases of cancer progression and accounts for the negative feedback which TC subpopulation has on the proliferation rates of CSCs and PC1s. With the introduction of this feedback, we have been able to model cell auto growth limitation which has been demonstrated not only on stem cells during organogenesis [32], but also in cancer cells during tumor growth [33]. Applying our methodology to this model we found that the probability of CSC symmetric proliferation is responsible for a switchinglike behavior. Specifically, P_{ sy } discriminates between two possible scenarios: tumorigenesis and unsustainable tumor growth. More precisely, if CSC symmetric proliferation probability has low values, then the system falls in the nongrowth scenario. Otherwise, for high P_{ sy } values, the only possibility is the tumor growth, and no other parameters considered individually are able to characterize the switch. Notice that the condition identifying the non growth scenario is a relaxed condition since the choice r = 4 did not allow us to unequivocally discriminate the threshold among tumorigenesis and unsustainable tumor growth. However, this initial analysis revealed the existence of such a threshold, that we had afterwards investigated in depth by means of bifurcation analysis. Moreover, we found that there is a transition zone in which it is necessary to consider together CSC symmetric proliferation, CSC proliferation rate, and CSC differentiation rate, in order to precisely characterize the tumor evolution. This supports the notion that CSC phenotypic plasticity is able to lead to functionally distinct cancer subpopulations that support and modulate the overall tumor growth and maintenance [2]. Moreover, our finding has been supported by experimental evidences suggesting that the deregulation of asymmetric proliferation contributes to cancer initiation [34, 35].
All dynamics considered in our phenomenological mathematical model are related to interactions among cell populations and are based on cancer stem cell fate, including cell proliferation, differentiation and death. These mechanisms are expressed through the phenotypic parameters of the model, which provides a global description of the tumor growth. However, a deeper characterization of this phenomenon might consider also: (i) the intrinsic noise present in biological data, and (i) the cellular processes within cancer stem cells which control cell fate, and which are tagged as hallmark of cancers. Examples of these processes are the unfolded protein response [36] and the autophagy [37], which are stress response phenomena and which modulate tumor microenvironment, leading to metabolic reprogramming and changes in cancer stem cells fate.
A more detailed investigation could be conducted  as a future work  integrating in our model both the noise and those internal cellular mechanisms which control cancer stem cell fate, such that it might be possible to better characterize which are the microenvironment stimuli that mostly influence symmetric proliferation decision. Results presented in this paper could hence be used to facilitate and improve this integration. Indeed, following the idea presented in our recent paper on multilevel modeling [38], our population model (1) could be corroborated by an additional level focusing on cellular internal dynamics.
Notes
Declarations
Declarations
Publication of the article was funded by grants from the Epigenomics Flagship Project EPIGEN, MIUR CNR, European Seventh framework program, Health.2012.1.21, NGSPTL grant n. 306242, and by project grant Nr. 10151432/HICI from the King Abdulaziz University of Saudi Arabia.
This article has been published as part of BMC Systems Biology Volume 9 Supplement 3, 2015: Proceedings of the Italian Society of Bioinformatics (BITS): Annual Meeting 2014: Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/9/S3.
Authors’ Affiliations
References
 Anderson A, Quaranta V: Integrative matheamtical oncology. Nat Reviews Cancer. 2008, 8: 227234. 10.1038/nrc2329.View ArticleGoogle Scholar
 Hanahan D, Weinberg R: Hallmarks of cancer: the next generation. Cell. 2011, 144: 646674. 10.1016/j.cell.2011.02.013.View ArticlePubMedGoogle Scholar
 Hanahan D, Weinberg R: The hallmarks of cancer. Cell. 2000, 100: 5770. 10.1016/S00928674(00)816839.View ArticlePubMedGoogle Scholar
 Conti L, Ruiu R, Barutello G, Macagno M, Bandini S, Cavallo F, Lanzardo S: Microenvironment, oncoantigens, and antitumor vaccination: lessons learned from balbneut mice. Hindawi Publishing Corporation BioMed Research International. 2014Google Scholar
 Frank N, Schatton T, Frank M: The therapeutic promise of the cancer stem cell concept. The journal of clinical investigation. 2010, 120:Google Scholar
 Cho R, Clarke M: Recent advances in cancer stem cells. Curr Opin Genet Dev. 2008, 18:Google Scholar
 Ho MM, Ng AV, Lam S, Hung JY: Side population in human lung cancer cell lines and tumors is enriched with stemlike cancer cells. Cancer Research. 2007, 67:Google Scholar
 Singh SK, Hawkins C, Clarke ID, Squire JA, Bayani J, Hide T, et al: Identification of human brain tumour initiating cells. Nature. 2004, 432:Google Scholar
 O'Brien CA, Pollett A, Gallinger S, Dick JE: A human colon cancer cell capable of initiating tumour growth in immunodeficient mice. Nature. 2007, 445Google Scholar
 MarciniakCzochra A, Stiehl T, Ho A, Jager W, Wagner W: Modeling of asymmetric cell division in hematopoietic stem cells  regulation of selfrenewals is essential for efficient repopulation. Stem Cells and development. 2009, 18:Google Scholar
 Mirams G, Fletcher A, Maini P, Byrne H: A theoretical investigation of the effect of proliferation and adhesion on monoclonal conversion in the colonic crypt. Journal of Theoretical Biology. 2013Google Scholar
 Dingli D, Traulsen A, Michor F: (a)symmetric stem cell replication and cancer. PLoS Comp Bio. 2007, 3:Google Scholar
 Gupta P, Fillmore C, Jiang G, Shapira S, Tao K, Kuperwasser C, Lander E: Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells. Cell. 2011, 146: 633644. 10.1016/j.cell.2011.07.026.View ArticlePubMedGoogle Scholar
 Nguyen L, Vanner R, Dirks P, Eeves C: Cancer stem cells: an evolving concept. Nature Review Cancer. 12: 2012Google Scholar
 Liu X, Johnson S, Liu S, Kanojia D, Yue W, Singh U, Wang Q, Wang Q, Nie Q, Chen H: Nonlinear growth kinetics of breast cancer stem cells: Implications for cancer stem cell targeted therapy. Sci Rep. 2013, 3:Google Scholar
 Scott J, Hjelmeland A, Chinnaiyan P, Anderson A, D B: Microenvironmental variables must influence intrinsic phenotypic parameters of cancer stem cells to affect tumourigenicity. PLOS Compultational Biology. 2014, 10:Google Scholar
 Youssefpour H, Li X, A L, Lowengrub J: Multispecies model of cell lineages and feedback control in solid tumors. Journal of theoretical biology. 2012, 304:Google Scholar
 Fornari C, Beccuti M, Lanzardo S, Conti L, Balbo G, Cavallo F, Calogero R, Cordero F: A mathematicalbiological joint effort to investigate the tumorinitiating ability of cancer stem cells. PLOSE One. 2014, 9:Google Scholar
 Marino S, Hogue I, Ray C, Kirschner D: A methodology for performing global uncertainity and sensitivity analysis in systems biology. Journal of Theoretical Biology. 2008, 254: 178196. 10.1016/j.jtbi.2008.04.011.PubMed CentralView ArticlePubMedGoogle Scholar
 Lo W, Chou C, Gokoffski FK, Wan , Lander A, Calof QA, Nie : Feedback regulation in multistage cell lineages. Mathematical Biosciences and Engineering. 2009, 6:Google Scholar
 Murray J: Mathematical Biology: I. An Introduction (Interdisciplinary Applied Mathematics). 2007, 3Google Scholar
 Saltelli A, Chan K, Scott EM: Sensitivity Analysis. 2009Google Scholar
 Helton J, Davis F: Samplingbased Methods for Uncertainty and Sensitivity Analysis. 2000View ArticleGoogle Scholar
 Iman R, Helton J: An investigation of uncertainty and sensitivity analysis techniques for computer models. Risk Anal. 1988, 8: 7190. 10.1111/j.15396924.1988.tb01155.x.View ArticleGoogle Scholar
 Mckay M, Beckman R, Conover W: Comparison of 3 methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. 1979, 21: 239245.Google Scholar
 Helton J, Davis F: Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliab Eng Syst Saf. 2003, 81: 2369. 10.1016/S09518320(03)000589.View ArticleGoogle Scholar
 Anderson T: An Introduction to Multivariate Statistical Analysis, third ed. edn. Wiley Series in Probability and Statistics. 2003, WileyInterscience, Hoboken, NJGoogle Scholar
 Tyson J, Baumann T, Chen C, Verdugo I, Tabassoly Y, Wang Y, Weiner L, Clarke R: Dynamic modelling of oestrogen signalling and cell fate in breast cancer cells. Nature Reviews Cancer. 2010, 11: 523532. 10.1038/nrn2850.View ArticleGoogle Scholar
 Turner C, Kohandel M: Investigating the link between epithelialmesenchymal transition and the cancer stem cell phenotype: A mathematical approach. Journal of Theoretical Biology. 2010, 265: 329335. 10.1016/j.jtbi.2010.05.024.View ArticlePubMedGoogle Scholar
 MolinaPena R, Álvarez M: A simple mathematical model based on the cancer stem cell hypothesis suggests kinetic commonalities in solid tumor growth. Plos One. 2012, 7:Google Scholar
 Conti L, Lanzardo S, Arigoni M, Antonazzo R, Radaelli E, Cantarella D, Calogero R, Cavallo F: The noninflammatory role of high mobility group box 1/tolllike receptor 2 axis in the selfrenewal of mammary cancer stem cells. FASEB Journal. 2013, 10:Google Scholar
 Johnston M, Edwards C, Bodmer W, Maini P, Chapman S: Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proc Natl Acad Sci USA. 2007, 104: 40084013. 10.1073/pnas.0611179104.PubMed CentralView ArticlePubMedGoogle Scholar
 RodriguezBrenes I, Komarova N, Wodarz D: Evolutionary dynamics of feedback escape and the development of stemcelldriven cancers. Proc Natl Acad Sci USA. 2011, 108: 1898318988. 10.1073/pnas.1107621108.PubMed CentralView ArticlePubMedGoogle Scholar
 Sugiarto S: symmetrydefective oligodendrocyte progenitors are glioma precursors. Cancer Cell. 2011, 20: 328340. 10.1016/j.ccr.2011.08.011.PubMed CentralView ArticlePubMedGoogle Scholar
 Ito T: Regulation of myeloid leukaemia by the cellfate determinant musashi. Nature. 2010, 466: 765768. 10.1038/nature09171.PubMed CentralView ArticlePubMedGoogle Scholar
 Erguler K, Pieri M, Deltas C: A mathematical model of the unfolded protein stress response reveals the decision mechanism for recovery, adaptation and apoptosis. BMC Systems Biology. 2013, 7:Google Scholar
 Tavassoly I, Shajahan A, Parmar J, Baumann W, Clarke R, Tyson J: Dynamical modeling of the interaction between autophagy and apoptosis in mammalian cells: a systems pharmacology framework. arXiv preprint. 2013Google Scholar
 Cordero F, Beccuti M, Fornari C, Lanzardo S, Conti L, Federica C, Balbo G, Calogero R: Multilevel model for the investigation of oncoantigendriven vaccination effect. BMC Bioinformatics. 2013, 14:Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 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.