## Abstract

Currently, most of the basic mechanisms governing tumour–immune system interactions, in combination with modulations of tumour-associated vasculature, are far from being completely understood. Here, we propose a mathematical model of vascularized tumour growth, where the main novelty is the modelling of the interplay between functional tumour vasculature and effector cell recruitment dynamics. Parameters are calibrated on the basis of different *in vivo* immunocompromised Rag1^{−/−} and wild-type (WT) BALB/c murine tumour growth experiments. The model analysis supports that tumour vasculature normalization can be a plausible and effective strategy to treat cancer when combined with appropriate immunostimulations. We find that improved levels of functional tumour vasculature, potentially mediated by normalization or stress alleviation strategies, can provide beneficial outcomes in terms of tumour burden reduction and growth control. Normalization of tumour blood vessels opens a therapeutic window of opportunity to augment the antitumour immune responses, as well as to reduce intratumoral immunosuppression and induced hypoxia due to vascular abnormalities. The potential success of normalizing tumour-associated vasculature closely depends on the effector cell recruitment dynamics and tumour sizes. Furthermore, an arbitrary increase in the initial effector cell concentration does not necessarily imply better tumour control. We evidence the existence of an optimal concentration range of effector cells for tumour shrinkage. Based on these findings, we suggest a theory-driven therapeutic proposal that optimally combines immuno- and vasomodulatory interventions.

## Major findings

We propose a tumour–effector cell recruitment model, where the main novelty is the low dimensional modelling of the complex interplay between the functional tumour-associated vasculature and the effector cell recruitment dynamics. Improved tumour vascularization, by means of either normalization or a microenvironmental stress alleviation strategy, is predicted to open a therapeutic window of opportunity in the treatment of cancer. Normalizing tumour vasculature as a potential therapeutic target closely depends on the immune recruitment dynamics and tumour sizes. Moreover, arbitrary low or high initial concentrations of effector cells can be detrimental to clinical outcomes, which evidences the existence of an optimal concentration range of effector cells that is crucial to obtain tumour control.

## 1. Introduction

Suppression of tumour angiogenesis has been recognized for more than four decades as a therapeutic target to treat solid tumours, as well as a complementary method of cancer prevention [1–3]. Traditionally, antiangiogenesis strategies attempt to inhibit the formation of new blood vessels within the tumour microenvironment, while at the same time reducing as much as possible the functionality of the existing tumour-associated vasculature. A wide array of drugs to inhibit tumour angiogenesis, such as bevacizumab (Avastin; Genentech/Roche), a ligand-trapping monoclonal antibody against vascular endothelial growth factor (VEGF), and the inhibitors of multiple protein kinase targets sorafenib (Nexavar, Bayer) and sunitinib (Sutent, Pfizer), have been successfully used in clinical practice to treat solid tumours [4]. However, different preclinical and clinical trials reveal that angiogenesis inhibitors alone have limited efficacy and only offer a modest survival benefit in a reduced cancer patient population [5,6]. The clinical benefits of antiangiogenesis treatments are not only limited in terms of tumour shrinkage, vasculature destruction and patient survival, but also are restricted by transient effects, insufficient efficacy or development of treatment resistance [2,7,8]. There is increasing evidence that the antitumour activity of most antiangiogenic drugs only becomes clinically significant in combination with conventional therapeutic modalities such as radiotherapy, chemotherapy or immunotherapy [9–11].

The tumour-associated vasculature exhibits opposing effects on tumour growth, invasion and progression [2,7,8]. On the one hand, blood vessels are required for delivery of oxygen and nutrients to the tumour, which promotes its growth, and they provide malignant cells with a way to spread, invade and metastasize into other healthy areas of the body [7,12–14]. On the other hand, blood vessels reoxygenate the tumour, which increases its radiosensitivity, and improve delivery of cytotoxic drugs to the tumour bulk as well as infiltration of immune cells into the tumour parenchyma [7,11,15,16]. In view of these opposing effects, the right treatment of a cancer patient relies on a mechanistic understanding and quantitative analysis of the cancer-mediated processes involved. This demands a mathematical model of tumour growth that considers the interplay between the functional tumour vasculature and the immune cell recruitment dynamics.

