Thierry Colin

Olivier Gallinato

Clair Poignard
email: clair.poignard@inria.fr

Olivier Saut

Tumor Growth Model for Ductal Carcinoma: from in Situ Phase to Stroma Invasion

Keywords: tumor growth model, nonlinear advection equations, Kedem-Katchalsky conditions, breast cancer, metalloproteinases modèle de croissance tumorale, équations d'advection non-linéaires, conditions de Kedem-Katchalsky, cancer du sein, métalloprotéinases

We propose a tumor growth model, based on partial differential equations, for breast ductal carcinoma. In particular, we focus on the transition from the in situ stage to the invasive phase, under the action of specific enzymes: matrix metalloproteinases (MMPs), that are produced by cancer cells. Advection equations describe the movement of biological cells, generated by cancer cell proliferation that exerts pressure on the surrounding tissues. In its early phase, ductal carcinoma remains confined to the galactophoric duct due to the basement membrane. However, as soon as MMP enzymes are produced, they are able to degrade the membrane, which leads to cancer invasion. The specificity of the model lies in the description of these phenomena and the modeling of the duct membrane degradation, thanks to nonlinear Kedem-Katchalsky type transmission conditions to describe the pressure.

Introduction

Ductal carcinomas, that occur in galactophoric ducts, are among the most common breast cancers. Under certain circumstances, these cancers may become invasive and agressive by crossing the duct membrane. In this paper, we propose a tumor growth model for two types of ductal carcinomas: the ductal carcinoma in situ (DCIS) and the invasive ductal carcinoma (IDC). DCIS remains confined to the duct due to the duct membrane which yields a barrier. But if some cancer cells produce enzymes called MMPs which degrade the membrane, DCIS then turns into an IDC, which is much more agressive.

We propose a continuous model based on partial differential equations, which describe the tumor growth by considering the cell densities, in the same way as in [START_REF] Ribba | A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents[END_REF][START_REF] Bresch | A viscoelastic model for avascular tumor growth. Discrete and continuous dynamical systems[END_REF][START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF]. The tumor growth is driven by an advection equation which takes into account the fact that cancer cell proliferation generates a pressure that pushes the healthy tissue. The new feature of our model lies in the modeling of the duct membrane, which separates the lumen from the breast tissue and which is degraded by the MMPs.

Biological background

Breast consists essentially of galactophoric ducts, lobes and fatty tissue [START_REF] Ramsay | Anatomy of the lactating human breast redefined with ultrasound imaging[END_REF]. A scheme of breast is given at Fig. 1. The ducts are surrounded by a connective tissue, the stroma, and consist of a basement membrane and an epithelium surrounding the central opening in which flows the lumen (Fig. 1,A). Ductal carcinoma derives from an epithelial cell inside the duct. It is a dysfunction of the epithelial cells that proliferate uncontrollably. Its origin and subsequent developments combine genetic and environmental factors, including many underlying mechanisms, which still Inria remain poorly understood. The development of a precancerous lesion in epithelium, as hyperplasia (Fig. 1, B andC), increases the risk of abnormalities. In particular, the disruption of interactions with the membrane may cause an accumulation of abnormal events in epithelial tissue, leading to cancer initiation and genetic instability [START_REF] Weaver | The importance of the microenvironment in breast cancer progression: recapitulation of mammary tumorigenesis using a unique human mammary epithelial cell model and a three-dimensional culture assay[END_REF][START_REF] Bissell | Why don't we get more cancer ? A proposed role of the microenvironment in restraining cancer progression[END_REF]. Among the consequences of genetic mutations, the increase in mitosis rate, the loss of cell polarity, the breakage of adhesion junctions and the ability to evade apoptosis lead to uncontrolled growth of the tissue and to tumor formation inside the duct. As long as tumor growth remains confined to the duct, carcinoma is called in situ (Fig. 1,D). The key point is then to determine whether or not the tumor will pass through the membrane (Fig. 1, E andF), since it makes the difference between a low-grade cancer, rather non-aggressive, and a high-grade cancer, invasive, progressing to a vascular stage and leading to metastasis and invasion of the body. MMP enzymes have the property to degrade the extracellular matrix (ECM) of biological tissues, enabling to reshape them under normal circumstances. In the case of ductal carcinomas, the production of MMP enzymes leads to the degradation of the basement membrane ECM. The membrane becomes porous, which makes possible the cancer cell invasion. The issue of MMP production, which is complex and not well understood yet, is crucial in this study.

Outline and positionning

The aim of the paper is to build a model that can reproduce qualitatively the different configurations observed in the medical experiments and the biopsy results. In Section 2, we derive our new model, based on the biology of ductal carcinoma. New developments, compared to other existing tumor growth models, can be classified into two types:

1. Specific aspects of the model as a model for breast cancer:

• the mechanics related to the duct are taken into account, • the vascularization is located outside the duct, therefore the growth phase in situ is an avascular phase, • the nutrients diffusion is artificially forced within the duct, which avoids introducing a nonlinear diffusion to obtain realistic rates of nutrients.

Non-specific aspects:

• the pressure is taken into account in the proliferation rate, as some studies show that it tends to limit cancer cell proliferation [START_REF]Externally applied forces can phenotypically revert malignant breast epithelial structures[END_REF][START_REF] Yang | To revert breast cancer cells, give them the squeeze[END_REF],

• the very low thickness of the basement membrane led us to model it as an interface with jump conditions to link the outer and inner media thanks to Kedem-Katchalsky type conditions.

Section 3 is devoted to 2D numerical results from single-cell or multi-cell models, depending on whether the cancer is initiated from one single cell or more, inside the duct. The results show the membrane degradation under the action of MMPs, and the progression of an invasive carcinoma through the membrane. However, whatever the fit of the parameters, the results are not very realistic, too homogeneous and symmetrical (see Fig. 18 at Section 4). The model is then refined in Section 4 to compare some possible causes of heterogeneity:

