# Proposal of a hybrid approach for tumor progression and tumor-induced angiogenesis

Patricio Cumsille, Aníbal Coronel, Carlos Conca, Cristóbal Quiñinao, Carlos Escudero

► To cite this version:

Patricio Cumsille, Aníbal Coronel, Carlos Conca, Cristóbal Quiñinao, Carlos Escudero. Proposal of a hybrid approach for tumor progression and tumor-induced angiogenesis. Theoretical Biology and Medical Modelling, 2015, 12, pp.13. 10.1186/s12976-015-0009-y. hal-01223150

HAL Id: hal-01223150

https://hal.sorbonne-universite.fr/hal-01223150

Submitted on 3 Nov 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Proposal of a hybrid approach for tumor progression and tumor-induced angiogenesis

Patricio Cumsille<fnref n="1" />,²,³,<fnref n="4" />*, Aníbal Coronel<fnref n="1" />,³, Carlos Conca⁵,², Cristóbal Quiñinao⁶ and Carlos Escudero<fnref n="1" />,<fnref n="4" />,⁷*

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

One of the main challenges in cancer modelling is to improve the knowledge of tumor progression in areas related to tumor growth, tumor-induced angiogenesis and targeted therapies efficacy. For this purpose, incorporate the expertise from applied mathematicians, biologists and physicians is highly desirable. Despite the existence of a very wide range of models, involving many stages in cancer progression, few models have been proposed to take into account all relevant processes in tumor progression, in particular the effect of systemic treatments and angiogenesis. Composite biological experiments, both *in vitro* and *in vivo*, in addition with mathematical modelling can provide a better understanding of these aspects. In this work we proposed that a rational experimental design associated with mathematical modelling could provide new insights into cancer progression. To accomplish this task, we reviewed mathematical models and cancer biology literature, describing in detail the basic principles of mathematical modelling. We also analyze how experimental data regarding tumor cells proliferation and angiogenesis *in vitro* may fit with mathematical modelling in order to reconstruct *in vivo* tumor evolution. Additionally, we explained the mathematical methodology in a comprehensible way in order to facilitate its future use by the scientific community.

**Keywords:** Hybrid approach, Tumor progression, Mathematical modelling, Parameter estimation, Experimental design

## Introduction

Cancer is a large group of diseases that could affect any part of the body, characterized by abnormal cell proliferation, and an increasing migration rate that could derive in invasion and organ spreading, becoming the leading cause of death all over the world [1]. Thus, cancer is a critical societal and scientific problem. Great amounts of human and material resources are yearly spent in developed countries in attempts to understand its root causes and to develop successful prevention and treatment strategies.

Most of human cancers have acquired six basic capabilities [2]: self-sufficiency in growth signals, insensitivity to growth-inhibitory signals, programmed cell death evolution, limitless replication potential, sustained angiogenesis, and tissue invasion, which could cause metastasis. In other words, the defense mechanism preventing each of these acquired capabilities must be thwarted before cells become a malignant and invasive tumor [2].

<footnote>*Correspondence:</footnote>

<footnote>pcumsille@ubiobio.cl;</footnote>

<footnote>escudero@ubiobio.cl</footnote>

<footnote n="1">¹Department of Basic Sciences,
Faculty of Sciences, Universidad del
Bío-Bío, Campus Fernando May, Av.
Andrés Bello s/n, Casilla 447 Chillán,
Chile</footnote>

<footnote n="4">⁴Group of Investigation in Tumor Angiogenesis (GIANT), Chillán, Chile
Full list of author information is
available at the end of the article</footnote>
Hence, it is necessary address new strategies for cancer understanding and treatment. For doing so, establishing correlation between observable phenomena, as well as, observing if a particular intervention produces a significant response could be very useful to state hypotheses postulating which physical processes are involved and how they interact. For example, tumor cell proliferation (observable phenomenon) can be correlated to oxygen concentration (another observable phenomenon). Experimentally speaking, if we would modify the oxygen availability to observe whether cell proliferation is significantly modified, then we could estimate the oxygen threshold over which cancer cells begin to proliferate. Biological experiments needed to test such hypotheses can be time-consuming, expensive and/or impossible with current technology. In these cases, mathematical modelling plays a key role, providing an independent check of the consistency of the hypothesis and can also improve experimental design by identifying which measurements are needed to test a particular theory, and additionally, whether new hypotheses can be established from experimental results (see the review paper by Byrne [3]). In this regard, mathematical modelling is a theoretical description of biological phenomena that may be calibrated by experimental data comparison. Moreover, by changing the parameter values of descriptive equations, the significance and functions of variables representing specific biological features can be easily tested. Then, practically all tumor growth features can be mathematically modelled, reducing the biological modelling complexity and offering a powerful tool to better understand tumor biology, facilitating drug development and also pre-clinical and clinical patient management.

In this work we describe the basic principles of mathematical modelling, as well as, suggest some experimental design in order to obtain relevant biological or clinical data necessary to estimate relevant parameters of mathematical modelling. We also define *hybrid approach* as the feedback capacity between mathematical modelling and biological experimental design which is required, as it has been shown in literature, for better understanding cancer, fitting parameters of a given model with a specific biological scenario, and for obtaining models with predictive capability. This work is not intended to deepen in the previous aspects, but it can serve as a first bridge of communication between applied mathematicians, biologist and physicians or even to be a tool to emphasize research in this field.

## Overview of tumor growth and tumor-induced angiogenesis

From a point of view of materials science, a cancerous tumor is a very complex, multiscaled material. The underlying genetic defects (nano-scale) affect cellular-level (micro-scale) properties, that in turn affect material properties at the tumor level (millimeter-scale). The larger scales may even feedback to the smaller scales. For example, a change in cellular stress may affect the type and occurrence of genetic mutations. Cancer can be studied at any of the three scales as suggested by Zheng et al. [2] and the references therein.

In this work we focus at the millimeter scale, i.e., tumor tissue level. Schematically, three successive stages can be identified in the growth of a solid tumor due to deregulation of cell division: initial *avascular growth*, *angiogenesis* and finally *vascular growth*, which in turn may induce metastatic spread. For a summary of models and main mathematical references discussed in this section, the reader is addressed to Table 1.
<table><caption>**Table 1** The aspects of cancer encompassed by the mathematical models discussed in this work</caption><thead><tr><th>Model</th><th>Avascular growth</th><th>Vascular growth</th><th>Angiogenesis delivery</th><th>Drug</th></tr></thead><tbody><tr><td>[24] first work in modelling avascular tumor growth</td><td></td><td></td><td></td><td></td></tr><tr><td>[32] first work in modelling tumor angiogenesis</td><td></td><td></td><td></td><td></td></tr><tr><td>[2, 10, 34, 35] works focused on tumor angiogenesis and invasion</td><td></td><td></td><td></td><td></td></tr><tr><td>[36–42] works focused in drug effectiveness</td><td></td><td></td><td></td><td></td></tr><tr><td>[31] review of avascular tumor growth</td><td></td><td></td><td></td><td></td></tr><tr><td>[3, 29, 30, 33] reviews of vascular tumor growth</td><td></td><td></td><td></td><td></td></tr></tbody></table>

## Overview of avascular tumor growth

During the avascular growth phase, oxygen and other nutrients are delivered to the tumor cells, and the waste products are removed from the tumor via diffusion from nearby blood vessels; under this condition the tumor cells proliferate rapidly consuming more oxygen than the host cells [4]. However, considering that healthy tissue has approximately 7% oxygen (53 mmHg) tension and the diffusion distance of oxygen in tissue is ~ 100 µm [5], tumor growth is limited in size [6] and maintained during a short period of time. Under these conditions, the neoplastic compartment rapidly exceeds the diffusion distance of oxygen and becomes hypoxic. Thus, oxygen tension in a tumor can range from physiological (7%) to severe hypoxic (< 1% oxygen) or even anoxic [5]. Moreover, the immediate molecular response to low oxygen is the hypoxia-induced factor (HIF) protein stabilization [7, 8], which in turn triggers the expression of target genes involved in hypoxia adaptation such as vascular endothelial growth factor (VEGF). In addition, under severe hypoxia, it is found surrounding areas of necrosis, which is a common characteristic of solid tumors. Nevertheless, quiescent cells generation in the periphery of culture dish was observed *in vitro* during a nutrient-deficient medium and hypoxia exposure in tumor cells [9]. Therefore schematically, in the hypoxic phase, a tumor grows within the limits of its local environment forming three characteristic layers: 1) cells towards the center, deprived of vital nutrients, will die and give rise to a necrotic core; 2) proliferating cells can be found in the outer cell layers; and 3) a layer of quiescent (or hypoxic) cells, which survive without dividing with slow metabolism, is found between the two others layers. In summary, tumor growth is severely restricted and in order to continue its developing needs to find additional nutrient sources.