There are several mathematical models of tumour growth in which some form of immune system dynamics is included [17–21]. Different immuno-oncology studies reveal that tumours can survive in a microscopic and undetectable dormant state [19,20,22], as inferred from clinical data [23]. This equilibrium can be disrupted by sudden events affecting the immune system. Indeed, the neoplasm develops several strategies to circumvent the action of the immune system [24,25], which could result in tumour evasion and regrowth [25]. Sustained oscillations by the immune system have been observed both in its healthy state and in pathological situations [22,26,27]. An elegant mathematical model of the dynamics of immunogenic tumours revealed that tumour–immune system interactions can produce cyclic behaviours [28]. In fact, the presence of an immune component in mathematical models has been described as crucial for reproducing observed phenomena such as tumour dormancy and oscillatory growth [22,28–30]. This was complemented by a detailed discussion of the importance of including an immune component in tumour growth models (see [29] and references therein), and was subject to several reviews [17,31–35]. Some of the models mainly focus on spatial aspects described either by partial differential equations or by cellular automatas or by realistic mechanical cell interactions [31,33,36,37], whereas others focus on non-spatial dynamics described by ordinary differential equations (ODEs) [34,38]. Mathematical modelling provides a valuable theoretical framework to perform *in silico* experiments, as well as to verify assumptions and suggest modifications of existing theories. Despite years of research, and the vast amount of reported mathematical models, there are still many unanswered questions regarding the critical components of the tumour microenvironment and their role in the associated immuno-oncological mechanisms. However, models have been rather successful in optimizing therapeutic protocols based on quantitative incorporation of experimental or patient data [29,39–45]. The work presented here follows this latter approach.

We intend to generate novel mechanistic insights not only into the interactions between the immune system and growing tumours, but also into the role of the functional tumour-associated vasculature. The main goal is to use this better understanding to suggest optimized treatment strategies. For that purpose, we combined a mathematical model of radially symmetric tumour growth [46–49] with an immune cell model [28] to propose a tumour–effector cell recruitment model. The main novelty is that the impact of the functional tumour vasculature on the effector cell recruitment dynamics is included. Previous models of tumour–immune system interactions did not directly incorporate the effect of tumour vascularization, and thus underestimated the relevance of this crucial mechanism for tumour growth and treatment outcomes. Based on well-supported biological assumptions, we derive that tumour vasculature normalization opens a therapeutic window of opportunity for tumour growth control. The potential benefit of a vascular normalization strategy relies on a specific immune recruitment dynamic into the tumour microenvironment. We also find that an unlimited increase in the initial effector cell concentration does not necessarily imply better tumour control. The proposed model is capable of providing a more complete view of vascularization effects on the immune responses to tumour growth, and it is therefore suited for analysis and prediction of treatment strategies.

## 2. Quick guide to model equations and assumptions

### 2.1. A radial tumour growth model

We consider a mathematical model of radial tumour growth which was originally proposed in [46]; see also [48,49]. The tumour is assumed to be an incompressible fluid flowing through a porous medium, where tissue elasticity is simplified. The tumour–host interface is considered to be sharp, and cell-to-cell adhesive forces are modelled as a surface tension at that interface. The tumour expands as a mass whose growth is governed by a balance between cell birth (mitosis) and death (apoptosis and necrosis). The mitotic rate within the tumour is assumed to be linearly dependent on the nutrient concentration (oxygen, glucose, etc.) and is characterized by its maximal value *λ*_{M} at the tumour–host interface. The death rate *λ*_{A} is considered uniform within the tumour and assumed constant in time. The concentration of nutrients obeys a reaction–diffusion equation in the tumour volume, where nutrients are supplied from the functional tumour vasculature and consumed by the tumour cells at a uniform consumption rate (see the derivation of the mathematical model provided in the electronic supplementary material for further details).

The resulting dimensional mathematical model for the tumour radius *R*(*t*) dynamics, under the assumption of radial symmetry, is given by
2.1

The model given by equation (2.1) can be interpreted as the temporal evolution of the average tumour radius, because radial symmetric growth is not a realistic behaviour. Note that multiplication of equation (2.1) by *R*^{2} describes the time evolution of the tumour volume (see the electronic supplementary material for more details). The non-negative and dimensionless parameter *B* represents the net effect of the tumour-associated vasculature on the tumour radius dynamic. Moreover, *L*_{D} is an intrinsic scale representing the average length of nutrient gradient, i.e. supply, diffusion and consumption.

### 2.2. A tumour–effector cell recruitment model considering the impact of the functional tumour-associated vasculature

