Benzekry Sébastien
email: benzekry@phare.normalesup.org

Mathematical and numerical analysis of a model for anti-angiogenic therapy in metastatic cancers

Keywords: AMS 2010 subject classification : 35F16, 65M25, 92C50 Anticancer therapy modelling, Angiogenesis, Structured population dynamics, Lagrangian scheme

We introduce and analyze a phenomenological model for anti-angiogenic therapy in the treatment of metastatic cancers. It is a structured transport equation with a nonlocal boundary condition describing the evolution of the density of metastasis that we analyze first at the continuous level. We present the numerical analysis of a lagrangian scheme based on the characteristics whose convergence establishes existence of solutions. Then we prove an error estimate and use the model to perform interesting simulations in view of clinical applications.

Nous introduisons et analysons un modèle phénoménologique pour les thérapies anti-angiogéniques dans le traitement des cancers métastatiques. C'est une équation de transport structurée munie d'une condition aux limites non-locale qui décrit l'évolution de la densité de métastases. Au niveau continu, des estimations a priori prouvent l'unicité. Nous présentons l'analyse numérique d'un schéma lagrangien basé sur les caractéristiques, dont la convergence nous permet d'établir l'existence de solutions. Nous démontrons ensuite une estimation d'erreur et utilisons le modèle pour produire des simulations intéressantes au regard de possibles applications cliniques.

Introduction

During the evolution of a cancer disease, a fundamental step for the tumor consists in provoking proliferation of the surrounding blood vessels and migration toward the tumour. This process, called tumoral neo-angiogenesis establishes a proper vascular network which ensures to the tumour supply of nutrients and allow the tumor to grow further than 2-3 mm diameter. It is also important in the metastatic process by making possible the spread of cancerous cells to the organism which then can develop in secondary tumors (metastases). Thus, an interesting therapeutic strategy first proposed by J. Folkman [START_REF] Folkman | Antiangiogenesis : new concept for therapy of solid tumors[END_REF] in the seventies consists in blocking angiogenesis with the goal to starve the primary tumor by depriving it from nutrient supply. This can be achieved by inhibiting the action of the Vascular Endothelial Growth Factor molecule either with monoclonal antibodies or tyrosine kinase inhibitors. Although the concept of the therapy seems perfectly clear, the practical use of the anti-angiogenic (AA) drugs leaves various open questions regarding to the best temporal administration protocols. Indeed, AA treatments lead to relatively poor efficacy and can even provoke deleterious effects, especially on metastases [START_REF] Paez-Ribes | Antiangiogenic therapy elicits malignant progression of tumors to increased local invasion and distant metastasis[END_REF]. Regarding to these therapeutic failures, it seems that the scheduling of the drug plays a major role. Indeed, as shown in the publication [START_REF] Ebos | Accelerated metastasis after short-term treatment with a potent inhibitor of tumor angiogenesis[END_REF], different schedules for the same drug can lead to completely different results. Moreover, AA drugs are never given in a monotherapy but always combined with cytotoxic agents (also named chemotherapy) which act directly on the cancerous cells. Again, the scheduling of the drugs seems to be highly relevant [START_REF] Riely | Randomized phase II study of pulse erlotinib before or after carboplatin and paclitaxel in current or former smokers with advanced non-small-cell lung cancer[END_REF] and the optimal combination schedule between these two types of drugs is still a clinical open question. Thus, the complex dynamics of tumoral growth and metastatic evolution have to be taken into account in the design of temporal administration protocols for anti-cancerous drugs.

In order to give answers to these questions, various mathematical models are being developed for tumoral growth including the angiogenic process. We can distinguish between two classes of models : mechanistic models (see for instance [START_REF] Billy | A pharmacologically based multiscale mathematical model of angiogenesis and its use in investigating the efficacy of a new cancer treatment strategy[END_REF][START_REF] Lignet | Identifying optimal combinations of anti-angiogenesis drugs and chemotherapies using a theoretical model of vascular tumour growth[END_REF]) try to integrate the whole biology of the processes and comprise a large number of parameters; on the other hand phenomenological models aim to describe the tumoral growth without taking into account all the complexity levels (see [START_REF] Swan | Applications of optimal control theory in biomedicine[END_REF] for a review and [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF][START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF][START_REF] Barbolosi | Optimizing drug regimens in cancer chemotherapy: a simulation study using a PK-PD model[END_REF]). Most of these models deal only with growth of the primary tumor but in 2000, Iwata et al. [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] proposed a simple model for the evolution of the population of metastases, which was then further studied in [START_REF] Barbolosi | Mathematical and numerical analysis for a model of growing metastatic tumors[END_REF][START_REF] Devys | A model describing the growth and the size distribution of multiple metastatic tumors[END_REF]. This model did not include the angiogenic process in the tumoral growth and thus could not integrate a description of the effect of an AA drug. We combined it with the tumoral model introduced by Hahnfeldt et al. [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] which takes into account for angiogenesis. The resulting partial differential equation is part of the so-called structured population dynamics (see [START_REF] Perthame | Transport equations in biology[END_REF] for an introduction to the theory) : it is a transport equation with a nonlocal boundary condition. Its mathematical analysis is not classical because the structuring variable is two-dimensional; as far as we know such models have only been studied in the case where one structuring variable is the age and thus has constant velocity (see [START_REF] Tucker | A nonlinear model of population dynamics containing an arbitrary number of continuous structure variables[END_REF][START_REF] Doumic | Analysis of a population model structured by the cells molecular content[END_REF]). This is not the case in our situation and the theoretical analysis of the model without treatment (autonomous case) was performed in [START_REF] Benzekry | Mathematical analysis of a two-dimensional population model of metastatic growth including angiogenesis[END_REF].

In this paper, we present some mathematical and numerical analysis of the model in the nonautonomous case that is, integrating both cytotoxic and AA treatments and with a general growth field G satisfying the hypothesis that there exists a positive constant δ such that G • ν ≥ δ > 0 where ν is the normal to the boundary. We first simplify the problem by straightening the characteristics of the equation. We perform some theoretical analysis first at the continuous level (uniqueness and a priori estimates) using the theory of renormalized solutions. Then we introduce an approximation scheme which follows the characteristics of the equation (lagrangian scheme). The introduction of such schemes in the area of size-structured population equations can be found in [START_REF] Angulo | Numerical schemes for size-structured population equations[END_REF] for one-dimensional models. Here, we go further in the lagrangian approach by doing the change of variables straightening the characteristics and discretizing the simple resulting equation, in the case of a general class of two-dimensional non-autonomous models. We prove existence of the weak solution to the continuous problem through the convergence of this scheme via discrete a priori L ∞ bounds and establish an error estimate in the case of more regular data.

Finally, we use this scheme to perform various simulations demonstrating the possible utility of the model. First, as a predictive tool for the number of metastases in order to refine the existing classifica-tions of cancers regarding to metastatic aggressiveness. Secondly, the model can be used to test various temporal administration protocols of AA drugs in monotherapy or combined with a cytotoxic agent.

Model

The model is based on the approach of [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF][START_REF] Barbolosi | Mathematical and numerical analysis for a model of growing metastatic tumors[END_REF][START_REF] Devys | A model describing the growth and the size distribution of multiple metastatic tumors[END_REF] to describe the evolution of a population of metastases represented by its density ρ(t, X) with X being the structuring variable, here two-dimensional X = (x, θ) with x the size (=number of cells) and θ the so-called angiogenic capacity. It is a partial differential equation of transport type. The behavior of each individual of the population (metastasis), that is the growth rate G(t, X) of each tumor is taken from [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] and is designed to take into account for the angiogenic process, as well as the effect of both anti-angiogenic (AA) and cytotoxic drugs (CT). Its expression will be established in the following subsection. The model writes

(1.1)    ∂ t ρ(t, X) + div(ρ(t, X)G(t, X)) = 0 ∀(t, X) ∈]0, T [×Ω -G • ν(t, σ)ρ(t, σ) = N (σ) Ω β(X)ρ(t, X)dX + f (t, σ) ∀(t, σ) ∈]0, T [×∂Ω ρ(0, X) = ρ 0 (X) ∀X ∈ Ω.
where Ω, the birth rate β(X), the repartition along the boundary N (σ) and the source term f (t, σ) will be specified in the sequel, T is a positive time and ν is the unit external normal vector to the boundary ∂Ω.

The model of tumoral growth under angiogenic control (Hahnfeldt et al. [16])

Let x(t) denote the size (number of cells) of a given tumor at time t. The growth of the tumor is modeled by a gompertzian growth rate modified by a death term describing the action of a CT. The equation is :

(1.2) dx dt = g 1 (t, x) = ax ln θ x -hγ C (t)H(x -x min ),
where a is a parameter representing the velocity of the growth, θ the carrying capacity of the environment, and the term hγ C (t)H(x-x min ) stands for the effect of a cytotoxic drug, where γ C is the concentration of the CT, x min is a minimal size for the drug to be effective (x min ≥ 1) and the function H is a regularization of the Heaviside function (for example H(t) = 1/2 + 1/2 tanh(t/K), with K being a parameter controlling the slope at 0), in order to avoid regularity issues in the analysis. The idea is now to take θ as a variable of the time, representing the degree of vascularization of the tumor and called "angiogenic capacity". The variation rate for θ derived in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] is :

(1.3) dθ dt = g 2 (t, x, θ) = cx -dθx 2 3 -eγ A (t)H(θ -θ min ),
where the terms cx and -dθx 2/3 represent respectively the endogenous stimulation and inhibition of the vasculature and eγ A (t)H(θ -θ min ) is the effect of an anti-angiogenic drug. The factor 2/3 comes from the analysis of [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] which concluded that the ratio of the stimulation rate over the inhibition one should be homogeneous to the tumoral radius to the square. In the figure 1, we present some numerical simulations of the phase plan of the system (1.2), (1.3). Following [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF], we assume a one compartmental pharmacokinetic for the AA and do the same for the CT (in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] there is no CT). We also assume that the drugs are administered as boli. This gives

γ A (t) = N i=1 D A e -clrA(t-t A i ) H(t -t A i )
where the t A i are the administration times of the AA, D A is the administered dose and clr A the clearance. The expression for the CT is the same, with C instead of A.

