Derivation of continuum models from discrete models of mechanical forces in cell populations

In certain discrete models of populations of biological cells, the mechanical forces between the cells are center based or vertex based on the microscopic level where each cell is individually represented. The cells are circular or spherical in a center based model and polygonal or polyhedral in a vertex based model. On a higher, macroscopic level, the time evolution of the density of the cells is described by partial differential equations (PDEs). We derive relations between the modelling on the micro and macro levels in one, two, and three dimensions by regarding the micro model as a discretization of a PDE for conservation of mass on the macro level. The forces in the micro model correspond on the macro level to a gradient of the pressure scaled by quantities depending on the cell geometry. The two levels of modelling are compared in numerical experiments in one and two dimensions.


Introduction
Mathematical modelling for simulation of cell populations is a tool complementing experimental studies to understand the complex biochemical and mechanical interactions between the cells in aggregations of unicellular organisms in bacterial cell colonies (Hellweger et al. 2016) and in multicellular organisms forming growing and developing tissues (Liedekerke et al. 2015;Lowengrub et al. 2010).
The mathematical models are either continuum models or discrete models. In continuum models, the evolution of the cell population can be modelled by the solution of time dependent, nonlinear partial differential equations (PDEs) for the cellular densities and concentrations of chemical compounds (Brodland et al. 2006;Frieboes et al. 2010;Humphrey 2003). In contrast, each cell, its internal state, and its motion are followed in time dependent discrete models (Liedekerke et al. 2015). There are techniques for analysis of a macroscale, continuum PDE model, for inference of parameters from experimental data on the population level, and established methods for the numerical solution of it. The continuum model is computationally efficient but it does not resolve any cellular details and heterogeneities in the population. A microscale, discrete model incorporates details of the cell behavior such as local chemical and mechanical interactions between cells, cell proliferation, and cell death. The number of cells can be a few to billions (Byrne and Drasdo 2009;Kang et al. 2014). Large numbers of cells may be necessary to bridge the gap between the individual cell level at the μm scale and the population level at the mm or cm scale. The problem is that the computational effort to simulate billions of individual cells may be prohibitive. By coarse-graining (or upscaling) the microscopic forces to the macroscale, large savings in computing time are possible.
In this paper, we derive novel continuum models from discrete models in one, two, and three dimensions (1D, 2D, 3D) for the biomechanical forces between the cells using the similarities between certain discrete models and discretizations of PDEs by a finite volume method (FVM) and of a pressure gradient. The force terms in the systems of ordinary differential equations (ODEs) governing the discrete models are identified as spatial discretizations of PDEs by the FVM. For a complete cell model we also need models for the motion in response to chemical gradients as in chemotaxis and haptotaxis, equations for the evolution of chemical species internally in the cells e.g. for the gene expression and cell metabolism, and externally e.g. for the signalling between cells, and transport of nutrients and oxygen. These models can also be continuous or discrete but are not addressed here.
Two different kinds of discrete models for the biomechanics of cells can be distinguished: on-lattice and off-lattice models. The cells are constrained in space to a lattice in on-lattice models. Each lattice point in a cellular automata (CA) model can host one or more cells. The cells move by stochastic jumps between the lattice points. A cell occupies several lattice sites in a cellular Potts (CP) model allowing greater geometric flexibility (Liedekerke et al. 2015). The cells move off-lattice and continuously in space and time in an agent based model (ABM) (Liedekerke et al. 2015;Osborne et al. 2017). Each individual cell can be modelled as an agent in an ABM and the motion of the cell satisfies Newton's equations of motion. The mechanical properties of the cells can be described in more detail in an off-lattice model with separate moving cell entities compared to the on-lattice CA and CP models but the latter are simpler and faster to simulate. Other advantages and disadvantages of ABM compared to CA and CP are discussed in Hellweger et al. (2016) and comparisons are made between CA, CP, and ABM models in Osborne et al. (2017).
A center based model (CBM) and a vertex based model (VBM) are ABMs where each cell is represented individually. The geometry of a cell in a CBM is a circle in 2D or a sphere in 3D in the overlapping spheres (OS) model (Drasdo and Höhme 2005) or a polygonal Voronoi cell defined by Voronoi tesselation (VT) (Kennedy et al. 2016). The forces due to the interactions between neighboring cells and external forces are applied at the cell center in OS and VT models. The adhesion and repulsion forces depend on the distance between the cell centers and parameters in various ways in the models (Ghaffarizadeh et al. 2018;Liedekerke et al. 2019). A summary and an evaluation of different CBM forces are found in Mathias et al. (2020). The cell membrane is modelled by a polygon (2D) or a polyhedron (3D) obtained by Voronoi tesselation in a VBM (Farhadifar et al. 2007;Fletcher et al. 2013Fletcher et al. , 2014Honda et al. 2004;Murisic et al. 2015;Nagai and Honda 2001;Staple et al. 2010;Weliky and Oster 1990). A motivation for the method is that the Voronoi polygon or polyhedron resembles the cell shape in certain tissues (Schaller and Meyer-Hermann 2005). The forces are applied at each vertex of the membrane and depend on the cell volumes and the perimeters of the cells adjacent to the vertex. Some vertex models are derived from an energy potential (Farhadifar et al. 2007;Murisic et al. 2015;Staple et al. 2010) which is minimized for the stationary solution. A local minimum of the potential will be reached depending on the initial state. All models have constant parameters and the sensitivity to these parameters in a VBM is investigated in Kursawe et al. (2017). The CBM and VBM are implemented in software such as Biocellion (Kang et al. 2014) , CellSys (Hoehme andDrasdo 2010), Chaste (Cooper et al. 2020;Mirams et al. 2013), MecaGen (Delile et al. 2017), and PhysiCell (Ghaffarizadeh et al. 2018). CP models are implemented in CompuCell3D (Swat et al. 2012) and Morpheus (Starruß et al. 2014).
The motion of the cells in CBM and VBM is overdamped in a viscous environment with small intertial terms compared to the disspative terms (Danuser et al. 2013;Fletcher et al. 2013) and the second derivatives in time in Newton's equations can be neglected. Thus, the CBM and the VBM are governed by a system of ODEs of first order in time for the coordinates of the cell centers (CBM) or the cell vertices (VBM).
The above deterministic models can be extended in different directions. A stochastic term models the micro-motility of a free cell as Brownian motion in the overview in Earnest et al. (2018) and also in Buttenschön et al. (2018), Middleton et al. (2014). The cells are ellipsoidal in Jin and Marshall (2020) with forces due to adhesion, lubrication, and hair-like appendages and bacterial biofilms are simulated with the model. Problems with growing domains are treated in Baker et al. (2010). Non-local effects between the cells correspond to integro-differential equations on the macro level (Buttenschön et al. 2018;Giniūnaite et al. 2020;Middleton et al. 2014). The relation between the micro level ABMs and macro level ODEs and PDEs has been studied in An et al. (2017), Byrne and Drasdo (2009), Drasdo (2005), Fozard et al. (2010, Liedekerke et al. (2015), Osborne et al. (2010) aiming at finding a direct correspondence between the models and their parameters at the two levels of modelling. If an expensive micro model is replaced by a cheaper macro model, then computational time and space can be saved. One approach is to first obtain an equation for the probability density of finding a cell at a certain position and then from that derive a PDE for the mean field of the cell density. Another approach is to regard the cell model as a finite difference approximation of a PDE which will be the macro model. The discrete micro level and the continuum macro level have been coupled in a number of papers, mostly in 1D. A nonlinear diffusion PDE in 1D is derived from a spring model for the cell forces in Murray et al. (2009Murray et al. ( , 2012. The diffusion coefficient is proportional to the spring constant and the viscosity coefficient and inversely quadratically proportional to the cell density in Murray et al. (2009). Another 1D continuum model is derived in Fozard et al. (2010) from an ABM with linear spring forces between the cells and drag due to the ambient substrate. The linear spring parameters in the discrete model in 1D in Murphy et al. (2019) are space dependent. The analysis shows how the parameters on the macro and micro levels are directly related. A mean field PDE in 1D for the motion and proliferation of cells is derived from a stochastic individual-based model in Chaplain et al. (2020). The resulting PDE is a diffusion equation with a source term. The CP model with stochasticity is the micro model in 1D and 2D in Lushnikov et al. (2008). A PDE with a nonlinear diffusion coefficient is derived there for the cellular density. A vertex model for an epithelium is homogenized to arrive at a linear elastic thin plate model in 2D at the PDE level which is valid on long length scales in Murisic et al. (2015). By adding a term in the potential of the vertex model, deformations in 3D are possible. A more elaborate cell geometry than that defined by Voronoi tesselation for a VBM is proposed in Jensen et al. (2020) to obtain a planar stress model at the macro scale.
The tissue level and the cell level are merged in computational hybrid models making them suitable for multiscale simulations of large aggregations of cells. Some examples follow. A PDE model for elasticity on the macro level is combined with a micro model on the cell level in Ghysels et al. (2009). The model for individual cells has spring forces and cell volume preservation. The PDE is discretized by a finite element method and the micro model is used at discrete points in space as in the heterogeneous multiscale method (Weinan et al. 2007). The forces in the stochastic on-lattice method in Engblom et al. (2018) are given by the solution of a PDE. A computational hybrid model of a tumour with individual cells in a proliferating outer layer and a continuum model in the interior of the tumour is proposed in Kim et al. (2007). In this way, expensive cell simulations are avoided in the quiescent center of the tumour.

Outline of paper
Here we establish relations between the mechanical properties of biological cells modelled by the CBM and the VBM on the discrete micro level and PDEs on the continuum macro level in 1D, 2D (CBM, VBM), and 3D (CBM). A center of cell i in a CBM with coordinates x i on the micro level is advanced in time by the ODE where v i is the velocity, f i j is the force between cell i and cell j, and the sum is over all other cells in the neighborhood of cell i. The same type of equation is solved in the VBM for the vertex coordinates. Then the sum in (1) is taken over all other vertices in the neighborhood. The equations for the coordinates of the cell centers and vertices are written in Lagrange coordinates following the flow. The spatial domain with the cells is tiled by Voronoi tesselation with one cell center in each Voronoi element defining the geometry of the cell as in the CBM-VT model.
The PDE for the evolution of the cell density ρ on the macro level is with a relation between the velocity and the pressure p according to Darcy's law for flow in a porous medium. It is written in Euler coordinates fixed in space. For closure of (2), a constitutive relation between p and ρ is needed. The space derivatives in (2) are discretized by a FVM on the Voronoi mesh. There is one cell per interval (1D), area (2D), or volume (3D) element in the discretized PDE and consequently, the cell density ρ is the reciprocal of the element size. By comparing the expressions for the macroscale velocity of the discretization of the PDE in (2) with the microscale velocity in (1) for the CBM, we find that the appropriately scaled pressure agrees with the CBM forces. Similarly for the VBM, the pressure gradient in (2) can be identified with microscale force terms in (1) after scaling. The scaling for both CBM and VBM depends on the cell geometry. Since it is not known in detail at the PDE level in 2D and 3D, some geometrical assumptions are necessary to determine the scaling of the forces between individual cells to obtain the pressure. The main contribution of the paper is the derivation of the relation between the CBM and VBM forces f i j on the micro level in (1) and the pressure p in the PDE (2) on the macro level.
When the PDE has been obtained, it can be solved numerically on any mesh, not only on the one defined by the geometry of the cells. If the variation of ρ is small, then a much coarser mesh than the Voronoi mesh for the cells is possible with savings in computational work and memory. In numerical experiments, ABM simulations of distributions of cells of high and low density are compared with PDE solutions on equidistant grids in 1D and Cartesian grids in 2D.
The paper is organized as follows. In the next section, the geometry of the cells and the forces between them are specified. In Sect. 3, macro level pressure is derived from the micro level descriptions by comparison with discretizations by FVM and a pressure gradient on Voronoi meshes. The PDEs are discretized by FVM on an arbitrary mesh in Sect. 4 using the pressure formulas. Section 5 contains the numerical experiments with the discrete and continuum models and conclusions are drawn in the final section.

Cell geometry and forces
The geometry of a biological cell is determined by the Voronoi tesselation based on the cell centers. The cell is an interval in 1D, a polygon in 2D, and a polyhedron in 3D. The cell centers are connected by the edges of the triangles (2D) or tetrahedra (3D) in the associated Delaunay triangulation of the cell centers. The forces in the CBM act along the Delaunay edges and in the VBM the forces act on the vertices of the polygons along the Voronoi edges. This section is a summary of the notation, the cell geometry, and the CBM and VBM forces. Vectors are written in boldface A and x andẋ denotes a time derivative dx/dt of x.

Cell geometry
The total number of cells in the system is N .
. . , N , of the center of the cell i. The position of cell center j relative to cell center i is r i j , i, j = 1, . . . , N , the distance between the centers is r i j , and the direction between the neighboring centers isr i j The full coordinate vector of all cells is denoted by . The coordinates of the nodes or vertices α, β, . . . on the perimeter of the cell are x α , x β , . . . and their relative positions are as in (3) After Voronoi tesselation there are a primal mesh and a dual mesh in Fig. 1. The primal mesh consists of the Voronoi elements representing the biological cells. The Delaunay triangulation defines the dual mesh composed of triangles in 2D. There are four primal elements with centers i, j, k, l, and two dual elements with centers α, β in the figure. The vertices β, γ , δ are connected to node α in the primal mesh. The three directionsr αβ ,r αγ ,r αδ start in vertex α and go to vertex β, γ , and δ along edges in the primal mesh. The three edgesr i j ,r ik ,r il in the dual mesh connect the element centers in the primal mesh, see Fig. 1. In an expanded mesh, γ, δ, ε, and ζ will also be triangle centers in the dual mesh.
The interval length in 1D, the element area in 2D, and the element volume in 3D of element i are denoted by V i . The perimeter length of an element in 2D is a i and a i j is the length of the edge separating elements i and j in 2D. Similarly, the area of the with element centers at i,j,k, and boundaries at α, β, γ , δ surface between elements i and j in 3D is denoted by a i j and the perimeter area is a i . In Fig. 1, a i j = r αβ . Since each Voronoi element corresponds to one biological cell, the cell density ρ is ρ i = 1/V i in element i and V i is the specific volume. Let J i be the index set of all adjacent elements to i. Then the perimeter a i of element i in 2D and 3D is The elements are intervals in 1D. The element indices are i, j, k, in Fig. 2 and the node indices are α, β, γ , and δ. The midpoint of cell i is x i and its length is

Center based models
The cell defined by a Voronoi element in a CBM can be associated with a circle in 2D or a sphere in 3D. Suppose that two circular (2D) or spherical (3D) cells are partly overlapping each other in the OS model. Then there is a force between them repelling the cell centers. There is an attraction force between two cells that are in the vicinity of each other without touching. The force depends on the distance r between the cell centers.
The strength of the force is denoted by g(r ). Let g(r ) be defined for r ≥ 0 with the properties where the parameters r A and s are positive. The cell centers are repelling each other when r < s and attracting each other when s < r < r A . The force vanishes at r = s and for r > r A . The following example of g satisfying (6) is from Mirams et al. (2013): There is a discontinuity in g(r ) at r = r A . The parameter μ determines the strength of the force and c > 0 the decay rate of the adhesion for r > s. Another example is found in Ghaffarizadeh et al. (2018); Macklin et al. (2012): This g(r ) is continuous for r ≥ 0. The adhesion or attraction parameter is c a and the repulsion parameter is c r . The force model in (8) also satisfies (6) if c a < c r and s is In the Hertz model, the repulsion force depends on the cell overlap r − s with g(r ) ∼ −(s − r ) 3/2 while the adhesion force has a more complicated dependence on r − s (Liedekerke et al. 2015). The Johnson-Kendall-Roberts (JKR) model also depends on r − s and requires the solution of nonlinear equations for g (Carpick et al. 1999;Drasdo and Höhme 2005). Other examples of CBM forces are found in Mathias et al. (2020). The force on cell i caused by cell j has the magnitude g(r i j ) and the directionr i ĵ Sincer i j = −r ji and g(r i j ) = g(r ji ) the forcing on the cell centers i and j is of equal strength but in opposite directions. Then the system of d ODEs for the center coordinates of cell i is as in (1) given by Newton's equations of motion with the viscosity η, neglecting the accelerations because they are small (Danuser et al. 2013;Fletcher et al. 2013),ẋ and summarized for all cellsẋ The Lagrange coordinates x follow each cell center or how the mass moves in space. This model has the same definition of the force in 1D, 2D, and 3D. There are a limited number of g(r i j ) = 0 in (10), see Sect. 3.4. Depending on r A and the cell size only indices in J i may contribute to the sum.

Vertex based models
The forces are applied in the corners of the cell polygon in 2D in a VBM according to Weliky and Oster (1990). Let the areas of the cells be V = (V 1 , V 2 , . . . , V N ) T and their perimeters be a = (a 1 , a 2 , . . . , a N ) T . Then the force acting on node α due to node β isr with a directionr αβ and a magnitude f (V, a) depending on the surrounding V i and a i . The velocity of node α satifieṡ where J α = {β, γ , δ} in Fig. 1. This equation corresponds to (10) for the CBM. The Lagrange coordinates x α for the nodes define the shape of the 2D cell. The generalization to 3D and its implementation are more complicated (Honda et al. 2004). The Weliky-Oster phenomenological force model Weliky and Oster 1990) at x α depends on the area and perimeter of cell k as follows (see Fig. 1) The area force coefficient is ς and κ is the perimeter force coefficient. Sum over all three cells with a corner at x α for the velocity of node α With the force function in (13), the velocity is Another example is the Nagai-Honda model , App. A), (Nagai and Honda 2001), which is derived by taking the gradient of an energy function. The force at x α due to V k and a k is almost a linearization of (14) Compared to Nagai and Honda (2001), the force in (18) is slightly simplified here in the factor r αβ in the first term. The parameter λ is a deformation energy coefficient, β is a membrane surface energy coefficient, and γ is a cell-cell adhesion energy coefficient. The ideal area is V 0k and the ideal perimeter length is a 0k . They are reached in an equilibrium configuration when the forces vanish. As in Nagai and Honda (2001), the forces in the VBMs in Farhadifar et al. (2007), Murisic et al. (2015), Staple et al. (2010) resemble the force in (18) and are derived from a potential.
Consider a small perturbation from equilibrium in cell k (or i or j) in (14) and (18). The perturbations in v α are equal with both models if

Comparison of the microscale models with a macroscale PDE
A macroscale PDE is derived for the cell density ρ in this section. The density is transported by a pressure gradient. The PDE is discretized by a FVM on the Voronoi mesh defining the biological cells. The discretized pressure gradient on the macro level corresponds to the forces in the CBM and the VBM on the micro level. The geometrical detail of the cells used in the CBM and the VBM on the microscale is missing on the macroscale. The cell size is available as 1/ρ but, for example, the cell perimeter is not known. By assuming a regular shape of the cell, the size of the edges and the distance between the center and the vertices can be obtained as functions of the cell size. The pressure will then be a function of ρ determined by the specific CBM or VBM. The viscosity coefficient in (10) and (13) is assumed to be η = 1 and is ignored in the formulas. It can always be removed by including it in the force parameters.

The conservation law
A fluid parcel has volume V and surface S. The initial position at t = 0 is at x 0 and its velocity field is v. Lagrange coordinates x(t) follow the parcel. This parcel is small on the macroscale and will later be identified as the specific volume V of a cell. Use Gauss's formula to arrive at a geometric conservation law for the time evolution of V at In Euler coordinates with fixed x in space we have with the total derivative D/Dt Introduce the density ρ = 1/V in (21) to obtain a PDE for conservation of mass in a domain Ω ∂ρ ∂t On the boundary of Ω with normal n, the boundary condition is n · ∇ρ = 0. Then there is no flux of mass across the boundary and the total mass Ω ρ dV is constant in time. A 1D model in Fozard et al. (2010) adds an extra term in (22) for internal viscosity in the cells resulting in a linear spring model for the force.
The PDE (20) is discretized in space by a FVM for cell i with volume V i and center at x i in the Voronoi mesh defining the cell geometries, see Fig. 1. The time derivative is approximated by the Euler forward method at discrete time points t n , n = 1, 2, . . . , The sum is taken over J i , the set of all indices of cells sharing a common boundary with cell i (at least j, k, l in Fig. 1). The outward normal on the edge (2D) or surface (3D) between cells i and j is n i j and v i j is also evaluated on that edge or surface. For the conservation law (22), we consider any fixed Euler mesh with element (or computational cell) area or volume ω i . The mesh here is fixed in space contrary to the Voronoi mesh for the biological cells which is continuously deformed in time. The elements i and j ∈ J i have a common perimeter of size σ i j with an outward normal n i j from cell i to j. The time derivative is approximated by the Euler forward method at discrete time points t n , n = 1, 2, . . . , starting at t 0 = 0 with timestep Δt = t n+1 − t n . With FVM discretization in space, the equation for ρ at t n+1 is where v i j is evaluated on the edge or surface σ i j . The sum is taken over J i , the set of all indices of elements sharing a common boundary with element i. The density at σ i j is approxmated by 0.5(ρ i + ρ j ). Now there are two alternatives: compute the cell density with the continuous model (24) or first advance the cell centers with the discrete model (10) or (13) and then compute V i using (23) or the present geometry of the cell and ρ i for each cell i. If the elements can be chosen much larger than the cells, ω i V i , then the first alternative will save computing time and memory. This is possible if ρ varies slowly in space.
We need a relation between v and ρ for closure of (22). In Darcy's law, v is proportional to the gradient of a pressure p v = −μ∇ p.
Inserted into (22) the PDE is The equation is closed with a constitutive relation between p and ρ: p = p(ρ). This is the macroscopic PDE for the evolution of ρ in Chaplain et al. (2020) where a source term for injection of mass on the right hand side models cell proliferation on the macro level (see also Byrne and Drasdo (2009)). One cell splits into two in cell division at the micro level thus doubling the cell density gradually or instantaneously and very locally in space. Then the numerical solution of (22) has to resolve a density peak in space and time requiring refinements in the temporal and spatial discretization. If the proliferation occurs frequently in the domain, then an averaged, smooth source term is preferred to model cell division. The pressure p expressed in V or ρ will be determined by the CBM and the VBM in the next two sections.

Center based models
Apply FVM to (25) in a cell i. Then the discretized law is where p i j is computed at the midpoint of the interface between cell i and j. With the cell indices in the neighborhood of cell i in J i , the equality (10) for the CBM is rewritten in A g with the properties (6) is nonzero for a limited number of j ∈ J i . Since n i j =r i j , the same velocity is obtained in (27) and (28) if p i j is chosen as This is the bridge between the macro and the micro levels but geometric relations between V i , a i j , and r i j are missing to obtain a macroscale pressure p(V ).

Vertex based models
The relation (25) is discretized in 2D on the dual triangular mesh, see Fig. 1. The triangle with center at x α has three edges between triangle α and three neighbors β, γ , and δ. The normals of the edges are n αβ , n αγ , and n αδ and are equal tor αβ ,r αγ , and r αδ . By (25), the pressure and the velocity at the midpoint of the edge between triangle α and β are assumed to satisfy where F α and F β are the VBM forces at x α and x β in (15). After discretization of (30), the relation between the forces and the pressures p α and p β at x α and x β using (17) is The pressures p α and p β are functions of the microscopic forces in (31) and the inner products involvingr, V i , and a i . They depend on properties of the geometry of the cell polygons.

Geometry of the cells
We need relations for the cells in (29) and (31) between their area or volume V , their perimeter a, a distance r , and the inner products such asr αβ ·r αγ to be able to find the macroscopic pressure by coarse-graining the microscopic forces. This is easy in 1D but in 2D and 3D complete geometric information is not available and assumptions on the cell geometry are necessary.

One dimension
Let r i be the distance between the center of cell i and the cell boundary in 1D in Fig. 2. Then the length of the cell is There is no perimeter but let a i j in (5) be 1 and J = 2. The geometric relations between V i , a i j , and r i j in (29) are The inner products are either −1 or 1.

Two and three dimensions
The PDE approximation of the forces is expected to be the most accurate when there is a smooth and small variation of the cell density. By assuming that the cell shapes are close to regular polygons in 2D and regular polyhedra in 3D, we can derive approximate expressions for the dependence on V of r i j and a i j in (29) and a i , r αβ , and the inner products in (31). The shortest distance between the cell center and the perimeter in cell i is denoted by r i which is equal to the radius of the inscribed circle. Hence, r i j = r i + r j . Table 1 The factor ξ 2 in the relation r = ξ 2 (J )V 1/2 in regular polygons J 2 D n a m e ξ 2 3 Triangle (3 tan(π/3)) −1/2 = 0.4387 4 Square (4 tan(π/4)) −1/2 = 0.5000 5 Pentagon (5 tan(π/5)) −1/2 = 0.5247 6 Hexagon (6 tan(π/6)) −1/2 = 0.5373 A regular polygon in 2D with J corners and J edges fulfills Introduce ξ 2 (J ) = 1/ √ J tan(π/J ). Then by (34) we find that The polygon becomes a circle when J → ∞ and ξ 2 → 1/ √ π . A triangle yields the lower bound on ξ 2 . The ξ 2 -values for different J are found in Table 1. The maximum number of inscribed circles of equal radius in other cells touching circle i is 6 (the kissing number (Conway and Sloane 1993)). In a densly packed cell colony, J will be close to 6. Choosing 5 or 6 for J in ξ 2 does not alter its value very much in the table and ξ 2 appears to vary slowly when J → ∞. The inner products betweenr αβ and otherr-vectors in (31) are equal to ± cos(2π/J ) if the cell is a regular polygon with J edges.
The growth of epithelia are simulated in Staple et al. (2010) with a planar vertex model. The statistics of the cell geometry shows that most cells have five or six neighbors. This is also the conclusion in Farhadifar et al. (2007) Lewis (1943), it is found that the average number of edges of a cell is six with a hexagonal shape.
A regular polyhedron in 3D with J faces or neighbors satisfies The quotient between r i and V 1/3 is written ξ 3 (J ). Then it follows from (36) that We obtain a sphere when J → ∞. The ξ 3 -values can be calculated for the five Platonic solids (Weisstein 2003). Four of them are tabulated in Table 2. As in Table 1, the value of ξ 3 increases slowly when J increases from 6 to ∞. The highest density packing of spherical cells is achieved by J = 12 and the kissing number for a sphere is 12 (Conway and Sloane 1993). Therefore, the icosahedron with J = 20 is excluded. The formulas in dimension d, d = 1, 2, 3, in (33), (35), and (37) are summarized by

Pressure in center and vertex based models
The geometric relations for the cells from Sect. 3.4 are introduced in the pressure formulas from Sects. 3.2 and 3.3.

Center based models
Insert r i j = r i + r j and V i /a i j from (38) into the pressure formula for CBM (29) to arrive at for the pressure between cell i and j. It can be interpreted as a discretization of at a point 1 2 (x i + x j ). In d dimensions, the pressure formulas for p d are with the number of neighbors J = 2, 5, 10, for d = 1, 2, 3. Let g be as in (7) and let s = 2ξ d V 1/d 0 . If V < V 0 then p is positive and decreases when V increases. This is the usual behavior in fluids: the pressure decreases when the density (∼ 1/V ) decreases.

Vertex based models
With the VBM force in (15) in 2D, f in (31) can be written as a sum By the geometric relations in 2D we have Use (38) with d = 2 and r αβ = a j /J and identify p α and p β in (43) with to satisfy the equality. These are interpreted as discrete versions of the pressure in Thus, the pressure corresponding to the Weliky-Oster model in (15) is In summary, we have derived macroscopic constitutive laws for the pressure from the microscopic CBM and VBM in (40) and (46) by comparing a discretization of a pressure gradient with the microscopic forces. Only the cell area or volume V is available on the macro level but the CBM and VBM forces depend on more details of the cell geometry such as the perimeter a and the distance r . The missing information is obtained in (38) by assuming that the cells are regular polygons or polyhedra.

One dimension
A comparison between the CBM and the VBM for small perturbations in 1D shows how the parameters in the methods are related. A perturbation δV i from the equilibrium r = s in CBM and V = V 0 in VBM is introduced in the force formulas in Sects. 2.2 and 2.3. The notation is as in Fig. 2.
The CBM velocity of cell i according to (28) is Perturb r i , r j , and r k about r i j0 = s. Then the change in velocity is The VBM due to Weliky-Oster in (15) depends on the perimeter a. It has no immediate interpretation in 1D and is ignored here. The velocity at Perturb the cell lengths as in (48) δv to obtain Comparing δv i in (48) and (51) yields g (s) = ς/V 2 0 . With g(r ) in (7), g (s) = μ = ς/V 2 0 .

Finite volume method for the PDE
The equation for conservation of mass (22) is advanced in time by time stepping with the Euler forward method and a FVM discretization in space on an arbitrary mesh as in (24) in this section. A finite difference method or a finite element method would be a possible alternative for the spatial discretization. The pressure in Darcy's law (25) is defined by the CBM or the VBM on the micro level as in Sect. 3.5. The conservation law for transport of density is a nonlinear parabolic PDE for ρ in stationary Euler coordinates in (26). Insert p(ρ) into the right hand side of (26) The expression ρ ∂ p ∂ρ is a diffusion coefficient in (26). This is observed by Murray et al. (2012) and the coefficient is derived in 1D in a coordinate system different from ours. For ρ to remain stable in the diffusion equation, we require ∂ p ∂ρ > 0. The diffusion coefficient for the CBM in (40) In particular the coefficient in 1D is The diffusion coefficient is negative if ρ is such that g + 2ξ d ρ −1/d g < 0 in (53). This is possible e.g. when 2ξ d ρ −1/d is close to r A in (7) since g > 0 but small and g < 0.
A continuous g in (6) with g > 0 for s < r < r A and g = 0 for r > r A necessarily has this property. The diffusion coefficient for the VBM in (46) is The coefficient is positive when J ≥ 4 with stable solutions of (26). Compressible fluids are modelled by constitutive laws for the relation between p and ρ or V similar to (40) and (46) (Anderson 1990). A perfect gas satisfies where R is the specific gas constant and T is the temperature. The constitutive law for an isentropic gas such as air (γ = 1.4) is The assumption in Chaplain et al. (2020) is that p(ρ) = K γ ρ γ , K γ > 0, γ > 1, as in (57). The space derivatives in (26) are discretized by a finite volume method on an arbitrary stationary mesh with element areas (or volumes) ω i , ω j . The separating line (or area) between ω i and ω j is denoted by σ i j and the distance between the element centers is i j . Then the equation for ρ i in element i is The equation (58) is integrated in time from t n to t n+1 by the forward Euler method with the time step Δt resulting in The right hand side is evaluated at t n . On a Cartesian equidistant spatial grid with grid size Δx, we have σ i j = Δx d−1 , i j = Δx, and ω i = Δx d .

Multiply (59) by ω i and sum over all M elements in the mesh
The double sum over i and j vanishes since each term in the interior of the domain appears twice with opposite sign and a term on the boundary is zero due to the discretized Neumann boundary condition ρ i = ρ j implying p i = p(ρ i ) = p(ρ j ) = p j . Therefore, the total mass i ω i ρ i is conserved also in the discretization, cf. ρ in the PDE in (22). The pressure at i and j in (59) is approximated by the CBM in dD in (40) with the specific volume V k resulting in The constant J in (62) is chosen to be fixed with different values in 2D and 3D. It is 6 in the numerical experiments in 2D in the next section assuming that the cell environment is dense and slowly changing. The PDE model is not suitable in rapidly changing domains. There the distance r may not depend on V as smoothly as in (38) where a change in V i will change r i equally in all directions from the cell center. The approximation of p at x k with the Weliky-Oster model (46) in 2D is This p is inserted into (59) to advance the cell system in time with the VBM. The PDE is discretized on an equidistant grid in 1D with grid size Δx with CBM forces (40) and indices as in Fig. 2. Let ρ i j = 1 2 (ρ i + ρ j ) and update ρ in time by Apply FVM to (22) for a comparison of the PDE in 1D with the discrete CBM model. (22) is discretized by The difference between the PDE solution in (64) and the CBM solution in (65) is in the evaluation of the pressure and the variable cell size in (65). The pressure in the Weliky-Oster model in 1D is derived from (46). The constant 2/(J ξ 2 ) → 1/(2ξ 1 ) = 1 in 1D and f a = 0. Multiplying the force function by V results in a constant p This is avoided by letting a discrete pressure p j at x j to update element i be Then the integration of ρ n i is achieved by cf. (64) and (65). With a general macroscopic definition of the pressure as in (40), it is possible to integrate the vertex coordinates x α in (13) with the CBM force by taking The macroscopic pressure in the PDE is the same with this f . In the same manner, p in (45) can be used to advance the center coordinates in (10) with a g derived from (68).

Numerical results
The discrete methods and their PDE approximations in Sect. 3 are compared in this section in numerical experiments in 1D and 2D. The methods are implemented in Matlab in a straightforward manner using their voronoin for the Voronoi tesselation without any attempts to optimize the efficiency of the code. The parameters in the methods are displayed in Table 3. The distance r between two cell centers with no force between them is denoted by s in the CBM. The force is scaled by μ. The exponential decay of the adhesion force is given by c and the force vanishes when r > r A . In the VBM, ς and κ are the scalings of the area and the perimeter dependent force components, respectively  Table 3

One dimension
The density of the particles is computed in 1D using the CBM and VBM forces in the micro model and the corresponding macro level PDEs in a comparison of the methods in Figs. 4, 5, 6, 7. The parameters in the methods are found in Table 3. The relation between the parameters ς and κ in the VBM is important but not the scaling of them. It follows from (10) and (13) that a different parameter scale of μ, ς, and κ will only change the time scale of the evolution of the system.
A system consisting of N particles is simulated for t ≥ 0 in an interval A in space. The spatial cells or intervals are computed by a Voronoi tesselation of the particle system as in (32). Then the cell density in the interval is given by ρ i = 1/V i . The original positions of the cells at t = 0 are such that ρ varies smoothly. The initial density distribution is interpolated to a grid with constant step size Δx for use as initial conditions to the solution of the PDE discretized by FVM in the interval B. The cells at the boundaries of A are fixed and ρ satisfies Neumann conditions at the boundaries of B. These boundary conditions guarantee that mass is conserved in A and B. A stationary or equilibrium solution for the discrete model and the continuum model have a constant density with equal cell size.
The repulsion and adhesion between the cells is modelled by the CBM in (7) in the first example. The behavior of p(ρ) based on CBM with the parameters in Table 3 is depicted in Fig. 3. The particle simulation with N = 40 and the solution of the PDE in (64)  When ρ > 1 in Fig. 4 then there is a repulsive force between the cells in the discrete model and ∂ p/∂ρ > 0 in (54) in the PDE model (to the right in Fig. 3). The solution diffuses toward the constant state. An adhesive force acts between the cells when 1/r A = 0.667 < ρ < 1 and the force disappears when ρ < 0.667. The diffusion in the PDE is negative and unstable when ρ < 0.901 and vanishes when ρ < 0.667 (to the left in Fig. 3). For small perturbations about r i j = s, the effect is given by g (s) in (48) which is equal to μ with the force in (7). In Fig. 5, the solutions approach a steady state with two ρ-levels: one low with separated cells and one high with cells touching each other. Because of the instability in the PDE for low ρ, some of the ρ values in the solution decrease as time increases. The solutions are close to steady state at t = 3 and the high density solutions agree between the CBM and the PDE. The instability in the PDE causes oscillations in the numerical solution for the low density values. This phenomenon is not ameliorated by refining the mesh or reducing the timestep. Since the diffusion vanishes for low density values, the oscillations remain in the stationary solution.
The  is compared with the corresponding PDE in Figs. 6 and 7. The PDE is stable with a positive diffusion coefficient. The initial data and the boundary conditions are the same as in the CBM example above and the force coefficient is ς = 50. The vertices of 39 cells are advected in the micro model and the macro model has 20 grid points. As observed in Sect. 4, there is no directly corresponding continuous pressure in 1D for the VBM but a discretization of the PDE is still possible with the pressure (66) in (67) which yields a fair result. The PDE solutions diffuse toward a constant steady state as time progresses which is different from the steady state with CBM in Fig. 5. The discrete and the PDE solutions are close except for the peak in ρ. Peaks cannot be represented very well in a PDE discretization on a coarse mesh.

Two dimensions
The density of the particles is computed in 2D in this section using the CBM and VBM forces and the corresponding PDEs in a comparison of the methods. The parameters in Table 3 are used in the methods. A system consisting of 1777 particles in A = [0, 40] × [0, 40] is simulated in a time interval with t ≥ 0. The density of the particle system is computed by a Voronoi tesselation. The area V i of the Voronoi cell is first determined and then the density by ρ i = 1/V i . The original positions of the cells are located in a hexagonal pattern. It is disturbed in the center of the domain such that the density varies radially from there. An example is found in Fig. 8 with ρ ∈ [1.11, 1.97] in the central parts of the figure and ρ = 1.15 in the outer parts. The cells will be close to hexagons also in the steady state. Therefore, J is chosen to be 6 in ξ 2 (J  (59) with the grid sizes Δx = Δy = 1.263, see Fig. 8. The initial ρ in the PDE solution at t = 0 is interpolated to the grid from the initial CBM or VBM solution. The solutions of the discrete system and the PDEs are compared along a line with a constant y value through the center of the domain along the x axis.
The system with CBM forces is simulated for t ∈ [0, 0.45] with the timesteps Δt = 2.5 · 10 −4 . The CBM solution and the PDE solution of (62) are compared at five time points in Fig. 9. The PDE solution decays at a slower rate than the CBM solution. The explanation is that r in the force formula varies more rapidly in a neighborhood of large variations in ρ than V does in the PDE approximation r = ξ 2 V 1/2 thus inducing weaker forces and slower dissipation in the PDE solution. The cells at the center of the domain are less regular than assumed to obtain ξ 2 . Refining the PDE mesh does not change the result very much. The steady state solution is a constant density and the boundary conditions are such that the total mass is conserved. Close to the steady state The solutions with the VBM forces due to Weliky and Oster are compared at five time points in the interval [0, 0.6] in Figs. 10, 11, 12, and 13 for different parameters ς and κ. The timestep is Δt = 10 −3 . The difference between the particle simulation and the PDE solution of (59) with the pressure (63) is largest at the peak of the density. There the approximations of a length scale r by ξ 2 V 1/2 and the perimeter a by 2V 1/2 /ξ 2 are the least accurate to determine a pressure for the PDE.
In the figures, the solutions computed by the CBM and the VBM and the PDE approximations are compared. The computing time to obtain the solutions is at least an order of magnitude faster with the PDEs. For example, take the simulation in Fig. 10. Initialization of data structures necessary for a VBM simulation required 84.9 s on a 1.

Conclusions
The CBM and the VBM are models on a micro level for the mechanical forces between neighboring biological cells in time dependent motion. Each cell is treated as an entity whose motion is determined by the forces. PDEs are suitable models on a macro level for the motion of large aggregations of cells. We have derived these PDEs from two established CBM and VBM in Sect. 3 and compared the two levels of modelling in numerical examples in Sect. 5. The advantage with PDEs is that tissues with billions of cells are easier to simulate than the detailed particle models with the same number of cells. A disadvantage with a PDE is that local events such as cell proliferation cannot be well represented. The PDE on the macro level, corresponding to the CBM and the VBM on the micro level, is nonlinear and parabolic. It is an equation for transport of cell density ρ = 1/V where V is the area or volume of a cell. The forces in CBM and VBM are transformed to a pressure p as a function of V , p(V ), in the PDE and there is a direct translation of the parameters in the microscopic and the macroscopic models. By the relation between the forces and the pressure, a pressure obtained from a CBM force can be used to define a VBM force and vice versa. The microscopic forces depend on the cell perimeter a in 2D and the cell radius r in 2D and 3D. They are involved in the specification of the macroscopic pressure but are not readily available on the macro level. This is usually the case in a multiscale model that information on the micro level is missing on the macro level. Instead, a and r are approximated in the PDE assuming that the cell is a regular polygon in 2D and a regular polyhedron in 3D. This approximation is not necessary in 1D. The PDE derived from the CBM is unstable in a density interval because the diffusion coefficient is negative there. The motion of the cells is assumed to obey Darcy's law, which governs the flow of a fluid through a porous medium. An alternative would be to assume that the cell aggregation behaves as a solid tissue with shear forces such as a tumour described by a nonlinear elastic model ) and then relate a given microscale model to those equations or choose a micro model suitable for those equations.
There are discrepancies between the solutions in 2D computed with the discrete particle models and the PDE approximations in the numerical experiments in Sect. 5. These are explained by the approximations of a and r in the PDE which are less accurate at peaks in the density distribution where the cell geometry deviates the most from a regular polygon. In 1D, the agreement between the fine micro level and the coarser macro level solutions is good. A possible solution in dimensions higher than one is to derive transport equations for a and r in the same spirit as for V (or ρ) and use these a, r , and V to define p. This would be more complicated than it is to arrive at an equation for ρ and would require detailed knowledge on the macro level of the initial cell geometry to provide initial data for a and r . Another solution is to let the forces between the cells on the micro level depend only on the cell volume V to simplify the transfer to the macro level. In this way, the macroscopic model would define the forces in the microscopic model to which additional, detailed features could be added to motivate and refine a cell based approach.
The advantage with a PDE model is that the spatial discretization of it for numerical solution is chosen to resolve the variation in the density without regard to the size of the biological cell. If the cell size is almost constant and the density varies smoothly, then major reductions in computational work are possible by a coarse mesh on the macro level. With a micro model and a macro model one can treat parts of the cell domain by the CBM or VBM where detailed modelling is required and more quiescent parts by the PDE model and glue them together at an interface in a multiscale model as in Kim et al. (2007). The PDE would define conditions at the boundary of an embedded domain with discrete, microscopic modelling. That appears to be possible already with the PDEs in Sect. 5 if the variation in density is low in the PDE domain.