The proposed tumour–effector cell recruitment model is a combination of the radial tumour growth model given by equation (2.1) and an effector cell model initially introduced in [28]; see also [20,50]. In addition, we consider the impact of the functional tumour vasculature on the tumour growth and immune response dynamics. The system variables are the average tumour radius *R*(*t*) and the effector cell concentration *E*(*t*) in the tumour vicinity. The model is formulated as a system of ODEs given by
2.2and
2.3where the time coordinate *t* on the system variables is omitted for notational simplicity. The non-dimensionalized version of the model is provided in the electronic supplementary material. In words, equations (2.2) and (2.3) can be respectively read as follows:
The expression for is a phenomenological scaling function that models the infiltration of effector cells into the tumour bulk through the existing functional tumour vasculature. This function modulates the terms related to tumour–effector cell interactions, such as killing of tumour cells owing to effectors and inactivation of effectors owing to their antitumour activity. As a particular case, we assume that *α* = *B*, where represents the net effect of functional levels of the tumour-associated vasculature. In the limit of avascular tumour growth *B* = 0, such tumour–effector cell interactions take place only at the tumour surface. At the other extreme, for *B* = 1, effector cells can potentially interact with any cancer cell within the tumour bulk. Note that such scaling functions have also been considered in the classical von Bertalanffy approach and more recently in allometric models [51]. Alternatively, we can view the encountering of effector and tumour cells mediated by fractal vasculature as a percolation process within the tumour mass. In this context, the term *f*(*R*, *α*) scales the average effector–tumour cell encountering rate in a vascularized tumour. In particular, the scaling term should be proportional to the corresponding percolation cluster mass of system size *R* and fractal dimension 3*α*, where the latter depends on the structure of the tumour-associated vasculature [52].

### 2.3. Model assumptions

(1) The temporal evolution of the average tumour radius is considered, where invasive and diffusive tumour mechanisms are not taken into account.

(2) The death rate

*λ*_{A}of tumour cells reflects the lump effect of apoptotic and necrotic processes.(3) Innate immunity or base immune surveillance is represented as a minimum presence of active effector cells at any time, even in the absence of cancer cells.

(4) Effector cells are recruited depending on the concentration of tumour cells according to Michaelis–Menten dynamics.

(5) The efficacy of immune killing depends on the ability of effector cells to penetrate the tumour bulk via the functional tumour-associated vasculature [10,11]. With improved vascularization, the effectors kill tumour cells not only on the surface of the tumour, but also further inside.

(6) Effector cells die at a constant rate and are inactivated in accordance with their antitumour activity.

## 3. Material and methods

### 3.1. Experimental data of tumour growth in mice

We now describe the murine experiments of the *in vivo* tumour growth considered. BALB/c mice were purchased (Janvier, France). Rag1^{−/−} mice were purchased from The Jackson Laboratory (C.129S7(B6)-Rag1^{tm1Mom}/J). Ig*α*^{−/−} mice were obtained from Michael Reth, Freiburg, Germany (Cd79a^{tm2Pld}). All recombinant mice were bred at the Helmholtz Center for Infection Research (HZI), Braunschweig, Germany, and all experiments were performed with 8–12 week old female mice, unless otherwise stated, and under approval of LAVES (Niedersächsisches Landesamt für Verbraucherschutz und Lebensmittelsicherheit); permission: 33.9-42502-04-12/0173. Unlike wild-type (WT) BALB/c mice, immunocompromised Rag1^{−/−} mice cannot produce functional T and B cells, i.e. they have no adaptive immune responses.

CT26 (ATCC CRL-2638) cells were cultured in Iscove's modified Dulbecco's medium supplemented with 10% fetal calf serum, 100 U ml^{−1} penicillin, 100 µg ml^{−1} streptomycin, 50 µM 2-mercaptoethanol and 2 µM l-glutamine and were maintained at 37°C and 5% CO_{2}. The F1A11 (H-2* ^{d}*) cell line is a murine fibrosarcoma that expresses β-galactosidase (β-gal) and was obtained by transduction of the spontaneously transformed BALB/c fibroblast cell line F1 with the LBSN retroviral vector [53]. Tumours CT26 and F1A11 were set by injecting 5 × 10

^{5}cells in 100 µl phosphate-buffered saline subcutaneously (s.c.). Growth was monitored by caliper, and the mean tumour radius evolution between days 4 and 9 was observed. Volume was calculated as where

*h*= height and

*w*= width.

### 3.2. Estimation of model parameters

The experimental data corresponding to immunocompromised Rag1^{−/−} and WT BALB/c mice were analysed in two different phases. In the first phase (days 0–4), exponential avascular tumour growth was assumed [54]. Using the tumour radius at these early times of growth, the net proliferation rate of cancer cells was estimated. The tumour radius evolutions in the second phase of growth (days 4–9) were then used to determine the remaining model parameter values. Because of uncertainties owing to the finite number of observations, consisting in 15 and 20 individual experimental trials, respectively, we implemented a bootstrapping procedure for each dataset considered. The goal was to obtain reliable estimates of the summary statistics, where a resampling procedure was performed. A set of bootstrap samples was artificially generated and represented the natural variations in the experimental datasets [55]. More precisely, we considered non-parametric bootstrap samples, i.e. samples that had been generated from the original experimental datasets by drawing with replacements. The default size managed, also called the cloud size, was 100 000 bootstrap samples per dataset, where the resulting mean and standard deviation of such samples was considered for parameter calibration. Another crucial factor to obtain precise model fitting is to consider different initial parameter values. Consequently, we started the fitting procedure from a large number of different random initial conditions, as well as by randomly perturbing the model parameter set to be estimated. Then, the solution with the lowest residual variance was selected, which thus provides the best fit.

