Scalable population-level modeling of biological cells incorporating mechanics and kinetics in continuous time

TL;DR

This paper introduces a scalable, multiscale computational framework for modeling biological cell populations, integrating mechanics and kinetics to simulate interactions and responses at the tissue level efficiently.

Contribution

The novel framework combines biomechanical properties with continuous-time biochemical signaling, enabling scalable and flexible population-level cell modeling.

Findings

Framework is highly scalable and efficient.

Can be coupled with biochemical signaling models.

Demonstrated flexibility in biological simulations.

Abstract

The processes taking place inside the living cell are now understood to the point where predictive computational models can be used to gain detailed understanding of important biological phenomena. A key challenge is to extrapolate this detailed knowledge of the individual cell to be able to explain at the population level how cells interact and respond with each other and their environment. In particular, the goal is to understand how organisms develop, maintain and repair functional tissues and organs. In this paper we propose a novel computational framework for modeling populations of interacting cells. Our framework incorporates mechanistic, constitutive descriptions of biomechanical properties of the cell population, and uses a coarse graining approach to derive individual rate laws that enable propagation of the population through time. Thanks to its multiscale nature, the resulting simulation algorithm is extremely scalable and highly efficient. As highlighted in our computational examples, the framework is also very flexible and may straightforwardly be coupled with continuous-time descriptions of biochemical signalling within, and between, individual cells.

Equations

$$\frac{\partial u}{\partial t}+\nabla \cdot I$$

$$I$$

$$\epsilon \frac{\partial p}{\partial t}$$

$$-\Delta p$$

$$\displaystyle p\bigr{|}_{\partial\Omega_{D}}$$

$$\displaystyle(\partial p/\partial n)\bigr{|}_{\partial\Omega_{N}}$$

$$\displaystyle\left[\alpha p+(\partial p/\partial n)\right]\bigr{|}_{\partial\Omega_{R}}$$

$$-Lp$$

$$p_i$$

$$I(i\to j)\propto -{\int }_{v_i\cap v_j}\nabla p(x)\text{d}\vec{S}=\frac{e_{ij}}{d_{ij}}(p_i-p_j)=:I_{ij},$$

$$\displaystyle R(i\to j;\;\mbox{$u_{i}\geq 1$, $u_{j}=0$, $v_{j}$ never visited})$$

$$\displaystyle R(i\to j;\;\mbox{$u_{i}\geq 1$, $u_{j}=0$, $v_{j}$ visited previously})\;$$

$$\displaystyle R(i\to j;\;\mbox{$u_{i}>1$, $u_{j}=1$})$$

$$-{\int }_{\partial \Omega }(\nabla p\cdot n)\text{d}\vec{S}$$

$$a_{ijk}$$

$$\delta I_{\text{mbox}adhesion}(i\to j)$$

$${\epsilon }_p\frac{\partial p}{\partial t}$$

$${\epsilon }_S\frac{\partial S}{\partial t}$$

$$-\Delta p$$

$$0$$

$$S(x_1,x_2)\to 0,$$

$$\frac{\partial S}{\partial x_1}\to 0,$$

$$S(x_1,x_2)=\frac{Q}{2\sqrt{k{D}_S}}{e}^{-\sqrt{k/D_S}|x_1-X_s|}.$$

$$I(i\to j)=-{\chi }_2{\int }_{v_i\cap v_j}\nabla S(\vec{x})\text{d}\vec{\mathrm{S}}={\chi }_2\frac{e_{ij}}{d_{ij}}(S_i-S_j),$$

$$\displaystyle\begin{array}[]{rcl}n^{\prime}_{i}&=&f(\bar{d}_{i})-n_{i},\\ d^{\prime}_{i}&=&v\left(g(n_{i})-d_{i}\right),\\ \end{array}$$

$${\bar{d}}_i$$

$$-Lc$$

$$c_i$$

$$D_1=0.01\text{leavevmode}{\text{s}}^{-1},D_2=25\text{leavevmode}{\text{s}}^{-1},D_3=0.01\text{leavevmode}{\text{s}}^{-1},\lambda =0.0015\text{leavevmode}{\text{s}}^{-1},$$

$${\kappa }_{\text{mbox}prol}=0.65,{\rho }_{\text{mbox}prol}=0.125\text{leavevmode}{\text{s}}^{-1},{\kappa }_{\text{mbox}death}=0.55,{\rho }_{\text{mbox}death}=0.125\text{leavevmode}{\text{s}}^{-1},{\rho }_{\text{mbox}deg}=0.01\text{leavevmode}{\text{s}}^{-1},$$

Full text

In this paper we propose a novel computational framework for modeling populations of interacting cells. Our framework incorporates mechanistic, constitutive descriptions of biomechanical properties of the cell population, and uses a coarse graining approach to derive individual rate laws that enable propagation of the population through time. Thanks to its multiscale nature, the resulting simulation algorithm is extremely scalable and highly efficient. As highlighted in our computational examples, the framework is also very flexible and may straightforwardly be coupled with continuous-time descriptions of biochemical signalling within, and between, individual cells.

