A mathematical system for human implantable wound model studies

In this work, we present a mathematical model, which accounts for two fundamental processes involved in the repair of an acute dermal wound. These processes include the inflammatory response and fibroplasia. Our system describes each of these events through the time evolution of four primary species or variables. These include the density of initial damage, inflammatory cells, fibroblasts and deposition of new collagen matrix. Since it is difficult to populate the equations of our model with coefficients that have been empirically derived, we fit these constants by carrying out a large number of simulations until there is reasonable agreement between the time response of the variables of our system and those reported by the literature for normal healing. Once a suitable choice of parameters has been made, we then compare simulation results with data obtained from clinical investigations. While more data is desired, we have a promising first step towards describing the primary events of wound repair within the confines of an implantable system.


Introduction
All dermal wounds primarily heal via three fundamental mechanisms. These include the deposition of new connective tissue, epithielization and contraction (Diegelmann & Evans, 2004). Although these actions constitute a majority of the healing process, they are by no means the only factors involved in the repair. In order for the body to effectively mitigate the damage caused by an injury, it must also mount an immune response, and it must replenish the supply of blood to the wound through a process known as angiogenesis (Monaco & Lawrence, 2003). Together these mechanisms are critical for removing debris and invasive organisms from the wound as well as providing basic nutrients which support the various aspects of wound repair. Each of these fundamental processes takes place during four highly regulated yet overlapping phases. These phases are referred to as hemostasis, inflammation, proliferation and remodelling (Broughton, Janis, & Attinger, 2006;Monaco & Lawrence, 2003).

Mathematical models of wound repair
There are a wide variety of mathematical systems that have been proposed to describe the fundamental processes associated with the repair of a cutaneous wound. These models range in their level of complexity based upon the degree to which they incorporate key mechanisms associated with a given aspect of repair. For example, many models only account for a subset of the mechanisms involved with a particular wound healing process, such as fibroplasia or epithielization. Often these models require mathematical apparatus to account for the underlying mechanisms, such as cell mitosis and death, capillary tip formation, oxygenation, growth factor induced chemotaxis and mitogenesis, collagen matrix deposition and/or wound contraction. We provide a brief account of mathematical systems that have been proposed to account for the various aspects of wound repair. This is by no means a comprehensive list, but it provides a context of work that has already been established towards mathematically modelling cutaneous wounds. Thackham, McElwain, and Turner (2006) investigated several in silico approaches for solving different types of advection-dominated models. These models attempt to describe a subset of the events involved in angiogenesis. Although their work primarily focuses on the suitability of numerical techniques for solving the equations that arise from wound healing, they present several systems of partial differential equations (PDEs) which account for angiogenesis. These PDE systems model neovasularization by accounting for the development of either endothelial cell or capillary tip density in response to chemotactic mediators, such as VEGF. In their investigation, they analyze several wound healing scenarios in both one and two dimensions; however, they do not attempt to describe the repair process beyond the confines of neovascularization. Menke et al. (2010) have focused on a more complete repair model using an integrated system approach to describe the repair process in hypoxic environments. Their work focuses on the extension of a system of ordinary differential equations (ODEs) originally proposed by Reynolds et al. (2006). Reynolds et al. originally created their system to describe the acute inflammatory response at a systemic level. This system was modified by Menke et al. to include fibroblast activity and incorporate tissue oxygenation producing a wound repair system. The result consists of a set of ODEs which account for tissue damage, inflammation, fibroblasts and pathogen level. Although the model is able to account for the effects of low tissue oxygenation in the presence of varying degrees of contamination, it lacks a spatial component which describes the evolution of each species within the wound geometry during the repair process. Sherratt and Murray (1990), (1991) created a wound healing model consistent with the repair of the epidermis. They implement a Fickian diffusion model which incorporates stimulus for increase in epidermal mitosis through an activator/inhibitor mechanism. Their work illustrates that the biochemical regulation of mitosis through chemical activation or inhibition plays a key role in the repair process. Although simulations with their model compare well with empirical results, the model does not account for events which occur in full depth wounds. Vermolen, van Baaren, and Adam (2006) model a particularly problematic aspect of wound repair, namely wound contraction. Their work focuses on closure of the wound mediated through growth factor controlled tissue regeneration. In their work, a wound healing rate equation is combined with a reaction transport equation within an active layer adjacent to the wound where growth factor is produced. Tissue regeneration and expansion of the active layer proceeds if the concentration of growth factor exceeds a certain threshold value. Their work builds upon a set of mathematical models for wound healing first presented by Adam (1999). These models Adam (1999), Adam (1999), (2001) primarily focus on the Critical Size Defect (CSD) associated with a wound. CSD is defined as the smallest wound that does not heal within the lifetime of a test subject. Although their model accounts for this difficult aspect of repair, they do not account for other fundamental repair processes such as inflammation or tissue oxygenation which directly affect the rate of healing.
Wound healing's dependence on oxygen environments has also been studied using PDE models (Schugart, Friedman, Zhao, & Sen, 2008;Xue, Friedman, & Sen, 2009;Friedman and Xue, 2011). These systems account for several aspects of wound repair, including wound contraction. Their system of equations consists of a set of coupled of PDEs in the partially healed region modelling the wound edge as a free boundary. The free boundary moves with the velocity of the extracellular matrix (ECM) at the wound edge. The system of equations models the concentration oxygen and cytokines, such as PDGF and VEFG. The model also contains terms for macrophage, fibroblasts, capillary tip density as well as ECM velocity. Simulations of their model demonstrate that ischemic conditions limit macrophage recruitment to the wound site and impair wound closure.

