framework for cell-population dynamics

UPTEC F 18036

Examensarbete 30 hp Juni 2018

Computational modeling of avascular tumours using a hybrid on-lattice

framework for cell-population dynamics

Lina Viklund

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Computational modeling of avascular tumours using a hybrid on-lattice framework for cell-population

dynamics

Lina Viklund

This master thesis consists of integrating an advanced model for avascular tumour growth into a multiscale computational framework for cell-population dynamics. In order to demonstrate this framework, the authors of

the framework developed a simple model for avascular tumour growth, governed by the Laplace equation and a set of rules relating to the oxygen concentration, that they integrated with the framework.

The realisation of this model closely followed previously established growth patterns of tumour cells up until the last phase, where predictions by some authors say the volume of the tumour should plateau. However, an

additional phase was observed where asymmetric protrusions appear on the tumour's surface, causing the volume to continue to increase. The authors speculate that this could be due to their model lacking the inclusion

of cell-cell adhesion, which restrains the migration of cells.

The advanced tumour model implemented in this thesis include several additional mechanisms that the simple model lacks, and thus the driving research question for the thesis is whether integrating this advanced model

into the framework will produce results that do not exhibit protrusions. Results suggest that despite this more advanced model, the outcome will still be an asymmetric tumour with unrestricted growth.

ISSN: 1401-5757, UPTEC F18036

Examinator: Tomas Nyberg

Ämnesgranskare: Per Lötstedt

Handledare: Stefan Engblom

Populärvetenskaplig sammanfattning

Cellbaserade beräkningmodeller är ett verktyg där beräkningar används för att undersöka hur grup- per av celler beter sig över tid. Då det är enkelt att upprepa experiment med hjälp av program kan dessa beräkningsmodeller användas för att på ett snabbt och smidigt sätt testa olika ideér kring cellers beteende. Den här sortens modeller har inom forskning använts för att modellera många olika biologiska processer, exempelvis hur blodkärl bildas, sår läker, embryon utvecklas eller hur cancertumörer växer.

I detta examensarbete implementeras en sådan cellbaserad modell som beskriver tillväxt av tumörer och integreras i ett beräkningsramverk för celldynamik, ett skalprogram som beräknar cellers rörelse utifrån en fix mekanik men som också kan modifieras utifrån vilken sorts cellpopulation man vill modellera. Den sorts tumör som modellerats har varit en tumör i dess allra första tillväxtstadium.

I detta stadium finns det inga blodkärl till tumören som kan förse den med syre, utan tumören är beroende av syre som tränger in i den från omkringliggande vävnad.

Det finns två karaktäristiska egenskaper för det första tillväxtstadiet hos tumören. För det första uppkommer vanligen tre zoner av olika sorters celler - längst ut på tumörens rand finns levande, växande celler, i mitten finns levande men vilande celler, och längst in i tumören finns en kärna av döda celler. Vidare så ser tillväxten i detta stadie ofta ut på ett specifikt sätt - först växer tumören mycket snabbt, vilket sedan åtföljs av en långsammare växtperiod. Experiment med cellsfäroider, som härmar tumörer i detta tillväxtstadium mycket väl, har även visat att tillväxten till slut mät- tas och tumören får en finit storlek.

När det ramverk som examensarbetet använt sig av skulle demonstreras valde skaparna av ramver-

ket att implementera en sådan, väldigt enkel, tumörmodell, men den tumör som modellen produc-

erade gav ingen finit tumör, utan uppvisade en fas där utskjutande "armar" uppkom på tumörens

kant. Den drivande forskningsfrågan i detta examensarbete har varit att undersöka om det med

en mer avancerad modell fortfarande uppkommer sådana utsprång eller inte.

Contents

1 Introduction 5

2 The DLCM framework 7

2.1 Description of the framework . . . . 7

2.2 Simple avascular tumour model . . . . 10

3 The Jagiella model 12 3.1 Individual-based description of single-cell dynamics . . . . 12

3.1.1 Proliferating cells . . . . 12

3.1.2 Quiescent cells . . . . 13

3.1.3 Necrotic cells . . . . 13

3.2 Continuum-based description of the dynamics of extracellular substances . . . . 14

3.2.1 Glucose and oxygen . . . . 14

3.3 Lactate . . . . 14

3.4 Extracellular matrix . . . . 14

3.5 Waste materials . . . . 15

4 Interpreting the Jagiella model into the DLCM framework 16 4.1 Event rates . . . . 16

5 Implementation 18 5.1 Methods for solving the continuous dynamics . . . . 18

5.1.1 Diffusion-reaction equations . . . . 18

5.1.2 ODEs . . . . 18

5.2 Tracking the state of cells . . . . 19

5.3 Boundary conditions . . . . 19

5.4 Updated algorithm . . . . 20

6 Results 22 6.1 Boundary conditions on far away boundary . . . . 22

6.2 Boundary conditions close to tumour . . . . 22

6.3 Effect of varying the conversion factors D 1 and D 2 . . . . 23

7 Discussion 25 8 Future work 25 8.1 Reproducability . . . . 26

A Parameters 29

1 Introduction

Cell-based compuational modeling is a tool that allows for easy exploration of how populations of cells behave under various circumstances, lending complimentary knowledge to biological experi- ments. They represent cells as discrete agents in a domain, where their state is updated over time with regards to a set of rules, often associated with the state of environmental conditions.

Cell-based models fall into one of the categories of being either on-lattice or off-lattice, this be- ing in regard to how they represent the position of individual cells within the modeling context.

On-lattice models constrain the domain by a fixed grid of polygonal shapes, where the cells reside within the lattice sites. Typical examples of on-lattice models include the cellular automaton, where the state of individual cells are determined by the their own and their neighbours present state, and cellular Potts [7], which uses energy minimization to determine the cells’ changes of state over time. Off-lattice models, on the other hand, do not constrain cells to be located within a discretised grid. In this case, a common way to describe the cell’s position is cell-centre models, where cells are seen as single points where energy potentials make sure cells stay within proximity of one another. Another example is vertex models, where cell geometry is approximated by a polygon but where the geometry does not remain static in the domain as in the case of on-lattice models, but the vertices are able to move according to deterministic rules of motion [8]. Compre- hensive reviews of on- vs. off-lattice models, including discussion of their respective strengths and weaknesses, can be found in [9] and [10].

Environmental conditions impacting the state change of cells, e.g. concentrations of extracellu- lar substances such as metabolites, are often governed by continuous-time processes modeled by ordinary differential equations (ODEs) or partial differential equations (PDEs). This leads to a so called hybrid model, merging discrete variables with continuous ones. This property means that hybrid models are well suited for modeling the interaction between individual cells and their environment, and the behaviour arising thereof. Two distinct categories of hybrid models arise, one that handles large populations of cells but with simplified geometry, and one where the geom- etry can be very complex but the population size stays small. The former allows for insights into collective population behaviour, while the latter is a suitable tool for exploring the relationship between cells and their neighbours, and cells and their environment [9].

Hybrid model or not, cell-based computational models have been used to model a wide variety biological processes, a few examples being blood vessel formation [13, 14], wound healing [11, 23], embryonic development [12, 22], and most notably in the context of this thesis - tumour growth [3, 9, 15, 19, 21].

