## Abstract

Traumatic and chronic wounds are a considerable medical challenge that affects many populations and their treatment is a monetary and time-consuming burden in an ageing society to the medical systems. Because wounds are very common and their treatment is so costly, approaches to reveal the responses of a specific wound type to different medical procedures and treatments could accelerate healing and reduce patient suffering. The effects of treatments can be forecast using mathematical modelling that has the predictive power to quantify the effects of induced changes to the wound-healing process. Wound healing involves a diverse and complex combination of biophysical and biomechanical processes. We review a wide variety of contemporary approaches of mathematical modelling of gap closure and wound-healing-related processes, such as angiogenesis. We provide examples of the understanding and insights that may be garnered using those models, and how those relate to experimental evidence. Mathematical modelling-based simulations can provide an important visualization tool that can be used for illustrational purposes for physicians, patients and researchers.

## 1. Introduction

Wounds, induced either by trauma (e.g. surgery or burn) or by chronic diseases (e.g. diabetic ulcers or immobility-related pressure injuries), continue to be a considerable medical challenge. They seriously not only compromise the quality of life of affected individuals, but, sometimes, also endanger life. Pressure injuries, for example, are a huge problem, with prevalence that is estimated to be 5.7–12.8% in hospitals, 6.2–12.8% in surgical intensive care units and as high as 25–33% in community settings [1,2]. The cost of each serious pressure injury case is estimated at over $10 000. Because wounds are very common and their treatment is so costly, approaches to reveal the responses of a specific wound type to different medical procedures and treatments could accelerate healing and reduce patient suffering. The effects of treatments can be forecast using mathematical modelling that has the predictive power to quantify the effects of induced changes to the wound-healing process. Wound healing involves a complex combination of fundamental processes, including cell migration, cell proliferation and apoptosis, cell differentiation, and the eventual regeneration of new tissue among many others. Different types of wounds involve diverse biological mechanisms and should be modelled accordingly. One can distinguish between small superficial wounds, rapidly formed deep wounds (such as third-degree burns) which involve wound contraction, and prolonged damage involved in pressure injuries.

Wounds are repaired by motile cells (e.g. fibroblasts or fibroblast-like cells) that migrate into the wound site to close the formed gaps. Cell migration includes various modes and involves changes in the cell structure and morphology, its alignment, mechanics and applied forces for cell–cell and cell–matrix interactions [3–6]. The recruitment of cells to the repair site and their migration to that site determine the extent and rate of the repair process, in filling tissue gaps by cells in the short time scale, and via synthesis of new extracellular matrix (ECM) and its density, directionality and overall structural organization at longer times. Consequently, evaluating the parameters controlling and affecting the dynamics of migrating cell populations (*en masse* migration) during damage repair is critical as a scientific basis for improving clinical wound treatment and management.

Different aspects of cell migration in the context of wound healing have been studied experimentally *in vivo* and *in vitro*. To link those diverse experimental, biological results with proposed theories one may formulate hypotheses in a quantitative approach. The quantification of the biological subprocesses then provides the backbone for mathematical formalisms that may be used to simulate (parts of) the wound-healing process. Mathematical-modelling-based simulations can provide an important visualization tool that can be used for illustrational purposes for physicians, patients and researchers.

Here, we provide a review of contemporary approaches to mathematically model and simulate gap closure in the context of wound healing. We focus on two-dimensional settings, such as superficial wounds, and also provide perspectives for three-dimensional models, which are also much more computationally expensive, and extensions that have been and can be developed; two-dimensional wounds can, for example, be on the skin or on the retina. We have worked on a breadth of problems related to modelling wound healing. In this review, we highlight our work through a set of problems and approaches that we believe would be interesting and useful for the community. To do so, we have focused on three different types of wound models, namely models related to wound closure, healing and contraction.

We begin by summarizing examples of approaches to experimentally obtain quantitative measures of gap closure kinematics, where such parameters serve as the basis for model design. We then provide a wide review of the different approaches that we used to mathematically model wound healing. We start from general approaches looking at the gap closure process from a bird's eye view of a migrating edge, then adding complexity by looking at the motion of the individual cells within the migrating monolayer, and eventually adding mechanics (deformability) of those cells. Thus, we provide a review of the different models that can be used and what may be learned from each of them.

