D Barbolosi

A Benabdallah

S Benzekry
email: benzekry@cmi.univ-mrs.fr

J Ciccolini
email: jo.ciccolini@pharmacie.univ-mrs.fr

C Faivre
email: faivre@cmi.univ-mrs.fr

F Hubert
email: fhubert@cmi.univ-mrs.fr

F Verga

B You

A mathematical model for growing metastases on oncologists's service

The dual classification of cancer as localized or metastatic disease is one of the key point in the elaboration of the best therapy for each patient. Nevertheless, many studies reveal that part of these localized diseases are already metastatic. The presence of undetectable or micro-metastases explains the necessity of adjuvant chemotherapies after resection of the primary tumor even for some T1N0M0 cancer. There is probably a continuum between these two stages.

We expose here how a mathematical model of growing metastases could reflect this continuum of the disease and how such a model could help the oncologists in the choice of the treatment. This phenomenological model is based on a structured transport equations with non local boundary condition describing the evolution of the density of metastasis. Thanks to this model, we forge a new numerical index, that we call Metastatic Index, able to reveal either the micro-metastatic state of the patient, or the visible metastatic one. Numerical illustrations show how this new index can be used.

1 A family of model of growing metastases

Tumor growth models with treatment

We assume that the dynamic of the primary tumor is given by some ODE equations. The simplest one is the Gompertz equation

x ′ p (t) = ax p ln θ x p (1)
where a quantifies the growth rate of the tumor and θ denotes the carrying capacity of the media or the maximal size that the tumor can reach. It expresses that the growth rate of the tumor

x ′ p
x p follows some exponential law. The influence of a chemotherapy can then be seen as a shift of this rate:

x ′ p = ax p ln θ x p -γ c (t, x p ) := g(t, x p ) (2)
where γ c reflects the pharmacokinetics and the pharmacodynamics of the chemotherapy treatment. The determination of this function γ c is discussed in [START_REF] Meille | An interface model for dosage adjustment connects hematotoxicity to pharmacokinetics[END_REF][START_REF] Meille | Modélisation du risque des toxicités induites par des substances thérapeutiques[END_REF].

To handle tumoral neo-angiogenesis of the tumor, Folkman & col. propose in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling : a dynamical theory of tumor growth, treatment, response and postvascular dormancy[END_REF] to see θ as a time variable function characteristic of the vascularization of tumor.

x ′ p = ax p ln θ x p -x p γ c (t, x p ) := g 1 (t, x p , θ ), (3)
θ ′ (t) = cx p -dθ x 2 3 p -γ aa (t, θ ) := g 2 (t, x p , θ ) (4)
where the terms cx p and dθ x

2 3
p represent, respectively, the endogenous stimulation and inhibition of the vascularisation and γ aa the action of an anti-angiogenic drug. We refer to these models as

X ′ = G(t, X(t)) := (g 1 (t, x p , θ ), g 2 (t, x p , θ )) (5)

Metastatic growth models

Our first approach consists in deriving a transport equation for the metastatic colony distribution ranged by their size x. We suppose that the primary tumor is generated by a single cell at t = 0 and grows with the rate g(t, x) per unit time where x is the tumor size represented by the number of cells in the tumor. The growing tumor emits metastatic cells with the rate β (x). Each metastatic cell develops into a new tumor, which also grows at rate g(t, x) and emits new nuclei of metastasis just as the primary tumor does. Let v(x,t) represent the colony size distribution with cell number x at time t, actually v(x,t)dx means the number of metastatic tumors whose sizes range from x and x + dx at time t. In the case considered here we suppose that the nuclei of colonization are located far enough from each other so that their ranges do not overlap for a long time period. The dynamics of the colony size distribution is given by the following Mac-Kendrick Von Foerster equation:

   ∂ ∂t v(t, x) + ∂ ∂ x [g(t, x)v(t, x)] = 0, x ∈ [1, b), t ≥ 0 g(1)v(t, 1) = b 1 β (x)v(t, x)dx + β (x p (t)) v(0, x) = 0, (6)
with g and x p defined in [START_REF] Benzekry | Mathematical analysis of a two-dimensional population model of metastatic growth including angiogenesis[END_REF].