Keywords: Continuous-time Markov chain, Computational cell biology, Cell population modeling, Notch signalling pathway, Avascular tumour model.

AMS subject classification: 60J28, 92-08, 65C40.

1 Introduction

Development, disease and repair all require the tightly coordinated action of populations of cells. In almost every case a plethora of interacting components, acting on a range of spatial and temporal scales, combines to drive the observed tissue-level behaviours. Research in the biosciences is now advanced to the stage where we have sufficient understanding of intercellular processes to build relatively sophisticated models of a wide range of cellular behaviours. A promising approach to generate, test and refine hypotheses as to the relative contributions of various mechanisms to tissue-level behaviours is that of cell-based computational modeling: individual cells are explicitly represented, each cell has a position that updates over time, and cells may also have an internal state or program determining their behaviour. Cell-based models hold great potential in this regard because they can naturally capture both stochastic effects and cell-cell heterogeneity, and they can be used to explore tissue-level behaviours when complex hypotheses on the cellular scale prevent straightforward continuum approximations at the tissue level. Recent applications of cell-based models to study population-level behaviours include embryonic development [13, 29, 52, 32, 4], wound healing [57, 56, 63] and tumour growth [3, 1, 39, 40, 33].

Multiple cell-based modeling approaches exist, and they can be categorised according to their approaches to representing cell positions as being either on- and off-lattice. In the on-lattice approach, space is divided up into a discrete grid of lattice sites. A common type of on-lattice model is the cellular automaton, in which each cell occupies a single lattice site and attempts to move to a new site at each time step according to a set of update rules. This volume exclusion rule can be relaxed to allow multiple cells per lattice site, depending on the level of spatial description required by the problem under consideration. Position update rules typically take into account the number and type of neighbouring cells, and can also depend on other information associated with lattice sites, such as nutrient or signalling factor concentrations [20]. Hybrid models often use systems of ordinary differential equations (ODEs) or partial differential equation (PDEs) to model the evolution of biochemical concentrations (see, for example, [1, 39, 41, 45]). Position update rules can also be stochastic, so that cells move according to, e.g., a biased random walk [36, 37, 5, 61, 44, 50]. In most cases a regular lattice is used, for example square or hexagonal, but unstructured lattices have also been used in several cases (see [26, 38] and references therein). A different on-lattice technique is the cellular Potts model, wherein each cell is allowed to occupy multiple lattice sites, and energy minimisation is used to propagate the shape of each cell over time. The cellular Potts model has been used to study biological processes ranging from cell sorting [17, 62] and morphogenesis [55, 21] to tumour growth [47, 43], and methods to provide a macroscopic limit to these models have been provided [2, 28].

A basic off-lattice approach to cell-based modeling is the cell-centre-based model which assumes cells are, in effect, point particles that interact with each other via some specified potential function [11, 31, 39]. Meanwhile, in vertex models [27, 58, 14] cell populations are modelled as a tessellation of polygons or polyhedra, whose vertices move due to forces originating from the cells, whilst in immersed boundary models [9, 40, 48, 49] cell boundaries are represented as a set of points that move like elastic membranes immersed in a fluid. Other cell-based models take the subcellular composition of cells into account, such as the subcellular element model [34] or the finite element vertex model [6].

Each cell-based modeling approach has specific advantages and disadvantages, and incorporate different representations of both biochemical signalling and biomechanics [12, 45, 48, 23, 54, 35]. For example, an immersed boundary model allows a detailed representation of cell shapes, but this benefit comes at an increased computational cost in comparison to other methods. Cellular Potts models are very versatile in the possible effects that can be modelled but, because updates are controlled via a Metropolis-type algorithm together with a user-specified Hamiltonian, time can only be measured in terms of Monte Carlo steps. When modeling a specific application, it is necessary to weigh the benefits of the existing cell-based models against each other in the context of the specific application. A comprehensive review of mechanocellular models, including a table that compares the different modelling approaches, is provided in [12, 54, 15] and a comparison of different model output is detailed in [54, 35].

In this work we will focus on lattice-based approaches that allow a user-specified maximum number of cells to occupy each lattice site [26, 10, 20]. On a given structured or unstructured tessellation of space we develop constitutive equations governing the dynamics of the cell population. As in the cellular Potts model, our update rules are stochastic and are established from global calculations. An important difference, however, is that our simulations take place in continuous time, thus allowing for a meaningful coupling to other continuous-time models. Our rationale for developing this approach is a desire to be able to quickly and efficiently simulate three-dimensional tissues consisting of large numbers of cells, with inclusion of biomechanical effects, intercellular signalling and external inputs. Since our approach is to base the modeling of cell biomechanics on the Laplace operator over the spatial tessellation, we refer to our method as Discrete Laplacian Cell Mechanics (DLCM). Relying on the discrete Laplace operator is advantageous because one may make use of highly developed and scalable numerical methods to evolve the biomechanical details of the cell population over time. An important advantage provided by our highly efficient framework is then the potential ability to conduct parameter sensitivity analysis, parameter inference and model selection for cell-based models; these are increasingly important research tools in this era of quantitative, interdisciplinary biology.