## 2. Experiments of gap closure

### 2.1. *In vitro* experiments

Clinical tests are subject to an enormous complexity of biological subprocesses, and because understanding of biomedical processes is often obtained by the systematic isolation of biophysical mechanisms, well-designed *in vitro* experiments are desirable. The *in vitro* experiments to monitor gap closure, either in response to injury or as part of normal development provide specific kinematic parameters that can be used in the mathematical models. The work by Anon *et al.* [7] experimentally considers the closure of unwounded epithelial gaps through collective cell migration by cell crawling. In contrast, to mimic wound healing *in vitro*, Gefen and co-workers have experimentally studied *en masse* cell migration after inflicting localized damage to monolayer cultures of fibroblasts, myoblasts and pre-adipocytes (3T3, 3T3-L1 and C2C12, respectively, all cell lines from ATCC) [8,9]. The monolayers were damaged in a robust approach, generating reproducible areas denuded of cells, where the migration was monitored following removal of cell debris; debris and dead cells may affect the migration in gap closure.

In the works by Gefen *et al.*, the damage site was monitored using time-lapse phase contrast microscopy, where the kinematic parameters of the *en masse* migration process post-damage infliction were determined by image processing. Using a custom image-processing-based algorithm, they continuously determined the time-dependent gap area, through local texture homogeneity measures. The damage area versus time plots were fitted with Richards-type asymmetrical sigmoid functions, which provided kinematic parameters of the *en masse* cell migration. Specifically, the Richards functions provided, for example, the maximum and average migration rates; the time for onset of mass cell migration, when 10% of the damage area is covered by the cells; and the time for end of mass cell migration, when 95% of that damage area is covered. In the context of wound healing, the above methodology was shown to be useful in quantifying the effects of ischemic factors such as low glucose, acidic pH and low temperature on cell migration in multiple relevant cell types (fibroblasts, pre-adipocytes and myoblasts) as related to wound healing [8], and was further used for testing the efficacy of natural medications used in different cultures for wound healing, such as curcumin, aloe vera and ginger [9].

Recent *in vitro* experimental studies have shown that cells, especially fibroblasts, respond to mechanical cues from their environment. For example, Moreno-Arotzena *et al.* [10] have observed that fibroblast migration in three-dimensional hydrogels is controlled by haptotaxis with a key role for myosin II; haptotaxis is migration up a biochemical gradient of adhesion. Similarly, Ladoux & Nicholas [11] have shown that fibroblasts under specific conditions preferentially migrate towards stiffer regions of the substrate, in a process known as durotaxis.

### 2.2. *In vivo* experiments

Numerous *in vivo* and clinical experimental studies on wound healing or on specific aspects of healing have been reported. In these studies, different therapies to enhance wound healing were investigated experimentally. We list just a few of them here. In the work by Fontana *et al.* [12], the enhancement of wound healing through inclusion of modified porous silicon microparticles is investigated; the microparticles increase cell proliferation and thus aid in healing. Botticelli & Berglundh [13] have investigated the influence of biomaterials on the closure of hard tissue (bony) defects, adjacent to implants in a canine animal model. Healing of burns, which also form an important class of wounds, was investigated by Knox *et al.* [14], including the effect of fibres and dye degradation products on improving wound healing.

## 3. Mathematical models in wound healing

The mathematical models for wound healing that exist in the literature range from very simple descriptions to formalisms with a very high degree of complexity. The simple models have the advantage of being very tractable with only a few parameters that are needed as input for the simulations, and herewith regression with respect to experimental outcomes is straightforward. In contrast, the more complicated models incorporate more biological detail making them more realistic; however, in parallel, those models are frequently intractable. In many cases, it is difficult or even impossible to obtain realistic values for the parameters that are required for simulations, because often those are inseparable from other parameters, are not specifically measured, or are available only with large ambiguity. Thus, one of the greatest challenges is to simplify the models as much as possible in such a way that the relevant cell and process characteristics are simulated adequately with the proposed formalism.