Implantable systems
Although a wide variety of wound repair models currently exist, researchers have not proposed a wound healing model that describes the various aspects of repair within the context of a human implantable model. Clinical investigators have devised these systems in order to gain insight into mechanisms associated with abnormal repair in humans. This is because of animal models, which are often used in clinical investigations, lack insight into abnormal healing states, which plague humans. For example, Seoka et. al. have shown that acute inflammatory stresses in humans produced from different causes correlate poorly with mouse models (Seoka et al., 2013). Humans have a tendency to exhibit either excessive or deficient healing (Alaish et al., 2003). These extremes pronounce themselves in the form of keloids and hypertrophic scar formations or chronic ulcers, respectively. Implantable systems are porous cylindrical structures that are surgically placed within the cutaneous layer of the skin. In clinical investigations, multiple implantable structures are positioned at distinct locations throughout a test subject's appendage (Jorgensen, 2003). The surgical procedures used to position an implant are carried out under aseptic conditions in order to prevent infection. Hence, the immune response associated with these injuries is geared entirely towards the removal of cellular debris and damaged matrix components produced during implantation.
Once surgically positioned, implants are removed at distinct times points, at which time they are assayed for biological markers that characterize the phase of repair as well as the processes involved in the healing response (Alaish et al., 2003). Examinations often include differential cell counts, which determine the presence and quantity of neutrophils, macrophages and fibroblasts, as well as biochemical assays, which determine the amount of collagen deposited in the wound. Researchers have developed a variety of these in vivo models. The design of each of these systems has evolved over time. Although these implantable systems can be employed to examine different aspects of the healing process, they all conform to an overall cylindrical design as illustrated by Figure 1. Here, the inner region indicated by r a may be empty as in the case of the Schilling-Hunt steel mesh chamber, or it may be filled with viscose cellulose sponge (Cellstick model) (Alaish et al., 2003). In other synthetic implantable wound healing systems, the inner region may represent a contiguously integrated matrix composed entirely of polyvinyl alcohol or polytetrafluoroethylene (ePTFE) sponge (Alaish et al., 2003;Jorgensen, 2003). The typical dimensions of a tube are 1 mm in diameter by 6 cm in length.

The implantable PDE model
An acute dermal wound generally proceeds through a well-orchestrated sequence of events all of which lead to the restoration of normal anatomic structure and function. In this work, we seek to account only for a subset of these events. We primarily focus on the inflammatory response and events that occur during fibroplasia. In order to describe these events, our model consists of a set of four conservation equations, which describe the time rate of change of the density of damaged tissue, inflammatory cells, fibroblasts and newly deposited collagen matrix.
We use a number of well-established mathematical techniques to describe both the physical and biological processes involved in the repair process. The underlying mathematical building blocks we use include (1) The Logistic Equation -models population growth, (2) First order kinetic models -account for population death or decay rates, (3) The Mass Action Law -models interactions between reactive species, (4) The Keller-Segel model -accounts for cell chemotaxis.
Using these well-established modelling techniques, our wound repair system is as follows Here, m( x, t), n( x, t), f ( x, t) and c( x, t) represent the density of damaged tissue, activated inflammatory cells, fibroblasts and new collagen deposition, respectively, where x ∈ R 3 . In the following sections, we discuss the rationale behind each of these equations.

Damage marker equation
Damage may be interpreted as the adulteration of the wound by foreign material and the destruction of functional tissue, such as collagen matrix. We view tissue damage as fixed in space or having a timescale that is relatively slower than the chemotaxis of inflammatory cells and fibroblasts. Hence, we do not account for the diffusion of damaged tissue or foreign debris in the wound. In addition, we make the assumption that the injury is a result of a surgical insertion of the ePTFE tube which occurs under sterile conditions. Therefore, we do not account for the presence of pathogen in this model. It is the immune system's responsibility to remove debris from the wound during the healing process. Inflammatory cells, such as neutrophils and macrophages, carry out this function by engulfing debris and digesting it in lysosomes through a process known as phagocytosis. Inflammation also leads to harmful side effects. As inflammatory cells enter the wound to remove damaged tissue, they also damage healthy functional tissue. This side effect can become pronounced during abnormal repair where inflammation lasts longer. In our work, we do not model this side effect.
The damage equation consists of a conservation equation for initial tissue damage, m( x, t), stemming from the initial injury. From a conservation of mass perspective, the damage equation in our wound healing model may be expressed as follows Rate of removal of initial tissue damage = Rate of removal by phagocytosis .
We account for phagocytosis through a reaction-based mechanism. This process may be described by the following reaction m + n → n * , where m, n and n * represent the density of initial damaged tissue, inflammatory cells and inactive inflammatory cells. Essentially, this equation depicts the process whereby an activated inflammatory cell removes damaged tissue from the wound. The result is an inactive inflammatory cell. When we refer to the term n * , we mean an inflammatory cell that has succumbed to proteolytic processes that were vital in removing non-functional tissue. The Law of Mass Action states that the rate of removal of debris is proportional to the density of damaged tissue and inflammatory cells within the wound. This is represented in mathematical terms using the following differential equation: Here, μ mn is the consumption rate of damaged tissue by inflammatory cells. Because the rate of removal of m( x, t) has no spatial dependence with respect to the wound healing process, we can express this as ∂m ∂t = −μ mn mn.