As with all cells, tumour cells require oxygen to live and grow. The initial growth phase of tumours is, however, always avascular, meaning the tumour is cut off from a blood vessel that can supply it with oxygen. Later in the growth phase the tumour can begin creating its own blood vessels and thus cancer cells can escape into the blood stream of the body, causing metastasis where new tumours appear in places unrelated to the location of the initial tumour [15]. There has been a lot of research put into the modeling and understanding of the mechanisms of all these stages of tumour growth, but this thesis is concerned with the avascular growth phase of the tumours.

Characteristic of this growth phase is the emergence of three distinct zones of cell types, with pro- liferating cells on the rim of the tumour and a necrotic core in the middle, with a band of quiescent cells in between, a behaviour dependent on e.g. nutrient supply [16]. Another property that in vitro experiments with multicellular spheroids, which mimics the properties of avascular tumours very well, have shown is that after an initial exponential growth phase followed by a slower, linear growth, the size of the tumour saturates [3, 16, 17].

In [1] the authors present an on-lattice, cell-based computational framework for modeling the

dynamics of cell populations. The framework relies on the Laplace operator on a discrete grid to

inform cell migration, and thus is referred to as Discrete Laplacian Cell Mechanics (DLCM). When

demonstrating the functionality of the framework, they developed a simple model for avascular tu-

mour growth which when implemented results in a tumour exhibiting the three regions described

above. It did not, however, stop growing in the expected fashion, but instead the authors witnessed

another growth phase where asymmetrical protrusions appeared on the tumour, which caused the volume of the tumour to continue to increase. The authors speculate that this might be due to their model not including a mechanism for cell-cell adhesion, which regulates cell migration. The objective of this thesis is to further investigate the capabilities of the DLCM framework by in- tegrating a more advanced model for avascular tumour growth, specifically the one given in [2], which has been developed over time in [4, 3], into the DLCM framework. While the tumour model presented in [1] only took into account the oxygen levels in the tumour, this model considers the impact of far more extracellular substances and the results in [2, 3] do indeed exhibit both the cell heterogeneity and growth saturation described above.

Section 2 presents the DLCM framework, along with the simple avascular tumour model. The

model described in [2] is presented in Section 3, and Section 4 described how the model in [2] was

interpreted mathematically into the DLCM framework. Details on the implementation is described

in Section 5, with results of the implementation presented in Section 6. Finally, the results are

discussed in Section 7, followed by conclusions and thoughts on future work.

2 The DLCM framework

This section introduces the DLCM framework, modified for the task at hand. For a full and general description, see [1]. Section 2.2 presents the avascular tumour model implemented in [1], along with the results motivating the implementation of a more advanced tumour model.

2.1 Description of the framework

Consider a computational grid, v i discretised by i = 1, 2 . . . N vox voxels, which is populated by N cells cells. For the purposes of this thesis, this grid is cartesian. At any time, the voxels are either empty or contains a number of cells, u i . We constrain the number of cells allowed in each voxel to u i ∈ {−1, 0, 1, 2}, where u i = −1 denotes a dead cell. The mechanism at the core of the DLCM framework is that if the number of cells contained in a voxel is above the carrying capacity of the voxel, in this case 1, this will cause the cells to exert a pressure upon other cells. This pressure gradient will in turn induce a force that can cause cells to move between voxels. Likewise, if all voxels contain a number of cells below the carrying capacity, equilibrium is assumed. An overview of the disretised grid can be seen in Figure 1.

Figure 1: A schematic view of the numerical model, using an unstructured voronoi discretisation, with green voxels occupied by one cell and a red voxel occupied by two cells (left). One can think of the pressure being distributed evenly through springs connecting the voxel centres (right). The same derivation is used for a Cartesian grid. Adapted from [1].

The discretisation of the domain, along with a discrete counting of cells, gives a discrete state space where changes occur due to events that can be tracked over continuous time. To accomodate this state space a continous-time Markov chain is used to describe the governing physics, where the memoryless Markov property simulates the stochastic element of biological processes. The cell population model is built from three equations and three assumptions, which will be presented next.

Let u(t, x) represent the cell density at time t and point x. As mentioned above, the model is defined on a discrete grid of voxels, v i , with an integer number of cells u i . All the same, a continous notation is used throughout for ease of presentation. Consider the continuity equation

∂u

∂t + ∇ · I = 0, (1)

where I is a current, or flux. This equation will be used to derive the rates for the discrete events occuring in the Markov chain. To describe this current we now make our first assumption.

Assumption 2.1. The tissue is in mechanical equilibrium when all cells are placed in a voxel of their own.

This assumption allows us to neglect small Brownian-type movements performed by the cell around its center, while not excluding any active movements of cells induced by e.g. chemical stimulus.

This means that only voxels containing two cells will give rise to a rate to move. This increased

rate is described as a pressure source, which is set to unity identically as no other units are present.

Let p = p(t, x) denote the cellular pressure at time t and position x. The current I can be interpreted as the result of a pressure gradient, and we describe it with the phenomenological model

I = −D∇p, (2)

where we interpret D ≡ κ/µ, κ being the permeability and µ the viscosity. This can be interpreted as a fluid flow, but since the cells retain their shape, the current I will instead be taken as the rate of the discrete event that a cell will move into a new voxel.

Finally, we need an equation to describe the relationship between u and p. Considering that changes in pressure are induced by overcrowded voxels, we use the heat equation to model pressure evolution

ε ∂p

∂t = ∆p + s(u), (3)

where s(u) is a source function, described below. Eq. (3) follows by assuming an isotropic medium and a scalar pressure. For voxels populated at or below the carrying capacity, there is no net flux of the potential and the divergence theorem implies the Laplace operator. For an overpopulated voxel, there is instead a net outward flux, then captured via the divergence theorem as a source term. We simplify eq. (3) by making an additional assumption.

Assumption 2.2. The cellular pressure of the tissue relaxes rapidly to equilibrium in comparison with any other mechanical processes of the system.

This assumption gives that any processes affecting the mechanics of the individual cells, apart from the propagation of the cellular pressure, must happen on a slower time-scale. This allows the simplification of eq. (3) into the Laplacian equilibrium

− ∆p = s(u). (4)

Boundary conditions can be chosen according to what biological problem one wants to model. For this thesis, Dirichlet boundary conditions are applied,

p  

 _∂Ω = 0, (5)

The domain consists of the bounded subset of R ² which is populated by cells, along with its bound- ary ∂Ω.

As we are working on a discrete grid v i , we need to adapt our model onto this grid. Let Ω h

be the subset of voxels v i for which u i 6= 0, i.e. all voxels that are not empty. Further, let ∂Ω h

denote the discrete boundary. This is the set of empty voxels that share one edge with one voxel in Ω _h . We thus consider the discrete version of eq. (4),

−Lp = s(u), i ∈ Ω _h , (6)

p i = 0, i ∈ ∂Ω h . (7)