The mathematical models described here range from (systems of) partial differential equations (PDEs) to cellular models where individual cells are treated. The cellular models can even be focused to the microscale, where single-cell deformation is considered or processes are modelled at a subcellular level. We also touch upon hybrid models where the solution of PDEs provides the mechanical and chemical signals in individually modelled cells. We review some of the existing mathematical formalisms, where we start with the class of models that is only based on PDEs.

### 3.1. Models based on partial differential equations

One of the simplest formalisms that has been used is the model that was originally proposed by Adam [15]. In this model, a linear diffusion–reaction equation was formulated for a generic growth factor that was secreted by the cells located near the rim of the two-dimensional wound. This cell layer at the rim of the wound is therefore referred to as the *active layer*. The wound was assumed to start healing (closing the gap) when and if the concentration of the growth factor on its rim exceeds a certain threshold value. The PDE reads as in equation (3.1)
3.1where *c*, *D, λ* and *P* represent, respectively, the concentration of the generic growth factor, its diffusivity, the decay rate coefficient and the regeneration rate in the active layer around the rim of the wound. If on the rim of the wound, then the wound will start healing. For this model, it is possible to derive analytic expressions for the concentration and therewith the analysis of this model is very straightforward. Vermolen *et al.* [16] extended this model to dynamically model the wound closure of superficial (two-dimensional) wounds or wounds on a retina, by replacing equation (3.1) with its transient counterpart so that *c* *=* *c*(*t, x, y*) and by adding a dynamic equation of motion. The added equation of motion, equation (3.2), provides the time, *t*, and position, (*x, y*) dependence of the normal component of the speed of motion of the rim of the wound, *v*. Equation (3.2) depends on the local curvature and on whether the threshold concentration of the growth factor is exceeded at each specific location on the rim of the wound,
3.2where *A* and *B* are rate constants, and *κ* and *H* denote, respectively, the local curvature (at time *t* and two-dimensional position (*x, y*)) and a Heaviside function. This model was used for simulating soft-tissue wound healing as well as bone fracture healing. Important features of this model are its simplicity, and the fact that it is a free boundary problem where the interface between damaged and undamaged tissue is tracked as a part of the solution. The resulting free boundary model equations can be solved by techniques for moving boundary problems, such as moving mesh methods [17–19] or level-set methods [20–22]. Several years before the construction of Adam's model, Sherratt & Murray [23] constructed a reaction–diffusion model for epidermal closure, with PDEs based on the density of epithelial cells and a generic growth factor that serves either as an inhibitor or as an activator; the authors analysed the model using travelling wave solutions and finite difference methods for polar coordinates. This model did not, however, include the explicit tracking of the interface between undamaged and damaged tissue as in the previous model.

The structure and forces at a wound's edge and the characteristics of the cells at that location strongly affect the migration during gap closure, and different models were developed to account for aspects of this. An important biological mechanism that occurs during wound healing is angiogenesis, the process of sprout formation on pre-existing blood vessels, which requires restructuring of the vessel's edge. This was modelled by Gaffney *et al.* [24], Maggelakis [25,26] and Vermolen & Javierre [27] using PDE frameworks. The first two groups modelled angiogenesis as an isolated process where cell densities are used. Specifically, Gaffney *et al.* [24] distinguished between the cell densities at the tips and the sprouts of the capillary network using a reaction–diffusion model augmented with a snail-trail formalism for the simulation of the sprouts that follow the tips. Later on, Flegg *et al.* [28] used a similar model for angiogenesis and wound healing, focusing also on experimental validation of the modelling. Maggelakis [25,26] modelled angiogenesis in conjunction with hypoxia, where the secretion of vascular endothelial growth factor by macrophages under hypoxic conditions is modelled as the trigger for formation and migration of the capillary network at the wound region using a reaction–diffusion model. The third group, Vermolen & Javierre [27], coupled angiogenesis to other cellular processes such as structural contraction at the wound site and migration during wound closure. The last-mentioned model is based on a system of reaction–diffusion equations augmented with passive convection and linear elasticity to incorporate strains and stresses to represent the contraction of the wound. We also note the interesting work by Carlier *et al.* [29,30] in the context of angiogenesis coupled with healing of bone fractures. In the aforementioned work, a diffusion–reaction model is used in combination with a lattice-based model for angiogenesis. In [29], four different computational strategies were assessed for various kinds of bone fractures, focusing on many cases where fractures did not heal because of hypoxia.