• introduction of a mutation on certain cancer cells, which causes the MMP production, resulting in an anisotropic tumor invasion,

• modeling the role of cells of the microenvironment that produce MMPs or TNFs (tumor necrosis factors), with satisfactory results in terms of tumor heterogeneity.

In particular, the consideration of the microenvironment influence on the MMP production leads to one of the most significant result of the study, given by Fig. 2: numerical results (Fig. 2a) show the invasion of the breast tissue due to the membrane disruption, and are qualitatively similar to biopsies (Fig. 2b). Finally, we conclude by pointing out the pros and the cons of our modeling, and we raise some questions that are still opened in breast cancer modeling.

Tumor growth model

The structure of the model is outlined in Fig. 3. For the sake of clarity, we subdivide our model into three submodels. Submodel 1, presented in Subsection 2.2, is at the heart of the model: it describes the general movement of cell species caused by the pressure due to cell proliferation. The specificity lies in the modeling of the very thin membrane as an interface across which Kedem-Katchalsky type conditions are imposed to describe the pressure jump. Submodel 2, given at Subsection 2.3, describes the action of MMP enzymes, leading to cell invasion, which is crucial in this study. Finally, submodel 3, dealing with nutrients, will be discussed in Subsection 2.4. It considers that the supply of nutrients comes from the membrane vasculature and from the neo-vasculature resulting from the angiogenesis.

Notations and hypotheses

We focus on four cell species: at the initial time, healthy cells are outside the duct (stroma), the lumen and the tumoral tissues are inside. Initially there is no necrotic tissue. We denote Inria by P , N , S and L the respective fractions of proliferating, necrotic, healthy and luminal cell populations, such that, in the whole domain we have

P + N + S + L = 1 ∀t, x. (1)
These variables depend on the time t and the position x = (x, y). The total cell density, denoted by đ 0 , is assumed to be constant in time and in the whole domain. We assume that cellular species behave like an incompressible, passive and porous fluid, driven by a homogeneous movement.

The healthy tissue and the lumen cell activity are neglected. The basement membrane, whose porosity (or permeability) is very low when it is not degraded, is an insulating barrier between the interior of the duct and the breast tissue. We denote by Ω the considered domain composed of the duct (Ω 0 ), the membrane (Γ) and the outer medium (Ω 1 ).

Ω0

Ω1 Γ

∂Ω

Cell movement and pressure equations

The main steps of tumor growth are outlined in Fig. 5: after cancer initiation, tumor cells either turn into necrotic cells or proliferate, depending on the nutrient concentration. The growth of volume due to cellular division induces the movement of the whole tissue.

Advection

We assume that the cells move according to the same velocity field v, due to the pressure generated by the proliferation. Applying the principle of mass balance, we obtain that the general movement is governed by the following advection equations that hold in Ω:

∂ t P + ∇ • (vP ) = (α P -α N )P, ( 2a
)
∂ t N + ∇ • (vN ) = α N P, ( 2b
)
∂ t S + ∇ • (vS) = 0, ( 2c
)
∂ t L + ∇ • (vL) = 0, (2d)
where α P and α N denote the rates of proliferation and necrosis, respectively. The necrosis rate is a transition rate between proliferating and necrotic cells. The velocity field v arises from the division of proliferating tumor cells, at the rate α P , which are responsible for the volume growth. By summing the equations (2a)-(2d) and using the saturation condition (1), we obtain the divergence of the velocity in Ω:

∇ • v = α P P.
(3)

Proliferation and necrosis rates

The proliferation and necrosis rates are driven by genetic internal factors and by environmental factors: the nutrient concentration Θ and the pressure Π. Actually, among the environmental influences, the pressure exerted on tumor cells has a significant impact in reducing their rate of proliferation. This intuitive idea seems to be supported by the ongoing work about the role of the microenvironment on tumor growth, led by the teams of D.A. Fletcher and M.J. Bissell [START_REF]Externally applied forces can phenotypically revert malignant breast epithelial structures[END_REF][START_REF] Yang | To revert breast cancer cells, give them the squeeze[END_REF].

For the sake of simplicity, we shall generically call hypoxia the lack of nutrients. We therefore express the proliferation and necrosis rates by:

α P = α f P (4a) α N = α f N (4b)
where α is a generic rate which is weighted by the functions f P and f N . The functions f P and f N are dimensionless weight functions for proliferating and necrosis rates, respectively. They contain the genetic factors, the pressure, the nutrient influences, and are given by

f P (Θ, Π) = α G 1 + tanh [Λ P (Θ -Θ H )] 2 Π max -Π Π max , (5a) f N (Θ) = 1 -tanh [Λ N (Θ -Θ N )] 2 , ( 5b
)
where Inria • Π max is the maximal pressure that cancer cells can withstand before stopping their mitoses. Θ H and Θ N are the hypoxia and necrosis thresholds for nutrient concentration Θ. We impose 0 < Θ N < Θ H , such that for Θ smaller than Θ N cancer cells P turn into necrotic cells N , between Θ N and Θ H the cells P do not proliferate, but are not necrotic: they are in a quiescent phase.

Above Θ H the proliferating cells P proliferate and the tumor grows.

• Λ P and Λ N are strictly positive constants, homogeneous to the inverse of a concentration. They give the slopes of the curves representing f P and f N in the neighborhood of the thresholds Θ H and Θ N ,

• α G is a weight function whose value is between 0 and 1. It represents the quantitative influence of oncogenic factors, such as the overexpression of estrogen receptors or other growth factor receptors like HER-2 [START_REF] Osborne | Oncogenes and tumor suppressor genes in breast cancer: potential diagnostic and therapeutic applications[END_REF][START_REF] Thanoon | Computational framework for local breast cancer treatment[END_REF]. In this study, for the sake of simplicity, α G is supposed to be homogeneous and constant.

