Particle methods to solve modelling problems in wound healing and tumor growth

The paper deals with a compilation of several of our modelling studies on particle methods used for simulation of wound healing and tumor growth processes. The paper serves as an introduction of our modelling approaches to researchers with interest in biological cell-based models that use particle-based modelling approaches. The particles that we consider in the present models mimic either cells or points on cell boundaries that are allowed to migrate as a result of several chemical and mechanical factors. A distinct feature of our modelling frameworks with respect to conven-tional particle models, is that cells, mimicked by particles, are allowed to divide, differentiate and to die as a result of apoptosis or any causes for cell death. The paper is merely descriptive, rather than written in full mathematical rigor, and will show some of the potentials of the applied modelling.


Introduction
Wound healing and tumor growth are both biological mechanisms that are very common in mammalian tissues and organisms, where the first process is indispensable for the survival of the organism or human beings, whereas the latter one is often life-threatening to humans. A good understanding of the biological mechanisms is crucially important for an improved treatment related to both processes. In societies with ageing populations, elderly people are more sensitive to B F. J. Vermolen F.J.Vermolen@tudelft.nl 1 Delft Institute of Applied Mathematics, Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherland tumor development, which makes cancer treatment more and more important. Besides cancer, elderly people are also more likely to become bed-bound patients and/or to suffer from diabetes, which makes them more likely to develop pressure ulcers and/or to develop ischemic wounds that will heal very slowly in the more fortunate cases. Since the optimisation of treatment requires a thorough understanding of the biological processes involved, scientific research is crucially important. This scientific research is often based on in vivo and in vitro experiments where one either analyses clinical data or data from modelling samples which for instance could originate from mono-layer cell cultures. Scientific research aims at revealing links between various processes and parameters, along with the development of hypotheses. Since the validation of the hypotheses with the experimental data needs a quantification of the relations between respective processes and parameters, mathematical relations and descriptions of the processes and parameters come into play. These mathematical relations and processes are the backbone of the mathematical models that are used to simulate the various biological (sub-) processes.
The mathematical models for tumor growth and wound healing often consider various aspects of the overall process due to the enormous complexity of the processes involved. Furthermore, the models can have different natures and be applicable for a specific scale only, such as -Space-Free level, where one can model DNA of cells, often stochastic, signalling, systems biology and network models are used here. Besides this one can also use this modelling type by using ordinary differential equations on very large scales for processes like wound healing or tumor growth. On DNA level, we report [26,30] as references on tumor sensitivity with respect to DNA. A review on macroscopic modelling combined with ther-apies related to wound healing is given in [18]. Further examples are [26,31]. An example of a system of ordinary differential equations is given in [55] where angiogenesis in the context of solid cancer growth is modelled; -Compartmental level, where one models the dynamics in several components of one cell, this approach could be partial differential equations based where one models diffusion over or through cell membranes, such as in [32]; -Cellular-Based Models -Cell-shape evolutionary models, where one models cell deformation and migration or other processes of one single cell or at most a low number of cells, here one uses either partial differential equations, large coupled systems of ordinary differential equations, or cellular automata models (e.g. Cellular Potts Modelling), see for instance the work in [4,33,34,37,48]; -Constant cell-shape models, where one models a large number of cells for migration and sometimes including differentiation, division and death, using (stochastic) ordinary differential equations sometimes combined with partial differential equations for the distribution of chemicals interacting with the cells. Here also cellular automata models are used. Examples are in [7,20,41,49], for tumor-induced angiogenesis we refer to [2]; -Tissue-level, where one models cells using upscaled densities, simulated by the use of partial differential equations, where in some models moving boundary formulations or spinodal decomposition theory have been adopted. Examples are in [1,2,12] for tumors and [25,43,45] for wound contraction and healing.
In this paper, we will describe some of the results obtained using modelling on a cellular level where deformation of cells is taken into account, and on a cell colony level where cells are assumed to have a predefined geometry that does not change during the process that is being considered. Similar modelling principles can be found in, among others, [4,7,20,37,41,49]. Cellular automata models, which form an important and interesting class of modelling approaches, see for instance [19,34], will not be dealt with in the present paper. Advantages of the small scales is that they are close to in vitro experimental studies in the sense that one is able to directly incorporate measured observations such as cell deformation rates, division rates, and cell migration velocities. One can use observations on one cell only under laboratory circumstances, or one is able to use cell-colony experiments with mono or multi layers. A limitation is the interpretation to clinical cases, where the cell environment is barely known due to its vast chemical and mechanical complexity. This complexity is analogous to the complexity in weather forecasting, and where predictions only face a lim-ited reliability. Despite the claim of many biologists that the dynamics regarding cell differentiation and division (proliferation) is exactly known in cases where the history path of the cell is known, we use probabilistic principles to simulate the aforementioned processes due to lack of knowledge of the history paths of all the cells involved. The larger-scale models on tissues are predominantly based on continuum modelling using systems of partial differential equations, see for instance the studies in [16,24,25], that could make use of the results of the smaller-scale models by the use of averaging principles. The paper will commence with the description of a cell deformation and migration model. Subsequently, we will deal with the migration of cells in colonies. The description of the models will be rather rough since the models were already presented in earlier papers. Furthermore, some modelling results will be presented and described in terms of simulation and biological interpretations. We will end up with a discussion and some conclusions concerning the applicability and limitations of the modelling approaches described here.

Mathematical modelling approaches for cellular and colony modelling
The models that we currently consider are based on cellbased formalisms, where models in which the cell shape evolves over time and models in which the cell-shape is fixed, are considered. The first modelling approach allows more biology in the sense that the actual cell deformation is taken into account. The deformability of the cells, in particular of cells that are active in the immune system response where hazardous pathogens have to be engulfed, is of importance for the efficiency of the performance of the immune system response. Using this modelling approach, one is able to quantify the influence of for instance the cell stiffness on the efficiency to engulf detrimental pathogens. In the case of modelling colonies with large numbers of cells, the inclusion of cell deformation provides a higher model-resolution, which makes the computations more expensive in terms of longer computation times, and here the second modelling approach, where the cell shape is fixed, becomes more attractive from a computational point of view. Both approaches are based on semi-stochastic principles in the sense that random processes are incorporated in the modelling. The inclusion of random processes is motivated by the heterogeneity of the extracellular matrix and the uncertainties in cellular parameters in terms of age, condition and maturation levels (if it comes to differentiation, death or proliferation processes). The main differences in the modelling approaches of this paper concern the scale of modelling. The cell-deformity models are suitable for a lower number of cells, and hence for smaller areas than the modelling approaches where the cell shape is fixed at all times. First a description of a celldeformation model will be given for cell deformation and cell migration. This description is followed by a model description for the colony-based models where the cell-shape is fixed at all times that we used in our studies.