The issue of wound contraction has also gained reasonable attention in models. Wound contraction results from the pulling forces that are exerted by adherent fibroblasts and myofibroblasts. Wound contraction is characteristic in deep cutaneous wounds and in burns, and is typically a defensive mechanism to minimize the wound area and reduce the influx of hazardous chemicals and pathogens into the organism. However, in some cases, such as in burns, wound contraction is undesirable, and may even cause permanent contractures, resulting in loss of patient mobility. Very relevant to this issue are the modern studies by Green *et al.* [31], Javierre *et al.* [32] and Valero *et al.* [33] that present modelling of wound contraction in the dermis, where fibroblasts differentiate to myofibroblasts that then pull harder on their extracellular environment. In the study by Valero *et al.* [33], anisotropy of the dermis was taken into account, which, to the best of our knowledge, has never been done in this context. Further, Valero *et al.* [34] present the main challenges in mechanistic models in soft tissues. Green *et al.* [31] analyse the function that characterizes the contraction forces exerted by the cells in an *in vitro* gel. They present one-dimensional simulations as well as a rigorous mathematical analysis of the well-posedness (existence, uniqueness and continuity) of solutions to the PDEs. Recently, Koppenol *et al.* [35] have considered the formation of excess scar tissue (hypertrophic scar) after the occurrence of a burn, where the dynamics of fibroblasts and myofibroblasts has been taken into account. In that work, the tissue is mechanically represented as a neo-Hookean material, and the migration of cells is obtained by solving Keller–Segel-like equations with reactive terms, and hence these equations fall within the category of transport (predominantly convective-like and diffusion-type)-reaction equations. Examples of chemotaxis Keller–Segel-like systems can be found in the study by Byrne & Owen [36], for the various cell types involved. Nardini *et al.* [37] consider a nonlinear diffusion equation for modelling of wound healing, where the leading edge is followed over time. A square-root-like behaviour is shown in the simulations, which is, in a time-region after possible transients, compared with experimental outcomes. The authors quantify cell–cell adhesion and in their model backward diffusion occurs if the amount of cell–cell adhesion is large compared with diffusion.

Two recent reviews discuss continuous-scale models. Tartarini *et al.* [38] present a review on continuous-scale models, as well as models of cell populations with several cell phenotypes with an application to stem cell therapies in wounds. Jorgensen & Sanders [39] present a review about continuum-based mathematical models for consecutive phases during wound healing, where they emphasize that, for instance, models combining wound closure with suturing (stitches for wound closure) techniques are lacking in the present literature.

### 3.2. Cell-based models

Next to the PDE-based models to simulate processes that are related to wound healing, one can distinguish several smaller scale, cell-based models, which treat cells as individual entities. These models can be divided into *cellular automata* (or in particular *the cellular Potts models* as an important subclass) [34–39] or into the so-called *particle methods* [40–47]*.* The cellular automata models are often used in various physical and biological applications, such as wound healing and angiogenesis. The cellular automata models include a lattice in which each lattice point contains a flag denoting whether that location is occupied by a single cell or other element and what is the element's state. The elements and their states may, for example, be blood vessels versus extracellular matrix in a biological context, indicated as wounded or unwounded, or in a metallurgical context may denote phase transformations, such as liquid versus solid; such models were applied to phase transformations in metals and alloys. This class of models often contains a physical or biological law that dictates the change of state at a specific lattice point. In some cases, the change of state is determined by a random process where its probability at a specific lattice point is affected by various parameters, such as chemical concentrations or the state of the neighbouring lattice points. Some examples of cellular automata models are the studies by Chaplain & Anderson [40] and by Hatzikirou & Deutsch [48]. Chaplain & Anderson [40] developed a deterministic cellular automata model to track individual endothelial cells during tumour-related angiogenesis; several experimental set-ups were also suggested in this work. Hatzikirou & Deutsch [48] use a so-called Fermi-gas cellular automata approach to represent cell migration in heterogeneous environments. In their formulation, they incorporate the feedback between cells and their environment.