We have to deal with a transport equation with a nonlocal boundary condition and an initial data equals to zero which supposes that there is no metastatic tumor at time t = 0. The boundary condition means that the number of metastatic cells newly created per unit time at time t is the total rate of occurences of metastases due to metastatic tumors (corresponding to the integral term) and the primary tumor. The colonization rate β (x) is choosen as β (x) = mx α , where m is the colonization coefficient and α is the fractal dimension of blood vessels infiltrating the tumor. The parameter α expresses how the blood vessels geometrically distribute in or on a tumor. If the vascularity is superficial the fractal dimension α is assigned to be 2/3 because we suppose that the tumor has the shape of a sphere; hence, the surface area is proportional to 2/3. Otherwise, in the case of a vascularization homogeneously distributed in the whole tumor, α is supposed to be equal to 1. This model was introduced by Iwata et al. [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF] for non-treated tumors. An efficient discretization of the problem is proposed in [START_REF] Barbolosi | Mathematical and numerical analysis for a model of growing metastatic tumors[END_REF], the extension of the model to tumor treated by chemotherapy can be found in [START_REF] Verga | Modélisation du risque d'évolution métastatique chez les patients supposés avoir une maladie localisée[END_REF][START_REF] Verga | Modélisation mathématique du processus métastatique[END_REF]. Let us define the metastatic index

MI θ max θ min = θ max θ min v(t, x) dx. (7)
with 1 ≤ θ min ≤ θ max ≤ θ . We recall that θ represents the maximal size of the tumor, its value classicaly taken between 10 11 and 10 12 corresponds to a mass tumor of 100g to 1kg. The choice θ min = 1, θ max < 10 8 leads to the evaluation of micrometases, whereas the choice θ min = 10 8 , θ max = θ corresponds to visible ones. The metastatic index MI is in model ( 6) determined through the four parameters a, θ , m, α. Each of these parameters has to be determined for each patient. The first two parameters can be estimated thanks to two values of the size of the primary tumors. A population study is needed to obtain mean values of parameters of m, α as well as their distribution.

Remark 1 (How to traduce in the model the resection of the primary tumor?). The resection of the primary tumor can be interpreted in our model as sudden diminution of the size of the primary tumor. The function x p becomes close to 0 at T 1 time corresponding to the surgery. Note that Peeters & al have also suggested in [START_REF] Peeters | Decrease in circulating antiangiogenic factors after emoval of primary colorectal carcinoma coincides with increased metabolic acticity of liver metastases[END_REF] that removal of a primary colorectal tumor resulted in an increase in metabolic activity in its liver metastasis (see also [START_REF] Elias | Rationnel de la chirurgie oncologique au sein du traitement multimodal des cancers[END_REF]). They indicate that the primary tumor suppressed angiogenesis in its distant metastasis, and that removal of the primary lesion caused a flare-up in vessel neoformation and, thus, enhanced metabolic activity in its liver metastasis. This stimulation of angiogenic activity could be integrated in our model , by a stiff change in the value of m.

Remark 2 (How to take into account for tumoral angiogenesis?). We assume that the metastatic colony distribution v depends not only on the time t and the size x, but also on the carrying capacity θ . The balance equation then reads

∂ ∂t v(t, X) + div[G(t, X)v(t, X)] = 0 with X = (x, θ ) and G given in (5). ( 8
)
We refer to Benzekry [START_REF] Benzekry | Mathematical and numerical analysis of the anti-angiogenic therapy in metastatic cancers[END_REF][START_REF] Benzekry | Mathematical analysis of a two-dimensional population model of metastatic growth including angiogenesis[END_REF] for more details on the model and its discretization.

2 How can this mathematical model be helpful for oncologist?

Calibration of the parameters of the model

We first exhibit the good behavior of the model [START_REF] Iwata | A dynamical model for the Growth and Size Distribution of Multiple Metastatic Tumors[END_REF]. Koscielny & al made in Institut Gustave Roussy [START_REF] Koscielny | Breast cancer. Relationship between the size of the primary tumour and the probability of metastatic dissemination[END_REF] a retrospective study, over more than twenty years, on 2648 patients with a breast cancer. They reported the percentage of patients (without adjuvant chemotherapy) with at least one metastase (detectable or not) with respect to the size of the initial tumor. Table 1 presents a comparison of these results with percentage from our model. To simulate a panel of patients, we choose the parameters m and α with a normal distribution (the respective parameters of the distribution are µ m = 2.5×10 -8 , σ m = 10 -8 and µ α = 0.48, σ α = 0.3). The similarity between observed and predicted data is emphasized by the value r 2 = 0.98 of correlation coefficient.