As p = 0 for steady state, the source term s(u) is given by ( s(u i ) = 0 u i ≤ 1,

s(u _i ) = 1 u _i = 2. (8)

L in eq. (6) is the discrete Laplacian operator over the grid Ω _h . For this thesis, L is defined using finite elements with linear basis functions and a lumped mass matrix, which simplifies to traditional finite difference stencils on structured grids. It follows that the continuous interpretation of the source term is formally consistent with a Dirac delta function, which is either one or zero depending on the whether the voxel in question is above or below carrying capacity; in continuous language

s(u) = s(u(x)) = X

u

>c

δ v

(x), (9)

where v i here denotes the center point of the ith voxel and c the carrying capacity.

In the discrete setting we have now introduced, we define the current I induced by pressure changes. Eq. (2) gives that I ∝ −∇p. Let I(i → j) denote the current from voxel v i to neighbour v j . The current is calculated by integrating the pressure gradient across the edge between two voxels,

I(i → j) ∝ − Z

v

∩v

∇p(x)d~ S = e _ij d ij

(p i − p j ) := I ij , (10) where d ij is the distance between voxel centres and e ij the shared edge length. The simplification considers the discrete nature of the model, and thus the gradient is simply taken as a scaled difference of pressure in between voxels. Eq. (10) giving a negative result means the current is reversed, I(j → i). This means that if there is a pressure source active in the domain, there will be a non-zero pressure gradient everywhere. We now go on to specify rules for the movement induced by the pressure gradient, along with D in eq. (2).

Assumption 2.3. The cells in a voxel occupied with n cells may only move into a neighbouring voxel if it is occupied with less than n cells.

Let R(e) be the rate of the event e, and use the notation i → j for the particular event that one cell moves from voxel i to j. In this case there are two movement rates; one for cells moving into an empty space, and one for cells moving into already occupied voxels.

R(i → j; u i ≥ 1, u j = 0) = D 1 I ij , (11) R(i → j; u _i > 1, u _j = 1) = D ₂ I _ij , (12) where D ₁ and D ₂ are conversion factors from a unit pressure gradient into a movement rate for the three cases.

To understand the scale of D 2 we note that from the divergence theorem we have

− Z

∂Ω

(∇p · n) d ~ S = Z

Ω

s(u) dV, (13)

i.e. the total pressure gradient over any closed surface ∂Ω is equal to the enclosed sources. With

∂Ω the boundary surrounding the whole cell population, the ratio D ₁ /D ₂ expresses how likely it is

that an event occurs on the boundary, with cells entering empty voxels, rather than in the internal

regions of the population, with cells in doubly occupied voxels moving into an already occupied

voxel. For example, if we assume a single internal source voxel, and that D 1 /D 2 ≡ 1, then the

probability of a cell in the source voxel to move is the same as the total probability of a cell to

move into any of the boundary voxels in ∂Ω.

Algorithm 1 Outline of the DLCM simulation methodology

1: Intialise: at time t = 0, given a state (u _i ), u _i ∈ {0, 1, 2} over (a subset of) the mesh of voxels (v _i ), i = 1, . . . , N _v ox.

2: while t < T do

3: Solve for the cellular pressure (p _i ), eqs. (6)-(7). Compute all movement rates (r _j ) for the subset of voxels where the cells may move, eqs. (11)-(12).

4: Determine also the rates for any other processes taking place in the model such as prolifer- ation and death events, or active migration.

5: Compute the sum λ of all transition rates thus computed.

6: Sample the next event time by τ = − log(U 1 )/λ, for U 1 a uniformly distributed random variable in 0, 1.

7: Determine which event happened by inverse sampling: find n such that P n−1

j=1 r j < U 2 λ ≤ P n

j=1 r j , for U 2 a second U (0, 1)-distributed random number.

8: Update the state of all cells with respect to any other continuous-time processes taking place in [t, t + τ ), e.g. intracellular kinetics or cell-to-cell communication.

9: Update the state (u i ) by executing the state transition associated with the determined event.

10: Set t = t + τ .

11: end while

12: Result: a sampled outcome of the system consisting of states observed at discrete times (t k ) ∈ [0, T ].

Algorithm 1 summarises the simulation method in pseudo code. For any given state of the cell population, the cellular pressure is calculated and the rates of all possible events are determined (lines 3-4). Here the sampling procedure of Gillespie [18] may be used; the sum of all the rates decides the time for the next event, and a proportional sampling next determines the event that happens (lines 5-7). Until the time of this next event, any other processes local to each voxel may be simulated in an independent fashion (line 8). Finally, as the event is processed, a new cell population state is obtained and the loop starts anew.

In addition to the rates in eqs. (11) -(12), one can easily extend the basic framework to include active cell movements as prescribed via additional rates and to be handled in line 4 of Algorithm 1. Depending on the details of the application, cell-cell adhesion and surface tension effects could also be added in this direct way. Another option to include similar effects could be to modify the pressure source function and/or the boundary conditions.

For this thesis, several rates relating to e.g. proliferation and death events will be added di- rectly at line 4 in Algorithm 1. Further, additional continuous dynamics relating to extracellular substances will be performed at line 8. All details about these rates and dynamics are described in Section 3, and a more detailed version of this algorithm is presented in Algorithm 2.

2.2 Simple avascular tumour model

In [1] the authors demonstrate the DLCM framework by presenting results from an implementation of a simple model for avascular tumour growth, governed by the Laplace equation and a set of rules regarding the oxygen concentration. This section presents details on this model, along with examples of the growth pattern it produces.

We imagine the tumour as a population made up of proliferating tumour cells, i.e. cells that are growing and ultimately divides into two new cells causing the population to grow, and qui- escent cells that are dormant within the population. Over the course of the simlation these cells can die and become necrotic. As mentioned above, in the case of an avascular tumour there is no direct supply of oxygen to the tumour itself. Instead, oxygen in the surrounding tissue diffuses into the tumour. The concentration of oxygen is thus modeled using the Laplace equation, with a boundary condition on the far boundary ∂Ω _ext , modeling the external oxygen supply,

−Lc = −λa(u), (14)

c i = 1, i ∈ ∂Ω ext , (15)

where λ is the rate of consumption of oxygen for a single cell and a(u i ) is the number of alive cells (u i ∈ {0, 1, 2}) in the i th voxel. Effectively, this means doubly occupied voxels consume twice as much oxygen while empty voxels or those containing dead cells consume none. The cells at the boundary of the tumour, ∂Ω h , can consume oxygen and proliferate if the concentration of oxygen is above some threshold value. Further from the rim of the tumour, the oxygen concentration will be lower and will not allow for proliferation. However, the cells will still be alive and consume oxygen in a quiescent state. Near the core of the tumour the oxygen will be low enough that the cells become necrotic and then disintegrate and are removed from the population.

Figure 2: Snapshots of the time evolution of the tumour at times a) 0 hours, b) 150 hours, c) 300 hours, d) 700 hours and e) 1000 hours, along with a graph of cell type count over time (bottom right). Adapted from [1].

From eqs. (14)-(15) we can calculate the rates for all possible events, which are as follows. A cell that is alone in its voxel, u _i = 1, will proliferate at ρ _prol if c _i > κ _prol , where κ _prol is the minimum oxygen concentration for proliferation to occur. A living cell will die at a rate ρ death if c i < κ death , where κ death is the threshold oxygen concentration for cell survival. In the case of cell death, u i = 1 is replaced with u i = −1. A dead cell, u i = −1, can degrade at a constant rate ρ deg to free up the voxel (u i = 0) which other cells can then move into again. Along with these event rates we have the following movement rates,

R(i → j; u i ≥ 1, u j = 0, v j never visited) = D 1 I ij , (16) R(i → j; u i ≥ 1, u j = 0, v j visited previously) = D 2 I ij , (17)

R(i → j); u i > 1, u j = 1) = D 3 I ij , (18)

where D 1 and D 2 have different values, as the extracellular matrix, a substance that surrounds

cells, of unvisited voxels is thought to be stiffer than if a voxel has previously been occupied. Figure

2 shows snapshots from the simulation of the simple model. Figures 2a-2d all exhibit a population

in accordance with the characteristic behaviour of avascular tumours, while Figure 2e shows the

