A MATHEMATICAL MODEL FOR CHRONIC WOUNDS

Chronic wounds are often associated with ischemic conditions whereby the blood vascular system is damaged. A mathematical model which accounts for these conditions is developed and computational results are described in the two-dimensional radially symmetric case. Preliminary results for the threedimensional axially symmetric case are also included.

1. Introduction.Chronic wounds represent a major public health problem affecting 6.5 million individuals in the United States.It is estimated that $25 billion is spend annually on the treatment of chronic wounds [1].Wound healing under normal conditions is a process consisting of four overlapping stages: haemostasis, inflammation, proliferation and remodeling [2,3,4].During haemostasis, which occurs immediately after injury, clotting factors are delivered by platelets to the injured site to stop bleeding.At the wound-site, platelets also release chemokines, such as platelet-derived growth factor (PDGF), which recruits blood-borne cells to the wound.During the inflammatory phase, mast cells release granules that contain enzymes promoting vascular leakiness.This enables neutrophils to migrate from the blood vessels into the wound site.Macrophages, differentiated from monocytes, also migrate into the wound, and, together with neutrophils, degrade and remove necrotic tissue and kill infectious pathogens.Macrophages also enhance the production of growth factors secreted by platelets, and other growth factors such as vascular endothelial growth factors (VEGFs), to attract fibroblasts and endothelial cells.The proliferative phase is characterized by the production of extracellular matrix (ECM) by fibroblasts, and by the directed growth and movement of new blood vessels (angiogenesis) into the wound.The newly deposited ECM on one hand serves as a bed for tissue repair, and on the other hand contributes to scar formation.During the remodeling phase, which may last several years, fibroblasts and other cells interact to increase the tensile strength of the ECM.Chronic wounds are those that fail to proceed through the above four stages timely and have persistent inflammatory stage, primarily due to venous insufficiency [1,5].In this paper

AVNER FRIEDMAN AND CHUAN XUE
we consider the first three stages of normal and chronic wound healing, i.e., the wound closure; the remodeling stage was considered by Dale et al. [6].
Oxygen plays a critical role in the healing process, and it depends on the formation of new blood vessels that move into the wound.There are several mathematical models of wound healing which incorporates the effect of angiogenesis [7,8,9,10].Mathematical models of angiogenic networks, such as through the induction of vascular networks by VEGFs [11,12], were developed by McDougall and coworkers [13,14,15], in connection with chemotherapeutic strategies.The role of oxygen in wound healing was explicitly incorporated in [9] and [10].In particular, it was demonstrated in [10] that enhanced healing can be achieved by moderate hyperoxic treatment.A more recent model by Xue, Friedman and Sen [16], which builds on [10], includes also the velocity of the ECM and treats the wound boundary as a free boundary.This model was developed for a two-dimensional radially symmetric wound.In the present paper, we extend the model to a general three-dimensional geometry, and provide initial simulations for the three-dimensional axially symmetric wounds.(m), oxygen (w), PDGF (p), VEGF (e), fibroblast (f ), extracellular matrix density (ρ) and velocity (v), capillary tips (n) and capillary sprouts (b).As illustrated in Figure 2, the wound occupies a region and the healed, or partially healed, region is Ω t = D\W t , where H 0 and L are such that W 0 lies above the plane z = −H 0 and inside the cylinder The ECM in Ω t is a growing collagen matrix which is elastic on a short time scale and viscous on a long time scale.We model it as upper convected Maxwell (viscoelastic) fluid with pressure depending on its density.The continuity equation for the matrix density ρ is where and k ρ , K wρ , ρ m , λ ρ are positive constants.The momentum equation is where σ is the total stress.We can write σ = −P I + τ where P is the isotropic pressure and τ is the deviatoric stress.Since healing is a slow, or quasi-stationary, process with negligible inertia, the last equation can be approximated by For compressible material the isotropic pressure is a function of the density, i.e., P = P (ρ), and we take where β, ρ 0 are positive constants.
For an upper convected Maxwell fluid, the stress-strain relationship is given by, where η is the shear viscosity, and, as shown in Xue et al. [16], the left-hand side is very small, so after dropping it we obtain, Hence, (3) is equivalent to, In summary, the functions ρ and v satisfy the system (1) -( 5) in Ω t .The (moving) wound's boundary is We denote the velocity of Γ t in the direction of the inward normal ν (pointing into In addition to (1) -( 2), the variables listed above satisfy a system of partial differential equations in Ω t : with specified structural functions B i and χ i .The explicit form of the B i 's and the chemotactic/haptotactic functions χ i are given below: The level of oxygen in the wound is a critical factor in the healing process.Moderate hypoxia improves healing; it stimulates macrophages to produce growth factors.Severe hypoxia impairs healing, since there is not enough oxygen for cells to grow and proliferate.Moderate hyperoxia improves healing, as it enables cells to proliferate faster.Extreme hyperoxia is toxic, and thus impairs healing.These facts are accounted for in the model equations.For example, in Equation (9), where w 0 is the oxygen concentration in a healthy tissue, k e and λ ei are constants, and G e has the profile shown in Figure 3. Ischemia is a condition where blood supply to organ or tissue is decreased as a result of constriction or obstruction of blood vessels.We denote the boundary of D in {z < 0} by ∂ 1 D, and suppose that blood supply is cut off from a portion ∂ 10 D of ∂ 1 D. Then on ∂ 10 D and on the boundary of Ω t which lies on z = 0, all the fluxes in Equations ( 7) -( 13) in the normal direction is zero.On the remaining boundary ∂ 11 D = ∂ 1 D\∂ 10 D the functions w, p, e, m, f , n, b take the same values as in normal healthy tissue.Under the above ischemic condition, the oxygen level w b in the vasculature (which appears in one of the terms of the function B w (w, f, m, p)) needs to be adjusted accordingly.
To complete the model we need to prescribe boundary conditions on the free boundary Γ t , and initial conditions.At the wound's boundary PDGF is secreted by platelets, but as the wound closes the secretion is diminished.We take where |W t,z | is the area of the cross section W t,z of the open wound (W t ) with the z-plane at time t.The function g is a monotone decreasing function of |W t,z | and g(0) = 0.For the remaining variables w, e, m, f , n, b, the flux in normal direction on Γ t is taken to be zero.Finally, the initial conditions on D\W 0 are the same as for a healthy tissue.For the parameter values for the system (1) -( 13), we refer to [16].
3. The radially symmetric 2-D case.We shall now specialize the model to the two-dimensional case, and furthermore, assume radial symmetry.Since the ischemic conditions as described above will break the symmetry, we shall implement the ischemic setup in a different way.
Suppose a function u satisfies where D δ consists of arcs of length δ on r = L, the distance between each adjacent δ-arcs is ε, and E ε = {r = L}\D δ .Then, by [17], if ε ∼ exp(−c/δ), δ → 0, while c is a positive constant, the functions u ε will converge to a solution of where α ∈ (0, 1) is a function of the constant c.Thus α represents the degree of ischemia: α = 1 means total cutoff of blood supply and α = 0 means no ischemia.Using this observation, which can be extended to the system (1) -( 13), we can now specialize this system to the two dimensional radially symmetric geometry, with wound {0 ≤ r < R(t)}, where we replace the conditions on ∂ 10 D and ∂ 11 D by a single mixed condition on all of ∂ 1 D. For example, w will satisfy In [18], the significance of ischemia on wound healing was experimentally addressed by a novel pre-clinical experimental model.Ischemic wounds were developed on a full-thickness bipedicle dermal flap where blood supply was isolated from underneath the flap and from the two long edges, as shown in Figure 4.One circular wound was then developed in the center of the flap (ischemic wound) and another on the normal skin (pair-matched non-ischemic wound) of the same animal served as control.
Figure 5 compares the experimental results obtained in [18] with the simulation results of our model.We see a very good fit after the initial 2-3 days, thus suggesting that since the radius of the wound in the experiments in [18] is small compared to the distance between the two parallel cuts (see Figure 4), a suitable choice of α can equivalently represent the level of ischemia in the three-dimensional geometry of the in vivo experiments.We note that since the wound contraction that occurs in vivo for the first few days is not included in our model, we cannot expect to get a good fit during this initial period of time.It was recently proved by Friedman, Hu and Xue [19] that (for any parameter values) in the radially symmetric two dimensional geometry the system (1) -( 13) has a unique global solution for any α ∈ [0, 1], and that if 1 − α is small (extreme ischemia), then R(t) = const.> 0 for all t ≥ T 1 for some T 1 = T 1 (α) > 0. Thus, in the case of extreme ischemia the wound does not heal.For the parameter values used in [16] it was further shown numerically in [19] that the radius R(t) ≡ R α (t) is monotone increasing in α for any time t.