Influence of the surgery

In this simulation, the values of the parameters a, θ , α are fixed to a = 4.71 × 10 -4 , θ = 10 12 , α = 0.5. Initialy, the value of the parameter m is 2.5 × 10 -8 . The bigger the primary tumor is, the bigger is m. We assume that the resection is almost radical (10 cells left). To simulate different sizes of tumors, we choose in table 2 m = 2.5 × 10 -7 , m = 2.5 × 10 -6 , m = 2.5 × 10 -7 . Remark that the model is able to quantify the impact of the surgery on the possible recidive of the disease.

Size of the initial tumor MI after 5 years MI after 10 years 10 9 0 1 5 × 10 10 3 13 10 11 5 19

Table 2 Evolution of the metastatic index after surgery

Inter individualization of an adjuvant treatment

Let us now consider a classical treatment of breast cancer [START_REF] Viens | Epirubicin-docetaxel combination in first line chemotherapy for patients with metastatic breast cancer[END_REF]. The standard protocol of 6 cycles of 21 days consists in 75mg of Docetaxel and 100mg of Epirubicin.

The pharmacokinetics and pharmacodynamics of the drugs follow the model developped in [START_REF] Meille | An interface model for dosage adjustment connects hematotoxicity to pharmacokinetics[END_REF][START_REF] Meille | Modélisation du risque des toxicités induites par des substances thérapeutiques[END_REF]. All the patients are assumed to be metastatic at diagnosis with a primary tumor of 1g. Table 3 quantifies the variability of the efficacity of the treatment with respect to m. The model allows to rationalize the number of cycle of chemotherapy to choose for each patient in order to avoid the longer as possible the recidive risk. Table 3 Estimation of the number of visible metastasis MI θ 10 8 with respect to the aggressivity of the tumor. The cases presenting no metastatic risk are in bold.

Combined cytotoxic and antiangiogenic agents

We use the model [START_REF] Meille | Modélisation du risque des toxicités induites par des substances thérapeutiques[END_REF] to exhibit the importance of the scheduling in a combined treatment. We simulate in Figure 1 a 24h injection of 500mg/m2 Etoposide followed or foreruned by a 90 minutes injection of 5mg/kg of Bevacizumab.

Conclusion

The developpment of such metastatic model is fundamental to complete the usual classifications (TNM, SBR, ...) of the disease and to optimize the synergy between the different available treatments (surgery, radiotherapy, drugs...). Besides, the scheduling plays a crucial role in the efficacity of the treatment (see [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] You | Inhibiteurs de tyrosine kinase et chimiothrapie : et si ctait seulement une question de scheduling ?[END_REF]).

Our models are particularly flexible and can describe a large spectrum of therapies, as combined cytotoxic and antiangiogenic agents (with the use of model ( 8)) and then optimize the organization of all treatment for each patient.

We are currently validating in vivo such models on little animal. The choice of phenomenological model is crucial in view of clinical applications. One of the wellknown difficulty of using mathematical approach in clinical application is obviously the identification of the parameters involved in the models. The advantage of these phenomenological models is to involve little parameters that we can expect to estimate properly in a finite time.

Fig. 1

 1 Fig. 1 Influence of the scheduling in a combined therapy. (a) Tumor growth, (b) visible metastasis, (c) carrying capacity, and (d) total number of metastasis

Table 1

 1 Percentage (computed vs observed) of patients having at least one metastase with respect to the size of the primary tumor

	Size of the initial tumor	Computed percentage	Percentage observed in [7]
	1-2.5cm	25.5%	27%
	2.5-3.5cm	44.25%	42%
	3.5-4.5cm	60.5%	56.7%
	4.5-5.5cm	68.6%	66.5%
	5.5-6.5cm	75.5%	72.8%
	6.5-7.5cm	78.25%	83.8%
	7.5-8.5cm	83.25%	81.3%
	>8.5cm	89.25%	92%

Acknowledgements

The authors were partially supported by l'Agence Nationale de la recherche under Grant ANR JC07-07-JCJC-0139-01.