unexpected asymmetric protrusions that cause an increase in tumour size.

3 The Jagiella model

This section introduces the model to be implemented, presented by Jagiella et al. in [2]. It is an on-lattice multi-scale model describing in-vitro tumour growth, meaning the growth takes place outside of its biological context. It is a hybrid model, so it uses an individual-based description of tumour cells, with PDEs modeling the extracellular substances in the medium surrounding the cells. The cells are modeled so that they can sense their environment, move, divide and die. The cells can move to an empty neighbouring lattice site with rate k mig , or by being pushed by grow- ing cells, as described in Section 3.1.1. The communication between single cells depend on direct cell-cell contact and uptake/secretion of the extracellular substances. As the model presented here is based on previous work by Jagiella in [4], the model will henceforth be referred to as the Jagiella model.

Some modifications of this model had to be done in order to accommodate the characteristics of the DLCM framework. However, this section will present the Jagiella model in its entirety, while the modifications will be addressed in Section 4.

3.1 Individual-based description of single-cell dynamics

The cells considered can take on three states: proliferating, quiescent and necrotic, and they live in a static, unstructured lattice. A lattice site x can only be occupied by one cell at any time t.

The dynamics of the cell located at site x at time t is governed by the time-dependent local con- centrations of the extracellular substances: oxygen o(t, x), glucose g(t, x), adenosine triphosphate (ATP) k a (t, x), lactate l(t, x), extracellular matrix (ECM) e(t, x), waste material from debris of necrotic cells w(t, x), and the distance L(t, x) to the next vacant lattice site. Figure 3 shows an overview of what events can happen to the three types of cells presented in this section.

Initially, the domain is populated by cells within a sphere of radius L init , located in the cen- ter of the lattice. Cells are either quiescent with probability q init , or otherwise proliferating and in the first stage of the cell cycle.

Notation: Dynamics in this section makes use of the Heaviside function, defined as

H[x] =