The outline of this work is as follows: in §2 we detail the model and its computational implementation; in §3 we describe four examples that showcase the utility of our modeling framework; and in §4 we summarise our results and discuss avenues for future exploration.

2 Methods

The method we propose is developed by distributing the cells onto a grid of voxels and defining a suitable physics over this discrete space. The Laplace operator emerges as a convenient and basic choice to describe evolution of the biomechanics of the population, but more involved alternatives could also be employed in its place. We enforce a bound on the number of cells per voxel such that processes at the scale of individual cells may be meaningfully described on a voxel-local basis. For the simulations performed in this paper the voxels contain a maximum of two cells, but larger carrying capacities than this can also be supported. The choice of discretisation (and so the maximum number of cells that can be accommodated in any voxel) should be made on a case-by-case basis, taking into account the need to balance computational complexity with the extent to which data on individual-cell-level processes are available. By evolving the individual cells via discrete PDE operators, e.g. the discrete Laplacian, processes at the population level are connected in an efficient and scalable way to those taking place inside the individual cells. In §2.1 we offer an intuitive algorithmic description of our framework, and a more formal development is found in §2.2.

2.1 Informal overview of the modeling framework

We consider a computational grid consisting of voxels $v_{i}$, $i=1,\ldots,N_{\mbox{{\tiny vox}}}$, in two or three spatial dimensions. We make no assumptions as to whether the grid is structured or not, however, we do require that a consistent Laplace operator may be formed over the grid. Each voxel $v_{i}$ shares an edge with a neighbour set $N_{i}$ of other voxels. In two dimensions, each voxel in a Cartesian grid has four neighbours and on a regular hexagonal lattice, each voxel has six neighbours. On a general unstructured triangulation, each vertex of the grid has a varying number of neighbour vertices and, in this general and flexible case, the voxels themselves can be constructed as the polygonal compartments of the corresponding dual Voronoi diagram (Figure 2.1).

At any given point in time the voxels are either empty or may contain a certain number of cells. If the number of cells is at or below the carrying capacity, the system is assumed to be at steady-state. Hence in the absence of any other processes, the cell population is then completely static. If the number of cells in one or more voxels exceeds the carrying capacity, the cells push each other and exert a cellular pressure. Eventually, this non-equilibrium state is changed by an event, for example one of the cells moves into a neighbouring voxel, and the pressure is redistributed. This goes on until, possibly, the system relaxes into a steady-state.

What is the relevant constitutive equation for this cellular pressure? We make a more detailed investigation of this in §2.2 below, but chiefly, for an isotropic medium and a scalar potential, thus essentially assuming the pressure to be spread evenly as in Figure 1(b), the answer is that the pressure is distributed according to the negative Laplacian, with source terms for all voxels where the carrying capacity is exceeded.

At any instant in time, a pressure gradient between two neighbouring voxels induces a force that, in turn, can cause the cells within the voxels to move. The rate of this movement is proportional to the pressure gradient, with a conversion factor that may depend on the nature of the associated movement. For example, it may be reasonable to assume that a cell may move into an empty voxel or into an already occupied voxel with different rates per unit of pressure gradient.

The simulation method is event-based and takes the form of an outer loop over successive events, see Algorithm 1, lines 2–11. Because the dynamics is driven by a discrete numerical PDE operator, i.e. the Laplacian in our context, we call the simulation method Discrete Laplacian Cell Mechanics (DLCM). It should be noted that it might be beneficial in certain situations to use a different driving PDE operator: for a concrete example, see §3.2.

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 [16] 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 summary, our new modelling framework is lattice-based, allows for a range of physical and biomechanical phenomena to be accurately described, and facilitates explicit modelling of cell size/excluded volume effects through the pressure-dependence of cell movements. The incorporation of stochasticity via the Gillespie algorithm renders it naturally able to represent the noise observed in most biological systems. Moreover, it is simple to simultaneously integrate either deterministic or stochastic models of biochemical reaction networks. Through the choice of sensible numerical linear algebra techniques, the framework is efficient and even complex models can be simulated in short times. This means that it will be possible to perform statistical inference of model parameters using quantitative experimental data together with techniques such as approximate Bayesian computation [51, 25].

2.2 Formal description

We now develop the details of the DLCM framework. At some instant in time, let the grid $(v_{i})$, $i=1,\ldots,N_{\mbox{{\tiny vox}}}$, be populated with $N_{\mbox{{\tiny cells}}}$ cells. In the interest of a transparent presentation we only allow each voxel to be populated with $u_{i}\in\{0,1,2\}$ cells; other arrangements may also be useful, but a suitable physics for voxels populated below the carrying capacity should then depend on biological details such as the tendency of the cells to stay in close proximity to each other.