Cell-based modelling
In this approach, one considers individual cells in reasonable detail in the sense that the geometry of the cell is not fixed during the simulation. Here, we model the process using a set of coupled ordinary differential equations, and we present the approach as in [48], where only the cell boundary is divided into mesh points, that are all connected to the central point of the cell, which mimics the cell nucleus, see Fig. 1. The boundary points are allowed to move depending on an external stimulus, such as the gradient of a chemical, to model chemotaxis or haptotaxis. Besides this movement, we also take into account the equilibrium cell shape, which could change in the course of time as a result of cell differentiation to a different phenotype. This tendency towards a cell shape equilibrium models can be physically interpreted as a trick to model the stiffness of the cell where partial differential equations modelling is avoided. This approach allows cheap computations in the sense of computation speeds. For the sake of illustration, we consider the case of a cell that experiences a concentration gradient and which will migrate and deform as a result of this concentration gradient. This case could biologically model chemotaxis or haptotaxis, with applications to the immune system response during wound healing for instance, or the directed movement of white blood cells towards cancer cells as a result of a mechanical or chemical stimulus. Let c(t, x) denote the concentration that makes the points on the cell boundary move. Note this quantity could, without loss of generality, correspond to a mechanical stimulus as well. Here t and x, respectively, denote time Fig. 1 The division of the cell boundary into mesh points and its connection to the neighbours and to the cell centre as in [48] and position within the domain of computation. Further, let X k j (t) denote the position of nodal point j at time t of cell with index k, then movement as a result of a gradient of c is determined by where we assume that M k (t) := {X k j (t)} m k j=1 denotes the set of discrete points on the boundary (being a closed curve or surface in R 2 and R 3 , respectively) of cell with index k, where this cell boundary is divided into m k mesh nodes. The situation is as has been plotted in Fig. 1. The above equation models a cell migrating and deforming in a medium such as an extracellular matrix as a result of a chemical gradient. Several experimental studies evidence that chemical gradients bias the migration of cells, such as [15] or other studies on the Boyden assay [6], where cells are placed on a porous filter with a chemical distribution. Byrne and Owen [8] mathematically investigate the relation between the random cell motility, chemotaxis and chemokinetic cell responses. On the basis of the behaviour of cell migration regarding concentration gradients, we postulate in our model that the position of each point on the cell boundary moves according to the gradient of concentration. This is motivated by the fact that the Keller-Segel model is a special case of the Fokker-Planck equation, where the chemotactic part represents the drift term in a stochastic differential equation used for single particles or cells. We note that the current formulation is based on a quasi-equilibrium hypothesis. In [11], a kinetic equation for the directional migration of individual cells according to a chemical gradient has been taken into account. In Eq. (1), the cell stiffness has not been incorporated yet. The above equation would allow the cell to deform arbitrarily, which is not realistic from a biological point of view since cells contain a large portion of water, which is incompressible. Since cells also contain a considerable solid fraction, which is compressible, cells do possess a limited compressibility. To mimic the cell stiffness, equation (1) is extended as follows where α > 0 can be interpreted a cell-stiffness parameter andX k j (t) denotes the equilibrium position of cell k that is subject to equilibrium phenotype formally dictates mathematical stability of the equilibrium points. Once the positions of the nodal points on the cell boundary, that is X k j (t) have been determined, the cell centre is componentwisely determined using The equilibrium points,X k j (t) on the cell boundary are determined over a time-step as follows: After having displaced the boundary points, we determine the orientation of the cell, which in two dimensions, represents the angle between the line connecting the rear-to the front position of the cell and the horizontal axis, as well as the translation by using the centre point of the cell. This orientation is added as an angular rotation to the reference position of the nodes on the cell that is determined by the equilibrium shape of the cell (being circular, spherical or for instance elliptic). In order to describe the equilibrium shape of the cell, we use a given parametrisation. For the sake of illustration, this gives for a circular cell k with radius R c , central position X k c (t) and rotation θ k (t): In such a way, the cell is inclined to become circular again once the stimulus for migration has disappeared. The α-term could represent a penalisation if the α-value is set to a very large value, since then X k j (t) −→X k j (t) if the components of the concentration gradient stay bounded. For small and intermediate values, the α-term models stability of the equilibrium positions of the nodes on the cell boundary, which implies a stable cell equilibrium shape. This provides us with a nonlinear problem at each time-step, which is solved using an implicit-explicit time-integration method where all nonlinear terms were evaluated at the previous time-step. It is noted that larger values of α make the cell deform less easily. We note that inertial effects have been neglected. In [48], the stiffness parameter is allowed to depend on the volume of the cell by which the amount of compressibility of the cell can be enlarged if desirable. Another option is to change the mathematical problem to an optimality principle where the cells follows an external signal as well as possible under the condition that the volume of the cell is not allowed to change or only allowed to change up to a maximal extent. The lastmentioned option is currently under investigation. The above differential Eq. (2) can be adjusted such that it also accounts for perturbations in the extracellular matrix. In [48], this is done by the incorporation of a Wiener process, which is in line with the formalism adopted by [20], hence Here dW(t) = [dW 1 (t) dW 2 (t)] T denotes a vector-Wiener process, which we defined here for two dimensions for the sake of illustration over time interval dt, The Wiener processes dW j (t) corresponding to the coordinate axes, and they are independent stochastic events that are normally distributed with zero mean and variance of dt, i.e. dW l (t) ∼ N (0, dt). The aforementioned (stochastic) differential equations need the initial positions for the boundary positions. The central position is computed by averaging over all the boundary points. Note that the above implementation of the Wiener Process only models an stochastic aberration of the migration of the nodal points on the cell boundary, and that this formalism cannot be used to incorporate the effect of random walk by the cells since the averaged displacement of the cell will be zero by the nature of the Wiener process. Hence the above σ parameter accounts for a stochastic variations to the cell shape. If one wants to include random cellular motion, which models cellular diffusion, then one should add a √ 2D c dW(t)-term, with D c denoting the diffusion rate of the cell, to the position of the centre of the cell. Next to c, other signals may determine the dynamics of the cellular points, such as hard impingement with other cells. In [48], a potential function is used that was inspired on a Lennard-Jones formulation in molecular dynamics where a repulsive term is taken into account that incorporates an additional term depending on the distance between the current node of the current cell and points on the boundaries of the other cells. The term used in [48], for intercellular repulsion is given by the addition of a phenomenological term Here ε > 0 represents a measure for the repulsive force between two adjacent points from different cells. Having all the boundary points, one determines the cell centre, cell volume/area and boundary area/length easily as post-processing parameters. The details are given in [48]. These parameters also allow the determination of the cell geometry parameters such as the cell shape index. As mentioned earlier, this modelling can be used to simulate the immune system response on infections in wounds, or on the clearance of cancer cells by the immune cells. Note that this model is presented in more detail in [48] and in [51] for the immune system response.

