## 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:

$$\mathit{pr}(a,c)=\{\begin{array}{cc}10c{e}^{-10(a-1)}{(2a-1)}^{5}\hfill & \text{if}\phantom{\rule{0.25em}{0ex}}0.5<a<6\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$ | (1) |

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}_{pq}$ (from proliferation to quiescence) and ${t}_{qp}$ (vice versa). As such, the transition rate ${t}_{pq}$ is assumed to roughly decay with oxygen concentration c, whilst ${t}_{qp}$ 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:

$$\sigma ={D}_{c}(i/2+1)\sqrt{dt}$$ |

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 $\rho (i/2+1)$. As before, $i=0,\dots ,4$ is the mutation type, and $\rho (i/2+1)$ stands for the haptotaxis intensity for the corresponding mutation, assumed to grow with the mutation type. The constant $\rho >0$ is the intensity of the wild-type haptotaxis. In the IBM, at every time step, each cell jumps a length

$$l(x,y)=\rho (i/2+1)\Vert \mathrm{\nabla}f(x,y)\Vert dt$$ |

$${\mu}^{\ast}(c)=\mathrm{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 $\eta m\phantom{\rule{0.2em}{0ex}}dt$, where $\eta >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(\text{new MDE particle})=[\kappa (i/2+1)\sum _{i=0}^{4}({n}_{ip}+{n}_{iq})]\phantom{\rule{0.2em}{0ex}}dt$$ | (4) |

#### 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 $\omega (i/2+1)$, with $\omega >0$ a constant, the constant decay rate is $\varphi >0$ and the constant production rate is $\nu >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.

**Table 1**

The IBM computational scheme

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.

#### 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:

$$5G(\mu =(5.0,\phantom{\rule{0.25em}{0ex}}5.0);\phantom{\rule{0.25em}{0ex}}\sigma =1.0)$$ |

$$\frac{1}{2\pi {\sigma}^{2}}\mathrm{exp}\{-\frac{1}{2{\sigma}^{2}}[{(x-{\mu}_{x})}^{2}+{(y-{\mu}_{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.3G((2.5,\phantom{\rule{0.25em}{0ex}}2.5);\phantom{\rule{0.25em}{0ex}}1.0)+0.3G((7.5,\phantom{\rule{0.25em}{0ex}}7.5);\phantom{\rule{0.25em}{0ex}}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.

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 ${\mathrm{\Delta}}_{x},{\mathrm{\Delta}}_{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)+{\mathrm{\Delta}}_{x}]\times [{y}^{k}(t)-{\mathrm{\Delta}}_{y},{y}^{k}(t)+{\mathrm{\Delta}}_{y}]$$ |

$$[{x}^{k}(t)-{\mathrm{\Delta}}_{x},{y}^{k}(t)]\times [{y}^{k}(t)-{\mathrm{\Delta}}_{y},{y}^{k}(t)+{\mathrm{\Delta}}_{y}]$$ |

$${D}^{k}(t):=\frac{{F}_{L}^{k}+{F}_{R}^{k}}{4{\mathrm{\Delta}}_{x}{\mathrm{\Delta}}_{y}}$$ |

$${G}^{k}(t):=\frac{{F}_{R}^{k}-{F}_{L}^{k}}{{\mathrm{\Delta}}_{x}}$$ |

Several values of ${\mathrm{\Delta}}_{x}$ and ${\mathrm{\Delta}}_{y}$ were tested, but in the end we found the value ${\mathrm{\Delta}}_{x}={\mathrm{\Delta}}_{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)),\phantom{\rule{1em}{0ex}}k=1,\dots ,N$$ |

$$d{x}^{k}={a}_{1}(x,y,t)\phantom{\rule{0.2em}{0ex}}dt+{b}_{1}(x,y,t)\phantom{\rule{0.2em}{0ex}}d{W}_{x}^{k}(t)$$ |

$$d{y}^{k}={a}_{2}(x,y,t)\phantom{\rule{0.2em}{0ex}}dt+{b}_{2}(x,y,t)\phantom{\rule{0.2em}{0ex}}d{W}_{y}^{k}(t)$$ |

We define the empirical measure of the system by:

$${P}_{N}(t):=\frac{1}{N}\sum _{k=1}^{N}{\delta}_{{P}^{k}(t)}$$ |

$$\underset{N\to \infty}{\mathrm{lim}}\langle {P}_{N}(t),f(\cdot ,t)\rangle =\underset{{\mathbb{R}}^{2}}{\int}f(x,t)\rho (x,t)\phantom{\rule{0.2em}{0ex}}dx$$ |

$$\frac{\partial \rho}{\partial t}=\frac{1}{2}\frac{{\partial}^{2}}{\partial {x}^{2}}\left({b}_{1}^{2}\rho \right)+\frac{1}{2}\frac{{\partial}^{2}}{\partial {y}^{2}}\left({b}_{2}^{2}\rho \right)-\frac{\partial}{\partial x}({a}_{1}\rho )-\frac{\partial}{\partial y}({a}_{2}\rho )$$ |

**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}_{iq}(P,t)$;
- • pass to a different class particle, at a rate given by ${t}_{iq}(P,t)$, for $i\ne q$, with ${t}_{ii}=0$.

We define an empirical measure

$${P}_{i,{N}_{i}}(t):=\frac{1}{{N}_{i}}\sum _{k=1}^{{N}_{i}}{\delta}_{{P}_{i}^{k}(t)}$$ |

$$\frac{\partial {\rho}_{i}}{\partial t}=\frac{1}{2}\frac{{\partial}^{2}}{\partial {x}^{2}}\left({b}_{1i}^{2}{\rho}_{i}\right)+\frac{1}{2}\frac{{\partial}^{2}}{\partial {y}^{2}}\left({b}_{2i}^{2}{\rho}_{i}\right)-\frac{\partial}{\partial x}({a}_{1i}{\rho}_{i})-\frac{\partial}{\partial y}({a}_{2i}{\rho}_{i})+(-\sum _{q=1,\phantom{\rule{0.2em}{0ex}}q\ne i}^{R}{t}_{qi}+\sum _{q=1,\phantom{\rule{0.2em}{0ex}}q\ne i}^{R}{t}_{iq}+\sum _{q=1}^{R}{b}_{iq}-\sum _{q=1}^{R}{d}_{i}){\rho}_{i}$$ |

### 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\u2a7e0$. By saying $N(t)$ is distributed as a Poisson process with variable intensity $\lambda (t)$ we mean its probability density is given by

$$P(N(t)=j):=\frac{1}{j!}(\underset{0}{\overset{t}{\int}}\lambda (u)\phantom{\rule{0.2em}{0ex}}du){e}^{-{\int}_{0}^{t}\lambda (u)\phantom{\rule{0.2em}{0ex}}du}\phantom{\rule{1em}{0ex}}j=0,1,\dots $$ |

$$P(N(t)=1)=(\underset{0}{\overset{t}{\int}}\lambda (u)\phantom{\rule{0.2em}{0ex}}du){e}^{-{\int}_{0}^{t}\lambda (u)\phantom{\rule{0.2em}{0ex}}du}$$ |

$$P(\text{death})=(\underset{{t}_{0}}{\overset{t}{\int}}d(\sigma )\phantom{\rule{0.2em}{0ex}}d\sigma ){e}^{-{\int}_{{t}_{0}}^{t}d(\sigma )\phantom{\rule{0.2em}{0ex}}d\sigma}$$ |

##### 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:

where ${W}^{k}$ are independent Brownian motions for each cell k. The components $\rho (i/2+1)(\partial f/\partial x),\rho (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.$$d{P}^{k}(t)=(\rho (1+\frac{i}{2})\frac{\partial f}{\partial x},\rho (1+\frac{i}{2})\frac{\partial f}{\partial y},1)\phantom{\rule{0.2em}{0ex}}dt+{D}_{c}(1+\frac{i}{2})\phantom{\rule{0.2em}{0ex}}d{W}^{k}$$ - •
**Demography.**The model includes the following Poisson processes:

#### 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

where ${D}_{m}>0$ stands for the diffusion coefficient and ${W}^{k}$ are independent Brownian motions for each MDE particle k.$$d{P}^{k}(t):={D}_{m}\phantom{\rule{0.2em}{0ex}}d{W}^{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 ${\sigma}^{\ast}>0$.

#### 4.2.4 Oxygen

- •
**Motion.**Oxygen particles experience Brownian motion given by$$d{P}^{k}(t):={D}_{o}\phantom{\rule{0.2em}{0ex}}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 $\nu >0$.
- – Uptake: oxygen particles are uptaken by cells of type i at a rate $\omega (i/2+1)$.
- – Decay: oxygen particles disappear naturally at a constant rate $\varphi >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}_{ip}}{\partial t}+\frac{\partial {n}_{ip}}{\partial a}={D}_{c}(\frac{i}{2}+1)\mathrm{\Delta}{n}_{ip}-\rho (\frac{i}{2}+1)\mathrm{\nabla}({n}_{ip}\mathrm{\nabla}f)-{\mu}^{\ast}(c){n}_{ip}-pr(a,c){n}_{ip}+{t}_{qp}(c){n}_{iq}$$ - –
**Quiescent state**:$$\frac{\partial {n}_{iq}}{\partial t}+\frac{\partial {n}_{iq}}{\partial a}=-{\mu}^{\ast}(c){n}_{iq}+pr(a,c){t}_{pq}(c){n}_{ip}-{t}_{qp}(c){n}_{iq}$$ - –
**Boundary conditions**:

$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].$${n}_{ip}(a=0)=2\underset{0}{\overset{{a}_{M}}{\int}}{n}_{ip}pr(a,c)(1-{t}_{pq}(c))(1-{t}_{m})\phantom{\rule{0.2em}{0ex}}da+2\underset{0}{\overset{{a}_{M}}{\int}}{n}_{i-1,p}pr(a,c)(1-{t}_{pq}(c)){t}_{m}\phantom{\rule{0.2em}{0ex}}da$$

- –
- •
**MM.**$$\frac{\partial f}{\partial t}=-\eta mf$$ - •
**MDE.**$$\frac{\partial m}{\partial t}={D}_{m}\mathrm{\Delta}m+\kappa \sum _{i}({n}_{i,p}+{n}_{i,q})-\sigma m$$ - •
**Oxygen.**$$\frac{\partial c}{\partial t}={D}_{c}\mathrm{\Delta}c+\nu f-\omega \sum _{i}({n}_{i,p}+{n}_{i,q})-\varphi 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 $\mathrm{\Delta}x=\mathrm{\Delta}y=1/8$ and $\mathrm{\Delta}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 $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z=1/12$ and $\mathrm{\Delta}t=0.005$. Moreover, the finite-difference procedures are well-suited for parallel implementation.

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).

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.

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.

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}^{\alpha}$, with $\alpha \u2a7e2$ 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.