Due to the spatial discretisation and the discrete counting of cells, the task is to track changes over this chosen state space. In continuous time this amounts to figuring out which cell will move to what voxel, and when it will move. This requires a governing physics defined over the discrete state. A continuous-time Markov chain respects the “memoryless” Markov property and stand out as a promising approach, requiring only movement rates in order to be fully defined. Our model of the population of cells follows from three equations (2.1), (2.2), and (2.3) below, understood and simplified under three assumptions, Assumptions 2.1–2.3. We present each in turn as follows.

Let $u=u(t,x)$ represent the cell density at time $t$ and at the point $x$. We need to keep in mind that our numerical model is to be formulated on an existing grid of voxels $(v_{i})$ containing a bounded (integer) number of cells $u_{i}$. Hence the continuum limit $h\to 0$ (of some suitable measure of the voxel size going to zero) is not meaningful. However, we shall use a continuous notation initially for ease of presentation. Our starting point is the continuity equation

[TABLE]

where $I$ is the current, or flux. Since we are aiming at an event-based simulation we will later use eq. (2.1) to derive rates for discrete events in a continuous-time Markov chain. To prescribe the current $I$ we now make a starting assumption.

Assumption 2.1**.**

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

Assumption 2.1 expresses the idea that small Brownian-type movements of each cell about its (voxel-) center can be ignored. It does not exclude an additional description of any active movements, such as chemotaxis or haptotaxis. With sufficient conditions for equilibrium specified, it follows from Assumption 2.1 that only doubly occupied voxels will give rise to a rate to move, and we will describe this increased rate as a pressure source. In the absence of any other units we can set this pressure source to unity identically.

Let $p=p(t,x)$ denote the cellular pressure at time $t$ and position $x$, again using a continuous notation for variables which will be implemented on a discrete grid. Interpreting the current $I$ as the result of a pressure gradient, we take the simple phenomenological model

[TABLE]

which, with the interpretation $D\equiv\kappa/\mu$, the quotient between the permeability $\kappa$ and the viscosity $\mu$, is a form of Darcy’s law in porous fluid flow. Notably, Darcy’s law can be obtained from first principles using homogenisation arguments [59]. Our use of eq. (2.2), however, will be different as cells preserve their geometrical integrity and do not flow freely like a fluid. Rather, the current in eq. (2.2) will be understood as a rate of the discrete event that a cell changes voxel.

To complete the model we require a constitutive equation relating $u$ and $p$. With the cellular pressure driven by sources in the form of overcrowded voxels, we provide a constitutive model of pressure evolution using the heat equation,

[TABLE]

with $s(u)$ a source function that will be prescribed below. Specifically, eq. (2.3) follows by assuming an isotropic medium and a scalar potential pressure. For each voxel 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.

Eq. (2.3) is a time-dependent PDE and would be complicated to handle within the current context. To move forward we therefore need to bring in 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.

In cases where the biochemical kinetics of the individual cells also affect their mechanical behaviour, for example via signal transduction or proliferation, Assumption 2.2 entails that these processes must occur on a slower time-scale than the propagation of the cellular pressure. Importantly, Assumption 2.2 simplifies eq. (2.3) into the non-singular $\varepsilon\to 0$ limit, the Laplacian equilibrium,

[TABLE]

The domain $\Omega$ understood here generally consists of the bounded subset of $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$ which is populated by the cells. Its boundary $\partial\Omega$ can be written as $\partial\Omega=\partial\Omega_{D}\cup\partial\Omega_{N}$, the Dirichlet and the Neumann boundaries, respectively, which can be chosen according to the specific biological problem under consideration. It is also possible to interpolate between the two boundary conditions,

[TABLE]

to model, for example, an increasingly impenetrable cellular matrix as $\alpha\to 0$ by a homogeneous Robin condition. To simplify the presentation here, we employ homogeneous Dirichlet conditions ($\partial\Omega=\partial\Omega_{D}$) throughout the computational examples in §3.

We now re-interpret the developed model onto the given tessellation $(v_{i})$, aiming specifically for a time-continuous and event-driven simulation. Denote by $\Omega_{h}$ the subset of voxels $v_{i}$ for which $u_{i}\not=0$. Similarly, let $\partial\Omega_{h}$ denote the discrete boundary; this is the set of unpopulated voxels that are connected (i.e. share an edge) with a voxel in $\Omega_{h}$. See Figure 2.2 for an illustration. We first consider the discrete version of eq. (2.4),

[TABLE]

where the source term is given by, in non-dimensionalized form, $s(u_{i})=0$ for $u_{i}\leq 1$ and $s(u_{i})=1$ whenever $u_{i}=2$, in view of the normalisation $p=0$ for the static case, following Assumption 2.1. In eq. (2.8), $L$ is a discrete Laplacian over the currently active grid $\Omega_{h}$. The precise choice of (consistent) numerical method chosen to define $L$ is not likely to have a strong influence on the model output since the Laplacian operator is quite forgiving to such details. In our experiments we used a finite-element-based discrete Laplacian operator 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; in continuous language, $s(u)=s(u(x))$ is just $\sum_{u_{i}>c}\delta_{v_{i}}(x)$, where $v_{i}$ here denotes the center point of the $i$th voxel and $c$ the carrying capacity.

