# Spatial Modeling of Tumor Drug Resistance: the case of GIST Liver Metastases

Guillaume Lefebvre, François Cornelis, Patricio Cumsille, Thierry Colin, Clair

Poignard, Olivier Saut

► To cite this version:

Guillaume Lefebvre, François Cornelis, Patricio Cumsille, Thierry Colin, Clair Poignard, et al. Spatial Modeling of Tumor Drug Resistance: the case of GIST Liver Metastases. Mathematical Medicine and Biology, 2014, pp.26. 10.1093/imammb/dqw002. hal-01089452v3

HAL Id: hal-01089452

https://inria.hal.science/hal-01089452v3

Submitted on 26 Jan 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
informatics mathematics
Inria

Spatial Modeling of
Tumor Drug Resistance:
the case of GIST Liver
Metastases

Guillaume Lefebvre, François Cornelis, Patricio Cumsille, Thierry Colin, Clair Poignard, Olivier Saut

RESEARCH
REPORT

Nº 8642

December 2014

Project-Team MC2
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Spatial Modeling of Tumor Drug Resistance: the case of GIST Liver Metastases

Guillaume Lefebvre*†, François Cornelis*‡, Patricio Cumsille§ ¶,
Thierry Colin*†, Clair Poignard†‖, Olivier Saut‖†

Project-Team MC2

Research Report n° 8642 — December 2014 — 26 pages

**Abstract:** This work is devoted to modeling gastrointestinal stromal tumor (GIST) metastases to the liver, their growth and resistance to therapies. More precisely, resistance to two standard treatments based on tyrosine kinase inhibitors (imatinib and sunitinib) is observed clinically. Using observations from medical images, we build a spatial model consisting in a set of nonlinear partial differential equations. After calibration of its parameters with clinical data, this model reproduces qualitatively and quantitatively the spatial tumor evolution of one specific patient. Important features of the growth such as the appearance of spatial heterogeneities and the therapeutical failures may be explained by our model. We then investigate numerically the possibility of optimizing the treatment in terms of progression free survival time and minimum tumor size reachable by varying the dose of the first treatment. We find that according to our model, the progression free survival time reaches a plateau with respect to this dose. We also demonstrate numerically that the spatial structure of the tumor may provide much more insights on the cancer cell activities than the standard RECIST criteria, which only consists in the measurement of the tumor diameter.

**Key-words:** mathematical modeling, numerical simulations, GIST, metastases, resistance and relapse treatments

Contact: guillaume.lefebvre@math.u-bordeaux1.fr

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

<footnote>* Université de Bordeaux, IMB, UMR CNRS 5251, Talence, France</footnote>

<footnote>† Inria Bordeaux-Sud-Ouest, Talence, France</footnote>

<footnote>‡ Service d'imagerie diagnostique et interventionnelle de l'adulte, Hôpital Pellegrin, CHU de Bordeaux, Bordeaux, France</footnote>

<footnote>§ Group of Applied Mathematics (GMA) and Group of Tumor Angiogenesis (GIANT), Basic Sciences Department, Faculty of Sciences, Universidad del Bío-Bío, Campus Fernando May, Chillán, Chile</footnote>

<footnote>¶ Centre for Biotechnology and Bioengineering, University of Chile, Beauchef 850, Santiago, Chile</footnote>

<footnote>‖ CNRS, IMB, UMR 5251, Talence, France</footnote>
# Modélisation spatiale de résistance aux traitements tumoraux :
cas des métastases hepatiques de GIST

**Résumé** : Ce travail est consacré à la modélisation de la croissance et de la résistance aux thérapies de métastases hépatiques en provenance du stroma gastrointestinal (GIST). Plus précisément, la résistance à deux traitements standard à base d'inhibiteur de tyrosine kinase (imatinib et sunitinib) est observée cliniquement. En utilisant les observations provenant des images médiacles, nous construisons un modèle spatial qui est un système non linéaire d'équation aux dérivées partielles. Après calibration des paramètres avec les données cliniques, le modèle reproduit qualitativement et quantitativement l'évolution spatiale de la tumeur pour un patient spécifique. D'importantes caractéristiques de la croissance, comme l'apparition d'hétérogénéités spatiales ou les échecs thérapeutiques, peuvent être expliquées par notre modèle. Nous avons ensuite étudié numériquement la possibilité d'optimiser le traitement en termes de survie en progression libre et de taille minimale atteinte, en faisant varier la dose du premier traitement. Selon notre modèle, le temps de survie en progression libre atteint un plateau en fonction de la dose. Nous démontrons également numériquement que la structure spatiale de la tumeur pourrait fournir plus d'informations sur l'activité des cellules cancéreuse que le critère RECIST standard qui consiste seulement en la mesure du diamètre de la tumeur.

**Mots-clés** : modélisation mathématique, simulations numériques, GIST, métastases, résistance et rechute aux traitements
# 1 Introduction

Gastrointestinal stromal tumors (GIST) are the most common mesenchymal tumors of the gastrointestinal tract, with an incidence of 9-14 cases per million people per year (see [16]). In 25% of cases (see [11]), this type of cancer spreads to the liver. Even though GISTs resist to most of conventional cancer chemotherapies, the discovery of activating mutations of the KIT as well as the role of PDGFR and the new subsequent therapeutic development have revolutionized GIST treatments. Thanks to the availability of these highly active targeted therapeutic agents, GISTs have become typical models of personalized treatment of cancer [6]. In particular, the survival of patients with GIST has been improved with the use of tyrosine kinase inhibitors such as imatinib in first-line setting and a multi-targeted receptor tyrosine kinase inhibitor, such as sunitinib or sorafenib, that inhibits PDGFRs, VEGFRs and KIT, as second-line treatment. However several limitations in terms of diagnosis and outcomes still remain.

First, an important variability exists in the molecular and genetic characteristics that drive the pathogenesis of these tumors. Hirota et al. have proved that molecular alterations in KIT gene with gain-of-function mutations occurred in these tumors (see [12]). Moreover, in addition to the primary mutation, secondary mutations have been identified in patients with advanced GIST pretreated with tyrosine kinase inhibitor. At this time, 10 different molecular subsets of GIST with different molecular alterations have been reported. For patients with KIT mutations, an imatinib resistance is frequently observed, as reported in [5]. For other patients, imatinib controls the metastatic disease during a more or less long period, around 20-24 months in 85% of cases. Physicians have then to switch to another molecule or use an alternative therapy. Since the prognosis and the sensitivity to the targeted treatments have been reported to be patient-dependent, we aim at developing a patient-dependent mathematical model based on medical images of liver metastases. We focus on locally advanced GIST in order to determine, for each patient, the time of emergence of mutations in cancer cells, the relapse times after the first-line and the second-line treatments, as well as the geometric features of tumor growth.

Second, the new anticancer agents with targeted mechanisms of action as used for the treatment of GIST have demonstrated the inherent limitation and unsuitability of the usual anatomic tumor evaluation, that only considers the largest diameter of the lesion (i.e. the RECIST criteria, see [20]). For clinicians, the challenge consists in optimizing these cancer treatments and in particular to determine the more adequate time to switch from the first-line to the second-line treatment, in order to increase the overall survival. The estimate of the relapse time is therefore crucial.

Clinical follow-up to monitor the disease evolution is mainly performed with CT-scans. We emphasize that the effect of these new drugs changes the paradigm according to which the tumor sensitivity to the treatment is measured (see [19]), since CT-scans have reported other information such as tumor heterogeneity: the RECIST criteria seems no more sufficient.

The aim of this work is to provide a spatial model of standard treated GISTs in order to compare the model with the images, and possibly to highlight the peculiarities of the tumor growth or regrowth. It is worth noting that this paper is a first step in the modeling of tumor drug resistance based on clinical images.

We provide a model, which consists in a non-linear system of partial differential equations (PDEs), in order to account for the spatial aspect of the tumor growth. Actually, the models based on ordinary differential equations (ODEs) as the models of Mendelsohn, Gompertz or Bertalanffy make it possible to track the tumor area growth but they do not consider the spatial aspects of the growth. We refer to the review by Benzekry et al. for more details on such 1D models [3] Our model is derived in the same vein as Ribba et al. [17] – we also refer to [9,10] – in order to describe the evolution of the disease. The main novelty of the model lies in the description of the treatments. Two treatments are considered: the first treatment consists in a cytotoxic effect while the second-line treatment has both cytotoxic and anti-angiogenic effects. 3 different proliferative cells are used to describe the resistance to the treatment: one cell population is sensitive to the two treatments, another one is only sensitive to the second treatment while the third cell population is resistant to both treatments. We also provide a simple model of angiogenesis, which is
crucial in the metastasis growth.

Once the model is built, we provide an appropriate numerical scheme that makes it possible to solve the PDEs. In particular, we present in Section 3 a new WENO5 type scheme that stabilizes the computation by using a combination of the classical WENO5 stencil and a twisted stencil. We then extensively compare our model to the clinical data of one specific patient, for whom we have the whole clinical protocol, in Section 4. Let us mention that it has been necessary to introduce a numerical reconstitution of the CT-scan from the numerical results, in order to compare the grey-levels of CT-scans with our simulations. Once the fitting has been obtained, we then investigate numerically the effect of the treatment dose on the tumor growth progression. The counter-intuitive result lies in the fact that according to the parameters we use for the fitting with the data, the increase of the dose of the first treatment does not improve the progression free survival rate. This result is explained in subsection 4.2.2. We then conclude on the consistency of our model by presenting the different behaviors of the tumor evolution that can be obtained. We also fit the numerical tumor area with another patient, whose tumor is close to the liver boundary. It is worth noting that for such patients, the shape of the tumor cannot be recovered. Our model does not take the corresponding mechanical constraints into account. However, the tumor area seems to be well reproduced by our simulations.