( 0, if x < 0, 1, if x ≥ 0.

3.1.1 Proliferating cells

Proliferating cells progress in a discretised cell cycle with m d stages, where m = 1, . . . , m d , and m d

denotes the stage where the cell will divide. Cells transition from stage m → m + 1 with propensity k div,m (t, x) = k ^div _max m d

 1 − 1

2 H[o ^div − o(t, x)] 

1 − 1

2 H[w(t, x) − w ^div ] 

, (19)

where k _max ^div is the maximal division rate, o ^div the oxygen division threshold and w ^div the waste division threshold. This propensity increases with the availability of oxygen, while the presence of waste decreases it. When the cell reaches stage m = m _g , the cell grows so that it occupies an additional lattice site, next to its original site. If no such site is available, the neighbouring cells are pushed along the shortest path toward the closest vacant lattice site.

Once a cell reaches m = m d it divides into two daughter cells, which each decides to either proliferate or become quiescent according to the probability

p re (t, x) = exp h

− L(t, x) L ^div

i

H[e(t, x) − e ^div ]H[k a (t, x) − k _a ^div ]H[l(t, x) − l ^div ]H[n w,o,max − n w,o (t, x)],

(20)

where L ^div is the division depth, e ^div the ECM division threshold, l ^div the lactate division threshold,

k _a ^div the ATP synthesis division threshold, n w,o,max the maximum time under waste exposure/oxy-

gen deprivation that allows proliferation, and n _w,o (t, x) is the time of oxygen deprivation and waste

exposure, defined as

n _w,o (t, x) = Z t

0

1 − H[w ^div − w(τ, x(τ ))]H[o(τ, x(τ )) − o ^div ] dτ, (21) where x(τ ) denotes the time-dependent spatial location of the cell located at site x at time t.

Note that this formulation of eq. (21) has been corrected, since the equation in [2] had the ar- gument in the second Heaviside function flipped. That formulation would not be consistent with what it represents, and also, in [3] the authors state the above formulation for calculating n w,o (t, x).

The ATP synthesis rate k a (t, x) is defined as

k a (t, x) = 2q g (t, x) + 17

3 q o (t, x), (22)

where q _g (t, x), q _o (t, x) is glucose and oxygen consumption respectively, defined in Section 3.2.1. If a daughter cell becomes proliferating, it begins at cell cycle stage m = 1.

3.1.2 Quiescent cells

Quiescent cells are cells that have become dormant in the cell cycle, meaning they do no longer grow. They can, however, reenter the cell cycle, which they try to do with propensity

k re (t, x) = k _max ^re  1 − 1

2 H[w(t, x) − w ^div ] 

1 − 1

2 H[o ^div − o(t, x)] 

, (23)

and succeeds with probability p _re (t, x). If it succeeds it once again becomes proliferating with m = 1. k ^re _max denotes the maximum reentry rate.

3.1.3 Necrotic cells

Proliferating and quiescent cells become necrotic with propensity k _nec (t, x) = k _max ^nec H[k _a ^nec − k a (t, x)] l(t, x) ²

l(t, x) ² + (l ^nec ) ² , (24) where k _a ^nec is the ATP synthesis necrosis threshold and l ^nec is the lactate necrosis threshold. This propensity increases with local lactate concentration, if the ATP synthesis rate is low. Necrotic cells are subject to lysis, i.e. disintegration, according to the constant propensity k lys . If lysed, a cell is then removed from its lattice site.

Living cells

Proliferating cells progress

division

m = 1, 2, …, m

1-p

p

Quiescent cells reentry

migration

death Necrotic cells

lysis k

k

p

k

k

k

Figure 3: Overview of what events can happen to a certain type of cell in the Jagiella model.

3.2 Continuum-based description of the dynamics of extracellular sub- stances

This section provides details on the dynamics of the extracellular substances. They are governed by PDEs, and all depend on the state of the cell occupying lattice site x at time t.

3.2.1 Glucose and oxygen

Glucose and oxygen are the cell’s primary energy sources, and are modeled by diffusion-reaction equations

∂g(t, x)

∂t = D g ∆g(t, x) − q g (t, x)c(t, x), (25)

∂o(t, x)

∂t = D _o ∆o(t, x) − q _o (t, x)c(t, x), (26) where D g , D o are diffusion coefficients, c(t, x) is the indicator function which is 1 if a proliferating or quiescent cell occupies x at time t, and q g , q o are the local consumptions of glucose and oxygen respectively, given by

q g (t, x) = V m,g (t, x) g(t, x)

g(t, x) + k _m,g (27)

q o (t, x) = V m,o (t, x) o(t, x) o(t, x) + k m,o

, (28)

where k m,g , k m,o are the Michaelis-Menten constants. V m,g and V m,o are the maximum consump- tion rates, given by

V m,g = q ^max _g  1 − 

1 − q ^min _g q ^max _g

 o(t, x) o(t, x) + k o



, (29)

V m,o = q ^max _o  1 − 

1 − q _o ^min q _o ^max

 g(t, x) g(t, x) + k _o



. (30)

As can be seen, these rates account for cross-dependencies between glucose and oxygen. It has been observed that cells lacking in one metabolite will compensate for this by consuming more of the other one, so that the net production of ATP molecules remain constant. q _g ^min , q _g ^max , k g , q ^min _o , q _o ^max and k o are consumption parameters.

Glucose and oxygen enter the simulation domain Ω from the surrounding medium, and we as- sume Dirichlet boundary conditions, g(t, x) = g 0 and o(t, x) = o 0 ∀ x ∈ ∂Ω. Further, the initial conditions are equivalent to the boundary conditions, g(0, x) = g 0 and o(0, x) = o 0 ∀ x ∈ Ω.

3.3 Lactate

Lactate is produced by living, i.e. proliferating and quiescent, cells and is a by-product of the anaerobic energy metabolism. It is modeled by

∂l(t, x)

∂t = D _l ∆l(t, x) + 2 

q _g (t, x) + min n

q _g (t, x), 1

6 q _o (t, x) o

c(t, x), (31)

where D l is the lactate diffusion coefficient. We assume that lactate dilutes, and zero Dirichlet boundary conditions, l(t, x) = 0 ∀ x ∈ ∂Ω, along with zero initial concentration, l(0, x) = 0 ∀ x ∈ Ω.

3.4 Extracellular matrix

Extracellular matrix is involved in many processes regarding the structural integrity of the cells, such as cell adhesion, and cell-to-cell communication. The extracellular matrix is produced by living cells and can be degraded, which yields the following model

∂e(t, x)

∂t = k ^pro _e c(t, x) − k ^deg _e e(t, x), (32)

where k ^prog _e , k _e ^deg is the production rate and degradation rate, respectively. These are assumed

to be constant, and are expressed in units of intensity, since the absolute extracellular matrix

concentration cannot be assessed experimentally. Zero Dirichlet boundary conditions are assumed,

e(t, x) = 0 ∀ x ∈ ∂Ω, and zero initial concentration e(0, x) = 0 ∀ x ∈ Ω.

3.5 Waste materials

Waste materials are produced by necrotic cells and absorbed by living cells with a constant rate, which yields the following model

∂w(t, x)

∂t = k _w ^pro c ^nec (t, x) − k _w ^upt w(t, x)c(t, x), (33)

where k ^pro _w , k ^upt _w is the waste production rate and waste uptake rate, respectively. c ^nec (t, x) is the

necrotic indicator function, which is 1 if a necrotic cell occupies x at time t, and 0 otherwise. As no

necrotic cells exist initially, the initial waste concentration is zero, w(0, x) = 0 ∀ x ∈ Ω. Further,

as waste is not transported and as there are no cells on the boundary, zero Dirichlet boundary

conditions are assumed, w(t, x) = 0 ∀ x ∈ ∂Ω.

4 Interpreting the Jagiella model into the DLCM framework

This section describes the mathematical model that comes from merging the framework with the Jagiella model, and what we arrive at is a 2D model of in vitro tumour growth. Incorporating the Jagiella model into the framework was straightforward in the case of event rates, Section 4.1, but there were some exceptions.

Firstly, as stated in Section 2, we allow the occupancy of voxels to be u i ∈ {−1, 0, 1, 2}, in contrast with the constraint set in [2] where doubly occupied voxels are not allowed. Allowing this means we do not need a growth state, m = m _g , in the cell cycle where the cell grows to occupy more than one lattice site, before it then divides into two cells at m = m _d . However, we only allow cells with m = m _d in singly occupied voxels to divide, which mimics this constraint and the growth stage, as this can be seen as one cell "taking up" two sites where cells might reside.

This change also propagates to the indicator functions c(t, x) and c nec (t, x). Since these func- tions decide what changes will happen in which voxels, it makes sense that if a voxel contains two cells this change will be doubled, as is the case in the model described in 2.2. That is, if a voxel contains two cells, twice the amount of e.g. oxygen will be consumed in this voxel. Likewise, if a voxel contains two dead cells, twice the amount of waste will be produced in this voxel. Thus we have c(t, x) ∈ {0, 1, 2} and c nec (t, x) ∈ {0, 1, 2}.

Secondly, the dependence on L(t, x), and L div in eq. (20) is substituted for a dependence on pressure,

p re (t, x) = exp

"

−p(t, x) p ^div

#

H[e(t, x) − e ^div ]H[k a (t, x) − k ^div _a ]H[l(t, x) − l ^div ]H[n w,o,max − n w,o (t, x)], (34) where p ^div is the pressure division threshold, and p(t, x) is the pressure found in voxel x at the time t. The value of the threshold value p _div is decided using the value provided for L _div . By creating a population of cells with a radius the size of L _div , with a single source voxel in the middle, we solve for the pressure. The resulting pressure in the source voxel is then taken as p div , as this is what the pressure would be in the voxel if there are free voxels at a distance L div away from the source vox- els.Thus, if p(t, x) is higher than this value, the probability that a cell will be proliferating decreases.

This change is motivated by differences in pressure being, as proposed in [1], a natural way for cells to sense their surroundings. Also, the DLCM framework is built around the idea that changes in pressure lead to cell migration, to not take advantage of this functionality would be counter- intuitive.

No changes were made to the spatio-temporal equations, (25)-(33), describing the evolution of the extracellular substances. They were taken as they are described in Section 3.2 and implemented accordingly.

4.1 Event rates

As the Jagiella model also uses the Gillespie algorithm to simulate the individual cell dynamics, the propensities given by equations (19), (23), (24), along with k lys , were simply taken as event rates. In addition to the movement rates, eqs. (11) and (12), five more rates arise which are as follows.

The rate with which a cell progresses in the cell cycle is given by

R(m → m + 1) = k _div,m , (35)

for cells that have 0 < m < m _d , otherwise zero.

The rate with which a mother cell, c m , divides into two new daughter cells, c d , is given by

R(c m → c d + c d ) = k div,m , (36)

for cells that have m = m d and reside in singly occupied voxels, otherwise zero. This rate is also k div,m since division is just the final step in the cell cycle progression.

The rate with which a quiescent cell successfully reenters the cell cycle is given by

R(m = 0 → m = 1) = p re k re , (37)

for cells that have m = 0.

The rate with which a cell dies and becomes necrotic is given by

R(m ≥ 0 → m = −1) = k _nec , (38)

for any living cell, either proliferating at any stage in the cell cycle or quiescent.

Finally, the rate with which a necrotic cell is removed from the population is given by

R(u i = −1 → u i = 0) = k lys . (39)

5 Implementation

The rates presented in Section 4.1 were implemented the same way as the movement rates given by eqs. (11)-(12), and events sampled using the Gillespie algorithm. A "null"-event was also added, to ensure that if for some reason all event rates would be zero, e.g. if the population consists of only quiescent cells and the process is waiting for one of them to reenter the cell cycle, the time step does not become infinite. Instead the simulation continues with nothing else happening except the continuous dynamics being recalculated and a new event sampled. The implementation has been done using Matlab and URDME [5, 6], and the parameters used are tabulated in Appendix A.

Details on the methods used for solving for the continous dynamics are presented in Section 5.1.

Section 5.2 describes the method used for tracking the state of individual cells. Some rearranging of the simulation process had to be done to incorporate continuous calculations of the extracellular substances. Section 5.4 presents an updated version of Algorithm 1, taking into account the details of the Jagiella model.

5.1 Methods for solving the continuous dynamics

The continuous dynamics were modeled using finite difference schemes in the case of the diffusion- reaction equations, and for the ODEs analytic solutions were produced. For each new event, the time-discretised equations were time stepped up until the next event time.

5.1.1 Diffusion-reaction equations

Eqs. (25)-(26) and (31) are all diffusion-reaction equations, where eqs. (25)-(26) were solved using a semi-implicit Euler method, where the diffusion term is implicit, while the reaction term is explicit.

This was motivated by a need for stability, along with a wish for simplicity. As explicit Euler gave way to instability, and we wished to avoid implicit Euler as it would be too time-consuming, we settled on a method in between. For an arbitrary function u(t, x) with u ⁿ ≈ u(t n , x) we can write the scheme as

u ⁿ⁺¹ − u ⁿ = ∆t(Du ⁿ⁺¹ − R(u ⁿ )) (40)

where ∆t is the time step, D denotes a diffusion operator on matrix form and R(u ⁿ ) models reaction. This semi-implicit scheme can be rearranged to give a linear relationship on the form Au = b, which can be solved easily using LU-decomposition and Matlab’s backslash operator.

In this case, D is the discrete Laplacian L described in Section 2, multiplied with an arbitrary diffusion constant. The time step was decided according to

∆t = 1

λ _D , (41)

where λ _D ≈ ||Du ⁿ ||/||u ⁿ ||, where the norm in this case is the Euclidean norm. Eq. (31) was solved by assuming equilibrium, in the same way we solve for pressure in the framework,

0 = D _l ∆l(t, x) + 2 

q _g + min n q _g , 1

6 q _o o

. (42)

This could be done for eqs. (25)-(26), but since these equations are non-linear, it was deemed a simpler solution to instead discretise them in the manner described above.

5.1.2 ODEs

Eqs. (32) and (33) are ODEs that were solved analytically. The solution to eq. (32) is given by e(t) = e T exp[−k ^deg _e (t − T )] + k _e ^pro

k ^deg e

c(t, x) 

1 − exp[−k ^deg _e (t − T )] 

, (43)

where

e _T = e _init exp[−k _e ^deg T ] + k _e ^pro k e ^deg

c(t, x)(1 − exp[−k ^deg _e T ]), (44)

and T is the time for the previous event, and e init is the previous solution. Eq. (33) can be solved in a similar fashion. However, doing this would give a solution that would give division by zero when c(t, x) = 0. Instead, a limit of eq. (32) was chosen so that the solution instead becomes

( w(t, x) = w(t − 1, x) + k ^pro _w t, if c nec (t, x) = 1,

w(t, x) = w(t − 1, x) exp[−k _w ^upt t], if c nec (t, x) = 0, (45) where w(t − 1, x) denotes the previous solution.

5.2 Tracking the state of cells

As all propensities and continous dynamics described in Section 3 depends on the state of cells, a system for tracking this was needed. For this purpose, a N vox × 2 identificaton matrix, ID, was created, where the second column is only used if a voxel contains two cells. The matrix entries can be categorised into four cases



 

 

 

 

ID _xk = 1, ID xk = m + 1 ID _xk = −1, ID xk = 0,

(46)

with k ∈ {1, 2}, and where the first case signifies quiescent cells, the second proliferating cells with m = 2, . . . , m d + 1 being which stage (minus one) in the cell cycle they are in, the third dead cells and the fourth signifies an empty cell. These were chosen so to make use of Matlab’s sparse function in the best possible way.

To mimic a natural flow in the dynamics of the cells, movement is always prioritised so that the cell in the first column of ID is the first to move, in a first-in-first-out approach. This also means that when the cell in the first column moves it is then placed in the second column for the new voxel. In turn, in the voxel the first cell moved from, the cell represented in the second column will then be moved to the first column.

5.3 Boundary conditions

For the simple model presented in 2.2, Dirichlet boundary conditions were placed on a "far away"

circle surrounding the tumour. What far away means depends on the resolution of the grid, with the side of the mesh measuring N c d , where N = √

N vox and c d is the cell diameter. When placing the boundary conditions far away, they are always placed on a circle with radius c d N/2. Through correspondence with the authors of [2] we learned that they place their Dirichlet boundary condi- tions much closer to the tumour, on a radius that allows for the measured size of the final tumour, r _{f inal} ≈ 450µm, along with some margin.

To determine the impact the placement of boundary conditions has on the concentrations de- scribed in eqs. (25) and (26), a small numerical experiment was performed. The following initial value boundary problem,

u t (t, x) = u xx (t, x) − 0.5u(t, x) if x ∈ [0, 0.5], (47)

u _t (t, x) = u _xx (t, x) otherwise, (48)

u(0, x) = 1, (49)

u _x (t, 0) = 0, (50)

u(t, L) = 1, (51)

where u(t, x) models some arbitrary concentration, was implemented and solved for two different

intervals, I = [0, L 1,2 ], where L ₁ = 1 and L ₂ = 10, using explicit Euler for time stepping and

a central difference for the second derivative. Figure 4 shows the result for a run with end time

T _end = 10 ⁴ .

Figure 4: Result of numerical test to determine the significance of boundary condition placement, where a 1D diffusion-reaction equation was solved on two different sized intervals, I = [0, 1] and I = [0, 10].

We see that at x = 1 there is a difference in concentration for the two cases. For the case I = [0, 1], x(T end , 1) = 1, while for the case I = [0, 10], x(T end , 1) ≈ 0.89. This means that there is indeed a correlation between how far boundary conditions are placed from the area of reaction, in this case x ∈ [0, 0.5], and the concentration surrounding this area, which is what matters in this case (note that the concentration at x = 0 seems to be the same for both cases). As the tumour relies on uptake of nutrients on and close to its rim in order to grow and proliferate, it seems the placement of the boundary condition is indeed very relevant for the growth pattern of the tumour. Both the

"far away" and the "close" way of placing the boundary conditions were implemented. For the close case, the boundary conditions were placed at r = 512µm by decreasing the disretisation so that c _d N/2 = 512µm.

5.4 Updated algorithm

Algorithm 2 presents an updated outline of the simulation process, taking into account the details of Jagiella’s model. As stated in Section 2, changes in the cellular pressure are assumed to occur instantaneously. Eqs. (25)-(26) however, need to be time-stepped, lines 7-10 in Algorithm 2.

They are time stepped with ∆t up until the next event time τ . As eqs. (32)-(33) have analytical

solutions, these do not need to be time stepped, but are instead solved at time t + τ . The same

goes for eq. (31), as it is solved assuming equilibrium. The event is then executed, and new event

rates are calculated before the process starts anew.

Algorithm 2 Updated simulation outline

1: Intialise: at time t = 0, given a state (u _i ), u _i ∈ {0, 1, 2} over (a subset of) the mesh of voxels (v _i ), i = 1, . . . , N _vox and initial conditions for extracellular substances.

2: Solve for cellular pressure (p _i ), eqs. (6)-(7). Compute all initial event rates: movement, birth, progress, reentry, death and degradation.

3: while t < T do

4: Compute the sum λ of all event rates.

5: Sample the next event time by τ = − log(U 1 )/λ, for U 1 a uniformly distributed random variable in (0, 1).

6: Determine which event happened by inverse sampling: find n such that P n−1

j=1 r _j < U ₂ λ ≤ P n

j=1 r _j , for U ₂ a second U (0, 1)-distributed random number.

7: for k < τ do

8: Simulate continuous-time processes taking place in [t, t + τ ), eqs. (25)-(26).

9: Set k = k + ∆t.

10: end for

11: Simulate continuous-time processes taking place in [t, t + τ ), eqs. (31)-(33).

12: Update the state (u _i ) by executing the state transition associated with the determined event.

13: Solve for the cellular pressure. Compute all event rates: movement, birth, progress, reentry, death and degradation.

14: Set t = t + τ .

15: end while

16: Result: a sampled outcome of the system consisting of states observed at discrete times (t k ) ∈

[0, T ].

6 Results

This section presents the results of simulations using the implemented model for two different placements of Dirichlet boundary conditions, Sections 6.1, 6.2. Section 6.3 presents a comparison of final tumour size when varying the parameters D 1 and D 2 presented in eqs. (11)-(12). The tumours are simulated for 960 hours, as this is what is done in [2, 3], and using the parameters for what is termed "Scenario II" in [2]. In [2, 3] they observe saturation at approximately 580 hours.

6.1 Boundary conditions on far away boundary

Figure 5 show snapshots of the tumour’s evolution over time when the boundary conditions are placed far away. We can see that the three zones of cell types appear in some fashion, with necrotic cells appearing in the middle at around 100 hours and proliferating cells on the rim. However, much like the results presented in Figure 2, the tumour becomes asymmetrical. In this case, one protrusion on the left side appears, which then splits into two branches, with no proliferation anywhere else on the rim. Further, comparing the bottom right graph of Figure 2 and Figure 5f we see that they both show the same behaviour - fast growth in the beginning that slows down as proliferating cells decrease and necrotic cells appear, but there is no indication of the growth decreasing. The final maximum radius of the tumour is r max = 387µm.

(a) (b) (c)

(d) (e) (f)

Figure 5: Snapshots of the time evolution of the tumour when boundary conditions are placed far away, at times a) 0 hours, b) 48 hours, c) 278 hours, d) 662 hours and e) 960 hours, along with f) a graph of cell type count over time. Green signifies a single quiescent cell, orange double quiescent cells, yellow a single proliferating cell, pink double proliferating cells and black necrotic cells. Note that in Figure 5f all proliferating cells are included in orange, and all quiescent cells are included in green.