### 3.3. Early exponential tumour growth phase (days 0–4)

For a precise understanding of the early tumour growth phase predicted by our mathematical model, we used the small tumour radius assumption and Taylor expanded the r.h.s. of equation (2.1) accordingly. Keeping the terms up to the first nonlinear term, we obtain the following approximation: 3.1

The very early tumour growth phase, dictated by the first term of this approximation, is always of an exponential nature and does not depend on the parameter *B*, i.e. the tumour vascularization effects, as reported also in [48,49,54]. Therefore, this initial phase depends exclusively on the net proliferation rate of cancer cells, Considering the initial tumour radius in the immunocompromised Rag1^{−/−} mouse experiments of 0.40 mm, which corresponds to a tumour spheroid of 5 × 10^{5} cells of about 10 µm in diameter, and the tumour radius at day 4 of 2.0 mm (figure 1*a*), we find that . This value is in agreement with previous estimates of growth rate (early stages) for murine CT26 spheroids [56]. Later, the initial exponential tumour growth is affected by tumour-associated vascularization, which starts playing a significant role when the two terms of the r.h.s. of equation (3.1) have similar orders of magnitude. This occurs around a critical tumour radius *R _{B}* derived from the condition
3.2

Because is necessary for the initial tumour growth, equation (3.2) requires where *B* = 0 represents an avascular tumour and *B* → 1 a fully vascularized tumour. This yields
3.3

Equation (3.3) provides an estimate of the tumour size at which nonlinearities, involving the vascularization mechanisms, start to significantly influence the initial exponential tumour growth. The effect of the nonlinear terms is to saturate the growth, i.e. d*R*/d*t* = 0, when the tumour radius reaches *R _{B}*. Because only exponential growth is observed for small tumour radii, angiogenesis should be activated before the critical tumour radius is reached and suggests that

*R*is an upper bound for the initiation of angiogenic processes. The value of

_{B}*R*can be predicted provided the values of the parameters in equation (3.3) are known. Assuming time-invariant values of the characteristic mitotic

_{B}*λ*

_{M}and death

*λ*

_{A}rates of cancer cells, and the relation together with the physiological plausible value for the mouse colon carcinoma cell line (or CT26 murine tumour cells) [56,57], yields The parameters

*B*and

*L*

_{D}, i.e. the functional tumour vasculature and intrinsic scale representing the average length of nutrient gradient, remain to be determined; these are estimated from the growth phase ranging between days 4 and 9.

### 3.4. Linear tumour growth phase (days 4–9)