## Main Insights

On the left part of Figure 1 we provide the sequence of the CT-scans, and on the right part we give the tumor area evolution of the GIST liver metastasis of one specific patient called patient A. On each CT-scan of patient A, we have depicted the darkest region, which corresponds mainly to the necrotic cells. During the tumor evolution, one can see that the heterogeneity of the tumor changes (for instance in Figure 1b the tumor is homogeneous while in Figure 1f is heterogeneous). We distinguish in the tumor area evolution the points corresponding to a homogeneous tumor (in fulfilled circle) from the point corresponding to heterogeneous tumor.

Interestingly, just before the tumor regrowth at Day 776 (Figure 1c) and Day 1116 (Figure 1f), one can see on the tumor a rim clearer than the dark core, while the response to the treatment is followed by a darkening of the tumor (see Figure 1b and 1e). Such successive stages of tumor homogeneities/heterogeneities, are particularly pronounced in GIST liver metastases, and our goal is to provide an explanation of such behavior.

The first results of this paper lies in the fact that it makes it possible to describe the tumor evolution in terms of tumor area compared with the CT-scans measurements, as presented by Figure 1g: the continuous line corresponds to the numerical results. It is worth noting that we do not provide a 1D-model that describes the tumor volume. We deal with a complex non-linear PDE model, which is phenomenological and which describes the behaviour of cancer cells with respect to the space and the time variables. Therefore this first fit with the tumor area is a non trivial insight.

The second main insight is that our model brings new information on the tumor structure, that seems corroborated with the CT-scans. Actually, as it will be presented in the following, we link the tumor heterogeneity to an increase in the cellular activity, meaning that a resistant phenotype is emerging in the clearer region. According to our modeling, such a behavior can be seen before the treatment failure, while RECIST criteria has not changed. For instance in Figure 2, we compare the structure of the tumor at two different days: the heterogeneity of the tumors seems to be well captured by the numerical simulations. Therefore our paper can be seen as a first step in developing new tools to evaluate the tumor response to treatment based on tyrosine kinase inhibitors.
<figure><img src="image_4.png" /></figure>

(a) Sept 16, 2008 –
Day 119

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

(b) June 30, 2009 –
Day 406

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

(c) July 5, 2010 –
Day 776

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

(d) Oct 25, 2010 –
Day 888

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

(e) Jan 7, 2011 –
Day 962

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

(f) June 10, 2011 –
Day 1116

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

(g) Tumor area. Each point represents the tumor area measured on CT-scans and the line stands for our numerical results. The letter refers to the CT-scans shown on the left. The symbol ⊕ stands for heterogeneous tumor, ● stands for rather homogeneous tumor and ○ stands for CT-scan on which it is difficult to detect heterogeneous aspect.

The value of parameters used in the numerical simulation are given in Table 2.

<figure><img src="image_11.png" /><figcaption>**Figure 2.** Comparison between CT-scans of patient A and the numerical simulations.</figcaption></figure>
## 2 Description of the model

Our model is a system of PDEs similar to the model of Bresch et al. [8]. Roughly speaking, the tumor is described thanks to densities of proliferative and necrotic cancer cells denoted by $P$ and $N$ respectively. Cell proliferation leads to an area increase in the tissue, which creates a pressure whose gradient transports the surrounding healthy cells $S$ away from the core of the tumor. The cells located at the center of the tumor turn into a necrotic phase, because of a lack of oxygen for instance, except if angiogenesis has occured to provide them nutrient supply. Angiogenesis and nutrient supply are taken into account with a simplistic description similar to [4, 18]: since the vascularization drives the nutrient concentration towards the tumor, we introduce a variable $M$ that describes both vascularization, neovascularization and nutrients brought to the tumor thanks to an advection-diffusion equation. We also introduce growth factor effects through a variable $\xi$ that modulates $M$. The set of quantities used in our model is resumed in Table 1.

The main novelty lies in the modeling of the treatments. Two treatments are considered (one can see in [15], the recent works by Lorz et al. for more elaborated models on drug resistance). The first one is a specific tyrosine kinase inhibitor, such as imatinib, that has a cytotoxic effect on proliferative cells. The second treatment is a multitargeted kinase inhibitor, such as sunitinib or sorafenib, that has both cytotoxic and anti-angiogenic effects, meaning that in addition to the cytotoxic effect, it blocks the production of growth factors such as vascular endothelial growth factors (VEGF) and thus decreases the nutrient supply brought to the tumor. It is well-known that the cytotoxic drugs do not impact similarly all the metastatic cancer cells since resistant phenotype can appear in the proliferative cell population. Moreover it is well-known that cancer cells can resist differently to hypoxia. Therefore we split, as in [7], the density $P$ of proliferative cells ($P$-cells) into 3 subpopulations $P_1$, $P_2$ and $P_3$, such that $P = P_1 + P_2 + P_3$, where

* $P_1$ denotes the fraction of proliferative cells that are sensitive to the first-line treatment $\mathcal{T}_1$, based on imatinib molecule and also to the second-line treatment $\mathcal{T}_2$, based on sunitinib or sorafenib, that has both cytotoxic and antiangiogenic effect,

* $P_2$ describes the density of proliferative cells that are resistant to the first-line treatment $\mathcal{T}_1$ and sensitive to treatment $\mathcal{T}_2$,

* $P_3$ stands for proliferative cells that are resistant to both treatments.

It is worth noting that we do not aim at describing the evolution of the tumor from the very beginning of the GIST cancer, but we only focus on the evolution of the metastasis located at the liver. Therefore, according to the clinical observations, it seems relevant to consider that the three cell subpopulations are present when the GIST metastasis is detected.

<table><caption>**Table 1.** List of quantities computed by the model – $d = \text{day}$</caption><thead><tr><th>Name</th><th>Meaning</th><th>Unit</th></tr></thead><tbody><tr><td>P<sub>1</sub>(t, x)</td><td>Fraction of cells that are both sensitive to treatments T<sub>1</sub> and T<sub>2</sub></td><td>-</td></tr><tr><td>P<sub>2</sub>(t, x)</td><td>Fraction of cells that are resistant to treatment T<sub>1</sub> and sensitive to treatment T<sub>2</sub></td><td>-</td></tr><tr><td>P<sub>3</sub>(t, x)</td><td>Fraction of cells that are both resistant to treatments T<sub>1</sub> and T<sub>2</sub></td><td>-</td></tr><tr><td>N(t, x)</td><td>Fraction of necrotic cells</td><td>-</td></tr><tr><td>S(t, x)</td><td>Fraction of healthy cells</td><td>-</td></tr><tr><td>M(t, x)</td><td>Fraction of nutrients // Vascularization</td><td>-</td></tr><tr><td>ξ(t)</td><td>Average velocity of nutrients transport in direction of the tumor</td><td>cm.d<sup>-1</sup></td></tr><tr><td>v(t, x)</td><td>Velocity of the passive movement of the tumor under the pressure</td><td>cm.d<sup>-1</sup></td></tr><tr><td>Π(t, x)</td><td>Medium pressure<sup>1</sup></td><td>kg.cm<sup>-1</sup>.d<sup>-2</sup></td></tr></tbody></table>

<footnote n="1">¹ The mass unit in the pressure $\Pi$ and in the permeability $k$ has no importance, since only $k\nabla\Pi$ is relevant, and this term is homogeneous to $cm.d^{-1}$. Thus $k$ and $\Pi$ have just to be in the same arbitrary mass unit.</footnote>
## 2.1 The mathematical model

Let $\Omega$ be the domain where the tumor grows, and $\partial\Omega$ denotes its boundary. Note that the following model is valid for tumors that never reach the boundary $\partial\Omega$.

### 2.1.1 PDEs on proliferative, necrotic and healthy cells

The proliferative cell densities are driven by the following transport equations:

$$ \partial_t P_1 + \nabla \cdot (\mathbf{v} P_1) = (\gamma_{pp}(M) - \gamma_{pd}(M)) P_1 - (\mu_1 \chi_1(t) + \mu_2 \chi_2(t))(1+M)P_1 \quad \text{in } \Omega, \quad (2.1) $$

$$ \partial_t P_2 + \nabla \cdot (\mathbf{v} P_2) = (\gamma_{pp}(M) - \gamma_{pd}(M)) P_2 - \mu_2 \chi_2(t)(1+M)P_2 \quad \text{in } \Omega, \quad (2.2) $$

$$ \partial_t P_3 + \nabla \cdot (\mathbf{v} P_3) = (\gamma_{pp}(M) - \gamma_{pd}(M)) P_3 \quad \text{in } \Omega, \quad (2.3) $$

where $\chi_i(t) = 1_{[T_{\text{ini}}^i, T_{\text{end}}^i]}(t)$ is the time-characteristic functions of treatment $T_i$ and $\mu_i$ stands for the rate of the death<fnref n="2" /> due to $T_i$ of the proliferative cells, for $i \in \{1, 2\}$. The term $\mathbf{v}(t, \mathbf{x})$ denotes the velocity due to the tumor area changes and $M(t, \mathbf{x})$ stands for the vascularization and the nutrient supply. The rate of proliferation (resp. death) of P-cells, denoted by $\gamma_{pp}$ (resp. $\gamma_{pd}$), depends on M as follows:

$$ \gamma_{pp}(M) = \gamma_0 \frac{1 + \tanh(R(M - M_{th}))}{2}, \quad (2.4) $$

$$ \gamma_{pd}(M) = \gamma_1 \frac{1 - \tanh(R(M - M_{th}))}{2}, \quad (2.5) $$

where R is a numerical smoothing parameter<fnref n="3" />, $\gamma_0$ and $\gamma_1$ are respectively proliferative/decay parameters and $M_{th}$ is the hypoxia threshold.