6.2 Boundary conditions close to tumour

Figure 6 show snapshots of the tumour’s evolution over time when the boundary conditions are placed close to the tumour. Here we can clearly see the three zones of cell types in Figure 6d.

However, we should already be seeing the growth saturating at this point, but instead we see a sharp increase in the number of cells, with no indication of saturation at the end of the simulation.

The final maximum radius is r max = 451µm.

(a) (b) (c)

(d) (e) (f)

Figure 6: Snapshots of the time evolution of the tumour when boundary conditions are placed close to the tumour, at times a) 0 hours, b) 48 hours, c) 278 hours, d) 662 hours and e) 960 hours, along with f) a graph of cell type count over time. Green signifies a single quiescent cell, orange double quiescent cells, yellow a single proliferating cell, pink double proliferating cells and black necrotic cells. Note that in Figure 6f all proliferating cells are included in orange, and all quiescent cells are included in green.

6.3 Effect of varying the conversion factors D ₁ and D ₂

Tests showed that varying the size of the conversion factors D 1 and D 2 in equations (11) and (12) had a significant impact on how quickly the tumour grows. This is not very surprising, since D 1

and D 2 directly affect the size of the migration rate. However, it should be noted that for the tumour not to grow outside of the domain when the boundary conditions were placed close to the tumour, D ₁ and D ₂ had to be set to values that would in many cases give a migration rate that is lower than the one reported in the supplementary information of [3].