Renewal equation for the density of metastasis

We denote X = (x, θ) and

G(t, X) = (g 1 (t, x, θ), g 2 (t, x, θ)). We define b = c d 3 2 and Ω = (1, b) × (1, b)
where b is the maximal reachable size and angiogenic capacity for (x, θ) solving the system (1.2),(1.3) with initial size 1 (see [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by Hahnfeldt et al[END_REF] for a study of this system without the CT term). We consider that each tumor is a particle evolving in Ω with the velocity G. Writing a balance law for the density ρ(t, X) we have

∂ t ρ + div(ρG) = 0, ∀(t, X) ∈]0, T [×Ω that we endow with an initial condition ρ 0 ∈ L ∞ (Ω).
Metastasis do not only grow in size and angiogenic capacity, they are also able to emit new metastasis. We denote by b(σ, x, θ) the birth rate of new metastasis with size and angiogenic capacity σ ∈ ∂Ω by metastasis of size x and angiogenic capacity θ, and by f (t, σ) the term corresponding to metastasis produced by the primary tumor. Expressing the equality between the number of metastasis arriving in Ω per unit time (l.h.s in the following equality) and the total rate of new metastasis created by both the primary tumor and metastasis themselves (r.h.s.), we should have for all t > 0

(1.4) - ∂Ω ρ(t, σ)G(t, σ) • νdσ = ∂Ω Ω b(σ, X)ρ(t, X)dX + f (t, σ)dσ.
We assume that the emission rate of the primary and secondary tumors are equal and thus take f (t, σ) = b(σ, X p (t)) where X p (t) represents the primary tumor and solves the ODE system (1.2)- (1.3). We also assume that the new metastasis created have size x = 1 and that there is no metastasis of maximal size b nor maximal or minimal angiogenic capacity because they should come from metastasis outside of Ω since G points inward all along ∂Ω. An important feature of the model is to assume that the vasculature of the neo-metastasis is independent from the one which emitted it. This means that b(σ, X) = N (σ)β(x, θ) with N (σ) having its support in {σ ∈ ∂Ω; σ = (1, θ), 1 ≤ θ ≤ b} and describing the angiogenic distribution of the metastasis at birth. We assume it to be uniformly centered around a mean value θ 0 , thus we take N (1, θ) = 1 2∆θ 1 θ∈[θ0-∆θ,θ0+∆θ] , with ∆θ a parameter of dispersion of the new metastasis around θ 0 . Following the modeling of [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] for the colonization rate β we take β(x, θ) = mx α , with m the colonization coefficient and α the so-called fractal dimension of blood vessels infiltrating the tumor. The parameter α expresses the geometrical distribution of the vessels in the tumor. For example, if the vasculature is superficial then α is assigned to 2/3 thus making x α proportional to the area of the surface of the tumor (assumed to be spheroidal). Else if the tumor is homogeneously vascularised, then α is supposed to be equal to 1. Assuming the equality of the integrands in (1.4) in order to have the equality of the integrals, we obtain the boundary condition of (1.1).

Analysis at the continuous level

In the autonomous case, that is when G depends only on X and there is no treatment, the analysis of the equation (1.1) has been performed in [START_REF] Benzekry | Mathematical analysis of a two-dimensional population model of metastatic growth including angiogenesis[END_REF]. It was proven the existence, uniqueness, regularity and asymptotic behavior of solutions. We present now some analysis on the equation (1.1) with a more general growth field G than the one defined in the section 1.2.

Let Ω be a bounded domain in R 2 , with ∂Ω being piecewise C 1 except in a finite number of points. Let G : R × Ω → R 2 be a C 1 vector field. We make the following assumption on G :

(2.1) ∃ δ > 0, G • ν(t, σ) ≥ δ > 0 ∀ 0 ≤ t ≤ T, σ ∈ ∂Ω.
We do the following assumptions on the data :

(2.2) ρ 0 ∈ L ∞ (Ω), β ∈ L ∞ (Ω), N ∈ L ∞ (∂Ω), N ≥ 0, ∂Ω N (σ)dσ = 1, f ∈ L ∞ (]0, T [×∂Ω).
Remark 1. In the case of G being the one of the section 1.2 if there is no treatment (that is, if e = h = 0, or t ≤ t 1 ) then we don't have G • ν(t, σ) ≥ m > 0 all along the boundary since G vanishes at the point (b, b). But since the problem was solved in this case (see [START_REF] Benzekry | Mathematical analysis of a two-dimensional population model of metastatic growth including angiogenesis[END_REF]) we consider that the time 0 is the starting time of the treatment and that e or h is positive, which makes the assumption (2.1) true.

Definition 1 (Weak solution). We say that ρ

∈ L ∞ (]0, T [×Ω) is a weak solution of the problem (1.1) if for all test function φ in C 1 ([0, T ] × Ω) with φ(T, •) = 0 T 0 Ω ρ(t, X) [∂ t φ(t, X) + G(t, X) • ∇φ(t, X)] dXdt + Ω ρ 0 (X)φ(0, X)dX + T 0 ∂Ω {N (σ)B(t, ρ) + f (t, σ)} φ(t, σ) = 0 (2.3)
where we denoted B(t, ρ) := Ω β(X)ρ(t, X)dX. Remark 2. By approximating a Lipschitz function by C 1 ones, it is possible to prove that the definition of weak solutions would be equivalent with test functions in W 1,∞ ([0, T ] × Ω) vanishing at time T .

Change of variables

Let X(t; τ, σ) be the solution of the differential equation

(2.4) d dt X = G(t, X) X(τ ; τ, σ) = σ .
For each time t > 0, we define the entrance time τ t (X) and entrance point σ t (X) for a point X ∈ Ω :

τ t (X) := inf{0 ≤ τ ≤ t; X(τ ; t, X) ∈ Ω}, σ t (X) := X(τ t (X); t, X).
We consider the sets

Ω t 1 = {X ∈ Ω; τ t (X) > 0}, Ω t 2 = {X ∈ Ω; τ t (X) = 0} T 1 b Σ 2 ]0, T [×Ω Q 1 Σ 1 Q 1 Q 2 Σ 2 Q 1 T Σ 1 Σ 1 X 1 0 1 τ x x t t t b T ]0, T [×Ω Σ 2 Σ 2 X 2 σ θ θ Figure 2:
The two changes of variables X 1 and X 2 (represented only on the plane θ = 1).

and

Q 1 := {(t, X) ∈ [0, T ] × Ω; X ∈ Ω t 1 }, Q 2 := {(t, X) ∈ [0, T ] × Ω; X ∈ Ω t 2 }.
We also define

Q 1 := {(t, τ, σ); 0 ≤ τ ≤ t ≤ T, σ ∈ ∂Ω} = X -1 (Q 1 ) and notice that Σ 1 := [0, T ]×∂Ω = {(t, X); τ t (X) = 0}, and Σ 2 = {(t, X(t; 0, σ)); 0 ≤ t ≤ T, σ ∈ ∂Ω} = {(t, X); τ t (X) = 0}.
See figure 2 for an illustration. We can now introduce the changes of variables that we will constantly use in the sequel.

Proposition 1 (Change of variables). The maps

X 1 : Q 1 → Q 1 (t, τ, σ) → X(t; τ, σ)
and

X 2 : [0, T ] × Ω → Q 2 (t, Y ) → X(t; 0, Y ) are bilipschitz. The inverse of X 1 is (t, X) → (t, τ t (X), σ t (X)) and the inverse of X 2 is (t, X) → (t, Y (X))
with Y (X) = X(0; t, X). Denoting J 1 (t; τ, σ) = |det(DX 1 )| and J 2 (t; Y ) = |det(DX 2 )|, we have :

(2.5) J 1 (t; τ, σ) = |G(τ, σ) • -→ ν (σ)|e t τ div G(u,X(u;τ,σ))du and J 2 (t; Y ) = e t 0 div G(u,X(u;0,Y ))du
We refer to the appendix for the proof of this result and to the figure 2 for an illustration.

Using these changes of variables we can write for a function f

∈ L 1 (]0, T [×Ω) T 0 Ω f (X)dX = T 0 t 0 ∂Ω f (X 1 (t; τ, σ))J 1 (t; τ, σ)dσdτ + T 0 Ω f (X 2 (t; 0, Y ))J 2 (t; Y )dY.
We want to decompose the equation (1.1) into two subequations : one for the contribution of the boundary term and one for the contribution of the initial condition since they are "independent". Defining (2.6) ρ 1 (t, τ, σ) := ρ(t, X(t; τ, σ))J 1 (t; τ, σ) and ρ 2 (t, y) := ρ(t, X(t; 0, Y ))J 2 (t; Y )

we have, when the solution is regular : ∂ t ρ 1 = (∂ t ρ + div (ρG)) J 1 = 0 and the same for ρ 2 . It is thus natural to introduce the following equations

(2.7) ∂ t ρ 1 (t, τ, σ) = 0 0 < τ ≤ t < T, σ ∈ ∂Ω ρ 1 (τ, τ, σ) = N (σ) B(t, ρ 1 , ρ 2 ) + f (t, σ) 0 < τ < T, σ ∈ ∂Ω
where we denoted

B(t, ρ 1 , ρ 2 ) = t 0 ∂Ω β(X(t; τ, σ)) ρ 1 (t, τ, σ)dσdτ + Ω β(X(t; 0, Y )) ρ 2 (t, Y )dY,

and

(2.8)

∂ t ρ 2 = 0 t > 0, Y ∈ Ω ρ 2 (0; Y ) = ρ 0 (Y ) Y ∈ Ω.
We precise the definition of weak solutions to these equations.

Definition 2. We say that a couple

( ρ 1 , ρ 2 ) ∈ L ∞ ( Q 1 )×L ∞ (]0, T [×Ω) is a weak solution of the equations (2.7)-(2.8) if for all φ 1 ∈ C 1 ( Q 1 ) with φ 1 (T, •) = 0 we have : T 0 t 0 ∂Ω ρ 1 (t, τ, σ)∂ t φ 1 (t, τ, σ)dσdτ dt + T 0 ∂Ω N (σ) B(t, ρ 1 , ρ 2 ) + f (t, σ) φ 1 (t, t, σ) = 0, (2.9)
and for all

φ 2 ∈ C 1 ([0, T ] × Ω) with φ 2 (T, •) = 0 we have (2.10) T 0 Ω ρ 2 (t, Y )∂ t φ 2 (t, Y )dt + Ω ρ 0 (Y ) φ 2 (0, Y )dY = 0.
Remark 3. If ρ 1 is a regular function which solves (2.7), then the weak formulation is satisfied since we have :

T 0 t 0 ∂Ω ρ 1 (t, τ, σ)∂ t φ 1 (t, τ, σ)dσdτ dt = T 0 ∂Ω φ 1 (T, τ, σ) ρ 1 (T, τ, σ)dτ dσ =0 - T 0 t 0 φ 1 (t, τ, σ)∂ t ρ 1 (t, τ, σ)dσdτ dt - T 0 ∂Ω φ 1 (t, t, σ) ρ 1 (t, t, σ)dσdt.
We prove now the following theorem, establishing the equivalence between the problem (1.

(2.11) ρ(t, X) := ρ 1 (t, τ t (X), σ t (X))J -1 1 (t, τ t (X), σ t (X))1 X∈Ω t 1 + ρ 2 (t, Y (X))J -1 2 (t, Y (X))1 X∈Ω t 2
is a weak solution of (1.1).

Proof.

• Direct implication. Let ρ be a weak solution of the equation (1.1). We will prove that ρ 2 defined by (2.6) solves (2.8)

. Let φ 2 ∈ C 1 ([0, T ] × Ω) with φ 2 (T, •) = 0. We define for X ∈ Q 2 φ 2 (t, X) := φ 2 (t, Y (X)) ∈ W 1,∞ (Q 2 )
and we intend to extend it in a Lipschitz function of [0, T ] × Ω so that we can use it as a test function in the weak formulation for ρ (see remark 2). We define, for (t, τ, σ)

∈ Q 1 , φ ε 1 (t, τ, σ) = φ 2 (t, σ)ζ ε (τ ) with ζ ε (τ ) being a truncature function in C 1 ([0, +∞[) such that 0 ≤ ζ ε ≤ 1, ζ ε (0) = 1, ζ ε (τ ) = 0 for τ ≥ ε. Then φ ε 1 ∈ W 1,∞ ( Q 1 ) and we set φ ε 1 (t, X) := φ ε 1 (t, τ t (X), σ t (X)) ∈ W 1,∞ (Q 1 )
since τ t (X) and σ t (X) are Lipschitz from proposition 1. We define then

φ ε := φ ε 1 on Q 1 φ 2 on Q 2 .
The function

φ ε is Lipschitz on Q 1 , Lipschitz on Q 2 and φ ε ∈ C([0, T ]×Ω) since Q 1 ∩Q 2 = {(t, X); τ t (X) = 0}. Thus φ ε ∈ W 1,∞ ([0, T ] × Ω) with φ ε (T, •) = 0. Using φ ε as a test function in (2.3), we have Q1 ρ[∂ t φ ε 1 + G • ∇φ ε 1 ]dXdt + T 0 ∂Ω {N (σ)B(t, ρ) + f (t, σ)} φ ε 1 (t, σ)dtdσ + Q2 ρ[∂ t φ 2 + G • ∇φ 2 ]dXdt + Ω ρ 0 (X)φ 2 (0, X)dX = 0 = I 1 ε + I 2 .
By doing the change of variables X 1 in the term

I ε 1 and noticing that φ ε 1 (t, σ) = φ ε 1 (t, t, σ) = φ 2 (t, σ)ζ ε (t), we obtain I ε 1 = T 0 t 0 ρ 1 (t, τ, σ)∂ t φ 1 (t, σ)ζ ε (τ )dσdτ dt + T 0 ∂Ω B(t, ρ) φ 2 (t, σ)ζ ε (t)dσ ---→ ε→0 0.
Now doing the change of variables X 2 in the second term I 2 and noticing that

∂ t φ 2 (t, Y ) = ∂ t (φ 2 (t, X(t; 0, Y ))) = ∂ t φ 2 (t, X(t; 0, Y )) + G(t, X(t; 0, Y )) • ∇φ 2 (t, X(t; 0, Y ))
gives the result. The equation on ρ 1 is proved in the same way.

• Reverse implication. Let ρ 1 and ρ 2 be solutions of (2.7) and (2.8) respectively. Define ρ(t, X) by (2.11), and consider a test function

φ ∈ C 1 ([0, T ] × Ω) with φ(T, •) = 0. Then φ 1 := φ |Q1 ∈ C 1 (Q 1 ), with φ 1 (T, •) = 0, thus φ 1 (t, τ, σ) := φ 1 (t, X 1 (τ, σ)
) is valid as a test function in the weak formulation of (2.7). In the same way φ 2 (t, Y ) := φ 2 (t, X 2 (Y )) with φ 2 := φ |Q2 is valid as a test function for (2.8). Thus we have

Q1 ρ 1 (t, τ, σ)∂ t φ 1 (t, τ, σ)dσdτ dt + T 0 ∂Ω B(t, ρ 1 , ρ 2 ) φ 1 (t, t, σ)dσdt + T 0 Ω ρ 2 (t, y)∂ t φ 2 (t, y)dtdy + Ω ρ 0 (y) φ 2 (0, y)dy = 0
Doing the changes of variables gives the weak formulation of (1.1).

This theorem simplifies the structure of the problem (1.1). In some sense, it formalizes the method of characteristics in the framework of weak solutions for our problem. The characteristics are straightened (see figure 2) and the directional derivative along the field (t, G) is transformed in only a time derivative. Moreover, integrating the jacobians (which contains the transformation of areas) in the definitions of ρ 1 and ρ 2 , these functions are constant in time. The continuous analysis and discrete approximation of the problem (1.1) is thus simplified.

A priori continuous estimates and uniqueness

In order to obtain a priori properties on the solutions of the equation, we will use the theory of renormalized solutions first initiated by DiPerna-Lions [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] in the case of R n and further developed by Boyer [START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF] in the case of a bounded domain. Let us first recall the result that we will use, which can be found in [START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF]. We need to introduce the following measure on ]0, T [×∂Ω :

dµ G = (G • ν)dtdσ Proposition 2 (Renormalization property). Let ρ ∈ L ∞ (]0, T [×Ω)
be a solution, in the distribution sense, to the equation :

(2.12)

∂ t ρ + div(ρG) = 0. (i) The function ρ lies in C([0, T ]; L p (Ω)), for any 1 ≤ p < ∞. Furthermore, ρ is continuous in time with values in L ∞ (Ω) weak- * . (ii) There exists a function γρ ∈ L ∞ (]0, T [×∂Ω; |dµ G |) such that for any h ∈ C 1 (R), for any test function φ ∈ C 1 ([0, T ] × Ω), and for any [t 0 , t 1 ] ⊂ [0, T ], we have t1 t0 Ω h(ρ)(∂ t φ + div(Gφ))dtdX + Ω h(ρ(t 0 ))φ(t 0 )dX - Ω h(ρ(t 1 ))φ(t 1 )dX - t1 t0 ∂Ω h(γρ)φG • νdtdσ - t1 t0 Ω h ′ (ρ)ρdiv(G)φdtdX = 0 (2.13) Remark 4.
• By approximating the function s → |s| by C 1 functions, it is possible to show that the formula (2.13) stands with h(s) = |s|.

• The second point of the proposition implies in particular that h(ρ) has a trace which is h(γρ).

• In [START_REF] Boyer | Trace theorems and spatial continuity properties for the solutions of the transport equation[END_REF], this proposition is proved in the case of a much less regular field G but with the technical assumption that divG = 0, which is not the case here. Though, the proof can be extended to our case.

Thanks to this result, we can prove the following proposition.

Proposition 3 (Continuous a priori estimates). Let ρ ∈ L ∞ (]0, T [×Ω) be a weak solution of the equation (1.1). The following estimates stand

(2.14) ||ρ(t, •)|| L 1 (Ω) ≤ e t||β|| L ∞ ||ρ 0 || L 1 (Ω) + t 0 e (t-s)||β|| L ∞ ∂Ω |f (s, σ)|dσds and ||ρ|| L ∞ (]0,T [×Ω) ≤ C ∞ (2.15) with C ∞ = (||N || L ∞ ||β|| L ∞ ||ρ|| L ∞ + ||f || L ∞ ) ||G|| L ∞ e T ||div G|| L ∞ + ||ρ 0 || L ∞ e T ||div G|| L ∞ Proof. • Estimate in L 1 .
Let ρ be a weak solution of the equation (1.1). Then in particular it solves (2.12) in the sense of distributions. Thus the proposition 2 applies and gives a trace γρ ∈ L ∞ (]0, T [×∂Ω; |d µG |). Now, by using (2.13) with h(s) = s and the definition of weak solutions to the equation (1.1) we have that for all

φ ∈ C 1 c ([0, T [×Ω) T 0 ∂Ω γρ(t, σ)φ(t, σ)G(t, σ) • νdσdt = T 0 ∂Ω N (σ) Ω β(X)ρ(t, X)dX + f (t, σ) φ(t, σ)dσdt which gives (2.16) -γρ(t, σ)G(t, σ) • ν = N (σ) Ω β(X)ρ(t, X)dX + f (t, σ), a.e.
In view of the remark 4, we know that |ρ| is also a weak solution to the equation (1.1), with initial data |ρ 0 | and boundary data |N (σ)B(t, ρ) + f (t, σ)|. By integrating this equation on Ω and using the divergence formula, we obtain in the distribution sense :

d dt Ω |ρ(t, X)|dX = - ∂Ω G(t, σ) • ν|γρ(t, σ)|dσ = ∂Ω |N (σ)B(t, ρ) + f (t, σ)| dσ and thus d dt Ω |ρ(t, X)|dX ≤ ||β|| ∞ Ω |ρ(t, X)|dX + |f (t, σ)|.
A Gronwall lemma concludes.

• Estimate in L ∞ . Using the proposition 1, we have ρ 1 and ρ 2 solving (2.7) and (2.8). By doing the changes of variables, using the definitions of ρ 1 and ρ 2 and the formulas (2.5), we see that

||ρ(t, •)|| L 1 (Ω) = || ρ 1 (t, •)|| L 1 (]0,t[×∂Ω) + || ρ 2 (t, •)|| L 1 (Ω) , ∀ t > 0 ||ρ|| L ∞ (]0,T [×Ω) ≤ || ρ 1 || L ∞ ( Q1) ||G|| L ∞ (∂Ω) e T ||div G||∞ + || ρ 2 || L ∞ (]0,T [×Ω) e T ||div G||∞
But solving explicitely the equation (2.7), we have

| ρ 1 (t, τ, σ)| = | ρ 1 (τ, τ, σ)| = N (σ) B(t, ρ 1 , ρ 2 ) + f (t, σ) ≤ ||N || ∞ ||β|| ∞ (|| ρ 1 (τ, •)|| L 1 + || ρ 2 (τ, •)|| L 1 ) + ||f || L ∞ ≤ ||N || ∞ ||β|| ∞ ||ρ(τ, •)|| L 1 + ||f || L ∞ On the other hand, for ρ 2 we have || ρ 2 || L ∞ (]0,T [×Ω) = || ρ 2 (0)|| L ∞ (Ω) = ||ρ 0 || L ∞ (Ω) .
Remark 5. The expression (2.16) shows that in the case of a zero boundary data f , the trace γρ has some extra regularity, namely it is C([0, T ]; L 1 (∂Ω)).

Corollary 1 (Uniqueness). If ρ and ρ ′ are two weak solutions of the problem (1.1), then ρ = ρ ′ almost everywhere.

Construction of approximated solutions and application to the existence

In this section, we build a weak solution to the equation (1.1) which, in view of the previous considerations, can be achieved by building a couple (ρ 1 , ρ 2 ) of solutions to the equations (2.7)-(2.8) (recall proposition 1). We will achieve the existence by convergence of an approximation scheme to the problem (2.7)-(2.8) where the difficulty is restricted to the approximation of the boundary condition. Then we establish an error estimate in the case of more regular data. In order to avoid heavy notations, we forget about the tilda when referring to the problem (2.7)-(2.8). We place ourselves in the case where Ω = (1, b) 2 .

Construction of approximated solutions of the problem (2.7)-(2.8)

Let 0 = t 0 < ... < t k < .. < t K+1 = T be a uniform subdivision of [0, T ] with t k+1 -t k = δt. For the equation (2.8), we give ourself uniform subdivisions 1 = x 1 < ... < x l < ... < x L+1 = b and 1 = θ 1 < ... < θ m < ... < θ L+1 = b, with x l+1 -x l = θ m+1 -θ m = δx.
The scheme for the equation (2.8) is then given by :

ρ 0 2 (l, m) = 1 (δx) 2 x l+1 x l θm+1 θm ρ 0 (x, θ)dxdθ 1 ≤ l, m ≤ L ρ k+1 2 (l, m) = ρ k 2 (l, m) 1 ≤ k ≤ K, 1 ≤ l, m ≤ L . (3.1) That is, ρ k 2 (l, m) = ρ 0 2 (l, m) for all k, l, m. For the discretization of the equation (2.7), for each k let 0 = τ 0 < ... < τ i < ... < τ k = t k with τ i+1 -τ i = δt. Let σ : [0, 1] → ∂Ω s → σ(s)
be a parametrization of ∂Ω with |σ ′ (s)| = 1 a.e., so that for g ∈ L 1 (∂Ω) we have ∂Ω g(σ)dσ = 1 0 g(σ(s))ds. Let 0 = s 1 < ... < s j < ... < s M+1 = 1 be an uniform subdivision with s j+1 -s j = δσ. The scheme is given by

   ρ 0 1 (0, j) = N j B 0 ((ρ 0 2 ) l,m ) + f 0 j 1 ≤ j ≤ M ρ k+1 1 (i, j) = ρ k 1 (i, j) 1 ≤ k ≤ K, 0 ≤ i ≤ k, 1 ≤ j ≤ M ρ k+1 1 (k + 1, j) = N j B k+1 (ρ k+1 1 , ρ k+1 2 ) + f k+1 j 1 ≤ j ≤ M (3.2) with B k (ρ k 1 , ρ k 2 ) = k-1 i=1 M j=1 β 1 i,j ρ k 1 (i, j)δtδσ + L l,m=1 β 2 l,m ρ k 2 (l, m) (δx) 2 ≃ t k 0 ∂Ω β(X(t k ; τ, σ))ρ 1 (t k , τ, σ)dτ dσ + Ω β(X(t k ; 0, Y ))ρ 2 (t k , Y )dY and (3.3) β 1 i,j := 1 δtδσ τi+1 τi σj+1 σj β(X(t k ; τ, σ))dσdτ, β 2 l,m := x l+1 x l θm+1 θm β(X(t k ; 0, (x, θ)))dxdθ f k j := 1 δtδσ t k+1 t k σj+1 σj f (t, σ)dσdt, N j := 1 δσ σj+1 σj N (σ)dσ.
Notice that the schemes (3.1) and (3.2) are well-posed since the definition of ρ k+1

1 (k + 1, j) involves values of ρ k+1 1 (i, j) only with 0 ≤ i ≤ k.
We denote by h = δt + δσ + δx and define now the piecewise constant functions ρ 1,h and ρ 2,h on Q 1 and [0, T [×Ω by, for 0

≤ k ≤ K, 1 ≤ i ≤ k, 1 ≤ j ≤ M and 1 ≤ l, m ≤ L (3.4) ρ 1,h (t, τ, σ(s)) = ρ k 1 (i, j) for t ∈ [t k , t k+1 [, τ ∈]τ i-1 , τ i ], s ∈ [s j , s j+1 [ ρ 1,h (t, τ, σ(s)) = 0 for t ∈ [t k , t k+1 [, τ ∈]t k , t], s ∈ [s j , s j+1 [ ρ 2,h (t, x, θ) = ρ k 2 (l, m) for t ∈ [t k , t k+1 [, x ∈ [x l , x l+1 [, θ ∈ [θ m , θ m+1 [.
See the figure 3 for an illustration. Notice that we have (3.5)

t Q 1 0 0 0 0 T k = K + 1 k = 0 k = 1 t k i = 0 i = 1 τ i-1 τ i i = K + 1 τ σ t k+1
||ρ 1,h (t k , •)|| L 1 (]0,t k [×∂Ω) = k i=1 M j=1 ρ k 1 (i, j) δtδσ, ||ρ 2,h (t k , •)|| L 1 (Ω) = M l,m=1 ρ k 2 (l, m) (δx) 2 .
Remark 6.

• We take the same discretization step in x and θ for ρ 2 but it would work the same with two different steps.

• For more regular data, we could take point values instead of (3.3).

• It will be clear from the following that the scheme would converge the same regardless to the value that we give to ρ 0 1 (0, j).

Discrete a priori estimates

We prove the equivalent of the proposition 3 in the discrete case. Notice that there exists a constant C σ such that

M j=1 N j δσ ≤ ∂Ω N (σ)dσ + C σ δσ = 1 + C σ δσ := ||N || h and M j=1 f k+1 j δσ ≤ ||f || L ∞ (]0,T [;L 1 (∂Ω)) + C σ δσ := ||f || h .
Proposition 4 (Discrete a priori estimates). Let ρ k 1 (i, j) k,i,j and ρ k 2 (l, m) k,l,m being given by (3.1) and (3.2) respectively. Then for all k

||ρ 2,h (t k , •)|| L 1 (Ω) = ||ρ 0 || L 1 (Ω) , ||ρ 2,h || L ∞ (]0,T [×Ω) = ||ρ 0 || L ∞ (Ω) (3.6) ||ρ 1,h (t k , •)|| L 1 (]0,t k [×∂Ω) ≤ e t k ||β|| L ∞ ||N || h ||ρ 0 || L 1 (Ω) + ||f || h ||β|| L ∞ ||N || h , (3.7) ||ρ 1,h || L ∞ ( Q1) ≤ ||N || L ∞ ||β|| L ∞ max k ||ρ 1,h (t k , •)|| L 1 + ||ρ 0 || L 1 + ||f || L ∞ .
Moreover, if ρ 0 ≥ 0 then ρ k 1 (i, j), ρ k 2 (l, m) ≥ 0 for all k, i, j, l, m.

Proof. The non-negativity of the scheme is straightforward from the definition. The estimate for ρ 2,h follows directly from the scheme (3.2). For the L 1 estimate on ρ 1,h we compute, using the scheme (3.1)

||ρ 1,h (t k+1 , •)|| L 1 (]0,t k+1 [×∂Ω) = k+1 i=1 M j=1 ρ k+1 (i, j) δtδσ = k i=1 M j=1 ρ k (i, j) δtδσ + B k+1 (ρ k+1 1 , ρ k+1 2 ) δt M j=1 N j δσ + δt M j=1 f k+1 j δσ ≤ ||ρ 1,h (t k )|| L 1 (]0,t k [×∂Ω) + B k+1 (ρ k+1 1 , ρ k+1 2 ) δt||N || h + δt||f || h Now from the expression of B k+1 (ρ k+1 1 , ρ k+1 2 ) B k+1 (ρ k+1 1 , ρ k+1 2 ) ≤ ||β|| L ∞ ||ρ 1,h (t k , •)|| L 1 + ||β|| L ∞ ||ρ 2,h (t k , •)|| L 1 .
Thus we obtain

||ρ 1,h (t k+1 , •)|| L 1 ≤ (1 + ||β|| L ∞ δt||N || h ) ||ρ 1,h (t k , •)|| L 1 + ||β|| L ∞ δt||N || h ||ρ 2,h (t k , •)|| L 1 + δt||f || h
Now using a discrete Gronwall lemma we obtain

||ρ 1,h (t k+1 , •)|| L 1 ≤ e ||β|| L ∞ ||N || h t k ||ρ 1,h (t 0 , •)|| L 1 + ||β|| L ∞ ||N || h ||ρ 2,h (t k , •)|| L 1 + ||f || h ||β|| L ∞ ||N || h Using ||ρ 1,h (t 0 , •)|| L 1 = 0 and ||ρ 2,h (t k , •)|| L 1 (Ω) = ||ρ 0 || L 1 (Ω) ends the proof of the L 1 estimate.
For the L ∞ estimate, we remark that

||ρ 1,h || L ∞ ( Q1) = max k max i,j |ρ k 1 (i, j)| = max k max j B k (ρ k 1 , ρ k 2 )N j + f k j ≤ ||N || L ∞ max k B k (ρ k 1 , ρ k 2 ) + ||f || L ∞ ≤ ||N || L ∞ ||β|| L ∞ max k (||ρ 1,h (t k , •)|| L 1 + ||ρ 2,h (t k , •)|| L 1 ) + ||f || L ∞ .

Application

T 0 Ω ρ 2,hn (t, Y )∂ t φ 2 (t, Y )dY dt = K k=0 L l,m=1 ρ k 2 (l, m) t k+1 t k x l+1 x l θm+1 θm ∂ t φ 2 (t, x, θ)dθdxdt = K k=0 L l,m=1 ρ k 2 (l, m)Φ 2 (t k+1 , l, m)(δx) 2 - K k=0 L l,m=1 ρ k 2 (l, m)Φ 2 (t k , l, m)(δx) 2
where we denoted Φ 2 (t k , l, m) := 1 (δx) 2

x l+1 x l θm+1 θm φ 2 (t k , x, θ)dθdx. Using the scheme (ρ k 2 (l, m) is constant in k) and Φ 2 (t K+1 , l, m) = 0 since t K+1 = T , we obtain

T 0 Ω ρ 2,hn (t, Y )∂ t φ 2 (t, Y )dY dt = L l,m=1 ρ K 2 (l, m)Φ 2 (T, l, m)(δx) 2 - L l,m=1 ρ 0 2 (l, m)Φ 2 (0, l, m)(δx) 2 = - L l,m=1 ρ 0 2 (l, m)Φ 2 (0, l, m)(δx) 2 = - Ω ρ 0 2,hn (Y )φ(0, Y ) ----→ hn→0 - Ω ρ 0 (Y )φ(0, Y )dY since ρ 0 2,hn L 1 ----→ hn→0 ρ 0 .
Observing that the left hand side converges to

T 0 Ω ρ 2 ∂ t φ 2 (t, Y
)dY dt gives the result.

•The function ρ 1 is a weak solution of (2.7). Let φ 1 be a test function for (2.7). Then the same calculation as above shows, with Φ 1 (t k , i, j) := 1 δtδσ τi τi-1 σj+1 σj φ 1 (t k , τ, σ)dσdτ and using that Φ 1 (t K+1 , i, j) = 0 as well as ρ k+1

1 (i, j) = ρ k 1 (i, j) for 1 ≤ i ≤ k and 1 ≤ j ≤ M Q1 ρ 1,hn (t, τ, σ)∂ t φ 1 (t, τ, σ)dσdτ dt = K k=1 k i=1 M j=1 ρ k 1 (i, j)Φ 1 (t k+1 , i, j)δtδσ - K k=1 k i=1 M j=1 ρ k 1 (i, j)Φ 1 (t k , i, j)δtδσ = K i=1 M j=1 ρ K 1 (i, j)Φ 1 (t K+1 , i, j)δtδσ + K-1 k=1 k i=1 M j=1 ρ k 1 (i, j)Φ 1 (t k+1 , i, j)δtδσ - K-1 k=1 k+1 i=1 M j=1 ρ k+1 1 (i, j)Φ 1 (t k+1 , i, j)δtδσ - M j=1 ρ 1 1 (1, j)Φ 1 (t 1 , 1, j)δtδσ = - K-1 k=1 M j=1 ρ k+1 1 (k + 1, j)Φ 1 (t k+1 , k + 1, j)δtδσ - M j=1 ρ 1 1 (1, j)Φ 1 (t 1 , 1, j)δtδσ = - K k=1 M j=1 N j B k (ρ k 1 , ρ k 2 ) + f k j Φ 1 (t k , k, j)δtδσ
Defining the following piecewise constant functions : We need the following lemma in order to conclude.

B h (t, ρ 1,h , ρ 2,h ) = B k (ρ k 1 , ρ k 2 ), N h (σ(s)) = N j , f h (t, σ(s)) = f k j and Φ 1,h (t, σ(s)) = Φ 1 (t k , k, j) on [t k , t k+1 [×[s j ,
Lemma 1. We have

B hn (t, ρ 1,hn , ρ 2,hn ) ⇀ hn→0 B(t, ρ 1 , ρ 2 ) * -L ∞ (]0, T [).
Proof. We define the piecewise constant function β 1 h (τ, σ) as for N h and f h and

β 2 h (X) = β 2 l,m for X ∈ [x l , x l+1 [×[θ m , θ m+1 [. Let t ∈ [t k , t k+1 [, then B h (t, ρ 1,h , ρ 2,h ) = B k (ρ k 1 , ρ k 2 ) = t 0 ∂Ω β 1 h (τ, σ)ρ 1,h (t, τ, σ)dτ dσ- M j=1 β 1 k,j ρ k 1 (k, j)δtδσ+ L l,m=1 β 2 l,m ρ k 2 (l, m)(δx) 2
since we defined ρ h (t, τ, σ) = 0 for τ ∈]t k , t]. Thus, for ψ ∈ L 1 (]0, T [) we have

T 0 B h (t, ρ 1,h , ρ 2,h )ψ(t)dt = T 0 t 0 ∂Ω β 1 h (τ, σ)ρ 1,h (t, τ, σ)ψ(t)dσdτ dt -δt K k=0 M j=1 β 1 k,j ρ k 1 (k, j) t k+1 t k ψ(t)dtδσ + T 0 Ω β 2 h (X)ρ 2,h (t, X)ψ(t)dXdt
and we obtain the result by using ρ 1,hn ⇀ hn→0

ρ 1 * -L ∞ , ρ 2,hn ⇀ hn→0 ρ 2 * -L ∞ , β hn L 1 ----→ hn→0 β, ||β hn || L ∞ ≤ C
and noticing that the second term goes to zero in view of the L ∞ bounds on ρ 1,h (proposition 4) and β.

Using the lemma as well as

N hn , f hn ⇀ hn→0 N, f * -L ∞ , ||N hn || L ∞ ≤ C and Φ 1,hn C([0,T ]×∂Ω) --------→ hn→0 φ(t, t, σ),
the previous calculations give

Q1 ρ 1,hn (t, τ, σ)∂ t φ 1 (t, τ, σ)dσdτ dt ---→ h→0 - T 0 ∂Ω N (σ) B(t, ρ 1 , ρ 2 ) + f (t, σ) φ(t, t, σ)dσdt.
On the other hand the left hand side also goes to Q1 ρ 1 (t, τ, σ)∂ t φ 1 (t, τ, σ)dσdτ dt. This proves that ρ 1 verifies the definition 2 and ends the proof.

Error estimate

We establish now an error estimate for the approximation of the equations (2.7)-(2.8). For this section, we make the following assumptions on the data :

(3.8) ρ 0 ∈ W 1,∞ (Ω), β ∈ W 1,∞ (Ω), N ∈ W 1,∞ (∂Ω), N ≥ 0, ∂Ω N (σ)dσ = 1, f ∈ W 1,∞ (]0, T [×∂Ω)
It can be noticed that in order to perform the weak convergence of the approximated solutions and establish theoretical existence to the continuous problem, we did not need to approximate the characteristics X(t; τ, σ) of the equation. In view of the error estimate though, we need to use another approximation of the data than (3.3). For β(X(t; τ, σ)) we have to introduce an approximation X h (t; τ, σ) of the characteristics given by a numerical integrator of the ODE system (1.2)-(1.3). Then we define (3.9)

β 1 i,j := β(X h (t k ; τ i , σ j )), β 2 l,m := β(X h (t k ; 0, (x l , θ m )) f k j := f (t k , σ j ), N j := N (σ j ).
For g 1 and g 2 being two continuous functions on Q 1 and ]0, T [×Ω respectively, we define

P 1 g 1 (t, τ, σ(s)) = g 1 (t k , τ i , σ j ) for t ∈ [t k , t k+1 [, τ ∈]τ i-1 , τ i ], s ∈ [s j , s j+1 [ P 1 g 1 (t, τ, σ(s)) = 0 for t ∈ [t k , t k+1 [, τ ∈]t k , t], s ∈ [s j , s j+1 [ P 2 g 2 (t, x, θ) = g 2 (t k , x l , θ m ) for t ∈ [t k , t k+1 [, x ∈]x l , x l+1 ], θ ∈ [θ m , θ m+1 [ . Lemma 2 (Projection error). Let (g 1 , g 2 ) ∈ W 1,∞ ( Q 1 ) × W 1,∞ (]0, T [×Ω). Then there exists C P1 and C P2 such that ||g 1 (t k , •) -P 1 g 1 (t k , •)|| L ∞ (]0,t k [) ≤ C P1 h, ||g 2 (t k , •) -P 2 g 2 (t k , •)|| L ∞ (Ω) ≤ C P2 h. (3.10)
We don't give the proof of this lemma since it is classical. We define e 1,h := ρ 1,h -P 1 ρ 1 and e 2,h := ρ 2,h -P 2 ρ 2 the errors of the schemes, with ( ρ 1 , ρ 2 ) solving the problem (2.7)-(2.8). From the equation (2.7) we have

ρ 1 (t k+1 , τ i , σ j ) = ρ 1 (t k , τ i , σ j ), 0 ≤ k ≤ K, 0 ≤ i ≤ k, 1 ≤ j ≤ M ρ 1 (t k+1 , τ k+1 , σ j ) = N (σ j ) B(t k+1 , ρ 1 , ρ 2 ) + f (t k , σ j )
and thus, subtracting this to (3.1) and denoting e k 1 (i, j) = e 1,h (t k , τ i , σ j ) we obtain

(3.11) e k+1 1 (i, j) = e k 1 (i, j), 0 ≤ k ≤ K, 0 ≤ i ≤ k, 1 ≤ j ≤ M e k+1 1 (k + 1, j) = N j E k+1 + r k+1 j with E k+1 = k i=1 M j=1 β 1 i,j e k+1 1 (i, j)δtδσ + L l,m=1 β 2 l,m e k+1 2 (l, m)(δx) 2 r k+1 j = N j B k+1 ( ρ 1 (t k+1 , τ i , σ j )) i,j , ( ρ 2 (t k+1 , x l , θ m )) l,m -B(t k+1 , ρ 1 , ρ 2 ) .
Hence the truncation error of the scheme r k+1 j comes only from the quadrature error coming from the approximation of the integral in B(t k , ρ 1 , ρ 2 ).

Lemma 3 (Truncation error). Assume (3.8) and that

(β•X 1 ) ρ 1 ∈ W 1,∞ ( Q 1 ), (β•X 2 ) ρ 2 ∈ W 1,∞ (]0, T [×Ω).
Then there exists C r such that max k,j

|r k j | ≤ C r h.
Proof. We have

r k j = N j [ k-1 i=1 M j=1 β 1 i,j -β(X 1 (t k ; τ i , σ j )) ρ 1 (t k , τ i , σ j )δtδσ + L l,m=1 β 2 l,m -β(X 2 (t k ; x l , θ m )) ρ 2 (t k , x l , θ m ) (δx) 2 + k-1 i=1 M j=1 β(X 1 (t k ; τ i , σ j )) ρ 1 (t k , τ, σ)δtδσ + L l,m=1 β(X 2 (t k ; x l , θ m )) ρ 2 (t k , x l , θ m ) (δx) 2 - t k-1 0 ∂Ω β(X 1 (t k ; τ, σ)) ρ 1 (t k , τ, σ)dτ dσ - Ω β(X 2 (t k , Y )) ρ 2 (t k , Y )dY - t k t k-1 ∂Ω β(X 1 (t k ; τ, σ)) ρ 1 (t k , τ, σ)dτ dσ].
Thus

r k j ≤||N || L ∞ {||β|| W 1,∞ (||X 1,h -P 1 X 1 || L ∞ ||P 1 ρ 1 || L 1 + ||X 2,h -P 2 X 2 || L ∞ ||P 2 ρ 2 || L 1 ) + k-1 i=1 M j=1 τi τi-1 σj+1 σj |P 1 [(β • X 1 ) ρ 1 ] (t k , τ, σ) -(β • X 1 ) ρ 1 (t k , τ, σ)| dτ dσ + L l,m=1 x l+1 x l xm+1 xm |P 2 [(β • X 2 ) ρ 2 ] (t k , Y ) -(β • X 2 ) ρ 2 (t k , x, θ)| dxdθ + ||(β • X 1 ) ρ 1 || L ∞ h}.
Using the lemma 2 and the L 1 a priori estimate of proposition 3 gives the result.

Remark 7 (Order of the truncation error). In order to have a better order for the truncation error we could use a more sophisticated quadrature method like for instance the trapezoid method on Ω for ρ 2 and on [0, t k-1 [×∂Ω for ρ 1 (completed by a left rectangle method on [t k-1 , t k [×∂Ω). Adapting the previous proof shows that if the numerical integrator used for the characteristics has order larger than 2, then the truncation error would have order 2 (order of the trapezoid method).

Proposition 5 (Error estimate). Assume (3.8) and that ( ρ

1 , ρ 2 ) ∈ W 1,∞ ( Q 1 ) × W 1,∞ (]0, T [×Ω)
||ρ 1,h (t k , •) -ρ 1 (t k , •)|| L 1 (]0,t k [) ≤ C 1 h, ||ρ 2,h (t k , •) -ρ 2 (t k , •)|| L 1 (Ω) ≤ C 2 h (3.12)
Proof. In view of the lemma 2, it is sufficient to prove the proposition with P s ρ s (t k , •) instead of ρ s (t k , •) (with s = 1, 2). For the second estimate, we notice that

||ρ 2,h (t k , •)-P 2 ρ 2 (t k , •)|| L 1 (Ω) = ||e 2,h (t k , •)|| L 1 (Ω) = l,m e k 2 (l, m) (δx) 2 = l,m ρ 0 2 (l, m) -ρ 0 (x l , θ m ) (δx) 2
and the result follows from the definition of ρ 0 2 (l, m). For the first one, we have

||ρ 1,h (t k , •) -P 1 ρ 1 (t k , •)|| L 1 (]0,t k [) = ||e 1,h (t k , •)|| L 1 (]0,t k [) = k i=1 M j=1 e k 1 (i, j) δtδσ.
We can compute, using (3.11)

||e 1,h (t k+1 , •)|| L 1 ≤ k i=1 M j=1 e k+1 1 (i, j) δtδσ + E k+1 δt M j=1 N j δσ + δt M j=1 r k+1 j δσ ≤ ||e 1,h (t k , •)|| L 1 + δt||β|| ∞ ||N || h {||e 1,h (t k , •)|| L 1 + ||e 2,h (t k+1 , •)|| L 1 } + C r hδt ≤ (1 + δt||β|| ∞ ||N || h )||e 1,h (t k , •)|| L 1 + C 2 ||β|| ∞ ||N || h h p δt + C r hδt
and conclude using a discrete Gronwall lemma.

Remark 8 (Order of the error).

• By looking more carefully at the propagation of errors in the proof, we see that if we set ρ 0 2 (l, m) = ρ 0 (l, m) (which is valid under (3.8)), the error on ρ 2 comes only from the projection error.

• If in addition, we follow the remark 7 for the approximation of the data, then the error between ρ 1,h and P 1 ρ 1 would be of order 2 if we had used a trapezoid method for the integral term in B(t k , ρ 1 , ρ 2 ).

Application to the approximation of the problem (1.1)

We explain now how we approximate the solution of (1.1) from the approximation of the solutions of equations (2.7)-(2.8) given by the schemes (3.1), (3.2). We translate formula (2.11) at the discrete level thanks to ρ 1,h , ρ 2,h given by (3.4) and the solutions ρ k 1 (i, j), ρ k 2 (i, j) of the schemes (3.1) and (3.2) :

ρ h (t, X) := ρ 1,h (t, τ t (X), σ t (X))J -1 1,h (t, τ t (X), σ t (X))1 X∈Ω t 1 :=ρ 1,h + ρ 2,h (t, Y (X))J -1 2,h (t, Y (X))1 X∈Ω t 2 :=ρ 2,h . (3.13)
The jacobians of the changes of variables

J 1 (t; τ, σ) = |G(τ, σ) • -→ ν (σ)|e t τ div G(u,X(u;τ,σ))du and J 2 (t; Y ) = e t 0
div G(u,X(u;0,Y ))du are approximated respectively by J 1,h and J 2,h , piecewise constant functions constructed similarly as in (3.4) through J k 1 (i, j) := e T1(k,i,j) and J k 2 (l, m) := e T2(k,l,m) , where T 1 and T 2 are one-dimensional quadrature methods such that T 1 (k, i, j) ≃ t k τi divG(X(s; τ i , σ j ))ds and T 2 (k, l, m) ≃ t k 0 divG(X(s; 0, (x l , θ m )))ds. The errors of these quadrature methods are denoted by r 1 , r 2 and are assumed to be of order α 1 , α 2 :

r 1 := max k,i,j |r 1 (k, i, j)| ≤ C q (δt) α1 , r 2 := max k,l,m |r 2 (k, l, m)| ≤ C q (δt) α2 .
Hence, we have

(3.14) J k 1 (i, j) = J 1 (t k , τ i , σ j )e -r1(k,i,j) , J k 2 (l, m) := e T2(k,l,m) = J 2 (t k , x l , θ m )e -r2(k,l,m) .
We define the following meshes :

V 1 (k, i, j) = {(t, X(t; τ, σ(s))); t ∈ [t k , t k+1 [, τ ∈]τ i-1 , τ i ], s ∈ [s j , s j+1 [} V 2 (k, l, m) = {(t, X(t; 0, (x l , θ m ))); t ∈ [t k , t k+1 [, x ∈ [x l , x l+1 [, θ ∈ [θ m , θ m+1 [} and, for a function g ∈ C([0, T ] × Ω) (3.15) Pg(t, X) = g(t k , X(t k ; τ i , σ j ))1 (t,X)∈V1(k,i,j) + g(t k , X(t k ; 0, (x l , θ m )))1 (t,X)∈V2(k,l,m) .
Remark 9. In the same way as the lemma 2, there exists a constant C P such that for all function

g ∈ W 1,∞ (]0, T [×Ω) ||g -Pg|| L 1 (]0,T [×Ω) ≤ C P h. Theorem 3. Suppose that ρ ∈ W 1,∞ (]0, T [×Ω)
is a regular solution of the equation (1.1) and let ρ h be defined by (3.13). Then there exists a constant C such that

sup t∈[0,T ] ||ρ h (t, •) -ρ(t, •)|| L 1 (Ω) ≤ Ch.
Proof. In view of the remark 9, it is again sufficient to prove the proposition with Pρ instead of ρ.

Let t ∈ [t k , t k+1 [, then ||ρ h (t, •) -Pρ(t, •)|| L 1 (Ω) = ||ρ 1,h (t k , •) -Pρ 1 (t k , •)|| L 1 (Ω t k 1 ) + ||ρ 2,h (t k , •) -Pρ 2 (t k , •)|| L 1 (Ω t k 2 )
with ρ s (t, X) := ρ(t, X)1 X∈Ω t s (s = 1, 2). We do the proof only for ρ 1 since it is similar for ρ 2 . We also don't write the dependency in σ in order to avoid heavy notations. To obtain the complete proof it suffices to add integrals with respect to σ in the following and σ in all the functions. Doing the change of variables X 1 we have, noticing that Pρ

1 (t k , X(t k ; τ )) = P 1 ρ 1 (t k , τ )P 1 J -1 1 (t k , τ ) ||ρ 1,h (t k , •) -Pρ 1 (t k , •)|| L 1 (Ω t k 1 ) = t k 0 ρ 1,h (t k , τ )J -1 1,h (t k , τ ) -P 1 ρ 1 (t k , τ )P 1 J -1 1 (t k , τ ) J 1 (t k , τ )dτ ≤ t k 0 | ρ 1,h (t k , τ )| J -1 1,h (t k , τ )J 1 (t k , τ ) -1 dτ + t k 0 | ρ 1,h (t k , τ ) -P 1 ρ 1 (t k , τ )| dτ + + t k 0 |P 1 ρ 1 (t k , τ )| 1 -P 1 J -1 1 (t k , τ )J 1 (t k , τ ) dτ.
Now we have, using the definition (3.14)

J -1 1,h J 1 -1 = PJ -1 1 e -r1 J 1 -1 ≤ e r1 1 |PJ 1 | |J 1 -PJ 1 | + e -r1 -1 Thus, since 1 J1 L ∞ < ∞ from formula (2.5
) and the fact that G • ν ≥ m > 0, and using |e -r1 -1| ≤ 2r 1 , there exists C J such that

||J -1 1,h (t k , τ )J 1 (t k , τ ) -1|| L ∞ ≤ C J h, and ||1 -P 1 J -1 1 (t k , τ )J 1 (t k , τ )|| L ∞ ≤ C J h.
The last inequality comes from the lemma 2 since J 1 ∈ W 1,1 ((0, +∞) ; L 1 (Γ)) 1,∞ from the formula (2.5).

Using then the continuous and discrete a priori L 1 estimates and the proposition 5 gives the result.

Remark 10. In the case of less regularity on the solution, we still have

ρ h ⇀ h→0 ρ, * -L ∞ (]0, T [×Ω).

Indeed, we write

ρ h = ρ 1,h J -1 1,h + ρ 2,h J -1 2,h = ρ 1,h J -1 1 + ρ 2,h J -1 2 + ρ 1,h (J 1,h -J 1 ) + ρ 2,h (J 2,h -J 2 ). Then we use that for s = 1, 2 J -1 s,h L 1 ---→ h→0 J -1 s as well as J -1 s,h L ∞ ≤ Ce rs with C a constant. Using the theorem 2
for the convergence of ρ 1,h and ρ 2,h gives the result.

Remark 11.

In practical situations we are often only interested in the number of metastases and not in the density ρ itself. Notice that thanks to the formula Ω ρ(t, X)dX = t 0 ∂Ω ρ 1 (t, τ, σ)dσdτ + Ω ρ 0 (X)dX, we don't have to compute the jacobians J 1 , J 2 to get the number of metastases. Yet, we still have to compute the characteristics since they are requested in the computation of the boundary condition (see formula (3.9)).

Numerical simulations 4.1 Simulation technique and parameters

Since our equation is two-dimensional, the computational cost can be relatively high because of the integral term in the boundary, especially for large-time simulations. In order to take into account that the metastases are born with a vasculature very close to a given value θ 0 , we examine replacing the function N (σ) by a dirac measure. In [START_REF] Benzekry | Passing to the limit 2D-1D in a model for metastatic growth[END_REF], we demonstrate that if we take

N (1, θ) = N ε (1, θ) = 1 2ε 1 θ∈[θ0-ε, θ0+ε
] and let ε go to zero the solution of the problem (1.1), in the case of an autonomous velocity field G and initial condition equal to zero, converges to the measure solution of a limit problem consisting in replacing N by a dirac measure in (1, θ 0 ). We use here these considerations to reduce the computational cost and simulate only along the characteristic passing through (1, θ 0 ), that is to say the scheme (3.1) with only one discretization point σ 0 on ∂Ω and N (σ) = 1 σ=σ0 . Moreover, we use a Runge-Kutta method of order 4 for the approximation of the characteristics and a trapezoid method for the approximation of the boundary condition. The values of the parameters for the tumoral growth are taken from [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF], where they were fitted to mice data. Following [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] and [START_REF] Barbolosi | Mathematical and numerical analysis for a model of growing metastatic tumors[END_REF], we take α = 2/3 and fix the value of m arbitrarily. The values of the parameters (without the treatment) are gathered in the table 1. The size (= volume) is expressed in mm 3 though until now it was thought as the number of cells.

Simulations without treatment

A very important issue for clinicians is to determine the number of metastases which are not visible with medical imaging techniques (micro-metastases). Having a model for the density of metastases structured in size allows us to compute the number of visible and non-visible metastases. We took as threshold for a metastasis to be visible a size of 10 8 cells, that is 100 mm 3 by using the conversion 1 mm 3 ≃ 10 6 cells.

In the figure 4, we plotted the result of a simulation showing both the total number of metastases as well as only the visible ones. We observe that at day 20 the model predicts approximately one metastasis though it is not visible. At the end of the simulation, the total number of metastases is much bigger than the number of visible ones. Thus, an interesting application of the model would be to help designing a predictive tool for the total number of metastases present in the organism of the patient. In this perspective, we define a metastatic index as the integral of ρ on Ω :

M I(t) := Ω ρ(t, X)dX.
Of course, this index depends on the values of the parameters, for example on the parameter m, as shown in the table 2. The larger m, the larger the metastatic index. In this table, we remark that at least for M I(1.5) M I(7.5) M I( 15) m = 10 -4 5.80 × 10 -3 6.60 × 10 -2 2.79 × 10 -1 m = 10 -3 5.80 × 10 -2 6.60 × 10 -1 2.81 m = 10 -2 5.80 × 10 -1 6.62 30.1

Table 2: Variation of the number of metastases with respect to m.

times less than 15 days, it seems that the metastatic index is linear in m. Indeed, this can be explained by the fact that at the beginning, most of the metastases come from the primary tumor and not by the metastases themselves (see figure 5.A). This means that the renewal term in the boundary condition of (1.1) could be neglected for small times and that the solution of (1.1) is close to the one of   

∂ t ρ + div(ρG) = 0 -G • νρ(t, σ) = N (σ)β(X p (t)) ρ 0 (X) = 0.
But then, integrating the equation on Ω gives M I(t) = t 0 β(X p (s))ds = m t 0 x p (s) α ds, where X p (s) = (x p (s), θ p (s)) represents the primary tumor and solves the system (1.2)-(1.3) with initial condition (x 0 , θ 0 ). The figure 5.B shows that for larger times metastases emitted by the metastases themselves are more important than the ones emitted by the primary tumor. The metastatic index for large time is then not anymore linear in m (result not shown).

Simulations with treatment

Anti-angiogenic drug alone

We present various simulations of treatments, first involving an anti-angiogenic drug (AA) alone, in order to investigate the difference in effectiveness of various drugs regarding to their pharmacokinetic/pharmacodynamic parameters. The first result shown in figure 6 takes the three drugs which were used in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] where only the effect on tumor growth was investigated, and simulates the effect on the metastases. The three drugs are TNP-470, endostatine and angiostatine and each drug is characterized by two parameters in the model : its efficacy e and its clearance rate clr A . These parameters were retrieved in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] by fitting the model to mice data. The administration protocols are the same for endostatine and angiostatine (20 mg every day) but for TNP-470 the drug is administered with a dose of 30 mg every two days. We observe that TNP-470 seems to have the poorest efficacy due to its large clearance. As noticed in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF], the ratio e/clr A should govern the efficacy of the drug and its value is 0.13 for TNP-470 and 0.39 for both endostatine and angiostatine. The model we developed is now able to simulate efficacy of the drugs on the metastatic evolution (figure 6.C). Interestingly, the drug which seems to be more efficient regarding to the tumor size at the end of the simulation (day 15), namely angiostatine, is not the one which gives the best result on the metastases. Indeed, the lower efficacy of endostatine regarding to ultimate size is due to a relatively high clearance provoking a quite fast rebound of the angiogenic capacity once the treatment stops. But since the tumor size was lower for longer time, the metastatic evolution was better contained. This shows that the model could be a helpful tool for the clinician since the response to a treatment can differ from the primary tumor to metastases, but the clinician has no data about micro-metastases which are not visible with imagery techniques. One of our main postulate in the treatment of cancer is that for a given drug, the effect can vary regarding to the temporal administration protocol of the drug, due to the combination of the pharmacokinetic of the drug and the intrinsic dynamic of tumoral and metastatic growth. To investigate the effect of varying the administration schedule of the drug, we simulated various administration protocols for the same drug (endostatine). The results are presented in figure 7. We gave the same dose and the same number of administrations of the drug but either uniformly distributed during 10 days (endostatine 2), concentrated in 5 days (endostatine 1) or in 2 days and a half (endostatine 3). We observe that the tumor is better stabilized with a uniform administration of the drug (endostatine 2) but the number of metastases is better reduced with the intermediate protocol (endostatine 1). It is interesting to notice that again if we look at the effects at the end of the simulation, the results are different for the tumor size and for the metastases. The two protocols endostatine 1 and endostatine 2 give the same size at the end, but not the same number of metastases. Moreover, the best protocol regarding to minimization of the final number of metastases (endostatine 1) is neither the one which provoked the largest regression of the tumor during the treatment (endostatine 3) nor the one with the most stable tumor dynamic (endostatine 2). In the figure 8, we investigate the influence of the AA dose (parameter D A ) on tumoral, vascular and metastatic evolution. We observe that the model is consistent since it exhibits a monotonous response to variation of the dose.

Combination of anti-angiogenic and cytotoxic drug

An important question in clinical oncology is to determinate how to combine a cytotoxic drug (CT) that kills the proliferative cells and an anti-angiogenic (AA) drug which acts on the angiogenic process, either by blocking the angiogenic factors like VEGF (monoclonal antibodies, e.g. Bevacizumab) or by inhibiting the receptors to this molecule. The AA drugs are classified as part of the cytostatic drugs as they aim to stabilize the disease. For instance, in the treatment of breast cancer, patients which express the receptor HER receive a combination of Docetaxel (CT) and Herceptine (a tyrosine kinase inhibitor, AA). Two questions are still open : which drug should come before the other and then what is the best temporal repartition for each drug? Here, we perform a brief in silico study of the first question. Since we don't have real parameters for the cytotoxic drug we fix arbitrarily the value of each parameter h, clr C and D C to 1, and perform simulations of the model to investigate combination of the CT and the AA. In the figure 9, we present the results of two simulations : one giving the AA before the CT (fig. 9.A) and the other one doing the opposite (fig. 9.B). Although in both cases the effect on the metastases is very good since the growth seems stopped (fig. 9.D), it appears that the qualitative behaviors of the tumoral and metastatic responses are different regarding to the order of administration of the drugs (fig. 9.C and 9.D). According to the model, it would be better to administrate first the CT in order to reduce the tumor burden and then use the AA to stabilize the disease. Indeed, the number of metastasis at the end of the simulation is lower when the CT is applied first than in the opposite case. Of course, this conclusion depends on the tumoral growth and drugs parameters but this simulation shows that the model is able to exhibit different responses regarding to the order of administration between CT and AA drugs.

Conclusion

In this paper, we combined the models of [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] and [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] to obtain a model aiming at describing the effect of anti-angiogenic drugs on the metastatic growth. We established the well-posedness of the model and developed an efficient numerical scheme to perform simulations, which could be adapted to similar models in higher dimensions. The model can now be used in order to rationalize the temporal administration of the anti-angiogenic drugs. To achieve this, we have to implement the various pharmacokinetic models of the different AA drugs and then compare the in silico predictions to real patient data. An important open problem in this direction is the mathematical parameter identifiability of the model, that is to say the inverse problem of uniqueness of the parameters resulting in a given observation. It is also important to develop efficient numerical methods able to achieve the parameter identification from the data. Indeed, identifying the parameters m and α in a given patient could determine the metastatic aggressiveness of its cancer, through the metastatic index. This could lead to interesting clinical applications such as a refinement of the existing classifications like TNM or SBR, which deal only with the visible metastases. As shown in [START_REF] Ebos | Accelerated metastasis after short-term treatment with a potent inhibitor of tumor angiogenesis[END_REF], the metastatic response to AA treatment depends on the time schedule of the drug. The results of the simulations are encouraging in the perspective of using the model as a tool able to test various real temporal administration protocols of the drugs and to perform predictions of the mathematically optimized schedule for a given drug. Moreover, AA are never used in a monotherapy but rather combined with a cytotoxic drug, and determining the best way to combine both drugs is still a clinical open question [START_REF] Riely | Randomized phase II study of pulse erlotinib before or after carboplatin and paclitaxel in current or former smokers with advanced non-small-cell lung cancer[END_REF]. As shown in the figure 9, the model could help in this direction, regarding both to tumor regression and metastatic evolution of the disease. We should also develop further the modeling in order to take into account for the competition effects between CT and AA. Indeed, by reducing the vasculature AA drugs should induce worse supply of both drugs and on the contrary some arguments are expressed in favor of a normalization effect on the tumor vasculature by AA therapy [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combination therapy[END_REF], at least at the beginning of the treatment. These elements should be incorporated to the model via nonlinear terms involving the drugs concentrations in the equations (1.2)-(1.3). The relative simplicity of the model (6 parameters without treatment) is a great advantage in view of concrete applications since we have to be able to fit the model to patients' data in order to retrieve their parameters and then perform predictions about the optimal schedule. A fundamental problem that we have to integrate in our model is the one of toxicities which have to be dynamically controlled to optimize the scheduling of the drug. In the case of CT and on the tumoral growth, a model dealing with hematological toxicities is used to drive phase I clinical trials [START_REF] You | A mechanistic model predicting hematopoiesis and tumor growth to optimize docetaxel + epirubicin (ET) administration in metastatic breast cancer (MBC): Phase I trial[END_REF][START_REF] Barbolosi | Optimizing drug regimens in cancer chemotherapy: a simulation study using a PK-PD model[END_REF]. In

Figure 1 : 3 2

 13 Figure 1: Two phase plans of the system (1.2)-(1.3), for different values of the parameters, without treatment (h = e = 0). In green, the nullclines. In both, b := c d

Figure 3 :

 3 Figure 3: Description of the discretization grid for Q 1 , only in the (τ, t) plane. The arrows indicate the index used in assigning values to ρ 1,h in each mesh (formula (3.4)).

  initial x) θ 0 (initial θ) m α (day -1 ) (day -1 ) (day

Figure 4 :

 4 Figure4: Evolution of the total number of metastases and of the number of visible metastases, that is whose size is bigger than 100 mm 3 (≃ 10 8 cells).

Figure 5 :

 5 Figure 5: Number of metastases emitted by the primary tumor and by the metastases themselves. A. T=50. B. T=100

Figure 6 :

 6 Figure 6: Effect of the three drugs from [16]. The treatment is administered from days 5 to 10. Endostatine (e = 0.66, clr A = 1.7) 20 mg every day, TNP-470 (e = 1.3, clr A = 10.1) 30 mg every two days and Angiostatine (e = 0.15, clr A = 0.38) 20 mg every day. A : Tumor size. B : Angiogenic capacity. C : Number of metastasis.

CFigure 7 :

 7 Figure 7: Three different temporal administration protocols for the same drug (Endostatine). Same dose (20 mg) and number of administrations (6) but more or less concentrated at the beginning of the treatment. Endostatine 1 : each day from day 5 to 10. Endostatine 2 : every two days from day 5 to 15. Endostatine 3 : twice a day from day 5 to 7.5. A : Tumor size. B : Angiogenic capacity. C : Number of metastasis.

Figure 9 :

 9 Figure 9: Combination of an anti-angiogenic drug (AA) : endostatine (e = 0.66, clr A = 1.7 and D A = 20mg) and a cytotoxic one (CT). The parameters for the CT are : h = 1, clr C = 1, and D C = 1. A. AA from day 5 to 10 then CT from day 10 to 15, every day. Tumor growth and angiogenic capacity. B. CT from day 5 to 10 then AA from day 10 to 15, every day. Tumor growth and angiogenic capacity. C : Comparison between both combinations on the tumor growth. D : Comparison between both combinations on the metastatic evolution.

to existence of solutions to the continuous problem

  Theorem 2 (Existence). Under the assumptions (2.2), there existsρ 1 ∈ L ∞ ( Q 1 ) and ρ 2 ∈ L ∞ (]0, T [×Ω) such that ρ 1,h ⇀ 2 for the weak- * topology of L ∞ . Furthermore, (ρ 1 , ρ 2 ) is the unique weak solution of (2.7)-(2.8).Proof. Uniqueness of the solution is straightforward for the problem (2.8) and follows from the L 1 estimate on ρ 1 which can be derived following the proof of the proposition 3. The proof for the existence is rather classical and consists in passing to the limit in discrete weak formulations of (2.7) and(2.8). From the previous proposition, we obtain that the families {ρ 1,h } δt, δσ and {ρ 2,h } δt, δx are bounded in L

	(2.7)-(2.8)

h→0 ρ 1 and ρ 2,h ⇀ h→0 ρ ∞ and thus there exist

ρ 1 ∈ L ∞ ( Q 1 ), ρ 2 ∈ L ∞ (]0, T [×Ω)

and some subsequences ρ 1,hn and ρ 2,hn such that ρ 1,hn ⇀ hn→0 ρ 1 and ρ 2,hn ⇀ hn→0 ρ 2 for the weak- * topology of L ∞ . We have to prove now that (ρ 1 , ρ 2 ) is a weak solution of (2.7)-(2.8). The uniqueness of solutions to the equation implies then by standard argument that the whole sequence converges. It remains to prove that (ρ 1 , ρ 2 ) solves (2.7)-(2.8).

•The function ρ 2 is a weak solution of (2.8). Let φ 2 be a test function for (2.8). We have

Table 1 :

 1 Values of the parameters without treatment.

			-1 vol -2/3 )	(vol)	(vol)	(N b of meta)(day -1 )(vol -α )
	0.192	5.85	8.73 × 10 -3	10 -6	625	10 -3	2/3

our case, we also have to integrate a module to control the toxicity and address the resulting problem of optimization under constraints. Eventually, our model can be used to run in silico tests about the paradigm of metronomic chemotherapies which consists in delivering a cytotoxic drug at low doses and uniformly distributed in the treatment cycle rather than administrating the maximum tolerate dose (MTD) at the beginning of the cycle. Indeed, these metronomic protocols seem to have a dynamical anti-angiogenic effect [START_REF] Hahnfeldt | Minimizing long-term tumor burden : the logic for metronomic chemotherapeutic dosing and its antiangiogenic basis[END_REF][START_REF] Barbolosi | Mathematical modeling of MTD and metronomic temozolomide[END_REF] that can be integrated in the model of [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] for the tumour growth and in our model for the effect on metastases.

A Proof of the proposition 1

The result for the second map is classical. For the first one, we have to deal with irregular points of the boundary ∂Ω. We denote by χ the set of such points and set χ t := {X(t; τ, ξ); ξ ∈ χ, 0 ≤ τ ≤ t}. In order to prove the result, it is sufficient to prove that the map

1 is bilipschitz and that its inverse is X → (τ t (X), σ t (X)). For the first point, since we avoid the irregular points of the boundary by excluding the set χ, we have the C 1 regularity. It remains to prove that X t 1 (τ, σ) is one-to-one and onto, and that its inverse is C 1 .

• The map X t 1 is one-to-one and onto. Let t > 0 and X ∈ Ω t 1 . We have

For the injectivity, we remark that if we have X(t; τ, σ) = X(t; τ ′ , σ ′ ) with for instance τ ′ < τ , then σ = X(τ ; τ ′ , σ ′ ) which is prohibited by the assumption that G • ν(τ, σ) > 0. Thus X t 1 is one-to-one and we have, for (τ, σ)

1 is a diffeomorphism. We will prove the formula (2.5) for J 1 which will conclude the proof by using the local inversion theorem. We have

the derivative in Y of X(t; τ, Y ) viewed as the flow on Ω. We compute

We compute now directly the value of J 1 (t; t, σ). We define

and now notice that we can write

Solving the differential equation between times τ and t and taking the absolute value then gives the formula (2.5).

• Globally,

On the other hand, using the formula (DX t 1 ) -1 = J -1

Acknowledgment

The author would like to express its gratitude to the following people for great support and helpful discussions : D. Barbolosi, A. Benabdallah, F. Hubert and F. Boyer. This work was partially supported by ANR project MEMOREX.