We assume that healthy cells are only sensitive to hypoxia, and they are passively transported by the tumor area changes:

$$ \partial_t S + \nabla \cdot (\mathbf{v} S) = -\gamma_{sd}(M) S, \quad (2.6) $$

where $\gamma_{sd}$ is the rate of death of healthy cells due to hypoxia:

$$ \gamma_{sd}(M) = C_S \gamma_1 \max(0, -\tanh(R(M - M_{th}))). \quad (2.7) $$

Note that $\gamma_{sd}$ vanishes exactly if $M \ge M_{th}$ in order to ensure that the $S = 1$ on the outer boundary at any time. Finally, necrotic cell density satisfies

$$ \partial_t N + \nabla \cdot (\mathbf{v} N) = \gamma_{pd}(M) P + \gamma_{sd}(M) S + (\mu_1 \chi_1(t) P_1 + \mu_2 \chi_2(t) (P_1 + P_2))(1+M) - \delta(1+M)N, \quad (2.8) $$

where

$$ P = P_1 + P_2 + P_3, \quad (2.9) $$

and $\delta$ is a parameter that controls the elimination rate of the necrotics cells by the immune system.

The following Dirichlet conditions are used on the boundary if the velocity is incoming:

$$ P_1 = P_2 = P_3 = N = 1 - S = 0 \quad \text{for } \mathbf{x} \in \partial\Omega, \text{ if } \mathbf{v} \cdot \mathbf{n} < 0, \quad (2.10) $$

where $\mathbf{n}$ is the outgoing normal vector of the domain $\Omega$.

<footnote n="2">²The death rate due to the treatment is clearly linked to the dose of drug delivered to the patient, but not only. For instance, the sensitivity of the patient and the dose that really reaches the tumor are also involved.</footnote>

<footnote n="3">³Note that the functions $\gamma_{pp}$ and $\gamma_{pd}$ are merely regularized Heaviside functions. For the numerical simulations, we arbitrarily set R to 5.</footnote>
## 2.1.2 Mechanics

By using the following saturation condition (as in [1])

$$P + N + S = 1, \qquad (2.11)$$

and summing (2.1), (2.2), (2.3), (2.6) and (2.8), we obtain

$$\nabla \cdot \mathbf{v} = \gamma_{pp} P - \delta(1 + M)N. \qquad (2.12)$$

Darcy's law ensures the solvability of the system, similarly to [17]:

$$\left\{ \begin{array}{ll} \mathbf{v}(t, \mathbf{x}) = -k \nabla \Pi(t, \mathbf{x}) & \text{in } \Omega, \\ \Pi(t, \mathbf{x}) = 0 & \text{on } \partial \Omega, \end{array} \right. \qquad (2.13)$$

where $\Pi$ is the pressure (or potential) of the medium and $k$ its permeability. The homogeneous Dirichlet condition is used since we consider that the domain of interest is not isolated, and the outer medium does not impose a pressure on the tumor. This assumption is valid for small tumors that are not mechanically constrained by the extratumoral region.

## 2.1.3 Vascularization, Nutrient Supply and Angiogenesis

It remains to describe the vascularization/nutrient supply $M$ and the impact of treatment $\mathcal{T}_2$ on it. It is worth noting that the second-line treatment does not impact directly $M$, but it blocks the production of growth factors that drive the quantity $M$.

We thus introduce a scalar variable $\xi$, which is related to the mean concentration of growth factors. It has been reported by [13] that hypoxic cells increase their production of growth factors, while highly proliferative cells do not need additional nutrient supply. Therefore, if $M$ is below $M_{th}$ then $\xi$ should increase. Note also that the anti-angiogenic effect of treatment $\mathcal{T}_2$ decreases the production of $\xi$, but only for the cells $P_1$ and $P_2$ since $P_3$ is the density of cells that are sensitive neither to $\mathcal{T}_1$ nor to $\mathcal{T}_2$. We thus describe the evolution of $\xi$ as

$$\partial_t \xi = \alpha \int_{\Omega} \left(1 + \epsilon_{\xi} - \gamma_{pp}(M)/\gamma_0\right) \left( (1 - \nu_2 \chi_2(t))(P_1 + P_2) + P_3 \right) d\mathbf{x} - \lambda \xi. \qquad (2.14)$$

The dimensionless parameter $\nu_2 \in (0, 1)$ stands for the anti-angiogenic effect of $\mathcal{T}_2$, assumed to be similar for $P_1$ and $P_2$, while $\epsilon_{\xi}$ reflects the ground production of growth factors by cancer cells.

The quantity $M$ is then driven by the following equation:

$$\left\{ \begin{array}{ll} \partial_t M - \xi \frac{\nabla S}{\|\nabla S\|} \nabla M = C_0 S \left(1 - \frac{M}{2M_{th}}\right) - \eta PM + \psi \Delta M & \text{in } \Omega, \\ M(t, \mathbf{x}) = 2M_{th} & \text{on } \partial \Omega, \end{array} \right. \qquad (2.15)$$

where $C_0$ is the angiogenic capacity of healthy cells, $\eta$ denotes for the destruction of the vascularization by proliferative cells, and $\psi$ is a diffusion parameter. The diffusive term describes the infiltration of blood vessels into the tumor. From the numerical point of view, this term has a regularizing effect on the vascularization $M$, and thus stabilizes the numerical scheme.

Note that if we have initially

$$0 \le \|M|_{t=0}\|_{L^{\infty}} \le 2M_{th} \qquad (2.16)$$

then at any time $t$, $0 \le M(t) \le 2M_{th}$. This reflects the fact that the healthy tissue surrounding the tumor is well-vascularized and supplied with enough nutrients.
## 2.2 Vector formulation of the equations on cancer cells

To simplify the notations, let us denote the vector of cell populations by $\mathbf{W} = {}^t(P_1, P_2, P_3, N)$. We define the vector-function

$$ \mathbf{G}(M, \mathbf{W}) = \begin{pmatrix} (\gamma_{pp}(M) - \gamma_{pd}(M) - (\mu_1\chi_1(t) + \mu_2\chi_2(t))(1+M))P_1 \\ (\gamma_{pp}(M) - \gamma_{pd}(M) - \mu_2\chi_2(t)(1+M))P_2 \\ (\gamma_{pp}(M) - \gamma_{pd}(M))P_3 \\ \gamma_{pd}(M)P + \gamma_{sd}(M)(1-P-N) + (\mu_1\chi_1 P_1 + \mu_2\chi_2(P_1 + P_2))(1+M) - \delta(1+M)N \end{pmatrix}. $$

It is also convenient to define $F(M, \mathbf{W})$ as

$$ F(M, \mathbf{W}) := \gamma_{pp}(M) \sum_{i=1}^{3} W_i - \delta(1+M)W_4, \qquad (2.17) $$

so that the set of equations (2.1)-(2.3) and (2.8) and (2.12) closed by the Darcy law read

$$ \left\{ \begin{array}{ll} \partial_t \mathbf{W} + (\nabla \mathbf{W}) \cdot \mathbf{v} + \mathbf{W}(\nabla \cdot \mathbf{v}) = \mathbf{G}(M, \mathbf{W}) & \text{on } \Omega, \\ \mathbf{W} = 0 & \text{on } \partial\Omega, \text{ if } \mathbf{v}.\mathbf{n} < 0, \end{array} \right. \qquad (2.18) $$

and

$$ \left\{ \begin{array}{ll} -\nabla \cdot (k\nabla\Pi) = F(M, \mathbf{W}), & \text{in } \Omega, \\ \Pi(t, \mathbf{x}) = 0 & \text{on } \partial\Omega. \end{array} \right. \qquad (2.19) $$

Note that the density of healthy cells $S$ is then given by (2.11).

## 3 Numerical methods

We use a 2D cartesian staggered grid with a finite volume method. For the numerical calculations, the domain $\Omega$ is the rectangle $[0, L] \times [0, D]$. The domain is meshed by a cartesian grid with $N_x$ points along the x-axis and $N_y$ points along the y-axis.

The cancer cell densities are discretized at the center of the cell grid and the velocities are discretized at the middle of each edge as shown in Figure 3.

<figure><img src="image_12.png" /><figcaption>Figure 3. Discretization of unknown variables on one cell grid.</figcaption></figure>

Note that the equality (2.11) gives straightforwardly $S$, without solving equation (2.6). We split the computation into the following steps.
• Given<fnref n="4" /> $\mathbf{W}^n$, and $M^n$ at the time $t^n$, we infer $F^n = F(M^n, \mathbf{W}^n)$ and $\mathbf{G}^n = \mathbf{G}(M^n, \mathbf{W}^n)$.

• We first compute the pressure $\Pi^n$ solution to (2.19) with $F^n$ as right hand side, from which we infer the velocity $\mathbf{v}^n$ thanks to the equation (2.13).

• Then the new computation time $t^{n+1} = t^n + \Delta t$ is determined using the equation (3.17).

• Afterwards we compute $\mathbf{W}^{n+1}$ from (2.18), from which we infer $S^{n+1}$.

• We end by computing $\xi^{n+1}$ and $M^{n+1}$ thanks to (2.14)-(2.15).

Let present precisely the schemes used in the numerical simulations.

## 3.1 Computation of pressure and velocity

According to (2.12) and (2.13), the pressure $\Pi^n$ is given by

$$ \begin{cases} -\nabla \cdot (k\nabla\Pi^n) = F^n := F(M^n, \mathbf{W}^n) & \text{on } \Omega, \\ \Pi^n = 0 & \text{on } \partial\Omega, \end{cases} \qquad (3.1) $$

where F is defined by (2.17). We solve this equation thanks to a classical 5-points scheme. The component along $\mathbf{e}_x$ of the velocity $v_{i+1/2j}^{x,n}$ at $\mathbf{x}_{i+1/2j}$ (resp. the component along $\mathbf{e}_y$, $v_{i,j+1/2}^{y,n}$ at the point $\mathbf{x}_{i,j+1/2}$) are given thanks to

$$ v_{i+1/2,j}^{x,n} = -k \frac{\Pi_{i+1,j}^n - \Pi_{i,j}^n}{\Delta x}, \quad v_{i,j+1/2}^{y,n} = -k \frac{\Pi_{i,j+1}^n - \Pi_{i,j}^n}{\Delta y} $$

and the velocity $\mathbf{v}_{ij}^n$ at the point $\mathbf{x}_{ij}$ is approached by

$$ \mathbf{v}_{ij}^n = v_{ij}^{x,n} \mathbf{e}_x + v_{ij}^{y,n} \mathbf{e}_y = \frac{1}{2} \left( v_{i+1/2,j}^{x,n} + v_{i-1/2,j}^{x,n} \right) \mathbf{e}_x + \frac{1}{2} \left( v_{i,j+1/2}^{y,n} + v_{i,j-1/2}^{y,n} \right) \mathbf{e}_y. \quad (3.2) $$

## 3.2 Advection equation

By definition of F given by (2.17), the equation (2.12) leads to $\nabla \cdot \mathbf{v} = F$. Thus, the equation (2.18) can read also as a non-conservative form

$$ \partial_t \mathbf{W} + (\nabla \mathbf{W}) \cdot \mathbf{v} = G(M, \mathbf{W}) - (\nabla \cdot \mathbf{v}) \mathbf{W} = G(M, \mathbf{W}) - F(M, \mathbf{W}) \mathbf{W}, \quad (3.3) $$

solved thanks to the following time-splitting scheme

$$ \frac{\mathbf{W}^* - \mathbf{W}^n}{\Delta t/2} = \mathbf{G}^n - F^n \mathbf{W}^n, \quad (3.4) $$

$$ \frac{\mathbf{W}^\# - \mathbf{W}^*}{\Delta t} + (\nabla \mathbf{W}^*) \cdot \mathbf{v}^n = 0, \quad (3.5) $$

$$ \frac{\mathbf{W}^{n+1} - \mathbf{W}^\#}{\Delta t/2} = \mathbf{G}^n - F^n \mathbf{W}^n. \quad (3.6) $$

A WENO5 type method as given by [14] is used to approach the gradient $\nabla \mathbf{W}$ involved in (3.5).

<footnote n="4">⁴The superscript $n$ stands for the discrete time $t^n$ of the quantity (for instance $S^n$ is the density of healty cells at the time $t^n$)</footnote>
## 3.3 Computation of the quantities $\xi$ and $M$

The scalar equation (2.14) is solved thanks to the classical forward Euler method:

$$ \frac{\xi^{n+1} - \xi^n}{\Delta t} = \alpha \int_{\Omega} \left(1 + \epsilon\xi - \gamma_{pp}(M^n)/\gamma_0\right) \left((P_1^n + P_2^n)(1 - \nu_2\chi_2(t^n)) + P_3^n\right) dx - \lambda\xi^n, \quad (3.7) $$

where we use a standard rectangle rule to approach the integral in the right hand side.

In order to solve the advection-diffusion equation (2.15) satisfied by $M$, we use the following time-splitting schemes:

$$ \frac{M^* - M^n}{\Delta t/2} = C_0 S^n \left(1 - \frac{M^n}{2M_{th}}\right) - \eta P^n M^n, \quad (3.8) $$

$$ \frac{M^n - M^*}{\Delta t} - \psi \frac{\Delta M^n + \Delta M^*}{2} = \xi^n \frac{\nabla S^n}{\|\nabla S^n\|} \nabla M^n, \quad (3.9) $$

$$ \frac{M^{n+1} - M^n}{\Delta t/2} = C_0 S^n \left(1 - \frac{M^n}{2M_{th}}\right) - \eta P^n M^n. \quad (3.10) $$

Equation (3.9) is computed as a heat diffusion equation with a standard 5-points scheme on the grid. The right hand side is approached by a WENO5 type scheme.

## 3.4 Modified WENO5 scheme

The standard WENO5 scheme as given by [14] is accurate in most of the cases, however some sets of parameters<fnref n="5" /> make us face numerical instabilities. More precisely, starting from irrotational initial data, the simulation can generate a clover-like structure, as reported by Figure 4, while the circular shape should be preserved.

Note that such errors on the shape provide error on the evolution of the area of the lesion, since the clover shape increases the contact surface and thus modifies the interaction between the vascularization and the tumor. These instabilities have to be fixed. The problem is due to the WENO5 stencil that tends to favor the directions of the grid where changes in the velocity direction occur. As we can see on Figure 4d, on the center of the tumor (around $x = 5$ cm in Figure 4d), there is a compression point: the velocity has centripetal directions, since $v_x$ is positive on the right and negative on the left. Moreover, around 1.5 cm from the center of the tumor, there is a rim of proliferative cells that induce a spreading movement: the velocity is centripetal close to the center but centrifugal far from it.

More precisely, for the standard WENO5 scheme, at any point $x_{ij}$ of the grid, the numerical approximation $\mathbf{W}_{ij}^{n+1}$ to equation (3.5) at the time $t^{n+1}$ is given by

$$ \mathbf{W}_{i,j}^{n+1} = \mathbf{W}_{i,j}^{n} + \Delta t \left( v_{i,j}^{x,n} \mathcal{F}(\Delta x, (\mathbf{W}_{i+k,j})_{k=-3, \dots, 3}) + v_{i,j}^{y,n} \mathcal{F}(\Delta y, (\mathbf{W}_{i,j+k})_{k=-3, \dots, 3}) \right), \quad (3.11) $$

where $v_{i,j}^{x,n}$ and $v_{i,j}^{y,n}$ are defined by (3.2) and where $\mathcal{F}$ is the WENO5 functional given by [14]. In order to avoid the numerical instabilities, we introduce the following twin-WENO5 scheme, which is a combination of the standard WENO5 stencil and a rotation at the angle $\alpha$ of the WENO5 stencil (see Figure 5), where $\alpha$ is defined by the grid steps $\Delta x$ and $\Delta y$ as

$$ \alpha = \arctan(\Delta y/\Delta x) \in (0, \pi/2). $$

<footnote n="5">⁵The set of parameters have been found incidently by fitting the tumor area evolution of patient B, see Section 5.1.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Numerical simulations with the standard WENO5 stencil for the specific set of parameters. Starting from circular initial data, clover-like structures appear.</figcaption></figure>

We introduce the coefficients ($v_{i,j}^{r,n}$, $v_{i,j}^{\theta,n}$) and $\Delta r$ defined by<fnref n="6" />

$$ \begin{pmatrix} v_{i,j}^{r,n} \\ v_{i,j}^{\theta,n} \end{pmatrix} = \begin{pmatrix} \cos \alpha & \sin \alpha \\ -\sin \alpha & \cos \alpha \end{pmatrix} \begin{pmatrix} v_{i,j}^{x,n} \\ v_{i,j}^{y,n} \end{pmatrix}, \quad \Delta r = \sqrt{\Delta x^2 + \Delta y^2}, \qquad (3.12) $$

and we discretize equation (3.5) thanks to our twin-WENO5 scheme:

$$ \mathbf{W}_{i,j}^{n+1} = \mathbf{W}_{i,j}^{n} + (1-\beta)\Delta t \left( v_{i,j}^{x,n} \mathcal{F}(\Delta x, (\mathbf{W}_{i+k,j})_{k=-3, \dots, 3}) + v_{i,j}^{y,n} \mathcal{F}(\Delta y, (\mathbf{W}_{i,j+k})_{k=-3, \dots, 3}) \right) \\ + \beta \Delta t \left( v_{i,j}^{r,n} \mathcal{F}(\Delta r, (\mathbf{W}_{i+k,j+k})_{k=-3, \dots, 3}) + v_{i,j}^{\theta,n} \mathcal{F}(\Delta r, (\mathbf{W}_{i-k,j+k})_{k=-3, \dots, 3}) \right), \qquad (3.13) $$

where $\beta \in (0, 1)$ is a numerical parameter that has to be chosen. In particular, the standard WENO5 scheme holds for $\beta = 0$. As we can see on Figure 6, our new scheme keeps the irrotational property in cases where WENO5 does not.

$$ \mathbf{v}_{i,j}^{n} = v_{i,j}^{x,n} \mathbf{e}_{x} + v_{i,j}^{y,n} \mathbf{e}_{y} = v_{i,j}^{r,n} \mathbf{e}_{r} + v_{i,j}^{\theta,n} \mathbf{e}_{\theta}, \quad \text{with} \quad \mathbf{e}_{r} = \cos \alpha \mathbf{e}_{x} + \sin \alpha \mathbf{e}_{y}, \quad \mathbf{e}_{\theta} = -\sin \alpha \mathbf{e}_{x} + \cos \alpha \mathbf{e}_{y}. $$

<footnote n="6">⁶The coefficients $v_{i,j}^{r,n}$ and $v_{i,j}^{\theta,n}$ are defined such that</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 5. Stencil of the twin-WENO5 scheme for a uniform grid (left) and a non-uniform grid (right).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 6. Numerical simulations with the twin-WENO5 scheme ($\beta = 0.26$). Comparing in Figure 4, the conservation of irrotational invariance is very clearly improved.</figcaption></figure>

## 3.5 CFL condition

In addition, CFL-type restriction condition is required to preserve numerical stability. First, the WENO5 type scheme leads to a CFL condition that writes :

$$ \Delta t < \min \left( \frac{\Delta x}{\max |v_x|}, \frac{\Delta y}{\max |v_y|}, \frac{\min(\Delta x, \Delta y)}{\xi} \right) := \Delta t_{adv}. \qquad (3.14) $$

The forward Euler scheme on equation (3.4) leads to the following inequality (coordinate by coordinate)

$$ \mathbf{W}^* = \mathbf{W}^n + \frac{\Delta t}{2}(\mathbf{G}^n - \mathbf{F}^n \mathbf{W}^n) \geq \left(1 + \frac{\Delta t}{2}(\bar{\mathbf{G}}^n - \mathbf{F}^n \mathbf{1})\right) \circ \mathbf{W}^n, $$
where $\mathbf{1} = t(1, 1, 1, 1)$, the symbol ◦ stands for the Hadamard product (the pointwise product of the vectors, coordinate by coordinate) and where $\bar{\mathbf{G}}^n$ reads:

$$ \bar{\mathbf{G}}^n = \begin{pmatrix} \gamma_{pp}(M^n) - \gamma_{pd}(M^n) - (\mu_1\chi_1^n + \mu_2\chi_2^n)(1 + M^n) \\ \gamma_{pp}(M^n) - \gamma_{pd}(M^n) - \mu_2\chi_2^n(1 + M^n) \\ \gamma_{pp}(M^n) - \gamma_{pd}(M^n) \\ -\delta(1 + M^n) \end{pmatrix}. $$

A similar relation between $\mathbf{W}^\#$ and $\mathbf{W}^{n+1}$ can be read from equation (3.6). Thus, assuming that $\mathbf{W}^n \ge 0$ and $S^n = 1 - \sum_i W_i^n \ge 0$ at the time $t^n$, the following restriction of the time step

$$ \Delta t < \min \left( \frac{1}{\max_i \|\bar{\mathbf{G}}_i^n - \mathbf{F}^n\|_{\infty}}, \frac{1}{\|\gamma_{sd}(M^n) + \mathbf{F}^n\|_{\infty}} \right) := \Delta t_W. \quad (3.15) $$

ensures that $\mathbf{W}^{n+1} \ge 0$ and $S^{n+1} \ge 0$. Similarly, the forward Euler scheme in equation (3.7) and in equation (3.10) on the vascularization, leads to

$$ \Delta t < \min \left( \frac{1}{\eta}, \frac{1}{\lambda} \right) := \Delta t_{angio}. \quad (3.16) $$

Finally, since the velocity might be very small, to prevent too large $\Delta t$, we arbitrarily choose a velocity $\bar{\mathbf{v}}$ and our CFL condition reads<fnref n="7" />

$$ \Delta t = C_{CFL} \min \left( \frac{1}{\bar{\mathbf{v}}} \min(\Delta x, \Delta y), \Delta t_{adv}, \Delta t_W, \Delta t_{angio} \right), \quad (3.17) $$

for a given constant $C_{CFL} < 1$.

# 4 Numerical Results

## 4.1 Numerical tools to compare the results with the CT-scans

In order to compare the numerical results to the CT-scans, we have to define the appropriate quantities of interest as well as to develop a numerical tool that reproduces the grey scale.

### 4.1.1 Numerical determination of the tumor area, the necrotic part and the tumor mass

Let the threshold $\epsilon_{th}$ be the minimal fraction of tumor cells above which we define numerically the tumor. The tumor area, that is numerically measured reads

$$ \mathcal{A}(t) = \int_{\Omega} \mathbf{1}_{\{\mathbf{x}: P(t,\mathbf{x}) + N(t,\mathbf{x}) > \epsilon_{th}\}}(\mathbf{x}) \, d\mathbf{x}. \quad (4.1) $$

We also define the area of each cancer cell population as

$$ \mathcal{A}_J(t) = \int_{\Omega} \mathbf{1}_{\{\mathbf{x}: J(t,\mathbf{x}) > \epsilon_{th}\}}(\mathbf{x}) \, d\mathbf{x}, \quad \text{for } J \in \{P_1, P_2, P_3, N\}. \quad (4.2) $$

We also define the mass of each population, and the total proliferative mass at any time as

$$ \mathcal{M}_J(t) = \int_{\Omega} J(t, \mathbf{x}) \, d\mathbf{x}, \quad J \in \{P_1, P_2, P_3, N\}, \quad (4.3) $$

$$ \mathcal{M}(t) = \int_{\Omega} P(t, \mathbf{x}) \, d\mathbf{x}. \quad (4.4) $$

<footnote n="7">⁷In the simulations, we choose $\bar{\mathbf{v}} = 1 \text{ cm/month}$, considering 30 days in one month and $C_{CFL} = 0.4$.</footnote>
### 4.1.2 Numerical reconstitution of CT-scans

The comparison of the numerical simulations with the CT-scans raises an important difficulty. Actually, unlike the numerics that provide the numerical values at any point of each quantity S, P and N, the CT-scans provide grey-levels that are related to local radiodensities thanks to the Hounsfield scale. Since Hounsfield unit (HU) makes it possible to quantify the tumor area and to detect its location on CT-scans, we introduce a numerical HU, which is a linear combination of the numerical results. More precisely, we consider a linear grey scale ranging from black to white. To each species (P-cells, healthy cells and necrotic cells), we allocate a coefficient $\tau_P$, $\tau_S$ and $\tau_N$, we then plot the quantity

$$ \tau_P P + \tau_N N + \tau_S S, \qquad (4.5) $$

which is a kind of numerical grey-level. Since for abdominal CT-scans, the Hounsfield scale is limited from -200 to +200, we arbitrarily fix the above coefficients to

$$ \tau_P = 60, \quad \tau_S = 120, \quad \tau_N = -140, $$

and we set the value -200 to black color and +200 to white color.

## 4.2 Extensive study of one specific patient

### 4.2.1 Comparison of the numerical results with to clinical data

We focus on patient A for whom we have the whole clinical protocol, as well as the clinical data of the tumor area evolution and a sequence of CT-scans. The numerical simulations are performed in a square of side $L = D = 6$ cm with 120 points in each direction. The time step $\Delta t$ is computed by using equation (3.17).

We choose the numerical parameters in order to reproduce the evolution of the tumor area. In particular, we did not try to fit with the images, we only verify that the spatial evolution is plausible compared to the images. The parameters are summarized in Table 2.

The numerical tumor area is then compared to the measured areas on Figure 1g (the circles represent the real data and the solid lines represent the numerical simulation). It is worth noting that according to this figure, the evolution of the tumor area is well reproduced. We emphasize that the initial conditions are crucial in the tumor growth. In order to match qualitatively with the shape of the lesion at the initial time (see Figure 1a), the initial condition is chosen as a perturbed ellipse. More precisely, given 3 parameters $r_1$, $r_2$ and $e$, we define in the domain $\Omega$, $d(\mathbf{x})$ as

$$ d(\mathbf{x}) = \sqrt{\left(\frac{x - L/2}{e(x - L/2) + cr_1}\right)^2 + \left(\frac{y - D/2}{e(y - D/2) + cr_2}\right)^2}, \qquad (4.6) $$

with $\mathbf{x} = (x, y)$ and where $c = \frac{2\pi}{2\pi - \arccos(1 - 2\epsilon_{th})}$. We then use the function Y defined by<fnref n="8" />

$$ Y(\mathbf{x}) = \begin{cases} 1 & \text{if } d(\mathbf{x}) \le 0.5, \\ 0 & \text{if } d(\mathbf{x}) \ge 1, \\ \frac{1}{2}(1 - \cos(2\pi d(\mathbf{x}))) & \text{else,} \end{cases} \qquad (4.7) $$

in order to impose the initial conditions:

$$ P_1|_{t=0} = (1 - \Sigma_{ini})Y, \quad P_2|_{t=0} = \frac{\Sigma_{ini}}{1 + q_{ini}}Y, \quad P_3|_{t=0} = q_{ini}P_2(t=0), \quad N|_{t=0} = 0, \qquad (4.8) $$

$$ M|_{t=0} = 2M_{th}, $$

<footnote n="8">⁸Note that if $e = 0$, then Y is rotationally invariant.</footnote>
<table><caption>**Table 2.** List of parameters of the models and their values for the two patients considered – *d* = day</caption><thead><tr><th>Name</th><th>Meaning</th><th>Unit</th><th>Value for fit patient A (Figure 1g)</th><th>Value for fit patient B (Figure 10g)</th></tr></thead><tbody><tr><td>γ<sub>0</sub></td><td>Tumor cells growth rate</td><td>d<sup>-1</sup></td><td>2.0e-2</td><td>6.33e-3</td></tr><tr><td>γ<sub>1</sub></td><td>Tumor cells apoptosis rate</td><td>d<sup>-1</sup></td><td>8.0e-3</td><td>4.46e-2</td></tr><tr><td>C<sub>S</sub></td><td>Healthy tissue apoptosis rate compared to γ<sub>1</sub></td><td>-</td><td>10</td><td>10</td></tr><tr><td>M<sub>th</sub></td><td>Hypoxia threshold</td><td>-</td><td>2</td><td>2</td></tr><tr><td>δ</td><td>Elimination rate of the necrotic tissue by the immune system</td><td>d<sup>-1</sup></td><td>1.33e-2</td><td>8.19e-2</td></tr><tr><td>ψ</td><td>Diffusion rate of the oxygen</td><td>cm<sup>2</sup>.d<sup>-1</sup></td><td>1.33e-2</td><td>3.33e-3</td></tr><tr><td>η</td><td>Consumption rate of tumor cells</td><td>d<sup>-1</sup></td><td>6.67e-2</td><td>8.05e-3</td></tr><tr><td>α</td><td>Angiogenic excitability</td><td>cm<sup>-1</sup>.d<sup>-2</sup></td><td>1.11e-3</td><td>8.0e-3</td></tr><tr><td>λ</td><td>Elimination rate of angiogenic growth factor signal</td><td>d<sup>-1</sup></td><td>2.0e-2</td><td>0.68</td></tr><tr><td>C<sub>0</sub></td><td>Angiogenic capacity of healthy tissue</td><td>d<sup>-1</sup></td><td>3.33e-2</td><td>3.33e-2</td></tr><tr><td>k</td><td>Tissue permeability</td><td>kg<sup>-1</sup>.cm<sup>3</sup>.d</td><td>1</td><td>1</td></tr><tr><td>T<sup>1</sup><sub>ini</sub></td><td>Beginning (in days) of the treatment T<sub>1</sub> administration</td><td>d</td><td>119</td><td>0</td></tr><tr><td>T<sup>1</sup><sub>end</sub></td><td>Ending (in days) of the treatment T<sub>1</sub> administration</td><td>d</td><td>867</td><td>845</td></tr><tr><td>T<sup>2</sup><sub>ini</sub></td><td>Beginning (in days) of the treatment T<sub>2</sub> administration</td><td>d</td><td>867</td><td>1049</td></tr><tr><td>T<sup>2</sup><sub>end</sub></td><td>Ending (in days) of the treatment T<sub>2</sub> administration</td><td>d</td><td>1298</td><td>1600</td></tr><tr><td>μ<sub>1</sub></td><td>Proliferative cells death rate due to treatment T<sub>1</sub></td><td>d<sup>-1</sup></td><td>7.17e-3</td><td>3.45e-3</td></tr><tr><td>ν<sub>2</sub></td><td>Inhibition rate of the angiogenesis by treatment T<sub>2</sub></td><td>-</td><td>0.8</td><td>0.90</td></tr><tr><td>μ<sub>2</sub></td><td>Proliferative cells death rate due to treatment T<sub>2</sub></td><td>d<sup>-1</sup></td><td>4.27e-3</td><td>3.0e-4</td></tr><tr><td>ε<sub>th</sub></td><td>Minimal proportion of tumor cells that can be detected on scans – Minimal threshold for the numerical location of the tumor</td><td>-</td><td>1.0e-2</td><td>0.1</td></tr><tr><td>Σ<sub>ini</sub></td><td>Proportion of cells that are resistant to imatinib at the time t = 0 – Equivalent to (P<sub>2</sub> + P<sub>3</sub>)<sub>t=0</sub></td><td>-</td><td>3e-06</td><td>0.10</td></tr><tr><td>q<sub>ini</sub></td><td>Proportion of imatinib resistant cells that also resist to sunitinib at time t = 0 – Equivalent to (P<sub>3</sub>/P<sub>2</sub>)<sub>t=0</sub></td><td>-</td><td>7.5e-3</td><td>0.41</td></tr><tr><td>ξ<sub>ini</sub></td><td>Growth factor signal at time t = 0</td><td>cm.d<sup>-1</sup></td><td>3.33e-3</td><td>0</td></tr><tr><td>ε<sub>ξ</sub></td><td>Residual production of growth factor</td><td>-</td><td>0.1</td><td>0.1</td></tr><tr><td>L, D</td><td>Dimensions of the computational domain</td><td>cm</td><td>6</td><td>12</td></tr><tr><td>N<sub>x</sub>, N<sub>y</sub></td><td>Number of point for each dimension of the computational domain</td><td>-</td><td>120</td><td>132</td></tr><tr><td>r<sub>1</sub></td><td>Radius along x-axis of the initial condition</td><td>cm</td><td>0.47</td><td>0.5</td></tr><tr><td>r<sub>2</sub></td><td>Radius along y-axis of the initial condition</td><td>cm</td><td>0.36</td><td>0.5</td></tr><tr><td>e</td><td>Kind of eccentricity of the initial condition</td><td>-</td><td>0.35</td><td>0</td></tr><tr><td>β</td><td>Twin-WENO5 weight</td><td>-</td><td>0</td><td>0.3</td></tr></tbody></table>
where $\Sigma_{ini}$ denotes the proportion of cells that are resistant to treatment $\mathcal{T}_1$ and $q_{ini}$ is the proportion of $P_2$ cells that resist also to treatment $\mathcal{T}_2$. For patient A, we choose $e = 0.35$ and $r_1 = 0.47$; $r_2 = 0.36$.

<figure><img src="image_16.png" /><figcaption>**Figure 7.** Numerical simulations for patient A: spatial evolution of the lesion with numerical reconstitution of CT-scans.</figcaption></figure>

The spatial aspects of the numerical simulations presented in Figure 7, make appear the following facts:

i) During the phase without treatment, from Day 0 to Day 119, the tumor grows. Since necrotic and proliferative cells are present, the numerical tumor is heterogeneous, as reported by Figure 7a.

ii) Then treatment $\mathcal{T}_1$ is delivered from Day 119 to Day 867. $P_1$-cells are killed and necrotic cells become predominant. Due to the choice of the coefficient $\tau_N$, the numerical tumor becomes homogeneous and darker as shown by Figure 7b.

iii) The rebound of the proliferative activity just before the regrowth of the tumor at Day 776 is characterized by an increase in tumor heterogeneity: a proliferative rim appears and gradually fulfill the necrotic interior of the tumor as illustrated by Figure 7c. It is worth noting that treatment $\mathcal{T}_1$ is still delivered and thus resistant cells start to be predominant.

iv) Then treatment $\mathcal{T}_2$ is delivered from Day 867 to Day 1298. Once again, the necrotic population increases, and the numerical tumor is darker, as shown by Figure 7d-7e.

v) Finally, at Day 1116, new therapeutic failure is getting ready. It is characterized once again by a proliferative rim on the tumor boundary (see Figure 7f).
Interestingly, the numerical spatial results are in accordance with the data, at least up to the last relapse around Day 1000. Actually, the five steps that have been numerically observed are reported by CT-scans (see Figure 1a-1f).

After Day 1116, our model is probably no longer valid since the numerical structure of the tumor is very different from the CT-scan. Maybe other phenomena that are not accounted for by our model occur, such as micro-environment changes or cell mutations. Some interactions in the 3rd direction (not accounted for here) can be also involved.

We emphasize that the model seems to provide important information that clinicians cannot have access to with imaging devices. More precisely, Figure 8, which presents the area and mass evolution of each cell population, makes it possible to state that according to our model:

* During the first shrinkage of the tumor caused by treatment $T_1$, from Day 119 to Day 406, we observe that

i) Treatment $T_1$ kills $P_1$-cells that become necrotic tissue.