In Figure 7 we see the results for three cases of varying D 1 and D 2 , with boundary conditions close to the tumour. It is evident that very small changes in these parameters cause significant changes in the final size of the tumour. The difference in D 1 , D 2 for the cases presented in Figures 7a and 7b is only 0.02/h, and yet the final radius in the case of Figure 7b has almost double the final radius that the tumour in Figure 7a has, r max,b = 410µm while r max,a = 228µm. The tumour in Figure 7c has D 1 , D 2 ten times larger than for the case in Figure 7a, and here we see that the tumour grows beyond the domain. It is thus not a valid result, but is shown here to emphasise the impact varying D 1 and D 2 has on the growth of the tumour.

Further, we note that for none of these cases do we observe an indication of saturating growth,

Figures 7d-7f. Even for the case with the smallest value for D 1 , D 2 , Figure 7d, we see a perpetual

increase in the number of cells.

(a) (b) (c)

(d) (e) (f)

Figure 7: Final population appearance and change in number of cells over time for a) + d)

D 1 = D 2 = 0.01, b) + e) D 1 = D 2 = 0.03, and c) + f) D 1 = D 2 = 0.1.

7 Discussion

We have implemented a model for avascular tumour growth by merging an existing model with a computational framework for cell population dynamics. The goal of the work has been to in- vestigate whether the tumours produced by this merged model would be symmetric and obtain a finite size. We have presented results from simulations with boundary conditions both far away and close to the tumour, along with results displaying the sensitivity of the parameters D 1 and D 2 . As the results presented in Section 6 show, we could not reproduce the results presented in [2, 3], as growth does not saturate in any case and in all cases the tumour becomes asymmetric.

The final maximum radius of the tumour presented Section 6.2 is in accordance with results in [2], but still asymmetric. Several other radii for the placement of boundary conditions were explored, but none yielded a satisfactory solution. Placing the boundary conditions closer than what was presented in Section 6.2 causes the tumour to grow larger than the domain it resides in.

According to [3], the mechanism that controls saturation of growth is eq. (21), which says that cells that have been exposed to waste or deprived of oxygen for longer than some time n w,o,max

will become quiescent and not reenter the cell cycle again. No fault could be found in our imple- mentation of this mechanism. Instead we pose the hypothesis that what is causing this difference in results lies in how cell migration is handled by the DLCM framework as opposed to how it is done in the Jagiella model, as this is the major difference in the two models.

As it is, whether the rate for migration is non-zero or not is much more restricted in the Jagiella model. A cell can move in only two cases - if there is a free voxel adjacent to it, or if a neighbouring cell has grown and needs to occupy a voxel next to it. In the former case, the cell has rate k mig to move. In the latter case however, the rate to move is rather something akin to being infinite, as it is required that cells neighbouring the growing cell must move in order to make room for that cell.

This means that movement of cells that are not on the boundary of the tumour and have empty spaces adjacent to them, will be limited to having a movement rate only when their neighbouring cells grow.

In the case of the DLCM framework, whenever there is a doubly occupied voxel present, all cells that are sufficiently close to this voxel to have a resulting pressure gradient and an empty voxel, or in the case of doubly occupied voxel - a single voxel, next to it, will have a rate to move.

Thus, movement is allowed freely in the framework, which might be the reason why the tumours produced by the framework become asymmetric. A possible explanation to why the results in [2]

are more symmetric might be the fact that when a cell is growing in the Jagiella model and there is no empty neighbouring voxel, it pushes cells along the shortest path towards an empty voxel. This would promote a more symmetric growth, as this would prevent large folds on the tumour rim, like the ones we see in Figure 6, to appear by filling up "holes" on the rim. In the framework, there exists no such prioritisation of path when cells move. The probability that a cell will move into the hole on the boundary is the same as for the cell moving into a voxel next to the hole that might be further away from the growing cell. Further, this might cause a snowball effect, where once cells have started to move and grow on a very small protrusion, the environmental conditions further out in the domain will be more advantageous and promote growth causing the protrusion to grow larger.

Considering this, we propose that the protrusions appearing on the rim of tumours produced by the framework is a natural consequence of the mechanics of the framework. And, since we consider the physics of the framework, where cells move according to pressure gradients, to be a natural way to describe how migration occurs in a population, we suggest that these protrusions are not "wrong". However, it might be that we need to add more mechanisms to the model in order to have the model produce symmetric and finite sized tumours.

