Abstract
The first function of the skin is to serve as a protective barrier against the environment. Its loss of integrity as a result of injury or illness may lead to a major disability and the first goal of healing is wound closure involving many biological processes for repair and tissue regeneration. In vivo wound healing has four phases, one of them being the migration of the healthy epithelium surrounding the wound in the direction of the injury in order to cover it. Here, we present a theoretical model of the reepithelialization phase driven by chemotaxis for a circular wound. This model takes into account the diffusion of chemoattractants both in the wound and the neighbouring tissue, the uptake of these molecules by the surface receptors of epithelial cells, the migration of the neighbour epithelium, the tension and proliferation at the wound border. Using a simple Darcy's law for cell migration transforms our biological model into a freeboundary problem, which is analysed in the simplified circular geometry leading to explicit solutions for the closure and making stability analysis possible. It turns out that for realistic wound sizes of the order of centimetres and from experimental data, the reepithelialization is always an unstable process and the perfect circle cannot be observed, a result confirmed by fully nonlinear simulations and in agreement with experimental observations.
1. Introduction
Adult skin is made up of three layers: the epidermis and the dermis separated by the basement membrane. When a profound injury occurs which destroys a part of the dermis, it has to be mended rapidly to restore the protective barrier [1,2]. In vivo wound healing is a complex repairing process orchestrated by intra and intercellular pathways. To recover the integrity of the skin, this process is accomplished by four successive but overlapping stages: clotting, inflammation, reepithelialization and remodelling. Immediately after an injury, a clot composed mainly by fibrin fibres and platelets is formed to plug the wound (stage 1). Subsequently in the next 2–10 days (stage 2) the clot is continuously infiltrated by inflammatory cells which clear up the debris and release chemical factors, such as vascular endothelial growth factor and transforming growth factorbeta. Once the chemical gradient is established, different cells are recruited to fabricate the new tissue (stage 3). On the one hand, neovasculature [3] and fibroblasts (depositing collagen) are recruited in the dermis, reforming the clot into a granulation tissue, which nutritiously and physically support the repair of upper layers. On the other hand, keratinocytes migrate and proliferate at the edge of the wound to extend the newly formed epithelial carpet made of several layers of cells in the epidermis. This process is called reepithelialization and it lasts for two to three weeks. At the end of stage 3, myofibroblasts transformed from fibroblasts contract and try to bring the wound edge together. They disappear by apoptosis in the dermis. The remodelling (stage 4) continues for months or even for years in order to restore the homoeostasis of the normal skin. However, the normal anatomical structure is not truly recovered and a scar is formed from the granulation tissue. While this is true for large wounds for adult humans, those in human embryos [4] may be perfectly closed, the reasons being not fully understood. Interestingly, the scar formation is believed to be an evolutionary sacrifice in order to achieve rapid wound closure [5] and this is indicated by the intense and redundant inflammatory response [6] which mediates the systematic cell behaviours [7] during reepithelialization through releasing chemical signals. In all cases, the final scar is an outcome of the interplay between cells and the microenvironment [8] through physical and chemical factors during the four stages.
Here, we focus on reepithelialization when keratinocytes, the epidermal predominant cells, detach from the basal membrane, migrate towards the wound margin and proliferate, driven by morphogens in the granulation tissue [9]. Very careful in vitro epithelium migration experiments have been performed recently by several groups [10–13]. Realized on a solid substrate, they involve an advancing monolayer motivated by the collective behaviour of cells through migration, generation of forces and reorganization of the cell cytoskeleton at the margin. The findings indicate a nonintuitive manner of cells to achieve robustness during the repair where chemotaxis seems not to be decisive for the global cell migration in vitro. Nevertheless, owing to the abundant morphogenic signals released in the wound in vivo compared with in vitro experiments, chemotaxis dominates the collective behaviour of cells during the woundhealing process. Closer to in vivo wound healing (figure 1a(i)(ii)), an artificial cornea has been reconstituted [15] (figure 1b(i)(ii)) by the introduction of epithelial cells, fibroblasts and growth factors in a threedimensional Petri dish. Interestingly, all these works [11,12,15] concern the circular geometry rather than the linear cut done in earlier experiments [16,17]. Although well controlled, it is unclear that these systems present the whole complexity of in vivo wound healing [14]. Nevertheless, both in vivo and in vitro experiments present disordered border evolution, which indicates a possible universal irregularity owing to chemotaxis (figure 1).
In this work, we describe the reepithelialization by considering a circular hole in a tissue, adapting a recent dynamical model of chemotactic migration [18] driven by morphogenetic gradients. Being continuous, the model describes the tissue at a scale larger than the cell size but smaller than the size of the hole. Here, we do not distinguish cell species or the chemoattractant categories. The latter is continuously supplied by an incoming flux from the third dimension normal to the epithelial layer and uptaken by cell receptors. As shown in [18], the thinness of the moving layer transforms the threedimensional process into a twodimensional model, where the localized thickness variation at the border contributes to a tension T. In addition, we assume that a strong viscous friction exists between the moving layer and the substrate [19–23], which allows writing a Darcy's law for the average velocity field inside the tissue, while cell–cell interactions and mitosis are transformed into interface boundary conditions in the sharp interface limit. We first present the physical model, then the analytical treatment giving an explicit solution for the circular geometry and a study of its stability. Indeed in the case of a linear border, it has been shown that a long wavelength instability occurs at low velocities owing to a Goldstone mode (translational invariance) that surface tension alone cannot succeed to prohibit [18]. Intuitively, one may expect that it disappears for small holes as the surface tension effect accentuates. On the contrary, as for experiments in vitro, larger holes may exhibit contour instabilities leading to dynamical behaviours that we aim to predict. To go beyond the stability analysis, the timedependent freeboundary problem requires numerical methods able to solve equations on evolving domains. We choose to tackle our problem by levelset methods, first developed in the 1980s by Osher & Sethian [24], which are able to track moving interfaces and topological changes automatically and have been successfully applied to problems, such as liquid–gas interactions [25], image processing [26] and tumour growth [27–30]. Owing to the nature of our theoretical problem, we select the methodology developed for tumour growth [27–30]. In the following, we first present the model, introducing the most relevant physical parameters, then the analysis for weak amplitude instability, and finally numerical methods and simulations for the ultimate state of the closure.
2. The model
We consider a circular hole Ω in a twodimensional cellular population of density ρ immersed in a morphogenetic environment (figure 2). Cells migrate up the morphogen gradient, the chemotactic flux being proportional to the concentration gradient of morphogens , with Λ_{c} being a mobility constant. The morphogen repartition inside and outside the hole Ω is alimented through sources coming from below. Mitosis happens only at the border [15,31] and the closure is mostly achieved by the migration. The mass balance for the cell population is then 2.1
Neglecting volumetric growth (γ = 0) in the tissue except at the periphery, the cell density is constant and ρ = ρ_{0}. Equation (2.1) simplifies to and gives that the normal velocity of the front is directly proportional to the normal concentration gradient. At extremely low velocities, the cell migration satisfies a Darcy's law , M_{p} being a porosity coefficient equal to the square of the epithelium height divided by a friction coefficient. As shown in [18], this law is deduced in tissues when friction between phases or friction with a substrate balances the hydrostatic part of the elastic stress acting on the cells [18,20]. The wound perturbs the homoeostatic state of the surrounding and a source of morphogens will try to restore an optimal equilibrium value in the aperture c_{0}. Taking this concentration as unit for the morphogen concentration (giving ), τ_{c} the uptake time as time units, L_{e} as length unit (with , D_{e} being the diffusion coefficient inside the epithelium), we get 2.2with the index h or e referring to the inner (hole) or outer (epithelium) domain. δ = D_{h}/D_{e} represents the ratio between the diffusion coefficient in the hole D_{h} and the tissue D_{e}, it is bigger than 1, and α gives the strength of the transverse flux which maintains the morphogen level inside the hole. Owing to the relative slowness of the cell migration, we neglect the time dependence in equation (2.2). Taking D_{e}/M_{p} as pressure unit simplifies Darcy's law into , where for simplicity we keep the same notation for p and v. From mass balance equation (2.1), we get the following Laplace equation coupling the unknown pressure p to the chemoattractant concentration c: 2.3So, we get , where ϕ is a holomorphic function which satisfies Δϕ = 0. Finally, we must pay special attention to the interface boundary conditions. For equation (2.2), they concern the concentration continuity and the flux discontinuity owing to the morphogen consumption by mitosis at the border (for demonstration, see [18]): 2.4where N is the outer normal. Γ_{2} is the uptake rate constrained at the border discussed below together with the mitosis rates Γ_{1} and at the border. Capillarity fixes the pressure jump at the interface 2.5which gives the freeboundary condition considering the geometric effect where is the local radius of curvature. is equal to the radius if the geometry is a circle. However, it considers the local effect at small length scale when the boundary is perturbed. At the same time, the velocity of the interface differs from that of the epithelium v by the cell proliferation at the border and we get 2.6where Γ_{1} is the mitosis rate and σ the capillary number related to the tension T at the interface so Wound healing in vivo is dominated by twodimensional migration and the morphogens are probably not the nutrients, so Γ_{2} vanishes and Γ_{1}c must be replaced by if we also consider mitosis without the quantitative effect of the nutrient concentrations. The capillary number σ may involve the activities of actin cable (bundled microfilaments) in the cells. During embryonic wound healing, the leading edge cells can coordinate their actin cables globally which generate an effect on the macroscopic level. However, this coordinated behaviour is lost in adults, so the effect of actin cables may act only locally. This resembles σ in our model which is not effective for large wounds. Considering the set of equations (2.2) and (2.3) in conjunction with the set of boundary conditions (2.4), (2.5) and (2.6) shows that chemotaticdriven migration is indeed a freeboundary problem [32,33] involving several parameters to represent the biological complexity and the study of simple cases, for example the circular closure, may help in understanding their corresponding roles. Thus, we present first the analytical results for a hole remaining circular at all times, then its stability.
3. The results
3.1. Regular circular closure
In the quasistatic approximation, equation (2.2) can be solved analytically and gives 3.1
I_{0} (resp. K_{0}) being the modified Bessel function of zero order, regular at r = 0 (resp. r → ∞). Equation (3.1) takes into account the continuity at the interface, A being given by the flux continuity and reads 3.2with the following definition for (ratio of two successive Bessel functions) and the equivalent for The pressure P_{0} inside the epithelium becomes 3.3We use the Laplace law to fix the unknown degree of freedom. The holomorphic function ϕ, proportional to log(r), represents a possible driving force which appears at the interface (i.e. the socalled kenotaxis defined in [13]). Indeed, even in the absence of morphogens, in vitro epithelium migration on solid substrate is observed, perhaps owing to the reorganization of the cellular cytoskeleton at the wound margin. Note that the model remains valid if we discard the chemotaxis provided the definition of length and time units is modified. Furthermore, C_{0} and P_{0} are dependent on time via R(t). For simplicity, we drop the time dependence of The velocity of the closure deduced from equation (3.5) is then 3.4
The closing velocity is constant for large holes and is exponentially small when the radius R becomes tiny, if one restricts on chemotatic migration and proliferation. So the radius of a large hole begins with a linear decrease in time but the total closure will take an infinite time for a complete achievement. For small holes, chemotaxis becomes subdominant compared with kenotaxis A_{in} at the interface which controls the closure dynamics and the radius satisfies a diffusive law in t^{1/2} [12]. Figure 3 shows as a function of R for different values of α, δ and A_{in}/Λ and §4 contains a discussion of parameters.
3.2. Loss of circularity
However, equation (3.4) is valid only if the circular contour is maintained. It is why we perform a linear stability analysis for large wounds (A_{in} ∼ 0) assuming a small harmonic perturbation for the radius as inducing variations on the pressure P and concentration field C of the same order ε as follows: 3.5where the subscript i indicates either a quantity relative to the hole (h, 0 < r < R) or to the epidermis (e, R < r < ∞). Although all perturbative quantities depend on the selected mode n, the linear perturbation analysis treats these modes independently. So we drop the n index and calculate the perturbative concentration fields c_{i}(r) from equation (2.3) 3.6while the perturbed pressure field is given by where we take into account the harmonic modes of the holomorphic function ϕ with B fixed by the Laplace law Owing to the weakness of ε, our system of equations, once linearized, can be solved analytically and the growth rate of the mode n reads 3.7where Equation (3.7) represents an implicit relationship for Ω_{n}, solved by iterative techniques and positive values indicate the modes responsible for the destabilization of the circular border as the migration proceeds. The results are presented for different values of Λ, α and δ, where Ω_{n} is displayed in figure 4. See §4 for the discussion of parameters. The results indicate an instability leading to a deviation from a circle in short time, for n up to a critical mode n_{c} (fixed by the capillarity). It is why numerical simulations are necessary to go beyond the linear analysis and fully consider the nonlinearities.
3.3. Full dynamics and numerical methods
We discretize c in equation (2.2) and pressure p in equation (2.3) on a Cartesian mesh in space and implicitly in time, using a nonlinear adaptive Gauss–Seidel iterative method [29,30]. At the domain boundary, we impose both vanishing values of c and normal gradient of p. Then the noise is added on C = c + χ with Beginning with a perfect circle with R = 90, the wound closing is tracked by the levelset method developed in [27,28] where a scalar function Φ with describes the wound (Φ < 0), the epithelium (Φ > 0) and the interface (Φ = 0). The normal and curvature in equation (2.5) are calculated by standard differential geometry: and Φ is updated by the relation , where V_{ext} is the constant extension of V_{int} on the interface from equation (2.5) following the gradient Note that this method simulates a perfect circular closing without noise. According to our results, the wound deviates from a perfect circle soon after the growth as shown in figures 5 and 6 and in agreement with the stability analysis.
4. Discussion and conclusion
There are several independent parameters in this model. We can fix some of them with the published experimental data (table 1). Our velocity unit is compatible with the closing velocity of the cornea (3 days for a hole of radius approximately 3 mm [15]) giving Λ ∼ 1 for this experiment according to equation (3.2). In figure 4, we vary the parameters α and δ, which are more difficult to estimate, and also Λ. Owing to the large size of wounds in practice, the linear stability analysis gives always an instability and this conclusion is robust to parameter changes. This finding is confirmed by fully nonlinear simulations under the same parameter range (figures 5 and 6). In the simulations, large Λ (approx. 10, total time approx. 60) contributes to faster closing compared with small Λ (approx. 1, total time approx. 600). A typical example of the closing process with α = 1 and δ = 2 is shown in figure 5a. The advancing interface becomes wavy as the wound heals; however, when the wound becomes small, surface tension restabilizes the wound to a potato shape consistent with the linear stability analysis. This surface tension may also pinch wounds off to smaller pieces as shown in figure 5b(ii), figure 6 and in [12]. This event requires at first a more irregular closing (red arrows in figure 6), when the concentration of chemoattractant is less homogeneous given inadequate transverse flux (α = 0.1). The pinchoff can happen both at the intermediate or the end stage of the closing (yellow arrows), consistent with observations at the late stage of wound healing on a monolayer of MDCK cells (figure 6, right), suggesting the same behaviour in more complex woundhealing processes.
In this work, we have shown theoretically that wound reepithelialization gives a border instability under clinically realistic parameters. Driven by chemotaxis, our model does not introduce other cell activities which close the wound at the ultimate stage. However, we conjecture that this chemotactic instability at the border of the wound may affect the quality of the final repair. During embryonic wound healing where perfect reconstruction is observed, the wound is closed by a ‘purse string’ which represents a coordination of leading edge cells facilitated by the actin cable [35]. This mechanism is lost in the adult skin, where the wound is closed by the crawling of cells (several layers) on the granulation tissue [36]. Compared with a static substrate in vitro, the granulation tissue undergoes contraction at the end of the reepithelialization by myofibroblasts [37]. Indeed, this purpose is to bring the edge together which resembles the ‘purse string’, but the contraction by those cells needs to be compatible with the synchronous material constitution, which is mechanically and systematically a challenge [7].
At the end, we discuss our model in the context of the skin wound healing in vivo. As the wound goes deeply into the dermis which is the case for most wounds, the two layers behave differently during reepithelialization. The depth of a wound contributes to the thickness in the lower layer where the chemoattractant is released. Considering the reepithelialization where only several layers of cells become motile above the lower stacks where the chemoattractant is released, the depth of the wound contributes mathematically to the transverse chemoattractant flux strength α. Besides α, which may involve the depth of the wound, the geometry, the size of the wound as well as the diffusion coefficient ratio δ are considered. Comparing with experiments in vivo, the human woundhealing process may be different from that in experimental mice owing to different skin biomechanical properties. However, if one intends to consider the problem in vivo accounting for the bioelasticity, the challenge does not only come from different components of the two layers, but also the fact that the junction between the dermal and epidermal layer is highly disordered even for the intact skin [38], as shown in our previous work. A regular pattern should be defined in the third dimension which we consider as ‘normal’. Also, as the wound tissue undergoes remodelling during the healing, the mechanical parameters cannot be considered as constants and the relaxation of the tissue probably evolves on a larger time scale than the reepithelialization. Indeed, these relaxation processes should be also considered in predicting the quality of the final repair. Before that, the irregularity of the wound border driven by chemotaxis should be considered as an ingredient.
Funding statement
This work was supported in part by AAP Physique Cancer 2012.
Acknowledgements
We acknowledge Pascal Silberzan and Olivier CochetEscartin for discussion and for providing photographs from their experiments. We also thank Patrick Carrier for discussion.
Appendix A. Numerical implementation
The simulation is done on a 200 × 200 Cartesian lattice with a grid size of 10 μm. This freeboundary problem is easy to be implemented with the levelset method owing to its dependence on the curvature in solving the pressure in equation (2.5). The full method is adapted from [28,29] and has an order of accuracy of more than 1.5. While the boundary can be implicitly given by the zero contour of the levelset function Φ, the curvature can also be readily calculated. At each time step, the domain from Ω_{t} to Ω_{t}_{+Δt} is updated as follows:

— Solve equation (2.2) with the fixed domain Ω_{t} in the steady state with the freeboundary condition equation (2.4). At the boundary of the computational lattice, the Neumann boundary condition is imposed.

— Update the solution of concentration C by C = c + χ pointwisely, where χ is randomly generated out of uniform distribution between and is used for the presented simulations.

— Solve equation (2.3) with known C under the interfacial boundary condition from the first part of equation (2.5), where the local curvature is calculated from the levelset function by The Neumann boundary condition is imposed at the boundary of the computational lattice.

— Calculate the velocity V of the moving boundary from the second part of equation (2.4), where is calculated from the levelset function by

— Find the appropriate Δt given by the CFL condition by and update the domain Ω_{t} to Ω_{t} _{+} _{Δt}.

— Start with the newly updated domain Ω_{t}_{+Δt} and repeat the last five steps.
 Received November 8, 2013.
 Accepted January 2, 2014.
© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.