From the experimental data corresponding to immunocompromised Rag1^{−/−} and (WT) BALB/c mice at the stages of growth ranging between days 4 and 9, we observe an approximately linear evolution of the tumour radius (figure 1). This observation is consistent with the long-term behaviour of growth, where linear radial evolution is prominent, as reported in [54]. In equations (2.2) and (2.3), this assumption is equivalent to assigning *B* = *λ*_{A}/*λ*_{M}. In particular, linear radial tumour growth in immunocompromised Rag1^{−/−} mice requires *c* = 0 in equation (2.2), because Rag1^{−/−} mice have no adaptive anti-cancer immune responses. Thus, the first and last terms on the r.h.s. of equation (2.2) vanish and this allows us to determine This value is consistent with the well-known characteristic nutrient diffusion length that ranges between 0.2 and 0.3 mm [58,59]. In figure 1*a*, the linear fitting results for the Rag1^{−/−} mouse dataset are represented. The estimated parameter values allow for the calculation of >*L*_{D} from equation (3.3) for a tumour in the early phase of growth which begins to be influenced by the tumour-associated vasculature. The following parameter values are taken from [20,28,50,60,61]: *d*_{0} = 0.37 d^{−1}, *d*_{1} = 0.01 mm^{−3} d^{−1}, *K* = 2.72 mm^{3} and *σ* = 0.13 × 10^{5} cells d^{−1}.

Radial tumour growth can also be assumed linear for the (WT) BALB/c mouse dataset; from figure 1*b* we can infer ^{−1} constantly in the time interval between days 4 and 9, and in agreement with the growth velocity measured previously for murine CT26 tumours [62]. Then, for equation (2.2) becomes
3.4where and for it follows that because Considering that equation (3.4) provides the following evolution of the effector cell concentration *E*(*t*) required for linear radial tumour growth
3.5

Formula (3.5) allows for the exact quantification of the initial concentration of effector cells *E*_{0}, as well as for the reduction in the number of free model parameters that enter in the fitting procedure, because the net proliferation *λ*_{p} and death *λ*_{A} rates of cancer cells, as well as *L*_{D}, have been estimated from the experimental data of immunocompromised Rag1^{−/−} mouse tumour growth. Accordingly, only three free model parameters are left in equations (2.2) and (2.3), i.e. the functional tumour vasculature *B*, immune cell recruitment rate *r* and the death rate of tumour cells owing to effectors *c*, which are estimated from the (WT) BALB/c mouse experiments imposing the expression (3.5) during parameter calibration. Considering variations in the initial conditions, and randomly perturbing the parameter set to be estimated, the resulting parameter values are provided in table 1. In figure 1*b*, the linear fitting results for the experimental data of tumour growth in (WT) BALB/c mice are represented.

## 4. Results

### 4.1. The dual effect of functional tumour-associated vasculature

The proposed mathematical model given by equations (2.2) and (2.3) is used to investigate the effect of the functional vasculature on the gross tumour growth. Tumour-associated vascularization enhances the supply of oxygen and nutrients, and therefore favours tumour growth and progression [12,14]. At the same time, blood vessels facilitate the infiltration of effector cells into the tumour bulk, which potentially exterminate cancer cells [15,16]. These opposing effects of tumour vascularization suggest the existence of an optimal tumour growth control, such that a quantitative modulation of functional vasculature can be considered as a suitable therapeutic target against cancer.

In order to gain insight into the impact of the functional vasculature on the short-term tumour evolution, we vary the model parameter *B* in equation (2.2). The radial growth rate d*R*/d*t*, which indicates whether the tumour radius grows (positive), decays (negative) or remains stable (zero) is monitored. Starting from different initial concentrations of effector cells *E*_{0}, figure 2 shows the rate of tumour growth in relation to the functional levels of vascularization *B* and initial tumour sizes *R*_{0}. Above a specific threshold concentration of effector cells, a degree of vascularization exists that maximizes tumour growth. Moreover, above this threshold increasing vascular functionality reduces the tumour growth rate, because it allows a high number of effector cells to penetrate into the tumour parenchyma; see diagrams for cells in figure 2. Tumour shrinkage d*R*/d*t* < 0 is possible for very low functional vascularization and large tumour radii, irrespective of the number of effector cells in the tumour vicinity. However, for small tumours and low initial concentrations of effector cells, no reduction in the tumour burden can be obtained for deficient vascularization. Moreover, there exist critical *E*_{0} values at which tumour shrinkage also becomes possible for well-vascularized tumours, irrespective of their sizes; see diagram for *E*_{0} = 30 × 10^{5} cells in figure 2.

A more systematic investigation into the dependence of the long-term model dynamics given by equations (2.2) and (2.3) on the set of initial conditions (*E*_{0}, *R*_{0}) is conducted through a tumour–immune phase portrait analysis (figure 3). For functional vascular levels of *B* = 0.50, long-term tumour control is possible, in the sense of damped oscillations, for all values of *R*_{0} considered (figure 3*a*). Tumour dormancy is described by decaying oscillations converging to a steady state that, once obtained, is maintained for a long period of time. However, at larger values of *R*_{0}, the amplitude of the oscillations in *E*(*t*) and *R*(*t*) increases and reaches uncontrolled growth above a threshold for *R*_{0} which determines the tumour fate. For a higher tumour-associated vascularization, e.g. *B* = 0.70, long-term tumour control cannot be achieved for any set of values of the model parameters considered (figure 3*b*). In contrast, for well-vascularized tumours, e.g. *B* = 0.95, a window of opportunity for tumour control is predicted (figure 3*c*). As long as *R*_{0} stays below a threshold value, tumour evolution can be controlled. At high and low vascularization levels, tumour eradication is not possible for sufficiently large tumours, revealing the relation between tumour size and the potential therapeutic benefit of targeting the tumour vascular network. A bifurcation analysis on model parameters *B* and *r*, i.e. the functional tumour vasculature and immune recruitment rate, is provided in the electronic supplementary material.

### 4.2. Critical radius and control range of effector cell concentration determine tumour fate