We now go back to considering the induced movement of the cells according to the continuity equation (2.1). By eq. (2.2), a pressure gradient gives rise to a proportional current; $I\propto-\nabla p$. Let us write $I(i\to j)$ to denote the current from voxel $v_{i}$ to neighbour voxel $v_{j}$. In our discrete setting this current is calculated by integrating the pressure gradient across the edge between the two voxels,

[TABLE]

with $d_{ij}$ the distance between voxel centres, $e_{ij}$ the shared edge length, and where the simplification takes the discrete nature of the model into account by simply substituting a scaled difference for the gradient. In eq. (2.10), the current is reversed $(j\to i)$ when the result is negative. It follows that, in the presence of a pressure source, a non-zero pressure gradient results everywhere except on solid (Neumann) boundaries. It remains to specify $D$ in eq. (2.2), including the subset of cells that are free to move as a result of this pressure gradient.

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 us write $R(e)$ for the rate of the event $e$ and use the notation $i\to j$ for the particular event that one cell moves from voxel $i$ to $j$. Under our present framework, $u_{i}$ is constrained to $\{0,1,2\}$ and a total of three distinct cases emerges for the movement rates:

[TABLE]

Here, $D_{1}$ is the conversion factor from a unit pressure gradient into a movement rate to a voxel $v_{j}$ never visited before, $D_{2}$ similarly, but into a voxel previously occupied, and $D_{3}$ covers the “crowding” case where a cell in a doubly occupied voxel enters a singly occupied voxel.

Although dependent on the specific application under consideration, we generally expect to have $D_{1}\ll D_{2}$, since the extracellular matrix contained in previously unvisited voxels is usually thought to be less penetrable than that of previously occupied voxels. To understand the scale of $D_{3}$ we note that from the divergence theorem

[TABLE]

In other words, the total pressure gradient over any closed surface $\partial\Omega$ is equal to the enclosed sources. With $\partial\Omega$ the boundary surrounding the whole cell population, it follows that the ratio $D_{2}/D_{3}$ expresses the preference for events at the tissue boundaries (cells entering empty voxels) to events internal to the region (cells in doubly occupied voxels moving to an already occupied neighbouring voxel). For example, assuming the presence of a single internal source voxel and that $D_{2}/D_{3}\equiv 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 $\partial\Omega$. The rates in eqs. (2.11)–(2.13) are illustrated on a common scale with $D_{1}=D_{2}=D_{3}$ in Figure 2.3. Note that, on this Cartesian grid, by the integration over an edge in eq. (2.10), the current is non-zero only along the cardinal directions up/down and right/left, respectively.

In addition to the rates in eqs. (2.11)–(2.13), one can easily extend the basic framework to include active cell movements and/or behaviours, as prescribed via additional rates and to be handled in lines 3–4 of Algorithm 1. For example, a possible approach to include the effects of cell-cell adhesion is to define the resistance force for a cell moving from voxel $i$ to $j$ due to cell-cell adhesion with cells in the non-empty voxel $k$ as

[TABLE]

for some adhesion constant $\alpha$ and for neighbor voxel pairs $(i,j)$ and $(i,k)$. Here $e_{ik}$ is the edge length over which adhesion bonds are active and $(\vec{r}_{ij},\vec{r}_{ik})$ denote unit vectors pointing from voxel $i$ to $j$ and, respectively, from voxel $i$ to $k$. The $\min$-construct in eq. (2.15) ensures that adhesion acts as a resistance (i.e. negative) force. Summing up the adhesion forces due to all populated neighbor voxels of voxel $i$ and scaling by the inter-voxel distance $d_{ij}$ we thus evaluate the net adhesion current contribution as

[TABLE]

which is to be added to the current in eq. (2.10). The resistance force works against the pressure-driven current in eq. (2.10) to reduce cell movements.

Other options to capture similar effects could be to modify the pressure source function and/or the boundary conditions; the best choice is clearly dependent on the details of the modeling situation.

3 Results

We now present some computational results for the proposed modeling framework. In doing so we focus on highlighting the three distinct advantages of our approach: (1) the mechanics is based on constitutive equations obeying a well-understood physics; (2) the framework is fast and scalable, and therefore convenient to experiment with; and (3), it is also flexible and allows for a seamless coupling to other processes of interest.

Throughout the examples below we constructed our computational grid over a rectangular base geometry; we used structured Cartesian and hexagonal grids in 2D and in one case a Cartesian grid in 3D. The grid was initially populated by cells in the middle of the domain, cf. line 1 of Algorithm 1. As previously mentioned, we relied upon a finite-element-based discrete Laplacian operator defined using linear basis functions and a lumped mass-matrix, implemented in the Matlab add-on package PDE Toolbox [30]. The discrete Laplacian was inverted using Matlab’s built-in LU-decomposition (line 3 of Algorithm 1); this was fast enough for our purposes and in all tests performed here. Larger models may well benefit from more advanced linear solver techniques. All code and the required data is available, and we refer the reader to §4.1 for details.

