A fully continuous individual-based model of tumor cell evolution

Pablo Gómez-Mourelo; Eva Sánchez; Luis Casasús; Glenn F. Webb

doi:10.1016/j.crvi.2008.08.010

Biological modelling / Biomodélisation

A fully continuous individual-based model of tumor cell evolution

Présenté par : Pierre Auger

Pablo Gómez-Mourelo ¹ ; Eva Sánchez ¹ ; Luis Casasús ¹ ; Glenn F. Webb ²

¹ Dpto. Matemática Aplicada. ETSI Industriales (Universidad Politécnica de Madrid), c/ José Gutiérrez Abascal, 2, 28006 Madrid, Spain
² Department of Mathematics, Vanderbilt University, Nashville, TN 37240, USA

Comptes Rendus. Biologies, Volume 331 (2008) no. 11, pp. 823-836.

Résumé

The aim of this work is to develop and study a fully continuous individual-based model (IBM) for cancer tumor invasion into a spatial environment of surrounding tissue. The IBM improves previous spatially discrete models, because it is continuous in all variables (including spatial variables), and thus not constrained to lattice frameworks. The IBM includes four types of individual elements: tumor cells, extracellular macromolecules (MM), a matrix degradative enzyme (MDE), and oxygen. The algorithm underlying the IBM is based on the dynamic interaction of these four elements in the spatial environment, with special consideration of mutation phenotypes. A set of stochastic differential equations is formulated to describe the evolution of the IBM in an equivalent way. The IBM is scaled up to a system of partial differential equations (PDE) representing the limiting behavior of the IBM as the number of cells and molecules approaches infinity. Both models (IBM and PDE) are numerically simulated with two kinds of initial conditions: homogeneous MM distribution and heterogeneous MM distribution. With both kinds of initial MM distributions spatial fingering patterns appear in the tumor growth. The output of both simulations is quite similar.

Métadonnées

Reçu le : 2008-04-22
Accepté le : 2008-08-12
Publié le : 2008-09-27

PMID

DOI : 10.1016/j.crvi.2008.08.010

Mots clés : Individual-based models, Computational mathematics

Affiliations des auteurs :

Pablo Gómez-Mourelo ¹ ; Eva Sánchez ¹ ; Luis Casasús ¹ ; Glenn F. Webb ²