Figure 3 reveals information about the existence of different dynamical regimes of the tumour–effector cell recruitment system described by equations (2.2) and (2.3) for the set of parameter values considered. Based on the analysis of the critical fixed points and their positions in the phase space (see the phase space and critical fixed points analysis provided in the electronic supplementary material), we obtain that a spiral and a saddle point are present. At the spiral point, the tumour radius oscillates around a small tumour size, i.e. a dormancy-associated equilibrium point. This implies a concentration range of effector cells that keep the tumour under control for a period of time that depends on the point stability. Moreover, there exists a critical tumour radius *R*_{S}, which coincides with the saddle point, above which tumours of equal or larger sizes always grow uncontrollably and the immune system is completely suppressed. However, below this threshold, tumours are always controlled, or at least for some period of time in a transient dormant state, provided the immune responses are in a specific range. The questions that arise are to (i) analytically evaluate the critical tumour size *R*_{S} that for implies uncontrollable growth and (ii) estimate the range of initial effector cell concentration *E*_{0} that allows for long-term tumour control.

#### 4.2.1. Critical radius for tumour dormancy

These results demonstrate that the saddle point position in the phase space separates uncontrolled tumour growth and immune-induced dormant states, as well as delineates the potential and limits of therapeutic benefits of immuno- and vasomodulatory interventions. For large well-vascularized tumours, we derive an analytical expression that approximates the critical tumour radius *R*_{S} corresponding to an upper bound estimation of the saddle point (see the estimation of the critical radius for tumour dormancy provided in the electronic supplementary material for further details), which is given by
4.1

We demonstrate that, for values of the tumour radius equal to or higher than *R*_{S}, the eigenvalues corresponding to the radial time evolution are positive. This implies that, for radii beyond *R*_{S}, the tumour grows uncontrollably, evading immune responses. The nonlinear system of equations (2.2) and (2.3) was numerically solved using the trust-region dogleg method [63]. Figure 4*a* shows the dependence of *R*_{S} on the immune recruitment rate *r* obtained by means of model simulations. *R*_{S} represents an upper bound of the critical tumour size that determines the long-term tumour fate, and as *r* increases the simulated and estimated values get closer. For *R* < *R*_{S}, long-term or at least transient tumour dormancy can be obtained, whereas results in tumour escape from the immune system surveillance. In the situation of tumours with low functional vascular levels, *R*_{S} can be directly estimated by means of model simulations.

#### 4.2.2. Control range of effector cell concentration regulates tumour evolution

The efficacy of infiltrating effector cells to control the tumour growth is derived from the properties of the spiral point, which is always obtained for tumours with size less than the critical threshold *R*_{S}. Tumour dormancy, and then neoplasia control, takes place when the tumour radius stays close to the spiral point for a long period of time. Therefore, the attraction range of the spiral point with respect to the *E*-axis provides useful information about the number of effector cells required for long-term tumour control.

We claim that, for arbitrary high or low initial concentrations of effector cells *E*_{0}, immune-induced tumour dormancy is not guaranteed. To illustrate the counterintuitive loss of tumour control for high or low initial numbers of effector cells, we focus on the phase portrait diagram in figure 3*b*. There we observe that, when the amplitude of the *E*(*t*) oscillations exceeds a certain range, the tumour escapes from the immune system surveillance. Mathematically, this is translated to the existence of a concentration range of effector cells Δ*E* that is crucial for the control of tumours with size smaller than the critical radius *R*_{S}. Details of the Δ*E* calculation are provided in the electronic supplementary material.

Figure 4*b* shows the dependence of the control range of effector cell concentration Δ*E* on the functional tumour-associated vasculature *B* and immune recruitment rate *r*. In particular, we observe that Δ*E* increases for high enough immune recruitment rates and increasing tumour vascularization. However, at low values of Δ*E* becomes smaller as *B* increases. The latter is intuitive considering that low immune recruitment rates could not provide enough effector cells to exploit the enhanced potential infiltration offered by increasing functional tumour vasculature. Therefore, the pro-tumoral effect of the functional vascular network overcomes its antitumoral action.

### 4.3. Therapeutic potential of immuno- and vasomodulatory interventions

Immune-mediated tumour dormancy depends not only on the effectiveness of the immune system responses, but also on the functional tumour-associated vasculature. In consequence, we investigate the impact of model parameters *B* and *r*, i.e. the functional tumour vasculature and immune recruitment rate, on the overall therapeutic potential for long-term (5 years) model predictions. The therapeutic potential is interpreted as the average tumour control, i.e. the tumour radius stays bounded, with respect to a certain regime of model parameter given by
4.2where *t** = 1825 days (5 years) is the target time, *R*_{S}(*p*) is the tumour radius at the saddle point and *H*(·) is a Heaviside step function. The effects of the initial tumour radius *R*_{0} and effector cell concentration *E*_{0} have been also investigated. It should be noted that a patient who has been tumour-free for 5 years is typically considered to be in a complete remission stage, i.e. permanent absence of disease, and thus cured.