In §3.1 we look at the possible dependency of results on the underlying computational grid, and in §3.2 we investigate how to couple the cellular behaviour with signals in the local external environment. In §3.3 we couple the modeling framework with intracellular continuous-time processes, including cell-to-cell signalling and, finally, in §3.4 we demonstrate how the complex behaviours of a tumour model may be investigated.

3.1 Simulated cell population is free of grid artefacts

A disadvantage with grid-based methods and local update rules that has been highlighted in the literature is that grid artefacts may appear unless some care is exercised [60]. We test for grid artefacts by artificially forcing a number of cells into a square configuration and then relaxing the system to equilibrium (Figure 3.1). In the absence of any other mechanical forces, clearly, the equilibrium state should be approximately circular.

In Figure 3.2 we quantify the average cellular density across $N=100$ independent runs of the model and find that, although a very small memory effect from the initial data is still visible (left/right and up/down), there is no dependence on the underlying grid except for the natural distortions due to discreteness. This comes from the fact that the physics is based on well-understood constitutive equations and that we employ a grid-consistent discretisation of the Laplacian (eqs. (2.4) and (2.8)).

3.2 Population-level behaviour via sensatory control

A notable extension of this basic framework is to model a population of cells that respond to external stimuli in their local environment. An advantage of the DLCM framework is that we are not restricted to considering a Laplacian operator but can freely choose other PDE operators, depending on the mechanics of the problem under consideration. As a test case we take the two-dimensional migratory model of neurons responding to an inhibiting chemical known as Slit [7]. The cellular pressure, $p$, of the neurons is now described by the chemotaxis equation

[TABLE]

where $\chi_{1}$ is the chemotactic sensitivity of the neurons to the chemical $S$, i.e., the concentration of Slit. The parameter $\varepsilon_{p}$ is taken to be small and represents our assumption that the pressure equilibrates on a faster timescale than the movement of the neurons. The source term, $s(u)$, is as defined previously. Eq. (3.1) represents the movement of cells down the gradients of the Slit concentration. We set the domain for the experiment to be $\vec{x}=(x_{1},x_{2})\in\mathbb{R}^{2}$ and the equation governing the concentration of Slit is given by

[TABLE]

where $D_{S}$ is the constant diffusivity of Slit, $k$ is the rate of degradation and $Q$ is the strength of the Slit source on the line $x_{1}=X_{s}$. The parameter $\varepsilon_{S}$ is again taken to be small and represents the assumption that the diffusion of Slit occurs on a faster timescale than the movement of the cells. Imposing Assumption 2.2 thus allows us to solve the steady-state problem

[TABLE]

where the boundary conditions for the Slit chemical are

[TABLE]

We note that eqs. (3.4)–(3.6) can be solved analytically to give

[TABLE]

We also include a pressure-independent movement mechanism, representing the active motion of the cells due to the chemo-repellent (see Assumption 2.1). To include the chemotaxis-driven movement of a cell from voxel $i$ to voxel $j$ we derive the following current

[TABLE]

where $\chi_{2}$ is a measure of the affinity of the cells to Slit for this type of pressure-independent movement. We initialise a circular region, with radius $r=10$, of doubly occupied sites in the centre of the domain, where the voxel size is taken to be $h=1$. We specify the parameters $[\chi_{1},\chi_{2},D_{S},k,Q,X_{s}]=[100,5,50,0.1,2,50]$, and run $100$ simulations until terminal time $T=1000$. We present the average cell density profile in Figure 3.3. As expected, on average the cell population moves away from the source of Slit, consistent with it being a chemo-repellent.

3.3 Continuous-time mechanics allows for seamless coupling to

cellular signalling processes

A great flexibility with the proposed framework comes from the fact that additional processes may be simulated concurrently with the overall mechanics of the cell population (see Algorithm 1, line 8). Notably, such processes may include cell-to-cell communication driving cell fate differentiation. To illustrate this we consider the classical Delta-Notch intracellular signalling model [8], which produces pattern formation via a simple lateral inhibition feedback loop. With $(n_{i},d_{i})$, respectively, the Notch and Delta concentrations within cell $i$, this model takes the dimensionless form

[TABLE]

We take parameters $[a,b,v,k,h]=[0.01,100,1,2,2]$ and the average in eq. (3.12) is taken over the set of neighbours $N_{i}$ of the cell $i$. Additional stochastic noise terms may also be added to eq. (3.11), but for simplicity we let the model be fully deterministic.

To produce a dynamic population we let the cells in the middle of the region proliferate at a constant rate and we terminate the simulation when $N_{\mbox{{\tiny cells}}}=1000$. The average in eq. (3.12) is appropriately modified to account for any doubly occupied voxels. To get a more realistic contact pattern a hexagonal grid was used; hence each voxel has six neighbours.

In Figure 3.4 a typical time-series of this model is summarised. Here the time-scale of the Delta-Notch dynamics, eq. (3.11), is on a par with that of the proliferation process, so that their dynamics equilibrates after the tissue stops growing. A second simulation is summarised in Figure 3.5, where the right-hand side of eq. (3.11) has been scaled by a factor of 50 so that the Delta-Notch model is in quasi-steady-state as the tissue grows. The significant differences between the patterns could potentially be used, together with experimental data, to better understand how the timescales for patterning compare with those of tissue growth.