ii) The reduction of the area is due to the elimination of necrotic cells by the immune system.

iii) $P_2$ and $P_3$-cells that are not sensitive to treatment $T_1$ keep dividing.

* During few months, from Day 406 to Day 778, the tumor area continues to slowly decrease due to the death of $P_1$-cells, however $P_2$ and $P_3$ cells keep growing and replace progressively the eliminated necrotic cells. Even though the cellular activity of $P_2$ and $P_3$ is not affected by the treatment, this leads in a first time, to a stabilization of the tumor area, before a regrowth of the tumor at Day 778. Actually, when $P_2 + P_3$ becomes too high, the tumor growth reoccurs, governed by cells that are resistant to treatment $T_1$.

* During treatment $T_2$, from Day 867 to Day 1298, we can notice that

i) $P_1$-cells are still sensitive to the treatment.

ii) $P_2$-cells become necrotic since they are sensitive to treatment $T_2$.

iii) $P_3$-cells that are resistant to the two treatments keep growing.

It is worth noting that for each relapse, the proliferative activity occurs on the tumor boundary, where there are nutrients. Moreover, our model produces differences in the evolution of the tumor area $\mathcal{A}$ and the tumor mass evolution $\mathcal{M}$ given respectively by (4.1) and (4.4).

In particular, the mass of cancer cells decreases right after the drug delivery, while the tumor area decreases with a delay in each case. This delay may be due to the fact that the killed proliferative cells turn into the necrotic phase. Therefore the total area is still the same until the necrotic cells are eliminated by the immune system. In addition, our model distinguishes the effect on angiogenesis of the two treatments. In both cases the angiogenic signal decreases but for different reasons. Actually, while treatment $T_2$ inhibits directly the angiogenic signal, treatment $T_1$ kills $P_1$-cells, which implies indirectly a decrease of the production of this signal.