Colony-based modelling
In this section, it is assumed that the geometry of the cells is fixed and determined by the specific phenotype. At this stage, we assume the cells to be spherical or hemi-spherical in the three-and two dimensional cases. For the two-dimensional case, it is assumed that the projection of all cells onto the substrate is circular. This substrate could mimic an extracellular matrix. In [49,50], it is assumed that the cells migrate as a result of communication based on exerting strains on the substrate or the extracellular matrix. One can easily apply a finite-element method to compute the displacements of the substrate or extracellular matrix. In [49], it has been chosen to approximate the displacement on the basis of a decaying exponential function, which is maximal but bounded on the cell centre. The adopted formalism is based on the strain energy density, which has the advantage of being additive. To this extent, using Hooke's Law, the strain energy density is given by here M 0 i denotes the strain energy density exerted by cell i, which is located on position x i (t), which exerts a force F i , on a substrate or extracellular matrix with elasticity E s on the spatial location of cell i. In the present modelling studies, the values of the cell force F i and the extracellular stiffness E s are assumed to be known and hence do not have to be computed as part of the solution. Hence, it is assumed that the cells do not adjust the stiffness of their environment. This assumption is reasonable for epithelial cells and immune cells. For cells like fibroblasts and myofibroblasts that are involved in wound healing or wound contraction, this assumption is questionable since these cells are indeed able to change the structure and mechanical stiffness of their environment. The cell has radius R. Furthermore, σ and ε, respectively, denote the stress and strain of the extracellular matrix. Based on the analytic solution, where we want to incorporate a finite strain energy density at the cell centre, we approximate the attenuation of the strain energy density signal away from the cell by As far as it is known at this stage, Reinhart-King and coworkers [9,39] were the first to experimentally evidence cells communicate through mechanical signalling as well as that there exists a threshold value for the signal under that can be detected by the cells. Using the fact that energy is a scalar physical parameter and therewith additive, the above equation is extended for a colony of cells by Under the assumption that cellular motion is fully deterministic, the displacement vector of cell i over a time interval dt is a linear combination of all unit vectors connecting the cell i with the other cells, of which the contribution of the strain energy density exceeds a certain threshold. To have larger influences to the migration direction from neighbouring cells that give high values of the local strain energy density sensed by cell i, the weight factor from cell j is given by the contribution to the strain energy density as a result of cell j. The strain energy density of each cell-pair determines the weightfactor. Henceforth, cell i migrates in the following direction Terms for which ||x i − x j || = 0 are set equal to zero. We normalise by usinĝ In [49], the displacement of cell i reads as Here the κ-parameter, dimensionally representing the ratio between cell velocity and exerted pressure (hence being proportional to the inverse of a stiffness (therefore a compliance)), is given by where μ denotes the dimensionless resistance parameter of the substrate friction or connectivity index in a phenomenological extent, where f = μF i defines the friction force that the cell experiences from its surroundings as a result of migration. Further γ i stands for the mobility of the cell per unit of time, andF denotes the pulling force exerted by the cell if it is entirely viable. From the above relation, it is clear that the friction and cell pulling forces inhibit the migration of the cell. Furthermore, if the cell is not viable, then F i = 0 and then migration is halted as well. The φ-function accounts for modelling the influence of the chemical environment on cell mobility. The model is enriched with chemotaxis or haptotaxis by adding a term with the gradient of a chemical. Furthermore, diffusion (random walk / Brownian motion) is added through a Wiener Process in all coordinate directions. These additions result into where D c denotes the diffusivity of the cell type under consideration. Here c denotes the concentration of a certain chemo-attractor. This field can be computed using different methods, by finite-element like discretisations or by the use of superpositions of Green's Fundamental Solutions. The latter formalism is characterised by the advantage of having to compute the solution only at these positions where one actually needs the solution values, whereas finite-element discretisations need to compute the solution over the entire domain of computation, and hence large systems of (linear) algebraic equations need to be solved at each time-step. A disadvantage of using the Green's Fundamental Solutions is the need for simplifying assumptions for the model parameters and geometry of the domain of computation, whereas finite-element approaches allow a full geometrical flexibility. Examples, where the approach with Green's Fundamental Solutions has been used in the context of the immune system response in ischemic wounds and cancer, can be found in [51] and [52], respectively. Regarding the stochastic part in the migration of cells, one could use a uniform probability density where a cell is allowed not to move further than over a certain distance as it is done in [27,42,48]. Another possibility is to use a Wiener process, as has been presented here. The introduction of the Wiener process (being a normal distribution with zero mean and the time-interval as the variance) leaves the possibility that a cell can move an arbitrarily large distance within a predefined time-interval although the likelihood tends to zero. This latter consideration could be an undesired side-effect. However, in order to be consistent with the Keller-Segel model for chemotaxis, which is a special case of the Fokker-Planck equation, the modelling approach using the Wiener process has been chosen. Probably the differences in the results from the two modelling approaches will not be spectacularly large. Cells will also possibly come into mechanical contact, which is hard-impingement. To this extent, we use the contact mechanics by Hertz for the contact of two cilinders or spheres in the two-dimensional and three-dimensional case, respectively. For the purpose of modelling this phenomenon, we employ the invagination model that was presented in [17]. The invagination forces are responsible for pulling cells from one-another whenever they are in mechanical contact. We derived in [48] the following contribution to the strain energy density Since this contribution, being nonzero if and only if cells are in physical contact, will cause the cell to move away from the neighbouring cell by the pushing force due to the neighbouring cell (and itself), it has to be subtracted from the contribution to the strain energy density resulting from long-distance communication (making the cells converge). Hence, the effective adjustment to the strain energy density function is then given bŷ whereM i and M i j , respectively, represent the total strain energy density and the contribution to the strain energy density from the elastic interaction between neighbouring cells. Furthermore, the set N i denotes the collection of indices of cells that are in mechanical contact with cell i, that is, the distance is less than 2R. This set can be written as The negative sign results from treating the strain energy density as a potential, where the long-distance and near-by contributions give opposite migration directions and hence contributions to the potential with opposite signs. For more details, we refer to [49,51]. Next to cell migration, cells are also subject to division (proliferation) and death. To this extent, we model these processes using stochastic processes. In [52], we take the cell-cycle into account for the modelling of cell division. A simplified, and less elaborate procedure, is to include a simple exponential distribution for the probability of cell division and death for t > 0 (after the appearance of the cell): Note that at t = 0, the probability that a cell dies or divides within the interval (τ, τ + Δt) is given by f (τ ; p e )Δt for an infinitesimal time-interval. This distribution is 'memoryless', and hence where In [50], the probabilities for cell division and death were assumed to be independent of the infectious agent that is released by bacteria. In [52], we assumed that the probabilities were effected by the strain energy density as a result of contact forces. In the case of tumor cells, the probabilities for division are generally set to higher values to mimic the increased proliferation rate of the tumor cells relative to the other constituent cells.