3.4 Non-trivial dynamics emerges from the combination of grid-based physics and local behaviour

Finally, we utilise our proposed framework to build a model for avascular tumour growth. It is known that an important factor in the growth of tumours is the availability of oxygen to the tumour. In the case of an avascular tumour there is no direct supply of oxygen to the tumour itself, but instead oxygen in the surrounding tissue diffuses into the tumour. We consider the tumour to be a growing population of cells that can potentially proliferate and/or die over the course of the simulation. We extend the occupancy of voxels such that $u_{i}\in\{-1,0,1,2\}$, where $u_{i}=-1$ corresponds to a dead cell occupying a voxel. The cells at the boundary of the tumour can consume oxygen and proliferate if the concentration of oxygen is sufficient. Further into the tumour, where the concentration of oxygen is generally lower, the cells continue to consume oxygen but can no longer proliferate; these are known as quiescent cells. Near the centre of the tumour the oxygen concentration is low enough that the cells die and eventually degrade; this region is known as the necrotic core.

As before we solve the discrete Laplace equation (2.8) for the pressure in each voxel, with homogeneous Dirichlet boundary conditions (2.9) on the boundary $\partial\Omega_{h}$, the collection of empty voxels that are adjacent to occupied voxels. We also have an equation for the oxygen concentration, $c$, which diffuses through the domain with sources at the fixed external boundary, $\partial\Omega_{\mbox{{\tiny ext}}}$, thus,

[TABLE]