## Overview of vascular tumor growth

As a result of hypoxic pressure some tumor cells secrete a number of diffusive chemical substances - called *tumor angiogenic factors* (TAF), such as VEGF, into the surrounding tissue, which encourage the body to vascularize the tumor and therefore provide new nutrients. Indeed, transition from the avascular to the vascular state, depends on the tumor ability to induce new blood vessels formation from the surrounding tissue [10]. These blood vessels sprout towards the tumor and then gradually surround and penetrate it, providing an adequate micro-circulation and blood supply. Tumor-induced angiogenesis, the process by which new blood vessels develop from an existing vasculature, through
sprouting, proliferation and fusion of endothelial cells, therefore is a critical step in solid tumor growth.

Under hypoxia, cancer stem cells (CSC), that reside inside the tumor, could differentiate toward endothelial progenitor cells and mature endothelium, which in turn generates new blood vessels inside the tumor [11]. This last process is called neovasculogenesis, and in the past was thought to happen only intrauterine, but nowadays it is known that this event happens also in adulthood. In cancer field, presence of CSC has been associated with tumor recurrence, resistance to chemotherapy, tumor metastasis, and in general with poor clinical prognosis [12, 13]. Tumor-induced angiogenesis is believed to start when a small avascular tumor exceeds a critical diameter (~ 2 mm), above which normal tissue vasculature is no longer able to support its growth. By the time a tumor has grown to a size whereby it can be detected by clinical means, there is a strong likelihood that it has already reached the vascular growth phase (see [10] and the references therein).

Tumor-induced angiogenesis is characterized by a chaotic tumor vessels development associated with both angiogenesis and vasculogenesis. Thus, endothelial cells proliferate and capillaries are rapidly formed allowing tumor growth, but in cancer framework theses processes are deregulated. Indeed, vessel diameter is five times bigger than in normal tissue [14], rarely differentiated into arterioles or venules, with frequently blind endings, and incomplete and abnormal endothelial cell lining [15–17]. All these abnormalities generate an irregular blood flow [5], perpetuate the intermittent low oxygen delivery, increase HIF activity and promote pro-angiogenic signals generation [18–20].

Over the past decade much work has been performed to understand the angiogenic process. For instance, it is well-known that increased density of blood vessels (the so-called “hot spots”), and high VEGF plasma levels are a powerful prognosis tool in many human tumor types [21]. In fact, the capacity to modify this process has been considered as a keystone for cancer treatment, which includes some molecules that reduce the ability of TAF, and in particular of the VEGF, to provide blood supply toward the tumor and thus controlling tumor growth. Indeed, a recent meta-analysis [22] including 24 randomized trials with 8 different types of cancers, in which the synthetic antagonist of VEGF (Bevacizumab®) was used in combination with chemotherapy, shown a statistically significant improvement in the overall survival and progression of free cancer survival in patients who received Bevacizumab® compared with those who did not receive this drug. These beneficial effects were more evident in patients with colon cancer and renal cell carcinoma, but less evident in those who had breast, pancreatic or prostate cancer. This suggests that despite the overall benefit of Bevacizumab®, some patients and types of cancer are more resistant to antiangiogenic therapy. A possible explanation to this phenomenon is the presence of CSC, whose differentiation could occur independently of VEGF [23].

Altogether these last evidences indicate that angiogenesis during tumor growth is a complex process which study demands the use of many experimental approaches, and actually experimental analysis is in general expensive and/or difficult to be carried out and moreover no experiment can fully explain this process. For these reasons, mathematical modelling might theoretically combine a broad range of biological events for giving a better vision of the overall tumor progression including angiogenesis.

Since the seminal work of Greenspan [24], the mathematical modelling of avascular solid tumor growth has been rapidly expanding. Most models in this area consist of
nonlinear partial differential equation systems (e.g. see References [25–28]), and may
be described as a macroscopic approach. In this work we will focus on tumor-induced
angiogenesis and on some of the key aspects of vascular tumor growth, so we will omit a
bibliographical discussion about avascular tumor growth in this manuscript, but we refer
to exceptional work [3, 29–31] and the references therein.

Modelling tumor-induced angiogenesis has a well-established history beginning with the work of Balding and McElwain [32]. These authors proposed a simple model of tumor angiogenesis to describe experiments in which tumor cells implanted in the rabbit cornea stimulated formation, growth and migration of new blood vessels from the corneal limbus to the tumor [3]. Since the pioneer work of these authors [32], much of the mathematical modelling has focused on the way in which TAF initiate and coordinate capillary growth [33]. The next phase in tumor development, namely vascular tumor growth, has received less attention than avascular growth and angiogenesis in mathematical modelling literature. See also [34, 35] for early work on vascular tumor growth and invasion.

In order to show some examples in which mathematical modelling combined with
biological experiments have contributed to a better cancer understanding, we will
highlight some selected results centered on targeted therapies efficacy. In the papers
[36–38], the authors studied multiscale models to address questions related to prediction
of chemotherapy and radiotherapy efficacy. Likewise, Lignet et al. [39] and Panovska et
al. [40] developed models to investigate the combined effect of anti-proliferative therapy
with anti-TAF or anti-vascular therapy. Additionally, Billy et al. [41] studied the efficacy
of a new anti-angiogenesis treatment and provide some indications about the best way
to optimize this cancer treatment strategy. In our very recent work [42], we dealt with
tumor drug resistance modelling in hepatic gastrointestinal stromal tumor (GIST) metas-
tases, which exhibit resistance to two standard treatments: imatinib and sunitinib. Based
on an accurate analysis of medical images, we provided a patient-dependent model that
reproduces qualitatively and quantitatively spatial tumor evolution, as shown the follow-
up clinical data that we have carried out. Interestingly, specific aspects of tumor growth
as spatial heterogeneity and treatment failures could be explained by our model. See
Section Looking for a hybrid approach.

Additionally, also concerning the combination of mathematical modelling with clinical
data, another recent work from our group described the basic principles of mathemat-
ical modelling, as well as, the advantages, limitations and future prospects using both
oncology imaging and modelling [43]. Combination of imaging and modelling can resolve
complex problems and describe many aspects of tumor growth or response to treat-
ment, and therefore nowadays is possible to consider its clinical use in the medium term.
Accordingly, in the present work we also describe some aspects of tumor growth mod-
elling giving more details and, at the same time, we suggest some biological experiments,
as well as, the way of combining them to provide new insights in tumor progression,
angiogenesis and response to treatments. To do this, in the next section we focus on
deterministic models which analyze tumor growth as a macroscopic mechanical process
involving the whole tumor tissue rather than analyzing individual cells, in the same way
as models presented so far in this section. In addition, in the Section Looking for a hybrid
approach, we give some insights about how mathematical modelling should be modified
for fitting with the biological scenario of a specific cancer.
## Basic principles of mathematical modelling of tumor growth

Taking into account the complexity of the processes involved in all the stages of tumor growth, it is not difficult to understand that mathematical models are to a large extent basically phenomenological and simplified compared to what is happening in a biological context. Roughly speaking, from a biological point of view two main types of modelling can be distinguished: those that only consider the tumor as a whole, and those that also consider the spatial distribution of tumor components [43].

The key challenge is to develop mathematical models for tumor progression including the aspects mentioned above, in which optimal combination of drugs is simulated in order to suggest the best protocols for improving clinical outcomes. In addition, these models may be used to shed light about the existence of new mechanisms that could explain phenomena clinically observed and may have predictive capability. To do this, parameter estimation of these models is a key issue for which performing experiments or collecting clinical data are required. Then, in this work we want to emphasize two main aspects; the first one is the mathematical modelling of tumor progression; and the second one is how data coming from experimental biology or clinical data may contribute to parameter estimation.

## Non spatially-structured models of tumor growth based on ordinary differential equations

We present below some simple models considering only the change in volume at which tumor growth occurs, without taking into account tumor environment. We focus on models based on ordinary differential equations (ODE). To do this, we borrow some ideas by Byrne [44]. See Appendix A for a more detailed discussion.

One of the simplest models that can be used to describe the way in which the number of cells $N(t)$ within a solid tumor changes over time is the exponential growth law. In this model, there are no constraints on cell growth: all nutrients and other vital growth factors are assumed to be available in abundance. In consequence, the model predicts that the population will increase exponentially, without limit. See Fig. 1.