The numerical methods
In this section, the following components in the numerical implementation are described: time-integration for cell displacement, fundamental solutions for chemical signals, and finite-element approaches for the mechanical balance.

Time integration for cell displacement
In both modelling cases, coupled systems of ordinary differential equations have to be solved. Since analytic solutions are impossible to generate, numerical solution techniques are used based on time-integration methods. In this sense, one better uses explicit time-integration methods instead of implicit methods since the contact mechanics pose a discontinuity in the acceleration once cells that did not have mechanical contact earlier just come into mechanical contact. Implicit methods need nonlinear solution techniques for the solution of the resulting algebraic systems where the object functions may no longer satisfy continuity and smoothness requirements. This loss of continuity is often detrimental to convergence behaviour. Another reason for the use of explicit time-integrators is that upon contact mechanics, one wants to limit the cell migration distance over a time-interval to avoid cells being located such that their centres move over one-another. Despite the fact that explicit numerical time-integration methods are well-known for exhibiting conditional stability and hence limiting the maximum allowable size of the time-step, in the current modelling, the maximum size of the time-step is always limited so that the cellular displacement is of the order of the cell-diameter in any time-integration method that is used. This limitation of the time-step warrants a sufficiently small value for the time-step so that numerical stability is attained for the explicit timeintegration method. Note that this limitation relates for the case that the cell displacement is at most one fourth of the diameter of a cell. Hence the time-step is related to the cell velocity and cell diameter. Runge-Kutta schemes can be used to obtain a high order accuracy if the solution remains smooth, that is, contact forces are not suddenly switched on or off. As soon as contact forces are switched on or off, then, smoothness of higher order derivatives does not apply any longer, which makes the higher-order methods useless under these circumstances. Therefore, here we most often rely on the Euler Forward method for the time-integration of the deterministic part of the equations. Besides the deterministic part, the stochastic part of the stochastic differential equations are integrated using the Maruyama method.
Regarding the cell-deformation and migration model, the nonlinear problem is solved using an Implicit-Explicit (IMEX) method, where all factors that make the problem nonlinear are evaluated at the previous time-step.

Fundamental solutions for chemical concentrations
To determine the concentrations that drive the migration of cells or its points on the boundary, fundamental solutions are used as an alternative to approaches that are based on discretisation techniques like the finite-element method. The chemokines are assumed to be subject to diffusion and regeneration. Regeneration taking place place as a result of the competition between pathogens and constituent cells for oxygen and nutrients is modeled by the use of point sources that coincide with the centres of the pathogens. The same approach is used for modelling the concentration of tumor-cell derived chemokines. Let X S (t) be the point source position (for instance coinciding with a tumor cell or a pathogen) at time t where a cytokine is secreted, then its concentration, c S (t, x) is modeled by the following initial value problem It is assumed that initially the amount of cytokines is zero.
Here D represents the diffusivity of the cytokine, and ν S (t) represents the amount of chemicals that are produced per unit of time by the source. Note that δ represents the Dirac Delta Distribution, which is zero at each point except in the origin, and satisfies δd = 1, if the origin is fully contained in . Since the above partial differential equation is linear, to obtain the solution for multiple point sources (that is multiple pathogens or multiple tumor cells), the concentration, c(t, x), is obtained by Superposition, i.e. by summing over all contributions, to get The solution components c S are obtained as fundamental Green's Functions (Duhamel Principle), which in the case of unbounded domain with dimensionality d read as In the above formulae, so far no correction has been carried out to deal with sources that are appearing due to proliferation or disappear due to engulfment or (programmed) cell death.
In the case of newly appearing cells, one integrates from the time at which the new-born cell becomes active ν S > 0. For dead cells, the source function ν S is set equal to zero, while the history path has still to be taken into account. As time proceeds, the contribution of a dead cell long after its death, is no longer significant, and hence, one can truncate these values.
The main advantage of this approach compared to finiteelement like approaches is that the concentration is obtained immediately only at those points where they are of interest. It is not necessary like in finite-element techniques to evaluate the solution on all the nodes in the domain of computation.
Another advantage of the current approach is that it is not necessary to map the concentration from the finite-element meshes onto the positions of cells where the concentration or its gradient is actually needed. Hence one only evaluates the solution where the solution is actually desired without the use of mapping operations or gradient-recovery techniques to retain the accuracy of the solution. As an alternative to the current approach, based on an unbounded domain, one could use Fundamental solutions based on Fourier expansions for bounded domains. Disadvantages of the current approach are of course the assumptions regarding isotropy and homogeneity of the medium. We used the approach that was presented here in [52] for modelling early stages of tumor development and in [48,50,51] for modelling infections and the immune system response.