8 Future work

In order to test the hypothesis that cells pushed along the shortest path cause a symmetric tumour,

one could perhaps modify this implementation so as to force cells to move along the shortest path,

and see if this makes a difference in symmetry and saturation. Further, one could implement active

migration towards the core as is done in e.g. [19], or a mechanism that causes the tumour to shed

cells on its rim. In [16] they suggest that at some stage of the tumour growth, the thickness of the proliferating rim becomes so thin that all new cells appearing in this zone is lost by either shedding or inwards migration, and thus leads to saturation.

As is evident from results presented in Sections 5.3, 6.1, and 6.2, the placement of boundary conditions have a significant impact on the final size and appearance of the tumour. However, it is desirable to place the boundary conditions far away in vitro, because if the tumour gets too close to the boundary the growth will become too large and not give a valid result, as in the case of Figure 7c. Using Robin boundary conditions,

∂u

∂n = −κ(u − u _D ), I

∂Ω

∂u

∂n dS < M,

would allow us to control the flow of oxygen and glucose into the domain, and restrict it so we only allow the inflow to be some maximum value M .

A modification to this implementation that one might want to explore in order to make it more similar to the Jagiella model, is to let the source function in eq. (8) be incremental, s(u) = [0.1, 0.2, . . . , 0.9, 1], to create a sense of the cell expanding and starting to exert pressure on its surroundings as it progresses through the cell cycle. This might cause proliferating cells to remain in more favourable conditions for a longer time and thus proliferate longer. As the model is im- plemented now, we might "lose" some proliferating cells, as in the Jagiella model cells are pushed towards the rim whenever a neighbouring cell reaches its growth stage, which is m = m g = 2. In this implementation, cells are not pushed until m = m d = 10, meaning proliferating cells might remain in conditions unfavourable for proliferation longer, causing its daughter cells to become quiescent once it divides.

Another modification to consider if one wants to strive for similarity with the Jagiella model is to reflect on the dependency on p(t, x) in eq. (34). In the Jagiella model, this equation depends on the distance to the closest empty voxel, eq. (20), which is linearly increasing away from the considered voxel. The distribution of pressure in the domain is probably not linear however, and one might want to rethink how to scale this variable in this equation. Perhaps something along the line of

exp h

− p(t, x) p _div

i → exp h

− p(t, x) ^1/2 p _div

i would better mimic the dependence in eq. (20).

8.1 Reproducability

These results can be reproduced using the 1.4 release of the URDME simulation framework, which is open-source and can be downloaded at www.urdme.org, and the parameters tabulated in Appendix A.

Acknowledgements

Many thanks to Jan Hasenauer, corresponding author of [2], for his quick response and clarifying

answers. A very sincere thank you to my supervisor Stefan Engblom for his help and suggestions

throughout the course of this work.

References

[1] S. Engblom, D. B. Wilson and R. E. Baker. Scalable population-level modeling of biological cells incorporating mechanics and kinetics in continuous time. Accepted for publication in Roy.

Soc. Open Sci., 2018.

[2] N. Jagiella, D. Rickert, F. J. Theis and J. Hasenauer. Parallelization and high-performance computing enables automated statistical inference of multi-scale models. Cell Systems, 4(2):194206.e9, 2017.

[3] N. Jagiella, B. Müller, M. Müller, I. E. Vignon-Clementel and D. Drasdo. Inferring growth control mechanisms in growing multi-cellular spheroids of NSCLC cells from spatial-temporal image data. PLoS Comput. Biol. 12:e1004412, 2016.

[4] N. Jagiella. Parameterization of lattice-based tumor models from data. PhD thesis, Université Pierre et Marie Curie, Paris, France, 2012.

[5] B. Drawert, S. Engblom, and A. Hellander. URDME: a modular framework for stochastic simulation of reaction-transport processes in complex geometries. BMC. Syst. Biol., 6(76):1–17, 2012.

[6] S. Engblom and others. URDME: Unstructured Reaction-Diffusion Master Equation. Multiple versions exist. Available at www.urdme.org. 2008-2018.

[7] F. Graner and J. A. Glazier. Simulation of biological cell sorting using a two-dimensional extended Potts model. Phys. Rev. Lett., 69(13), 1992.

[8] A. G. Fletcher, M. Osterfield, R. E. Baker, and S. Y. Shvartsman. Vertex Models of Epithelial Morphogenesis. Biophys. J., 106: 2291 – 2304, 2014.

[9] K. A. Rejniak and A. R. Anderson. Hybrid models of tumor growth. Wiley Interdiscip. Rev.

Syst. Biol. Med., 3: 115–125, 2011.

[10] A. G. Fletcher, F. Cooper, and R. E. Baker. Mechanocellular models of epithelial morphogen- esis.Phil. Trans. R. Soc. B, 372: 20150519, 2017

[11] D. Weihs, A. Gefen, and F.J. Vermolen Review on experiment-based two- and three- dimensional models for wound healing. Interface Focus, 6: 20160038, 2016.

[12] J. Delile, M. Herrmann, N. Peyriéras, and R. Doursat. A cell-based computational model of early embryogenesis coupling mechanical behaviour and gene regulation. Nat. Commun., 8:

13929, 2017.

[13] T. Jackson and X. Zheng. A Cell-based Model of Endothelial Cell Migration, Proliferation and Maturation During Corneal Angiogenesis. Bull. Math. Biol., 72(4): 830-868, 2010.

[14] A. Shirinifard, J. A. Glazier, M. Swat, J. S. Gens, F. Family, Y. Jiang, and H. E. Grosskniklaus.

Adhesion Failures Determine the Pattern of Choroidal Neovascularization in the Eye: A Com- puter Simulation Study. PLoS Comput. Biol., 8(5): e1002440, 2012.

[15] T. Roose, S. J. Chapman, and P. K. Maini. Mathematical Models of Avascular Tumor Growth.

SIAM Rev., 49(2): 179-208, 2007.

[16] J. P. Freyer and R. M. Sutherland. Regulation of Growth Saturation and Development of Necrosis in EMT6/Ro Multicellular Spheroids by the Glucose and Oxygen Supply. Cancer Research, 46(7):3504-3512, 1986.

[17] J. Folkman, and M. Hochberg. Self-regulation of growth in three dimensions. J. Exp. Med., 138(4): 745-753, 1973.

[18] D. T. Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys., 22(4): 403-434, 1976.

[19] S. Dormann and A. Deutsch. Modeling of Self-Organized Avascular Tumor Growth with a

Hybrid Cellular Automaton. In Silico Biol., 2(3): 393-406, 2002.

[20] Y. Jiang, J. Pjesivac-Grbovic, C. Cantrell, and J. P. Freyer. A Multiscale Model for Avascular Tumor Growth. Biophys. J., 89(6): 3884-3894, 2005.

[21] K-A. Norton and A. S. Popel. An agent-based model of cancer stem cell initiated avascular tumour growth and metastasis: the effect of seeding frequency and location. J. R. Soc. Interface, 11:20140640, 2014.

[22] M. S. Hutson, M. C. K. Leung, N. C. Baker, R. M. Spencer, and T. B. Knudsen. Computational Model of Secondary Palate Fusion and Disruption. Chem. Res. Toxicol., 30(4): 965-979, 2017.

[23] W. M. Boon, D. C. Koppenol, and F. J. Vermolen. A multi-agent cell-based model for wound

contraction. J. Biomech., 49(8): 1388-1401, 2016.