Whilst the exponential growth law is not realistic enough, it provides an accurate description of the early stages of a tumor development. In particular, reduced growth and eventual saturation are observed when avascular tumors are grown *in vitro* or when vascular tumors develop *in vivo*. This discrepancy arises because as the tumor increases in size, competition for nutrients and other vital resources, such as space, can no longer be neglected. A simple modification of the exponential growth law which takes account of competition for resources (without specifying what those resources are) is the logistic growth law.

Whilst the logistic growth law predicts almost exponential growth of small tumors and growth saturation when the tumor reaches its carrying capacity (see Fig. 1), the symmetry of $N(t)$ about its point of inflection means that it is not particularly flexible to fit or describe experimental data. A more general family of curves, which, depending on the choice of a parameter $\alpha$ can saturate more or less fastly than the logistic growth law, is given by

$$ \frac{dN}{dt} = \frac{k}{\alpha} N \left[ 1 - \left( \frac{N}{\theta} \right)^{\alpha} \right], \quad \text{with } N(t=0) = N_0, \qquad (1) $$
<figure><img src="image_4.png" /></figure>

where $k > 0$ represents the net rate at which tumor cells proliferate and $\theta > 0$ denotes the carrying capacity.

We remark that the logistic growth law is a special case of (1) (set $\alpha = 1$) and that the Gompertzian growth law is recovered in the limit as $\alpha \to 0^+$, which states that the net proliferation rate of tumor cells exponentially decreases with time.

In order to compare the four models presented above, we plot in Fig. 1 growth curves for each model for fixed values of the proliferation rate $k$ and the carrying capacity $\theta$.

These models are too simple to provide useful, reliable and reproducible information about the tumor from the phenomenological change in its volume, which varies between different tumors and patients. They may, however, be made more complex to incorporate tumor cell heterogeneity (coming from the cell cycle) or to take into account different processes such as angiogenesis; see Reference [43].

In order to get models with predictive capability an important issue is parameter estimation, which has to be performed for fitting modelling with experimental or clinical data. For example, estimating parameters $k$ and $\theta$ for the logistic model (see (13) in appenA) is roughly made by minimizing the error

$$E(k, \theta) = \sum_{i=1}^{M} (\hat{N}_i - N(t_i))^2, \qquad (2)$$

where $\hat{N}_i$, $i = 1, \dots, M$ are $M$ experimental or clinical measurements collected at different times, say $t_i$ ($i = 1, \dots, M$), during the tumor evolution, corresponding to the observed
numbers of tumor cells at the time instants $t_i$, and $N(t_i)$ are the numbers of tumor cells predicted by the model at time instants $t_i$. The minimization of the error $E(k, \theta)$ (2) is made in this case with respect to the parameters $k$ and $\theta$ in such a way that the responses of the model $N(t_i)$ approximately match the data $\hat{N}_i$. When the minimum of (2) is reached for some $(\hat{k}, \hat{\theta})$, the error $E(\hat{k}, \hat{\theta})$ is called the least-squares error whereas $(\hat{k}, \hat{\theta})$ is known as the least-squares estimate of the parameter $(k, \theta)$. It is worth noting that the optimization process can be quite complicated depending on the model, the number of parameters and of the availability of data; the interested reader is referred to the textbook [45].

Classically the mathematical models used for clinical applications are based on systems of non-linear ODE, which do not consider the spatial aspect of tumor growth. Despite of it, these models have great interest in biological applications and typically they are parametrized using statistical methods and may provide a prognosis of tumor volume, but neither shape nor location of tumor can be estimated. See references [46, 47].

In order to appreciate the impact of parameter estimation in the clinical context, for example we refer to a recent meta-analysis [48] which describes natural development of meningiomas, a kind of neurological tumor. In [48], 22 studies reporting 675 patients with untreated meningiomas were found, followed by serial magnetic resonance imaging (MRI) during 5 years. From the analysis, authors show that tumors which initial diameter was > 2.5 cm exhibit a linear growth rate > 10% per year, which leads to the highest risk for developing progressive symptoms. By contrast, untreated meningiomas which initial diameter was < 2.5 cm do not show tumor growth over a follow-up period of 4.6 years. Then, the finding of the parameter 2.5 cm as a threshold for meningiomas in vivo allows to describe tumor biologically more aggressive and classify patients who might require or not surgery, which may avoid excessive intervention.

## Spatially-structured models of tumor growth based on partial differential equations

We describe below generic modelling, based on partial differential equations (PDE), which consider the effect that changes in the composition of the medium surrounding the tumors have on their growth. We focus on the study of three key issues: vascular tumor growth, tumor-induced angiogenesis and invasion, and efficacy of treatments such as anti-proliferative and anti-angiogenic therapies. See Appendix B for a more detailed discussion.

The most basic principle for all quantitative models is conservation of mass [49]. Conservation of mass of a component in a dynamic and open system states that:

$$ \begin{pmatrix} \text{Net rate of} \\ \text{change} \\ \text{of mass} \\ \text{of component} \\ \text{in the system} \end{pmatrix} = \begin{pmatrix} \text{Mass flow} \\ \text{of the} \\ \text{component} \\ \text{into} \\ \text{the system} \end{pmatrix} - \begin{pmatrix} \text{Mass flow} \\ \text{of the} \\ \text{component} \\ \text{out of} \\ \text{the system} \end{pmatrix} + \begin{pmatrix} \text{Rate of} \\ \text{production} \\ \text{of the} \\ \text{component by} \\ \text{transformations} \end{pmatrix} - \begin{pmatrix} \text{Rate of} \\ \text{consumption} \\ \text{of the} \\ \text{component by} \\ \text{transformations} \end{pmatrix} $$

The local mass balances are the mathematical form of equality, which can be written as

$$ \frac{\partial C}{\partial t} = -\nabla \cdot \mathbf{J} + r, \qquad (3) $$

where $t$ is time; $C$ is the concentration of the component in the system; $\mathbf{J}$ is the mass flux of the component; and $r$ is the net production (or growth) rate of the component. This is the equation of continuity for a component, either tumor cell densities or chemical
concentration of a growth rate limiting (e.g. oxygen or glucose), or of TAF, or of drugs against cancer.

The transport processes that regularly are considered in cancer models are advection (for cell densities) and molecular diffusion (for chemical concentrations). The general expression to model the specific mass flux of a component is

$$J = vC - D\nabla C, \qquad (4)$$

where v is the velocity field due to the mechanical forces which act on the system (e.g. the pressure within the tumor) and D is the diffusion coefficient.

We assume that transport of tumor cells is driven only by advection. Consequently, applying the principle of mass balance, the spatiotemporal dynamics of cell populations are formulated as the generic PDE:

$$\frac{\partial f}{\partial t} = -\nabla \cdot (v f) + r(f), \qquad (5)$$

where r(f) is the net growth rate of the cell density f, which depends on the cell population considered, and $\frac{\partial f}{\partial t}$ stands for the partial derivative with respect to time (representing the rate of change of the mass density f). For example, if f represents the proliferating tumor cells density, designed as f = P, then

$$r(P) = \Gamma P - \delta P. \qquad (6)$$

In replacing r(P) in (5), one arrives to

$$\frac{\partial P}{\partial t} + \nabla \cdot (vP) = \Gamma P - \delta P. \qquad (7)$$

This formula is called a reaction-advection equation, because of the presence of the advection term $\nabla \cdot (vP)$ and the reaction term $\Gamma P - \delta P$. The term $\Gamma P$ represents the net tumor growth rate, which is determined by local oxygen concentration. For example

$$\Gamma = \gamma_0 \frac{1 + \tanh(R(C - C_{hyp}))}{2}, \qquad (8)$$

where C represents the local oxygen concentration, $C_{hyp}$ is a parameter which describes the sensitivity of cells to hypoxia and $\gamma_0$ stands for the maximal proliferation rate. $\Gamma$ is chosen according to (8) in such a way that it satisfies the following properties: $\Gamma \to \gamma_0$ as $C \gg C_{hyp}$, i.e. cells proliferate as the oxygen concentration is high enough, and $\Gamma \to 0$ as $C \ll C_{hyp}$, i.e. cells do not proliferate as the oxygen concentration is low enough (hypoxia). Moreover $R > 0$ is a numerical smoothing parameter, in such a way that $\Gamma$ be a smooth version of a Heaviside function.

The term $\delta P$ represents the necrosis due to a generic anti-proliferative drug, where the parameter $\delta$ represents the doses of such drug. Eq. 7 is a basic form for modelling tumor proliferation and obviously there are others terms which we could add, for example, a term taking into account the passage to quiescence (see Appendix B).

Similarly, we assume that transport of chemical concentrations is driven only by diffusion. Consequently, applying the principle of mass balance, the spatiotemporal dynamics of chemical concentrations are formulated as the generic PDE:

$$\frac{\partial g}{\partial t} = \nabla \cdot (D \nabla g) + r(g), \qquad (9)$$
where $r(g)$ is the reaction term of the chemical concentration $g$. This formula is called reaction-diffusion equation. For example, the local oxygen concentration is given by:

$$ \frac{\partial C}{\partial t} - \nabla \cdot (D_C \nabla C) = -g_C + h_C, \qquad (10) $$

where $D_C$ stands for the diffusion coefficient of local oxygen in the blood flow, $h_C$ is the production of oxygen in the functional blood vessels, and $g_C$ is the consumption of oxygen mainly due to proliferation tumor cells.

Likewise, the anti-proliferative drug concentration is given by

$$ \frac{\partial M}{\partial t} - \nabla \cdot (D_M \nabla M) = -\lambda M - aMP + \rho \delta_M, \qquad (11) $$

where $D_M$ is the diffusion coefficient of anti-proliferative drug in the blood flow, $\lambda$ is the mean-life of such drug, $a$ is the consumption rate due mainly to proliferation tumor cells, $\rho$ represents the blood flow, and finally $\delta_M$ denotes the doses of the drug administered at different time intervals. Similar equations also hold for the blood flow $\rho$, the anti-angiogenic drug concentration and the TAF concentration, particularly the VEGF.

Mathematical modelling described above relies on systems of non-linear PDE, in which a set of parameters takes into account the complexity of the underlying biological phenomena. In order to apply such models in practical situations, these parameters need to be identified, that is to say, biologically meaningful values have to be estimated. One way to determine their values is by means of inverse problems theory, exploiting data coming from medical images, as achieved in References [50, 51]. One of the main difficulties is that the amount of data for system identification is scarce. Although medical scans allow an accurate localization of tumors in space, little information can be inferred regarding cellular nature or nutrient distribution inside tumors. In addition, usually only two scans are available before treatment, which makes estimation of tumor evolution a challenging problem.

On the other hand, retrieving tumor shape evolution may provide useful information since PDE based models are spatially distributed. We will return to this important issue in the next section. One possible approach to formulate the inverse problem is by optimal control theory, as was done for instance in [50]. In this approach a PDE-constrained optimization problem has to be solved. For instance, one can consider a functional to be minimized by matching spatiotemporal evolution of tumor density predicted by mathematical modelling with the corresponding tumor density maps estimated from serial scans for one particular patient. In order to carry out the minimization process an adjoint-based algorithm is used for evaluating the gradient of the functional. This algorithm is expensive, since each optimization iteration requires solving a number of forward problems equal to number of variables. A different approach, as was done for instance in [51], consists in using the difference between a variable predicted by the model and its corresponding observed value (the residuals of the model) within a Newton method to solve the inverse problem. This identification procedure is based on proper orthogonal decomposition (POD), a way to identify complex natural processes with models that are intrinsically much simpler than biological scenarios. In this approach the solution is sought in a given low-dimensional functional space, which basis gives an optimal representation of a sufficiently large number of solution samples. In particular, the parameter space is sampled in such a way that all the possible different biological behaviors are represented.
We finish this section by showing some results which are part of a work in progress. We have carried out some numerical simulations of a tumor growth model under normal conditions, i.e., in absence of treatments. These simulation results are depicted in Fig. 2, which shows the distribution of endothelial cells (right), as well as, the distribution of the sum of the proliferative and quiescent tumor cells (left). The time unit is 12 hours. Because of the configuration (four blood vessels initially placed far from tumor cells), the proliferative and quiescent tumor cells initially decrease, which can be clearly seen in Fig. 2. At the same time, this would produce an increment of the VEGF in the region where the quiescent cells are concentrated. This would trigger the proliferation of endothelial cells making them to move by chemotaxis towards the tumor (see Fig. 2). This process takes place between $t = 30[12h]$ and $t = 140[12h]$, generating four sources of nutrients which are enough to increase the oxygen concentration over the hypoxic level. Further, as a result a lack of production in VEGF may be observed which it is thought to be consumed by endothelial cells. At this stage, the main effect is the reactivation of the proliferation of peripheral tumor cells due to high levels of nutrients in a neighborhood of the new blood vessels. However, after some time the oxygen concentration would be consumed by proliferative cells to levels that produce the usual distribution in three layers of tumor cells (see Fig. 3) described in the schedules of fully developed tumors [44]. A diminution in the spread of endothelial cells can also be observed in Fig. 2 as a consequence of tumor growth, which is reflected as a reduction in density of the four blood vessels between $t = 159[12h]$ and $t = 250[12h]$. Accordingly vascular collapse and regression in tumors has been suggested [52], and it is thought that is a consequence of the biomechanical stresses and of the action of the interstitial pressure [53].

## Looking for a hybrid approach

Parameter estimation is highly required in order to improve models for tumor progression and provide new clinical insights. In this regard, we have defined hybrid approach as the feedback capacity between mathematical modelling and biological experimental design. This approach may lead to a better understanding cancer, by fitting mathematical parameters with biological processes, in order to achieve accurate biological predictive models. In this section we will give two examples which illustrate the hybrid approach.

The first example is concerned with tumor cell proliferation. As mentioned in Section Introduction, tumor cell proliferation can be correlated to oxygen concentration, and indeed one could think that the higher the oxygen concentration the higher the tumor cell proliferation. However, this is not always like that. In particular, elevation of anaerobic metabolism in tumor cells during avascular or even vascular tumor growth has been described, a phenomenon called Warburg effect [54], which involves a shift in metabolism away from oxidative phosphorylation (i.e., aerobic) towards anaerobic glycolysis. Therefore, tumor cells are resistant to hypoxia which indeed, depending on the threshold, may stimulate cell proliferation. In this regard, unpublished results from our lab using two ovarian cancer cell lines, HEY and UCI, were cultured under two different oxygen concentrations, 21% (normoxia) and 5% oxygen (hypoxia) respectively. Under these conditions, UCI cells were significantly more sensitive to oxygen variations than HEY cells. Indeed under 21% oxygen UCI cell proliferation was accelerated about 3 times compared to 5% oxygen, whereas in the case of HEY cells, the proliferation did not significantly change from both different oxygen concentrations (see Fig. 4). Considering these last results into
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Four snapshots of tumor evolution in absence of treatment at times $t = 10, 70, 159$ and $250[12 h]$. (Left) Spatial distribution of proliferative plus quiescent tumor cells. (Right) Spatial distribution of endothelial cells</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Three layers fully developed distribution of the tumor at time $t = 270[12 h]$. Spatial distribution of proliferative cells (up), quiescent cells (middle) and necrotic cells (down)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Different cell proliferation in ovarian cancer cells. Two cell lines derived from human ovarian cancers were used (HEY and UCI). Cells were cultured under 21% and 5% oxygen during 0, 3, 6, 12 and 24 hours. Cell proliferation was analyzed by bromouridine incorporation as previously reported [58]. Data is presented as the logarithm of $t_{50}$ (replication time) ± SEM. *N* = 4 per group and analyzed time. \*p < 0.05 vs HEY at 21% oxygen. \*\*p < 0.05 vs UCI at 21% oxygen</figcaption></figure>

the mathematical modelling described before, net tumor growth rate $\Gamma$ in Eq. 8 need to be remodelled. In Fig. 4, it is clearly seen that UCI cells behave as one would expect, i.e., the higher the oxygen concentration the higher the tumor cell proliferation; in contrast, no matter the oxygen concentration, HEY cells undergo the same proliferation rate. Consequently, whereas net proliferation rate for UCI cells behaves as function $\Gamma$ defined by (8), one should modify the equation for HEY cells. In addition, parameters $C_{hyp}$ and $\gamma_0$ could be estimated using the same experimental setting.

The second example is concerned with tumor drug resistance in gastrointestinal stromal tumors (GIST) liver metastases. In our very recent work [42], we have modelled and simulated resistance to two standard treatments: imatinib and sunitinib. In this regard, the specific tyrosine kinase inhibitors, imatinib, is used as standard first-line treatment. In most of the cases this drug leads to tumor reduction during several months, but then, most of the patients relapse and the tumor is no longer responding to treatment. Then the standard care switches to a second-line treatment using a multi-targeted tyrosine kinase inhibitor (sunitinib), which has both cytotoxic and antiangiogenic effects. After use this second drug, metastasis is controlled for some additional time before a new therapeutic failure occurs. Considering this clinical data, we have split the density *P* of proliferative cells (*P*-cells) into 3 subpopulations $P_1$, $P_2$ and $P_3$, such that $P = P_1 + P_2 + P_3$, where

* $P_1$ denotes the fraction of proliferative cells that are sensitive to the first-line treatment, based on imatinib molecule, and also to the second-line treatment, based on sunitinib;