Some cellular automata models involve a state of change determined by energy considerations. These energetic principles are the basis of the cellular Potts models used by Merks & Koolwijk [49], Merks *et al.* [50], Van Oers *et al.* [41] in the context of angiogenesis. The cellular Potts models were introduced to the angiogenesis modelling community by Graner & Glazier [42], and the models originate from theoretical physics. Even cell deformation processes are modelled using this approach and combinations with solutions of PDEs describing mechanical deformation and chemical content have been developed, as in [41].

The other small-scale modelling formalisms are the so-called *particle methods*. Vermolen & Gefen [43] formulated a model based on the strain–energy density field where cells are attracted to each other yet repel if they begin to overlap. Overlapping is prevented by forces, modelled by the Hertz contact model, which was used for the first time in this context by Gefen [44] to model invagination of viruses into cells. The model was based on the *in vitro* observations by Reinhart-King *et al.* [45,46], where it was found that individual cells communicate remotely with each other by mechanical signals. Another cell-based model based on mechanical signals was formulated by Rey & Garcia-Aznar [47]. The equation of motion of the attracting–repelling cells roughly reads as
3.3where ** x**(

*t*)

*,*(

**v***t,*(

**x***t*))

*, σ,*(

**W***t*)

*,*respectively, denote the cell's position over time

*t*, the deterministic component of the cell migration speed (by for instance chemotaxis, haptotaxis, durotaxis etc.), the motion owing to random walk, where (in which

*D*is the cell types diffusivity), and the vector Wiener process, in which each component is distributed according to the normal distribution with zero mean and variance of d

*t*. Other cell-based models were developed by Byrne & Drasdo [51], Drasdo & Höhme [52], Höhme & Drasdo [53], Vermolen

*et al.*[54] and Schluter

*et al.*[55], where applications beyond wound healing (such as cancer) were dealt with.

In the context of wound healing, Vermolen *et al.* [56] extended the cell-based model to incorporate the effects of bacterial infections and the response of the immune system. Figure 1 shows an *in silico* simulation of the process of gap closure in a monolayer of cells; the simulation is performed in two dimensions for simplicity, yet the cells are rounded as is typical in a three-dimensional environment, thus the simulation could also represent a single slice of a three-dimensional gap. The model run represents a superficial wound with a blood vessel network underneath. The simulation includes vital, tissue-constituent cells, such as epithelial cells, pathogens such as bacteria that accumulate at the gap, immune cells that migrate to the gap site and blood venules in the tissue that facilitate immune cell arrival. In the simulation, Vermolen *et al.* have taken into account that the constituent cells' motion will be inhibited as a result of the biotic lactates that are released from the pathogens, and also incorporate the process of engulfment of the pathogens by the immune cells (figure 1*a*). Using this model, the impact of the pathogen proliferation rate on the gap closure rate was quantified (figure 1*b*); the effects of other parameters can similarly be evaluated. Similarly, it was found that under the conditions that were employed in the modelling, the immune system operated optimally if the immune cells had a lifespan of approximately 7 h [56].

Vermolen & Gefen [57] and Boon *et al.* [58] have further extended the cell-based models to the simulation of wound contractures that occur in deep burns. The models were modified by combining mechanical balances with the cell-based models. In addition, differentiation of fibroblasts to myofibroblasts was added. These models that combine cell-based formulations with PDEs for chemicals or for mechanical signals fall within the class of the so-called *hybrid models*. These models can further be extended to simulate, for example, effects of externally applied forces that change the material balance in cells, affecting their homeostasis and their viability [59].