### 4.2.2 Numerical study of the influence of treatment $T_1$ efficacy

We focus now on the numerical study of different outcomes of treatment $T_1$, to investigate their influence on the tumor growth. We take the parameters of Table 2, except that we let vary $\mu_1$. Let us define two characteristic durations:

* $T_{PFS}$, which is the progression free survival time. It is the time duration for which the tumor is smaller (in term of area) than at the beginning of the treatment.
<figure><img src="image_17.png" /><figcaption>**Figure 8.** Evolution of the mass (integral of grey levels, in arbitrary unit) and area (mm²) of each cellular population and angiogenic signal evolution (cm.d⁻¹) given by the numerical simulations.</figcaption></figure>

• T<sub>double</sub>, which is the time duration for which the tumor area is below twice its initial value at the beginning of the treatment.

Note that this study is purely theoretical and cannot be used for treatment optimization since the whole evolution of the disease (including the relapse phase) is needed in order to obtain the parameters that are used for the simulation. It is therefore clear that our approach, for the time being, cannot lead to the determination of an optimal protocol, but this numerical investigation is a crucial step in addressing this important challenge.

In Figure 9a, we show the progression free survival time T<sub>PFS</sub> with respect to μ<sub>1</sub>. If μ<sub>1</sub> is below a threshold value μ<sub>th</sub> (μ<sub>th</sub> ~ 0.0047 for patient A), then the tumor growth is not stopped. For μ<sub>1</sub> above this threshold, T<sub>PFS</sub> increases rapidly and reaches a plateau, which means that it is not necessary to increase the dose, since it has no effect on T<sub>PFS</sub>.