Finite-element approach for mechanical balance
In some other works, one uses discretisation techniques to evaluate the concentrations or mechanical displacements. In our studies, we used the finite-element to model mechanical displacements. Here we consider a bounded domain ⊂ R d , where we solve the following partial differential equation that describes mechanical equilibrium: where σ represents the stress tensor and F represents internal forces. Constitutive relations describe the relation between the stress and strain tensors. In our modelling, we used Hook's Law. For each cell that pulls on the surrounding extracellular matrix, we consider pulling forces on the cell boundary. For this purpose, the cell boundary is treated as a polygon or a polyhedron in the respective two-and three dimensional cases. To this extent, the cell boundary is divided into boundary elements and on each of this element a point force is applied with a magnitude of the cell traction per unit area/length of the cell boundary multiplied by the size of its area/length in the direction normal to the cell boundary. The total force F is obtained by summing over all N boundary elements, as in Here P(t, x k ), n(x k ), and k , respectively, denote the cellular normal force per unit of length (or area in the threedimensional case), the unit normal vector pointing into the cell and length (or area) measure of a boundary element. The limit transition N −→ ∞, i.e., Δ = max j Δ j −→ 0, gives if C is a piecewise, smooth curve. In the case of multiple, say N c , cells, one sums again over all cells in the domain of computation, to get Spring force boundary conditions are used. In the finiteelement method, the entire domain of computation, , is divided into mesh points and elements that have the mesh points as vertices. Using the solution for the displacements and strain energy densities, one can determine the differentiation rate of cells, or the displacement of the boundaries of the domain of computation. In our modelling the boundary of the domain coincided with the wound boundary and herewith the important process of contraction in burns could be modelled. The interested reader is able to read [53] to get to know more about the integration of the finite-element method into the cell-based modelling.

Alternative modelling approaches
A widely used alternative modelling strategy is the use of cellular automata formalisms, where the domain of computation is divided into a discrete set of points that form a lattice, where each lattice point is either occupied by a cell or not determined by some mathematical and/or probabilistic processes typically based on optimality principles from (virtual) energy functionals. This approach has been chosen by, among others, Merks et al [34]. Merks predominantly applies this principle to modelling angiogenesis. It is noted that cellular automata (or cellular Potts) models can also be used both on tissue scale for processes like wound healing or the development of tumors. In this scale, each lattice point con-tains an identifier saying whether or not that this point is part of the wound/tumor or the surrounding (undamaged) tissue. The work by [46] provides a Cellular-Potts framework where cell-matrix mechanical interactions relate to collective endothelial durotaxis, in which cells move according to the gradient of extracellular rigidity. The mechanical model that they use, was developed by [28], where cellular tractions are modelled by means of a grid that coincides with the Cellular-Potts lattice. Here each mesh point that is occupied by a cell is equipped with a body-force pointing towards or away from the cell centre. Our approach is different in the sense that we only assign forces on the cell boundary. In this sense, our modelling approach has more similarities with the immersed boundary method that was employed in ICCells by [40], except that in our approach we only consider interactions between the cell boundary and solid extracellular matrix.
Alternatively to the cellular-automata models, the evolution of the cell geometry can also be modeled using partial differential equations where one solves free and moving boundary problems for each cell. Madzvamuse [33] uses visco-elasticity theory for the deformation and displacement of cells. This requires the approximation of the solution of a evolutionary partial equations for the evolution of the mechanical balance. This problem poses a challenging moving surface and interface problem, which contains more physics than in the present approaches, however, on the other hand, the solution to the resulting problem using finiteelement like discretisation becomes more time-consuming.

Numerical simulation results
In this section, we show some simulation results to demonstrate the potentials of the models. First, we deal with cell deformation and migration, and subsequently, we deal with the colony-based models.

Cell-deformation modelling
We start with the deformation of a three-dimensional cell under a chemical signal and subsequently, we consider an application of many cells in parallel applied to the immune system.

The basic one-cell model
We show the deformation of a cell in R 3 under the influence of three chemical sources in Fig. 2. The three red dots mimic pathogen sources, which do not move in this simulation. The three pathogens release a chemical that attracts the cell, being the red surface. The concentration gradient was determined by the superposition of three Green's Functions for each cell. The Green's Functions hold for an infinite domain of computation. We note that one could also have computed the concentrations by the use of a finite-element method. It can be seen that the cell deforms such that three branches appear, which engulf the pathogen sources. After engulfment of the pathogens, the concentration gradient disappears and therewith the cell deforms back to its equilibrium position. This simulation was also presented in [48].

Application to many cells leaving a small blood vessel
The second illustration is shown in Fig. 3, where a small blood vessel is shown containing a number of white blood cells, see the white areas in the vessel. Blood flows through the vessel and the flow is modelled by the use of Poissieulle flow, of which the velocity profile is added to the velocity of all the nodal points of the boundary of all the cells. At a certain time, it is assumed that pathogens appear at a location next to the small blood vessel. These pathogens compete with the constituent cells for oxygen and nutrients and thereby the acidity increases locally. The acidity is represented by the biotic lactates, which diffuse over the entire tissue and eventually they will reach the small blood vessel. In the modelling, the pathogens are considered as migrating point sources for the biotic lactates. These lactates trigger a chemical signal in the small blood vessel, which make the white blood cells move towards the vessel walls, and subsequently, they transmigrate through the vessel and enter the extracellular matrix area around the blood vessel. The white blood cells are assumed to appear randomly at the inflow boundary of the portion of the small blood vessel that we modeled. Proliferation of the white blood cells is not taken into account since white blood cells are produced from multi potent cells in the bone marrow and they are found in the lymphatic system. Furthermore, cell death of the white blood cells is not taken into account in this modelling. Having arrived there, they move, driven by the gradient of the lactate concentration, towards the pathogens and subsequently, they engulf the pathogens. This process can be seen in Fig. 3. We realise that the current picture is too simple to reflect reality. In [50], the parameters have been changed and adopted such that drug treatments can be simulated. We note in addition that in [50], we studied the influence of the vessel wall transmittivity on the white blood cell motility in the simulation results to be able to model diseases that inflict changes of the small vessel structure. This is realised by locally changing the μ and σ -values in Eqs. (2,5,6). Similar adaptations can be made to include cell motility, cell stiffness or leukocytes counts that can be worsened in certain diseases like diabetes, dengue or AIDS.