* $P_2$ describes the density of proliferative cells that are resistant to the first-line treatment and sensitive to the second-line treatment; and

* $P_3$ stands for proliferative cells that are resistant to both treatments.
According to the clinical observations, it seemed relevant to consider that the three cell sub-populations are present when the GIST metastasis is detected. In the mathematical models, we describe tumor growth evolution in terms of tumor area, as well as, its spatial structure and compared with the CT-scans measurements. We believe that intra-tumoral spatial heterogeneity may be related to an increase in the cellular activity, which would mean that a resistant phenotype should be emerging. Indeed according to our modeling, such heterogeneity may be seen as the first stage of the treatment failure, while Response Evaluation Criteria In Solid Tumors (RECIST) [55] do not bring any information about this fact. Therefore our work can be seen as a first step in developing new tools to evaluate tumor response to treatment based on tyrosine kinase inhibitors. Consequently, as mentioned in the previous section, retrieving the evolution of the tumor shape may provide useful information for clinical applications.

According to the discussion above, data coming from clinical observations were crucial in order to develop a complex model that takes into account tumor heterogeneity and treatments resistance. In addition, our model parameters were calibrated in a patient-dependent way, in particular, estimating parameters for two patients and fitted with their specific tumor evolution. Consequently, in this example, mathematical modelling discussed before (in particular Eq. 7 among others) was adapted from cross-talk between physicians and mathematicians based on clinical observations, leading to a new modelling capable to provide better insights for enhancing clinical meaning of mathematical modelling.

## Conclusions and future directions

In this work we described and discussed mathematical modelling of tumor progression, including some aspects as vascular tumor growth, angiogenesis and effectiveness and/or resistance to therapies. We also emphasized what we have called *hybrid approach* in which inverse and progressive confirmation by biological experimentation and use of data from patients with cancer is suggested for enhancing clinical meaning of mathematical modelling. We believe that hybrid approach is a key tool for better understanding cancer, which has been illustrated by two specific examples carry out by our group.

Hybrid approach also involves a big challenge for researchers in cancer field, since despite cancer is deeply analyzed from the particular perspective in each disciplines, when mathematical modelling and biological research face each other, an easier and common language is required in order to go beyond understanding of the pathology. This is not easy to deal with, because there are many differences in the expertise, nomenclature, language and even in the way of thinking that could slow down the communication. Then, one of the main challenges to be overcame is listening and learning each other, maintaining academic discussions where we could bridge the gap on these differences. With hope we have seen a growing number of publications and working groups which are seeking to combine science, but it is still insufficient. In addition, at the same time another step that can be done is the inclusion of physicians in this discussion/analysis in order to look forward potential clinical interest and application of the generated knowledge from basic sciences. These interactions not only will enhance knowledge of cancer, but also improve clinical meaning of mathematical modelling.

Parameter estimation in mathematical modelling may ask very specific questions that need to be designed and tested *in vivo* or *in vitro*. These questions might constitute a
challenge for biologist due to limitation in available technology, and also because some parameters are not easy to determine or even not necessarily represent an *in vivo* system, but they are only associated with. For instance, reaction-diffusion Eq. 10 models the balance between variation, diffusion, production and consumption of local oxygen. Despite in mathematical modelling constitutes a usual equation, in biological terms series of experimentation are required in order to analyze how oxygen is delivered towards tumor, how is diffused from vessels to tumor, how tumor cells may control blood flow and oxygen delivery, among others. In particular, future directions should focus in a suitable design of *in vitro* or *in vivo* models for answering such mathematical questions. In this regard, another experimental setting we would like to suggest is isolation of endothelial cells derived from patient's tissue after biopsies, as it has been already previously developed [56, 57]. This particular cell type might be used for testing specific questions such as how and how much VEGF released from tumor cells may enhance tube formation capacity of that endothelium. Interestingly, anti-angiogenic drugs could also be tested using this experimental setting, which results might be used for parameter estimation.

Future directions in tumor growth modelling should include a cross-talk between biological and theoretical researchers in order to unify effort and generate an active interaction between them. A more common language and understanding each other will be helpful in order to enhance generation of high impact knowledge in tumor growth. In this regard, combination of *in vivo* and *in vitro* experimentation with theoretical analysis is a major challenge which hopefully is changing. In this manuscript we intend to contribute in this field, which of course is limited in comparison with what is needed to be done in order to achieve a better understanding of tumor growth, but it constitutes an example which may alert both mathematicians and clinicians.

## Appendix A: Tumor growth laws based on ordinary differential equations

In this appendix we give some details about of the non spatially-structured models of tumor growth based on ordinary differential equations.

One of the simplest models that can be used to describe the way in which the number of cells $N(t)$ within a solid tumor changes over time is the exponential growth law which states

$$ \frac{dN}{dt} = kN, \quad \text{with } N(t=0) = N_0, \qquad (12) $$

whose solution is $N(t) = N_0 \exp(kt)$. In Eq. 12, $k > 0$ represents the net rate at which the cells proliferate, and $N_0$ denotes the number of cells initially present within the tumor. In this model, there are no constraints on cell growth: all nutrients and other vital growth factors are assumed to be available in abundance. In consequence, the model predicts that the population will increase exponentially, without limit.

Whilst the exponential growth law is not realistic enough, it provides an accurate description of the early stages of a tumor development. In particular, reduced growth and eventual saturation are observed when avascular tumors are grown *in vitro* or when vascular tumors develop *in vivo*. This discrepancy arises because as the tumor increases in size, competition for nutrients and other vital resources, such as space, can no longer be neglected. A simple modification of (12) which takes account of competition for resources
(without specifying what those resources are) is the logistic growth law

$$ \frac{dN}{dt} = kN \left(1 - \frac{N}{\theta}\right), \quad \text{with } N(t=0) = N_0. \qquad (13) $$

By using elementary calculus one arrives to

$$ N(t) = \frac{\theta N_0}{N_0 + (\theta - N_0) \exp(-kt)} \rightarrow \theta \quad \text{as} \quad t \rightarrow \infty. \qquad (14) $$

In (13), $\theta > 0$ represents the carrying capacity of the tumor cells population. Whilst the logistic growth law predicts almost exponential growth of small tumors and growth saturation when the tumor reaches its carrying capacity ($N = \theta$), the symmetry of $N(t)$ about its point of inflection (where $\frac{d^2N}{dt^2} = 0$ and $N = \theta/2$) means that it is not particularly flexible to fit or describe experimental data. A more general family of curves, which, depending on the choice of a parameter $\alpha$ can saturate more or less fastly than (13), is given by (1), which solution is

$$ N(t) = \theta \left( \frac{N_0^{\alpha}}{N_0^{\alpha} + (\theta^{\alpha} - N_0^{\alpha}) \exp(-kt)} \right)^{1/\alpha}. \qquad (15) $$

We remark that the logistic growth law is a special case of (1) ($\alpha = 1$) and that the Gompertzian growth law is recovered in the limit as $\alpha \rightarrow 0^+$. Gompertzian growth law is given by

$$ \frac{dN}{dt} = kN \ln\left(\frac{\theta}{N}\right), \quad \text{with } N(t=0) = N_0, \qquad (16) $$

whose solution is

$$ N(t) = N_0 \exp\left(\frac{r}{k}(1 - \exp(-kt))\right), \quad \text{for } r = k \ln\left(\frac{\theta}{N_0}\right). \qquad (17) $$

The Gompertzian growth law states that the net proliferation rate of tumor cells exponentially decreases with time. This is apparent because of that

$$ \frac{dN}{dt}(t) = k \ln\left(\frac{\theta}{N(t)}\right) = r \exp(-kt), \qquad (18) $$

and noting that $k \exp(kt) \ln(\frac{\theta}{N(t)}) = k \ln(\frac{\theta}{N_0})$ (i.e. $k \exp(kt) \ln(\frac{\theta}{N(t)})$ remains constant in time).

In order to compare the four models presented above, we plot in Fig. 1 growth curves for each model for fixed values of the proliferation rate $k$, the carrying capacity $\theta$ and the initial condition $N_0$.

## Appendix B: Tumor growth laws based on partial differential equations

In this appendix we give some details about of the spatially-structured model of tumor growth based on partial differential equations, from which we have obtained the simulation results depicted in Figs. 2 and 3. It is worth noting that this is a work in progress, therefore the previous simulation results are not yet definitive. Moreover, we have not included yet in this model the influence of drugs.