In Figure 9b, we provide T<sub>double</sub>. As one can see, T<sub>double</sub> is not increased by the increase the dose above the threshold μ<sub>th</sub>, which means once again it is not necessary to increase the drug delivery above a certain threshold value.

In Figure 9c, we have represented the minimum size reached by the lesion with respect to the dose μ₁. Note that this curve is decreasing: therefore for high values of μ₁, the minimum of the lesion area is smaller. However, as shown in the Figure 9d, the relationship between the minimum size of the lesion and the doubling time is not monotonic. In particular, if the minimum size of the lesion is very small, then the doubling time can be smaller. This could be thought of as some Darwinian selection mechanism: P₁-cells, that are predominant when the lesion is detected, are killed faster by the treatment and therefore more room and more nutrients for P₂-cells are available for their growth. Thus, the control time becomes smaller. These curves show that there exists a threshold μth above which treatment T₁ is efficient. Beyond this threshold, the minimal tumor area decreases again, but the lifetime of the patient is not increased.
<figure><img src="image_18.png" /><figcaption>(a) Progression free survival time ($T_{PFS}$ in days) with respect to $\mu_1$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c) Minimal area reached ($A_{min}$ in mm²) with respect to $\mu_1$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) Time corresponding to the growth of tumor area by a factor 2 ($T_{double}$ in days) with respect to $\mu_1$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) Phase portrait.</figcaption></figure>

**Figure 9.** Efficacy of treatment $\mu_1$ for patient A. The star corresponds to the parameters used in Figure 1 for the fit of the tumor area.

# 5 Discussion

In the previous section, we have extensively studied our model on one specific patient, patient A, whose tumor lesion was followed up by a sequence of CT-scans. We found parameters that make it possible to compare qualitatively the images with the numerics, thanks to our numerical reconstitution of CT-scans, and we also fitted the tumor area. The lesion of this patient is interesting since it is confined in the inner of the liver, and thus the tumor evolution was not constrained by mechanics of the organ. However, in some cases, the tumor metastasis is close to the liver boundary. In such a case, there is no hope of providing quantitative results on the spatial evolution of the lesion, but in the next subsection, we show that the tumor area evolution can be well described.

## 5.1 Tumor area evolution for patient B

In this subsection, we focus on patient B, whose tumor evolution is quite different from patient A. Actually, the clinical protocol of this patient was the following:

i) Patient B was treated in a first time with a specific tyrosine kinase inhibitor (imatinib), which stabilizes the increase of the tumor area during more than 10 months before a relapse.

ii) The multi-targeted tyrosine kinase inhibitor sunitinib was started, but unlike for patient A, it was
totally inefficient.

iii) A third treatment, was then delivered. Sorafenib is another multi-targeted tyrosine kinase inhibitor, that has both antiangiogenic and cytotoxic effects. The tumor seemed to be sensitive to this drug during several months until the treatment failure occurs.

Unfortunately, as illustrated by Figure 10g, the tumor evolution was so quick between the CT-scan at Day 429 and the CT-scan at Day 845 that the metastasis has reached the boundary of the liver, and thus with our model we have no hope to provide a numerical tumor growth that would be spatially in accordance with the CT-scans, since mechanical effect of the liver membrane have to be accounted for. Therefore we focus on the tumor area as given by the clinicians, the main challenge being to capture such quick tumor growth.

<figure><img src="image_22.png" /><figcaption>(a) May 23, 2007 –
Day 0</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) July 25, 2008 –
Day 429</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) Sept 14, 2009 –
Day 845</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) April 06, 2010 –
Day 1049</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) Sept 28, 2010 –
Day 1224</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(f) May 20, 2011 –
Day 1458</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(g) Tumor area. Each point is the tumor area measured on CT-scans and the line stands for our numerical results.
The letter refers to the CT-scans present on the left. The symbol ⊕ stands for heterogeneous tumor, ● stands for rather homogeneous tumor and ○ stands for CT-scan on which it is difficult to detect the heterogeneous aspect.
The value of parameters used in the numerical simulation are given in Table 2.</figcaption></figure>

**Figure 10.** Spatial evolution of the liver metastasis of patient B on a series of CT-scans. The small round lesion in the top of Figure 10a is considered.

Since sunitinib is totally inefficient in this case, we consider that treatment $T_1$ is delivered from Day 0 to Day 845, while treatment $T_2$ consists in the sorefinib, delivered from Day 1049 to Day 1600. We find

<footnote n="8">⁸Do not confuse the metastasis with the gallbladder that are bigger on the first two CT-scans.</footnote>
parameters that make it possible to fit qualitatively the tumor area measured by the clinicians and our numerical tumor area. These parameters are given in Table 2. The initial data has been chosen rotationally invariant, in accordance with the data of Figure 10a. As we can see, the tumor evolution is very stiff: between Day 416 and Day 614 the tumor area has been multiplied by more than 9, and once sorafenib is delivered the tumor area decreases from 4 500 mm² to 2 850 mm². We manage, at least qualitatively, to obtain such behavior, even though the fit is not perfect.

We have also investigated the efficacy of the treatment $T_1$, as for patient A. It has been observed that the progression free survival time $T_{PFS}$, the doubling time $T_{double}$ and the minimum area reached by the lesion $A_{min}$ have the same profile as the patient studied in the previous section. In particular, there exists a threshold $\mu_{th}$ below which treatment has no effect on $T_{PFS}$. Then an increase of the dose does not improve the $T_{PFS}$ that reached, as for patient A, a plateau beyond $\mu_{th}$. The doubling time $T_{double}$ with respect to $\mu_1$ is also not monotonic, contrary to $A_{min}$. Thus, as previously, beyond the threshold $\mu_{th}$, the tumor area decreases again, but the overall survival time of the patient is not increased.

## 5.2 Consistency of the model

<figure><img src="image_29.png" /><figcaption>Figure 11. Different behaviors accounted for by the model.</figcaption></figure>

Our model reproduces the clinical data given for patient A and patient B. Moreover, it is possible to account for the following behaviors that have been reported by physicians as shown in the Figure 11 (see Table 3 in the supplementary informations section, for the different values of the used parameters):

a) The metastasis is controlled by the treatment $T_1$ (imatinib). In this case, there is no clinical need to change the treatment.

b) The metastasis is controlled by the treatment $T_1$ but then the tumor regrows. The treatment $T_2$ (sunitinib or sorafenib) is then delivered successfully, and the tumor area is controlled.

c) The metastasis is controlled by the treatment $T_1$ before a first relapse. Then, the treatment $T_2$ is efficient before a second relapse.

d) The metastasis is controlled by the treatment $T_1$ before a first relapse. Then, the treatment $T_2$ is totally inefficient.

e) The treatment $T_1$ is totally inefficient. Then, the treatment $T_2$ is efficient and the tumor area is controlled.

f) The treatment $T_1$ is totally inefficient. Then, the treatment $T_2$ is efficient before a relapse.
g) The treatment $T_1$ is totally inefficient. Then, the treatment $T_2$ is totally inefficient. This profile holds for a patient with a genetic mutation EXON, as reported by Andersson et al. [2], or Hirota et al. [12].

# 6 Conclusion

In this paper, we have provided a patient-dependent model, based on PDEs, that describes the global behavior of GIST metastasis to the liver during the different stages. We have presented the numerical methods used to solve the PDE system and we introduced a new WENO5 type scheme, called twin-WENO5. Then, this model has been numerically compared with clinical observations concerning patient A, who have been treated successively by imatinib and sunitinib. As presented by Figure 1g, our model provides results that are qualitatively in accordance with the clinical data. In particular, our model is able to describe not only the evolution of the size of the lesion, but also its structure, as illustrated by Figure 1 and 2.

Interestingly, it has been reported by our simulations that a rim of proliferative cells appears on the tumor boundary just before the relapse time. This seems to be corroborated in the CT-scans images by an increase of the tumor heterogeneity, in the sense of grey-level, before the regrowth. The more the metastases are heterogeneous, the quicker is the relapse. This result underlines the fact that the RECIST criteria is not sufficient to evaluate the efficiency of a treatment.

We also investigated numerically the effect of the parameter $\mu_1$, linked to the efficacy of treatment $T_1$. We have shown that, according to the numerical model, increasing the value of $\mu_1$, which can be seen as an increase in the drug delivery, does not provide better results in terms of progression free survival time as reported by Figure 9.

It is worth noting that our model fits well with the data, but it is not predictive. Indeed, Figure 11b-c-d shows that the knowledge of the first 400 days is not sufficient to determine uniquely the tumor growth. This means that more precise data such as functional imaging could be necessary for a better analysis of the inner structure of metastases.

In conclusion, we have provided a model that fits with the clinical CT-scans we had. The forthcoming work will consist in adding more biological information that cannot be obtained from the CT-scans, in order to provide a predictive model. Note that functional imaging data (TEP or MRI) or biopsies might be crucial in enriching the present model.

# Acknowledgement

This study has been carried out within the frame of the LABEX TRAIL, ANR-10-LABX-0057 with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx (ANR-10-IDEX-03-02).

Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr/)

The work of Patricio Cumsille (PC) was partially supported by Associative Research Program (PIA) from Conicyt under grant number FB0001. The work of PC was also partially supported by Universidad del Bío-Bío under grant DIUBB 121909 GI/C and DIUBB 122109 GI/EF.