Pressure equation

Cell division generates a pressure in the tissue and thus the cellular movement. Darcy law is imposed to link the velocity to the pressure [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF]. Since the membrane is a thin-layer of very low permeability (Fig. 6), we link the outer and inner media thanks to Kedem-Katchalsky conditions [START_REF] Kleinhans | Membrane permeability modeling: Kedem-Katchalsky vs a two-parameter formalism[END_REF][START_REF] Perrussel | Asymptotic expansion of steady-state potential in a high contrast medium with a thin resistive layer[END_REF], which are an asymptotic version of the conditions used in [START_REF] Ribba | A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents[END_REF]:

Ω0 Ω1 Ω1   +  -  -  +
-∇ • (k∇Π) = α P P in R + × (Ω 0 ∪ Ω 1 ), ( 6a
)
k∂ n Π = 0, ( 6b
)
κ m Π = k∂ n Π| Γ + , ( 6c
)
Π| ∂Ω = 0, (6d)
where

• Γ denotes the membrane,

• brackets • denote the jump from Ω 0 to Ω 1 across the membrane Γ,

RR n°8502

• n is the normal to Γ,

• κ m represents the surface membrane permeability, homogeneous to the membrane permeability divided by its thickness.

The medium permeability is denoted by k and as cellular species have different permeabilities, we set on Ω \ Γ:

k = k L L + k S S + k N N + k P P.
Actually we shall call permeability of the medium, by abuse of language, the ratio between the permeability and the viscosity, where the viscosity is assumed to be constant.

The homogeneous condition imposed on the boundary is a kind of transparent boundary condition: it reflects the hypothesis that tumor growth is negligible compared to the volume of tissues that make up the breast. Tumor pushes healthy tissue outside the computational domain without feedback on the internal pressure.

Role of MMPs

In our model, P cells are supposed to produce MMPs which are locally diffused. Then, they degrade the ECM of the biological tissues: membrane, stroma, tumor (Fig. 7).

Cell porosities kS, kP

MMPs

Action on membrane permeability

The pressure accumulated in the tumor that was previously confined to the duct, tends to limit and block the proliferation. But as soon as the membrane is damaged by the action of MMPs, the pressure is released, the cell divisions restart and tumor cells can cross the membrane to invade healthy tissue. MMPs make the basement membrane porous and can break it, by degrading the extracellular matrix, that is to say by breaking cellular adhesions. The surface membrane permeability as a function of MMP concentration is therefore essential in the model. Actually, the medium permeability reflects the degree of cell adhesion. While the membrane is intact, the membrane permeability is denoted by κ 0 which is very small. Then, relation (6c) shows that the pressure problem tends to a Neumann problem: cells cannot cross the membrane. Conversely, Inria a high permeability means that the junctions are fewer and the tissue is more permeable: relation (6c) then expresses the fact that the pressure flux is non-zero across the membrane, as well as the velocity (according to Darcy law). Thus, the action of MMPs, whose concentration is denoted by M , can be modeled by modulating the membrane permeability: tumor remains confined inside the duct if permeability tends to zero and invades if it is larger. We therefore set:

κ m (t, x) = κ 0 + (κ max -κ 0 ) sup s∈[0,t] 1 + tanh [Λ M (M (s, x) -M th )] 2 , ( 7
)
where

• κ 0 denotes the very low membrane permeability when it is intact,

• κ max denotes the membrane permeability when it is completely degraded,

• Λ M is strictly positive, homogeneous to the inverse of a concentration, and gives the slope of the curve,

• M th is the enzymatic threshold from which the MMPs start to degrade the membrane.

We neglect the regeneration of the basement membrane. In the expression [START_REF]Externally applied forces can phenotypically revert malignant breast epithelial structures[END_REF], the supremum in time therefore reflects the irreversibility of the degradation phenomenon.

Action on other tissues permeability

It is obvious that if MMPs act on the membrane extracellular matrix, they have a similar effect on other surrounding tissues. By modeling the ECM degradation of tumor tissue P and tissue S (of permeability k P and k S respectively), we simulate a greater mobility of cancer cells and a more rapid increase in the stroma, consistent with biological processes. k max denotes the maximum permeability of tissues P and S, and a couple (K X , M thX ) is assigned for each tissue X (X stands for P or S) whose permeability is denoted by k X :

k X (t, x) = k 0 X + (k max -k 0 X ) sup s∈[0,t] 1 + tanh [K X (M (s, x) -M thX )] 2 .
where Λ X is homogeneous to the inverse of a concentration. Finally, MMPs have no influence either on the necrotic tissue N or on lumen L. The first one is dead tissue and is subject to calcification that we consider as immediate and relatively low permeability. The lumen is a liquid whose permeability is assumed to be close to that of water, higher than k max .

Basic model of production and diffusion of MMPs

The key issue is the production of MMPs. This question has been extensively studied in biological literature but remains still unclear because of its complexity [START_REF] Duffy | Metalloproteinases: role in breast carcinogenesis, invasion and metastasis[END_REF][START_REF] Overall | Strategies for MMP inhibition in cancer: innovations for the post-trial era[END_REF][START_REF] Köhrmann | Expression of matrix metalloproteinases (MMPs) in primary human breast cancer and breast cancer cell lines: new findings and review of the literature[END_REF]. As a first step, we assume that MMPs are just produced by the proliferating cells P at the rate α M and from the carcinoma initiation. Then, they diffuse through the tissues and are partially inhibited at rate µ by another type of protein: the TIMPs (Tissue Inhibitors of MetalloProteinase). To simplify, we consider that the rate of production α M and the rate of degradation µ are constant. A reaction-diffusion equation is therefore used:

∂ t M -∇ • (D M ∇M ) = α M P đ 0 -µM, ( 8a
) M | t=0 = 0, ( 8b
) M | ∂Ω = 0. ( 8c
)
We assume that MMPs diffuse locally and do not reach the boundary, hence the homogeneous boundary conditions. The evolution equation reflects the fact that although the biological phenomenon of diffusion is faster than the advection on a small spatial scale (a few nanometers), it is much slower at larger scales (a few millimeters) [START_REF] Truskey | Transport Phenomena in Biological Systems[END_REF]. To simulate a local action, we assume that the diffusion coefficient D M has a higher value in the region where the molecule is produced. We therefore define the MMPs diffusion coefficient D M by

D M = D prod M P + D dif f M (1 -P ), 0 < D dif f M < D prod M < ∞.

Nutrients and angiogenesis

Cells need nutrients, in particular for mitosis. As a result, cancer cells require more nutrients than normal cells. Nutrients, whose concentration is denoted by Θ, are provided by the membrane vascularization and by the angiogenic process, assumed to occur during the invasive stage (Fig. 8).

S

Nutrients Θ

Membrane

Avascular stage

The vascularization of the membrane is the only source of nutrients at the avascular stage of tumor growth. We assume that this is an average vascularization on Γ, providing at each point of Γ a constant nutrient concentration equal to Θ max at the constant rate α θ . It is further assumed that the capillaries are damaged at a constant rate β, when tumor cells P invade and crush them.

Inria

We thus define this source of nutrients by

S memb = (α θ -β P ) Θ max 1 Γ . ( 9
)
where 1 denotes the indicator function.

Angiogenesis

The second source of nutrients is angiogenesis, at the stage of vascular growth. We use a model of angiogenesis similar to the one of F. Billy [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]. When the tumor infiltrates outside the duct, it needs more and more nutrients as the number of proliferating cells grows. Then, the hypoxic cells react by producing angiogenic factors, of concentration denoted by ξ, which are diffused:

-∇ • (D ξ ∇ξ) = α ξ f ξ P đ 0 , on Ω, ( 10a
)
ξ| ∂Ω = 0. ( 10b
)
The production rate of angiogenic factors, denoted by α ξ , is weighted by the function f ξ that contains the hypoxic condition:

f ξ (Θ) = 1 -tanh [Λ ξ (Θ -Θ H )(Θ -Θ N )] 2 (11)
where the parameter Λ ξ is homogeneous to the square of a concentration. The signal ξ implies the formation of new blood vessels, composed of endothelial cells, according to the well-known biological process of chemotaxis. Endothelial cells are a part of S cells. Their proportion relative to the cells S is denoted by ν. They are provided at the rate α ν by the initial blood network R 0 , which is assumed to be composed only of endothelial cells. They are then transported in the direction of the gradient of the angiogenic factors, leading to the formation of a new vasculature R. ν is thus driven by the following equations:

∂ t ν + ∇ • (χ ν ∇ξ) = α ν 1 R0 + (α ν -β P ) ν (1 -ν) 1 R , on Ω, ( 12a
)
ν| ∂Ω = 0 (12b) ν| R0 = 1 (12c) ν| t=0 = 0, on Ω R 0 . ( 12d
)
where

R| t=0 = Ø, R = {(t, x) : ν(t, x) ≥ ν min } .
The new vasculature R also provides endothelial cells at the rate α ν , and is subject to degradation by P cells at the rate β. The term (1ν) is a term of saturation to prevent overcrowding. The chemotaxis coefficient χ is built so that the endothelial cells are easily transported in lumen and healthy tissue and do not progress in the tumor tissue and through the basement membrane:

χ = χ L L + χ S S.
Finally, considering that the new vascularization diffuses a concentration of nutrients denoted by Θ max ν at the constant rate α θν , the second source of nutrients is given by

S angio = α θν Θ max ν 1 R0∪R . ( 13
)

Diffusion

The diffusion of nutrients from the two sources S memb and S angio , is assumed to be instantaneous and is described by a static diffusion equation:

-∇ • (D θ ∇Θ) = S angio + S memb -φP Θ, on Ω, ( 14a
) Θ| ∂Ω = 0. ( 14b
)
P cells are supposed to consume nutrients at rate φ = λα P , in proportion to their rate of division. D θ denotes the diffusion coefficient.

Correction

Inside the duct, at the avascular stage, the diffusion radius should be of the order of a few hundred microns for a duct diameter of 1 millimeter. It is therefore decided to truncate artificially the diffusion of nutrients in the center of the duct, otherwise the equation (14a) gives a uniform distribution of nutrients. Truncation avoids introducing a nonlinear diffusion coefficient, which would complicate the numerical resolution. Thereby, we bring an explicit correction to Θ by removing nutrients in the center of the duct according to an exponential law:

Θ corr = Θ exp(-c δ-l R ), if δ ≥ l, Θ, otherwise, (15)
where c is a positive constant, R the duct radius, l the diffusion radius from the membrane towards the center of the duct and δ the distance to the membrane inside the duct (estimated using a level set function).

3 Numerical results

Implementation

The following 2D-simulations of DCIS and IDC are implemented in C++ from a finite-volume discretization on cartesian grid. For the sake of simplicity, the duct is assumed to be cylindrical (Fig. 9). We consider the domain from two points of view: in longitudinal and cross sections. In this paper, results are mostly presented in cross section, in order to be compared with some medical biopsies. Simulations are carried out on a 100 × 100 mesh. Significant parameters used in the simulations are given in Table 1. They are empirically determined to reproduce biological observations. Note that the rate of MMP production α M is vanishing for DCIS simulations and equal to 1000 for IDC simulations.

Inria

Total cell density Pressure threshold

đ 0 1 Π max 5 Rates Diffusion coefficients α 8 D prod M 1 α M 0 -1000 D dif f M 0.1 µ 2 D θ 0.1 α θ 1 D ξ 2 α θν 1 Chemotaxis coefficients β 0.2 χ L 1 α ξ 10 χ S 1 α ν 1 Lengths Concentration thresholds l 0.05 Θ H 3 R 20 Θ N 2.

