Louise Viger
email: louise.viger@coria.fr

Fabrice Denis

Martin Rosalie

Christophe Letellier

A cancer model for the angiogenic switch

Keywords: Cancer model, angiogenic switch, endothelial cells, Chaos

The occurrence of metastasis is an important feature in cancer development. In order to have a one-site model taking into account the interactions between host, effector immune and tumor cells which is not only valid for the early stages of tumor growth, we developed in this paper a new model where are incorporated interactions of these three cell populations with endothelial cells. These latter cells are responsible for the neo-vascularization of the tumor site which allows the migration of tumor cells to distant sites. It is then shown that, for some parameter values, the resulting model for the four cell populations reproduces the angiogenic switch, that is, the transition from avascular to vascular tumor.

Introduction

Tumor growth is a complex process depending on various cell types as mutant (tumor) cells, host (normal-tissue) cells, immune cells (lymphocytes, macrophages), endothelial cells... In order to get real insights into critical parameters that control system dynamics, theoretical models are required [1]. They can also be used for designing new effective treatments without an extensive experimentation [START_REF] Byrne | Using mathematics to study solid tumour growth[END_REF]. Typically, there are three levels for describing interactions among these cells, namely the cell-, tissue-and organ-level. At the individual-cell level, the complexity of the model increases in proportion to the number of specific cells which contribute to the tumor growth. These models are well suited for investigating in details a particular mechanism arising in a specific cancer. For instance, such an approach was used to describe the role of the Adenosine TriPhosphate (ATP) metabolism as the cellular energy carrier in tumor angiogenesis [START_REF] Nagy | Evolution of uncontrolled proliferation and the angiogenic switch in cancer[END_REF]: the corresponding cell energetics model is made of three ordinary differential equations. The ATP is thus used to drive a second model -at the tissue level -for tumor growth made of three other ordinary differential equations governing tumor mass, immature vascular endothelial cell mass and total microvessel length; the model depends on twenty-three parameters. Such a model type is necessarily specific to a given type of cancer since based on a description of interactions at the cell level. In order to have a generic model, the tissue level was here retained for modelling the tumor growth. This level of description was required to reduce the model complexity, thus avoiding a model too difficult to parametrize and to investigate [START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF].

Indeed, most of the models used to investigate tumor growth are at the tissue level as reviewed by Arujo and McElwain [START_REF] Araujo | A history of the study of solid tumour growth: the contribution of mathematical modelling[END_REF] or by Eftimie and coworkers [START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF]: at this level of description unavoidable simplifications are made, focusing on some detailed mechanisms depending on the objectives. Most of the models are devoted to the interactions between tumor cells and immune cells [START_REF] Eftimie | Interactions between the immune system and cancer: a brief review of non-spatial mathematical models[END_REF]. For a generic description of interactions between tumor and immune cells, the number of ordinary differential equations is at least two [START_REF] Boon | Human tumor antigens recognized by T lymphocytes[END_REF][START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF][START_REF] Onofrio | Metamodeling tumor-immune system interaction, tumor evasion and immunotherapy[END_REF]. Some specifications are developed for designing therapies in particular contexts for which specific cell types can be taken into account, thus increasing the model dimension [START_REF] Byrne | Macrophage-tumor interactions: in vivo dynamics[END_REF][START_REF] Berner | IFN-γ mediates CD4 + T-cell loss and impairs secondary antitumor responses after successful initial immunotherapy[END_REF][START_REF] Bunimovich-Mendrazitsky | Mathematical model of BCG immunotherapy in superficial bladder cancer[END_REF]. Thus, a set of 11 ordinary differential equations can be obtained [START_REF] De Boer | Macrophage T lymphocyte interactions in the anti-tumor immune response: a mathematical model[END_REF]. In these cases, models are specific to a given type of cancer. One of them was for instance developed for malignant melanoma [START_REF] Eikenberry | Tumor-immune interaction, surgical treatment, and cancer recurrence in a mathematical model of melanoma[END_REF] by taking into account tumor cells, healthy cells, tumor angiogenic factor, blood vessel endothelial cells, necrotic debris, spatial pressure of oxygen and basement membrane: this model is made of partial differential equations for describing the spatial tumor growth as well as tumor angiogenesis. There are also some models combining a tissue-level description with some migra-tions among various organs to take into account the circulating endothelial progenitor cells [START_REF] Stoll | A mathematical model of the contribution of endothelial progenitor cells to angiogenesis in tumors: implications for antiangiogenic therapy[END_REF] since these cells may contribute to tumor angiogenesis [START_REF] Rafii | Vascular and haematopoietic stem cells: novel targets for anti-angiogenesis therapy?[END_REF].