The authors thank very warmly Dr H. Fathallah-Shaykh for discussions and advices that significantly contributed to improving the quality of the paper.
## References

1. D. Ambrosi and L. Preziosi. On the closure of mass balance models for tumor growth. *Mathematical Models and Methods in Applied Sciences*, 12(05):737-754, 2002.

2. J. Andersson, P. Bumming, J.M. Meis-Kindblom, H. Sihto, N. Nupponen, H. Joensuu, A. Oden, B. Gustavsson, L.G. Kindblom, and B. Nilsson. Gastrointestinal stromal tumors with kit exon 11 deletions are associated with poor prognosis. *Gastroenterology*, 130:1573-1581, 2006.

3. S. Benzekry, C. Lamont, A. Beheshti, A. Tracz, J.M.L. Ebos, L. Hlatky, and P. Hahnfeldt. Classical mathematical models for description and prediction of experimental tumor growth. *arXiv preprint arXiv:1406.1446*, 2014.

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

5. J.Y. Blay. A decade of tyrosine kinase inhibitor therapy: Historical and current perspectives on targeted therapy for GIST. *Cancer Treatment Reviews*, 37(5):373-384, 2011.

6. J.Y. Blay, A. Le Cesne, P.A. Cassier, and I.L. Ray-Coquard. Gastrointestinal stromal tumors (GIST): a rare entity, a tumor model for personalized therapy, and yet ten different molecular subtypes. *Discov Med*, 13(72):357-67, May 2012.

7. D. Bresch, T. Colin, E. Grenier, B. Ribba, and O. Saut. A viscoelastic model for avascular tumor growth. *Discrete And Continuous Dynamical Systems*, Volume 2009:101-108, 2009.

8. D. Bresch, T. Colin, E. Grenier, B. Ribba, and O. Saut. Computational modeling of solid tumor growth: The avascular stage. *SIAM Journal on Scientific Computing*, 32(4):2321-2344, 2010.

9. T. Colin, O. Gallinato, C. Poignard, and O. Saut. Tumor growth model for ductal carcinoma: from in situ phase to stroma invasion. Submitted Report RR-8502, INRIA, March 2014.

10. T. Colin, A. Iollo, D. Lombardi, and O. Saut. System identification in tumor growth modeling using semi-empirical eigenfunctions. *Mathematical Models and Methods in Applied Sciences*, 22(06):1250003, 2012.

11. Ronald P DeMatteo, Jonathan J Lewis, Denis Leung, Satvinder S Mudan, James M Woodruff, and Murray F Brennan. Two hundred gastrointestinal stromal tumors: recurrence patterns and prognostic factors for survival. *Annals of surgery*, 231(1):51, 2000.

12. S. Hirota, K. Isozaki, Y. Moriyama, K. Hashimoto, T. Nishida, S. Ishiguro, K. Kawano, M. Hanada, A. Kurata, G. Muhammad Tunio, Y. Matsuzawa, Kanakura Y., Y. Shinomura, and Y. Kitamura. Gain-of-function mutations of c-kit in human gastrointestinal stromal tumors. *Science*, 279(5350):577-80, Jan 1998.

13. K. Holmes, O. Ll. Roberts, A. M. Thomas, and M.J. Cross. Vascular endothelial growth factor receptor-2: Structure, function, intracellular signalling and therapeutic inhibition. *Cellular Signalling*, 19(10):2003-2012, 2007.

14. X.D. Liu, S. Osher, and Chan T. Weighted essentially non-oscillatory schemes. *Journal of Computational Physics*, 115(1):200-212, 1994.
15. A. Lorz, T. Lorenzi, J. Clairambault, A. Escargueil, and B. Perthame. Effects of space structure and combination therapies on phenotypic heterogeneity and drug resistance in solid tumors. *arXiv preprint arXiv:1312.6237*, 2013.

16. B. Nilsson, P. Bumming, J.M. Meis-Kindblom, A. Oden, A. Dortok, B. Gustavsson, K. Sablinska, and L.G. Kindblom. Gastrointestinal stromal tumors: the incidence, prevalence, clinical course, and prognostication in the preimatinib mesylate era—a population-based study in western sweden. *Cancer*, 103(4):821–9, Feb 2005.

17. B. Ribba, O. Saut, T. Colin, D. Bresch, E. Grenier, and J.P. Boissel. A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents. *Journal of Theoretical Biology*, 243(4):532–541, 2006.

18. O. Saut, J-B. Lagaert, T. Colin, and H. Fathallah-Shaykh. A multilayer grow-or-go model for GBM: Effects of invasive cells and anti-angiogenesis on growth. *Bulletin of Math. Biol.*, 76:2306–2333, 2014.

19. N. Schramm, E. Englhart, M. Schlemmer, M. Hittinger, C. Übleis, C.R. Becker, M.F. Reiser, and F. Berger. Tumor response and clinical outcome in metastatic gastrointestinal stromal tumors under sunitinib therapy: comparison of recist, choi and volumetric criteria. *European journal of radiology*, 82(6):951–958, 2013.

20. C. Suzuki, H. Jacobsson, T. Hatschek, M. R. Torkzad, K. Bodén, Y. Eriksson-Alm, E. Berg, H. Fujii, A. Kubo, and L. Blomqvist. Radiologic measurements of tumor response to treatment: Practical approaches and limitations. *Radiographics*, 28(2):329–344, 2008.

# 7 Supplementary informations

The set of parameters used to compute the numerical results given in Figure 11 of section 5.2 is presented in Table 3.
<table><caption>**Table 3.** Values of the parameters for curves presented in Figure 11.</caption><thead><tr><th>Name</th><th>a)</th><th>b)</th><th>c)</th><th>d)</th><th>e)</th><th>f)</th><th>g)</th></tr></thead><tbody><tr><td>γ<sub>0</sub></td><td>2.0e-2</td><td>2.03e-2</td><td>2.0e-2</td><td>1.97e-2</td><td>1.33e-2</td><td>1.33e-2</td><td>1.33e-2</td></tr><tr><td>γ<sub>1</sub></td><td>6.67e-3</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td></tr><tr><td>C<sub>S</sub></td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>M<sub>th</sub></td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><td>δ</td><td>2.67e-2</td><td>3.0e-2</td><td>5.0e-2</td><td>3.0e-2</td><td>3.0e-2</td><td>3.0e-2</td><td>3.0e-2</td></tr><tr><td>ψ</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td></tr><tr><td>η</td><td>6.67e-2</td><td>6.67e-2</td><td>6.67e-2</td><td>6.67e-2</td><td>6.67e-2</td><td>6.67e-2</td><td>6.67e-2</td></tr><tr><td>α</td><td>1.11e-3</td><td>1.11e-3</td><td>1.11e-3</td><td>1.11e-3</td><td>1.11e-3</td><td>1.11e-3</td><td>1.11e-3</td></tr><tr><td>λ</td><td>2.0e-2</td><td>2.0e-2</td><td>2.0e-2</td><td>2.0e-2</td><td>2.0e-2</td><td>2.0e-2</td><td>2.0e-2</td></tr><tr><td>C<sub>0</sub></td><td>3.33e-2</td><td>3.33e-2</td><td>3.33e-2</td><td>3.33e-2</td><td>3.33e-2</td><td>3.33e-2</td><td>3.33e-2</td></tr><tr><td>k</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>T<sup>1</sup><sub>ini</sub></td><td>119</td><td>119</td><td>119</td><td>119</td><td>119</td><td>119</td><td>119</td></tr><tr><td>T<sup>1</sup><sub>end</sub></td><td>3000</td><td>867</td><td>867</td><td>867</td><td>300</td><td>300</td><td>300</td></tr><tr><td>T<sup>2</sup><sub>ini</sub></td><td>3000</td><td>867</td><td>867</td><td>867</td><td>300</td><td>300</td><td>300</td></tr><tr><td>T<sup>2</sup><sub>end</sub></td><td>3000</td><td>1700</td><td>1298</td><td>1700</td><td>1700</td><td>1700</td><td>1700</td></tr><tr><td>μ<sub>1</sub></td><td>8.33e-3</td><td>8.33e-3</td><td>8.33e-3</td><td>8.33e-3</td><td>8.33e-3</td><td>8.33e-3</td><td>8.33e-3</td></tr><tr><td>ν<sub>2</sub></td><td>0.9</td><td>0.9</td><td>0.99</td><td>0.9</td><td>0.9</td><td>0.9</td><td>0.9</td></tr><tr><td>μ<sub>2</sub></td><td>6.0e-4</td><td>6.0e-4</td><td>6.6e-4</td><td>6.0e-4</td><td>6.0e-4</td><td>6.0e-4</td><td>6.0e-4</td></tr><tr><td>ε<sub>th</sub></td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td><td>1.0e-2</td></tr><tr><td>Σ<sub>ini</sub></td><td>0</td><td>4e-06</td><td>4e-06</td><td>2e-06</td><td>1</td><td>1</td><td>0.9</td></tr><tr><td>q<sub>ini</sub></td><td>0</td><td>0</td><td>4.5e-2</td><td>1</td><td>2e-07</td><td>3.0e-2</td><td>0.9</td></tr><tr><td>ξ<sub>ini</sub></td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td><td>3.33e-3</td></tr><tr><td>L, D</td><td>6</td><td>6</td><td>6</td><td>6</td><td>6</td><td>6</td><td>6</td></tr><tr><td>N<sub>x</sub>, N<sub>y</sub></td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td><td>120</td></tr><tr><td>r<sub>1</sub>, r<sub>2</sub></td><td>0.62</td><td>0.62</td><td>0.62</td><td>0.62</td><td>0.62</td><td>0.62</td><td>0.62</td></tr><tr><td>e</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>β</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>
informatics mathematics
Inria

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399