4.
The three-dimensional axially symmetric case.The 2-D radially symmetric model described in Section 3 predicts quite well the change of the radius R(t) of ischemic wound on the surface z = 0 of the 3-D cutaneous wound developed in the experiments of [18].In order to determine the radius R(z, t) at depth z of an axially symmetric 3-D wound, we apply the model ( 1) -( 13) to a wound region In this case the homogenized boundary condition on r = L and at z = −H 0 (cf.( , where e r is the unit outward normal on r = L, e z is the unit outward normal on z = H 0 and α ∈ [0, 1] is the ischemic parameter.Similar boundary conditions are imposed for the other variables in (8) - (13).
Since the radius of the wound in the experiments conducted in [18] is small compared to the distance between the two long edges (see Figure 4), we conjecture that the radius R(z, t), which has not been measured yet in experiments, and the radius R(z, t) that will be computed by the mathematical model, for a suitable parameter α, will be in good agreement.Figure 6 gives the shape and macrophage density of a normal wound (α = 0) at Day 2 and Day 13 solved from the model.In the simulation, radial symmetry in the x and y direction has been imposed.The initial wound is given by {(r, z) : √ r 2 + z 2 < 0.75 mm}. 5. Future directions.Oxygen treatment of chronic wounds, either topical or in hyperbaric chamber have not yet achieved the desired level of effectiveness.Oxygen treatment can be introduced into our model by adding a control term Φ(w, t) to B w in Equation (7).We can then use our model to investigate which hypothetical treatment Φ(w, t) will best promote wound closure under ischemic conditions and this will suggest a biologically testable optimal protocol.

Figure 2 .
Figure 2. The geometry of the wound.

B
e (w, m, e, b, n) = k e mG e w w 0 − (λ en n + λ eb b + λ e )e, B m (w, m, p, b) = k

Figure 3 .
Figure 3.The profile of G e (w).

Figure 4 .
Figure 4. Schematic view of the wound environment developed by Roy et al. [18].The blue objects represent circulation barriers created using a bipedicle flap approach, and the red arrows illustrate the blood circulation near the wounds.Figure reprinted with permission.Copyright 2009 National Academy of Sciences, U.S.A.

Figure 6 .
Figure 6.Normal wound healing.The boundary flux function of PDGF is given as g(z) = R(z, t)/R 0 with R 0 = 0.75 mm.The black curve indicates the initial position of the wound boundary, which is given by √ r 2 + z 2 = R 0 .The color of each plot gives the macrophage density.L = H 0 = 1.5 mm.