Any increase in the initial tumour radius *R*_{0} always reduces the therapeutic potential of the immune system to control tumour growth (figure 5). However, an optimal concentration of effector cells for long-term tumour control is predicted, irrespective of the degree of functional tumour vascularization *B*. This observation implies that an arbitrary increase in the initial concentration of effector cells does not necessarily result in a greater therapeutic benefit. It is counterintuitive that high values of *E*_{0} do not enhance long-term tumour control. A large initial number of effector cells can induce a great initial tumour burden reduction that results in the tumour escaping immune control by an undercritical immune stimulus. In mathematical terms, the tumour radius time evolution crosses the separatrix, escapes from the limit cycle and enters into the unlimited growth dynamic regime (see for instance the green curve in the phase portrait diagram of figure 3*b*, as well as the phase space and critical fixed points analysis provided in the electronic supplementary material).

In addition to the initial concentration of effector cells, the immune recruitment dynamics play an important role in long-term tumour control (figure 6). Large values of the immune recruitment rate *r* enhance the tumour control, irrespective of the tumour vascularization regime *B*. However, for high levels of functional tumour vasculature, i.e. *B* > 0.55, large tumours with *R*_{0} > 4 mm are almost uncontrollable at any immune recruitment rate *r* (figure 6*c,d*). We also demonstrate that the initial concentration of effector cells required for long-term tumour control mostly depends on the functional tumour vascularization (see the phase space and critical fixed points analysis provided in the electronic supplementary material).

### 4.4. A theory-driven therapeutic proposal

Model predictions provide new insights that can be used to suggest and design novel therapeutic strategies against tumour growth. The aim is to propose alternative or complementary ideas for cancer treatments that are less toxic than conventional therapeutic modalities, such as surgery, radiotherapy or chemotherapy. The latter cancer therapies are related to immunosuppressive effects that disallow the intention to exploit immune system responses.

An appropriate tumour profiling via medical imaging techniques and tissue biopsy is essential. Imaging allows for the estimation of the initial tumour size *R*_{0} and its comparison with the critical radius *R*_{S} that potentially determines the long-term tumour fate. For calculation of the precise *R*_{S} value, the net proliferation rate of tumour cells could be extrapolated from biopsy samples by following the methodology proposed in [64]. Moreover, from the same pathology data the initial concentration of effector cells *E*_{0} can be estimated, e.g. by counting the number of CD8^{+} antitumour effector T cells. Then, the range of the effector recruitment rate could be assumed to be comparable to the one estimated from our experimental data (table 1). For *R*_{0} < *R*_{S}, mid- or long-term-controlled tumour growth is possible by clinically influencing the interplay between the immune recruitment dynamics and the functional tumour-associated vasculature, as suggested by the therapeutic potential results in figures 5 and 6. In the situation of tumour control could not be obtained according to the findings in the proposed model, and complementary treatment modalities are initially required.

Positive therapeutic responses, d*R*/d*t* < 0, can be achieved by increasing or decreasing the functional tumour vasculature. The latter is possible via antiangiogenic drugs, such as bevacizumab, sorafenib or sunitinib. However, low levels of tumour vascularization trigger hypoxia-induced phenomena, and this treatment choice typically fails after some short-term tumour burden reduction benefits [5,6]. Hypoxic tumour regions promote either the secretion of alternative angiogenic factors, e.g. VEGF [65], or enhanced tumour cell migration [66], or suppressed immunoactivity allowing tumour cells to evade host immunosurveillance [10,11]. Moreover, hypoxia is well known to reduce tumour cell sensitivity to radiation and chemotherapy [67]. To avoid these adverse effects, and following our model predictions, an alternative strategy favouring high functional tumour vasculature could be adopted, i.e. by increasing the model parameter *B*. This can be realized via normalizing the tumour vasculature [2,7,8] or by means of a microenvironmental stress alleviation, i.e. decompression of tumour vessels (e.g. [68]). Improved functional tumour vascularization allows for a wider range of effector cell types to infiltrate the tumour bulk and further exterminate cancer cells [15,16]. This is supported by several preclinical studies demonstrating that vascular normalization increases infiltration of immune effector cells into tumours and substantially improves survival [69–71]. Indeed, vascular normalization improves tumour blood vessel perfusion, reduces tumour hypoxia, and enhances the delivery of cytotoxics and radiotherapy efficacy [11]. Counterintuitively, normalization of the tumour vasculature mitigates the metastatic risk because it is correlated to diminishing the epithelial–mesenchymal transition of cancer cells [72].