Our model consists in a system of partial differential equations of type advection-diffusion-reaction. This system describes the evolution of the cell densities and of the molecules (oxygen and VEGF), as well as, the interaction between them. The rules that regulate such interactions are, in general, non-linear expressions of the coefficients in the equations. All of these variables undergo temporal dependency (noted by $t$) and spatial
<table><caption>**Table 2** Notations used in the model</caption><thead><tr><th>Variable</th><th>Description</th></tr></thead><tbody><tr><td>Ω</td><td>Computational domain</td></tr><tr><td>P</td><td>Proliferating tumor cells density</td></tr><tr><td>Q</td><td>Quiescent tumor cells density</td></tr><tr><td>N</td><td>Necrotic tumor cells density</td></tr><tr><td>S</td><td>Host cells density</td></tr><tr><td>ρ</td><td>Endothelial cells density</td></tr><tr><td>α</td><td>VEGF concentration</td></tr><tr><td>C</td><td>Oxygen concentration</td></tr><tr><td>v</td><td>Advection velocity</td></tr><tr><td>k</td><td>Permeability of the medium</td></tr></tbody></table>

dependency (noted by $\vec{x}$). However, to simplify the notation, such a dependence will be omitted and by instance $P$ shall be written instead of $P(\vec{x}, t)$.

The notations used in the model are shown in Table 2, and the full model proposed for the evolution of tumor growth coupled with tumor angiogenesis is shown in Table 3.

In above equations, $\Gamma$ represents the specific tumor proliferation rate, given by Eq. 8; $f_{QP}$ the specific rate of transition from the quiescent to the proliferative state; conversely, $f_{PQ}$ the specific rate of transition from the proliferative to the quiescent state; $f_{PN}$ and $f_{QN}$ represent the specific rate of necrosis, from the proliferative and the quiescent state, respectively. These transition functions are given by:

$$f_{QP} = \gamma_0 \frac{1 + \tanh(R(C - C_{hyp}))}{2}, \qquad (19)$$

$$f_{PQ} = \gamma_0 \left( \frac{1 - \tanh(R(C - C_{hyp}))}{2} \right) \left( \frac{1 + \tanh(R(C - C_{hyp}^{sev}))}{2} \right), \qquad (20)$$

$$f_{PN} = \gamma_0 \left( \frac{1 - \tanh(R(C - C_{hyp}^{sev}))}{2} \right), \qquad (21)$$

$$f_{QN} = \gamma_0 \frac{1 - \tanh(R(C - C_{ nec} ))}{2}, \qquad (22)$$

where $C_{hyp}^{sev}$ is the severe hypoxia threshold, that is to say, the oxygen concentration under which proliferative tumor cells die, and similarly $C_{ nec}$ is the oxygen concentration under which quiescent tumor cells die. One obviously has that $0 < C_{ nec} < C_{hyp}^{sev}$ to take into account the fact that it is more difficult the quiescent tumor cells die.

<table><caption>**Table 3** Summary of the equations used for the numerical simulation depicted in Figs. 2 and 3</caption><thead><tr><th>Variable</th><th>Equation</th></tr></thead><tbody><tr><td>Proliferative tumor cells density</td><td>P<sub>t</sub> + ∇ · (vP) = ΓP + f<sub>QP</sub>Q - (f<sub>QP</sub> + f<sub>PN</sub>)P</td></tr><tr><td>Quiescent cells density</td><td>Q<sub>t</sub> + ∇ · (vQ) = f<sub>QP</sub>P - (f<sub>QP</sub> + f<sub>QN</sub>)Q</td></tr><tr><td>Necrotic tumor cells density</td><td>N<sub>t</sub> + ∇ · (vN) = f<sub>PN</sub>P + f<sub>QN</sub>Q</td></tr><tr><td>Host cells density</td><td>S<sub>t</sub> + ∇ · (vS) = 0</td></tr><tr><td>Endothelial cells density</td><td>ρ<sub>t</sub> + ∇ · (vρ) + ∇ · (g<sub>CH</sub>) = f<sub>ρρ</sub>ρ</td></tr><tr><td>[O<sub>2</sub>]</td><td>C<sub>t</sub> - ∇ · (D<sub>C</sub>∇C) = -g<sub>C</sub> + h<sub>C</sub></td></tr><tr><td>[VEGF]</td><td>α<sub>t</sub> - ∇ · (D<sub>α</sub>∇α) = -g<sub>α</sub> + h<sub>α</sub></td></tr><tr><td>Advection velocity</td><td>v = -k∇φ</td></tr><tr><td>Pressure</td><td>-∇ · (k∇φ) = ΓP + f<sub>ρρ</sub>ρ - ∇ · (g<sub>CH</sub>)</td></tr></tbody></table>
The transport of endothelial cells, which are responsible by tumor vascularization, is mainly driven by chemotaxis toward the source of VEGF. This transport term is modelled by

$$g_{CH} = \chi \rho \nabla \alpha, \qquad (23)$$

whereas the source of endothelial cells is given by

$$f_{\rho\rho} = \chi_{prol} \left( \frac{\alpha}{\alpha_{half} + \alpha} \right). \qquad (24)$$

In Eq. 23, $\chi$ denotes the sensitivity to the chemotaxis, which is given by

$$\chi = \chi_{chemo} (1 - \rho), \qquad (25)$$

with $\chi_{chemo}$ is the maximum effect of VEGF on endothelial chemotaxis; $\chi_{prol}$ is the maximum effect of VEGF on the proliferation of endothelial cells; $\alpha_{half}$ is a half-proliferation constant, i.e., it denotes the VEGF concentration such that the proliferation of endothelial cells is exactly half of its maximum (0.5 $\chi_{prol}$).

For the oxygen concentration $C$ the source and the consumption are given respectively by:

$$g_C = (\lambda_C + \varepsilon_C P + \beta_C Q)C, \qquad h_C = C_{sou} \rho \left(1 - \frac{C}{C_{max}}\right), \qquad (26)$$

where $\lambda_C$ is the oxygen half-life, $\varepsilon_C$ is the consumption rate by the proliferating cells, $\beta_C$ is the consumption rate by the quiescent cells, $C_{sou}$ the maximum oxygen level available in functional blood vessels and $C_{max}$ the maximum oxygen level in the host tissue.

In a similar way, the source and consumption terms for the VEGF concentration are given respectively by:

$$g_{\alpha} = (\lambda_{\alpha} + \pi \rho) \alpha, \qquad h_{\alpha} = \alpha_{sou} Q \left( \frac{1 - \tanh(R(C - C_{hyp}))}{2} \right) \left( 1 - \frac{\alpha}{\alpha_{max}} \right) \qquad (27)$$

where $\lambda_{\alpha}$ is the VEGF half-life, $\alpha_{sou}$ is the source of VEGF located at quiescent tumor cells, $\alpha_{max}$ the maximum VEGF concentration and $\pi$ denotes the binding rate of VEGF to the receptors of the endothelial cells.

We have assumed that the total number of cells per volume unit is constant (and without loss of generality equal to 1):

$$P + Q + S + N + \rho = 1. \qquad (28)$$

Adding the equations for the cell densities $P$, $Q$, $S$, $N$ and $\rho$ (see Table 3), we get:

$$\nabla \cdot \nu = \Gamma P - \nabla \cdot (g_{CH}) + f_{\rho\rho}\rho. \qquad (29)$$

Thus by imposing the Darcy's law $\nu = -k\nabla\varphi$ for the advection velocity, where $k$ is the permeability of the medium, we deduce the equation for the pressure $\varphi$ (see Table 3).

Finally, we have noted $D_j$ the diffusion coefficients for the different molecules $j = 1, 2$ (oxygen and VEGF):

$$D_j = D_{j,max}(1 - \varepsilon(P + Q + N)), \qquad (30)$$