Initial setup

Cancer initiation

Some images raise questions about the starting point of a tumor. Fig. 10 seems to suggest that the tumor does not necessarily start from only one defective cell but from many. As the bibliographic sources are unclear on the topic, we assume that a tumor may be derived from one or more cells. In the case of a single initial cancer cell, we set P 0 = 1 in a small area supported by the membrane at the bottom of the duct. In the case with many initial cancer cells, we consider that they have appeared along the epithelial tissue, and they are initially randomly distributed in space. We assume that along the duct membrane, on the inner part, one fifth of the cell density is defective. For each simulation, we will specify if it is a 'single-cell simulation' or a 'multi-cell simulation'.

Healthy tissue and nutrients

Initially, for all runs, we set

N 0 = 0, in Ω, S 0 = 1, in Ω 1 , L 0 = 1 -P 0 , in Ω 0 .
The pressure is assumed to be zero and nutrient initial conditions are set such that

f 0 P = α G , f 0 N = 0.
Moreover, R 0 , the initial blood network for angiogenesis, is localized far away from the duct.

Preliminary results: nutrient influence for a DCIS

We first consider a carcinoma in situ, without MMPs, by setting α M = 0. As expected, if nutrient conditions are set to be sufficient (Θ max = 1000 for example) the tumor fills the duct (solid carcinoma). Fig. 11, obtained by modulating the source Θ max (single-cell simulation), shows the formation of a quiescent core, no longer proliferating, when the nutrient concentration is below the hypoxia threshold. In Fig. 11b, the nutrient concentration at the center of the tumor is below the threshold of necrosis, resulting in cell death. There is a difference in growth rate due to the wider area of quiescence and necrosis in Fig. 11b. Necrosis can occur in large ducts in the center of which the nutrient concentration is low. But it can also appear in ducts of all sizes, depending on the growth rate of the tumor. The quicker the cells proliferate, the more nutrients they need and the more they deprive their counterparts at the center of the tumor, resulting in necrosis formation. Then, the size of the necrotic core gives an indication of the Inria aggressiveness of the cancer. This kind of tumor is called comedo-carcinoma (Fig. 12, multi-cell simulation, Θ max = 50, t=10), a type of breast cancer often high-graded, associated with a poor diagnosis.

The limitation of cross-sectional 2D simulations

In the cross-sectional 2D-simulations, L cells cannot be evacuated by the tumor in the transverse direction: they are enclosed by the membrane and they cannot move. Therefore, the saturation condition (1) is used to compute rather than the advection equation (2d):

L = 1 -(P + S + N ).
Yet the problem is not resolved: the pressure cannot normally be removed from the duct. Actually, the homogeneous Dirichlet condition (relation (6d)) plays the role of a transparent condition on the boundary ∂Ω, which makes it possible for the pressure to be evacuated from the computation domain. In contrast to the case in longitudinal section, the interior media in cross section is not connected to the boundary. Consequently, the pressure is not evacuated and builds up in the duct, leading to the proliferation rate decrease. The result is a difference in growth rate between the transverse and longitudinal sections (Fig. 13, single-cell simulation). Moreover, in the case of a tumor derived from several cancer cells (Fig. 14), the excessive pressure causes an artificial competition between them. The tumor foci, in the area where they are more numerous, create a stronger pressure gradient. Consequently, the velocity field is directed towards the area where the foci are the fewer. These seem to shrink and disappear (Fig. 14b). In fact, tumor cells are dispersed within the lumen, which can be checked by changing the scale of observation. We overcome this difficulty by artificially decreasing the pressure accumulated in the duct according to an exponential law:

Π corr = Π exp(- δ R 2 ),
where δ is the distance to the membrane, inside the duct, and R is the duct radius. In the remainder of the paper, the corrected pressure Π corr will be denoted by Π. Inria

Simulation of tumor invasion

We assume that P cells produce MMPs, by setting α M = 1000. The results in Fig. 15 (single-cell simulation) show the transition from the in situ stage to invasion of a carcinoma, starting from a single initial cancer cell. It can be seen in Fig. 15b and Fig. 15c that infiltration logically starts at the bottom, where the original tumor was. With the same parameters, the simulation in longitudinal section (Fig. 15d) shows the length of duct affected by membrane degradation and tumor invasion. Taking angiogenesis into account, we obtain results from the complete model with the multi-cell simulation given by Fig. 16. In addition to growth and tumor invasion, Fig. 16 shows the formation of a central necrosis (due to selected nutrient conditions) and the outbreak of the angiogenesis that will allow the tumor to continue growing in healthy tissue.

In previous results, carcinoma completely degrades the membrane but we can also simulate a micro-infiltrating carcinoma. According to [START_REF] Man | The signifiance of focal myoepithelial cell layer disruptions in human breast tumor invasion: a paradigm shift from the "protease-centered" hypothesis[END_REF], membrane damage is multifactorial and is not limited only to the enzymatic action: we should also consider that some parts of the membrane are weaker than others. We decided to simulate a more sensitive area by modulating the slope coefficient Λ M in expression ( 7) of the membrane permeability. By assigning a lower value on a given area of the membrane (Λ M = 0.2 on this area, 2 elsewhere), MMPs act faster and degradation only concerns that area. Simulation, as shown in Fig. 17 (multi-cell simulation), then yields a micro-invasive tumor: the membrane degradation and invasion are highly localized.

Different ways for describing heterogeneous tumors

Although the model is capable of reproducing the behavior of a carcinoma that infiltrates through the membrane, it remains relatively unsatisfactory when compared with biomedical results. In particular, if the carcinoma has filled out the duct, the entire membrane is rapidly degraded, due to a homogeneous distribution of tumor cells and therefore of MMPs they produce. Thus, despite random initial data, the tumor growth is still too uniform (Fig. 18a vs. Fig. 18b).

