Sergey Astanin
email: astanin@calvino.polito.it

Luigi Preziosi
email: luigi.preziosi@polito.it

Mathematical modelling of the warburg effect in tumour cords

Keywords: tumour growth, tumour metabolism, Warburg effect, mathematical model, ghost fluid method, population dynamics

The model proposed here links together two approaches to describe tumours: a continuous medium to describe the movement and the mechanical properties of the tissue, and a population dynamics approach to represent internal genetic inhomogeneity and instability of the tumour. In this way one can build models which cover several stages of tumour progression. In this paper we focus on describing transition from aerobic to purely glycolytic metabolism (the Warburg effect) in tumour cords. From the mathematical point of view this model leads to a free boundary problem where domains in contact are characterized by different sets of equations. Accurate stitching of the solution was possible with a modified ghost fluid method. Growth and death of the cells and uptake of the nutrients are related through ATP production and energy costs of the cellular processes. In the framework of the bi-population model this allowed to keep the number of model parameters relatively small.

Metabolic processes in normal tissues require oxygen. Specifically, 6 molecules of oxygen are consumed per oxidated molecule of glucose with the yield of approximately 32 molecules of ATP [START_REF] Nelson | A fully conservative ghost fluid method stiff detonation waves[END_REF] (36 according to [START_REF] Smallbone | Metabolic changes during carcinogenesis: potential impact on invasiveness[END_REF], 29.85 according to [START_REF] Rich | The molecular machinery of keilin's respiratory chain[END_REF]). In many cancers intensive proliferation exceeds available oxygen supply, which leads to hypoxia. In such hypoxic conditions cells may rely only on glycolysis, the first step of glucose oxydation, to cover their energy needs. This process gives a smaller amount of ATP, 2 molecules of ATP per molecule of glucose, but it is possible in hypoxic conditions. In fact, most tumours are known to rely on glycolytic metabolism even in non-hypoxic conditions. This effect is known as Warburg effect or aerobic glycolysis [START_REF] Warburg | On the origin of cancer cells[END_REF][START_REF] Kim | Cancer's molecular sweet tooth and the Warburg effect[END_REF] and is one of the hallmarks of cancer [START_REF] Hannahan | The hallmarks of cancer[END_REF]. Glycolytic catabolism has the important side effect of tissue acidification. Lower pH is toxic to most normal cells while altered tumour cells are likely to be resistant to it and achieve another invasion advantage [START_REF] Gatenby | Acid-meditated tumour invasion: a multidisciplinary study[END_REF]. In the same time, there are several therapeutic strategies which allow for targeting tumours with glycolytic metabolism [START_REF] Kim | Cancer's molecular sweet tooth and the Warburg effect[END_REF][START_REF] Mathupala | Lactate and malignant tumors: a therapeutic target at the end stage of glycolysis[END_REF].

We want to describe phenomenologically the transition of tumours from normal to glycolytic metabolism in the framework of spatio-temporal model of tumour growth. Well aware that there are several mechanisms which contribute to the Warburg effect, we assumed that this switch in metabolism happens as an all or nothing event, after which cells rely only on glycolytic metabolism even with adequate oxygen levels. The switch is assumed to happen in hypoxic conditions. This problem received a lot of attention from the mathematical modelling community in the recent years. Such works as [START_REF] Gatenby | Acid-meditated tumour invasion: a multidisciplinary study[END_REF][START_REF] Gatenby | Cellular adaptations to hypoxia and acidosis during somatic evolution of breast cancer[END_REF][START_REF] Smallbone | Mathematical modelling of tumour acidity[END_REF] have pointed to the possibility of the acid-mediated tumour invasion and studied the role of the cellular adaptivity in tumour invasion. Their discrete and continuous models and experiments demonstrate that emergence of the acidresistant, glycolytic cell line is an important step towards cancer invasiveness. Another recent work is [START_REF] Gerlee | A hybrid cellular automaton model of clonal evolution in cancer: The emergence of the glycolytic phenotype[END_REF], where the cellular automaton model demonstrates that an advantage of the glycolytic phenotype may be also conditioned by the density of the matrix. In [START_REF] Venkatasubramanian | Incorporating energy metabolism into a growth model of multicellular tumor spheroids[END_REF] ability of the cells to rely on the glycolytic metabolism was incorporated in the continuous model of tumour spheroid. In [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF] the rigorous analytical study of a similar model was undertaken, with a proof that the necrotic core is formed for any sufficiently large spheroid. Both in [START_REF] Venkatasubramanian | Incorporating energy metabolism into a growth model of multicellular tumor spheroids[END_REF] and in [START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF] the switch to glycolytic metabolism is assumed to be reversible, and the cells return to normal aerobic metabolism when the hypoxia finishes. However, the shift towards glycolytic metabolism is often irreversible in tumours, and they may remain glycolytic even in the presence of sufficient oxygen [START_REF] Kim | Cancer's molecular sweet tooth and the Warburg effect[END_REF]. In this case we have to consider at least two distinct cell lines with different (probably very different) metabolic behaviour, because after the switch the cells are not the same any more.

Having to construct a multi-population model certainly complicates the numerical method, but it also brings some advantages of the discrete models in the reign of the continuous A c c e p t e d m a n u s c r i p t models. In fact, possible applications of the proposed numerical approach go far beyond tumour metabolism, and it may be used to describe other genetic transformations in the course of tumour evolution.

We also tried to build the model around the idea used in (Astanin & Tosin, 2007), that the growth in the tissue strictly depends on the amount of the nutrients consumed (or more precisely, the quantity of ATP produced). To put it simply, growth and intensity of the metabolic processes are two sides of the same medal in this model. Insufficient ATP production, in turn, limits growth and leads to necrosis.

Another aim of the paper is to understand whether the metabolic switch gives rise to a spatial stratification in the bipopulation model. In fact, near the vessels oxygen is more abundant and cells have no reason to change metabolism there. This might lead to the conclusion that tumour cords are constituted by an inner core near the vessel with cells using an aerobic metabolism, and an outer layer of cells using glycolysis.

Counter intuitively, in the simulations we observed that the cell line with glycolytic metabolism "wins", and spreads in the entire cord, though in a non-homogeneous way. The intensity of metabolism affects the cord geometry and the moment when tumour starts suffering from hypoxia. This is due to the fact that cells which rely on anaerobic glycolysis move also to the inner layer, spreading in the whole tumour. They may achieve a significant mass even if aquisition of this trait is a very rare event.

This paper has the following structure. In the first section (Bi-population model) we define the equations of the mathematical model, paying a particular attention to the coupling of the multiphase model with the subcellular catabolic mechanisms generating the growth, death, and switch terms in the mass balance equations. The second section (Simulations) contains some 2D simulation results and their discussion. Appendices present then a nondimensional model (Appendix A) and describe a modification to the ghost fluid method which allows to stitch solutions between two different modelling domains (Appendix B).

Bi-population model

In this paper we develop the approach of earlier works (Ambrosi [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF]Astanin & Tosin, 2007;[START_REF] Preziosi | Multiphase modeling of tumor growth and extracellular matrix interaction: mathematical tools and applications[END_REF][START_REF] Astanin | Mathematical model of tumour cord growth along the source of nutrient[END_REF] but, in addition to describing tissue growth and mechanics, we consider the transition of tumour cells from one metabolic behaviour to another, and incorporate effects of glycolytic switch in tumour cells. We carefully consider the implications of altered glucose metabolism on the energy balance in tumour cells and examine in detail the spatial effects.

In the previous works we assume that tumour grows in a host enviroment. The model deals with only two cell populations, tumour cells with volume ratio φ t and host cells with volume ratio φ h initially occupying different domains Ω t (t = 0) and Ω \ Ω t (t = 0).

Two other components of the tissue are considered in addi-tion to cells: extracellular liquid (with volume ratio φ l ) and extracellular matrix (ECM, with volume ratio φ m ). ECM is rigid, uniform and does not remodel:

φ i + φ l = 1 -φ m = φ * ≡ const, i := t, h. (1)
An immediate consequence of this hypothesis is that from the macroscopic point of view the tissue would behave like a rigid porous medium, with cells and extracellular liquid moving inside a rigid scaffold. Any stress acting on the bulk tissue would be sustained by the ECM and cells in the core of the tissue would experience no stress deriving directly from the external actions. It is possible to relax this hypothesis following for instance [START_REF] Ambrosi | Cell adhesion mechanisms and elasto-viscoplastic mechanics of tumours[END_REF].

The interface ∂ Ω th between tumour and environment is a material surface moving with the common velocity of the cells

n • d x t dt = n • v t = n • v h , on ∂ Ω th (t) , ( 2
)
which implies that if the two populations occupy initially different interfacing domains, they will always occupy different interfacing time-dependent domains. However, the two populations mechanically and chemically interact with each other.

In the sense that, for instance, the growing tumour will mechanically compress the host tissue and will compete for resources.

In particular, it will induce hypoxia and subsequent death of the surrounding host tissue. Then at the interface in addition to imposing continuity of velocity and stress (see ( 7) and ( 8) below), continuity of oxygen and glucose concentrations and fluxes holds. Treating for the sake of simplicity the ensemble of cells as elastic fluids, neglecting liquid-cell interactions, we have the following free boundary problem, see [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF] for details:

∂ φ t ∂t + ∇ • (φ t v t ) = Γ t , in Ω t , ( 3
)
v t = -K∇(φ t S(φ t )) , in Ω t , (4) ∂ φ h ∂t + ∇ • (φ h v h ) = Γ h , in Ω \ Ω t , ( 5
)
v h = -K∇(φ h S(φ h )) , in Ω \ Ω t , ( 6
)
v t • n = v h • n , on ∂ Ω th , (7) φ t S(φ t ) = φ h S(φ h ) , on ∂ Ω th , ( 8
)
where K is the "motility" of the cells or rate of tissue relaxation, being inversely proportional to the drag force coefficient between cells and rigid ECM, S(φ i ) is a scalar function which describes elastic forces of the compressed tissue (φ i > φ 0 , φ 0 = const), and intercellular adhesion forces in the stretched tissue (φ i < φ 0 ), Γ t and Γ h are growth/death rates for tumour and host cells respectively.

In (Astanin & Tosin, 2007) we paid special attention to the energy balance in cells, supposing that the amount of oxygen and glucose consumed stays in direct correspondence with the intensity of tumour proliferation. We have shown that tumour growth in this model is limited by the formation of the hypoxic zone and by the availability of oxygen.

In this work we build the model upon the same idea. We assume that if ATP production is insufficient the cells stagnate

A c c e p t e d m a n u s c r i p t

and die (like in the old model), but also with a certain (small) probability may irreversibly switch their metabolic behaviour and survive. This switch corresponds to upregulation of the glycolysis, and is not present in (Astanin & Tosin, 2007). The growth also becomes limited by the availability of glucose as well as oxygen (glucose was assumed to be abundant in the old model).

So we introduce the second population of tumour cells which are mechanically identical to the others, but consume only glucose through glycolysis. We refer to this subpopulation as anaerobic or glycolytic subpopulation, while the initial oxygen-consuming subpopulation is refferred to as aerobic one. Transition from the aerobic to the anaerobic population is supposed to be irreversible. We also assume that cells of the host always rely on aerobic glucose oxidation.

Let φ 1 be the volume fraction of cells with normal aerobic catabolism, and φ 2 be the volume fraction of cells with altered, anaerobic, catabolism. On the right hand side of the mass balance equation related to φ 1 and φ 2 we will distinguish a growth rate Γ + i , i = 1, 2, a death rate Γ - i , i = 1, 2, and a switch rate Γ s , so that

Γ 1 = Γ + 1 -Γ - 1 -Γ s , (9)
and

Γ 2 = Γ + 2 -Γ - 2 + Γ s . ( 10
)
Therefore, for the tumour tissue as a whole

Γ t = Γ 1 + Γ 2 . ( 11
)
Looking more closely at the two catabolic pathways of glucose oxydation (shown in Figure 1), full glucose oxidation requires 6 molecules of oxygen per molecule of glucose and produces an equivalent of approximately 30-32 molecules of ATP [START_REF] Nelson | A fully conservative ghost fluid method stiff detonation waves[END_REF] (36 according to [START_REF] Smallbone | Metabolic changes during carcinogenesis: potential impact on invasiveness[END_REF] or 29.85 according to the evaluation in [START_REF] Rich | The molecular machinery of keilin's respiratory chain[END_REF]). If only glycolysis takes place, the energy yield is much smaller (only 2 molecules of ATP per molecule of glucose) but it gives an evolutionary advantage to the cells capable of this kind of catabolism, because they no longer depend on oxygen supply.

We have to consider at least two substances: c o , the concentration of oxygen, and c g the concentration of glucose. These substances are often referred to as nutrients, which is not precise, because of the special role of oxygen in respiration process. For our purpose it is important that both of these substances are critical resources for tissue growth. We refer to all such chemical substances collectively as nutrients, though they do not necessarily have nutritional value. We consider that both nutrients are spread by diffusion with diffusion coefficient of oxygen D o , and diffusion coefficient of glucose D g . Lactic acid concentration may be considered later to take into account of effects due to the lower pH as done in [START_REF] Gatenby | Acid-meditated tumour invasion: a multidisciplinary study[END_REF].

In the model, we represent two pathways as two concurrent reactions. In aerobic conditions:

Glc + 6 O 2 -→ 6 CO 2 + N ATP,
where N is the number of ATP molecules produced during the complete oxydation of glucose; and in anaerobic conditions: Glc -→ 2 Lactic acid + 2 ATP.

We assume that the rate of ATP production in the former reaction (in the aerobic population) is proportional to the amount of cells and the concentration of nutrients:

q (1) AT P = κ (1) φ φ 1 c o c g , ( 12
)
and the rate of ATP production in anaerobic population is q

(2)

AT P = κ (2) φ φ 2 c g . ( 13
)
We should keep in mind that the reaction rates are not necessarily proportional to the concentrations of oxygen and glucose. This is only a basic assumption, and the structure of the model will be valid also with different rates of ATP production.

The coefficients κ

(1) φ and κ

(2) φ also reflect inhibition of metabolism. Here we only consider contact inhibition and take

κ (1) φ = κ(φ * -φ 1 -φ 2 ) ( 14
)
and κ

(2)

φ = κ 2 (φ * -φ 1 -φ 2 ) (15)
meaning that metabolism is inhibited when the cells occupy all the available volume φ * . If inhibition occurs at cell volume fraction φ 1 + φ 2 = φ = const, this case requires only rescaling of volume fractions.

For sake of simplicity it is useful to introduce the ratio k of glucose oxidation so that

κ (2) φ = kκ (1) φ . ( 16
)
Now let us consider the rate of ATP production per cell. In the aerobic population it is

E (1) = κ (1) φ c o c g , ( 17
)
and in the anaerobic population it is

E (2) = kκ (1) φ c g . ( 18
)
We suppose that ATP is not accumulated, and the amount of ATP produced is the amount consumed. As ATP is an energy carrier molecule, these quantities (E (1) and E (2) ) define how much energy the cells may spend. We assume that some of this energy is always spent for life maintenance, and every cell spends ATP with rate θ . The rest is what is left for growth and proliferation.

We propose that if there is not enought ATP to maintain life, that is

E (1) -θ < 0 (19)
in the aerobic population, or

E (2) -θ < 0 (20)
in the anaerobic population, then some cells start dying.

We propose that only as many cells survive, as may be fed with the current rate of ATP production: [START_REF] Nelson | A fully conservative ghost fluid method stiff detonation waves[END_REF]. Anaerobic pathway is less energy efficient and in addition it acidifies the microenvironement.

φ i,surviving φ i = min E (i) θ , 1 , i = 1, 2, (21)
and the fraction of dying cells is

φ i,dying φ i = 1 - E (i) θ + , i = 1, 2, (22)
where (•) + is the positive part of (•).

Given that τ 1/2 is the "half-life" of the dying cells, the death rate in hypoxic conditions is

Γ - i = φ i ln 2 τ 1/2 1 - E (i) θ + , i = 1, 2. ( 23
)
Finally, in hypoxic conditions aerobic cells may suddenly switch to the glycolytic metabolism. We suppose that the rate of conversion from the aerobic to the anaerobic population is

Γ s = νφ 1 H( θ -E (1) ), ( 24
)
where ν is the small conversion rate, and H is Heaviside function.

If the production of ATP is excessive, some cells may grow and proliferate, and the growth rate is

Γ + i = φ i θ ln 2 Q M E (i) θ -1 + , i = 1, 2, ( 25
)
where Q M is the amount of ATP required for the full cell cycle (G1 growth, DNA synthesis, G2 growth and mitosis).

If we assume that the host tissue consumes oxygen and glucose like aerobic cells, but it does not proliferate, it cannot switch to the anaerobic metabolism and consumes only as much as it needs. Then the normal rate of ATP production in the host is: q

(h) AT P = θ φ h , ( 26
)
and the corresponding death rate is

Γ h = -φ h ln 2 τ 1/2 1 - E (h) θ + , ( 27
)
where

E (h) = κ(φ * -φ h )c o c g . ( 28
)
Now let us consider again the rates of ATP production q

(1)

AT P , q

(2) AT P and q (h) AT P . They are directly related to the rates of oxygen and glucose consumption in the tissue. In the aerobic population the rate of oxygen consumption is

q (1) o = - 6 N q
(1)

AT P , ( 29
)
and the rate of glucose consumption is

q (1) g = - 1 N q (1) AT P . ( 30
)
Also in the host tissue

q (h) o = - 6 N q (h) AT P , (31)
and

q (h) g = - 1 N q (h) AT P . ( 32
)

A c c e p t e d m a n u s c r i p t

In the anaerobic population oxygen is not consumed and the rate of glucose consumption is

q (2) g = - 1 2 q (2) AT P . ( 33
)
So the overall consumption rates for oxygen is:

q o = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ - 6 N q (1) AT P , in Ω t , - 6 N q (h) AT P , in Ω \ Ω t , (34)
and for glucose it is:

q g = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ - 1 N q
(1)

AT P - 1 2 q (2) AT P , in Ω t , - 1 N q (h) AT P , in Ω \ Ω t . ( 35
)
We substitute these reaction terms in the diffusion equations for oxygen and glucose respectively. Assuming that both oxygen and glucose may freely diffuse in all components of the mixture, we have

∂ c o ∂t = D o Δc o + q o , ( 36
) ∂ c g ∂t = D g Δc g + q g . ( 37
)
Summarising, we have the following equations:

∂ φ 1 ∂t -K∇ • (φ 1 ∇((φ 1 + φ 2 )S(φ 1 + φ 2 ))) = + φ 1 θ ln 2 Q M ⎛ ⎝ κ (1) φ θ c o c g -1 ⎞ ⎠ + -φ 1 ln 2 τ 1/2 ⎛ ⎝ 1 - κ (1) φ θ c o c g ⎞ ⎠ + -φ 1 νH ⎛ ⎝ 1 - κ (1) φ θ c o c g ⎞ ⎠ , in Ω t , (38) ∂ φ 2 ∂t -K∇ • (φ 2 ∇((φ 1 + φ 2 )S(φ 1 + φ 2 ))) = + φ 2 θ ln 2 Q M ⎛ ⎝ k κ (1) φ θ c g -1 ⎞ ⎠ + -φ 2 ln 2 τ 1/2 ⎛ ⎝ 1 -k κ (1) φ θ c g ⎞ ⎠ + + φ 1 νH ⎛ ⎝ 1 - κ (1) φ θ c o c g ⎞ ⎠ , in Ω t , (39) ∂ φ h ∂t -K∇ • (φ h ∇(φ h S(φ h ))) = -φ h ln 2 τ 1/2 1 - κ(φ * -φ h ) θ c o c g + , in Ω \ Ω t , (40) n • ∇ ((φ 1 + φ 2 )S(φ 1 + φ 2 )) = n • ∇ (φ h S(φ h )) , on ∂ Ω th , (41) (φ 1 + φ 2 )S(φ 1 + φ 2 ) = φ h S(φ h ) , on ∂ Ω th , (42) ∂ c o ∂t = D o Δc o -H(ψ) 6 N κ (1) φ φ 1 c o c g -H(-ψ) 6 N θ φ h , in Ω, (43) ∂ c g ∂t = D g Δc g -H(ψ) 1 N κ (1) φ φ 1 c o c g + 1 2 kκ (1) φ φ 2 c g -H(-ψ) c * o 6c * g αθ φ h , in Ω, ( 44
)
where ψ > 0 in Ω t , and ψ < 0 in Ω \ Ω t . For the sake of simplicity we also assume that

S(φ ) = φ -φ 0 . ( 45
)
In Appendix A we rewrite the model in dimensionless variables and find that the system is characterised by the dimensionless parameters listed in Table 1.

We would like to highlight that the model parameters still have physiological meaning and depend on measurable quantities. The growth rate is proportional to the intensity of metabolic processes κ and is inversely proportional to the cost of the whole cell cycle Q M (probably measurable but to our awareness not known). The death rate depends not only on the intensity of metabolism, but also on the cost of maintenance ( θ ), and tolerance of the cells to the extreme conditions ("half-life" τ 1/2 of the dying cells). The same κ affects both the consumption rates of the nutrients and the growth/death rates of the cells.

The proposed approach of incorporating tumour metabolism in the model is different from the approach taken in [START_REF] Venkatasubramanian | Incorporating energy metabolism into a growth model of multicellular tumor spheroids[END_REF][START_REF] Bertuzzi | ATP production and necrosis formation in a tumour spheroid model[END_REF]. Those works essentially assume that cells use glycolytic metabolism as a fallback in hypoxic conditions and return to normal metabolism when the level of oxygen increases. So the same cells use two metabolic pathways depending on the current conditions, there are no memory effects in the system. In contrast, in this model a new population of cells is born, with a different and sustainable behaviour. We think that such a modelling approach has an advantage as it allows more straightforward translation of genetic transformations in tumour into modelling entities (cell populations), and is not limited to metabolic switch only. A similar approach can also be used for reversible epigenetic transformations and in particular for the inclusion of epigenetic factors in metabolic switches.

Here in the description of the Warburg effect we consider only irreversible transitions to the glycolytic metabolism, which are in fact observed experimentally [START_REF] Kim | Cancer's molecular sweet tooth and the Warburg effect[END_REF].

A c c e p t e d m a n u s c r i p t 2 Simulations

In some previous works we have already studied some aspects of the growth of tumour cords (Astanin & Tosin, 2007;[START_REF] Preziosi | Multiphase modeling of tumor growth and extracellular matrix interaction: mathematical tools and applications[END_REF][START_REF] Astanin | Mathematical model of tumour cord growth along the source of nutrient[END_REF] in particular in the case of mono-population models where growth only depends on oxygen supply. This form of growth results in a long cylindrical tumour mass arranged around a blood vessel, which forms the axis of the cord. Cells mostly proliferate near the axis and the tip of the cord. The hypoxic zone is formed in the outer rim of the cord. We will see here that also the bipopulation cord grows similarly to the mono-population one, and therefore we study its growth in a similar geometry.

We consider two geometrical settings in this section. Both simulate the growth of the cord in a two-dimensional domain along the blood vessel. The first setting focuses on the part of the cord with stationary radius far away from its tip (Fig. 2 a). This setting is similar to known uni-dimensional simulations of tumour cords [START_REF] Bertuzzi | Cell kinetics in a tumour cord[END_REF][START_REF] Bertuzzi | Cell kinetics in tumour cords studied by a model with variable cell cycle length[END_REF][START_REF] Bertuzzi | Regression and regrowth of tumour cords following single-dose anticancer treatment[END_REF][START_REF] Bertuzzi | A mathematical model for tumor cords incorporating the flow of interstitial fluid[END_REF]. The second setting focuses on the axial evolution and on the tip of the expanding cord (Fig. 2 b).

In all simulations we neglected consumption of oxygen and glucose by the host tissue ( 26), thus disabling also the death term ( 27) in the equation for the host cells.

One may note, that we have actually two different model for two domains, and even the number of variables in two domains do not coincide. We discuss boundary conditions between two domains below. To impose these boundary conditions numerically, we used a custom modification of the ghost fluid method. This modification is described in Appendix B.

Both numerical simulations share the same code for the numerical method. The source code was published under an open source license (GPL) at http://code.google.com/p/ cord/.

Typical values of model parameters used in the simulations are given in Table 1. Estimates for some of the model parameters may be found in [START_REF] Venkatasubramanian | Incorporating energy metabolism into a growth model of multicellular tumor spheroids[END_REF]. Stressfree volume fraction of cells φ 0 = 0.75 was estimated from π/4 area coverage of the 2D plane with fixed size circles. Choice of the time and spatial scales is discussed in Appendix A.

Boundary conditions

There are five boundaries in the model: the four exterior boundaries

∂ Ω W , ∂ Ω E , ∂ Ω S , ∂ Ω N of
the model domain Ω, and the boundary between Ω t and Ω \ Ω t , or tumour-host interface ∂ Ω th (Fig. 2).

We distinguish three roles of an exterior boundary:

• a boundary coinciding with a blood vessel,

• a remote boundary,

• a symmetry axis.

For boundaries coinciding with a vessel, we assume that the cells do not penetrate the wall of the vessel, and nutrient supply is always sufficient to maintain the constant concentrations c 0 o , c 0 g . So, for x at vascular boundaries, we have:

⎧ ⎨ ⎩ ∂ (φ S(φ )) ∂ n ( x,t) = 0, c k ( x,t) = c 0 k , k = o, g (46)
where n is a unit exterior normal of Ω. These boundaries are the only source of oxygen for the region of tissue in consideration. These boundary conditions are written under the assumption that oxygen diffusion through the blood vessel wall is fast with respect to the growth of tumour. A more elaborate approach to describe tumour-vessel interaction is suggested in [START_REF] Breward | Modelling the interactions between tumour cells and a blood vessel in a microenvironment within a vascular tumour[END_REF].

For remote boundaries we assume that they stay undisturbed by the growth and the flux of nutrients through them is zero. For x at remote boundaries, then

⎧ ⎨ ⎩ S(φ )( x,t) = 0, ∂ c k ∂ n ( x,t) = 0, k = o, g (47)
In some simulations we also assume growth to be symmetrical with respect to one of the domain boundaries, and for that boundary

⎧ ⎪ ⎨ ⎪ ⎩ ∂ φS(φ ) ∂ n ( x,t) = 0, ∂ c k ∂ n ( x,t) = 0, k = o, g (48)
The tumour-host interface ∂ Ω th coincides with the material boundary of tumour cells fractions and moves with the velocity of the cells v. Stress and velocity are continuous across the tumour-host interface ∂ Ω th . If the function S(φ ) is the same for both tissues then the condition above implies continuity of φ 1 + φ 2 and φ h . Conditions on the free boundary are defined in non-dimensional form equations ( 51), (53-55).

We describe our way to impose these boundary conditions numerically in Appendix B. This is not trivial because we need to interface the solutions in two domains where the number of variables do not coincide. However, the fact that the stresses and the velocities of the components and of the interface depend on the sum of the variables makes the problem well formulated and allows to extend the ghost fluid method as discussed in Appendix B.

Body of the cord and effects of the model parameters

In order to put in evidence some characteristics of the glycolytic switch, let us initially consider the simplified initial configuration, shown in Figure 2 a. This corresponds to starting with a long cord, with all cells using aerobic pathway, and focusing on its central region. Even though the central region cannot expand along the vessel but may only adjust its width and its structure may change according to the local condition. This configuration is essentially one-dimensional, but it allows better understanding of the various model parameters and their effects. Simulations were run on a narrow rectangular domain [0, 0.01] × [0, 1] with finite difference grid 5×500. Initially the tumour-host interface was positioned at y = 0.01. The same 2D code was used in these simulations as in the subsequent ones, but we studied growth only in one dimension along the y-axis.

Referring to Appendix A for the deduction of the dimensionless model and to Table 1 for the definition of the dimensionless quantities, we focus on the effect of the parameters related to metabolism.

Figure 3 shows a parametric study of the model. In each series of simulations only one parameter was changed while the others were fixed on the values from Table 1. The first two columns of the figure show how percentage of the glycolytic

cells p = 100 • φ 2 dx (φ 1 + φ 2 )dx
changes in the tumour for various values of α, k and ν. For the sake of clarity we omit the tildes from now on. The right column shows how these parameters affect the stationary width of the cord. The parameter α is the oxygen consumption rate and reflects the general intensity of metabolic processes. We observed that tumour cords with more intensive metabolism suffer from hypoxia earlier (Figure 3 for various values of model parameters. We observed, that for long simulation times anaerobic cells spread through the whole tumour, but never completely extinguish the original aerobic population. The left column (a, d, g) shows emergence of the anaerobic population with time. We observe that the onset of anaerobic switch is independent of k (d). The middle column (b, e, h) shows the final volume fraction of the anaerobic cells p ∞ as a function of model parameters. Slow oxygen uptake (low α) delays emergence of the anaerobic cells (a), but favours them in the long run (b). Relative rate of the anaerobic metabolism k determines composition of the tumour, cells with faster metabolism dominate (d, e). Metabolic switch from aerobic to anaerobic metabolism does happen even for very low conversion rates ν ("mutation rate"), and the significant number of cells is converted even for very low ν (h). We may also observe an approximate dependence like w ∞ ∝ 1/ 3 √ α for the thickness of the cord. Cords where anaerobic metabolism is faster are smaller (f). This fact may be explained by higher fraction of anaerobic cells, and shorter diffusion distance of glucose with respect to that of oxygen. Probably for the same reason cords with faster transition to anaerobic metabolism are smaller (i). Curves are automatically fitted and serve for visual aid only. depends on the rate of oxygen consumption approximately as w ∞ ∝ 1/ 3 √ α. Now let us consider the effects of two parameters specific for the bi-population model. They are ν, the conversion rate from aerobic to anaerobic metabolism, and k, the relative rate of the anaerobic metabolism with respect to aerobic one.

These parameters have a clear physiological meaning: ν is related to the mutation rate and the genetic instability (the higher it is, the more probable the mutation is), and k is related to the time required for a molecule of glucose to pass through all the catabolic reactions' chain. While we do not have definitive real-world estimations for these parameters at this moment, we expect ν to be small and k > 1, because the anaerobic pathway is much shorter than the aerobic one.

Simulation results for various values of the parameter k are presented in Figure 3 Figures 3 g,h, i represent simulations with various values of the parameter ν. Even with very small ν the tumour is successfully converting to glycolysis. However, the value of ν is affecting the steepness of the initial exponential stage of growth of the anaerobic subpopulation, and the duration of this period (Fig. 3 g). In the extreme case, this conversion may be even so fast, that a further correction may happen (see the results for ν = 0.2 in Fig. 3 g, the upper curve is not monotonous).

Smaller ν means that more tumour cells will die immediately due to hypoxia, while larger ν means that more cells may quickly adapt to anaerobic catabolism and survive. Their quantity is still limited by availability of glucose.

An open question is whether in the long term just one subpopulation (clone) remains or both subpopulations will coexist forever. If both subpopulations remain then what will be the ratio between them? Figures 3 b,e, h suggest that both population coexist, and the ratio does depend on the conversion rate ν, as well as other model parameters. Interestingly, even for very low conversion rates ν the anaerobic population plays significant role in the tumour.

To sum up, population with faster metabolism "wins", intensity of metabolism affects cord geometry and the moment when the tumour starts suffering from hypoxia, anaerobic population may achieve significant mass even with very low conversion rates.

Tip of the cord and spatial dynamics

Now let us consider the full-featured two-dimensional simulation of the cord growth in the configuration presented in Figure 2 b. We considered growth of a single cord in the domain [0, 1] × [0, 0.5]. In the simulations the relative rate of the anaerobic glucose catabolism k was taken to be k = 2.0. This value reflects the assumption, that the anaerobic pathway is shorter and faster in tumour cells than the complete aerobic oxidation of glucose. The other model parameters were as in Table 1.

We focus our attention on the energy supply level per cell, which is e 1 = E (1) -θ for the aerobic population, and e 2 = E (2) -θ for the anaerobic population. When these quantities become negative, the cells suffer from energy starvation, which leads to death or to metabolic switch. When these quantities are positive, the cells may proliferate (see Figure 4).

We shall refer to the isolines e 1 ≡ 0 and e 2 ≡ 0 as starvation limits for the aerobic and anaerobic population respectively.

The results are presented in Figure 5. Qualitatively the growth is very similar to the growth of a genetically homogeneous cord (Astanin & Tosin, 2007), with the radial size of the cord limited by diffusion distance of the nutrients (including glucose in this case). The dynamics, however, shows relevant differences involving two sub-populations. Initially the cord is composed of aerobic cells only (Fig. 5 and intensively grows, until some cells start suffering from hypoxia (Fig. 5 b). Hypoxic region is coloured on the plots, and the starvation limit for the aerobic population is shown with a point-and-dashed red line. However, as soon as anaerobic cells appear, we should also consider another starvation limit, above which also anaerobic cells suffer (dashed blue line, Fig. 5 c,d).

∂ Ω N ∂ Ω E e 1 = 0 (for aerobic) e 2 = 0 (for anaerobic) φ 1 ↑, φ 2 ↑ φ 1 ↓, φ 2 ↑ φ 1 ↓, φ 2 ↓ 3 2 1 ∂ Ω th ∂ Ω W ∂ Ω S : vessel Ω \ Ω t Ω t x y
We would like to note, that the positions of the starvation limits do not coincide. Moreover, with k > 1, the position of the starvation limit for anaerobic cells is further from the vessel than the position of the starvation limit for the aerobic cells.

Figure 4 schematizes the positions of the starvation limits. There are three distinct zones. In the first zone 1 near the blood vessel both populations have enough energy supply, and can proliferate. The aerobic population (φ 1 ) suffers from hypoxia in the second 2 and in the third 3 zones. However, anaerobic population (φ 2 ) suffers only in the third zone. This means that in the zone 2 anaerobic cells have a clear advantage over those using an aerobic metabolism.

One may expect that anaerobic cells would occupy zones 2 and 3 , while the aerobic cells would stay in the zone 1 . However, we did not observe this in the simulations. Instead, while the anaerobic cells initially emerge only in the hypoxic region (Fig. 6 We attribute this fact to the broken equilibrium in the tissue when the anaerobic population suddenly emerges. In our model, the cells move as a whole, with the same velocity (51), which depends on the gradient of the total volume fraction of cells ∇(φ 1 + φ 2 ). When the anaerobic cells appear, they do not feel their proliferation constraints initially, and grow much faster than their aerobic neighbours. However, the local conditions for the neighbour aerobic cells do not change or even become more severe, because the supply of glucose also reduces. This means that given two parts of the tissue in contact, one with anaerobic cells, and the other with only aerobic cells, the total volume fraction of cells will be higher in the tissue with the anaerobic cells, and this will lead to flow of the cells from the region with anaerobic cells into the region without them. This effect is clearly visible in Figure 6 a, where the cells in region 2 , as defined in Figure 4, expand almost radially. Velocity of the axial expansion slightly increases (Fig. 7).

Interestingly, a similar spreading of the glycolytic cells was described in the models based on cellular automata [START_REF] Smallbone | Metabolic changes during carcinogenesis: potential impact on invasiveness[END_REF][START_REF] Gatenby | Cellular adaptations to hypoxia and acidosis during somatic evolution of breast cancer[END_REF]. And in fact, those models

A c c e p t e d m a n u s c r i p t

share some common properties with the model presented here: cells die when ATP production is insufficient, and cells divide only if ATP is produced in excess. In both cases glycolytic cells have an advantage. These models in [START_REF] Smallbone | Metabolic changes during carcinogenesis: potential impact on invasiveness[END_REF][START_REF] Gatenby | Cellular adaptations to hypoxia and acidosis during somatic evolution of breast cancer[END_REF], however, do not take into account flow of the tissue, but consider acidity induced cell death, so the direct comparison of the results is not possible.

Concluding Remarks

In our model we described emergence of one of the hallmarks of cancer, the ability of the tumour to rely exclusively on glycolytic metabolism. In the model we assumed that some cells may spontaneously acquire this property in hypoxic conditions. We focused on the transition of the tumour from one type of metabolism to the other on a model of a growing tumour cord.

In this model we managed to build a reasonable link between the deduction of growth-death processes and consumption of the nutrients. In bi-population model, this approach allowed to naturally describe alternative metabolic pathways and different regions of energy starvation for different cell populations.

We believe that the proper description of the cancer development requires taking into account of both spatial phenomena in the tissue and evolution of the tumour cells themselves. This means that the cell lines with different behaviour and attributes may emerge during tumour development, and their emergence should be addressed by tumour growth models. The incorporation of population dynamics with different sub-cellular behaviours into continuous models is one of the ways to achieve that.

However, even dealing with two cell lines in such models poses already many technical difficulties. In our case we were able to solve them. Our modification of the ghost fluid method allows to stitch solutions between two porous media domains when the number of the components in those domains do not coincide, using the fact that the velocities of the components and of the boundary between the domains depend on the sum of the variables. Applicability and reliability of this method needs to be studied.

We described how the parameters related to metabolism affect the composition and width of the cord. It is interesting to notice that the width of the cord seems to go like the inverse of the cubic root of the oxygen uptake coefficient α. We also showed that the metabolic switch induces a tip-core heterogeneity rather than an axysymmetric stratification, as also described phenomenologically. This is due to the fact that the anaerobic population has its proliferation limits offset, grows relatively faster and invades the entire cord. The inclusion of epigenetic factors of the metabolic transition would probably induce a slight radial stratification that would be stronger if only epigenetic factors are considered. This issue, however, deserves further investigation.

The presented model may certainly benefit from some other extensions. Firstly, switch to glycolytic metabolism implies acidification of the environment. Ability to sustain its toxi-city is known as another very important hallmark of cancer, and is fundamental in acidity mediated invasion. Secondly, appearance of anchorage-independent and invasive cells are other phenomena which may greatly extend the applicability of the model. Using qualitatively accurate model parameters in the simulation is another priority.

A Model non-dimensionalization

Following Section 1 we may write the complete model:

∂ φ 1 ∂t + ∇ • (φ 1 v t ) = Γ + 1 + Γ - 1 -Γ s , in Ω t , ( 49
) ∂ φ 2 ∂t + ∇ • (φ 2 v t ) = Γ + 2 + Γ - 2 + Γ s , in Ω t , ( 50
)
v t = -K∇((φ 1 + φ 2 )S(φ 1 + φ 2 )) , in Ω t , (51) ∂ φ h ∂t + ∇ • (φ h v h ) = Γ - h , in Ω \ Ω t , ( 52
)
v h = -K∇(φ h S(φ h )) , in Ω \ Ω t , ( 53
)
v t • n = v h • n , on ∂ Ω th , ( 54
) (φ 1 + φ 2 )S(φ 1 + φ 2 ) = φ h S(φ h ) , on ∂ Ω th , (55) ∂ c o ∂t = D o Δc o + q o , in Ω, (56)
∂ c g ∂t = D g Δc g + q g , in Ω. (57)
Now let us expand the reaction terms and rewrite the model in dimensionless variables. We use tildes to distinguish them from the non-scaled variables.

Given the maximum available volume fraction of cells is φ * , we rescale cells' volume fractions as

φ i = φi φ * , i := 1, 2, h, (58)
and for concentration scales c * o and c * g of oxygen and glucose respectively (mass or moles per volume) we have:

c o = co c * o , ( 59
)
c g = cg c * g . ( 60
)
We also define the scaled S as

S( φ ) = S( φ φ * ) φ * , ( 61
)
and define dimensionless κφ as κφ ( φ1 , φ2 ) = κ

(1)

φ ( φ1 φ * , φ2 φ * ) φ * κ = 1 -φ1 -φ2 . ( 62
)
We choose as time scale:

T = Q M κφ * c * o c * g ln 2 , ( 63
)
and as spatial scale:

L = Q M D o κφ * c * o c * g ln 2 . ( 64
)

A c c e p t e d m a n u s c r i p t

Then the equations may be rewritten as

∂ φ1 ∂ t -K ∇ • φ1 ∇(( φ1 + φ2 ) S( φ1 + φ2 )) = + φ1 κφ co cg -θ + -ε φ1 θ -κφ co cg + -ν φ1 H θ -κφ co cg , in Ω t , (65) ∂ φ2 ∂ t -K ∇ • φ2 ∇(( φ1 + φ2 ) S( φ1 + φ2 )) = + φ2 k κφ cg -θ + -ε φ2 θ -k κφ cg + + ν φ1 H θ -κφ co cg , in Ω t , (66) ∂ φh ∂ t -K ∇ • φh ∇( φh S( φh )) = -ε φh θ -(1 -φh ) co cg + , in Ω \ Ω t , (67) n • ∇ ( φ1 + φ2 ) S( φ1 + φ2 ) = n • ∇ φh S( φh ) , on ∂ Ω th , ( 68
)
( φ1 + φ2 ) S( φ1 + φ2 ) = φh S( φh ) , on ∂ Ω th , ( 69
)
∂ co ∂ t = Δ co -H(ψ) α κφ φ1 co cg -H(-ψ) α θ φh , in Ω, (70)
∂ cg ∂ t = Dg Δ cg -H(ψ) cog α κφ φ1 co + N 2 k φ2 cg -H(-ψ) cog α θ φh , in Ω, (71)
where ψ > 0 in Ω t , and ψ < 0 in Ω \ Ω t . For the sake of simplicity we also assume that

S( φ ) = φ -φ0 , ( 72
)
and the dimensionless coefficients are defined as follows:

• cell motility K = T Kφ 2 * /L 2 = Kφ 2 * /D o , • death rate ε = T φ * c * o c * g κ ln 2/ θ τ 1/2 = Q M /( θ τ 1/2 ), • conversion rate ν = νQ M /(κφ * c * o c * g ln 2), • upkeep cost per cell θ = θ /(κφ * c * o c * g ), • oxygen uptake rate α = 6T φ 2 * c * g κ/N = 6Q M φ * /(Nc * o ln 2), • relative rate of glycolysis k = k/c * o , • glucose diffusion coefficient Dg = D g T /L 2 = D g /D o ,
• ratio between absolute concentrations cog

= c * o 6c * g .
Once the model is defined, we omit tildes in (38-44) everywhere in the paper, except for Table 1.

B Modified ghost fluid method

In our model we have to deal with a free boundary separating two domains, where each of the phase variables is defines in only one domain.

There is a technique to solve similar kind of problem with moving interface and discontinuities, where the level set method is used in conjunction with the ghost fluid method [START_REF] Fedkiw | A nonoscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Fedkiw | The ghost fluid method for viscous flows[END_REF]Nguyen et al., 2002;[START_REF] Macklin | An improved geometry-aware curvature discretization for level set methods: Application to tumour growth[END_REF].

The level set method [START_REF] Sethian | The design of algorithms for hypersurfaces moving with curvature-dependent speed[END_REF] is used to distinguish between the domains. It consists in introducing an additional variable ψ that is initially positive in one domain and negative in the other. The ghost fluid method consists in extrapolating each variable in the neighbour domain across the interface in a way which enforces the conditions on the interface. If I is continuous across the interface, then it may be extrapolated in the neighbour domain with

∂ I ∂t ± n • ∇I = 0. ( 73
)
Just a few time steps are typically sufficient as only values of I near the interface are required. They are used to extrapolate the variable of interest φ i in the neighbour domain.

∂ Ω th tumour host

φ 1 + φ 2 φ 2 φ 1 φ h
Figure 8: Cell packing density near the tumour-host interface. Neither variable is defined in the complement domain.

The standard single-phase solvers may be used to solve the equations for both phases separately, using extrapolated "ghost" phase from the other domain. Then the actual noncontinuous solution may be reconstructed in each of the domains using the values of one of the variables.

However, most known applications of the ghost fluid method deal with one variable per domain and use the same equation in both domains. In our model we have to evaluate In particular, it is not trivial to enforce the conditions on the tumour-host interface, because neither variable is continuous across the interface (Fig. 8) nor even defined in the complement domain.

We propose an extension of the ghost fluid method for this particular case when the number of phases in the liquid on one side of the interface is different from the number of phases on the other side.

Let us note, that the method should maintain not only the continuity of stress but also the continuity of the velocity of the interface, which both depend on the combination of the phase variables φ 1 + φ 2 through ( 51), (53-55), but each phase is governed by a separate equation. It is important to apply restrictions on the interphase ratios to satisfy the requirements above when extrapolating the phase variables in the ghost domain.

Let us assume that there are two phases on one side of the interface, and one phase on the other. Now we can rewrite the continuity requirements in a much simpler form:

(φ 1 + φ 2 )| T = φ h | H , ( 74
) n • ∇(φ 1 + φ 2 )| T = n • ∇φ h | H , ( 75
)
where |T and |H denote values on different sides of the interface.

Referring to Figure 9, let us consider a point x on the interface and the line normal to the interface. Let us assume that the initial configuration is conforming, i.e. the conditions above are satisfied. We introduce a phase ratio variable

∂ Ω th tumour host φh = φ 1 + φ 2 φ 1 φ 2 φ h r = φ 1 /(φ 1 + φ 2 ) r φ1 = rφ h φ2 = (1 -r)φ h
r = φ 1 φ 1 + φ 2 . ( 76
)
A bar above the variable denotes extrapolated (ghost) values. An upper index 0 denotes values on the interface:

φ 0 := (φ 1 + φ 2 )| T = φ h | H , (77) r 0 := r| T . ( 78
)
We do not need to extrapolate the stress explicitly as it is uniquely identified by the sum of the phase variables when the function S is the same in tumour and host: S t (φ ) = S h (φ ). We use (73) to extrapolate r in the host domain. Practically, we advect the value of the phase ratio on the interface into the neighbouring domain, and it is approximately equal to the value on the interface: r

( x + ε n)| H = r 0 + O(|ε|) (79)
Then we define extrapolated (ghost) values of φ 1 and φ 2 in the host domain as φ1 := rφ h and φ2 :

= (1 -r)φ h , ( 80
)
while φ h in the tumour domain is defined as

φh := φ 1 + φ 2 . ( 81
)
Evidently, equation ( 73) transports the value on the interface along the normal to the interface into the other domain. This naturally raises the question how to extrapolate the values near the concave points of the interface where the normals to the interface intersect. Referring to Figure 10, let the interface be defined locally as y = f(x), the normal to the interface be (-f (x), 1) T . Let us consider two points of the interface: x 0 = (x 0 , y 0 ), x 1 = (x 1 , y 1 ), and the point x = (x, y) at the intersection of the normals to the interface. We may express y as

y = y 0 - x -x 0 f (x 0 ) , ( 82
)
and as

y = y 1 - x -x 1 f (x 1 ) . ( 83
)
Then where Δx = x 1x 0 , (87)

y 1 -y 0 = 1 f (x 1 ) - 1 f (x 0 ) x - x 1 f (x 1 ) + x 0 f (x 0 ) , (84)
we have at the zero-th order

x ≈ x 0 - [1 + (f (x 0 )) 2 ]f (x 0 ) f (x 0 ) = x 0 - f (x 0 ) κ(x 0 ) 1 + (f (x 0 )) 2 ,
(88) and

y ≈ y 0 - [1 + (f (x 0 )) 2 ] f (x 0 ) = y 0 - 1 κ(x 0 ) 1 + (f (x 0 )) 2 , ( 89
)
where κ(x 0 ) is the curvature of the interface in point x 0 .

Numerically we need at least one grid point at each of the segments ( x 0 , x) and ( x 1 , x) to avoid uncertainty in the solution of ( 73). This gives an estimation of the maximal spatial step of the numerical grid:

h x ≤ 1 2 f (x 0 ) κ(x 0 ) 1 + (f (x 0 )) 2 , ( 90
)
h y ≤ 1 2 1 κ(x 0 ) 1 + (f (x 0 )) 2 . ( 91
)
With our choice of ghost values the extrapolated phase variables are continuous, and the total density in both domain is preserved as well as the stress. However, it is easily seen that gradients of φ i and φi are not continuous. On the side of the host domain

n • ∇ φ1 | H = n • ∇(r 0 φ h )| H = r 0 n • ∇φ h | H = r 0 n • ∇(φ 1 + φ 2 )| T , (92)
while on the side of the tumour domain

n • ∇φ 1 | T = ((φ 1 + φ 2 ) n • ∇r + r n • ∇(φ 1 + φ 2 ))| T = φ 0 n∇r + n∇φ 1 | H . (93)
This gives the jump of the φ 1 gradient across the interface:

n • (∇φ 1 | T -∇φ 1 | H ) = φ 0 n • ∇r| T . ( 94
)
In the same way, for φ 2 we have

n • (∇φ 2 | T -∇φ 2 | H ) = φ 0 n • ∇r| T . ( 95
)
Fortunately, we do not need values of ∇φ 1 and ∇φ 2 separately, but only ∇(φ 1 + φ 2 ), which determines the velocity of the interface (51), (53)(54)(55). With this choice of φ1 , φ2 and φ , the velocity of the tissue is continuous across the interface.

Evidently, for a monophase domain r ≡ 1. This practically means that for two monophase domains with the same S we do not have to do anything: r ≡ 1 everywhere, and the ghost variables coincide with the real variables.

This approach may be generalized for two domains A and B with S A (φ ) = S B (φ ) and greater number of subphases, when all phases move with the same velocity, and velocity and stresses depend on the sum of the phase variables. Given n the number of phases in the domain A, and m be the number of phases in the domain B, the algorithm to define "ghost" variables is represented in Table 2. Once this procedure is defined the construction of the complete numerical method is relatively straightforward (Table 3

Figure 1 :

 1 Figure1: Glucose catabolic pathways, based on[START_REF] Nelson | A fully conservative ghost fluid method stiff detonation waves[END_REF]. Anaerobic pathway is less energy efficient and in addition it acidifies the microenvironement.

Figure 2 :

 2 Figure 2: Simulation domain and boundary conditions for the body of the cord (a) and for the tip of the cord (b)

  a), are smaller (Figure 3 c) and have lower fraction of glycolytic cells (Figure 3 b).It is interesting to note that the stationary cord thickness

Figure 3 :

 3 Figure 3: Percent of anaerobic cells p = 100 • φ 2 d x/ (φ 1 + φ 2 ) d x and stationary thickness of the bi-population tumour cord w ∞for various values of model parameters. We observed, that for long simulation times anaerobic cells spread through the whole tumour, but never completely extinguish the original aerobic population. The left column (a, d, g) shows emergence of the anaerobic population with time. We observe that the onset of anaerobic switch is independent of k (d). The middle column (b, e, h) shows the final volume fraction of the anaerobic cells p ∞ as a function of model parameters. Slow oxygen uptake (low α) delays emergence of the anaerobic cells (a), but favours them in the long run (b). Relative rate of the anaerobic metabolism k determines composition of the tumour, cells with faster metabolism dominate (d, e). Metabolic switch from aerobic to anaerobic metabolism does happen even for very low conversion rates ν ("mutation rate"), and the significant number of cells is converted even for very low ν (h). We may also observe an approximate dependence like w ∞ ∝ 1/ 3 √ α for the thickness of the cord. Cords where anaerobic metabolism is faster are smaller (f). This fact may be explained by higher fraction of anaerobic cells, and shorter diffusion distance of glucose with respect to that of oxygen. Probably for the same reason cords with faster transition to anaerobic metabolism are smaller (i). Curves are automatically fitted and serve for visual aid only.

  d, e, f (the second row). It appears that k is one of the most important parameters in the model. It has decisive influence on what population (metabolic pathway) dominates in the tumour. If glycolysis is faster (k > 1) than normal metabolism, the glycolytic cells may compose up to 85% of the On the opposite, if pure glycolytic metabolism were for some reason slower (k < 1), then glycolytic cells are only a marginal fraction of the tumour mass (as low as 5% in the simulations).This difference in composition of the tumour may explain slightly different sizes of the cord for various values of k (Figure3 f).

Figure 4 :Figure 7 :

 47 Figure 4: Energy deficit limits, growth (↑) and death (↓) in two populations. There are three distinct zones.

  a), later they spread throughout the tumour (Fig.6 b).

Figure 5 :Figure 6 :

 56 Figure5: Bi-populational tumour cord growth, starvation limits, colour/darkness shows negative part of e 1 = E (1) -θ (deficit of ATP for aerobic cells, i.e. severity of hypoxia)

  one domain and one in the other, and the equations are different in different domains.

Figure 9 :

 9 Figure 9: Extrapolating component densities across the tumour-host interface. Interphase ratio r = φ 1 /(φ 1 + φ 2 ) is extrapolated continuously (73) and the stresses φ h S(φ h ) and (φ 1 + φ 2 )S(φ 1 + φ 2 ) are preserved. Extrapolated values shown with dotted lines.

Figure 10 :

 10 Figure 10: Extrapolation near the concave points of the interface.

  1 )f (x 0 )(y 1y 0 ) + f (x 0 )x 1f (x 1 )x 0 f (x 0 )f (x 1 ) . (85) Substituting f (x 1 ) = f (x 0 ) + f (x 0 )(Δx) + o((Δx)),(86)

Table 1 :

 1 Typical values of the model parameters for the simulations of the metabolic switch in the tumour. See also model definition in Section 1 and non-dimensionalization of the model in Appendix A.

	Nondimensional	Definition	Description	Value
	S( φ ) φ0 K	S( φ φ * )/φ * φ0 φ * K φ 2 * /D o	partial stress in the tissue stress-free packing density cells' motility (rate of tissue relaxation)	S = φ -φ 0 0.75 0.01
	γ	-	growth rate	1.0
	ε ν θ α Do	death rate g ln 2) aerobic-to-anaerobic conversion rate Q M /( θ τ 1/2 ) ν Q M /(κφ * c * o c * θ /(κφ * c * o c * upkeep cost per cell g ) 6Q M φ * /(Nc * oxygen uptake rate o ln 2) -oxygen diffusion coefficient	0.8 0.01 0.15 200 1.0
	Dg k L	D g /D o k /c * o (64)	glucose diffusion coefficient relative rate of glycolytic metabolism characteristic length	0.1 2.0, 4.0 1.0
	T	(63)	characteristic time	1.0
	co cg φi	c o /c * o c g /c * g φ i /φ	scaled oxygen concentration scaled glucose concentration	--

* scaled volume fraction of the i-th population -

Table 2 :

 2 ).£ Given: £ domain A, components u 1 ,...u n , £ stress function S A (u 1 + ...+ u n ); £ domain B, components v 1 ,...v m , £ stress function S B (v 1 + ... + v m ); 1 for x ∈ A £ In domain A do 2 £ evaluate the sum of components s A = ∑ A u i s A ← u 1 + ... + u n 3 £ evaluate the stress T T ← S A (s A ) 4 £ evaluate the "ghost" value sB sB ← x, such that s A S A (s A ) = xS B (x) 5£ and phase ratios of each componentr i = u i /s A for i ← 1 to n do r i ← u i /s A 6 for x ∈ B £ Do the same in domain B do s B ← v 1 + ... + v m T ← S B (s B ) sA ← x, such that xS A (x) = s B S B (s B ) for j ← 1 to m do r j ← v j /s B 7 £ Extrapolate phase ratios £ in the complement domains for i ← 1 to n do ri ← solution of (73) in B for j ← 1 to m do r j ← solution of (73) in A 8 £ Define "ghost" values of each component £ in the complement domains for x ∈ A do for j ← 1 to m do v j ← r j sB for x ∈ Bdo for i ← 1 to n do ūi ← ri sA Multiphase ghost fluid method: defining "ghost" variables in the complement domains. £ define the "ghost" variables φ1 , φ2 , φh £ in complement domains (see Fig.2and Table2) 4£ re-define model variables globallyφ i ← φ i , x ∈ Ω t , φi , x ∈ Ω \ Ω t for i = 1, 2 φ h ← φh , x ∈ Ω t , φ h , x ∈ Ω \ Ω t 5 £ evaluate solutions of φ 1 , φ 2 , φ h , c o , c g in Ω £ after time step δt φ i ← φ i (t + δt, φ 1 , φ 2 , φ h , c o , c g ), i = 1, 2, h c g ← c g (t + δt, φ 1 , φ 2 , φ h , c o , c g ) c o ← c o (t + δt, φ 1 , φ 2 , φ h , c g ) 6£ reconstruct total density variable globallyφ ← φ 1 + φ 2 , x ∈ Ω t , φ h , x ∈ Ω \ Ω t 7£ evaluate collective velocity v ← -K∇ (φ S(φ )) 8£ update the solution for the level set function ψ £ (position of ∂ Ω th )

	1 while t < T
	2	do
	3
		A c c e p t e d	m	a n u s c r i p t

ψ ← ψ(t + δt) t ← t + δt

Table 3 :

 3 Numerical method with multiphase ghost fluid method