Colony-based modelling
We start with the closure of a wound under the influence of pathogens that paralyse the constituent (epithelial) cells. Subsequently, we deal with wound contraction by  they disappear and the cell retreats to its equilibrium shape, being a spherical shape in this case. These figures were presented earlier in [48]. (Color figure online) (myo)fibroblasts. We end up with treating the initiation of cancer.

Wound closure under the influence of pathogens
First we consider the closure of a wound under the presence of pathogens that release an acidic chemical, i.e. a lactate that is able to paralyse the constituent cells. In Fig. 4, it can be seen that at the early stage there is a circular gap, where the constituent cells have been indicated by the red circles. As time proceeds, the constituent cells re-arrange, and pathogens proliferate and move randomly over the domain. It is assumed that the pathogens easily enter the cells, which is also observed in experimental studies such as by Topman et al. [44]. The pathogens paralyse the cells by the secretion of the lactates. In these simulations, the concentration of the chemicals is monitored using a superposition of Green's Fundamental solutions. As time proceeds, it can be seen that the gap, which mimics a superficial wound, does not close in time. Furthermore, it can be seen that between the constituent cells, gaps appear in the undamaged region and that these gaps do not vanish as time proceeds due to the immobilisation of the cells as a result of the lactates secreted by the pathogens. The same simulation, not shown here but shown in [51], has been done for the case where there are no pathogens. In the simulation without the pathogens, the gap closes entirely and there are no gaps in the undamaged region. It is noted that despite the fact that the model possesses many stochastic components, one could search for a threshold of the initial amount of pathogens under which the wound gap is still able to fully close. One should use statistical methods here with Monte-Carlo methods where for each parameter value multiple runs have to be performed.

Wound contraction
The next case that we show mimics contraction of a wound. Contraction generally takes place in burns, in particular in deep wounds. The contraction takes place via a sequence of processes. First, the immune system response clears up aggressive chemicals, hazardous pathogens through the white blood cells. Furthermore, the white blood cells release chemicals that trigger the migration of fibroblasts, which are a kind of skin cells, towards the wound area. The fibroblasts produce new collagenous tissue that repair the wound site. Furthermore, the endothelial cells are triggered and migrate towards the wound area, by which re-vascularisation takes place by the re-establishment of a new vascular network. In the present simulations that we show in Fig. 5 the wound area is portrayed by the red area. The white blood cells, not shown in the figure, release a growth factor (VEGF-hormone)  Fig. 4 Snapshots at consecutive times of a simulation of a mono-cell culture with an initial gap. At the centre of the gap, bacteria migrate randomly and divide, and subsequently infect the constituent cells, which are paralysed so that the gap does not close and so that the surrounding initially undamaged tissue gets more microgaps and hence a less regular cell density and hence the quality becomes worse. These results were also given in [50]. (Color figure online) that triggers the surrounding fibroblast to migrate into the wound area. The fibroblasts are visualised by the white circles. As the fibroblasts move into the wound area, they pull the tissue and under unfavourable chemical and mechanical circumstances they are able to differentiate to myofibroblasts, which are cells that pull at a larger force on the surrounding tissue, produce more, excessive collagen, and furthermore shorten the chemical bonds of the long polymeric chains that constitute the collagen. These chemically and mechanically unfavourable conditions typically refer to the concentration of certain cytokines (growth factors) and/or to pulling forces that the fibroblasts experience, which make them differentiate to myofibroblasts, see [45] for instance for a modelling study. It can be seen that the wound area decreases over time and that the wound region changes shape from the initial rectangular shape to a more star-shape due to the internal forces exerted by the (myo)fibroblasts. The myofibroblasts are represented by the white diamonds. Next to the snapshots at different times in Fig. 5, we show the evolution of the wound area in the contraction phase over time for various values of the differentiation rate of the fibroblasts in Fig. 6. It can be seen that the model is capable of predicting the impact of the differentiation rate on the contraction behaviour. It can be seen that wounds where there is no differentiation recover to the initial area. This is beneficial in the case of burns since a final contraction, referred to as a contracture, is an undesired side-effect. The contracture of the skin in general reduces the mobility of the patient. If the fibroblasts do differentiate to a large extent to myofibroblasts, then contraction becomes more severe and moreover the contraction becomes partly irreversible and hence a contracture remains.
This second model for the contraction of a burn is based on cell-based modelling with semi-stochastic cell migration incorporating random walk (diffusion) and chemotaxis. The mechanical forces exerted by the cells have to be monitored on the boundaries of the wound and therefore, finite-element simulations have been used. The details are presented in [53].

Initiation of cancer
Finally, some simulations are shown for the case of the initiation of cancer. We consider a spherical tissue domain with epithelial cells, which divide and die. During the division process, errors in DNA copying can occur, by which tumor cells appear. The tumor cells divide and press onto the epithelial cells, that constitute dermal tissue but could also constitute lung tissue. The contact forces onto the epithelial cells make them more prone to die and less likely to divide and therewith the tissue deteriorates. The tumor cells are Fig. 6 The wound area versus time for several differentiation rates from fibroblasts to myofibroblasts, see also in [53] allowed to divide at a high proliferation rate. The details are given in [52]. From the surrounding vascular network, white blood cells enter the tissue area and these white blood cells chase the tumor cells on the basis of chemical signals that are released by the tumor cells. The white blood cells have the ability to neutralise the tumor cells to a certain extent. Some snapshots of a run are shown in Fig. 7. The green spheres denote the constituent epithelial cells, the blue spheres represent the immune cells, and the red cells visualise the tumor cells. More details are presented in [52]. Finally we show the number of cells as a function of time for a specific run where it is clear that the number of cancer cells grows at the expense of the epithelial cells. In this run, the T-cells, i.e. the white blood cells, are not able to beat the cancer cells. Furthermore, in the same figure, the fraction of the tumor cells is plotted over time, which also clearly shows that the tumor is going to develop further.

Discussion and conclusions
First the relation to other modelling studies in literature will be assessed and subsequently the current modelling approaches with possible extensions will be discussed.