Using the above hybrid approach, a continuum model, based on a time-dependent reaction–transport equation (solved by the use of the finite-element method) was incorporated to simulate the time evolution of the growth factor concentration over the wound area, where, at the same time, the behaviour of the immune cells and fibroblasts was simulated on the basis of a discrete cell (particle) model. Specifically, Boon *et al.* [58] have evaluated the concentration field of the cytokine named tissue plasminogen activator (tPA), which breaks down the clot that initially forms to prevent blood loss through a wound puncture. The simulation clearly demonstrates an increase in the concentration of the tPA around the gap perimeter (figure 2*a*). Concurrently, the simulation also highlights advancement of the immune cells to the gap edge as well as differentiation of fibroblasts into myofibroblasts (figure 2*b*). Myofibroblasts are involved in the inflammatory response to injury, and are also able to apply stronger forces to their surrounding microenvironment; they are an intermediate between fibroblasts and smooth muscle cells. The migration of the myofibroblasts is modelled using ordinary, stochastic differential equations solved by the Euler–Maruyama method. As has recently been observed, both experimentally and in modelling, when fibroblasts migrate to close a gap they apply strong traction forces and generate new collagen networks, which accelerate gap closure [60]. To apply the strong traction forces fibroblasts also differentiate into myofibroblasts. Thus, in the model, Boon *et al.* [58] have added increasing numbers of myofibroblasts (by differentiation) which may apply traction forces. The pulling forces exerted by the myofibroblasts are modelled through inward (normally directed) forces applied by each cell on its boundaries; forces are applied using a Dirac delta function approach in a similar way as in the immerse boundary method. The forces are summed (using a superposition principle) and their effect on the entire domain of computation is determined using a finite-element method. This provides the displacement field, of interest specifically around the wound area. The paths of migration of the myofibroblasts and their traction forces appear to be in line with the collagen network (figure 2*b*). Thus, the wound edges in figure 2*b* become slightly curved compared with the initial wound, which had piecewise straight boundaries. Thus, the simulations using the hybrid model are able to show that concurrently at the wound site: tPA concentration increases, immune cells accumulate and fibroblasts differentiate into myofibroblasts that cause wound contraction; all of these are stages that occur simultaneously during wound healing.

As many models based on particle methods have been presented in the literature, we refer to the recent review on the use of particle methods for simulating wound healing by Vermolen [61]. The review includes many references of different kinds of particle (i.e. cell-based) methods (including hybrid methods, cellular automata methods and purely particle methods) for various subprocesses that may occur in wound healing.

Cell mechanisms such as differentiation (e.g. from fibroblasts to myofibroblasts [58]), division and death are commonly modelled as stochastic processes in the cell-based approaches. Recently, a cell-based model that models angiogenesis, a crucial stage during wound healing, was proposed by Bookholt *et al.* [62]. Figure 3 shows a snapshot of a simulation of the development of angiogenic sprouts (see the topography of the top coloured surface). An alternative hybrid model, for angiogenesis under the influence of mechanical signals, was developed by Stephanou *et al.* [63].

### 3.3. Models including cell deformation

Recently, models have been introduced that are based on an even smaller scale, where both the migration and deformation of the cells are taken into account; those are relevant to wound healing, immune cell response, as well as in tumour cell migration. Cell deformation models were developed based on PDEs, such as by Marth & Voigt [64], who use a phase-field approach to simulate cell motility under the influence of external signals, or by Madzvamuse & George [65], who use a moving mesh approach in a viscoelastic setting to model cell migration and deformation. Here, also the work by Barreto *et al.* [10,66,67] should be mentioned where a structural finite-element method is applied to model cell mechanics, and in particular with respect to environmental variations. These elegant approaches model cell deformation and migration very well, and are based on sound physical principles; however, the computational approach is relatively expensive. To this extent, Vermolen & Gefen [43] developed a three-dimensional, phenomenological cell deformation and migration model. That model is based on equations of motion for a set of points on the cell outer surface (its membrane) that are subject to elastic forces from the cell centre and their neighbouring points, and external (chemotactic) signals. The model also incorporates cell–cell interactions. For example, the model has been used to evaluate the deformation of an immune cell as it encounters pathogens that it aims to engulf (figure 4). The simulation shows that when the immune cells senses a factor secreted by the pathogens, it extends itself to engulf them. In the simulations, the chemical is modelled by solving a diffusion equation with three point sources whose locations coincide with the spatial positions of the centres of the pathogens. The solution is represented in terms of Green fundamental solutions.