Inflammation equation
Inflammatory cells, such as neutrophils and macrophages, typically respond to damage via autocrine and paracrine mechanisms. These mechanisms are mediated through the release of peptides referred to as cytokines. Two cytokines involved in wound repair at the onset of injury include PDGF and TGF-β. Inflammatory cells travel along gradients produced by these cytokines through a process known as chemotaxis. As these cells arrive in the wound, they begin consuming non-functional host tissue and foreign debris through phagocytosis. Eventually, neutrophils are removed by the macrophages in the wound. We view the overall rate of change of inflammatory cells in the wound as the sum of these two processes. In general terms, we have Rate of change of inflammatory cell density = Inflammatory cell chemotaxis − Inflammatory cell death from phagocytosis .
In our model, we account for inflammatory cell chemotaxis; however, we do not account for chemotaxis by way of chemical cues, such as PDGF and TGF-β. Instead, we assume these mediators, initially released by platelets at the onset of injury, form a gradient that is proportional to the initial tissue damage, m( x, 0). We assume this proportionality holds as the cytokine gradient changes with respect to time. In other words, if φ( x, t) represents the density of cytokines, then m = Kφ for all t > 0 where K is a positive constant. Hence, in our model, inflammatory cells respond to the evolution of the damage gradient, m(x, t).
We draw inspiration from the Keller-Segel chemotaxis model Hillen and Painter (2007), Keller and Segel (1970), (1971a), (1971b, to account for inflammatory cell chemotaxis. This system consists of the following two PDEs: Here, u denotes the cell density on a given domain ⊂ R 3 and v describes the concentration of the chemical signal. D v , k and μ represent the diffusion rate, rate of production and degradation of chemical, respectively. In Equation (6), χ is the chemotactic coefficient; D u is sometimes referred to as the motility or diffusion constant. In our work, we couple the chemotaxis equation, Equation (6), to the damage equation, Equation (1), to model the immune response to tissue damage. The inflammatory chemotaxis system initially appears as where D n and χ n represent the diffusion and chemotaxis constants, respectively. During phagocytosis, neutrophils eventually succumb to the same proteolytic processes that are responsible for the breakdown of foreign material and cellular debris within the wound. Additionally, macrophages inactivate themselves by releasing cytokines as a result of their phagocytic activity (Lang, Rutschman, Greaves, & Murray, 2002). To account for the deactivation of inflammatory cells, we implement the Law of Mass Action once again. Consequently, the inflammatory cell death term due to deactivation resembles Equation (1) and is given as follows ∂n ∂t = −μ nm mn.
We should note that the death rate of inflammatory cells governed by the constant μ nm is different from the rate constant μ mn shown in Equation (1). Its units are the same; however, since a single inflammatory cell may engulf numerous quantities of cellular debris, we expect μ nm < μ mn . That is the rate of deactivation of inflammatory cells occurs at a slower rate than the removal of damaged tissue from the wound. The overall rate of change of inflammatory cells is the sum of Equations (8) and (9). Thus, we arrive at the overall time rate of change of inflammatory cells in the wound, which we restate below for convenience ∂n ∂t = ∇ · D n ∇n − nχ n ∇m − μ mn mn.

Fibroblast equation
As inflammation subsides and inflammatory cells reduce in number, fibroblasts begin to enter the site of injury. This marks the beginning of the proliferation phase. During this period, fibroblasts respond to the same cytokine gradients as inflammatory cells, and they react to cytokines signals in two ways. First, cytokines provide chemical signals for fibroblast chemotaxis. Hence, these chemical cues influence fibroblast cell migration to the wound. Second, cytokines serve as important mitogenic agents causing fibroblasts to increase their rate of mitotic cell division. We account for both of these mechanisms in our wound repair model as well as the natural loss of fibroblast cells. Thus, the mass balance of fibroblast cell density is expressed as As with the inflammatory equation, we do not directly account for fibroblast chemotaxis via cytokine gradients. Rather, we assume that cytokines, which are produced by neutrophils and macrophages during inflammation, form gradients that are proportional to the density of these inflammatory cell populations. Therefore, fibroblasts respond to the evolution of the inflammatory cell gradient in our model. Once again, we use the Keller-Segel system to describe fibroblast chemotaxis. Here, the time rate of change of fibroblasts due to chemotaxis becomes where D f is the diffusion or motility constant associated with fibroblasts, and χ f is the chemotaxis constant. In this equation, D f and χ f serve the analogous roles that D n and χ n provided in the inflammation equation.
To account for fibroblast cell growth, we use a logistic model. The Logistic equation has long been used to describe the growth rate of populations. It has the following general form where k is a constant that depicts the growth rate and U 0 is the population's maximum sustainable value. All cell populations such as fibroblasts eventual experience apoptosis or cell death (Greenhalgh, 1998). In order to account for the overall life cycle of fibroblasts within the wound, we add the first-order decay term, −μ f f , to account for cell death. Therefore, the time rate of change of fibroblasts associated with these two terms appears Here, μ f is a constant that describes the death rate of fibroblasts, and F 0 is the nominal density of fibroblasts in healthy tissue. There is a subtle difference between the logistic equation and its implementation in Equation (10). In Equation (10), we introduced the term (2−f /F 0 ) rather than (1−f /F 0 ). This factor was introduced to model the equilibrium condition of the fibroblast equation. That is we desire the mitotic generation rate and death rate of fibroblasts to balance as t → ∞. We will discuss this subtlety in more detail in Section 2.5. Since cytokines influence the mitogenic response rate of fibroblasts, the growth rate would normally be a function of cytokine concentration, or in our case, its proxy, inflammatory cell density. Hence, s f (n) serves as a variable mitogenic growth rate in Equation (10); however, for simplicity, we let s f (n) = k f . Combining both chemotaxis and cell life cycle terms produce the overall rate equation for fibroblasts during repair. This has the following form:

Repair equation
The density of newly deposited collagen and granulation tissue is a direct indication of the progress of the repair. Any system that attempts to explain dermal wound repair must account for these variables. Thus, we monitor the status of these two species with a single unknown, c( x, t). Prolonged inflammation can damage newly deposited collagen and granulation tissue. In our modelling approach, we account for damage to newly deposited tissue due to this mechanism. Consequently, our conservation equation results from the mass balance of these two terms. The overall time rate of change of the density of newly deposited tissue during wound healing is expressed as follows Overall rate of new tissue deposition = Fibroblast tissue deposition rate − Tissue damage rate from phagocytosis .
First, we consider the rate of new collagen deposition to be proportional to the density of fibroblasts and collagen damage, c, in the wound. This gives where k c is a positive constant which describes the rate of collagen deposition. Collagen damage, c, can be expressed as the difference between the existing collagen density, C 0 , in undamaged tissue and undamaged collagen, c, within the wound. In other words, we write Substituting Equation (12) into (11) provides the rate of collagen deposition in the wound. This is given as Prolonged inflammation leads to tissue damage. Here, we account for damage to newly deposited collagen and granulation tissue by this mechanism using a mass action term. In our repair equation, we view the rate of damage to new collagen matrix and granulation tissue as proportional to the product of the density of inflammatory cells and newly deposited tissue in the wound domain. This yields the following damage rate term: where μ cn represents a positive rate constant for tissue damage due to inflammation. Combining Equations (13) and (14) yields the overall rate of repair or Equation (15).

The repaired state
In humans, an acute dermal wound only regains approximately 80% of its original structural integrity with regard to tensile strength; however, in our model, we assume that the wound fully recovers. What does this state look like in terms of our model? If repair proceeds in an ideal manner, we expect all damaged tissue to be removed from the wound during the inflammatory phase. Inflammatory cells respond to specific biological markers or antigens which once removed, provide no stimuli for inflammatory cell functionality. Thus, activated inflammatory cell levels should drop as the density of damaged tissue disappears from the wound. Ideally, damaged matrix and inflammatory cell levels should both approach zero as the system approaches the final repaired state.
On the other hand, as repair proceeds, we expect the rate of newly deposited collagen matrix to slow as the newly deposited collagen and granulation tissue approaches their nominal values in healthy or undamaged tissue. In our wound repair system, we model this value as C 0 , the density of healthy tissue. As collagen density recovers in the wound, fibroblasts enter an equilibrium where their overall rate of production and death due to apoptosis balance around their nominal levels in normal tissue, F 0 (Greenhalgh, 1998). Thus, in terms of all variables, we assume that the repaired state represents the point (m,n,f ,c) = (0, 0, F 0 , C 0 ). With regard to our system of equations, this is a special point or solution to our model, if Under these assumptions, if we let (m, n, f , c) = (m,n,f ,c), we have In other words, the repaired state represents a homogeneous steady-state solution in our system of equations.

Simulations
Implantable models provide experimentally controlled environments which aid in the development of mathematical models. Our goal is to create a mathematical system which describes aspects of wound repair within the framework of the synthetic implant. Consequently, we constrict the domain of our mathematical system so that it conforms to the geometry of the implant. This is illustrated by Figure 2. In Figure 2, R b is a constant that represents the boundary of the implant and r is a dependent variable such that 0 < r < R b . The domain, , is defined as the set {(x, y) : x 2 + y 2 ≤ R b }. From symmetry, we expect ∂u ∂θ = 0 and ∂u ∂z = 0, where u is a place holder for m, n, f , and c. This condition reduces both the gradient and the divergence operators as follows ∇u =r ∂u ∂r Consequently, our system reduces to a set of time-dependent equations with one spatial component, which are given as follows ∂m ∂t = −μ mn mn, where we now have functions m(r, t), n(r, t), f (r, t), and c(r, t) rather than m( x, t), n( x, t), f ( x, t), and c( x, t), respectively. Simulations were perform using a combination of PDEONE (algorithm 494) (Sincovec & Madsen 1975a, 1975b and the PDEPE solver provided by MATLAB.

Boundary conditions
At the left end point, r = 0, we impose the following symmetry condition for all equations.
At the right end point, r = R b , the boundary conditions for each equation are quite different from one another. For the damage equation, we impose the same type of boundary condition, but the rationale is different. Since the damage marker does not diffuse throughout the wound, we impose the condition This condition represents a zero-flow condition due to the fixed nature of the tissue damage. The right boundary condition for the inflammatory cell equation is difficult to model. Inflammatory cells are produced in the spleen, and they are activated thorough transport mechanisms via blood stream, which are beyond the scope of the model's domain. In our model, we do not model cytokines directly. Rather, we assume that gradient's of PDGF and TGF-β are proportional to the species that release them, namely inflammatory and fibroblast cells. To arrive at a suitable boundary condition, we must consider several factors. First, we consider the pool of resting inflammatory cells to be very large and their production rate to be under a quasi-steady state condition. Also, we assume the activation and transport of inflammatory cells to be very fast as compared to the timescale of wound repair. Hence, in this work, we assume that at the right boundary the injection of inflammatory cells into the wound is proportional to the quantity of damage at the boundary. We express this by the following relationship Another suitable choice of boundary conditions would be to set the density of inflammatory cells at the boundary proportional to the gradient of damage at the boundary. In other words We will explore this option in future studies. The damage caused by the initial infraction dissipates as one approaches the edge of the wound. Therefore, fibroblasts should approach their nominal levels as r → R b . For this reason, we approximate the density of fibroblasts at right boundary point to be equivalent to the nominal density of fibroblasts in healthy tissue. Thus, the right boundary condition becomes Since new granulation tissue and collagen deposited by fibroblasts is fixed, we impose the same no-flow condition as we did with right boundary point of the initial damage. Thus, we have following boundary condition for new tissue deposition: at the right boundary point.

Initial conditions
The initial conditions of our PDE system account for the state of the variables at the onset of injury. We assume that the profile of the initial damage marker is Gaussian, and it achieves a maximum at the centre of the wound or implant. The initial condition for damage is as follows m(r, 0) = M 0 e −C m r 2 .
Here, M 0 is the maximum density of damage, C m is the Gaussian decay constant. For our other two equations, we assume that their values at the boundary slough off into the wound. We use exponential forms to express this phenomenon. Here, we have These equations represent the initial conditions for inflammatory and fibroblast cells, respectively. For the former equation, n(R b , t) is the boundary condition for inflammatory cells at t = 0, and C n is the decay constant. The latter equation, which represents the initial condition for fibroblasts, has a very similar form. Here, F 0 is the density of fibroblasts in healthy tissue, and C f is the decay constant for fibroblast density.
Since C 0 represents the density of healthy tissue, then the initial wound profile, c(r, 0), may be represented as c(r,0) = C 0 -m(r,0), where c is the profile of collagen in the wound. However, we must be careful how we treat the time evolution of such an equation in our model because the removal of non-functional matrix is distinct from the deposition of a new collagen matrix. This leads to the following initial condition for the initial undamaged collagen profile.

Parameter values
One of the difficulties with mathematically modelling any biological system is the procurement of parameters which describe the fundamental properties of the system. The issue of biological variability and complexity make it difficult for biologists to make quantitative measurements of these parameters in the same way a physicist might measure a fundamental constant of nature. Consequently, there is currently very little information with regard to the fundamental properties associated with wound healing. For this reason, we attempt to elucidate values for coefficients in our mathematical system from simulations performed with our PDE model. Before we can use our model in this way, we first must identify a set of values that characterize the behaviour of a given unknown during the repair process. In addition, we must also define criteria so that we can arbitrarily evaluate how well a given solution predicts its characteristic information. For a given variable, we assign three characteristic time periods. These refer to a given unknown's behaviour during the repair process. We define each of these characteristic intervals as follows (1) Initial response time, t in = t in − t 0 : the time period in which an unknown begins to initially accumulate at a given location in the wound starting from the onset of injury or implantation, (2) Peak response time, t max = t max − t 0 : the time period in which an unknown attains its maximum value at a given location in the wound starting from the onset of injury or implantation, (3) Equilibrium or steady-state response time, t fn = t fn − t 0 : the time period in which an unknown approaches or becomes arbitrarily close to its steady-state value at a given location in the wound starting from the onset of injury or implantation.
Here, t 0 represents the time of injury or implantation.  Li et al. (2003), Mast and Schultz (1996) 1-2 Li et al. (2003), Mercandetti and Cohen (2013) 3 Li et al. (2003), Mercandetti & Cohen (2013) Macrophages 2 Diegelmann and Evans (2004), Li et al. (2003) 2-7 Li et al. (2003), Monaco and Lawrence (2003) 7 Broughton et al. (2006), Li et al. (2003), Monaco and Lawrence (2003) Fibroblasts 2-3 Greenhalgh (1998), Stadelmann, Digenis, and Tobin (1998 7 Li et al. (2003), Stadelmann et al. (1998) 14 Li et al. (2003) Collagen 7-10 Gailit and Clark (1994) 7 -1 0 Stadelmann et al. (1998) 14-21 Greenhalgh (1998), Li et al. (2003), Monaco and Lawrence (2003)  upon information obtained from our literature search. All information displayed in this table is reported in days. In Table 1, we have cross referenced each of the characteristic values with the source from which they were obtained. Using the data presented Table 1, we can construct a notional timeline which depicts the natural progression associated with each of our unknowns during repair. The figure illustrates this for neutrophils, macrophages, fibroblasts and collagen deposition. In Figure 3, all curves have been drawn to reflect their respective values tabulated in Table 1. The immune response reflects the sum of neutrophil and macrophage activity. We assume that active neutrophils and macrophages approach a steady-state value of zero. This is because once the inflammatory phase comes to an end, these cells are no longer required to remove debris from the wound. Figure 3 also outlines characteristics associated with fibroblast infiltration and collagen deposition. Since it has been reported that fibroblasts begin apoptosis once collagen deposition levels off (Greenhalgh, 1998), we assume that fibroblasts and collagen reach their equilibrium points congruently.  Figure 3 provides a valuable check when examining potential solutions that stem from our system of equations; however, we use the information presented in Table 1 to fit the parameters of our model. Once we fit our model with a suitable set of parameters (see Table 2), we validate the model and numerical solution with data obtained from ePTFE implants.

Simulation results
We now discuss the overall simulation results in the context of a wound which undergoes repair. Figures 4, 5, and 6 display the numerical results when these parameters are used to simulate wound repair within the geometric constraints of a synthetic implant. In Figures  4(a) and 5(a), we plot the normalized density of an unknown against the radius of the wound domain. Here, 0 ≤r ≤ 1, where 0 is the wound centre and 1 is the outer edge of the implant. The subplots shown in each figure represent a specific snapshot in time. As an alternative, we provide a three-dimensional view for each variable as repair proceeds. These are illustrated in Figures 4(b) and 5(b). These figures were generated by performing a three-dimensional rotation of their respective one-dimensional plots along the axis of symmetry. Snapshots for Figures 4(a) and 5(a) are provided for each time point illustrated in Figures 4(b) and 5(b), respectively. Each figure separates the solution into three distinct periods to illustrate the evolution of each variable described in our model. For convenience, we refer to these periods as the initial, intermediate and steady-state phases. The initial phase covers the first 24 h after injury. This period primarily captures the initial transient response of each unknown, particularly the inflammatory response. By comparison, the intermediate phase covers the first two weeks after injury. During this phase, unknowns still behave transiently; however, the proliferation response, rather than the inflammatory response, becomes the dominant response. Finally, the steady-state phase includes time points where the solution converges on its homogeneous steady-state solution, or repaired state. Figure 4 shows the response of the system within the initial 24 h directly following injury. Specifically, this plot displays the response for 0, 7.2, 14.4, 16.8, 19.2, and 24 h after injury. The first plot in the upper left hand corner of Figure 4 shows the initial conditions at the time of injury. In this plot, damaged tissue,m(r, t), increases as we approach the centre of the wound where we assume the injury occurred. The normalized density of inflammatory cells and fibroblasts,n(r, t) andf (r, t), respectively, exponentially drop off   as we move inwards from the edge of the wound,r = 1, into the region where tissue damage increases. For inflammatory cells, the value at the wound boundary represents the level of activated inflammatory cells which respond to the intensity of damage at the boundary. This is our boundary conditionn(r = 1, t) = κ * m (r = 1, t). In the case of fibroblasts, their values at the boundary represent their nominal levels in health tissue. Their levels fall off exponentially due to the initial damage asr → 0. The profile of undamaged collagen, c(r, t), is simply the complement of the damaged tissue profile.
During the initial transient phase, inflammatory cells appear to be the sole responders to the injury. This behaviour is expected, and it represents the inflammatory response associated with wound repair. During this period, inflammatory cells first diffusive into the wound at approximately 7.2 h or .3 days. After this point, these cells begin to migrate into the wound via a mechanism other than diffusion. This is the Keller-Segel chemotaxis model in action. As a result, inflammatory cells begin to accumulate at the wound's centre within 14.4 h or .6 days. The aggregation of inflammatory cells continues and becomes more pronounced by the end of the first day. As inflammatory cells migrate into the wound, they remove damaged tissue. During the initial phase, the removal is very subtle, but by the end of the first day, it becomes clearly apparent, particularly at the wound centre where the normalized density of tissue damage becomes less than unity.
The intermediate phase associated with our solution is illustrated in Figure 5. Here, we provide the simulations results for a two-week period following the initial 24 h post injury. In this figure, we examine the response of our four unknowns for days 2, 5, 7, 9, 11, and14. The inflammatory response continues to strengthen during the initial time points of this phase. Inflammatory normalized cell density reaches it maximum value during the simulation at day 2. As a result, the tissue damage profile reduces to approximately 50% of its initial value at the wound centre. In addition, we begin to see the first indications of fibroblast infiltration. In our model, fibroblasts respond to inflammatory cell gradients via the Keller-Segel model. Because inflammatory cells release factors associated with fibroblast chemotaxis, we assume that cytokine gradients may be approximated using normalized inflammatory cell density. Consequently, fibroblasts appear to aggregate near the inflammatory cell curve's inflection point. This is where the inflammatory cell gradient reaches its maximum value and is approximately the midway point between the wound's centre and edge.
At Day 5, the inflammatory response begins to subside but spreads outwards shifting its gradient in the direct the wound boundary. Fibroblasts respond to this spatial shift. As fibroblasts continue to migrate into the wound, they follow the inflammatory cell gradient. Their levels peak in the vicinity of this region, and they slowly diffuse from this point towards the centre of the wound. This process continues into days 7 and 9 where fibroblast levels peak. During this time frame, fibroblasts deposit new collagen matrix. As this new collagen approaches the nominal level in healthy tissue, fibroblast activity subsides. This process continues as both of these unknowns converge towards unity, which is depicted in the final plot shown in Figure 5.

Model validation
Although the numeric solution presented in Figures 4-6 is in reasonable agreement with the characteristic values reported in Table 1, it is important to note that the set of parameters, which produced this solution, were selected based upon their evaluation criteria. In order words, we selected the solution based upon its performance. This simply verified that our model successfully describes the ideal behaviour associated with the repair of an acute dermal wound. Now we answer the question: how well does this particular solution describe a clinical investigation carried out with the implantable system? Oswal, Belle, Diegemann, and Najarian (2013) recently reported on a wound repair study which implemented an expanded polytetraflouroethylene, ePTFE, sponge. The purpose of their investigation served to validate the results of an entropy-based automated cell nuclei segmentation and quantification algorithm. This routine automates the process of manually counting cell nuclei that have been histochemically treated for biological markers. In their examination, they compared results obtained from their algorithm with manual counts preformed by a pathologist. Their data-set consisted of 21 immunohistochemically stained images belonging to a single patient. Images were treated with Hematoxylin & Eosin (H&E), cluster of differentiation 68 (CD-68) or alpha-smooth muscle actin (α-SMActin) stains. Samples were acquired from human tissue sections derived from the ePTFE implant. Tubes were implanted subcutaneously into the upper arms of a healthy volunteer subject, and were removed at 5, 7, and 14 days after surgical implantation.
In this study, macrophages were identified by the presence of a specific cytoplasmic granule referred to as CD-68, and fibroblasts were characterized by Hematoxylin & Eosin (H&E) stain. In addition to this data, we have acquired results from this study (not published) which quantify the degree of collagen deposition in ePTFE cross sections courtesy of Dr. Robert F. Deigelmann. These images where treated with a trichrome stain, which is a specific marker for collagen. Table 3 summarizes the results obtained from these two sources.
According to Oswal, macrophages and fibroblasts were expressed in terms of the total number of cell nuclei per cross section of ePTFE sponge. Alternatively, collagen deposition in these cross sections was quantified by the total colour intensity of pixels. These measurements reflect total values with no spatial dependence. For this reason, we computed the total value associated with a given unknown over the entire domain at each time point using composite trapezoidal integration.
We then normalized a given data-set from either the clinical investigation or our simulation with its respective maximum value to account for discrepancy between units. Table  3 compares the final results of process. Here, the normalized values, which correspond to a particular unknown, are reported for each data-set for days 5, 7, and 14. In addition, we computed the absolute error between the respective normalized data-sets. We report these values in the last two columns of Table 3. Examination of these calculations indicates that our model accurately predicts the maximal response associated with fibroblasts and collagen deposition; however, there is only a qualitative level of agreement with the offpeak responses associated with fibroblasts. The model breaks down in its description of the inflammatory response. Comparison between the macrophage levels for both days 5 and 14 show the largest points of disparity among our model and the empirical data-set. For this case, it appears that macrophage infiltration lasts substantially longer than the 7-day period we reported in Table 1. Since we used this point as a characteristic value for our evaluation criteria, it is understandable why our results do not agree with these empirical points. Finally, our model's best performance occurs in its depiction of new collagen deposition. Here, the greatest disparity occurs at day 5 where we report an absolute error of 11%.
As a final comparison, we plot each empirical data point against its respective value obtained from the numerical solution. Figure 6 displays the normalized values associated with each data-set over a 16-day period. Here, we report total normalized levels of each variable that was recorded in Table 3. Figure 6 reports these values for macrophages, fibroblasts and collagen deposition, respectively.

Discussion
We have created a basic PDE model which describes the repair of an acute dermal wound within the confines of human implantable systems. In our approach, we primarily focused on the fundamental mechanisms associated with the inflammatory and proliferation phases. Our system of equations accounts for the time rate of change in the density of the initial damage, inflammatory cells, fibroblasts and new collagen. The underlying mechanisms associated with unknowns in our system of equations included the rate of mitotic cell division, cell death, collagen deposition, phagocytosis and chemotaxis. We implemented fundamental mathematical techniques to account for a majority of these processes, such as the first-order decay, the Law of Mass Action and the logistic model. In our model, the Keller-Segel chemotactic model accounts for the reaction-diffusion processes associated with inflammatory and fibroblast cell populations during repair. The coefficients associated with our PDE model were obtained through simulations with our system equations. We carried out simulations until there was general agreement between a given solution and a set of characteristic values that we obtained from the literature. This strategy appears to have produced promising results for an initial attempt towards modelling implantable systems.
Although our wound repair system predicts the general progression for most of the unknowns it models, there were cases where the model deviated from data-sets obtained from both theoretical and empirical perspectives. The inflammatory response appears to be the major exception in this last regard for both cases.
From a theoretical venue, the inflammatory response does not match all characteristic values obtained from our literature review. In our system, this response initiates within a reasonable period; however, inflammatory cell activity last appreciably longer than reported. For a typical acute dermal wound, the inflammatory phase usually subsides approximately seven days after the initial infraction (Broughton et al., 2006;Li et al., 2003;Monaco & Lawrence, 2003;Pallister & Watson, 2010). In our system, it is difficult to determine exact intervals of this type for all unknowns. This is because of the continuous behaviour associated with all variables in our model; however, according to our verification criteria, inflammatory cells maintain 10% of their peak value 13.4 days after the start of the simulation. This compares poorly with data reported within the literature for an acute inflammatory response. The situation might be corrected by the inclusion of a decay term that represents inflammatory cell death. We did not include this term in the inflammatory cell equation because we hypothesized that the death of inflammatory cells from phagocytosis should dominate their background decay rate.
When we compared the results produced by our system with empirical data obtained from ePTFE implants (Oswal et al., 2013), similar trends emerged between experiment and numerical solution. Our model accounts for the general behaviour of fibroblasts, myofibroblasts and collagen deposition for 5, 7, and 14 days after implantation. In this regard, there is particularly good agreement between the peak responses for these data-sets. The closest agreement between any data-set occurs with the deposition of new collagen. We computed the relative errors for collagen deposition for days 5, 7, and 14 as .11, .07, and .01, respectively. Although we had a reasonable prognosis for these unknowns, our model's prediction for macrophages (i.e. inflammatory cells) proved problematic. Data obtained from Oswal et al. showed that macrophages remained at appreciable levels 14 days after implantation. Although there are only two points in their data-set for macrophages, this point represents the peak value. In our model, inflammatory cells attain their accumulated peak value on day three. We selected parameters for our system of equations that would yield outcome consistent with data collected from our literature review (Broughton et al., 2006;Diegelmann & Evans, 2004;Gailit & Clark, 1994;Greenhalgh, 1998;Li et al., 2003;Mast & Schultz, 1996;Mercandetti & Cohen, 2013;Monaco & Lawrence, 2003;Stadelmann et al., 1998). Hence, this empirical observation is also in disagreement with those sources.
Human implantable models present tremendous opportunities for applied mathematicians interested in applying their skills in the field of biomathematics. Implants greatly simplify the mathematics required to describe the repair process because they minimize the modelling of complex features associated with wound healing. Furthermore, these empirical systems provide valuable data, which is collected in a controlled environment. Consequently, the results of these studies can be used to validate any hypothetical mathematical system which attempts to model certain aspects of the repair process, such as the inflammatory response or fibroplasia. These mathematical systems can in turn play a valuable role in the study of wound healing. For instance, mathematical models offer a noninvasive approach to study the dynamics of healing in human studies. These system allow for the study of hypotheses and therapies to be examined before clinical implementation (Reynolds et al., 2006). A numerical model can increase the success of clinical trials and aid in the understanding of disorders associated with abnormal healing, such as hypertrophic scar formation and chronic ulcers.
In this work, we have taken an important first step towards developing a basic mathematical system that describes the repair of acute dermal wounds within the context of implantable systems. In future work, we hope to incorporate additional features into our PDE model. These features would include accounting for the effect of low tissue oxygenation along with the inclusion of a pathogen term. In addition, we should also account for the side effects of inflammation, such as tissue damage. Furthermore, it would be advisable to include a cytokine equation rather than assume that gradients associated PDGF and TGF-α maintain levels of proportionality with the cell types that produce these cytokines. Finally, since the inflammatory response in our model has proved questionable, it might be advantageous to divide this response into two terms. The remedy would consist of separate equations for neutrophils and macrophages.