The crosstalk between fibroblasts and keratinocytes is a vital component of the wound healing process, and involves the activity of a number of growth factors and cytokines. In this work, we develop a mathematical model of this crosstalk in order to elucidate the effects of these interactions on the regeneration of collagen in a wound that heals by second intention. We consider the role of four components that strongly affect this process: transforming growth factor-β, platelet-derived growth factor, interleukin-1 and keratinocyte growth factor. The impact of this network of interactions on the degradation of an initial fibrin clot, as well as its subsequent replacement by a matrix that is mainly composed of collagen, is described through an eight-component system of nonlinear partial differential equations. Numerical results, obtained in a two-dimensional domain, highlight key aspects of this multifarious process, such as re-epithelialization. The model is shown to reproduce many of the important features of normal wound healing. In addition, we use the model to simulate the treatment of two pathological cases: chronic hypoxia, which can lead to chronic wounds; and prolonged inflammation, which has been shown to lead to hypertrophic scarring. We find that our model predictions are qualitatively in agreement with previously reported observations and provide an alternative pathway for gaining insight into this complex biological process.
Wound healing is a complex, finely balanced process and, in the case of successful healing, is generally thought to comprise four interconnected and overlapping stages: haemostasis (clot formation), inflammation, proliferation and remodelling . The formation of an initial fibrin clot, as well as its degradation and subsequent replacement with a stronger, mature matrix that is mainly comprised of collagen, occurs as a result of a series of interactions, involving the activity of a number of cells, growth factors (GFs) and cytokines [2,3]. Two key cells involved in this process are fibroblasts and keratinocytes , whose interaction, via a paracrine loop , is essential to the outcome of successful dermal remodelling  and re-epithelialization . This dermal–epidermal crosstalk is currently an active area of biological research (see Werner et al.  for a comprehensive review). It is known that keratinocytes and fibroblasts cooperatively induce gel contraction in vitro . More recently, it was observed that keratinocytes can regulate the modulation of fibroblasts into a phenotype that is more suited to the remodelling, rather than the deposition, of the extracellular matrix (ECM) [6,10], thus precipitating the subsequent stage of wound healing. Disruptions to this process can lead to serious pathologies such as hypertrophic scars , keloid scars  and non-healing wounds . It is therefore of crucial importance to throw light on such networks. As it is difficult to conduct in vivo investigations into these processes in a non-invasive manner, realistic mathematical models that are based on known cell behaviours provide a useful framework for studying the wound healing process and its dysregulation.
One of the earliest mathematical models of wound healing was developed by Sherratt & Murray , who used two coupled reaction–diffusion equations to describe re-epithelialization, as well as the dependence of the rate of wound closure on the rate of cell proliferation and the migratory activity of the cells. Tranquillo & Murray  developed a model for the closure of a dermal wound that included the production of collagen by the fibroblasts and the effects of cell traction on the rate of wound closure. This model has provided the basis for most subsequent mechanochemical models of dermal wound healing, including those by Olsen et al. , Tranqui & Tracqui  and Vermolen & Javierre . Cruywagen & Murray  developed a continuum mechanical model of the interactions between the dermal–epidermal layers of a developing embryo, assuming that ‘signal morphogens’ are emitted in one layer and have their paracrine effect in the other. Several models have also been developed to study the actions of GFs during the wound healing process, including the effect of keratinocyte growth factor (KGF) signalling on dermal–epidermal interactions  and the chemotactic and mitotic effect of platelet-derived growth factor (PDGF) on fibroblasts in dermal wound healing , while Murphy et al. [22–24] recently investigated the effects of transforming growth factor-β (TGF-β) and tension on dermal wound closure. There have also been several models that describe the laying down of collagen fibres in a wound, including those by Dale et al.  and Dallon et al. , who investigated the effect of TGF-β on this process, and more recently by Cumming et al. , who used a hybrid multiscale discrete-continuum model. Despite these numerous advances, to our knowledge, there are no mathematical models of wound healing that consider the fibroblast–keratinocyte crosstalk that we argue is central to the wound healing process.
In this work, we simulate the behaviour of a wound that is not surgically closed and which therefore heals by second intention (figure 1). In §2, we develop a mathematical model that describes a set of interactions between fibroblasts and keratinocytes that are vital to the proliferative phase of healing, as well as several critical GFs and cytokines that mediate this crosstalk. As illustrated in figure 1, during this stage, enzymes produced by fibroblast cells degrade the initial fibrin clot and synthesize a collagen-based ECM in its place , while keratinocyte cells migrate across the wound surface leading to re-epithelialization . In §4, we show that our model captures essential features of the normal healing process, and can be used to describe abnormal wound healing pathologies. It is widely accepted that persistent, chronic hypoxia inhibits healing [30,31] and is one of the most common causes of chronic wounds  in which healing fails to produce anatomical and functional integrity . In addition to chronic hypoxia, we consider the case of an excessively long inflammatory response, following Murphy et al. . We find that when treatment is initiated within an appropriate time frame, our model predicts that a normal healing outcome is achievable in each pathological case.
2. Mathematical model
We model a wound of partial thickness in which both the dermal and epidermal layers are damaged. We present our model in terms of spatial variables (x, y) and use this model to simulate wound healing in a two-dimensional domain. In particular, we focus on the complex network of interactions that comprise the crosstalk between keratinocyte and fibroblast cells, whose densities we represent by k(x, y, t) and f(x, y, t), respectively. Following the approach of Cumming et al. , we assume that a fibrin clot containing platelets is formed during the initial haemostasis phase and, on being degraded by enzymes produced by fibroblasts, is replaced with a matrix mainly composed of collagen that allows keratinocyte migration. In order to describe this process, the matrix in our model is represented by two species, c(x, y, t) and ρ(x, y, t), corresponding to the fibrin clot and ECM densities, respectively. Wound healing involves the activity of a large number of GFs and cytokines (see  for a recent review). We restrict our attention to the activity of those species that strongly impact the proliferation and migration of fibroblasts and keratinocytes, as well as the degradation and renewal of the ECM. Specifically, we consider the concentrations of interleukin-1, i(x, y, t); KGF, G(x, y, t); PDGF, P(x, y, t); and TGF-β, β(x, y, t).
In figure 2, we present a schematic of what we believe to be the most important interactions between these eight species during the dermal wound healing process. In this work, we develop a two-dimensional model of nonlinear reaction–diffusion partial differential equations to describe these interactions. We now discuss the biological phenomena that underpin the form of the equations.
Re-epithelialization is driven by the proliferation and migration of keratinocytes. It has been observed that this migration is dependent on the presence of an appropriate proximate substrate, such as the newly formed collagen-based matrix [35,36], and is stimulated by both the presence of mature collagen  and fibronectin (which is present in the fibrin clot) . We note that in several models of angiogenesis (e.g. ) it is assumed that endothelial cells require an ECM in order to migrate into the wound bed. In our model, we assume that keratinocyte migration is dependent on the densities of the fibrin clot and mature collagen. Owing to the rapid degradation of fibrin, this migration effectively occurs along the fibrin–collagen interface. This is consistent with the observation that the keratinocytes form an invasive ‘tongue’ prior to re-epithelialization, as shown in figure 1. The proliferation of keratinocytes is upregulated by the presence of fibrin  as well as mature collagen , and we model this using a logistic form with an intrinsic proliferation rate constant λ1. This proliferation is upregulated by the presence of KGF [42,43] and is downregulated by the presence of TGF-β , which we model using the coefficients λ2 and λ3, respectively. Assuming that this species has a maximum density and random motility coefficient Dk, the fibrin clot has initial density , and the density of mature collagen in the healed state is , we thus have: 2.1
The proliferation and migration of fibroblasts are critically important to successful dermal wound healing, as they replace lost collagen . It has been observed that there is an increased migration of fibroblasts owing to a chemotactic response to both PDGF  and TGF-β , and we model this using terms with chemotactic coefficients χfP and χfβ, respectively. We model the proliferation of fibroblasts using a logistic form with an intrinsic proliferation rate constant λ4. Keratinocytes are observed to have a complex effect on the behaviour of fibroblasts: while they secrete interleukin-1, which upregulates fibroblast proliferation , they also physically control the growth of fibroblasts by forming colonies . We model these processes using the coefficients λ5 and λ7, respectively. Fibroblast proliferation is further upregulated by the presence of TGF-β , which we model using the coefficient λ6. Assuming that this species has a maximum density and random motility coefficient Df, we obtain: 2.2
(c) Initial provisional matrix (fibrin clot)
The initial provisional matrix is degraded by enzymes produced by a variety of cells, particularly fibroblasts [49,50]. As an approximation to this more complicated process, we model the degradation of the provisional matrix as being (linearly) dependent on the fibroblast cell density at a rate characterized by the constant λ8. We thus have: 2.3
(d) Extracellular matrix
Fibroblasts produce collagen, the main component of the new ECM . We model this dependence in line with Flegg et al.  and Tracqui et al. , where a kinetic term proportional to the fibroblast population is used. This production is dependent on the concentration of oxygen ; however, to avoid further complexity, we assume the dependence of collagen production on oxygen is incorporated into the coefficient λ9. Interleukin-1 is known to downregulate the production of collagen by fibroblasts  and TGF-β is known to upregulate this process . We model these processes using the coefficients λ10 and λ11, respectively. As fibroblasts are known to degrade newly synthesized collagen , following Dallon and co-workers [26,58] we assume that this process occurs at a rate, characterized by the constant λ12, proportional to the collagen density. We thus have: 2.4
An important cytokine in the wound healing process is interleukin-1 , whose production by keratinocytes  we model using a linear term with rate constant λ13. Cytokines and GFs tend to have a half-life of the order of a few minutes , and so this decay is assumed to decay linearly with rate constant λ14. Assuming this species diffuses at a rate characterized by the constant Di, we thus have: 2.5
(f) Keratinocyte growth factor
KGF is known to play an important role in epidermal wound healing . Its production by fibroblasts, which is upregulated by the presence of interleukin-1 [42,43], is modelled using a nonlinear term with coefficients λ15 and λ16. The decay of this species is assumed to be linear with rate constant λ17. Assuming this species diffuses at a rate characterized by the constant DG, we thus have: 2.6
(g) Platelet-derived growth factor
PDGF is produced by degranulating platelets present in the initial fibrin clot . The rate of production of PDGF decreases as the fibrin clot is degraded, and we model this production using a term with rate constant λ18. This GF decays naturally at a rate characterized by the constant λ19, and is further removed by fibroblasts through endocytosis  at a rate characterized by the constant λ20. Assuming this species diffuses at a rate characterized by DP, we thus have: 2.7
(h) Transforming growth factor-β
The protein TGF-β exists in a number of isoforms, each of which has different biochemical activities . In our model, we consider the activity of the TGF-β1, as it is the most prevalent isoform in wounds  and its activity in the formation of scars is known to be significant . TGF-β is secreted in a latent form and is then converted to an active form by enzymes  or by myofibroblasts . For this investigation, we assume that this conversion occurs rapidly, and we represent both active and latent forms of TGF-β by a single species, β. Fibroblasts are one of the main sources of TGF-β in wounds , and their production of this GF is inhibited if TGF-β is already present . We model this production using a saturation form, in line with the work of Murphy et al. , with the constant λ21 characterizing the rate of production by fibroblasts and the constant λ23 characterizing the half-maximal rate of production. It has been observed that this production by fibroblasts is downregulated by the presence of keratinocytes , and we model this using a downregulatory coefficient λ22. Modelling the decay of TGF-β using a linear term with rate constant λ24 and assuming that it diffuses at a rate characterized by the constant Dβ, we thus have: 2.8
(i) Geometry, initial and boundary conditions
We assume that (x, y) represents a plane perpendicular to the surface of the skin, which includes both epidermal and dermal regions. Note that in figure 1, the x-axis runs parallel to the surface of the skin and the y-axis is perpendicular to the surface. As illustrated in the schematic in the electronic supplementary material, figure S1, we model a wound of depth Ly and length 2Lx, with symmetry around x = 0, and where y = 0 represents the lower edge of the wound. We assume that x = Lx represents the interface where the wound meets healthy tissue and that represents the epidermal region. For our chosen values of Lx, Ly and (table 1), this assumption corresponds to the physical extent of the human epidermis. Our computational domain also contains some of the surrounding unwounded tissue and extends to the region , .
We assume that keratinocytes are initially present in the epidermal region outside the wounded area, , , while fibroblasts and the ECM are present in the region y < 0, as well as the dermal region outside the wounded area, , . In our model, interleukin-1 and KGF are produced by keratinocytes and fibroblasts, respectively, and so we assume that they are initially absent in the wounded area and present outside. As PDGF and TGF-β are produced by macrophages and platelets present in the fibrin clot , we assume that they are initially present in the wounded area and absent outside. For the GFs we use no-flux boundary conditions along each boundary, while for the keratinocytes and fibroblasts we use no-flux boundary conditions along x = 0 and y = Ly, and Dirichlet boundary conditions along and . These initial and boundary conditions are simulated using steep tanh functions and are detailed in the electronic supplementary material. We now present a set of numerical results obtained using this model that describe different wound healing outcomes. We include three movies in the electronic supplementary material that show the evolving densities of keratinocytes over the entire domain for each outcome, illustrating the development of the invasive ‘tongue’ and the subsequent migration of these cells across the wound.
3. Material and methods
The model equations were non-dimensionalized (this is detailed in the electronic supplementary material) and numerical results were obtained by integrating the resulting system using the NAG subroutine D03RAF  for different choices of the dimensionless parameters over a two-dimensional domain, with 121 × 121 grid points. A schematic of our non-dimensionalized model is shown in figure S2 of the electronic supplementary material. The dimensional and dimensionless values for the parameters of our model are listed in table 1. While our model is intended to be a caricature of the network of interactions that comprise the fibroblast–keratinocyte crosstalk, we have as far as possible used parameter values that correspond to reported experimental data. We now discuss how these parameters were determined.
(a) Scaling, diffusion and random motility coefficients
The carrying capacity of fibroblast cells is roughly , while the volume of each cell, which weighs the same as water, is around . We thus take . Noting that the molecular weight of KGF is around 27 kDa  and taking the value for the KGF concentration from Wearing & Sherratt , we have . Assuming that the density of PDGF is the same as that of KGF and noting that the molecular weight of PDGF is around 24 kDa , we take . The molecular weight of interleukin-1 was found to be in the range 12–70 kDa , and for our simulations we consider the higher value 70 kDa. Assuming that the density of interleukin-1 is also the same as that of KGF, we take . We assume that the densities of the ECM and the fibrin clot are comparable and take both and to be , which lies within the range of values previously used in simulations by Namy et al. . We also assume that the density of the epidermal layer is similar to that of the dermal layer, and take to be . Using our assumed value for the molecular weight of interleukin-1, we estimate the random motility coefficient by considering the molecular weight, m, and random motility coefficient, D, for another chemical using . Following Wearing & Sherratt , we consider the species N-formyl-methionyl-leucyl-phenylalanine, for which and . We therefore obtain the estimate . Finally, we assume that the dimensional values of the chemotactic coefficients and are g cm−1 s−1 and , respectively.
(b) Reaction kinetic coefficients
For our simulations, we assume that the intrinsic growth rate constants λ1 and λ4 have the dimensionless values λ1 = 14 and λ4 = 3.5, such that the wound closes within three weeks. It is known that KGF stimulates a 4.2-fold rise in keratinocyte proliferation ; the proliferation of keratinocytes will halve under the influence of TGF-β ; interleukin-1 stimulates a 2.5-fold rise in fibroblast proliferation , downregulates the production of collagen by fibroblasts by 50 per cent  and increases the production of KGF by fibroblasts by a factor of 6 ; and the presence of keratinocytes decreases the production of TGF-β by a factor of 1.6 . We thus set the dimensionless values λ2 = 3.2, λ3 = 1, λ5 = 1.5, λ11 = 1, λ16 = 5 and λ22 = 0.6. We set the dimensionless value λ7 = 10 to ensure that keratinocytes sufficiently restrict the proliferation of fibroblasts. The half-life of interleukin-1 mRNA was found to be in the range 60–90 min . For our simulations, we use the upper limit and take the dimensional value . The dimensionless values of and are chosen such that we obtain the steady state inside the wound, while the dimensionless value of λ18 = 2.4 is chosen such that there is a balance in the production and decay of PDGF in the absence of fibroblasts. For the degradation of PDGF by enzymes produced by fibroblasts, we use the upper limit of the range of values estimated by Haugh  and take the dimensional value . Finally, in the model by Murphy et al. , the rate of TGF-β synthesis by fibroblasts is estimated to be ng (cell d−1). This value was deduced from experiments by Wang et al. , where the maximum dilution of TGF-β was . Thus, we take the rate constant of TGF-β production for a given density, , of fibroblasts to be . The values of the parameters λ6, λ8, λ9, λ10, λ12, λ17, λ19, λ23 and λ24 were taken directly from the corresponding references listed in table 1.
(a) Normal wound healing
We simulate the normal healing case using the parameter values listed in table 1. Figure 3a shows how keratinocytes migrate into the wound bed at the interface between the intact ECM and the fibrin clot, reproducing the invasive ‘tongue’ seen in figure 1. These cells subsequently migrate across the wound site, leading to re-epithelialization, as seen in figure 3b,c. Although not shown here, we also find that fibroblasts migrate over the entire domain, completely degrading the fibrin clot and synthesizing mature collagen in the process. The large initial concentration of TGF-β in the wound is found to result in the production of an excessive amount of collagen by fibroblasts; however, the rapid drop in TGF-β concentration, coupled with increasing levels of interleukin-1, quickly downregulates this production, and the ECM density approaches its pre-wounded value. We now consider two pathologies that can arise owing to a disruption to the wound healing process.
(b) Prolonged inflammation
First, we consider the case of an excessively long inflammatory response. Murphy et al.  noted that poorly regulated inflammatory mediators can lead to pathological scarring and simulated this case by decreasing the decay rate of TGF-β. Such an assumption is justified by the observation that prolonged inflammation results in reduced levels of neutrophils, which in turn negatively affects the ability of chemoattractants to dissipate . Thus, we follow the approach used by Murphy et al.  and model this case by decreasing the value of the parameter that characterizes the decay of TGF-β, λ24, by a factor of 10. We find that although keratinocytes quickly migrate across the domain to re-epithelialize the wound (figure 4a), the ECM levels in this case are significantly increased.
The average ECM density (relative to normal) is calculated by averaging the net ECM density in the wound bed at each timestep over the total number of grid points in this region. A smooth function is then obtained using spline interpolation, and this is then divided by the corresponding function obtained for the normal case. Figure 5 shows that the ECM levels in the prolonged inflammation case, shown as a thick red curve, are consistently much higher than normal. This arises owing to the prolonged presence of TGF-β in the wound (prolonged inflammation), which stimulates fibroblasts to produce an excessive amount of collagen and often leads to hypertrophic scarring.
(c) Chronic hypoxia
We consider the case of chronic wound hypoxia, where a significantly reduced concentration of oxygen inhibits the healing process . Prolonged hypoxia is one of the most common causes of chronic wounds , where healing fails to produce anatomical and functional integrity . Siddiqui et al.  noted that fibroblast cells placed in chronically low oxygen levels over 72 h exhibited a 50 per cent decrease in population doublings, a 1.48-fold reduction in collagen production and a 3.1-fold decrease in TGF-β production relative to a standard oxygen system. These observations confirm the findings that chronic hypoxia reduces collagen production and that fibroblast proliferation is greatly inhibited in this case . We investigate this situation by changing the parameters λ4, λ9 and λ21 by these experimentally determined amounts. As seen in figure 4b, keratinocytes take much longer to migrate across the wound in this case.
The rate of wound closure is estimated by considering the extent of the non-epithelialized region in each case (see the electronic supplementary material for further details). Figure 6 shows the wound takes much longer to heal than in the normal case, which is in agreement with the conclusions of Oberringer et al. . Additionally, as seen in figure 5, our model predicts that collagen is produced at a slower rate than in the normal case.
(d) Implementation of treatment
Finally, we simulate the application of a treatment for the two pathologies under consideration. Prolonged inflammation can be treated by blocking the TGF-β pathway, for instance, by suppressing the secretion of TGF-β by fibroblasts (see Murphy et al.  for a recent discussion). We therefore model the treatment of prolonged inflammation by returning the rate of TGF-β decay to the corresponding value for the normal healing case. We note that early treatment can significantly impact the wound healing outcome. From figure 5, we see that a treatment applied at early times can cause the ECM density to quickly approach the corresponding value for the normal case. Furthermore, figure 6 shows that a treatment applied at days 2.5 and 5 can cause the wound to re-epithelialize at a rate closer to the normal rate, while treatments applied after 12.5 days will not affect this process. These findings are in agreement with the results of Murphy et al. , who noted that a healing strategy implemented after around two weeks is less effective in preventing hypertrophic scarring than an early treatment.
Tissue oxygen levels can be increased by the use of hyperbaric oxygen therapy, an approach that has been used to treat chronic hypoxia . We simulate this treatment at various times by returning the proliferation rate of fibroblasts and the rates of production of collagen and TGF-β to the corresponding values for the normal healing case. As seen in figure 5, we find that the application of a treatment causes the ECM density to approach the corresponding value for the normal case. Furthermore, figure 6 shows that on applying treatment, the wound re-epithelializes at a rate closer to that of a normal wound.
Wound healing has been the subject of many mathematical investigations over the last two decades. These models have tended to concentrate on either epidermal or dermal wound healing; however, it is known that the crosstalk between epidermal cells (keratinocytes) and the primary cells of the dermis (fibroblasts) is critical to successful wound healing. The model we have developed here is more general than earlier models, as it explicitly takes these interactions into account. For the sake of simplicity, we assume that the fibrin clot is present initially, and so we do not model its formation, which occurs rapidly after wounding. As human wounds (healing by second intention) generally tend to heal by infilling rather than contraction, we do not include mechanical effects in our modelling. We note, however, that mechanical phenomena may play a role in severe hypertrophic scarring.
Our predictions, based on the numerical solutions of a coupled system of governing reaction–advection–diffusion partial differential equations, capture many key aspects of the wound healing process. In addition to reproducing a healing timeline (including the removal of a fibrin clot and re-epithelialization) that is qualitatively similar to the observed wound healing process, our model can simulate healing pathologies. To this end, we provide two illustrative examples of sub-optimal healing that lead to excessive collagen production and to a delay in healing, respectively.
We simulate the case of an excessively long inflammatory response by assuming that the decay rate of TGF-β is small. Our model predicts that this situation can give rise to hypertrophic scarring via the excessive production of collagen by fibroblasts, and that this outcome can be most effectively prevented by implementing a healing strategy effected within two weeks.
Sufficient local oxygen levels are required in order to facilitate normal fibroblast proliferation and collagen production , and TGF-β production. We simulate hypoxic conditions by changing the values of the parameters that correspond to these processes from their level in the normal healing situation by experimentally determined amounts. The significant delay in healing that we observe in the case of chronic hypoxia is also in agreement with the literature [13,31]. Our model predicts that hypoxic wounds can be successfully treated with the application of supplemental oxygen, which is consistent with recent reviews .
This model could be extended in ways that may better allow us to capture the healing process. For instance, the sharp profiles of keratinocytes (see figure 3) arise owing to the initial rectangular geometry, and more realistic behaviour could be obtained by using initial conditions that more closely represent the geometry of a dermal wound. Additionally, in order to prevent potential scenarios in which the wound cannot re-epithelialize, the dynamic generation of the fibrin clot could be incorporated into the model. Future models of this process could also consider the competition for space that occurs between fibroblasts, keratinocytes and mature collagen, which we ignore in this model to avoid increasing its complexity. While every attempt has been made to obtain parameter values that are in agreement with previous experimental results, several parameters in our model have been estimated in order to obtain appropriate steady states, and to qualitatively simulate observed healing behaviour. A detailed parameter sensitivity analysis, while outside the scope of the current work, will help to further validate this model. Although we have not attempted to address the full complexity of the fibroblast–keratinocyte crosstalk in this work, we find that our model captures the significant consequences of this critical biological process.
The authors would like to acknowledge funding from the following sources for R.C.S. to visit Queensland University of Technology (QUT): the Institute of Health and Biomedical Innovation (under the Visiting Researcher Scheme), the School of Mathematical Sciences at QUT and Western Kentucky University. Computational resources and services used in this work were provided by the HPC and Research Support Unit, QUT. This research was supported under the Australian Research Council's Discovery Projects funding scheme (project number DP0878011).
- Received February 15, 2012.
- Accepted April 27, 2012.
- This journal is © 2012 The Royal Society