Obviously, a more irregular duct shape as well as stroma anisotropy could partially explain the tumor asymmetry. But we decided to overcome the lack of heterogeneity by first considering the effect of a genetic mutation within the tumor, and then the influence of the microenvironment on the production of MMPs and TNFs. In all cases, we obtain a greater inhomogeneity and more realistic results.

Influence of genetic mutations of tumor cells

It is well known that genetic factors are critical in tumor initiation and progression. As cancer is the result of genetic mutations that can be seen as random in space, a tumor is usually inhomogeneous in terms of growth and substances which it produces. We assumed in previous section that each cell produced MMPs from the initial time. We now consider that a random mutation is occuring in a few selected cells that makes them produce MMPs (Fig. 19). A genetic mutation means that a cell acquires a specific property which is cloned during cell division. It requires the definition of new species for each mutation. We then introduce two species of proliferative cells, P 1 and P 2 . MMPs are only produced by P 2 . They satisfy the same advection equation (2a) than P :

∂ t P 1 + ∇ • (vP 1 ) = (α P -α N )P 1 , ( 16a
)
∂ t P 2 + ∇ • (vP 2 ) = (α P -α N )P 2 , ( 16b
)
P = P 1 + P 2 , ( 16c
)
and relation (8a) becomes

∂ t M -∇ • (D M ∇M ) = α M P 2 đ 0 -µM. ( 16d
)
We assume that the mutation has occured before the initial time of simulations. We perform a multi-cell simulation, in which P 0 1 and P 0 2 are randomly distributed such that: P 0 = P 0 1 + P 0 2 , P 0 1 = 0.8 P 0 , P 0 2 = 0.2 P 0 As the production of MMP is quantitatively changed, we choose specific values for the MMP thresholds given by Table 2. As a result, MMP concentration and membrane degradation are less homogeneous (Fig. 20) and we obtain some less symmetrical and more realistic tumors. We could consider that genetic mutations occur later and obtain the same kind of results. Obviously, it is possible to introduce some other genetic mutations. For example, some onco-genetic mutations imply overexpression of cell receptors like HER-2 that, when activated, increase the mitosis rate. This phenomenon has already been modeled [START_REF] Thanoon | Computational framework for local breast cancer treatment[END_REF] and leads to a rate α G becoming variable in space and time, in relation (5a). Then, the proliferation rate would be inhomogeneous and could also lead to various shapes of tumor. What can be retained is that different genetic mutations can be simulated and give more realistic results, with inhomogeneities in the growth and shape of the tumor.

M th 10 M thS 8 M thP 5

Influence of the microenvironment on production of MMPs

However, genetic factors are not the only ones involved in cancer progression. More and more biological studies underline that the microenvironnement is also essential. We focus on the impact modeled [START_REF] Man | The signifiance of focal myoepithelial cell layer disruptions in human breast tumor invasion: a paradigm shift from the "protease-centered" hypothesis[END_REF]. Some studies mention the role played by cells of epithelial tissue (myoepithelial cells) and stromal cells (fibroblasts, macrophages of the immune system, etc.) in MMP production [START_REF] Brummer | Matrix-metallo-proteinases 1, 2 and 3 and their tissue inhibitors 1 and 2 in benign and malignant breast lesions: an in situ hybridization study[END_REF][START_REF] Cichon | Microenvironmental influences that drive progression from benign breast disease to invasive breast cancer[END_REF][START_REF] Hu | Microenvironmental regulation of cancer development[END_REF][START_REF] Rha | Sequential production and activation of matrix-metalloproteinase-9 (MMP-9) with breast cancer progression[END_REF]. In response to a signal from the tumor, stromal cells migrate towards the duct and produce MMPs as inactive proenzymes. These are then activated by forming complexes with MMPs produced by tumor cells (Fig. 22).

To model these processes simply, some stromal cells S are randomly labeled as fibroblasts denoted by S 2 , which avoids introducing a new population that would migrate to the tumor: the important point is that these cells can produce MMPs. Other stroma cells are denoted by S 1 . Then, proliferative cells diffuse a chemical signal, of concentration Z, to activate the proenzymes (Fig. 23). The local action of the signal is simulated with the same diffusion coefficient than for MMPs. Actually, we assume that the MMP production starts only if the signal concentration reach a certain threshold Z s . We therefore set:

∂ t S 1 + ∇ • (vS 1 ) = 0, ( 17a
)
∂ t S 2 + ∇ • (vS 2 ) = 0, (17b)
S = S 1 + S 2 . ( 17c
)
Relationship (8a) is replaced by

∂ t M -∇ • (D M ∇M ) = γS 2 đ 0 -µM, ( 17d
)
∂ t Z -∇ • (D M ∇Z) = α z P đ 0 -µZ, ( 17e
)
Z| t=0 = 0, (17f)
Inria where

γ = α M 1 + tanh [Λ z (Z -Z th )] 2 ,
and Λ z is homogeneous to the inverse of a concentration. We give the value of the new parameters in Table 3 and, as previously, we choose specific values for MMP parameters. Regarding the initial conditions, S 2 cells are useful only in a neighborhood of the duct and they are randomly distributed such that, in this neighborhood:

Λ z 30 Z th 1 α z 30 M th 6 M thS 4 M thP 2
S 0 = S 0 1 + S 0 2 , S 0 1 = 0.9 S 0 S 0 2 = 0.1 S 0 .
The results in Fig. 24 (multi-cell simulation) show a more local breaking point in the membrane than previously (Fig. 24a). This more localized gap is due to the random distribution of stromal cells that produce MMPs. Before the membrane is further deteriorated, the tumor starts to come out and invade healthy tissue, brimming over the duct. The resulting tumor is asymmetrical and an analogy is observed between tumor shapes in Fig. 24c and Fig. 24d. However, it should be noted that Fig. 24d does not exactly represent an invasive ductal carcinoma, but its metastases that created a secondary focus in a lymph node, which in turn became invasive.