Vermolen *et al.* [68] extended the deformation model to simulate the immune system response in tissue. This includes the deformation of white blood cells that migrate out of a blood vessel (extravasation), to pursue pathogens and as they chase the pathogens in the tissue. The simulations allow fast processing times and are based on linear combinations of Green fundamental solutions as well as on stochastic differential equations for the migration and deformations of the points on the boundaries of the immune cells. The flow of blood through the venules is modelled via Poiseuille flow. A snapshot of an *in silico* simulation is shown in figure 5. In another approach to reduce the model complexity, Borau *et al.* [69] formulated a three-dimensional probabilistic voxel-based model for the deformation and migration of a single cell. This model is phenomenological, yet still provides quite a lot of biological detail. It was used to investigate the impact of certain antibiotic treatments on the engulfment of all the pathogens.

## 4. Evaluation of the presented mathematical models

We presented several types of models that were either continuum-based (or macroscale formulated) in terms of PDEs or cell-based (or microscale) where individual cells are treated. In some of the models of the second class, the cell geometry is kept constant, whereas in others, the cell geometry and structure are allowed to change over time. Regarding the PDE-based models, where cell densities rather than individual cells are considered, in some cases, such as in Sherratt & Murray [23], Gaffney *et al.* [24] and Adam [15], (near) closed-form expressions could be obtained for the solution of the equations. This is only possible under severely limiting conditions, such as low dimensionality, homogeneity, elementary geometry (such circular or planar) and having a relatively simple PDE (system) to deal with. These (near) closed-form expressions provide very valuable information about the nature of the solutions and can predict trends of effects, which can then be used to test the reliability of numerical approximation techniques. Further, closed-form solutions can be used for regression procedures with experimental data. On the other hand, under more realistic circumstances, these solutions will provide and overly simplistic representation.

The contraction models and more complicated transport–reaction models represent a set of PDEs whose solution can only be approximated using discretization techniques such as the finite-element or finite-difference methods. Here, the parameters may be time- and space-dependent with possible nonlinearities to deal with. In particular, the contraction models, which also involve a first-order differentiation with respect to spatial coordinates owing to passive convection, pose additional numerical challenges if the diffusive (random walk) term is dominated by chemotactic or passive convection terms. In that case, the PDE becomes more hyperbolic. This necessitates the use of either upwind-like techniques, such as streamline upwind Petrov–Galerkin (SUPG) methods, or the use of limiters, in order to preserve positivity and monotonicity of the numerical approximations; in any case, it becomes attractive to use a high mesh resolution in regions that are characterized by sharp transitions. SUPG techniques are used for instance in the works of Javierre *et al.* [32] and Valero *et al.* [34], whereas slope and flux-limiting strategies have been used in Koppenol *et al.* [70]. Epshteyn [71] uses a discontinuous Galerkin (DG) method to solve the equations for chemotaxis. This could possibly be extended to passive convection as well, which has not been done yet as far as we know. The DG method, combined with filtering, as designed by Mirzaee *et al.* [72], seems to be a promising method for highly accurate approximations of solutions in the case of presence of chemotaxis, where also gradients of the solution are needed in structured and unstructured meshes. In particular, the modern filtering techniques to get high-order accuracy for spatial derivatives of the approximations, where classical methods reduce the accuracy, are very promising and provide a large extent of geometrical flexibility. A drawback of the DG method however, where the numerical approximation is discontinuous over the edges of the elements, is the large number of unknowns in the case of higher-order piecewise polynomial approximations. Furthermore, the bandwidth of the resulting discretization matrix is usually large, which imposes more requirements on computer memory usage and hence on computational resources.