The model predictions shown in figure 5 describe variations in the resulting therapeutic window of opportunity in terms of long-term tumour control for different functional levels of tumour-associated vasculature *B*, the initial number of effector cells *E*_{0} and tumour sizes *R*_{0}. A negative correlation between initial tumour sizes and mid- or long-term control is predicted. For sufficiently small tumours, a range of optimal *E*_{0} values with higher therapeutic potentials exist. An increase in the immune recruitment rate results in a higher therapeutic potential and the window of opportunity for high functional levels of tumour vasculature still persists (figure 6). Therefore, in addition to tumour vasculature normalization, applying a monitored immunotherapy protocol could further increase the expected therapeutic benefits. Figure 7 describes the complete theory-driven therapeutic proposal.

## 5. Conclusion

In this work, we have explored the possibilities offered by mathematical modelling and computer simulations to gain insights into the interplay between the functional tumour-associated vasculature and immune recruitment dynamics. We have combined a radially symmetric tumour growth and an effector cell model, including the impact of vascularization on the immune system efficacy in killing tumour cells. The model parameter calibration was based on murine experiments of *in vivo* tumour growth. Model analysis revealed the existence of a tumour control window of opportunity when high functional levels of tumour vasculature are present. This can be biologically understood, because high vascularization allows for the penetration of effector cells into the tumour bulk and diminishes the adverse effects of hypoxia caused by low oxygenation. Improving tumour vascularization can be mediated by a normalization of blood vessel functionality [2,7,8,67,68]. Such an approach is successful only if the tumour size is below a critical threshold that depends on tumour cell metabolic needs. Furthermore, we predict that an arbitrary increase in the initial effector cell concentration does not necessarily imply a direct enhancement of tumour control, thus optimal values are needed. According to our model predictions when the initial reduction in tumour size exceeds a specific threshold, the tumour escapes immune control and relapses faster than the associated tumour-induced immune responses. Taking into account the dependency of tumour control on these various mechanisms, we have proposed an alternative less toxic therapeutic approach based on personalized tumour profiling. This therapeutic proposal, based on normalizing the tumour vasculature, can also be potentially combined with chemotherapy and radiotherapy to enhance the therapeutic outcomes [2,7,67].

We conclude by pointing out a number of related future research directions, as well as the limitations of this work. The current results are based on murine experiments and therefore a parametrization on human data would be essential. Introducing tumour vascularization dynamics should be not only a natural, but also an insightful extension of the proposed mathematical model. Moreover, we need to further refine the modelling of the term *f*(*R*, *α*) and elaborate the dependence of the exponent *α* on the functional tumour vasculature *B*. Additionally, we could extend the notion of parameter *α* to reflect the dynamic feedback of effector cell activity on local tumour density fluctuations as a time-evolving fractal dimension. In this state, vascularization *B* is considered to be a static parameter in time. Making tumour vasculature dynamic would make sense to model further hypoxic effects, such as necrosis- [73] or hypoxia-induced invasion [48]. Moreover, this will make it possible to further investigate the dual effects of functional blood vessels on tumour growth dynamics when cytotoxic drugs are administered as part of chemotherapy. On the other hand, the immune system is much more complex than its description in the current model, involving a wide range of immune cell types and intricate interactions. In particular, the immune system is often regarded as acting to suppress tumour growth; however, it is now clear that it can be both stimulatory and inhibitory [74]. Recent advances have indicated that tumour-associated macrophages actively promote all aspects of tumour initiation, progression and development [75]. Therefore, modelling the interplay between macrophages and tumour cells, while taking into account the functional tumour-associated vasculature and effector cell recruitment dynamics, may highlight new targets to develop novel anti-cancer strategies.

## Authors' contributions

H.H. conceived of and developed the methodology; C.S. and S.W. conceived of and designed the experiments; H.H., J.C.L.A., S.M. and M.M.H. analysed and interpreted the data and model results; H.H. and J.C.L.A. provided administrative, technical and material support; S.W. and M.M.H. supervised the study; H.H., J.C.L.A., S.W. and M.M.H. wrote and revised the manuscript. All authors read and approved the final manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

H.H. and J.C.L.A. gratefully acknowledge the support of the German Excellence Initiative via the Cluster of Excellence EXC 1056 Center for Advancing Electronics Dresden (cfAED) at the Technische Universität Dresden and the German Federal Ministry of Education and Research (BMBF) for the eMED project SYSIMIT (01ZX1308D). H.H. and M.M.-H. were supported by the Human Frontier Science Programme. M.M.-H. was supported by Measures for the Establishment of Systems Medicine projects SYSIMIT and Sys-Stomach (01ZX1310C), by the Federal Ministry of Education and Research, Germany, and by iMed—the Helmholtz Initiative on Personalized Medicine. C.S. and S.W. gratefully acknowledge the support of the Deutsche Krebshilfe.

## Acknowledgements

We thank Professor Miguel A. Herrero and Professor Andreas Deutsch for the fruitful discussions.

- Received May 15, 2015.
- Accepted October 6, 2015.

- © 2015 The Author(s)