where $D_{j,max}$ is the maximum diffusion coefficient of the molecule $j$ and $\varepsilon > 0$ is a parameter that represents the percentage of diminution of the diffusion within the tumor, due to the higher density of this last one with respect to the host tissue.
<table><caption>**Table 4** Summary of the parameter values used in the simulation depicted in Figs. 2 and 3</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>C<sub>hyp</sub></td><td>Threshold of hypoxia</td><td>5.5</td><td>M</td></tr><tr><td>C<sup>sev</sup><sub>hyp</sub></td><td>Threshold of severe hypoxia</td><td>1.52</td><td>M</td></tr><tr><td>C<sub>nec</sub></td><td>Threshold of necrosis</td><td>0.02</td><td>M</td></tr><tr><td>D<sub>C,max</sub></td><td>Maximum diffusion of oxygen</td><td>1</td><td>mm<sup>2</sup> h<sup>-1</sup></td></tr><tr><td>D<sub>α,max</sub></td><td>Maximum diffusion of VEGF</td><td>0.1875</td><td>mm<sup>2</sup> h<sup>-1</sup></td></tr><tr><td>ε</td><td>Percentage of loss of diffusion</td><td>20</td><td>%</td></tr><tr><td>C<sub>sou</sub></td><td>Oxygen concentration in functional blood vessels</td><td>8</td><td>M h<sup>-1</sup></td></tr><tr><td>λ<sub>C</sub></td><td>Degradation rate of oxygen</td><td>0.01</td><td>h<sup>-1</sup></td></tr><tr><td>ε<sub>C</sub></td><td>Rate of oxygen consumption by proliferative cells</td><td>3</td><td>cells<sup>-1</sup> mm<sup>2</sup> h<sup>-1</sup></td></tr><tr><td>β<sub>C</sub></td><td>Rate of oxygen consumption by quiescent cells</td><td>1.5</td><td>cells<sup>-1</sup> mm<sup>2</sup> h<sup>-1</sup></td></tr><tr><td>C<sub>max</sub></td><td>Maximum concentration of oxygen</td><td>8</td><td>M</td></tr><tr><td>λ<sub>α</sub></td><td>Degradation rate of VEGF</td><td>1.25 × 10<sup>-4</sup></td><td>h<sup>-1</sup></td></tr><tr><td>π</td><td>Binding rate of VEGF to endothelial cells</td><td>0.791</td><td>cells<sup>-1</sup> mm<sup>2</sup> h<sup>-1</sup></td></tr><tr><td>α<sub>sou</sub></td><td>Production rate of VEGF</td><td>2.11</td><td>M cells<sup>-1</sup> mm<sup>2</sup> h<sup>-1</sup></td></tr><tr><td>α<sub>max</sub></td><td>Maximum concentration of VEGF</td><td>2.11</td><td>M</td></tr><tr><td>α<sub>half</sub></td><td>Shape parameter of f<sub>ρρ</sub></td><td>1.25 × 10<sup>-3</sup></td><td>M</td></tr><tr><td>χ<sub>chemo</sub></td><td>Maximum effect of VEGF on endothelial chemotaxis</td><td>1.25 × 10<sup>-2</sup></td><td>mm<sup>2</sup> M<sup>-1</sup> h<sup>-1</sup></td></tr><tr><td>χ<sub>prol</sub></td><td>Maximum effect of VEGF on endothelial proliferation</td><td>5 × 10<sup>-2</sup></td><td>h<sup>-1</sup></td></tr></tbody></table>

For the boundary conditions, we assumed that the molecules will not go out or in through of the boundary of the computational domain, that is to say, we assumed homogeneous Neumann boundary conditions for C and α.

For the pressure field $\varphi$ we imposed $\varphi|_{\partial\Omega} = 0$. This homogeneous Dirichlet condition is used since we consider that the domain of interest is not isolated, and the outer medium does not impose a pressure on the tumor. This assumption is valid for small tumors that are not mechanically constrained by the extratumoral region.

Finally, Table 4 below shows the parameter values for which simulation depicted in Figs. 2 and 3 was carried out. In Table 4 we use the following notations: h(hour), mm(millimeter), M(molarity).

**Competing interests**

The authors declare that they have no competing interests.

**Authors' contributions**

PC, AC, CC and CQ designed and simulated mathematical modelling. CE designed and conceived biological background. PC and CE designed and conceived the study. All authors approved the final manuscript.

**Acknowledgements**

The work of P. C. was founded by PIA-CONICYT grant PFBasal-01, DIUBB GI 153209/C and DIUBB GI 152920/EF. The work of A. C. was founded by DIUBB 121909 GI/C. C. E. thanks partial support from Ecos-Conicyt C13E05. He is also partially supported by PFBasal-01 and PFBasal-03 projects and by Fondecyt 1140773. The work of C. E. was founded by Fondecyt 1140586, Fondequip EQM140104, DIUBB GI153109/EF, GI 152920/EF. We would like to acknowledge Dr. Gareth Owen from Pontificia Universidad Católica de Chile for donating us ovarian cancer cells. Finally, we would like to thank you to Dr. Andrés Rodríguez for your invaluable comments on this work, as well as, your English corrections.

**Author details**

<footnote n="1">¹Department of Basic Sciences, Faculty of Sciences, Universidad del Bío-Bío, Campus Fernando May, Av. Andrés Bello s/n, Casilla 447 Chillán, Chile. ²Centre for Biotechnology and Bioengineering, University of Chile, Beaucheff 851, Santiago, Chile. ³Group of Applied Mathematics (GMA), Chillán, Chile. ⁴Group of Investigation in Tumor Angiogenesis (GIANT), Chillán, Chile. ⁵Department of Mathematical Engineering (DIM) and Center for Mathematical Modelling (CMM), University of Chile, (UMI CNRS 2807), Beaucheff 851, P.O. Box 170-3, Correo 3 Santiago, Chile. ⁶Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie and Mathematical Neuroscience Team, CIRB, Collège de France, (UMR CNRS 7598), 4 place de Jussieu, F-75005 Paris, France. ⁷Group of Research and Innovation in Vascular Health (GRIVAS Health), Chillán, Chile.</footnote>
Received: 16 March 2015 Accepted: 12 June 2015

Published online: 02 July 2015

References

1. Organization WH. Cancer. http://www.who.int/cancer/en/ (2015), Accessed 23 Jun 2015.

2. Zheng X, Wise SM, Cristini V. Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method. Bull. Math. Biol. 2005;67(2):211–59.

3. Byrne HM. Dissecting cancer through mathematics: from the cell to the animal model. Nat. Rev. Cancer. 2010;10(3): 221–30.

4. Casciari JJ, Sotirchos SV, Sutherland RM. Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH. J. Cell. Physiol. 1992;151(2):386–394.

5. Heddleston JM, Li Z, Lathia JD, Bao S, Hjelmeland AB, Rich JN. Hypoxia inducible factors in cancer stem cells. Br. J. Cancer. 2010;102(5):789–95.

6. Sutherland RM. Cell and environment interactions in tumor microregions: the multicell spheroid model. Science. 1988;240(4849):177–84.

7. Rey S, Semenza GL. Hypoxia-inducible factor-1-dependent mechanisms of vascularization and vascular remodelling. Cardiovasc Res. 2010;86(2):236–42.

8. Fox SB. Assessing tumor angiogenesis in histological samples. Methods Mol. Biol. 2009;467:55–78.

9. Kim BJ, Forbes NS. Single-cell analysis demonstrates how nutrient deprivation creates apoptotic and quiescent cell populations in tumor cylindroids. Biotechnol. Bioeng. 2008;101(4):797–810.

10. Macklin P, McDougall SR, Anderson ARA, Chaplain MAJ, Cristini V, Lowengrub JS. Multiscale modelling and nonlinear simulation of vascular tumour growth. J. Math. Biol. 2009;58(4-5):765–98.

11. Carmeliet P, de Almodovar CR. Vegf ligands and receptors: implications in neurodevelopment and neurodegeneration. Cell Mol Life Sci. 2013;70(10):1763–78.

12. Bao S, Wu Q, Sathornsumetee S, Hao Y, Li Z, Hjelmeland AB, Shi Q, McLendon RE, Bigner DD, Rich JN. Stem cell-like glioma cells promote tumor angiogenesis through vascular endothelial growth factor. Cancer Res. 2006;66(16):7843–8.

13. Folkins C, Shaked Y, Man S, Tang T, Lee CR, Zhu Z, Hoffman RM, Kerbel RS. Glioma tumor stem-like cells promote tumor angiogenesis and vasculogenesis via vascular endothelial growth factor and stromal-derived factor 1. Cancer Res. 2009;69(18):7243–51.

14. Algire GH, Chalkley HW, Earle WE, Legallais FY, Park HD, Shelton E, Schilling EL. Vascular reactions of normal and malignant tissues in vivo. iii. vascular reactions' of mice to fibroblasts treated in vitro with methylcholanthrene. J. Natl. Cancer. Inst. 1950;11(3):555–80.

15. Jain RK, Duda DG. Role of bone marrow-derived cells in tumor angiogenesis and treatment. Cancer Cell. 2003;3(6): 515–6.

16. Gerwins P, Sköldenberg E, Claesson-Welsh L. Function of fibroblast growth factors and vascular endothelial growth factors and their receptors in angiogenesis. Crit Rev Oncol Hematol. 2000;34(3):185–94.

17. Roskoski JR. Sunitinib: a vegf and pdgf receptor protein kinase and angiogenesis inhibitor. Biochem. Biophys. Res. Commun. 2007;356(2):323–8.

18. Kimura H, Braun RD, Ong ET, Hsu R, Secomb TW, Papahadjopoulos D, Hong K, Dewhirst MW. Fluctuations in red cell flux in tumor microvessels can lead to transient hypoxia and reoxygenation in tumor parenchyma. Cancer Res. 1996;56(23):5522–528.

