Modelling the Effect of Matrix Metalloproteinases in Dermal Wound Healing

With over 2 million people in the UK suffering from chronic wounds, understanding the biochemistry and pharmacology that underpins these wounds and wound healing is of high importance. Chronic wounds are characterised by high levels of matrix metalloproteinases (MMPs), which are necessary for the modification of healthy tissue in the healing process. Overexposure of MMPs, however, adversely affects healing of the wound by causing further destruction of the surrounding extracellular matrix. In this work, we propose a mathematical model that focuses on the interaction of MMPs with dermal cells using a system of partial differential equations. Using biologically realistic parameter values, this model gives rise to travelling waves corresponding to a front of healthy cells invading a wound. From the arising travelling wave analysis, we observe that deregulated apoptosis results in the emergence of chronic wounds, characterised by elevated MMP concentrations. We also observe hysteresis effects when both the apoptotic rate and MMP production rate are varied, providing further insight into the management (and potential reversal) of chronic wounds.


Introduction
Wound healing is a physiological response to injury of tissue involving the coordinated interactions of many cell types and biochemical agents (Wallace et al. 2021;Martin 1997). In individuals such as diabetes patients, so-called 'chronic wounds' may persist, which require medical treatments to promote healing (Nyugen et al. 2016;Fan et al. 2021). There are approximately 2.2 million people in the UK suffering from chronic wounds, costing the NHS over £5 billion per year (Department of Health and Social Care 2021). These wounds last, on average, 12-13 months and recur in 60%-70% of individuals, potentially leading to the loss of function (Frykberg and Banks 2015). Improved understanding of the biochemical mechanisms underpinning chronic wounds is therefore crucial to support the development of new treatments.
The healing of wounds involves the complex interplay of various cell types and their mediators and cytokines (Flegg et al. 2015). A wound is defined as damage to the skin which is comprised of two layers: the epidermis and the dermis. The epidermis is the outermost layer and is responsible for protecting against infection, while the dermis is the innermost layer and provides tensile strength for the skin by means of the dermal extracellular matrix (ECM) (Flegg et al. 2015;Frantz et al. 2010). The process of wound healing can be categorised into four successive stages: haemostasis, inflammation, proliferation and remodelling (Wallace et al. 2021). A wound that is able to complete these four stages in a well-coordinated manner is defined as an acute wound; conversely, a wound that spends a prolonged time in any of the four stages is defined as a chronic wound (Nyugen et al. 2016). The inflammation and proliferation stages of wound healing involve the production of matrix metalloproteinases (MMPs) by many cells types such as keratinocytes, fibroblasts, endothelial cells, and inflammatory cells (Guo and DiPietro 2010;Nyugen et al. 2016). For the remainder of this work, we consider only the production of MMPs by fibroblasts. MMPs are responsible for lysing protein components of the ECM, allowing for fibroblast migration within the ECM and, as a result, the proliferation of cells and modification of tissue (Nyugen et al. 2016). Experimental data has shown that there is increased expression of MMPs at the edge of a wound in all stages of the healing of an acute wound (Krejner et al. 2016). It is essential to have a well regulated concentration of MMPs during this process: if MMP levels are too low, this leads to uncontrolled ECM production and can cause issues such as hypertrophic scarring or dermal fibrosis (Li et al. 2009). If MMP levels are too high, chronic wounds are observed, as the overexposure of MMPs results in the increased degradation of the ECM (Sabino and Keller 2015). One such cause for these elevated MMP levels is deregulated apoptosis of dermal tissue, which may be caused by uncontrolled blood sugar levels in diabetes patients (Rai et al. 2005;Arya et al. 2014).
Mathematical modelling of the wound healing process can provide a framework to understand the behaviour of chronic wounds, potentially to direct treatments and improve the quality of life of patients. Many mathematical models have been developed to describe the wound healing process, each focusing on specific mechanisms of interest. For example, agent-based models (ABMs) can emulate the individual-level stochastic nature of the biological processes involved in wound healing (Sun et al. 2009;Mi et al. 2007;Walker et al. 2004;Ziraldo et al. 2013;An et al. 2009;Bankes 2002). Although ABMs are able to detail specific properties of individual cells, using them to model dynamics of wound healing on a tissue level may be computationally infeasible. We therefore will consider continuous models of wound healing, which typically use partial differential equations (PDEs) and present opportunities for a range of analytical techniques for further study (Flegg et al. 2015).
Continuous mathematical models of both dermal and epidermal wound healing typically take the form of reaction-diffusion systems, extending the work of Sherratt and Murray (1990) that examined the interaction of epidermal cells with a representative chemical acting as a regulator of mitosis. Their model has subsequently been developed to consider additional chemical agents involved in epidermal wound healing, such as the role of epidermal growth factor in corneal wound healing and the role of keratinocyte growth factor in epidermal wound healing (Sherratt and Murray 1992;Dale et al. 1994;Sheardown and Cheng 1996;Gaffney et al. 1999). Contrastingly, other models of dermal wound healing consider the role of fibroblasts in restoring the ECM in response to various mechanisms at different stages in the wound healing process. Many of these models consider the roles of angiogenesis and consider the role of oxygen transport in wound healing (Landen et al. 2016;Chaplain and Byrne 1996;Pettet et al. 1996a, b;Schugart et al. 2008;Flegg et al. 2015), while others consider the process of ECM restoration by fibroblasts in response to growth factors (Dale et al. 1997;Wearing and Sherratt 2000), the crosstalk of the epidermis and the dermis to simulate the simultaneous healing of both layers (Menon et al. 2012;Menon and Flegg 2021) and the incorporation of collagen fibre orientation due to cell movement (Dallon et al. 1999(Dallon et al. , 2000(Dallon et al. , 2001McDougall et al. 2006;Cumming et al. 2009). Travelling wave solutions are often observed in these mathematical models, agreeing with the experimental evidence of wound healing assays (Maini et al. 2004a, b). However, these models are often limited to the healing of acute wounds.
In this work, we consider the role that MMPs play in wound healing, with particular consideration given to chronic wounds. We propose a two-component reactiondiffusion model that describes the interaction of MMPs with dermal tissue, in which the restoration of the ECM is affected by MMP concentration levels. Using biologically realistic parameter values, this model gives rise to travelling waves corresponding to a front of healthy cells invading a wound. From the arising travelling wave analysis, we observe that deregulated apoptosis results in the emergence of chronic wounds, characterised by elevated MMP concentrations. We also observe hysteresis effects when both the apoptotic rate and MMP production rate are varied, providing further insight into the management (and potential reversal) of chronic wounds.
We develop our mathematical model in Sect. 2. Direct numerical simulations indicate the existence of travelling wave solutions, motivating a travelling wave analysis in Sect. 3. We then consider the effect of deregulated apoptosis on the healing of a wound via bifurcation and phase plane analysis in Sect. 4. Finally, we consider the effects of elevated MMP production levels on the healing of a wound in Sect. 6, before discussing our results in Sect. 7.

Model Development
We begin by presenting a mathematical formulation to describe MMP concentration dynamics in the presence of a wound. Elevated levels of MMPs play a key role in the persistence of chronic wounds; in light of this, we focus on constructing a reactiondiffusion model of wound healing that incorporates the interaction between MMPs with dermal cells during the wound healing process. For simplicity, we consider a wound with a one-dimensional Cartesian geometry, where x represents the longitudinal distance across a fixed domain, which evolves over time t. This model consists of two populations within a wound: the dermal cell density, denoted as n(x, t), and the MMP concentration, denoted as m(x, t). We take n(x, t) to incorporate the actions of both fibroblasts and ECM for the sake of simplicity, since fibroblasts are responsible for creating the components of the ECM.
We assume that dermal cells undergo mitosis with intrinsic growth rate σ and carrying capacity K n , noting that n = 0 corresponds to the complete absence of cell tissue. As mentioned previously, MMPs can assist in the healing of the wound up to a threshold concentration, denoted here as m thresh , above which they adversely affect healing of the wound by degrading the surrounding ECM. This effect is incorporated in the mitosis process by means of a function f (m). Furthermore, dermal cells undergo apoptosis with rate δ n and dermal cell motion within the wound is represented by linear diffusion with effective diffusivity D n .
MMP production by fibroblasts, denoted by a recruitment function g(n), occurs in response to inflammatory chemicals which are found in abundance in a wound. MMP production is therefore assumed here to be proportional to the absence of healthy tissue. MMPs undergo natural decay with rate δ m and diffuse through the wound with diffusivity D m . These aforementioned processes are shown in Fig. 1. Combining these modelling elements provides the following system of PDEs: We now discuss the functional forms of f (m) and g(n) appearing in the system (1), (2). The function f (m) in (1) represents the contribution to wound healing by MMPs and is defined to be an increasing function of m until m = m thresh , at which f (m) = f max , the maximum factor by which mitosis can be enhanced due to the presence of MMPs. Larger concentrations of MMPs hinder the wound healing process and thus, for m larger than m thresh , f (m) is a decreasing function of m. A suitable choice of f (m) is therefore a schematic of f (m) is shown in Fig. 2a. We note that f (m) and hence the parameter f max are dimensionless. For elevated levels of MMPs, f (m) decreases to zero and the mitosis term reduces to logistic growth with zero healing contribution due to MMPs.
One could also allow f (m) to take negative values for large m, representing a negative healing contribution to the ECM. We have considered a functional form of this type in Appendix A, where we show that such a choice gives rise to largely similar qualitative features to those produced by the functional form of f (m) given in (3). Given these similarities, we restrict attention to the choice of f (m) as in (3) as this choice results in simpler analytic computations. The recruitment function g(n) in (2) represents the production of MMPs by the cells in proportion to the absence of healthy tissue. Consequently, we can represent g(n) as where k m is the maximum reaction rate and γ m is a Michaelis-Menten constant (Michaelis and Menten 1913). A schematic of g(n) is shown in Fig. 2b.

Initial Conditions and Boundary Conditions
The initial conditions of the system (1), (2) are chosen as follows to emulate a wound starting at x = λ ν , i.e the tissue surrounding the wound occupies the spatial interval [0, λ ν ), and the wound is the region [ λ ν , L]. Furthermore, a small concentration of MMPs, η, is initially present at the wound edge λ ν : In (5), H is the Heaviside function while n = K n corresponds to the cells being at carrying capacity. Assuming the wound to be symmetric about its centre and that the tissue is fully healed far away from the wound, we adopt zero-flux boundary conditions on both endpoints of the finite domain of length L, i.e:

Nondimensionalisation
In this section, we non-dimensionalise the system (1)- (7). We are interested in the dynamics of the system on the timescale of mitotic rate of dermal cells and therefore introduce the following dimensionless variables: With these scalings, we obtain the following dimensionless equations: where In taking ν = σ D n in (5), (6) for simplicity, the dimensionless initial and boundary conditions are: and where η = η m thresh and L = L σ D n . We note that the dimensionless values λ, η and L are chosen for illustrative purposes. The corresponding dimensional parameter values are easily computed and are consistent with physiological values.
We define the dimensionless system (9)-(14) as the 'wound healing model' and we note that 0 ≤ N (X , T ) ≤ 1, since we do not expect the cell density to exceed its carrying capacity once healing is completed.
Suitable values for the dimensional parameters appearing in (1)-(4) are given in Table 1; we remark that these parameters correspond to 'normal', or healthy, wound healing and, unless otherwise stated, are employed in all simulations. These parameter values are used to calculate the dimensionless parameter values in (11)-(14) and are given in Table 2. We remark that where parameter values were unable to be obtained from literature; these were estimated heuristically.  Table 2. a Evolution of cell density N (X , T ) across the spatial domain, b Evolution of MMP concentration M(X , T ) across the spatial domain

Model Simulations
In order to obtain numerical solutions of the wound healing model, we use the method of lines and discretise on a one-dimensional spatial domain, using second order finite differences to approximate spatial derivatives. We then integrate the resulting system of time-dependent ODEs using MATLAB's ODE15s solver.
An example numerical simulation of the wound healing model is shown in Fig.  3. As seen in Fig. 3a, the cell density N (X , T ) evolves into a travelling wave where a front of tissue of density N ≈ 1 is 'invading' the wound (N = 0). This state of N ≈ 1 is close enough to unity to be considered as a healed state; we elaborate on the classification of a healed state further in Sect. 3. Since we observe a travelling wave where a healed front of tissue invades the wound, this scenario corresponds to a wound healing to completion. In Fig. 3b, we also observe the formation of a travelling wave for M(X , T ) in which a 'spike' in MMP concentration is found at the edge of the healing wound. These qualitative phenomena are in agreement with biological literature, which suggests the increased expression of MMPs at the edge of an acute would at all stages in the healing process (Krejner et al. 2016).

Steady State and Travelling Wave Analysis
From the previous section, we observe that the wound healing model gives rise to travelling wave solutions. Therefore, in this section we conduct a travelling wave analysis of the model. To determine the far-field states of the travelling waves, we first examine the uniform steady-states of (9), (10). In particular, the invading far-field state for N (X , T ) will indicate whether or not a wound is healing to completion. The uniform steady states (N , M) of (9), (10) satisfy the following: Combining (15), (16), we find that all non-trivial steady states may be written in terms of the single variable N : Using the fimplicit function in MATLAB, we determine a numerical solution for N and subsequently M using equation (15).
It is clear from (16) that we always have the trivial steady state, i.e. (N , M) = (0, 0) for all positive parameters. As discussed in Sect. 2.3, we characterise a wound as healing to completion when a far-field state of N ≈ 1 is achieved. This is because for δ N = 0, N = 1 satisfies (16). For non-zero δ N , however, N does not attain the value one, since apoptosis of cells is continually occurring. For small δ N , i.e. values corresponding to healthy biological functioning, (16) provides N close to unity, and we hence consider this state to describes healed tissue. Additionally, if an unhealed or partially healed steady-state of N , i.e. not qualitatively close to one, invades a state corresponding to damaged tissue, we characterise this as a chronic wound.
The stability of the uniform steady-states (N , M) is given by linear stability analysis of the analogous spatially-independent form of the wound healing model: This stability analysis (detailed in Appendix B) is used to determine the stability of the branches of bifurcation diagrams which we examine in Sects. 4 and 6.
In order to verify the existence of travelling solutions to the wound healing model, we consider a travelling wave analysis. Employing the travelling wave coordinates ξ = X − cT ∈ R, where c is the wave speed, the wound healing model (9)-(14) reduces to the boundary value problem (BVP): with boundary conditions We note that when considering solutions in travelling wave coordinates, the stability of N and M in the bifurcation diagrams in Figs. 5 and 9 are swapped, i.e. stable  Table 2 branches are unstable in travelling wave coordinates and vice versa, due to the change of variables ξ = X − cT .
Using MATLAB's BVP5c solver, we present numerical simulations of the BVP (20)-(22) in Fig. 4 with parameter values given in Table 2, which verifies the simulations of the wound healing model in Fig. 3. We are also able to use BVP5c to obtain a nondimensional wavespeed of c = 9.6, which translates to a dimensional wavespeed of c dim = 1.504 × 10 −2 mm h −1 . This wavespeed is of the correct order of magnitude observed from wound healing assays from (Maini et al. 2004a).

Deregulated Apoptosis
As mentioned previously, a potential cause for the emergence of chronic wounds is deregulated apoptosis. In this section, we examine the behaviour of the wound healing model under variation of the apoptotic rate of cells, δ N . Figure 5 shows bifurcation diagrams for N and M against δ N ; the stability of the branches is determined by linear stability analysis that is detailed in Appendix B.
Inspection of the bifurcation diagrams in Fig. 5 show that we have four bifurcation points that separate five parameter regimes, in which solution behaviours of the wound healing model differ qualitatively. The intervals for each regime are denoted as follows: The positions of the interval boundaries β 1 < β 2 < β 3 < β 4 and the method used to calculate them are given in Appendix C. While each interval is discussed in greater detail in subsequent subsections, we begin with a brief overview of the qualitative behaviours observed in each region.
For δ N ∈ I 1 , we have two steady-states for N : a state approaching unity and the trivial state. As discussed in Sect. 2.3, this parameter regime corresponds to a wound  Table 2. The bifurcation diagrams show four bifurcation points β 1 , β 2 , β 3 and β 4 separating five regimes I 1 , I 2 , I 3 , I 4 and I 5 . Dashed lines represent unstable branches and solid lines represent stable branches. The stable branches are colour-coded such that the green branch of steady states of N in a corresponds to the green branch of steady states of M in b, and similarly for blue branches healing to completion, with an increased expression of MMPs at the wound edge. For δ N ∈ I 2 , I 3 , it is unclear from the bifurcation diagram whether or not a wound will heal to completion, due to the existence of other non-trivial steady states. The bifurcation diagram for N indicates that we have 3 and 4 nontrivial states for I 2 and I 3 respectively: with one being a healed state (N ≈ 1), and the others being partially healed states. At δ N = β 3 , we observe a saddle node bifurcation point at which the healed state is destroyed. For δ N ∈ I 4 , there are two non-trivial states for N , both being unhealed states which indicates that for δ N ∈ I 4 , a wound will not heal to completion and hence a chronic wound will persist. Finally, for δ N ∈ I 5 only the trivial state exists.

ı N ∈ I 1
The bifurcation diagram in Fig. 5 indicates that we have two steady-state solutions, with one being the trivial state. The other state is N ∼ 1 at leading order, which can be verified using a regular perturbation analysis by considering the case where δ N , δ M 1. As discussed in Sect. 3, this regime describes the invasion of a wound by a healed front of cells and thus the healing of an acute wound. These features are reflected in direct simulation of the full system in Fig. 6a, for which δ N ∈ I 1 . We also note that travelling waves for M(X , T ) in Fig. 6b display an increased expression of MMP concentration at the wound edge, which is in agreement with the biological literature as discussed in Sect. 1. From the phase plane diagram in Fig. 6c, we see that regardless of initial conditions, trajectories in (N , M) space always arrive at the healed state, (N , M). Therefore, for δ N ∈ I 1 and other parameter values as in Table  2, a wound will always heal to completion.

ı N ∈ I 2 , I 3
For δ N ∈ I 2 , I 3 , the bifurcation diagram in Fig. 5 shows that we have 3 and 4 nontrivial states, respectively. The phase plane diagrams in Fig. 6c indicate that the initial state determines whether or not the wound heals to completion, with high initial concentrations of MMPs resulting in the emergence of a chronic wound. Contrastingly, if a small concentration of MMPs is initially present at the wound edge, a wound can heal to completion (Fig. 6a). This suggests that baseline levels of MMP must be kept low in order to minimise risk of the persistence of a chronic wound. The simulations for δ N ∈ I 2 , I 3 exhibit behaviour similar to that for which δ N ∈ I 1 , demonstrating the healing of an acute wound.

ı N ∈ I 4
For δ N ∈ I 4 , all trajectories arrive at one of two unhealed states for all initial conditions (Fig. 6c). These states have complex eigenvalues, characterised by spirals in the phase plane diagram and we therefore expect a travelling wave corresponding to a chronic wound with fluctuating MMP levels about its nonzero state. This is verified by the simulations in Fig. 6a, b, where we observe travelling waves with qualitatively different features to those observed in previously discussed intervals.
This change in the qualitative features of the travelling waves corresponds to a saddle-node bifurcation point at δ N = β 3 observed in the bifurcation diagrams in Fig. 5. At this bifurcation point, the branch of healed states N ∼ 1, shown in blue (and the corresponding branch for M) is destroyed. This loss of a healed state of N indicates that the invading state in the travelling wave simulations will transition to an unhealed state, which is verified by Fig. 6a where we also observe increased density of dermal tissue at the edge of the wound. As discussed in Sect. 3, we characterise this as a chronic wound. The emergence of a chronic wound also corresponds to the invading state of M in the travelling wave simulations transitioning to a larger value as shown in Fig. 6b, suggesting that elevated MMP concentrations prevent the complete  N (X , T ) and b M(X , T ). The phase plane diagrams in c correspond to the spatially independent ODE system (18), (19). Blue lines represent nullclines, green lines indicate representative trajectories of the ODE system and red dots represent the steady states. The yellow highlighted section of the domain represents the set of all the initial conditions (ICs) in the (N , M) phase space whose trajectories go to the healed state, the orange section represents the ICs whose trajectories go to an unhealed state and the red section represents the ICs whose trajectories go to the zero state. All other parameter values are as given in Table 2. We note that for I 2 and I 3 , only the fully healed travelling waves (obtained using initial conditions as in (12), (13)) are shown in a and b. We also note that the M-axis changes in each row healing of the wound. These results are consistent with the biological literature in that deregulated apoptosis is a prominent cause of the emergence of chronic wounds, which are characterised by having raised MMP concentrations.

ı N ∈ I 5
For δ N ∈ I 5 , we find that N and M both tend to zero for all initial conditions. This collapse of the travelling wave is expected in the case where apoptosis of dermal cells is very high, as cells die faster than they can reproduce.

Multistability and Hysteresis
As discussed in Sect. 1, a potential cause for the emergence of chronic wounds is deregulated apoptosis which may be caused by uncontrolled blood sugar levels in diabetes patients. In this section, we therefore examine how the apoptotic rate, δ N controls the transition from a healthy to a chronic state. We also investigate whether or not a chronic wound may be reversed if the apoptotic rate were to be regulated, for example, by regulating blood sugar levels. As discussed in Sect. 4.3, a saddle-node bifurcation point occurs at δ N = β 3 , which results in the emergence of a chronic wound, where high concentrations of MMPs prevent the wound healing to completion. The subcritical nature of this bifurcation has important implications for the solution behaviour, which we demonstrate in direct simulations of the wound healing model by allowing δ N in (9) to slowly increase in time: where δ N = 0.1 ∈ I 1 . As we see in Fig. 7, the travelling wave for N (X , T ) transitions from a completely healed wound to a chronic wound with elevated MMP concentrations, which verifies the results outlined in Sect. 4.3. As mentioned in Sect. 4.3, a saddle-node bifurcation point also occurs at δ N = β 1 . As such, we consider the case of decreasing δ N to investigate the implications for the solution behaviour upon decreasing δ N below β 1 from a regime in which a chronic wound persists. We take a value of δ N such that δ N > β 3 (to ensure the persistence of a chronic wound), and allow δ N in (9) to slowly decrease in time via (24) by changing with − . In Fig. 8, we take δ N = 2.5 ∈ I 4 and simulate the wound healing model with δ N (T ) applied to (9). From Fig. 8c, we observe the emergence of a chronic wound as we expect, displaying consistent qualitative features as in Fig. 7c when considering increasing δ N . For β 1 < δ N < β 3 , i.e. δ N ∈ I 2 , I 3 , however, we observe another partially-healed wound in the travelling wave simulations in Fig. 8b in contrast to Fig. 7b and demonstrates the existence of two qualitatively different travelling wave profiles, depending on the choice of initial conditions. We observe from Fig. 8a that for a chronic wound to heal, δ N (t) must be decreased to when δ N < β 1 , i.e δ N ∈ I 1 , demonstrating the subcritical nature of the bifurcation point β 1 , and the multistable nature of different travelling wave profiles for β 1 < δ N < β 3 . The multistability of Fig. 7 Simulations of the wound healing model with (24) applied to (9), i.e. δ N increasing in time with δ N = 0.1 and = 0.01. Other parameter values are as given in Table 2. Note that the X -axis has been trimmed to highlight the qualitative features of each travelling wave (dotted lines) but has not been scaled, i.e. the X -axis is the same length scale as those in Fig. 6a and b. The bifurcation diagram is as given in Fig.  5a for δ N ∈ [0, 3] and coloured asterisks represent the δ N (T ) values at the time points given in a, b and c. The green arrows represent the evolution of the invading states of the travelling waves as δ N is increased. Note that all travelling waves shown connect to the zero state. We also note that the M-axis changes at each time point the travelling wave profiles occurs as there is history dependence upon increasing and decreasing the apoptotic rate, i.e. a hysteresis loop affects the invading states of the travelling waves profiles. We conclude that chronic wounds may persist if apoptotic rates reach an unregulated state, but these chronic wounds may heal to completion if δ N were able to be regulated from a deregulated state, for example by controlling blood sugar levels.  (24) and → − applied to (9), i.e. δ N decreasing in time with δ N = 2.5 and = 0.005. Other parameter values are as given in Table 2. Note that the X -axis has been trimmed to highlight the qualitative features of each travelling wave (dotted lines) but has not been scaled, i.e. the X -axis is the same length scale as those in Fig. 6a, b. The bifurcation diagram is as given in Fig. 5a for δ N ∈ [0, 3] and coloured asterisks represent the δ N (T ) values at the time points given in a, b and c. The green arrows represent the evolution of the invading states of the travelling waves as δ N is increased. Note that all travelling waves shown connect to the zero state. We also note that the M-axis changes at each time point

Elevated MMP Production
As described in Sect. 1, chronic wounds are characterised by elevated levels of MMPs which prevent its healing. MMP production is a factor that may potentially be targeted with treatments for chronic wounds. In a similar fashion to the varying of δ N ,  Table 2. Dashed lines represent unstable branches and solid lines represent stable branches. The stable branches are colour-coded such that the green branch of steady states of N in a correspond to the green branch of steady states of M in b; equivalent for blue branches we now examine the qualitative changes to travelling wave profiles when we vary the production rate of MMPs, α.
We present bifurcation diagrams for N and M against α in Fig. 9; the stability of the branches is determined by linear stability analysis detailed in Appendix B. Inspection of the bifurcation diagrams show that we have two saddle-node bifurcation points (μ 1 and μ 2 ) separating three parameter regimes. We note that at μ 2 (at which α is large), the branch of healed states of N is destroyed and hence, similar to in Sect. 4.3, we expect the emergence of a chronic wound. We also note that the unhealed branch in Fig. 9 consists of states much higher than those in Fig. 5 upon varying δ N . This contrast is due to our choice of f (m) in (3). If we allowed f (m) to take negative values for large m, then the unhealed branch would take much lower values. Similar to Sect. 5, we now investigate the dynamics of the wound healing model for both increasing and decreasing α. Firstly, we allow α in (10) to slowly increase in time: where α = 0 < μ 1 . As we see in Fig. 10, the travelling wave for N (X , T ) transitions from a completely healed wound to a partially healed wound with elevated MMP concentrations (see the bifurcation diagram in Fig. 9b). These dynamics occur as α exceeds a value μ * ≈ 103.2 determined heuristically, where, despite the invading state of the travelling wave attaining a healed value, the wave visits the unhealed state, resulting in a partially healed wound. Increasing α further results in the invading state transitioning to an unhealed value, and thus the emergence of a chronic wound at μ 2 . We deduce that the threshold MMP production rate at which wounds cease to heal to completion is μ * , rather than the saddle-node bifurcation point μ 2 as we may have originally expected. We now consider decreasing α to investigate whether or not a chronic wound may be reversed if MMP production levels were able to be regulated. In Fig. 11, we take  Table 2. Note that the X -axis has been trimmed to highlight the qualitative features of each travelling wave (dotted lines) but has not been scaled, i.e. the X -axis is the same length scale as those in Fig. 6a and b. The bifurcation diagram is as given in Fig. 9a for α ∈ [0, 300] and N ∈ [0.8, 1] and coloured asterisks represent the α(T ) values at the time points given in a, b, c and d. The green arrows represent the evolution of the invading states of the travelling waves as α is increased. The maroon dotted line represents the value α = μ * . Note that all travelling waves shown connect to the zero state. We also note that the M-axis changes at each time point α = 300 > μ 2 and allow α in (10) to slowly decrease in time via (25) by changing with − . From Fig. 11a, we observe the emergence of a chronic wound as we expect, showing consistent qualitative features as in Fig. 10d when considering increasing α. For μ 1 < α < μ 2 however, we observe wounds that do not heal to completion in contrast to the case of increasing α, again demonstrating the existence of multistable travelling wave profiles depending on the choice of initial conditions. We also observe a change in qualitative features at α = μ * . For α < μ * , the travelling wave visits the healed state resulting in a partially healed wound. From Fig. 11d, we deduce that α(T ) must be reduced such that α < μ 1 for a chronic wound to heal.

Fig. 11
Simulations of the wound healing model with (25) and → − applied to (10), i.e. α decreasing in time with α = 300 and = 1. Other parameter values are as given in Table 2. Note that the X -axis has been trimmed to highlight the qualitative features of each travelling wave (dotted lines) but has not been scaled, i.e. the X -axis is the same length scale as those in Fig. 6a and b. The bifurcation diagram is as given in Fig. 9a for α ∈ [0, 300] and N ∈ [0.8, 1] and coloured asterisks represent the α(T ) values at the time points given in a, b, c and d. The green arrows represent the evolution of the invading states of the travelling waves as α is decreased. The maroon dotted line represents the value α = μ * . Note that all travelling waves shown connect to the zero state. We also note that the M-axis changes at each time point Similar to in Sect. 5, there is a hysteresis loop which results in history dependence of the travelling wave profiles on increasing and decreasing α which results in multistable travelling wave profiles for μ 1 < α < μ 2 . Unlike in Sect. 5 however, we observe an inconsistency in qualitative features of travelling waves for μ 1 < α < μ 2 due to the existence of the threshold value μ * . We conclude that that wounds will only partially heal beyond α = μ * . Furthermore, a chronic wound may be reversed to a healed state if MMP production rate is controlled from a deregulated state back to a regulated state.

Conclusions and Discussion
In this work, we develop a two-variable reaction-diffusion model, describing the interaction of matrix metalloproteinases (MMPs) with dermal cells in the wound healing process. In particular, we focus attention on the emergence of chronic wounds, since biological literature suggests that elevated levels of MMPs play a key role in their emergence. Our mathematical model gives rise to travelling wave solutions and in particular is able to emulate key qualitative features of the wound healing process reported in biological literature. One such property is that under parameter regimes representing healthy biological functioning, we observe acute wounds healing to completion with an increased expression in MMP concentration at the edge of the healing wound.
To contrast acute wounds with chronic wounds, we consider the effect of varying the apoptotic rate of dermal cells. Small apoptotic rates correspond to healthy biological functioning and thus we observe the complete healing of a wound. Provided initial concentrations of MMPs are kept low, an increase in apoptotic rate also results in the healing of a wound, which suggests that in order to avoid risk of the emergence of a chronic wound, baseline levels of MMPs must be kept minimal. Elevated apoptotic rate beyond a threshold value leads to the emergence of a chronic wound, with increased levels of MMPs invading the wound which prevents its healing. Moreover, we observe mulitstable travelling waves depending on initial conditions due to the existence of a hysteresis loop in the bifurcation diagrams for apoptotic rate. Therefore, in order to reverse a chronic wound to a state of complete healing, the apoptotic rate must be decreased below a threshold value. This gives insights into the regulation of apoptotic rate in the healing of chronic wounds which may be achieved by, for example, regulating blood sugar levels in diabetes patients.
We also consider the effect of varying MMP production rate on the healing of a wound. Similar to the analysis for apoptotic rate, we obtain threshold values at which chronic wounds persist and may be reversed. However, contrasting from the variation of the apoptotic rate, we observe a change in the the qualitative features of the travelling wave solutions at a different threshold value μ * . At this new threshold, occurring in the middle of the multistable regime, a wound transitions from healing completely to only being partially healed, before reaching a regime where chronic wounds persist.
Due to the observation of multistable travelling wave solutions with varying qualitative features, a natural extension to this work includes conducting a stability analysis of these travelling waves. Moreover, since we have considered the effects of varying apoptotic rate and MMP production rate separately, future work should consider the impact of concurrent defects in apoptotic rate and MMP production levels on the healing of a wound. We note that we have defined a chronic wound by changes to apoptotic rates and MMP production rates starting from a parameter set corresponding to acute wound healing. Parameter values corresponding to a chronic wound may however involve additional changes beyond these, hence further work should consider a full parameter sensitivity analysis.
Despite the wound healing model successfully being able to emulate qualitative features of both acute and chronic wounds, it has greatly simplified the process of wound healing by coupling the action of fibroblasts and the extracellular matrix, as well as assuming that MMPs directly contribute to healing of wounds. In reality, MMPs indirectly contribute to the healing by allowing the migration of other key elements involved in the wound healing process, such as fibroblasts, keratinocytes and endothelial cells (Nyugen et al. 2016). Future research should therefore take these factors into consideration. Furthermore, the role of tissue inhibitors of MMPs (TIMPs) in the wound healing process should be considered as these have an impact on the regulation of MMP concentration. This consideration would ultimately give further insight into the behaviour of MMPs during the wound healing process, as well as potentially direct biological therapies of chronic wounds, i.e. by determining the optimum physical and chemical composition of hydrogel therapies.  Table 2. a Evolution of cell density N (X , T ) across the spatial domain, b Evolution of MMP concentration M(X , T ) across the spatial domain which provides the following eigenvalues: Since δ M , δ N > 0, the zero state is an unstable saddle for δ N < 1 and is a stable node for δ N > 1. For all other steady states, we use MATLAB's fsolve function to determine (N , M) which are then substituted into (27) and whose eigenvalues are computed using MATLAB's eig function.