In our cases, we would like to develop a generic model (not specific to a given type of cancer) at the tissue level for a single tumor site, that is, ignoring for now the diffusion process. Our aim is to include the micro-environment (host or healthy cells) of the tumor cells as considered in [START_REF] Owen | Modelling the macrophage invasion of tumours: effects on growth and composition[END_REF] or in [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumor: modeling, applications and biological interpretations[END_REF]. To overcome the limitation of these latter models which are only for avascular tumor growth, we introduced interactions between tumor, immune, host and endothelial cells as recommended by Merlo and co-workers [START_REF] Merlo | Cancer as an evolutionary and ecological process[END_REF]. Since we limit ourselves our model to a single tumor site, cell migration by diffusion or by circulating through the blood vessels is not considered in this paper. Nevertheless, we would like to have a model for which the presence (or not) of angiogenesis depends in a simple way on some parameters since angiogenesis is a relevant requirement for an expanding growth of multiple solid tumors [START_REF] Folkman | What is the evidence that tumors are angiogenesisdependent?[END_REF]. The angiogenic switch is a fundamental step which allows a tumor whose size is less than 3 mm to switch from an avascular type (with a slow growth) to a vascular type through new blood vessels under the influence Tumor Angiogenesis Factor (TAF) like Vascular Endothelial Growth Factor (VEGF) produced by tumor cells. The new vessels allow the tumor to get a faster growth and to become more invasive [START_REF] Folkman | Role of angiogenesis in tumor growth and metastasis[END_REF]. The next step is the acquisition of a metastatic phenotype [START_REF] Folkman | Angiogenesis in cancer, vascular, rheumatoid and other disease[END_REF] according to which cancer cells use neo-vessels to migrate at a significant distance from the initial tumor site, most often in other organs, to create new tumor sites [START_REF] Weidner | Tumor angiogenesis and metastasis-correlation in invasive breast carcinoma[END_REF][START_REF] Malanchi | Interactions between cancer stem cells and their niche govern metastatic colonization[END_REF]. Thus by introducing endothelial cells in the model initially proposed by de Pillis and Radunskaya [START_REF] De Pillis | A mathematical tumor model with immune resistance and drug therapy: an optimal control approach[END_REF][START_REF] De Pillis | The dynamics of an optimally controlled tumor model: a case study[END_REF] we constructed a one-site cancer model which is also valid for vascular tumor growth. As performed in these works, we do not consider cellular interactions at the molecular level, such an approach necessarily leading to too complex models for simulating tumor growth at the organ level as we planed to do in further works.

Most of cancer models do not describe the interactions between tumor cells and the organism (the host cells). Among the very few models dealing with host cells [START_REF] Owen | Modelling the macrophage invasion of tumours: effects on growth and composition[END_REF][START_REF] Eisen | Mathematical models in cell biology and cancer chemotherapy[END_REF][START_REF] Dibrov | Mathematical model of cancer chemotherapy. Periodic schedules of phase-specific cytotoxic-agent administration increasing the selectivity of therapy[END_REF][START_REF] Knolle | Cell kinetic modelling and the chemotherapy of cancer[END_REF] is the three-dimensional cancer model proposed by de Pillis and Radunskaya [START_REF] De Pillis | A mathematical tumor model with immune resistance and drug therapy: an optimal control approach[END_REF][START_REF] De Pillis | The dynamics of an optimally controlled tumor model: a case study[END_REF] which retained our attention as explained below. It describes the interactions between tumor, immune and host cells. A chaotic regime was observed in this model [START_REF] Itik | Chaos in three-dimensional cancer model[END_REF] and some bifurcation diagrams were investigated in [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]. Some clinical evidences were well reproduced by this model [START_REF] De Pillis | A mathematical tumor model with immune resistance and drug therapy: an optimal control approach[END_REF][START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF][START_REF] Denis | Chaos Theory and radiotherapy: the tit and the butterfly[END_REF][START_REF] Denis | Chaos theory: a fascinating concept for oncologist[END_REF]. In this latter study, it was shown that this model produces stiff oscillations of the population of tumor cells, corresponding to a fast-growing cancer after a quite long period of latency: it corresponds to tumor dormancy as discussed in [START_REF] De Pillis | A mathematical tumor model with immune resistance and drug therapy: an optimal control approach[END_REF]. We investigated this model using an observability analysis which consists in determining whether a given variable provides all the required information to distinguish states which are different in the original phase space. This is performed by investigating the property of the jacobian matrix of the coordinate transformation between the original phase space and the space reconstructed by using successive time derivatives of a given variable [START_REF] Letellier | Investigating nonlinear dynamics from time series: the influence of symmetries and the choice of observables[END_REF][START_REF] Letellier | Relation between observability and differential embeddings for nonlinear dynamics[END_REF]. Such an observability analysis of this cancer model showed that if one would like to investigate the dynamics of the system "tumor + organism", it would be better to "measure" the number of host cells [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]. It is rather hard to imagine how one could measure the number of host cells since all cells in the body can be considered as host cells but transposed to the clinical point of view, the host dynamics can be evaluated at the body level by the tumor-induced symptoms. For instance, if the population of host cells decreases, the tumor grows in size and symptoms appear. In the case of patients who received treatments for a lung cancer, we thus designed a follow-up based on weekly self-assessed symptoms (lack of apetite, fatigue, pain, cough, breathlessness) and weight to "evaluate" the environment (host cells) of the tumor: the reliability of such a follow-up is equivalent to those of a routine imaging [START_REF] Denis | Detecting lung cancer relapse using self-evaluation forms weekly filled at home: the sentinel follow-up[END_REF][START_REF] Denis | Detection of lung cancer relapse using self-reported symptoms transmitted via an Internet Web-application : pilot study of the sentinel follow-up[END_REF].

Before considering spatio-temporal models for tumor growth, it is relevant to have a model which reproduces the angiogenic switch, an important feature inducing metastasis during tumor growth. It is therefore needed to introduce endothelial cells in the three-dimensional model proposed by de Pillis and Radunskaya. The subsequent part of this paper is organized as follows. Section 2 briefly introduces the three-dimensional model describing the interactions between host, immune and tumor cells. Section 3 briefly describes the dynamics in the "tumor-free-limit", that is, when there is no tumor cell in the site. Section 4 is devoted to the interactions between these first three types of cells with endothelial cells. A dynamical analysis of the resulting four-dimensional model is then performed. Section 5 gives some conclusions.

Three-dimensional cancer model

Mathematical cancer model taking into account normal (non malignant) cells interacting with immune and tumor cells are not numerous. There is one proposed by Owen & Sheratt [START_REF] Owen | Modelling the macrophage invasion of tumours: effects on growth and composition[END_REF] which remains mainly focused on tumormacrophage interactions, the normal cells being only considered for their ability to colonize the site studied. Some others only focused on tumor-host interactions [START_REF] Eisen | Mathematical models in cell biology and cancer chemotherapy[END_REF][START_REF] Dibrov | Mathematical model of cancer chemotherapy. Periodic schedules of phase-specific cytotoxic-agent administration increasing the selectivity of therapy[END_REF][START_REF] Knolle | Cell kinetic modelling and the chemotherapy of cancer[END_REF]]. An interesting model by de Pillis and Radunskaya [START_REF] De Pillis | A mathematical tumor model with immune resistance and drug therapy: an optimal control approach[END_REF][START_REF] De Pillis | The dynamics of an optimally controlled tumor model: a case study[END_REF] incorporates host (normal), immune and tumor cells to reproduce certain qualitative aspects as ocillations in tumor size (Jeff's phenomenon) [START_REF] Tholimson | Measurement and management of carcinoma of the breast[END_REF] or tumor dormancy [START_REF] Farrar | Cancer dormancy vii. A regulatory role for CD8+ T cells and IFN-gamma in establishing and maintaining the tumor-dormant state[END_REF]. This model is rather generic in the sense that it is not specific to a given type of cancer. It is indeed based on quite common interactions between host, immune and tumor cells. The model is

         Ṅ = ρ 2 N(1 -b 2 N) -c 4 T N Ṫ = ρ 1 T (1 -b 1 T ) -c 2 IT -c 3 T N İ = s + ρIT α + T -c 1 IT -d 1 I (1)
where N represents the population of normal (host) cells, T the population of tumor cells and I the population of immune cells. In this model, immune cells can be B-or cytotoxic T-lymphocytes, or even Natural Killer. This is a single tumor-site model where immune cells are considered as arising from outside at a constant rate s as in an immunotherapy [START_REF] Kirschner | Modeling immunotherapy of the tumorimmune interaction[END_REF][START_REF] Kronik | Improving alloreactive CTL immunotherapy for malignant gliomas using a simulation model of their interactive dynamics[END_REF]. Since we are not considering such a case, we set s = 0. System (1) can be thus rewritten in the dimensionless form

       ẋ = ρ 1 x(1 -x) -α 13 xz ẏ = ρ 2 yz 1 + z -α 23 yz -δ 2 y ż = ρ 3 z(1 -z) -α 31 zx -α 32 zy (2)
by using the coordinate transformation (N, T, I) → (x, z, y). Thus, x designates the normalized population of host cells, y being the population of effector immune cells and z the population of tumor cells. Model (2) describes interactions between these three populations as being mainly in competition (Fig. 1). Tumor cells compete for resources with the two other populations of cells. From that point of view, they are the "generalist" competitors while the two others are "specific competitors". System (2) was investigated for the set of parameter values as follows. A previous study investigated in details the dynamics of this three-dimensional system [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]. With the help of bifurcation diagrams, it was shown, for instance, that increasing the growth rate ρ 1 of host cells reduces the probability for observing a large population of tumor cells. Nevertheless, when ρ 1 was large, the maximum of the latter population was nearly 1 and, when it occurred, the tumor growth was very fast, as in these cancers with a very rapid progression that leaves treatment inefficient. Another interesting result was that the tumor cell killing rate α 32 by the effector immune cells had no influence on the dynamics (it only needed to be non-zero), explaining why most of the therapy acting on the immune system are not very efficient [START_REF] Choudhury | Clinical results of vaccine therapy for cancer: learning from history for improving the future[END_REF][START_REF] Chi | Vaccine therapy for metastatic melanoma: systematic review and meta-analysis of clinical trials[END_REF].

With the above parameter values, a chaotic attractor characterized by a smooth unimodal map was observed (Fig. 2). The smooth character of this map is a signature of the period-doubling cascade observed as the route to chaos [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]. A topological analysis lead to the template shown in Fig. 3 where each branch is associated with one of the monotonic branches in the first-return map. The branch without half-turn corresponds to the increasing branch labelled by "0" and the second branch of the template, with a negative half-turn, is associated with the decreasing branch of the first-return map (symbol "1") (see [START_REF] Letellier | Unstable periodic orbits and templates of the Rössler system: toward a systematic topological characterization[END_REF][START_REF] Gilmore | The topology of chaos[END_REF] for a detailed discussion of the topological analysis).

Dynamics in the "tumor-free" limit

The purpose of this section is to check that when there is no tumor cell, the dynamics of the model settles down onto a state which corresponds to a site only made of host cells. When there is no tumor cell, the three-dimensional model reduces to the two-dimensional system

ẋ = ρ 1 x(1 -x) ẏ = -δ 2 y (3)
which has two fixed points

S 0 = 0 0 and S 1 = 1 0 . (4)
The fixed point S 0 corresponds to an empty site (without biological meaning). Its eigenvalues are

Λ 0 = 0.518 -0.5 . (5)
It is therefore a saddle point, that is, an empty site is thus never observed (the corresponding fixed point being unstable). The second fixed point has two negative eigenvalues

Λ 1 = -0.518 -0.5 , (6)
inducing that this is a stable node point. Only initial conditions chosen in the stable manifold of fixed point S 0 induces a trajectory converging to that fixed point. Consequently, almost all initial conditions from the positive quadrant lead to a trajectory ending at this fixed point (Fig. 4). 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9

A four-dimensional model with endothelial cells

In most of cancers, tissue invasion and the occurrence of metastasis is a key step in tumor growth associated with a bad prognostic. Tumor cells become thus more aggressive once tumor hypoxia triggered neo-vessels. Indeed, beyond a given size, numerous tumor cells become in a state of hypoxia and therefore trigger the tumor neo-angiogenesis, that is, induces the production of new (blood and/or lymphatic) vessels issued from the existing ones. Fed with supplementary resources, tumor cells proliferate with an increased growth rate and neo-vascularization opens new routes for tissue invasion and metastasis. This is the so-called angiogenic switch corresponding to the switch from an avascular to a vascular tumor. Such an angiogenic switch cannot exist without the endothelial cells that are responsible for the production of vessels.

Our objective is thus to build a four-dimensional model able to reproduce the angiogenic switch by the means of endothelial cells. Three main interactions have to be considered. The first type of interactions are between tumor cells and endothelial cells. When the tumor diameter is less than 1 or 2 mm, tumor cell needs in oxygen and nutriments are satisfied by blood vessels at a distance not exceeding 50-100 µm, that is, within the oxygen diffusion distance. When the tumor size increases beyond this threshold, tumor needs become too important and can no longer be satisfied by these pre-existing blood vessels. Tumor cells must therefore drain oxygen and nutriments from more distant locations, being beyond the oxygen diffusion limit. Pushed to cellular hypoxia [START_REF] Hillen | Tumour vascularization : sprouting angiogenesis and beyond[END_REF], tumor cells produce Hypoxia Inducible Factor (HIF) proteins which stimulate in turn transcription of Vascular Endothelial Growth Factor genes. The newly produced VEGF protein binds to VEGF receptors on endothelial cells producing vessels. Such a bind induces a signal stimulating a proliferation of endothelial cells. It corresponds to arrow 1 drawn in the flow graph of our four-dimensional model (Fig. 5).

The second type of interactions involving the endothelial cells concerns those with tumor cells. Once their growth is stimulated and that their migration is triggered, endothelial cells structure new blood vessels (in this model we neglect the production of lymphatic vessels). New vessels thus reach the neighborhood of tumor cells and, consequently, increase the proliferation of the latter. In principle, new vasculature would not be selective in the types of cells it provides nourishment and, consequently, it should also benefit to host cells. Nevertheless, only tumor cells have the ability to migrate through neo- vessels and to proliferate more importantly than host cells without saturating the considered tumor site. For this reason, tumor cells are more dependent on endothelial cells than host cells; we therefore neglected interactions between endothelial cells and host cells. The third type of interactions considered in our model is between endothelial cells and effector immune cells. Once cancerogenesis is initiated, effector immune cells become the main actors of the organism reaction against tumor cells. To reach tumor cells, immune cells travel in blood vessels. Thus, the vascular development in the neighborhood of the tumor also increases the ability of immune cells to interact with tumor cells and, consequently, to improve the response of the immune system against the cancer. Endothelial cells secrete chemoattractant molecules attracting immune cells; thus an increase of chemoattractant concentration induced by a growth of endothelial cells induces an additional influx of immune cells to the tumor site [1,[START_REF] Byrne | Using mathematics to study solid tumour growth[END_REF][START_REF] Young | Endothelial cells in the eyes of an immunologist[END_REF][START_REF] Salsman | Crosstalk between medulloblastoma cells and endothelium triggers a strong chemotactic signal recruiting T lymphocytes to the tumor microenvironment[END_REF]. By the same time, endothelial cells may have a negative effect on immune cells [START_REF] Young | Endothelial cells in the eyes of an immunologist[END_REF], a feature that we also took into account in our model by the indirect effect endothelial cells have on immune cells via interactions these two cell populations have with tumor cells. In turn, this decrease of immune cells provokes an increase of the tumor cell population. Tumor can also inhibit immune response [START_REF] Topalian | Safety, activity, and immune correlates of anti-PD-1 antibody in cancer[END_REF][START_REF] Mellman | Cancer immunotherapy comes of age[END_REF]. Although an endogenous immune response to cancer is observed in preclinical models and patients, such response is ineffective because tumors develop multiple resistance mechanisms, including local immune suppression, induction of tolerance and, systemic dysfunction in T-cell signaling. Tumors can exploit different ways to actively evade immune destruction, including endogenous "immune checkpoints" that normally terminate immune response after antigen activation. Blockade of programmed death 1 (PD-1), an inhibitory receptor expressed by T cells, can overcome immune resistance. Thus, tumor-immune cell interactions (described by the term α 23 yz in our model) reinforce the negative impact of tumor cells on immune cells whose population thus decreases.

Since our main objective is to obtain a cancer model reproducing the tumor neo-angiogenesis, we focused our attention on the previous three types of interactions. Other interactions between endothelial cells and host cells or effector immune cells are not relevant in the context of the present model. For instance, we assume that immune cells have no (or a very limited) impact on the development of endothelial cells as done in [START_REF] Brazzoli | A mathematical model of immune competition related to cancer dynamics[END_REF]. More importantly, there is no positive feedback loop of endothelial cells to themselves because their development essentially results from the interactions they have with tumor cells. Consequently, proliferation of endothelial cells is due, in our model, to a coupling between tumor cells and endothelial cells associated with a small natural death rate represented by a negative feedback loop.

The equations and a fixed point stability analysis

According to the flow graph, and based on the three-dimensional model (2), we built a set of four differential equations describing the tumor growth taking into account endothelial cells. The 4D model reads

                 ẋ = ρ 1 x(1 -x) -α 13 xz ẏ = ρ 2 yz 1 + z -α 23 yz -δ 2 y + α 24 yw ż = ρ 3 z(1 -z) -α 31 zx -α 32 zy + α 34 zw 1 + w ẇ = ρ 4 wz 1 + z -δ 4 w (7)
where x represents the population of host cells, y is associated with the population of effector immune cells, z with tumor cells and w corresponds to endothelial cells. Lymphocytes of the immune system circulate continuously between lymph and blood. When an immune response is triggered, effector lymphocytes reach the pathogen site (the tumor) via the blood flow. In other words, effector immune cells use vessels for their migration. Consequently, immune cells are thus promoted by endothelial cells contributing to new vessels through which they can reach the tumor [START_REF] Young | Endothelial cells in the eyes of an immunologist[END_REF]. If immune cells do not require endothelial cells to proliferate, the immune response is obviously stronger when they are present. This positive interaction is represented by the sole term α 24 yw in the second equation of model [START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF]. We assumed that immune cells have no impact on endothelial cells as it was done in [START_REF] Brazzoli | A mathematical model of immune competition related to cancer dynamics[END_REF].

Tumor cells interact with endothelial cells in a reciprocal manner, that is, both types of cells receive benefits from their interactions. When tumor cells are in hypoxia they produce pro-angiogenic factors whose most important is the VEGF. These molecules bind to receptors at the endothelial cell surface, thus leading to proliferation and migration of endothelial cells according to the gradient of pro-angiogenic factor concentration ; endothelial cells are thus attracted by tumor cells. At the cellular level, the presence of tumor cells in the environment triggers the endothelial cell proliferation. At the tissue level, such a feature can be described by a Michaelis-Menten term such as ρ 4 wz 1+z in the fourth equation of model [START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF]; the population of endothelial cells thus saturates at the maximal value ρ 4 . When new vessels are formed, tumor cells are no longer in hypoxia and their proliferation is increased: this is the so-called angiogenic switch leading to the vascular phase. A term reflecting the proliferation of tumor cells due to the presence of endothelial cells is thus added to the logistic term modelling the tumor cell growth. In the third equation of model [START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF] we choose to represent this additional term by α 34 zw 1+w as done for describing the proliferation of endothelial cells.

Both are type ii Holling response functions whose saturation effects are more adequate to model these interactions than a type i Holling response function. To complete this cancer model, the term -δ 4 w in the fourth equation describes the natural death rate of endothelial cells.

Our 4D model ( 7) is investigated with parameter values as α 24 = 0.3 stimulation of effector immune cells by endothelial cells; α 34 = 0.75 tumor cell growth rate due to neo-angiogenesis; ρ 4 = 0.86 endothelial cell growth rate;

δ 4 = 1 11
endothelial cell natural death rate, the other parameters being kept as for the 3D model. These parameter values were chosen because they correspond to a chaotic attractor which is very often observed when parameter values are varied. As for most mathematical models at tissue level, their biological meaning remains uncertain. These non specific parameter values are useful for investigating the qualitative dynamics of tumor growth as performed in [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumor: modeling, applications and biological interpretations[END_REF]. What is relevant is how a change in parameter values affects the dynamics. This will be investigated using some bifurcation diagrams in the subsequent part of this paper.

With these parameter values, the 4D model ( 7) has seven fixed points with positive coordinates, namely

• S 0 = 0 0 0 0 , • S 1 = 1 0 0 0 , • S 2 = 0 0 1 0 , • S 3 = 0 y 3 = 0.394 δ 4 ρ 4 -δ 4 = 0.118 1 α 24 δ 2 - 1 ρ 4 -δ 4 δ 4 (ρ 2 -α 23 ) - ρ 2 δ 2 4 ρ 4 = 0.160 , • S 4 = 0 ρ 3 -ξ α 23 = 0.347 ξ = 0.133 0 , • S 5 = 1 - α 13 ρ 1 ξ = 0.616 ρ 1 ρ 3 (1 -ξ) -α 31 (ρ 1 -α 13 ξ) α 32 ρ 1 = 0.101 ξ = 0.133 0 , • S 6 = 1 - α 13 δ 4 ρ 1 (δ 4 -ρ 4 )
= 0.658

y 6 = 0.131 δ 4 ρ 4 -δ 4 = 0.118 1 α 24 δ 2 - δ 4 [ρ 2 (δ 4 -1)] δ 4 -ρ 4 + α 23 ρ 4 (δ 4 -ρ 4 ) = 0.160 , where ξ = (ρ 2 -α 23 -δ 2 ) -(ρ 2 -α 23 -δ 2 ) 2 -4α 23 δ 2 2α
23 , y 3 and y 6 being too complicated to be explicitely reported. The numerical values here provided correspond to the parameter values previously given.

The first two fixed points correspond to the two fixed points obtained in the tumor-free limit investigated in Section 3. There are also four fixed points which have at least one negative coordinates; they do not contribute to the dynamics since they are not located in the positive domain of the phase space (a population cannot be negative).

• Point S 0 , located at the origin of the phase space, corresponds to an empty site. As for the tumor-free limit, this point must be unstable (since it has not biological meaning, it should not be possible to observe it).

• Point S 1 is associated with a site only inhabited by host cells whose growth is governed by the logistic function ρ 1 x(1x) (also observed in the tumor-free limit). It should be stable, at least for healthy patients.

• Point S 2 corresponds to a site where only tumor cells are observed with the growth according to ρ 3 z(1z). This is thus a pathological state for which tumor cells are at hypoxia: it must be unstable by definition.

• Contrary to this, point S 3 corresponds to a site where tumor cells are at equilibrium with effector immune and endothelial cells: it can be stable, a property which would be of a worse prognostic for the patient than point S 2 since tumor cells can migrate to other sites due to new blood vessels produced by endothelial cells.

• Point S 4 is associated with a site inhabited by immune and tumor cells: a priori, this point corresponds to an avascular tumor which cannot induce metastasis since there is no endothelial cell involved.

The tumor should remain localized and, consequently, could be quite well-treated by a radiotherapy, for instance.

• Point S 5 corresponds to a site where host, immune and tumor cells co-exist: such a state could be associated with a tumor before the angiogenic switch occurs. Host cells are still dominant. The existence of such a point reveals that, in certain cases, rare if this point is unstable, a tumor does not necessarily increase its size up to the angiogenic switch.

• Point S 6 represents a site in which the four populations co-exist, meaning that the angiogenic switch already occurred: this is therefore a vascular tumor. At such a point, the prognostic for the patient would be uncertain since metastasis are very likely expected.

The stability analysis is performed using the jacobian matrix

J =              ρ 1 (1 -2x) -1.5z 0 -1.5x 0 0 ρ 2 z 1+z -0.2z ρ 2 y 1+z -ρ 2 yz (1+z) 2 0.3y -0.5 + 0.3w -0.2y -z -2.5z 1 -2z -x 0.75z 1-w -0.75zw (1+w) 2 -2.5y + 0.75w 1+w 0 0 ρ 4 w 1+z -ρ 4 wz (1+z) 2 ρ 4 z 1+z -1 11             
whose eigenvalues are where the dependence of the eigenvalues on the three parameters ρ 1 , ρ 2 and ρ 4 is explicitely provided when it is simple. λ 1 and λ 2 are solutions to a fourth-degree polynomial.

Λ 0 = -0.5 -0.09 ρ 1 = 0.52 1 , Λ 1 = -ρ 1 = -0.52 -0.5 -0.09 0 , Λ 2 = -1 ρ 1 -1.
The point S 0 is a saddle point (thus unstable). Point S 1 is nearly a stable node: only one eigenvalue is null, inducing a marginal stability. Point S 2 is a saddle and unstable, confirming that a site cannot remain too long without host cells. Fixed point S 3 , corresponding to a vascular tumor, is a saddlefocus SF -with a 1D unstable manifold. Fixed point S 4 associated with an avascular tumor is a saddle-focus SF -with a two-dimensional unstable manifold. Point S 5 , corresponding to a site where tumor cells do not interact with endothelial cells, is a saddle-focus SF + , the two complex conjugated eigenvalues spanning the 3D unstable manifold with the positive real eigenvalue. The last fixed point, S 6 , is also a saddle-focus SF + with a 2D unstable manifold. These four latter fixed points are connected by pairs, namely, S 3 -S 5 and S 4 -S 6 . This could be a signature of some heteroclinic connections, one between S 3 and S 5 , and one between S 4 and S 6 : such connections would be the main "axis" around which the attractor would be structured as shown in Fig. 6 (a feature which remains to further investigate).

We are now investigating how eigenvalues λ 1 and λ 2 , both being real, change their sign under limited parameter value variations. The first eigenvalue

λ 1 = -2δ 4 + ρ 4 2ρ 2 (ρ 2 + δ 2 -α 23 -ξ) (8)
only depends on parameters ρ 2 , ρ 4 , δ 2 , δ 4 and α 23 . The second eigenvalue λ 2 is too complicated to be explicitely provided (its expression exceeded the Maple R capacity). We thus restrict our investigations by varying these five parameter values, one by one, the others being kept as previsously given.

Our results are thus valid for a certain domain of the parameter space which remains to be bounded (postponed for future works). The sign of λ 1 changes when ρ 4 = ρ 4 = 0.777 or ρ 2 = ρ 2 = 4.9536, the other parameter values being kept as previously mentioned; it does not depend on ρ 1 . This bifurcation only affects the dimensionality of the stable (unstable) manifold by ±1, but not the type of the concerned fixed points which remain of the saddle type. When initial conditions are taken in the neighborhood of S 3 for instance, the population of endothelial cells remains very small if not zero for ρ 4 < ρ 4 ; contrary to this, the population quickly increases from any conditions in the neighborhood of S 3 for ρ 4 > ρ 4 . This bifurcation has thus a relevant impact on the dimensionality of the chaotic attractor as discussed in the subsequent part of this paper. The eigenvalues λ 2 change their sign when ρ 4 = 0.9627, that is, for ρ 4 > ρ ∞ = 0.92 for which the trajectory is ejected to infinity as discussed below. This second bifurcation has therefore no impact on tumor growth.

Dynamical analysis

In order to perform the dynamical analysis, we choose to use the four differential embeddings induced by each variable of our four-dimensional model. Let us designate by s the "measured" variable. The corresponding differential embedding is thus spanned by

       X = s Y = ṡ Z = s W = ... s . (9)
The four X-Y plane projections corresponding to the attractors shown in Figs. [START_REF] Boon | Human tumor antigens recognized by T lymphocytes[END_REF] are shown in Figs. [START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF]. As for the 3D model (see a detailed discussion in [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]), the differential embedding induced by variable x is the single one which does not present a domain of the phase space where the trajectories are strongly confined, meaning that host cells should still provide the best observability of the dynamics (Fig. 11a).

In order to browse various patient conditions, we choose to vary two parameter values, namely the growth rate ρ 1 of host cells since we observed that it allows to switch from a common to a dormant cancer [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF], and parameters ρ 2 , ρ 4 , δ 2 , δ 4 and α 23 as previously discussed. Bifurcation diagrams versus one of these parameters are computed as introduced in [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF], that is, using minima and maxima of a given variable to obtain its range of variability. This is easily done in a differentiable embedding by defining the "double" Poincaré section as

P D ≡ (X n , Z n , W n ) ∈ R 3 | Y n = 0, Z n ≷ 0 ( 10
)
where Z n > 0 corresponds to minima and Z n < 0 to maxima of the "measured" variable. The bifurcation diagrams versus the endothelial cell growth rate ρ 4 are shown in Figs. 8a-8c. They are terminated at ρ 4 = ρ ∞ ≈ 0.92, when the trajectory is ejected to infinity. The most remarkable feature is that there is a threshold value ρ 4 = 0.777 under which there is no bifurcation, thus meaning that this growth rate does not affect the dynamics (the tumor growth) when ρ 4 < ρ 4 : there is thus a range of "patient conditions" for which the tumor growth does not depend on the endothelial cell growth rate, and which corresponds to non metastatic patients. Beyond the threshold value ρ 4 , the population of endothelial cells starts to be significantly different from zero and presents chaotic oscillations as the others (Fig. 8c). A patient with an endothelial cell growth rate greater than ρ 4 could present a neo-vascularized tumor as discussed below. To this bifurcation corresponds a threshold value ρ 4 beyond which the angiogenic switch we wanted to describe with our model can be observed. As a consequence, the tumor site can be X=z (tumoral cells) - considered as being avascular when ρ 4 < ρ 4 and vascular otherwise. This means that when the vascularization is large enough (ρ 4 ≥ 0.9), tumor cells saturate the site and start to migrate toward other sites leading to metastasis which are therefore strongly governed by the endothelial cell growth rate. If other parameter values are kept constant, there are thus two different groups of patients, one with a growth rate ρ 4 less than the threshold value and for which patients do not present vascular tumor, and one with a growth rate greater than the threshold value, leading to vascular tumor. It is clear that this threshold value depends on the other parameter values, that is, on other patient conditions. From a dynamical point of view, this means that for ρ 4 < ρ 4 , the chaotic

Endothelial cell growth rate ρ 4 0

Immune cell growth rate ρ 2 0

Endothelial cells w n (c) Extrema of endothelial cells (f) Extrema of endothelial cells attractor can be embedded in the three-dimensional sub-space R 3 (x, y, z) since w = 0. The first-return map to a Poincaré section built from the minima of variable x (Fig. 11b) is a smooth unimodal map exactly as observed in the 3D model for the same parameters value (Fig. 2b). Consequently, this attractor is topologically equivalent to the attractor (Fig. 2a) solution to the 3D model. Similar conclusions can be adressed when the immune cell growth rate ρ 2 is decreased. For values larger than ρ 2 = 4, 9536, the population of endothelial cells remains very small. Below this threshold value, this population of cells start to grow and vascular tumor can be observed. When ρ 2 ≤ 4.21, the population of tumor cells becomes too large to remain in the same tumor site and necessarily starts to migrate toward other sites. When the growth of endothelial cells (immune cells) is too large (small), metastatic cancer can be most likely expected. As show in Figs. 9, the natural death rate δ i affects the dynamics in an opposite way than the corresponding growth rate ρ i (i = (2, 4)). The other parameter values being kept as specified in the beginning of this section, the bifurcation values are δ 2 = δ 2 = 0.4521 and δ 4 = δ 4 = 0.1006, respectively. When parameter values are chosen as in Fig. 7 and α 23 is varied, there is no possibility to obtain a bifurcation in such a way that the population of endothelial cells remains at very small values. For these conditions, only vascular tumors are observed. The population of tumor cells saturates the site for α 23 ≥ 0.38. These five parameters affect the nature of different fixed points. The type of points S 4 and S 5 is changed when ρ 2 , ρ 4 or δ 4 is varied; contrary to this, their eigenvalues and, in particular, λ 1 is nearly unchanged when δ 1 or α 23 is varied. Parameter δ 2 can induce a change in the sign of λ 2 , thus affecting the type of points S 3 and S 6 . No bifurcation was identified in these two latter points when α 23 was varied.

When the angiogenic switch can appear according to patient conditions (for instance, if ρ 2 < ρ 2 , ρ 4 > ρ 4 , δ 2 > δ 2 , or δ 4 < δ 4 , the other parameter values being kept as specified earlier), populations of host and tumor cells increase. As imposed to our model, the population of host cells only slighly increases, remaining below the carrying capacity of the site. Contrary to this, the population of tumor cells can become larger than the one allowed by the carrying capacity of the "isolated" site without supplementary blood vessels (Fig. 8b). This is due to the fact that our model is a single site model. With a site placed in larger environment, supplementary tumor cells would have migrated toward sites through the blood vessels built by endothelial cells. It is also important to note that the growth of tumor cells is not too damaging for host cells, thanks to the additional ressources carried by new blood vessels. All these features correspond to common biological evidences which are qualitatively validating our four-dimensional model. Clinically speaking, let us consider one patient characterized, from our model point of view, by a given set of parameter values. It seems reasonable to assume that a patient in stable conditions would present parameter values constant in time. It is thus possible to imagine various scenarios, only depending on patient conditions, that is, on its parameter values. Let us first consider a patient with an endothelial growth rate ρ 4 = 0.72 < ρ 4 , the other parameter values being those previously reported. Assume that he has a tumor site with very few tumor cells, thus corresponding to initial conditions as x 0 = 1, y 0 = 0.01, z 0 = 0.01 and, w 0 = 0.01. Time series of cell populations are shown in Fig. 12a. First, tumor cells quickly proliferate, soon followed by a rapid growth of the population of effector immune cells, in reaction against the cancer progression. During this first period of time, the population of host cells decreases near zero, and there is a very light increase of endothelial cells. Then, all these populations start to oscillate, excepted the population of endothelial cells which tends to zero. A patient with such a tumor site would never have metastasis as long as none of his parameters change in such a way that the previously discussed bifurcation occurs. In the absence of any therapy, all these parameters should not change significantly over the period of time commonly considered when a tumor growth is clinically observed, that is, over a few years. Let us now take a second case for which the endothelial growth rate ρ 4 = 0.92 > ρ 4 (Fig. 12b). From the same initial conditions as used in the previous case, the time series in the first short period of time (t < 200 arbitrary units of time) presented by populations of host, immune and tumor cells are quite similar to those observed for ρ 4 = 0.72: consequently, clinically speaking, the early tumor progression (t < 200 arbitrary units of time) would not be distinguished from the previous one. In fact, the main departure is in the evolution of endothelial cells which are progressively increasing their population up to be significant (the angiogenic switch already occured) and, then to mainly affect the dynamics of immune and tumor cells. The populations of immune and endothelial cells present nearly synchronous oscillations in response to those of tumor cells. Endothelial cells are now able to build new blood vessels to drain supplementary resources, mostly for tumor cells, in this site. As already mentioned, since our model kept the site isolated from a larger environment, populations of immune and tumor cells become greater than those allowed by the carrying capacity.

In the second case, the patient would have metastasis for t > 380 arbitrary units of time. From the clinical point of view, what would help to predict the occurrence of metastasis would be to investigate the ability of the patient to produce new endothelial cells in response to the demand of tumor cells, a parameter which is not clinically assessed and thus leaving the impression that cancer progression is only governed by stochastic laws.

For large values of ρ 4 , the topology of the attractor (Fig. 13a) is not so different than the one obtained with a 3D model when ρ 1 is increased [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]: a third branch occurs in the first-return map to a Poincaré section (Fig. 13b). The main departure is in the differential structure of the attractor as evidenced by the unusual shape of the map with its flat maximum and its very stiff decreasing branch. Such a map is mainly organized around a period-one orbit (corresponding to tumor oscillation with a moderate amplitude whose neighborhood is the most often visited two or three times before wandering in another neighborhood) and a period-two orbit with one large amplitude tumor oscillation followed by one small amplitude oscillation (Fig. 13c). Nevertheless, we should not forget that the dynamics observed in the isolated site described by our 4D model should be, in principle, affected by interactions with other sites, at least when some of the populations are beyond of the carrying capacity.

As a last investigation, we computed bifurcation diagrams versus the host cell growth rate ρ 1 in two situations: below (ρ 4 = 0.72) and beyond (ρ 4 = 0.92) the bifurcation (Figs. [START_REF] Stoll | A mathematical model of the contribution of endothelial progenitor cells to angiogenesis in tumors: implications for antiangiogenic therapy[END_REF]). As we previously explained, the diagrams we obtained for ρ 4 = 0.72 are similar to those obtained with the 3D model (see [START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]). For ρ 4 = 0.92, the diagrams are less developed than for ρ 4 < ρ 4 (compare the diagrams shown in Figs. 14b to those shown in Figs. 14a) and dormant cancer cannot be observed with such endothelial cell growth rate. The growth rate ρ 1 cannot be increased as done with ρ 4 < ρ 4 , the trajectory being ejected to infinity at ρ 1 ≈ 0.52. Migration of immune and tumor cells would appear for ρ 1 = 0.47, the parameter value at which tumor cells exceed the carrying capacity: such a migration thus explains why dormant cancer is no longer possible, even for largest ρ 1 -value.

Conclusion

The 3D model proposed by de Pillis and Radunskaya is very interesting in the fact that it takes into account interactions of immune and tumor cells (two very often considered populations in cancer models) with host cells (very rarely considered). Its study led us to observe behaviors which are consistent to some clinical observations as, for instance, the fact that the tumor cell killing rate by effector immune cells does not influence the dynamics until it is different from zero and, that some fast growing tumor after dormant cancer can be observed for some parameter values. Nevertheless, this was a single tumor site model which was therefore not able to reproduce relevant phenomena as the occurrence of metastasis. In order to overcome such a limitation, we introduced in this 3D model the population of endothelial cells. Our 4D resulting model reproduces the angiogenic switch, a key phenomenon for tissue invasion and the production of metastasis. This 4D model must be considered as a prerequisite before considering a spatial model for tumor growth.

Contrary to de Pillis and Radunskaya's model which was limited to early stage of tumor growth, our model thus spontaneously reproduces the angiogenic switch which characterizes some developed cancers. In the present one tumor-site version, the population of tumor cells increases beyond the carrying capacity: when it will be connected to other sites, this will be transformed into a tumor-cell migration. Such a migration would thus result from interactions between tumor cells and their microenvironment [START_REF] Liao | Cancer associated fibroblasts promote tumor growth and metastasis by modulating the tumor immune microenvironment in a 4T1 murine breast cancer model[END_REF][START_REF] Raz | An inflammatory vicious cycle: Fibroblasts and immune cell recruitment in cancer[END_REF]. Our model does not take into account the biological complexity of cancer (genomic instability, expression of a given inhibition factor, etc.) but focuses on generic interactions between the different cell populations. It thus allows to reproduce situations observed in vivo or in clinics, as for instance, more or less long latence phases without metastasis as well as strongly invasive tumors providing very quickly metastasis [START_REF] Naumov | Role of angiogenesis in human tumor dormancy: animal models of the angiogenic switch[END_REF][START_REF] Naumov | Tumor-vascular interactions and tumor dormancy[END_REF].

Figure 1 :

 1 Figure 1: Flow graph for the 3D model proposed by de Pillis and Radunskaya describing interactions between host, immune and tumor cells.

  ρ

Figure 2 :

 2 Figure2: Chaotic behavior solution to the 3D model as investigated in[START_REF] Letellier | What can be learned from a chaotic cancer model[END_REF]. Parameter values as in the main text.

Figure 3 :

 3 Figure 3: Template synthesizing the topology of the chaotic attractor (Fig. 2a) solution to the 3D model (2).

Figure 4 :

 4 Figure 4: Phase portrait of the tumor free model (3). Only the host cells survive.

Figure 5 :

 5 Figure 5: Flow graph of our four-dimensional model describing interactions between host, effector immune, tumor and endothelial cells.

6 Figure 6 :

 66 Figure6: Chaotic attractor solution to the four-dimensional model[START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF]. Fixed points structuring the attractor are also shown. Parameter values as in the main text. In the z-w plane projection, the axis joining S 5 to S 6 is shown to structure the folding of the attractor.

Figure 7 :

 7 Figure 7: Differential embeddings induced by each variable of the 4D model. ρ 4 = 0.86 and other parameter values as previously reported.

  Extrema of tumor cells (e) Extrema of tumor cells 0,7

Figure 8 :

 8 Figure 8: Bifurcation diagrams versus the endothelial cell growth rate ρ 4 (ρ 2 = 4.5) and versus the immune cell growth rate ρ 2 (ρ 4 = 0.86). Other parameters values as specified for the other figures.

Figure 9 :

 9 Figure 9: Bifurcation diagrams versus the immune cell death rate δ 2 (δ 4 = 1/11) and the endothelial cell death rate δ 4 (δ 2 = 0.5). Other parameters values as specified for the other figures.

Figure 10 :

 10 Figure 10: Bifurcation diagrams versus the immune cell inhibition rate by tumor cells α 23 . Other parameters values as specified for the other figures.

Figure 11 :

 11 Figure 11: Chaotic regime observed for ρ 4 = 0.72. Other parameter values as in Fig. 7.

Figure 12 :

 12 Figure12: Time series of the four populations of cells produced by model[START_REF] Khar | Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumor regression[END_REF]. Other parameter values as in Fig.7. Initial conditions x 0 = 1, y 0 = 0.01, z 0 = 0.01 and w 0 = 0.01.

  Tumor cells versus timeFigure 13: Three branch chaotic attractor solution to the 4D model. ρ 4 = 0.90 and other parameter values as in previous figures.

Figure 14 :

 14 Figure 14: Bifurcation diagrams versus the host cell growth rate ρ 1 . Other parameters values as specified for the other figures.

Acknowledgements

Louise Viger's Ph.D. thesis is supported by the company Hypérion in collaboration with the companies Takeda and Chugai, and with Roche Group. It is associated with a Cifre agreement.