Relation to other studies in literature
The methods that have been described in the present paper are based on cellular models where a distinction between models allowing for cell-shape evolution and cell colony models where the cells do not change their shape has been made. In literature, many other studies can be found concerning cellshape evolution, such as in the works by [4,21,29,33,36,47]. In [21], the cell moves through a complicated domain with obstacles as a result of a chemical stimulus. The modelling approach used there is based on the representation of the cell by a set of connected nodes on which protrusion, cortical tension, friction and a force to keep the cell volume constant work. In this sense, the approach differs from the present continuous approach. The modelling of the cell through a domain with obstacles is very interesting and in our own modelling approaches, we also have some preliminary results where cells even have to adjust their geometry in order to go through an obstacle. Another cell-shape evolutionary model is described in [36,37], where a set of reaction-transport equations is solved over the cell boundary using an ALE formulation within a finite-element method with linear elements. The actual movement of the cell boundary is determined using a level-set method. Besides this cellular movement, a partial differential equation for the chemical stimulus is solved outside the cell and coupled to the obtained receptor distribution from the boundary finiteelement method. The method differs from the approaches in the current paper in the sense that the present methods are simpler and based on Fundamental solutions as well as on a simpler representation of cell boundary movement. Further, our method does not require the discretisation of the partial differential equations for the level-set method over the entire domain of computation. Further, no partial differential equations have to be solved over the cell boundary in our methods opposed to the studies in [36,37]. The model in [47] treats migration of cells in larger quantities with applications to fibroblasts. The model is based on traffic rules that were observed experimentally in collective cell movement. The modelling physics is elegant and experimentally-based and is based on constant magnitude forces that are exerted from the pseudopods radially away from the nucleus. Furthermore, pseudopod formation is influenced by the current directional motion as well as by the space-time profile of the chemokine concentration which is released by all the cells. In this modelling work, cell collisions are rare and give rise to randomisation of the distribution of the pseudopods. The study in [29] models, like in [21], cells as a collection of particles including material particles like plasma membrane and cortical biomolecules. All particles move resulting from Newton's Law of motion under conservative, dissipative, actomyosin forces and membrane tension forces. Further the extracellular matrix is incorporated in terms of particles as well, which directly facilitates intracellular communication. The equations of motion are solved using the Verlet algorithm for the time-integration for the particle position in the two-dimensional domain. The main extension and result of the model is that cellular rotation is taken into account, for which the cause was not revealed earlier.
In our simulations, results are obtained with small computation times, although we sacrifice on the physics behind the cell migration phenomena. On the other hand, the parameter space in our models is much smaller, which also poses some merits for experimental validation and regression techniques. This is the clearest advantage of the current approach that is shared by the approach in [47]. The model in [47] is very interesting in the sense that it is also based on first principles and that it is probably also cheap in terms of computation time. The only shortcoming may be its applicability to high cell densities where cells collide often.
Next, we review some models without cell-shape changes over time. In [13], the dynamics of tumor growth is considered based on modelling the interaction between single cells and the distribution of nutrients and biomechanics forces. All cells are individual objects which is parametrised by biophysical and cell-kinetic parameters. The simulations show that the tumor exhibits sub-exponential growth at large sizes and that at larger sizes the necrotic core develops as a result of depletion of oxygen or glucose. The individual cells are allowed to divide and to die, and they displace each other by acting pushing pressures on each other whenever they divide. The (in)active mechanical displacement is taken into account by the use of the Metropolis algorithm. The diffusion-reaction equation for oxygen/glucose is solved using a square discretisation mesh, where the reaction represents the consumption rate by each cell. A mapping of the solution onto the cellular centres is needed to compute the concentration. In [22], similar modelling techniques are considered regarding the solution of concentration equations on a rectangular discretisation mesh and forces of the cells on the substrate or extracellular matrix. In the last-mentioned article the simulation toolbox is described. In the study of [23], cell alignment along micro vessels is modelled. The modelling approach is roughly based on the previous two approaches, however, now the daughter cells resulting after division align according to the closest sinusoid. The model is also tested against experiments on mice and the model is able to explain the role of hepatocyte-sinusoid-alignment on the restoration of the restored liver architecture. The study in [10] deals with tumor-stromal interactions that stimulate the migration of cancer cells where cell polarisation is taken into account. Besides the polarisation, the heterogeneities, which also enhance the migration of the cancer cells, are taken into account. Strandkvist et al. [42] model cell sorting based on local variations in cell motility. The modelling approaches are similar to the those employed in [49], where contact forces are applied if cells approach each other by distances shorter than the equilibrium distance from hexagonal packing. Random walk is incorporated by picking velocities from a uniform distribution and translated to a diffusion coefficient, which is an approximation from a probabilistic point of view. This principle is used to estimate diffusivities of larger cell colonies consisting of several cell types with various motilities as well as the direct multi-cell approach is used to estimate the effective diffusivity.
All these approaches are very interesting from a modelling point of view. The methods that have been presented in the present manuscript are also able to carry out most of the same biological processes such as alignment according to sinusoids or considering multiple cell types with varying cell motility after some changes in the implementation. The polarisation issue described by [10], is currently dealt with in a study where we want to simulate development of muscle and fat along with differentiation of mesenchymal stem cells to myocytes and adipocytes. Furthermore, the approach of tumor cells where the necrotic core appears could be dealt with in the present approach by considering (point) sinks coinciding with the tumor cells that consume oxygen or glucose. Another study that has just been carried out in our department involves modelling angiogenesis using chemi-cal signalling and stalk and tip cells, see [3], which will be submitted as a future journal paper. This approach could be combined with modelling the growth of a tumor which possibly halts due to depletion of oxygen and continues to grow as soon as oxygen levels increase due to the arrival of small blood vessels where angiogenesis is triggered by chemokines that are released by death of tumor cells in the necrotic core.
Here the modelling could largely benefit from the principles outlined especially in [13,22,23].
A cellular-Potts model (cellular automata) is presented in [35], where cells are represented as sets of neighbouring lattice points in which lattice points may transform to a certain cell at a probability defined from minimisation of a virtual energy. The approach allows for inclusion of chemotaxis as well as mechanical cell-structure interactions. Comparative remarks have been given in [53]. The model is applied to plant physiology but also extended to modelling angiogenesis. The process of angiogenesis contains some similarities with development of dendrites in plants. The modelling of angiogenesis is easily captured by the cellular Potts model, and somewhat harder by the presented models in this paper. However, as mentioned earlier, we do have the first results using the cell-based modelling of this paper for angiogenesis. Although more biological detail is dealt with in the cellular automata models, the parameter space tends to be much larger.
In [40] an immersed boundary method is used within the IBCell computational framework developed by Rejniak to simulate cell migration in fluid structures in two spatial dimensions. This issue has not been dealt with in the present modeling, except for the simplified approach where white blood cells move through a small blood vessel according to Poissieulle flow. In the current approach, the cells did not have any influence on the fluid flow. In [40], the cells do influence the fluid pattern. It would take the solution of the Navier-Stokes equations including cell interactions by the immersed boundary method. The study is performed to simulate the development of carcinoma. The computation time would be enlarged, though new interesting effects in the case of fluidic environments can be dealt with which cannot be simulated in the present modelling approaches. In [14] the packing geometry has been modeled under the influence of mechanical force and proliferation. Cells are modeled as polygons and a stable packing is modelled, where also the cell-shape is part of the solution. The simulations are based on a vertex model where cell elasticity and junctional forces are taken into account. A force balance is reached from a minimisation principle. It would be interesting to add dynamic processes to their modelling such that one can track the cell network configuration in terms of packing under cell proliferation, death and possible differentiation. In [38], a subcellular element model is used to simulate multicellular structures. The main principle is that the cells are represented as a set of, unlike in the cellular Potts modelling, dynamically, moving discrete particles that are connected to one-another. In this sense, the approach by [38] possesses some similarities with the approaches by [21] and [29]. Further, these particles interact with the particles that correspond to other cells, and in the computation of the positions of all the particles, an equation of motion is solved that contains intra and inter cellular communication (in terms of the gradient of a Morse potential). This approach allows to follow the evolution of cell geometry in the course of time. In the simulation results, cobble-stone structures are obtained, which also have been obtained in our modelling approach, see Figures 8,9 and 11 in [48]. Finally, we report the study in [54] where a discrete element model is used to simulate directional collective cell migration. The model takes into account the cell orientation, mechanical contact forces if cells are in physical contact, coattraction forces, rotational turning, self-propulsion, contact damping and contact polarisation. The cells are represented by elastic spheres. In order to be able to deal with larger cell numbers, the simulation tools have been implemented in CUDA, which allows to employ GPU's that facilitate mass parallel computation on many small processor subunits.