@article{CRBIOL_2008__331_11_823_0,
     author = {Pablo G\'omez-Mourelo and Eva S\'anchez and Luis Casas\'us and Glenn F. Webb},
     title = {A fully continuous individual-based model of tumor cell evolution},
     journal = {Comptes Rendus. Biologies},
     pages = {823--836},
     publisher = {Elsevier},
     volume = {331},
     number = {11},
     year = {2008},
     doi = {10.1016/j.crvi.2008.08.010},
     language = {en},
}

TY  - JOUR
AU  - Pablo Gómez-Mourelo
AU  - Eva Sánchez
AU  - Luis Casasús
AU  - Glenn F. Webb
TI  - A fully continuous individual-based model of tumor cell evolution
JO  - Comptes Rendus. Biologies
PY  - 2008
SP  - 823
EP  - 836
VL  - 331
IS  - 11
PB  - Elsevier
DO  - 10.1016/j.crvi.2008.08.010
LA  - en
ID  - CRBIOL_2008__331_11_823_0
ER  -

%0 Journal Article
%A Pablo Gómez-Mourelo
%A Eva Sánchez
%A Luis Casasús
%A Glenn F. Webb
%T A fully continuous individual-based model of tumor cell evolution
%J Comptes Rendus. Biologies
%D 2008
%P 823-836
%V 331
%N 11
%I Elsevier
%R 10.1016/j.crvi.2008.08.010
%G en
%F CRBIOL_2008__331_11_823_0

Pablo Gómez-Mourelo; Eva Sánchez; Luis Casasús; Glenn F. Webb. A fully continuous individual-based model of tumor cell evolution. Comptes Rendus. Biologies, Volume 331 (2008) no. 11, pp. 823-836. doi : 10.1016/j.crvi.2008.08.010. https://comptes-rendus.academie-sciences.fr/biologies/articles/10.1016/j.crvi.2008.08.010/

Version originale du texte intégral

1 Introduction

In this paper we construct and analyze a fully continuous individual-based model (IBM) for cancer invasion.There have been many recent studies of tumor growth in spatial settings, including [1–11]. In our model of spatial tumor growth our aim is two-fold: first we construct an IBM founded on stochastic events, which yields readily interpreted graphical visualizations of tumor growth. Second, we transform the IBM to a system of partial differential equations (PDE) for the continuum densities of the IBM elements (in the limit of large numbers of particles). Then we obtain the numerical solution of the PDE, and this solution is qualitatively compared to the IBM output, in order to validate the IBM approach.

Individual-based models in mathematical biology have been used extensively recently, because increased computational capability allows the possibility of emulating in virtual worlds the huge variability of individual behavior. The consideration of different kinds of individuals is natural in IBM technology, with heterogeneity of individuals embedded in object-oriented programming. A frequent assumption in differential equations models of population dynamics is the so-called mean-field description, which assumes that individuals are homogeneous in a spatial context. In such models the influence of environment and the interactions between individuals are averaged throughout space, and spatial dependence is simplified or omitted. For many biological processes, however, spatial heterogeneity of individuals is an essential feature in understanding population behavior, and IBM approaches are advantageous for inclusion of spatial heterogeneity (a basic reference for IBM is [12]).

Although the visual aspect of IBM simulations are appealing to non-experts, IBMs still lack a general theoretical framework to interpret their output, which are different every run. In order to bridge the gap between computer simulation and mathematical description of our tumor invasion IBM, we use a set of scaling techniques to derive large-scale descriptions from the individual behavior. General references for this approach are [13] and [14], and the seminal paper [15].

The biological aims of our study are to develop a mathematical description of tumor growth with emphasis on spatial invasion into surrounding extracellular matrix. The biological complexity of this process invokes multi-scale levels – sub-cellular, cellular, and tissue. Mathematical frameworks provide a medium for the biological understanding of tumor growth through integration of these levels. The larger biological aim of our study is to develop a reliable predictor of tumor growth amenable to computer simulations specific to individual patients.

In this paper we focus on a hybrid discrete-continuous model recently presented by A. Anderson in [16], which we use as starting point. Related developments and numerical treatments can be found in [17]. We refer the reader to these references for the biological considerations and features of the model. We build and study an improved model for the same phenomenon. The improvements to Anderson's model are:

(1) Our model does not have a rigid discrete lattice spatial structure, but instead provides a more biologically realistic continuum spatial setting for cell growth and proliferation;
(2) Our model includes a cell age structure connected to a PDE formulation to describe individual cell behavior, and in particular allows inclusion of the phases of the cell cycle correspondent to cell age;
(3) Our model is convertible to a completely PDE formulation (the model in [16] was only partially convertible) and is advantageous for computational simulations;
(4) Our model can be analyzed theoretically and the significance of the parameters can be quantified.

As in [16], we analyze conditions under which tumor growth exhibits fingering, and in order to study this phenomena, design simulations with different initial MM distributions: homogeneous and heterogeneous.

With the aim of validating the individual-based approach, we scale-up the IBM (Lagrangian framework) to obtain a density-based Eulerian procedure. We construct a set of PDEs equivalent to a large number of runs of the computer simulations. We obtain numerically the solutions to these PDEs and qualitatively compare the results of both approaches with special interest on fingering patterns. We have made extensive numerical simulations, which strongly indicate that both approaches show very similar results. A formal validation of the equivalence of the two approaches will be explored in future work.

The organization of the paper is as follows: Section 2 contains a brief description of the biological features which are taken into account for the construction of the model; in Sections 3 and 4 we describe the tumor invasion IBM and its equivalent formulation in terms of PDEs. Section 5 focuses on the numerical analysis of the PDE model. Section 6 concludes with a discussion of results and conclusions.

2 Biological background

Tumors develop as mutations occurring in key regulatory genes of cell proliferation. According to [16,18], one of the most important mutations involves the p53 gene, the so-called Guardian of Genome, which is found in mutated form in over half of all cancers. The p53 protein controls three cellular functions: proliferation, death and DNA repair, so that a loss of p53 function due to mutation allows for the propagation of damaged DNA to other cells.

The aim of this paper is to develop a model for the growth of a generic solid tumor for which we will assume a blood supply has been established. We are especially interested in the influence of cell adhesion (i.e., cell–cell and cell–matrix adhesion) in the tumor shape, as this is a key question in the invasive process. Molecules in charge of cell adhesion not only regulate cell migration, but proliferation and differentiation as well.

Following [16], we restrict the model to the behavior of four principal elements: tumor cells (both proliferating and quiescent, and classified by mutation phenotype), extracellular matrix macromolecules (MM), matrix degradative enzymes (MDE) (which break down the MM), and oxygen (which supplies essential nutrients from the surrounding vasculature). Tumor cells proliferate, enter quiescence, or die according to the concentration of nutrients provided by the capillary network, which is simulated by the diffusion of oxygen. Proliferating cells have the ability to migrate, but quiescent cells do not and there is experimental evidence that cells move towards high MM concentrations. Some regions may become more nutrient deficient, until eventually cells can no longer survive, degenerating into necrotic regions. Macromolecules disappear progressively with time as they are degraded by tumor cells, the diminishing rate being proportional to the MDE concentration. MM are non-diffusible since they are bound within the surrounding tissue. Matrix degradative enzymes undergo several processes, and such enzymes facilitate degradation of extracellular matrix allowing adhesion and attachment of tumor cells. They are produced by cells, diffuse, and decay naturally. Finally, oxygen diffuses, is produced by MM, is consumed by cells and decays naturally.

3 Description of the IBM for tumor invasion

In this section we describe the way in which the behavior of the four kinds of particles that we consider in our model have been taken into account to built the IBM. For simplicity we assume the tumor is contained in a 2-dimensional spatial region, although our methods are applicable to the 3-dimensional case.

3.1 Description of the particles

3.1.1 Tumor cells

The tumor consists of populations of proliferating and quiescent cells. Proliferating cells are capable of growing, dividing, entering quiescence and becoming necrotic. We consider the possibility for a cell to mutate only at division, and the probability of mutation is taken to be $t_{m}$ . We consider four types of mutations, so the number of actual phenotypes is five with the wild-type is also taken into account. This scheme of mutation is a simplified treatment (as in [16] and [17]). Our simplified framework of proliferating mutation classes is readily adapted to a greatly expanded number of phenotypes, but with greatly increased computational cost. Nonetheless, we believe this simplification to be representative of the effects of mutations.

We track progress through the cell cycle with cell age, and newly divided cells start with cell age 0. The scheme of mutation is linear, so a mutation results in a passage to the next mutation phase. The higher the mutation phenotype, the higher the cell's aggressiveness. This aggressiveness is due to the possibility of higher motility, lower proliferation age, higher oxygen uptake and greater MDE production.

Proliferating cells divide through a distribution of ages, and the probability for a cell to divide at a given age depends both on its age and the oxygen concentration nearby. Based on [17] we choose the proliferation rate in the IBM as:

pr (a, c) = {\begin{matrix} 10 c e^{- 10 (a - 1)} {(2 a - 1)}^{5} & if 0.5 < a < 6 \\ 0 & otherwise \end{matrix}

(1)

where

c : = c (x, y)

stands for oxygen concentration at position

(x, y)

, and a stands for cell age. The form of

pr (a, c)

is taken to be roughly proportional to oxygen concentration c, and roughly bell-shaped on age a. The intention was to develop an adequate mathematical function shape rather than to choose the exact real parameter values.

An absence of oxygen can also lead cells to the quiescence state, and the return to the proliferating state is only achieved if sufficient oxygen becomes available again. In the IBM, both changes are modelled with transition rates $t_{p q}$ (from proliferation to quiescence) and $t_{q p}$ (vice versa). As such, the transition rate $t_{p q}$ is assumed to roughly decay with oxygen concentration c, whilst $t_{q p}$ behaves the opposite way. With these ideas in mind, and using the formulas from [17]:

(2)

Experiments show that cells have random motility. A Brownian motion seems suitable to reproduce this movement and the uncertainty inherent in the process of modeling, as explained in [19]. In the IBM, at every time step, each cell jumps in a random direction. The jump length is distributed as a Gaussian with mean zero and standard deviation:

σ = D_{c} (i / 2 + 1) \sqrt{d t}

where

i = 0, \dots, 4

is the mutation type,

D_{c} (i / 2 + 1)

is the Brownian intensity for the corresponding mutation, and dt is the time step. In this way, the Brownian motion intensity is assumed to grow with mutation, proportionally to

(i / 2 + 1)

with the constant

D_{c} > 0

as the intensity of the wild-type

i = 0

Brownian motion. The term

\sqrt{d t}

is introduced due to technical reasons in order to reproduce a Wiener process (the interested reader can refer to [20]).

Cells are attracted by high MM concentrations, in a phenomenon called haptotaxis in which the intensity of this attraction is assumed to be proportional to the concentration gradient. Haptotaxis is, from the mathematical side, similar to chemotaxis (see [21]) and, in general, equivalent to any motion whose speed is proportional to the gradient of concentration of any substance, particle or individual type. In our model the intensity of the haptotaxis is described by a corresponding coefficient $ρ (i / 2 + 1)$ . As before, $i = 0, \dots, 4$ is the mutation type, and $ρ (i / 2 + 1)$ stands for the haptotaxis intensity for the corresponding mutation, assumed to grow with the mutation type. The constant $ρ > 0$ is the intensity of the wild-type haptotaxis. In the IBM, at every time step, each cell jumps a length

l (x, y) = ρ (i / 2 + 1) ‖ \nabla f (x, y) ‖ d t

in the direction of

\nabla f (x, y)

, where

f (x, y)

is the MM density at the position

(x, y)

and

‖ \cdot ‖

is the usual Euclidean norm. Finally, we assume the main cause of cell death is a low concentration of oxygen. Accordingly, the IBM computes at every time step, a death probability rate given by (see [17]):

μ^{*} (c) = \max (1.0 - c, 0.0)

(3)

3.1.2 Macromolecules (MM)

Macromolecules form the tissue surrounding tumor cells and are crucial in our model, since they induce the haptotaxis phenomenon resulting in tumor fingering. In particular, their initial spatial distribution determines the evolution of the tumor spatial pattern. Following [16,18], we will consider different types of initial MM spatial distributions, homogeneous and heterogeneous as explained below.

MMs undergo degradation produced by the MDE particles, and the intensity of their degradation is proportional to the product of concentrations of MDE and MM. This will produce a so-called reaction term in the equivalent PDE formulation (see [22,23]). To reproduce this phenomenon in the IBM, the probability for a MM particle to disappear at any time step is taken to be proportional to the concentration of MDE nearby. Another way to model this process is by considering the number of collisions per time unit between MDE and MM particles. We believe both ways are equivalent for sufficiently small time steps. More precisely, the probability for an MM particle at position $(x, y)$ to disappear at any small enough time step dt is $η m d t$ , where $η > 0$ is a constant related to the intensity of the reaction and $m : = m (x, y)$ is the concentration of MDE. MM particles are assumed bound, and thus have no motion at all, at least in comparison to tumor cells.

3.1.3 Matrix-degradative enzymes (MDE)

MDE have a natural tendency to spread throughout the tumor environment and experimentally exhibit a random motion. As for cells, a Brownian motion is implemented for MDE particles in our IBM. The deployment of this motion is exactly the same as with cells, the intensity coefficient now called $D_{m} > 0$ . MDE are generated by cells, and the production rate is assumed to be proportional to cell concentration. To properly emulate this in the IBM, at each time step the probability for an MDE particle to be created is proportional to the surrounding cell concentration, more precisely the probability is taken to be

P (new MDE particle) = [κ (i / 2 + 1) \sum_{i = 0}^{4} (n_{i p} + n_{i q})] d t

(4)

where

n_{i p}

n_{i q}

stand for the density of mutation type i, proliferating and quiescent cells, respectively and

κ > 0

is a proportionality constant. MDE also decays naturally at a fixed rate

σ^{*} > 0

, so a constant death probability

σ^{*} d t

is considered at any time step for every MDE particle.

3.1.4 Oxygen

Oxygen spreads throughout the environment, is produced at rates proportional to the MM concentration, is consumed by cell uptake, and decays naturally. All these processes are incorporated in the IBM in the previously described manner. The Brownian motion intensity is called $D_{o} > 0$ , cell uptake rate is $ω (i / 2 + 1)$ , with $ω > 0$ a constant, the constant decay rate is $ϕ > 0$ and the constant production rate is $ν > 0$ .

3.2 IBM computational scheme

We model the tissue as a $S : = [0, 10] \times [0, 10]$ square (we allow the particles to move outside this square, but we will focus on the inner region). Summing up, we model the evolution of 4 types of individuals: cells, MM, MDE and oxygen. We consider 5 phenotypes for cells: wild-type plus mutation types 1, 2, 3 and 4 and consider a linear mutation scheme, so each mutation makes a cell increasingly aggressive, according to the values given in Table 1.

Phenotype	Oxygen uptake	MDE production	Brownian coefficient	Haptotaxis coefficient
Wild type	ω	κ	D	ρ
Mutation 1	3ω/2	3κ/2	3D/2	3ρ/2
Mutation 2	2ω	2κ	2D	2ρ
Mutation 3	5ω/2	5κ/2	5D/2	5ρ/2
Mutation 4	3ω	3κ	3D	3ρ

As the cells evolve through mutations, they became increasingly aggressive, due to stronger motility, stronger oxygen uptake and stronger MDE production. We remark that our scheme is space continuous, avoiding limitations of discrete approaches. We think this is important in order to connect the discrete simulation to the PDE approach. See Fig. 1 for a schematic flow chart of the IBM algorithm.

Fig. 1
Summarized flowchart of the IBM algorithm.

3.2.1 Initialization stage

We are interested in the development of tumor fingers, which form as a consequence of the attraction cells experience for macromolecules. To discern the influence of initial MM spatial distributions on the development of tumor patterns, we carry out two different experiments, one with MM homogeneously distributed and the other heterogeneously distributed.

Cells.

In both experiments all cells are initially spatially distributed according to the density formula:

5 G (μ = (5.0, 5.0); σ = 1.0)

where

G (μ : = (μ_{x}, μ_{y}); σ)

is a 2D isotropic Gaussian (by isotropic we mean the standard deviation σ is equal in any direction, so the Gaussian is completely symmetric) centered at μ with standard deviation σ, whose density at

(x, y)

is given by:

\frac{1}{2 π σ^{2}} \exp {- \frac{1}{2 σ^{2}} [{(x - μ_{x})}^{2} + {(y - μ_{y})}^{2}]}

In this way, we model a pre-invasive stage, where cells are initially present in the neighbourhood of the domain center, with lower densities along any radius starting at the tumor center. Initially we assume that about 100 cells are present. Also, we assume that all cells are initially proliferating with an age distribution between 0 and 2 days.

Macromolecules.

As we previously indicated, we will carry out two experiments corresponding to homogeneous and heterogeneous MM distributions, respectively:

(A) Macromolecules homogeneously distributed along the $[0, 10] \times [0, 10]$ domain at $t = 0$ with a fixed density of 0.3;
(B) Macromolecules intentionally heterogeneously distributed over the same domain at $t = 0$ . This heterogeneous distribution has two maximum peaks in the upper-right and the lower-left zones, more exactly at points $(2.5, 2.5)$ and $(7.5, 7.5)$ . We choose this distribution to intensify the expected consequences of haptotaxis. The density is
$0.2 + 0.3 G ((2.5, 2.5); 1.0) + 0.3 G ((7.5, 7.5); 1.0)$

Matrix-degradative enzymes.

The density of MDE is taken to be half the proliferating cells density in both experiments.

Oxygen.

The density of oxygen is taken to be 10 times the MM density in both experiments.

Remark

Our parameters are chosen for mathematical illustration as well as for biological realism (as in [17]). Biological validation of these choices requires the experimental support of current efforts [24]. Simulations were run from $t = 0$ through $t = 20$ . In order to be reasonably sure that the qualitative results of our simulations were representative, we run the simulation 200 times. The results of these simulations were quite similar.

3.3 Routine to compute local densities and gradients

The local densities and gradients of oxygen, cells, MDE and MM are computed in several subroutines of the IBM code. In order to measure these quantities, we assume that particles have a limited perception of the surrounding neighbourhood, so they can measure, at least roughly, the local densities and gradients. We assume this perception to be limited to a rectangle whose axis lengths are $Δ_{x}, Δ_{y}$ units in x and y, respectively. This concept of limited perception is quite usual for IBMs in mathematical biology, [25], due to their description of local interactions.

For each particle k, which may be a cell, MM, MDE or an oxygen particle, let $P^{k} (t) : = (x^{k} (t), y^{k} (t))$ be its position in the 2D space. For example, let $F_{R}^{k}$ denote the number of MM particles, perceived by individual k on its right side, i.e., inside the rectangle

[x^{k} (t), x^{k} (t) + Δ_{x}] \times [y^{k} (t) - Δ_{y}, y^{k} (t) + Δ_{y}]

and by

F_{L}^{k}

the number of MM particles perceived by individual k on its left side, i.e., inside the rectangle

[x^{k} (t) - Δ_{x}, y^{k} (t)] \times [y^{k} (t) - Δ_{y}, y^{k} (t) + Δ_{y}]

We define the discrete approximations of the density and gradient respectively as:

D^{k} (t) : = \frac{F_{L}^{k} + F_{R}^{k}}{4 Δ_{x} Δ_{y}}

G^{k} (t) : = \frac{F_{R}^{k} - F_{L}^{k}}{Δ_{x}}

These approximations are assumed to represent the local densities and gradients for large numbers of particles, as happens in our simulation. These equations also hold for computing the densities of MDE, cells and oxygen, and analogous formulae are used to compute vertical gradients.

Several values of $Δ_{x}$ and $Δ_{y}$ were tested, but in the end we found the value $Δ_{x} = Δ_{y} = \frac{1}{8}$ to be reasonable, because this value is large enough to include a high number of particles inside, and low enough to be small as compared to the whole square. The choice is a trade-off between both circumstances. The high number of particles considered slows down the computation loop, because the previous measure of gradients and concentrations must consider each particle individually. The execution of our simulations took about 24 hours on a Pentium III processor. The results of the simulations can be seen in the figures. The interested reader can also check some videos (.avi format) of the simulations on http://dmaii.etsii.upm.es/~cells.

4 Continuous PDE formulation for the IBM

In order to describe the algorithm underlying the IBM in mathematical language, we will adapt the stochastic approach described in [26] and [15]: the evolution of each particle can be written as a stochastic differential equation for the motion together with Poisson processes for the birth and death. Then, under some assumptions, a passage to the limit when the number of particles tends to infinity leads to a PDE for the density of each kind of particles. We also refer to [14] and [27].

4.1 General description of the procedure

When dealing mathematically with the IBM, we must distinguish two kind of processes: motion and demography. Motion processes are only related to the spatial evolution of IBM particles, while demographic processes include births, deaths and transitions between different types of particles. Roughly speaking, motion processes will include transport and Brownian phenomena and will be described by stochastic differential equations (SDE), while demographic phenomena will be described by Poisson processes with certain intensity rates key to the model. For a detailed explanation on SDEs and Poisson processes we refer the reader to [28] and [29].

We describe the general procedure in two steps:

Step 1. First of all we consider a general IBM where only two-dimensional motion takes place and only one kind of particle appears. The position of the kth individual at time t is given by

P^{k} (t) : = (x^{k} (t), y^{k} (t)), k = 1, \dots, N

and we assume that the following SDE system is satisfied:

d x^{k} = a_{1} (x, y, t) d t + b_{1} (x, y, t) d W_{x}^{k} (t)

d y^{k} = a_{2} (x, y, t) d t + b_{2} (x, y, t) d W_{y}^{k} (t)

where

a_{i}

b_{i}

i = 1, 2

, are given functions, and

W_{x}^{k} (t), W_{y}^{k} (t)

are the components of independent 2D Brownian motions

W^{k}

k = 1, \dots, N

We define the empirical measure of the system by:

P_{N} (t) : = \frac{1}{N} \sum_{k = 1}^{N} δ_{P^{k} (t)}

where

δ_{a}

is the Dirac's delta at point a. The main point consists of assuming that as

N \to \infty

, the empirical process

P_{N} (t)

tends to a deterministic process with density a

ρ (x, t)

, i.e.,

\lim_{N \to \infty} 〈 P_{N} (t), f (\cdot, t) 〉 = \int_{R^{2}} f (x, t) ρ (x, t) d x

for any function test f (see [13,15,27] for the details). Then, it can be shown that the density

ρ (x, t)

satisfies (in a weak sense) the so-called Kolmogorov or Fokker–Planck PDE:

\frac{\partial ρ}{\partial t} = \frac{1}{2} \frac{\partial^{2}}{\partial x^{2}} (b_{1}^{2} ρ) + \frac{1}{2} \frac{\partial^{2}}{\partial y^{2}} (b_{2}^{2} ρ) - \frac{\partial}{\partial x} (a_{1} ρ) - \frac{\partial}{\partial y} (a_{2} ρ)

Step 2. Now we consider R different classes of particles and we denote by $N_{i}$ the number of particles in class i, $i = 1, \dots, R$ . The equations of movement for each particle are as above. In addition, we assume that each particle of the class i can

• die, at a rate given by $d_{i} (P, t,)$ , where P denotes the position of the particle;
• give birth to a new particle of the class q, at a rate given by $b_{i q} (P, t)$ ;
• pass to a different class particle, at a rate given by $t_{i q} (P, t)$ , for $i \neq q$ , with $t_{i i} = 0$ .

The rates

d_{i}, b_{i q}, t_{i q}

represent the intensity of Poisson random processes, and can be understood as the probabilities of death, birth and transition in a small time step.

We define an empirical measure

P_{i, N_{i}} (t) : = \frac{1}{N_{i}} \sum_{k = 1}^{N_{i}} δ_{P_{i}^{k} (t)}

for each class i of particles,

i = 1, \dots, R

. As before, we assume that, as the numbers of particles of each family

N_{i}

tends to infinity, the empirical process

P_{i N_{i}} (t)

tends to a deterministic process with density

ρ_{i} (x, t)

. Then it can be shown that (see [15]) the densities

ρ_{i}

satisfy the system of PDEs:

\frac{\partial ρ_{i}}{\partial t} = \frac{1}{2} \frac{\partial^{2}}{\partial x^{2}} (b_{1 i}^{2} ρ_{i}) + \frac{1}{2} \frac{\partial^{2}}{\partial y^{2}} (b_{2 i}^{2} ρ_{i}) - \frac{\partial}{\partial x} (a_{1 i} ρ_{i}) - \frac{\partial}{\partial y} (a_{2 i} ρ_{i}) + (- \sum_{q = 1, q \neq i}^{R} t_{q i} + \sum_{q = 1, q \neq i}^{R} t_{i q} + \sum_{q = 1}^{R} b_{i q} - \sum_{q = 1}^{R} d_{i}) ρ_{i}

i = 1, \dots, R

4.2 SDEs of the IBM

We start with a few words about what we mean by a Poisson process with variable intensity [30]. Let $N (t)$ be an integer-valued stochastic process defined for every time $t ⩾ 0$ . By saying $N (t)$ is distributed as a Poisson process with variable intensity $λ (t)$ we mean its probability density is given by

P (N (t) = j) : = \frac{1}{j!} (\int_{0}^{t} λ (u) d u) e^{- \int_{0}^{t} λ (u) d u} j = 0, 1, \dots

In our model we are only concerned by a much more restricted case, as long as only one jump is considered, its probability given by

P (N (t) = 1) = (\int_{0}^{t} λ (u) d u) e^{- \int_{0}^{t} λ (u) d u}

We do not consider further jumps (

j = 2, 3, \dots

), because we use Poisson processes as stopping-time alarm-clocks, which determine the instants for discontinuous changes (either births, deaths or transitions) to occur. For example, let us consider a certain cell born at time

t_{0}

, whose death intensity is given by

d (t)

for each time

t > t_{0}

. Then, the probability for this cell to suffer death in the interval

[t_{0}, t]

is given by

P (death) = (\int_{t_{0}}^{t} d (σ) d σ) e^{- \int_{t_{0}}^{t} d (σ) d σ}

Thus, the Poisson process is only used as an alarm-clock: when the clock goes off, that is, when the Poisson process take the value

N (t) = 1

, the cell dies. Birth and transition processes are modelled in exactly the same way.

Remark

Combining the motion and the demographic processes:

• Every SDE for motion of a particle is valid during the lifetime of that particle and only during its lifetime;
• Each lifetime of a particle is a time interval whose start and end are determined as consequence of the relevant birth, death and transition Poisson processes.

We assume, for sake of simplicity, that each integer number k is used only once throughout the simulation; this means any newborn particle will be assigned an unused index k, which will not be used again after its death. As such, each particle k will have a unique lifetime interval. The convenience of this assumption is notation-related and will be apparent in subsequent paragraphs. In the next section we will explain the SDEs which describe the evolution of the four types of particles considered in the IBM, in terms of the corresponding Poisson processes.

4.2.1 Cells

• Motion. The position and age of the kth cell at time t are represented by $P^{k} (t) = (x^{k} (t), y^{k} (t), a^{k} (t))$ where $x^{k}$ is the horizontal position, $y^{k}$ is the vertical position, and $a^{k}$ is the age. According to the explanation given before in Section 3.1.1, we have:
$d P^{k} (t) = (ρ (1 + \frac{i}{2}) \frac{\partial f}{\partial x}, ρ (1 + \frac{i}{2}) \frac{\partial f}{\partial y}, 1) d t + D_{c} (1 + \frac{i}{2}) d W^{k}$
where $W^{k}$ are independent Brownian motions for each cell k. The components $ρ (i / 2 + 1) (\partial f / \partial x), ρ (i / 2 + 1) (\partial f / \partial y)$ represent the haptotaxis, i.e., a displacement proportional to the gradient of MM concentration, while the third component is the aging speed, equal to one per time unit. $D_{c} (i / 2 + 1)$ is a term storing the intensity of the Brownian motion, which increases with the mutation type.
• Demography. The model includes the following Poisson processes:
- – Transition from proliferating state to quiescent state, and from quiescent state to proliferating state, at rates $t_{p q}$ and $t_{q p}$ , respectively, given by (2);
- – Proliferation with rate given by (1);
- – Death, with rate given by (3).

4.2.2 MM

• Motion. MM particles experience no motion, so their position remains fixed throughout the simulation and neither transport nor Brownian effects appear.
• Demography. Only one death Poisson process appears: degradation at a rate given by ηm (see Section 3.1.2).

4.2.3 MDE

• Motion. MDE particles experience Brownian motion, so the motion of the kth particle is given by
$d P^{k} (t) : = D_{m} d W^{k}$
where $D_{m} > 0$ stands for the diffusion coefficient and $W^{k}$ are independent Brownian motions for each MDE particle k.
• Demography. We consider the following Poisson processes:
- – Production: MDE particles are produced by cells at a rate given by (4);
- – Decay: MDE particles disappear at a constant rate $σ^{*} > 0$ .

4.2.4 Oxygen

• Motion. Oxygen particles experience Brownian motion given by
$d P^{k} (t) : = D_{o} d W^{k}$
• Demography. We consider the following Poisson processes (see Section 3.1.4):
- – Production: oxygen particles are produced by MM particles at constant rate $ν > 0$ .
- – Uptake: oxygen particles are uptaken by cells of type i at a rate $ω (i / 2 + 1)$ .
- – Decay: oxygen particles disappear naturally at a constant rate $ϕ > 0$ .

4.3 Obtaining the PDE system for the IBM

As explained in Section 4.1, we develop the following limiting PDE system for the IBM:

• Cells.
- – Proliferating state:
  $\frac{\partial n_{i p}}{\partial t} + \frac{\partial n_{i p}}{\partial a} = D_{c} (\frac{i}{2} + 1) Δ n_{i p} - ρ (\frac{i}{2} + 1) \nabla (n_{i p} \nabla f) - μ^{*} (c) n_{i p} - p r (a, c) n_{i p} + t_{q p} (c) n_{i q}$
- – Quiescent state:
  $\frac{\partial n_{i q}}{\partial t} + \frac{\partial n_{i q}}{\partial a} = - μ^{*} (c) n_{i q} + p r (a, c) t_{p q} (c) n_{i p} - t_{q p} (c) n_{i q}$
- – Boundary conditions:
  $n_{i p} (a = 0) = 2 \int_{0}^{a_{M}} n_{i p} p r (a, c) (1 - t_{p q} (c)) (1 - t_{m}) d a + 2 \int_{0}^{a_{M}} n_{i - 1, p} p r (a, c) (1 - t_{p q} (c)) t_{m} d a$
  $i = 0, 1, 2, 3, 4$ . The boundary conditions appear naturally from the fact that a proliferating cell gives birth to two daughters, whose mutation type can be either the same (in case no recent mutation has occurred) or the next sequential type (in case the last division led to a new mutation). In our simulation we have taken $t_{m} = 0.1$ as in [16].
• MM.
$\frac{\partial f}{\partial t} = - η m f$
• MDE.
$\frac{\partial m}{\partial t} = D_{m} Δ m + κ \sum_{i} (n_{i, p} + n_{i, q}) - σ m$
• Oxygen.
$\frac{\partial c}{\partial t} = D_{c} Δ c + ν f - ω \sum_{i} (n_{i, p} + n_{i, q}) - ϕ c$

5 Numerical analysis of the PDE

To analyze this model computationally we used a finite difference discretization in space that for the range of parameters of this work, allows large time steps. For the 2D version, we use $Δ x = Δ y = 1 / 8$ and $Δ t = 0.01$ , in the domain $S = [0, 10] \times [0, 10]$ .

The results of the simulations can be seen in the figures. The interested reader can also check some videos (.avi format) of the numerical resolution on http://dmaii.etsii.upm.es/~cells.

A fully implicit scheme is far too expensive for practical calculations. We succeeded in using an improved predictor-corrector scheme with a double purpose: to test the convergence of our straightforward explicit discretization in time and space and to develop a robust general method for accurately solving the system given in Section 4.3. Discretization of the age variable does not suffer the difficulties mentioned in [17] because our methods are not semi-discrete.

In order to obtain a more complete picture of the process, we have started extending the two-dimensional work to three space dimensions. The experiments reveal that the general morphology is similar to that of the 2D version. As a significative (though preliminary) example, Fig. 8 shows a 3D shell of type-4 cells, with an inner necrotic core represented as an area void of cells. In the 3D case, we have been able to efficiently compute different tumors with $Δ x = Δ y = Δ z = 1 / 12$ and $Δ t = 0.005$ . Moreover, the finite-difference procedures are well-suited for parallel implementation.

Fig. 8
A picture of a 3D simulation of the IBM. The inner core is quiescent and necrotic.

We are presently developing a splitting procedure and different approximations of the div-grad operator to allow the accurate modeling of morphologic instabilities, and in particular, more intricate and irregular tumor boundary evolution and satellite tumors around the main tumor mass.

6 Conclusion

In this work we have studied the invasion of tumor cells into surrounding tissue by using two approaches: one, a computational individual-based perspective; the other, a partial differential equations perspective reproducing the same phenomena. Two computer experiments were carried out for both approaches: (a) one starting with an initially homogeneous MM density (homogeneous case); and the other (b) starting with an initially heterogeneous MM density (heterogeneous case). We considered four types of cell phenotypes evolved by sequential mutation from the wild-type. The higher the mutation type in the sequence, the higher the cell's aggressiveness in the invasion.

Simulations show the expansion of increasingly aggressive cell phenotypes, as can be observed in Figs. 2 and 3. In the homogeneous MM case, cells evolve in a ring-shaped fashion (Fig. 2). We note that quiescent and necrotic cells are not drawn in the figures and the inner zone of the ring contains only quiescent and necrotic cells. Proliferating cells advance outwards, attracted by MM. In the heterogeneous case, cells do not form a ring, but advance toward higher concentration of MM (Fig. 3).

Fig. 2
Snapshots of the evolution of the IBM (homogeneous MM initial distribution). For clarity of representation, quiescent cells are not represented. The color indicates the mutation (blue = wild-type, cyan = type 1, green = type 2, orange = type 3, red = type 4).

Fig. 3
Snapshots of the evolution of the IBM (heterogeneous MM initial distribution). As in Fig. 2, for the sake of clarity, only proliferating cells are represented. The color indicates the mutation with the same code (blue = wild-type, cyan = type-1, green = type-2, orange = type-3, red = type-4).

An important biological feature of tumor invasion is fingering patterns (Fig. 7). Our simulations systematically reproduce such patterns, both in the homogeneous case (Fig. 2) and the heterogeneous case (Fig. 3). Fingering was less evident in the homogeneous case, nevertheless it appears even under fairly uniform densities of MM. We conclude that very slight oscillations of MM density may be enough to induce finger patterns. We hypothesize that in vivo MM concentrations may present slight density heterogeneities, which are eventually responsible for fingering patterns.

Fig. 7
Finger patterns appear due to MM density fluctuations, not only in the heterogeneous case (right), but in the homogeneous case (left) as well.

Our simulations are the results of spatially explicit, fully continuous, individual based computation, whose building blocks are tokens for cells, MM, MDE elements and oxygen particles (Fig. 4). We derived an equivalent (in a mathematical sense) PDE model to describe the same behavior. The PDE numerical simulations match remarkably well those of the IBM. The homogeneous case develops a crown-shaped cell density (Fig. 5), with inner zone density close to zero, as corresponding to the quiescent and necrotic zone. In the heterogeneous case (Fig. 6), two main directions of invasion appear, as a consequence of the initial spatial configuration of the MM density. The concordance of the IBM and PDE approaches is evident in Figs. 2 and 5 and Figs. 3 and 6.

Fig. 4
The three boxes represent the oxygen, MM, and MDE density profile, respectively, at a representative stage of the simulation.

Fig. 5
Snapshots of the type-4 cells density evolution, as provided by the numerical solution of the PDE (case of homogeneous MM initial distribution). Only type-4 cells are taken into account in order to make the figure interpretation easier.

Fig. 6
Snapshots of the type-4 cells density evolution, as provided by the numerical solution of the PDE (case of heterogeneous MM initial distribution). As in Fig. 5, only the most aggressive (type-4) cell density is represented.

IBM computational simulations are accessible and intuitive. They do not require deep mathematical knowledge and can be very helpful for non-experts. However, high computation costs appear when dealing with many particles; in particular, the computation of interactions between N particles increases (roughly) with $N^{α}$ , with $α ⩾ 2$ making such simulations hardly feasible for very large values of N. As a consequence, it is important to acquire a correspondent PDE approach that reliably reproduces equivalent – at least in the mathematical sense of Section 4 – results. We believe IBM modeling to be valuable, but its practical limitations must be known and appreciated.

From a paradigmatic point of view, it seems to us more natural to model large tumor cell populations as continuous mathematical densities rather than as frameworks of discrete units, even taking into account the variable nature of individual cells. Furthermore, we believe the continuous approach in an IBM overcomes some limitations of the discrete one, bridging the gap between many-particle systems and differential equations-based models. We find it particularly relevant and fruitful to establish comparisons between computational and PDE approaches. Not only can such comparisons give insight into complex population behavior in general (not only tumors), but can as well be tools to test models. The scheme in Section 4 to establish the comparison between IBMs and PDEs remains valid in general cases, and can be used by computational modelers seeking mathematical validation.

As for future work, we plan to extend our work to three dimensions. In fact, we have already run tentative 3D simulations, which showed similar results to the 2D case. Fig. 8 is a snapshot of a 3D simulation displaying a quasi-spherical tumor with a necrotic inner core. Fingering patterns were also observed in our 3D simulations, which are still under development.

Bibliographie

[1] Cancer Modelling and Simulation (L. Preziosi, ed.), Chapman & Hall/CRC, Boca Raton, FL, 2003

[2] R.P. Araujo; L.S. Mcelwain A history of the study of solid tumour growth: The contribution of mathematical modelling, Bull. Math. Biol., Volume 66 (2004), pp. 1039-1091

[3] A.R.A. Anderson; A. Weaver; P. Cummings; V. Quaranta Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment, Cell, Volume 127 (2006), pp. 905-915

[4] C. Walker; G.F. Webb Global existence of classical solutions for a haptotaxis model, SIAM J. Math. Anal., Volume 38 (2007) no. 5, pp. 1694-1713

[5] I. Ramis-Conde; M.A.J. Chaplain; A.R.A. Anderson Mathematical modelling of cancer cell invasion of tissue, Math. Comp. Mod., Volume 47 (2008), pp. 533-545

[6] M.A.J. Chaplain Modelling aspects of cancer growth: Insight from mathematical and numerical analysis and computational simulation (V. Capasso; M. Lachowicz, eds.), Multiscale Problems in the Life Sciences: From Microscopic to Macroscopic, Lecture Notes in Mathematics, Springer, Heidelberg, 2008

[7] N. Bellomo; N. Li; P. Maini On the foundations of cancer modelling selected topics, speculations, and perspectives, Math. Models Methods Appl. Sci., Volume 18 (2008), pp. 593-646

[8] A. Belloquid; M. Delitala Mathematical Modeling of Complex Biological Systems: A Kinetic Theory Approach, Birkhäuser, Berlin, 2006

[9] J. Dyson; R. Villella-Bressan; E. Sanchez; G.F. Webb An age and spatially structured model of tumor invasion with haptotaxis, Discr. Cont. Dyn. Sys. B, Volume 8 (2007) no. 1, pp. 45-60

[10] J. Dyson; R. Villella-Bressan; G.F. Webb An age and spatially structured model of tumor invasion with haptotaxis II, Math. Pop. Studies, Volume 15 (2008), pp. 1-23

[11] J. Dyson; R. Villella-Bressan; G.F. Webb A spatial model of tumor growth with cell age, cell size, and mutation of cell phenotypes, Math. Mod. Nat. Phen., Volume 2 (2008) no. 3, pp. 69-100

[12] Individual-Based Models and Approaches in Ecology: Populations, Communities and Ecosystems (D.L. DeAngelis; L.J. Gross, eds.), Chapman and Hall, New York, 1992

[13] V. Capasso, Stochastic differential equations and systems of interacting particles, in: 1st ESMTB Summer School Lecture Notes, 2000, pp. 133–149

[14] A. Stevens Derivation of chemotaxis-equations as limit dynamics of moderately interacting stochastic many particle systems, SIAM J. Appl. Math., Volume 61 (2000) no. 1, pp. 183-212

[15] K. Oelschlaeger On the derivation of reaction-diffusion equations as limit dynamics of systems of moderately interacting stochastic many particle processes, Probab. Theory and Related Fields, Volume 82 (1989), pp. 565-586

[16] A.R.A. Anderson A hybrid mathematical model of solid tumor invasion: The importance of cell adhesion, Math. Med. Biol., Volume 22 (2005), pp. 163-186

[17] B.P. Ayati; G.F. Webb; A.R.A. Anderson Computational methods and results for structured multiscale models of tumor invasion, Multiscale Model. Simul., Volume 5 (2006) no. 1, pp. 1-20

[18] A.R.A. Anderson; M.A.J. Chaplain Continuous and discrete mathematical models of tumour-induced angiogenesis, Bull. Math. Biol., Volume 60 (1998), pp. 857-899

[19] A. Okubo Diffusion and Ecological Problems: Mathematical Models, Springer-Verlag, New York, 1980

[20] I.I. Gikhman; A.V. Skorokhod Introduction to the Theory of Random Processes, Dover, New York, 1969

[21] H. Othmer; A. Stevens Aggregation, Blowup and Collapse: The ABC's of generalized taxis, SIAM J. Applied Math., Volume 57 (1997), pp. 1044-1081

[22] T. Czárán Spatiotemporal Models of Population and Community Dynamics, Chapman and Hall, London, 1998

[23] R.S. Cantrell; C. Cosner Spatial Ecology via Reaction-Diffusion Equations, John Wiley and Sons, UK, 2003

[24] A.R.A. Anderson; V. Quaranta Integrative mathematical oncology, Nature Rev. Cancer, Volume 8 (2008), pp. 227-234

[25] M. Adioui; J.P. Treuil; O. Arino Alignment in a fish school: a mixed Lagrangian–Eulerian approach, Ecol. Model., Volume 167 (2003), pp. 19-32

[26] D. Morale; V. Capasso; K. Oelschlaeger An interacting particle system modelling aggregation behavior: from individuals to populations, J. Math. Biol., Volume 50 (2005) no. 1, pp. 49-66

[27] P. Gomez-Mourelo From individual-based models to partial differential equations. An application to the upstream movement of elvers, Ecol. Model., Volume 188 (2005) no. 1, pp. 93-111

[28] B. Oksendal Stochastic Differential Equations, Springer-Verlag, Heidelberg, 2003

[29] L. Arnold Stochastic Differential Equations: Theory and Applications, Wiley, New York, 1974

[30] G.R. Grimmett; D.R. Stirzaker Probability and Random Processes, Clarendon Press, Oxford, 1992