The cell-based methods represent a large system of ordinary differential equations, possibly augmented with PDEs in the case of hybrid methods (for representing mechanical or chemical signals for instance). These models are in general computationally demanding and accordingly are typically feasible for the smaller scales. Monolayer cultures are perfect candidates for these models. Fortunately, the resulting equations for the migration and evolution of cell populations allow parallel implementation of the solution approximation method. Parallel computing or usage of a GPU computing environment is very helpful for these models if one wants to use large numbers of cells. An example of a GPU (CUDA) study on colony models is the work by Woods *et al.* [73]. This computational framework facilitates the treatment of large numbers of cells in a colony, so that possible clinical circumstances and larger computational areas and interactions can be dealt with. One of the advantages of the cell-based models is that these models are experiment-based, and that the influence of cellular processes, such as cell deformations, can be quantified. In this way, it is possible to evaluate, for instance, the influence of cell stiffness on the ability to penetrate into the surrounding tissues (as in cancer metastases) or on the engulfment of pathogens. Another advantage is the ability to visualize the computational results, so that physicians, nurses, other relevant medical professionals as well as students can actually see how the various biological processes (like wound closure or contraction) are happening, which indirectly contributes to the quality of healthcare. Another advantage of the use of cell-based models is the straightforward inclusion of (partly) stochastic biological processes owing to a lack of knowledge of the full time-history path of all the cells.

Where the continuum-scale PDEs only predict the mean of the possible outcomes, the stochastic nature of the cell-based models allows one to evaluate the spread of all possible outcomes. However, if one wishes to study wound healing, or any other biomedical process, in a larger framework, then the upscaled PDEs are more suitable. We also note that research on upscaling techniques from microscales to macroscales is highly relevant and is a promising direction to take. Steps in this direction were already taken indeed, mostly in the community of modelling flow and transport in porous media; see for instance, Bringedal *et al.* [74] and Van Noorden *et al.* [75]. Such upscaling techniques can provide a rigorous mathematical basis for the transition between micro- and macroscale models by the use of averaging techniques; they are also applicable to fields other than wound healing where cell migration is relevant, such as cancer, organogenesis and immune system response studies.

## 5. Summary

We have provided a review of the current approaches to mathematically model and simulate wound healing processes. We focused on three general types of models, including increasing levels of system detail. Specifically, we began with simplified models based on PDEs, which focus on the edge of a wounded cell monolayer considering its migration, restructuring and contraction. These models fall within the class of moving boundary problems. Then, we have discussed somewhat more complicated models in terms of reaction–transport equations, as well as mechanically oriented models being considered, where no interface was followed. Here, the edge of the wound coincides with a certain level curve or surface in three dimensions (which is a curve or surface where the solution is equal to a predefined value) of the solution of one of the PDEs. These models already contain some additional representation of the relevant biology, which makes them more complicated in the sense that they often contain various parameters that are difficult to measure or to obtain from regression techniques using experimental outcomes *per se*. We then progressed to cell-based models, specifically cellular automata and particle methods, which consider specific cells at different location near and away from the closing gap. Finally, we showed how to add elements in order to take into account cell deformation in addition to the migration. While there is much to learn from the bird's-eye simplified models, the more detailed models provide an improved understanding of the processes that each cell undergoes. A drawback of very detailed models is that they typically require many parameters that are generally difficult to estimate. Thus, to be able to forecast efficacy of specific physical conditions or treatments, it is likely that a combined approach, of *in vitro* experiments but also modelling on various levels is required.

## Competing interests

We declare we have no competing interests.

## Funding

The work was partially supported by the Israeli Ministry of Science and Technology Third Age Program (D.W. and A.G.).

## Acknowledgements

The authors thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme and conference on The New Mathematical and Computational Problems Involved in Cell Motility, Morphogenesis and Pattern Formation, supported by EPSRC grant number EP/K032208/1 (F.J.V. and A.G. were fellows in this programme).

## Footnotes

One contribution of 12 to a theme issue ‘Coupling geometric partial differential equations with physics for cell morphology, motility and pattern formation’.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.