where $\lambda$ is the rate of consumption of oxygen for a single cell and $a(u_{i})$ is the number of alive cells ($u_{i}\in\{0,1,2\}$) in the $i$th voxel (doubly occupied voxels consume twice as much oxygen and empty voxels or those containing dead cells consume no oxygen). From these discretised equations we can calculate all the different rates for the possible events. These events are as follows. A cell occupying its own voxel, $u_{i}=1$, will proliferate at rate $\rho_{\mbox{{\tiny prol}}}$ if $c_{i}>\kappa_{\mbox{{\tiny prol}}}$, where $\kappa_{\mbox{{\tiny prol}}}$ is the minimum oxygen concentration for proliferation to occur. A living cell will die at a rate $\rho_{\mbox{{\tiny death}}}$ if $c_{i}<\kappa_{\mbox{{\tiny death}}}$, where $\kappa_{\mbox{{\tiny 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 $\rho_{\mbox{{\tiny deg}}}$ to free up the voxel ($u_{i}=0$) for other cells to move in to it. We also include the rates for cell movement as shown in eqs. (2.11)–(2.13).

In Figure 3.6 we present snapshots from a realisation of the model with parameters

[TABLE]

[TABLE]

that demonstrates the evolution of the tumour into the classical regions of proliferating cells, quiescent cells and a necrotic core. We also present the evolution of the numbers of each cell type over the duration of the realisation (Figure 3.6 bottom right). The widely-held view in the field is that a limited oxygen supply leads to a stable, finite-sized tumour [42, 22, 18]. Recent models for avascular tumour growth have been utilised to successfully predict the growth of tumour cell lines that exhibit sigmoidal growth curves [19]. Sigmoidal curves have three phases: for an initially small population of cells (Figure 3.6 (a)) in an abundance of oxygen we see exponential growth (Figure 3.6 (b)), which is then followed by a quasi-linear growth of proliferating cells (Figure 3.6 (c)). Eventually the lack of oxygen causes three heterogeneous cell types to appear (Figure 3.6 (d)); a proliferating ring, a quiescent annulus and a necrotic core. It is here where the volume of the avascular tumour should begin to plateau, as observed in simulations (Figure 3.6 bottom right (d)).

However, due to a possible lack of surface adhesion in our model we observe a fourth phase emerge, where asymmetric protrusions appear on the surface (Figure 3.6 (e)) of the tumour which increases the total uptake of oxygen and hence a further incline in the total cell number (Figure 3.6 bottom right). We sought a parameterisation of our model that did not exhibit this fourth phase and instead produces a steady-state tumour, however we could not do so despite a systematic search of parameter space. We can navigate the parameter space as follows; firstly, select a user-specified finite size tumour radius and with it solve the steady-state oxygen equation. Then the parameters ($\lambda,\kappa_{\text{prol}},\kappa_{\text{death}}$) can be chosen to decide where the proliferating ring and necrotic core should lie. Next, initialise a population of cells at the required radius with some proliferation and death rates $\rho_{\text{prol}}$ and $\rho_{\text{death}}$, and temporarily let $\rho_{\text{deg}}=\infty$ such that dead cells instantly degrade. We can then alter the movement rates $D_{1}$, $D_{2}$ and $D_{3}$ such that the number of cells proliferating is roughly equal to the number of cells being instantly absorbed in the necrotic core; this means that, on average, the size of the tumour will be constant. Finally, returning to the full avascular tumour model, we can vary the remaining parameters ($\rho_{\text{prol}},\rho_{\text{death}},\rho_{\text{deg}}$) in an attempt to build a finite-size tumour. Despite this comprehensive exploration of the parameter space, which is only possible due to the efficiency and versatility of the DLCM framework, at best we could only slow down the growth of the asymmetric protrusions by increasing $D_{2}$ to be far greater than $D_{1}$. That cells invade new matrix at a slower rate appears to be a reasonable assumption here, but the exact separation in scales between $D_{1}$ and $D_{2}$ will need to be investigated on a case-by-case basis using cell trajectory data collected from relevant experiments. Thus the current set of modelling assumptions, which only use passive physics, are not capable of developing the avascular tumour spheroids seen in experiments. Our results are in line with previous works [24, 10] which both conclude that additional mechanisms are necessary in order for tumour growth to be stabilised. The suggested mechanisms involve the necrotic cells secreting chemicals that can either act as a growth inhibitor for healthy cells to become quiescent [24], or as a chemotactic signalling chemical to attract tumour cells to migrate back towards the necrotic core [10]. The inclusion of such mechanisms fits naturally into our proposed DLCM framework and will be implemented in future, more focussed studies.

4 Discussion

The goal of this work was to provide an efficient, hybrid and multiscale computational framework for modelling populations of cells at the tissue scale. Our framework can harness either a regular or unstructured lattice, as required by the biological problem under consideration, and it draws upon a constitutive, mechanistic description of cellular biomechanics to drive expansion and/or retraction of the cell population as cells move and undergo proliferation and death. By assuming that the timescale on which the tissue relaxes to mechanical equilibrium is much shorter than that of other mechanical processes, we are able to assume that the “cellular pressure” within the population is governed by the Laplace equation with specified boundary conditions and source terms. We then use the pressure gradient within the tissue to derive rate equations for the movement of cells within the population. This enables us to propagate the model in continuous time, using an event-driven algorithm such as that originally proposed by Gillespie [16]. A significant advantage of our approach in this regard is that both discrete and continuous models of biochemical signalling can easily and flexibly be incorporated into the framework. This enables the user to develop hybrid and multiscale models that include constitutive descriptions of cellular biomechanics, inter- and intra-cellular signalling in an efficient and scalable framework. Our approach is flexible, computationally efficient and provides a consistent description of cellular mechanics. However, it does not, in the form presented here, allow for detailed tracking of individual cells, or a resolution of their shape, over time; other cell-based modelling approaches may be more suitable if such a spatial resolution is required [54, 15].

We have demonstrated the utility of our approach using four examples. First, we demonstrated that, as expected, a simple colony of proliferating cells relaxes isotropically. Second, we demonstrated how to integrate signals from the local microenvironment into cell movement laws, using the migration of cells within a neuronal explant exposed to a gradient in the secreted protein Slit as an example. Our third example showed that it is easy to integrate cell-cell signalling models into the framework, using the delta-notch lateral inhibition model as a test case. Finally, we developed a model for avascular tumour growth, wherein cell proliferation and death is controlled by the local oxygen concentration, and the cell population generates an oxygen profile within the tumour as the cells consume oxygen.

In summary, the real advance provided by our approach is the ability to quickly and efficiently simulate the behaviour of large populations of cells, in both two and three spatial dimensions. Our model is both hybrid and multiscale in nature; it can incorporate both continuum and discrete models of biochemical signalling both within and between cells. For the purpose of experimenting we have developed a serial Matlab implementation of the model that employs a direct factorisation of the discrete Laplace operator and, as such, scales conveniently to about 20,000 cells. Using a compiled language and parallelisation one could likely improve upon this figure to some extent. However real savings in computing time will require optimal Laplace solvers that utilise, for example, algebraic or geometric multigrid techniques [46, 53] that can take advantage of the incremental nature of the computational process when the domain evolves slowly. We leave this aspect of the implementation for future work.

In the modern era of biology, where quantitative data describing the evolution of cell populations and tissues can be collected with relative ease, we are now in a position to test and validate experimentally generated hypotheses using biologically realistic mechanistic models. However, these models need to be calibrated against the available data, and the sensitivity of model predictions to changes in model parameters needs to be explored. All of this requires repeated simulation of multiscale and hybrid cell-based models, the computational demands of which have to-date provided a barrier to significant progress. The DLCM method we outline here results in significant advances in our ability to efficiently simulate hybrid and multiscale cell-based models in two and three spatial dimensions, and therefore provides a very real opportunity to test and validate mechanistic models using quantitative data.

4.1 Availability and reproducibility

The computational results can be reproduced within release 1.4 of the URDME open-source simulation framework, available for download at www.urdme.org.

Acknowledgement

REB would like to thank the Royal Society for Wolfson Research Merit Award and the Leverhulme Trust for a Leverhulme Research Fellowship. DBW would like to thank the UK’s Engineering and Physical Sciences Research Council (EPSRC) for funding through a studentship at the Systems Biology programme of The University of Oxford’s Doctoral Training Centre.