The current modelling approaches and extensions
A mathematical cell-based framework for the simulation of several subprocesses in wound closure, wound contraction and tumor initiation has been developed. Cells are allowed to deform, to migrate and to divide, to differentiate to other phenotypes. Using the formalism, one is able to model migration processes like haptotaxis, chemotaxis or mechanotaxis which are all migration mechanisms based on chemical signalling such as concentration gradients or on mechanical signals. Tensotaxis, defined in the sense of the collective cell migration towards the lowest point of an applied stretching force as in [5] will probably need a more careful assessment before it can be treated with similar approaches. A possible way to do this is to compute the stress-strain distribution and the strain energy density in the assay used in [5] and subsequently check whether the point of minimal stretching force is reached by following the gradient of the strain energy density. To monitor the signals, one can use Green's Fundamental Solutions or finite-element like discretisation techniques. Where the discretisation techniques require the signals to be computed everywhere, also at locations where the signal is not needed, the fundamental solutions can be determined on these locations where we want to have them, such as on the locations on cell boundaries or on cell centres depending on the model that we consider. The Green's Fundamental solutions have the demerit that they only hold for very simplified circumstances such as linearity of the partial differential equations, physical constants being con- Fig. 8 The number of cells regarding the various phenotypes as a function of time and the fraction of tumor cells versus time. These results were also given in [52] stant and not allowing boundaries to influence the solution depending on the nature of the fundamental solution that is employed. Chemokines that satisfy transient diffusion equations involve Green's functions that need time-integration of the magnitude and location of the sources of the chemicals that coincide with pathogen positions. This time-integration can nevertheless be a toilsome operation in the algorithm and therefore this issue certainly needs more investigation.
Considerable advantages of the modelling approaches are that the models are close to experimental settings in the sense that the model parameters involve directly measurable parameters such as cell velocities, doubling rates for the proliferation, death rates and differentiation rates. The sensitivity of migration and deformation with respect to chemotaxis is also taken into account. A real challenge is the translation of one scale into the other scale. For instance, knowing that the cell with a certain stiffness is able to deform such that its effective radius from the cell centre in which pathogens can be engulfed is actually larger than the actual cell radius. The difference between the effective cell radius and the actual cell radius is larger if the cell is less stiff. This issue, which can be quantified using the cell-based model, can be incorporated into the cell-colony models for many of the subprocesses where chemotaxis and haptotaxis play a role. We also note that the present modelling approaches are very suitable for visualisation purposes and hence they are suitable for educational purposes to physicians, students, and interested patients.
Though the first results and simulations seem to be promising, real interesting in vivo cases, such as very large colonies of epithelial cells in lungs or very large communities of cells in the region of a burn, cannot yet be simulated due to a lack on computational power. The simulations would greatly benefit from CPU/CUDA implementations, which is a current topic of investigation. This step is needed in a parallel computational environment as a result of the stochastic nature of several components in the models. The stochastic nature necessitates running multiple simulations in order to be able to determine intervals of confidence such that statistically sound claims can be issued. Some preliminary steps in this direction have been taken in the papers that describe our studies. It is referred to [54], where such implementation has turned out to be successful for the discrete element model applied to the simulation of migration of cells. Next to an increase of computational power such that more interesting three-dimensional simulations can be done, one could invest in increasing the complexity of the models in terms of more accurately modelling the underlying biology. Despite the more biological detail that can be dealt with, the model would suffer from a higher complexity in the sense that even more parameters that are hard to get in terms of measurements or regression techniques would be part of the model. It is questionable whether this extension would increase the applicability of the models.
We finally note that the current paper describes a compilation serving to introduce our modelling work to researchers with interest in the field of particle-based methods as well as in bio-medical applications. The interested reader can find more information in our references.