TNF action

In the same vein, we could take another microenvironmental influence into acount: the TNF action on necrosis core (Fig. 25). The lack of nutrient is not the only responsible for the development of a necrotic core: the tumor necrosis factors (TNFs) seem to be involved with conflicting effects [START_REF] Waterston | TNF and cancer: good or bad ?[END_REF][START_REF] Hamaï | Le facteur de nécrose tumorale: de la biologie à la thérapie oncologique[END_REF]. It is a chemical secreted by immune cells that infiltrated the tumor and that cause cell death. TNFs could also play a role in promoting cancer: initiation, proliferation of malignant cells, angiogenesis and metastases production. Motivated by the idea of modeling more realistic necroses, larger and with an irregular shape (Fig. 26c), we decide to take the TNF actions into account. TNFs are supposed to be produced and diffused by infiltrating macrophages, denoted P ′ , which are chosen among P cells, randomly in space and time,

-∇ • (D τ ∇τ ) = α τ P ′ đ 0 , ( 18
)
where τ denotes the TNF concentration with τ 0 = 0. The action of TNF on the rate of necrosis α N is introduced in the weight function f N and is assumed to be described by equation ( 19) instead of equation (5b):

f N = 1 - 1 4 [1 + tanh (Λ N (Θ -Θ N ))] • [1 -tanh (Λ τ (τ -τ th ))] , ( 19
)
where Λ τ is homogeneous to the inverse of a concentration. For a given concentration of nutrients, the rate of necrosis increases if the TNF concentration exceeds a certain threshold τ th .

The simulations reveal a more extensive and irregular necrosis (Fig. 26b) than if only hypoxia is involved (Fig. 26a). The wider necrosis is due to the increase in necrosis rate, and its irregularity to the random distribution of macrophages. Introducing the action of TNFs in the calculation of the weight function f P of the proliferation rate (equation (5a)), in the same way as in equation ( 19), we could also simulate a faster growth under the influence of TNFs, as they are supposed to promote tumor cells proliferation. Once again, despite a very partial modeling of biological processes involved, the results seem more realistic (Fig. 26b and Fig. 26c) when the microenvironment is taken into account.

Conclusion

The aim of the study was to model the transition of a breast cancer between its intraductal phase and the invasive stage, by accounting for the action of MMP enzymes. A continuous model of tumor growth was developed from advection equations. The use of Darcy law has made it possible to describe the movement of cell species, depending on the pressure generated by the division of tumor cells. The very low thickness of the membrane was taken into account by using a Kedem-Katchalsky model: the membrane is considered as an interface through which there is a pressure jump. The issue of cross section in 2D-simulations was solved specifically by the addition of an explicit correction to the pressure inside the duct. At this stage, the simulations yielded the growth of a tumor inside a galactophoric duct. Then, the role of MMPs has been included in order to model the membrane degradation. MMPs are the cornerstone in the transition from an intraductal carcinoma to the invasion stage. Considering that the permeability of the medium simulates correctly the degree of cell adhesion within tissues, the process of membrane Inria degradation was described by the permeability, which is a function of enzyme concentration. The simulations highlighted the degradation of the membrane, the spread of malignant cells outside the duct and the tumor invasion. However, despite adding environmental conditions such as nutrient diffusion and angiogenesis, the model remained rather inaccurate by assuming that MMPs was produced by tumor cells in a homogeneous way: the obtained tumors were quite regular and symmetrical, much different from those observed in the few available biomedical results. Then, the model was enriched in the purpose of modeling more complex processes involved in tumor growth: taking some genetic factors into account first, and secondly some influences of the microenvironment, has led us to more satisfactory results, one of which is given by Fig. 24. Maybe another specific model could lead to the same kind of result. But it should be noted that to obtain Fig. 24, we need most of the new developments of the model, which therefore appear to be important in order to describe reality: the mechanical specificity related to the duct, the pressure influence on proliferation and its correction, the forced diffusion of nutrients within the duct to simulate the avascular stage, the multi-cell simulation and the microenvironment role in MMP production. In particular, this result seems to support the idea that microenvironmental factors play a crucial role in tumor progression, in agreement with many recent studies in oncology [START_REF] Bissell | Why don't we get more cancer ? A proposed role of the microenvironment in restraining cancer progression[END_REF][START_REF] Cichon | Microenvironmental influences that drive progression from benign breast disease to invasive breast cancer[END_REF][START_REF] Hu | Regulation of in situ to invasive breast carcinoma transition[END_REF].

Despite the satisfactory behavior of the model, we can highlight some points for further study. The use of Darcy law, at the heart of the model, raises issues that indicate some limits of the model. First, when the tumor grows and fills the duct, it may tend to deform the duct membrane, implying a feedback on the rate of tumor growth due to the membrane elasticity. To take the membrane elasticity into account, we could use a Stokes law [START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF][START_REF] Franks | Modelling the early growth of ductal carcinoma in situ of the breast[END_REF] rather than Darcy law, which would significantly complexify the model. Secondly, carcinoma in situ may have very different architectures from one cancer to another. The model provides only two of the main types of breast carcinomas: solid and comedo carcinomas. Other types of breast intraductal carcinoma mentioned in the medical literature are:

• the micropapillary carcinoma develops along the epithelium, with tumor cells projections towards the center of the duct, without filling (Fig. 27a)

• the papillary carcinoma develops inside the duct with irregular shapes (Fig. 27b)

