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

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


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 2018 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

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 discretization (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.

Informal overview of the modelling framework
We consider a computational grid consisting of voxels v i , i = 1, . . . , N 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 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, but chiefly for an isotropic medium and a scalar potential, thus essentially assuming the pressure to be spread evenly as in figure 1b, 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 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 [49]  (b) (a) Figure 1. Schematic explanation of the numerical model. An unstructured Voronoi tessellation (a) with green voxels containing single cells and a red voxel containing two cells. The modelling physics for the cellular pressure can be thought of as if the pressure was spread evenly via linear springs connecting the voxel centres (b). The grid here is unstructured, but the same derivation is used for regular, e.g. Cartesian and hexagonal, grids. 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 manner (line 8). Finally, as the event is processed, a new cell population state is obtained and the loop starts anew.

Algorithm 1
Outline of the DLCM simulation methodology.
Solve for the cellular pressure (p i ), equations (2.8) and (2.9). Compute all movement rates (r j ) for the subset of voxels where the cells may move, equations (2.11)-(2.13).

4:
Determine also the rates for any other processes taking place in the model such as proliferation 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 n−1 j=1 r j < U 2 λ ≤ 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].
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 [50,51].

Formal description
We now develop the details of the DLCM framework. At some instant in time, let the grid (v i ), i = 1, . . . , N vox , be populated with N cells cells. In the interest of a transparent presentation, we only allow each voxel to be populated with u i ∈ {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.
Owing to the spatial discretization 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 stands 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.3), 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 → 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 where I is the current, or flux. Since we are aiming at an event-based simulation we will later use equation (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-) centre 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 which, with the interpretation D ≡ κ/μ, the quotient between the permeability κ and the viscosity μ, is a form of Darcy's Law in porous fluid flow. Notably, Darcy's Law can be obtained from first principles using homogenization arguments [52]. Our use of equation (2.2), however, will be different as cells preserve their geometrical integrity and do not flow freely like a fluid. Rather, the current in equation (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 with s(u) a source function that will be prescribed below. Specifically, equation (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. Equation (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. 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 equation (2.3) into the non-singular ε → 0 limit, the Laplacian equilibrium The most immediate boundary conditions from equation (2.3) are The domain Ω understood here generally consists of the bounded subset of R 2 or R 3 which is populated by the cells. Its boundary ∂Ω can be written as ∂Ω = ∂Ω D ∪ ∂Ω 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 to model, for example, an increasingly impenetrable cellular matrix as α → 0 by a homogeneous Robin condition. To simplify the presentation here, we employ homogeneous Dirichlet conditions (∂Ω = ∂Ω 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 Ω h the subset of voxels v i for which u i = 0. Similarly, let ∂Ω h denote the discrete boundary; this is the set of unpopulated voxels that are connected (i.e. share an edge) with a voxel in Ω h . See figure 2 for an illustration. We first consider the discrete version of equation (2.4): and where the source term is given by, in non-dimensionalized form, s(u i ) = 0 for u i ≤ 1 and s(u i ) = 1 whenever u i = 2, in view of the normalization p = 0 for the static case, following assumption 2.1. In equation (2.8), L is a discrete Laplacian over the currently active grid Ω 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, where v i here denotes the centre point of the ith 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 equation (2.2), a pressure gradient gives rise to a proportional current; I ∝ −∇p. Let us write I(i → 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: 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 equation (2.10), the current is reversed (j → 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 equation (2.2), including the subset of cells that are free to move as a result of this pressure gradient.   Let us write R(e) for 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. Under our present framework, u i is constrained to {0, 1, 2} and a total of three distinct cases emerges for the movement rates: and 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 D 2 , as 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 (2.14) In other words, the total pressure gradient over any closed surface ∂Ω is equal to the enclosed sources.
With ∂Ω 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 ≡ 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 ∂Ω. The rates in equations (2.11)-(2.13) are illustrated on a common scale with D 1 = D 2 = D 3 in figure 3. Note that, on this Cartesian grid, by the integration over an edge in equation (2.10), the current is non-zero only along the cardinal directions up/down and right/left, respectively.
In addition to the rates in equations (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 owing to cell-cell adhesion with cells in the nonempty voxel k as a ijk := αe ik min(0, r ij · r ik ),  Figure 3. Schematic of the model: as in figure 1, green voxels contain single cells and red voxels contain two cells, giving rise to pressure sources. This pressure is propagated by a Laplacian operator over the voxels and induces a force, and hence a rate to move, for the cells in the subset of voxels indicated by the arrows. Cells in boundary voxels may move into empty voxels and cells in doubly populated voxels may move into voxels containing fewer cells. All arrows are drawn on a common scale.
voxels of voxel i and scaling by the inter-voxel distance d ij we thus evaluate the net adhesion current contribution as which is to be added to the current in equation (2.10). The resistance force works against the pressuredriven current in equation (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 modelling situation.

Results
We now present some computational results for the proposed modelling framework. In doing so, we focus on highlighting the three distinct advantages of our approach: (i) the mechanics is based on constitutive equations obeying a well-understood physics; (ii) the framework is fast and scalable, and therefore convenient to experiment with; and (iii), 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 two dimensions and in one case a Cartesian grid in three dimensions. 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 [53]. 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 are available, and we refer the reader to 'Data accessibility' 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 modelling framework with intracellular continuous-time processes, including cellto-cell signalling and, finally, in §3.4 we demonstrate how the complex behaviours of a tumour model may be investigated.

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 [54]. We test for grid artefacts   by artificially forcing a number of cells into a square configuration and then relaxing the system to equilibrium (figure 4). In the absence of any other mechanical forces, clearly, the equilibrium state should be approximately circular.
In figure 5, 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 are still visible (left/right and up/down), there is no dependence on the underlying grid except for the natural distortions owing to discreteness. This comes from the fact that the physics is based on well-understood constitutive equations and that we employ a grid-consistent discretization of the Laplacian (equations (2.4) and (2.8)).

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 [55]. The cellular pressure, p, of the neurons is now described by the chemotaxis equation where χ 1 is the chemotactic sensitivity of the neurons to the chemical S, i.e. the concentration of Slit. The parameter ε p is taken to be small and represents our assumption that the pressure equilibrates on a faster time-scale than the movement of the neurons. The source term, s(u), is as defined previously. Equation (3.1) represents the movement of cells down the gradients of the Slit concentration. We set the domain for the experiment to be x = (x 1 , x 2 ) ∈ R 2 and the equation governing the concentration of Slit is given by 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 ε S is again taken to be small and represents the assumption that the diffusion of Slit occurs on a faster time-scale than the movement of the cells. Imposing assumption 2.2 thus allows us to solve the steady-state problem where the boundary conditions for the Slit chemical are We note that equations (3.4)-(3.6) can be solved analytically to give We also include a pressure-independent movement mechanism, representing the active motion of the cells owing 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: where χ 2 is a measure of the affinity of the cells to Slit for this type of pressure-independent movement. We initialize 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 [χ 1 , χ 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 6. As expected, on average the cell population moves away from the source of Slit, consistent with it being a chemo-repellent.

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 [56], which produces where denotes differentiation with respect to time and (3.10) We take parameters [a, b, v, k, h] = [0.01, 100, 1, 2, 2] and the average in equation (3.10) is taken over the set of neighbours N i of the cell i. Additional stochastic noise terms may also be added to equation (3.9), 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 cells = 1000. The average in equation (3.10) 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 7, a typical time-series of this model is summarized. Here the time-scale of the Delta-Notch dynamics, equation (3.9), is on a par with that of the proliferation process, so that their dynamics equilibrates after the tissue stops growing. A second simulation is summarized in figure 8, where the right-hand side of equation (3.9) 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 time-scales for patterning compare with those of tissue growth.

Non-trivial dynamics emerges from the combination of grid-based physics and local behaviour
Finally, we use 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  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 ∂Ω h , the collection of empty voxels that are adjacent to occupied voxels. We also have an equation for the oxygen concentration, c, which and c i = 1, i ∈ ∂Ω ext , (3.12) 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 ith voxel (doubly occupied voxels consume twice as much oxygen and empty voxels or those containing dead cells consume no oxygen). From these discretized 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 ρ 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) for other cells to move in to it. We also include the rates for cell movement as shown in equations (2.11)-(2.13).
In figure 9, we present snapshots from a realization of the model with parameters 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 realization (figure 9, bottom right). The widely held view in the field is that a limited oxygen supply leads to a stable, finite-sized tumour [57][58][59]. Recent models for avascular tumour growth have been used to successfully predict the growth of tumour cell lines that exhibit sigmoidal growth curves [60]. Sigmoidal curves have three phases: for an initially small population of cells (figure 9a) in an abundance of oxygen we see exponential growth (figure 9b), which is then followed by a quasilinear growth of proliferating cells ( figure 9c). Eventually, the lack of oxygen causes three heterogeneous cell types to appear (figure 9d); 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 9, bottom right (d)). However, owing to a possible lack of surface adhesion in our model we observe a fourth phase emerge, where asymmetric protrusions appear on the surface (figure 9e) of the tumour which increases the total uptake of oxygen and hence a further incline in the total cell number (figure 9, bottom right). We sought a parametrization 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 (λ, κ prol , κ death ) can be chosen to decide where the proliferating ring and necrotic core should lie. Next, initialize a population of cells at the required radius with some proliferation and death rates ρ prol and ρ death , and temporarily let ρ deg = ∞ 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 (ρ prol , ρ death , ρ deg ) in an attempt to build a finite-size tumour. Despite this comprehensive exploration of the parameter space, which is only possible owing 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 a 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 [48,61] which both conclude that additional mechanisms are necessary in order for tumour growth to be stabilized. The suggested mechanisms involve the necrotic cells secreting chemicals that can either act as a growth inhibitor for healthy cells to become quiescent [61], or as a chemotactic signalling chemical to attract tumour cells to migrate back towards the necrotic core [48]. The inclusion of such mechanisms fits naturally into our proposed DLCM framework and will be implemented in future, more focused studies.

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 time-scale 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 [49]. 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 [45,47]. We have demonstrated the use 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 factorization of the discrete Laplace operator and, as such, scales conveniently to about 20 000 cells. Using a compiled language and parallelization, one could probably improve upon this figure to some extent. However, real savings in computing time will require optimal Laplace solvers that use, for example, algebraic or geometric multigrid techniques [62,63] 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.
Data accessibility. The computational results can be reproduced within release 1.4 of the URDME open-source simulation framework, available for download at www.urdme.org.