19. Moeller BJ, Cao Y, Vujaskovic Z, Li CY, Haroon ZA, Dewhirst MW. The relationship between hypoxia and angiogenesis. Semin. Radiat. Oncol. 2004;14(3):215–21.

20. Dewhirst MW. Cao Y, Moeller B. Cycling hypoxia and free radicals regulate angiogenesis and radiotherapy response. Nat Rev Cancer. 2008;8(6):425–37.

21. Fox SB, Harris AL. Histological quantitation of tumour angiogenesis. APMIS. 2004;112(7-8):413–30.

22. Amit L, Ben-Aharon I, Vidal L, Leibovici L, Stemmer S. The impact of Bevacizumab (Avastin) on survival in metastatic solid tumors–a meta-analysis and systematic review. PLoS One. 2013;8(1):e51780.

23. Ricci-Vitiani L, Pallini R, Biffoni M, Todaro M, Invernici G, Cenci T, Maira G, Parati EA, Stassi G, Larocca LM, De Maria R. Tumour vascularization via endothelial differentiation of glioblastoma stem-like cells. Nature. 2010;468(7325):824–8.

24. Greenspan HP. On the growth and stability of cell cultures and solid tumors. J. Theoret. Biol. 1976;56(1):229–42.

25. Byrne HM, Chaplain MAJ. Growth of nonnecrotic tumors in the presence and absence of inhibitors. Math. Biosci. 1995;130(2):151–81.

26. Byrne HM, Chaplain MAJ. Growth of necrotic tumors in the presence and absence of inhibitors. Math. Biosci. 1996;135(2):187–216.

27. Byrne HM, Chaplain MAJ. Free boundary problems arising in models of tumour growth and development. Eur. J. Appl. Math. 1998;8:639–58.

28. Ribba B, Saut O, Colin T, Bresch D, Grenier E, Boissel JP. A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents. J. Theoret. Biol. 2006;243(4):532–41.

29. Araujo RP, McElwain DLS. A history of the study of solid tumour growth: The contribution of mathematical modelling. Bull. Math. Biol. 2004;66(5):1039–1091.

30. Byrne HM, Alarcón T, Owen MR, Webb SD, Maini PK. Modelling aspects of cancer dynamics: a review. Philos Trans A Math Phys Eng Sci. 2006;364(1843):1563–78.

31. Roosevelt T, Chapman S, Maini P. Mathematical models of avascular tumor growth. SIAM Review. 2007;49(2):179–208.

32. Balding D, McElwain DLS. A mathematical model of tumour-induced capillary growth. J. Theoret. Biol. 1985;114(1): 53–73.

33. Mantzaris NV, Webb S, Othmer HG. Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 2004;49(2):111–87.
34. Chaplain MAJ. The mathematical modelling of tumour angiogenesis and invasion. Acta Biotheoretica. 1995;43(4):387-402.

35. Orme ME, Chaplain MAJ. A mathematical model of vascular tumour growth and invasion. Math. Comput. Modelling. 1996;23(10):43-60.

36. Ribba B, Colin T, Schnell S. A multiscale mathematical model of cancer, and its use in analyzing irradiation therapies. Theor. Biol. Med. Model. 2006;3:7.

37. Ribba B, Marron K, Agur Z, Alarcón T, Maini PK. A mathematical model of Doxorubicin treatment efficacy for non-Hodgkin's lymphoma: investigation of the current protocol through theoretical modelling results. Bull. Math. Biol. 2005;67(1):79-99.

38. Ribba B, Saut O, Colin T, Bresch D, Grenier E, Boissel JP. A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents. J. Theoret. Biol. 2006;243(4):532-41.

39. Lignet F, Benzekry S, Wilson S, Billy F, Saut O, Tod M, You B, Berkane AA, Kassour S, Wei MX, Grenier E, Ribba B. Theoretical investigation of the efficacy of antiangiogenic drugs combined to chemotherapy in xenografted mice. J. Theoret. Biol. 2013;320:86-99.

40. Panovska J, Byrne HM, Maini PK. A theoretical study of the response of vascular tumours to different types of chemotherapy. Math. Comput. Modelling. 2008;47(5-6):560-79.

41. Billy F, Ribba B, Saut O, Morre-Trouilhet H, Colin T, Bresch D, Boissel JP, Grenier E, Flandrois JP. A pharmacologically based multiscale mathematical model of angiogenesis and its use in investigating the efficacy of a new cancer treatment strategy. J. Theoret. Biol. 2009;260(4):545-62.

42. Colin T, Cornelis F, Cumsille P, Lefebvre G, Poignard C, Saut O. Spatial modeling of tumor drug resistance: the case of gist liver metastases. Research report RR-8642 INRIA Bordeaux. 2014pp.26.

43. Cornelis F, Saut O, Cumsille P, Lombardi D, Iollo A, Palussiere J, Colin T. In vivo mathematical modeling of tumor growth from imaging data: Soon to come in the future? Diagn. Interv. Imaging. 2013;94(6):593-600.

44. Byrne H. Modelling avascular tumour growth In: Preziosi L, editor. Cancer Modelling and Simulation. Chapman & Hall/CRC Mathematical Biology and Medicine Series. Boca Raton, Fla, USA: Taylor & Francis; 2003. p. 75-120.

45. de Vries G, Hillen T, Lewis M, Müller J, Schönfisch B. A Course in Mathematical Biology: Quantitative Modeling with Mathematical and Computational Methods (Monographs on Mathematical Modeling and Computation). Philadelphia: Society for Industrial and Applied Mathematics; 2006.

46. Simeoni M, Magni P, Cammia C, De Nicolao G, Croci V, Pesenti E, Germani M, Poggesi I, Rocchetti M. Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after administration of anticancer agents. Cancer Res. 2004;64(3):1094-1101.

47. Benzekry S, Lamont C, Beheshti A, Tracz A, Ebos JML, Hlatky L, Hahnfeldt P. Classical mathematical models for description and prediction of experimental tumor growth. PLoS Comput. Biol. 2014;10(8):e1003800.

48. Sughrue ME, Rutkowski MJ, Aranda D, Barani IJ, McDermott MW, Parsa AT. Treatment decision making based on the published natural history and growth rate of small meningiomas. J. Neurosurg. 2010;113(5):1036-42.

49. Wanner O, Eberl H, Morgenroth E, Noguera D, Picioreanu C, Rittmann BE, van Loosdrecht MCM. Mathematical Modeling of Biofilms. Report of the IWA Biofilm Modeling Task Group, Scientific and Technical Report No. 18. London: IWA Publishing; 2006. See page 63.

50. Hogea C, Davatzikos C, Biros G. An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects. J. Math. Biol. 2008;56(6):793-825.

51. Colin T, Iollo A, Lombardi D, Saut O. System identification in tumor growth modeling using semi-empirical eigenfunctions. M3AS. 2012;22(06):1250003.

52. Goldacre R, Sylven B. On the access of blood-borne dyes to various tumour regions. Br. J. Cancer. 1962;16(2):306-22.

53. Boucher Y, Baxter L, Jain R. Interstitial pressure gradients in tissue-isolated and subcutaneous tumors: implications for therapy. Cancer Res. 1990;50(15):4478-84.

54. Ngo H, Tortorella SM, Ververis K, Karagiannis TC. The warburg effect: molecular aspects and therapeutic possibilities. Mol. Biol. Rep. 2015;42(4):825-34.

55. Suzuki C, Jacobsson H, Hatschek T, Torkzad MR, Bodén K, Eriksson-Alm Y, Berg E, Fujii H, Kubo A, Blomqvist L. Radiologic measurements of tumor response to treatment: Practical approaches and limitations. 1. Radiographics. 2008;28(2):329-44.

56. Luo W, Hu Q, Wang D, Deeb KK, Ma Y, Morrison CD, Liu S, Johnson CS, Trump DL. Isolation and genome-wide expression and methylation characterization of cd31+ cells from normal and malignant human prostate tissue. Oncotarget. 2013;4(9):1472-1483.

57. Miebach S, Grau S, Hummel V, Rieckmann P, Tonn JC, Goldbunner RH. Isolation and culture of microvascular endothelial cells from gliomas of different who grades. J. Neurooncol. 2006;76(1):39-48.

58. Acurio J, Troncoso F, Bertoglia P, Salomon C, Aguayo C, Sobrevia L, Escudero C. Potential role of A2B adenosine receptors on proliferation/migration of fetal endothelium derived from preeclamptic pregnancies. Biomed. Res. Int. 2014;2014:11.