• the cribriform carcinoma is caracterized by distinctive holes between the cancer cells (Fig. 27c). These other types of carcinomas are more difficult to model: in each case (Fig. 27), projections and spans of tumor cells appear to indicate that some of adhesion and polarity properties have been kept by the cells. In our model, the permeability coefficient, introduced with Darcy law, does not distinguish between the different types of adherent junctions. Above all, it simulates neither the polarization of proliferating cells nor the orientation of cell divisions. Thus, the continuous model finds its limits when it comes to simulate phenomena that are played at the cellular level. A discrete model, where each cell is modeled by an agent, like in [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF][START_REF] Lesart | On the importance of the submicrovascular network in a computational model of tumour growth[END_REF][START_REF] Macklin | Patient-calibrated agentbased modelling of ductal carcinoma in situ (DCIS): from microscopic measurements to macroscopic predictions of clinical progression[END_REF], seems more appropriate. Note that Norton et al. model those different types of carcinomas at the stage in situ, in [START_REF] Norton | A 2D mechanistic model of breast ductal carcinoma in situ (DCIS) morphology and progression[END_REF]. Coupling the continuous model with a discrete model could be considered as in [START_REF] Kim | The role of the microenvironment in tumor growth and invasion[END_REF] to simulate an invasive phase of these types of cancer. There still remains numerous unanswered questions regarding the biological phenomena that are involved: the impact of the pressure, or mechanical stress, on tumor growth seems to be ambiguous [START_REF] Yang | To revert breast cancer cells, give them the squeeze[END_REF][START_REF] Tse | Mechanical compression drives cancer toward invasive phenotype[END_REF]. It could even have an impact on genetic mutations and production of MMPs. The production of MMPs seems even more complex than what we simulated, involving different kinds of MMPs. The structure and activity of the membrane, here reduced to an interface, could also be a determining factor in the invasive process [START_REF] Hu | Regulation of in situ to invasive breast carcinoma transition[END_REF]. Finally, the migration of specific cells of the microenvironment (such as immune cells, macrophages, fibroblasts), randomly simulated in our study, should be further studied in the hope of bringing predictability to the model.

It is worth noting that according to our model, either oncogenic or microenvironmental factors make possible to obtain heterogeneous tumors, in accordance with biopsies. Our paper can be seen as a first step in the study of the influence of these factors, and in their modeling in order to provide realistic tumor shapes.

Figure 1 :

 1 Figure 1: Simplified anatomy of the breast and different stages of the evolution of a duct carcinoma.

  (a) IDC simulation (b) Biopsy [9]

Figure 2 :

 2 Figure 2: Main result : qualitative comparison of our simulation with biopsy.

Figure 3 :

 3 Figure 3: Schematic diagram of the model structure.

Figure 4 :

 4 Figure 4: Simple geometry: Ω denotes the whole domain, Ω0 the duct, Ω1 the stroma, Γ the duct membrane with zero thickness.

Figure 5 :

 5 Figure 5: Schematic diagram of submodel 1.

Figure 6 :

 6 Figure 6: Γ is considered as an interface.

Figure 7 :

 7 Figure 7: Schematic diagram of submodel 2.

Figure 8 :

 8 Figure 8: Schematic diagram of submodel 3.

Figure 9 :

 9 Figure 9: Domain Ω for simulations.

Figure 10 :

 10 Figure 10: Micropapillary tumor in situ. Emergence of several cancerous foci. Biopsy from [20].

( a )Figure 11 :

 a11 Figure 11: Nutrient influence: quiescent cells are in light blue and necrosis in dark blue (t=11).

( a )

 a Necrosis by simulating hypoxia (b) Tumor with necrosis[START_REF] Sakr | Carcinomes canalaires in situ du sein: rôle potentiel de la biologie moléculaire[END_REF]

Figure 12 :

 12 Figure 12: Tumors with a necrotic center.

Figure 13 :

 13 Figure13: Difference of growth between sections: with the same parameters and at the same time (t = 6), the tumor fills the duct in longitudinal section, but not in cross section where an area consisting only of lumen is observed (white area).

Figure 14 :

 14 Figure 14: Some of the tumor foci retract (t = 7 and t = 10).

= 9 Figure 15 :

 915 Figure 15: Evolution of an invasive carcinoma.

Figure 16 :

 16 Figure 16: Invasive carcinoma with angiogenesis (α = 5 and Θ max = 200). Some blood vessels resulting from angiogenesis are indicated by arrowheads.

Figure 17 :

 17 Figure 17: Comparison of microinvasion between biopsy and simulation with the parameter αM = 600 at t = 8.

Figure 18 :

 18 Figure 18: Invasive tumors (biopsy from [24]). The simulation gives a too symmetrical tumor.

Figure 19 :

 19 Figure 19: Genetic influence on MMP production.

Figure 20 :Figure 21 :

 2021 Figure 20: Invasive carcinoma with two populations of proliferating cells. P 1 cells are in blue. P 2 cells, in red, produce MMPs (α = 5).

Figure 22 :Figure 23 :

 2223 Figure 22: Role of stroma in membrane degradation from [27].

Figure 24 :

 24 Figure 24: Improved model simulating the evolution of an invasive carcinoma (α = 5).

Figure 25 :

 25 Figure 25: Schematic diagram of TNF production and action on necrosis and proliferation rates.

(

  a) ατ = 0 (b) ατ = 1000(c) Comedo-carcinoma from[START_REF]Breast Pathology on the Web[END_REF]

Figure 26 :

 26 Figure 26: Simulation of the TNF effect on tumor necrosis (t = 11, Λτ = 8, Dτ = 2, τ th = 0.5).

Figure 27 :

 27 Figure 27: Other types of intraductal carcinoma. Taken from [21].

Table 1 :

 1 5 Simulation parameters. Values are empirically determined.

	Weight function coefficients

Table 2 :

 2 Changes in the MMP thresholds.

Table 3 :

 3 New parameters and changes in the MMP thresholds for a model that takes the microenvironment into account.

RR n°8502

Inria

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

Acknowledgements

This study has been carried out 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 Bordeaux -CPU (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/).

Inria