A THREE DIMENSIONAL MODEL OF WOUND HEALING: ANALYSIS AND COMPUTATION

This paper is concerned with a three-dimensional model of wound healing. The boundary of the wound is a free boundary, and the region surrounding it is viewed as a partially healed tissue, satisfying a viscoelastic constitutive law for the velocity v. In the partially healed region the densities of several types of cells and the concentrations of several chemical species satisfy a coupled system of parabolic equations, whereas the tissue density satisfies a hyperbolic equation. The parabolic equations include advection by the velocity v and chemotaxis/haptotaxis terms. We prove existence and uniqueness of a smooth solution of the free boundary problem, for some time interval 0 ≤ t ≤ T , T > 0. We also simulate the model equations to demonstrate the difference in the healing rate between normal wounds and chronic (or ischemic) wounds.

1. Introduction.Wound healing under normal conditions is a process consisting of four overlapping stages: haemostasis, inflammation, proliferation and remodeling [5,13,19].During haemostasis, which occurs immediately after injury, clotting factors are delivered by platelets to the injured site to stop bleeding.Platelets also release chemokines, such as platelet-derived growth factor (PDGF), which recruits blood-borne cells to the wound.During the inflammatory phase, macrophages migrate into the wound, remove necrotic tissue and kill infectious pathogens.They also enhance the production of growth factors secreted by platelets, such as vascular endothelial growth factors (VEGFs), to attract fibroblasts and endothelial cells towards the wound.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 serves as a bed for tissue repair, and also contributes to scar formation.During the remodeling phase, which may take 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 due to venous insufficiency [18,21].Chronic wounds represent a major public health problem affecting 6.5 million individuals in the United States.It is estimated that $25 billion is spent annually on the treatment of chronic wounds [18].
Oxygen perfusion depends on the formation of new blood vessels that move into the wound, and it plays a critical role in the healing process.There are several mathematical models of wound healing which incorporate the effect of angiogenesis [15,14,3,17].Mathematical models of angiogenic networks, such as through the induction of vascular networks by VEGFs [6,7], were developed by McDougall and coworkers [12,20,2], in connection with chemotherapeutic strategies.The role of oxygen in wound healing was explicitly incorporated in [3] and [17].In particular, it was demonstrated in [17] that enhanced healing can be achieved by moderate hyperoxic treatment.A recent model by Xue, Friedman and Sen [22], which builds on [17], includes also the velocity of the ECM and treats the wound boundary as a free boundary.This model was initially developed for twodimensional radially symmetric wounds, and the numerical results were shown to be in agreement with the experimental measurements for cutaneous wounds [16].More recently the model was extended to three-dimensional wounds [11].In the present paper we consider this three-dimensional model.We establish the existence and uniqueness of a solution for some time interval 0 ≤ t ≤ T , T > 0. We then proceed to simulate the model for both normal and chronic wounds.The simulation results focus on two variables: the decreasing region of the open wound, and the macrophage density at the free boundary, that is, at the boundary of the open wound.As in the two-dimensional case [22], the three-dimensional simulations demonstrate the difference in the healing process between normal and chronic wounds.
In Section 2 we introduce the three-dimensional mathematical model developed in [11].The model consists of two coupled systems of PDEs.The first system is a strongly elliptic system for the three components of the velocity v in the evolving partially healed region Ω t .This system involves the isotropic pressure P , which is a function of the ECM density ρ.The second system involves various cell densities and concentrations of growth factors, oxygen, and ρ.It is a parabolic system for 7 variables and one hyperbolic equation for ρ.All these equations include advection terms by the velocity v.
In Sections 3-4 we consider the free boundary problem for the strongly elliptic system for the three components of the velocity v, assuming that the pressure P is a given function; we call this system the reduced problem.In Section 3 we establish some properties of this system, and in Section 4 we prove existence and uniqueness of a solution of the reduced problem for some time interval 0 ≤ t ≤ T , T > 0. In Section 5 we use the results of Section 4 to prove existence and uniqueness for the full wound healing problem, i.e., for the two coupled PDE systems.The special case of axially symmetric wounds is described in Section 6. this case is simulated in Section 7 for both normal and chronic wounds.
2. Model equations.We introduce the following variables: w(t, x): concentration of tissue oxygen e(t, x): concentration of Vascular Endothelial Growth Factor (VEGF) -p(t, x): concentration of Platelet-Derived Growth Factor (PDGF) -m(t, x): density of macrophages f (t, x): density of fibroblasts n(t, x): density of capillary tips b(t, x): density of capillary sprouts ρ(t, x): density of the ECM v(t, x): velocity of the ECM The geometry of the wound.
The geometry of the wound is shown schematically in Figure 1.The open wound, W t , is given by W t = {(x, y, z) : Z(t, x, y) < z < 0, (x, y) ∈ A t }.The wound's boundary is Γ t = ∂W t ∩ {z < 0}, and the partially healed region is Ω t = D\W t , where is a cylindrical domain whose boundaries x 2 + y 2 = L and z = −H 0 border the normal healthy tissue.
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 shall model it as upper convected Maxwell (viscoelastic) fluid with isotropic pressure depending on its density.The mass conservation law for the ECM density ρ is ∂ρ ∂t 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 ∇ • σ = 0, or For compressible material the isotropic pressure is a function of the density, i.e., P = P (ρ), and we take where β, ρ 0 are positive constants, and with ε/ρ 0 1 is a C 2+α approximation to x + for any α ∈ (0, 1).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. [22], the first term on the left-hand side is very small, so after dropping it we obtain, Denote the free boundary by φ(x, y, z, t) = 0.By continuity, the free boundary points move with the velocity v at these points.Hence, the equation for the free boundary is We assume that the external force at the wound surface is zero and there is no surface tension; hence, σ • ν = 0. Hence, the boundary condition for v is, or If we assume that the free boundary can be written as z = Z(t, x, y), then φ = Z(t, x, y) − z, and The boundary conditions for v at the fixed boundaries are All the boundaries are characteristic curves for (2.1), thus no boundary conditions are imposed for ρ.The initial condition for ρ is ρ(0, x) = ρ 0 . (2.12) In summary, the equations for ρ and v are (2.1),(2.6), (2.9), (2.11) and (2.12), while the free boundary equation is (2.8).
In addition to the equations for ECM, the remaining variables listed above satisfy the following system of partial differential equations in Ω t (see [22,11]), ) ) ) ) ) where the χ i 's are chemotactic/haptotactic coefficients, and H(s) is the smooth approximation of the Heaviside function, 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.These facts should be incorporated by taking the profile of G e (w) as in figure 2; a similar profile should be assumed for G p (w).For simplicity we approximate such profiles by taking, We also take The profile of G e (w).
To complete the model we need to prescribe boundary and initial conditions for the above variables.On the wound's boundary Γ t , PDGF is secreted by platelets, but secretion rate diminishes as the wound closes.We take where ν is the outward normal vector, and |W t,z | is the area of the cross section W t,z of the open wound with the horizontal plane at depth z and time t.The function g is a monotone decreasing function of |W t,z | and g(0) = 0.The remaining variables satisfy (2.21) The boundary conditions on the fixed boundaries for p, e, w, m, f, n, b are where f 0 , b 0 are the constant densities of f, b for normal healthy tissues.We impose the following initial conditions: where b 0 (x), p 0 (x) have three continuous derivatives and satisfy the boundary conditions (2.20) -(2.24), and Ischemia is a condition where blood supply to organ or tissue is decreased as a result of constriction or obstruction of blood vessels.We associate chronic wounds with ischemia, and, accordingly, we change the source function B w and the boundary conditions on { x 2 + y 2 = L} and {z = −H 0 } as follows.
Notice that if γ i = 1 for i = 1, 2, then the fluxes across { x 2 + y 2 = L} and {z = −H 0 } are zero, which means a total cutoff of blood supply.Moderate ischemia is represented by an intermediate value of γ i .The above ischemic boundary conditions can be formally derived by homogenization [22,9].Definition 2.1.We shall refer to the system of (2.1), (2.6), (2.8), (2.13) -(2.19), with the boundary conditions (2.9), (2.11),(2.20)-(2.24), and initial condition (2.12), (2.25) as the normal wound model.The corresponding model when B w is given by (2.26) and the boundary conditions on the fixed boundaries are replaced by (2.27) for w and similar Robin conditions for p, e, m, f, n, b will be called the ischemic wound model.
In the following three sections we prove the existence and uniqueness of the solution for the normal wound model for a small time interval 0 ≤ t ≤ T , T > 0. The proof for the ischemic wound model is similar.

The reduced problem and its properties.
We first consider the subsystem (2.6), (2.9), and (2.11) (or, equivalently, (2.7), (2.10), and (2.11)) with the free boundary condition (2.8), and prove local existence and uniqueness.Later on we shall extend the proof to the complete normal wound model.Without loss of generality, we assume from now on that η = 1.
Definition 3.1.The system (2.6) with the boundary conditions (2.8), (2.9), (2.11), and the initial domain Ω 0 will be called the reduced problem.Definition 3.2.The system (2.6) in Ω t with the boundary conditions (2.9) and (2.11) when t is a fixed time will be called the static reduced problem.

3.1.
Properties of the static reduced problem.We first consider the static reduced problem and show that the system is strongly elliptic, satisfies the Supplementary Condition (see [1], p.39), the Complementing Boundary Condition [1, p.42], and that the only solution of the homogeneous system is zero.This will ensure (by [1]) that for any smooth boundary z = Z(t, x, y) and a smooth function P (t, x) (t fixed) the static reduced problem has a unique smooth solution.
Supplementary condition.For any linearly independent vectors Ξ = (ξ 1 , ξ 2 , ξ 3 ), Ξ = (ξ 1 , ξ 2 , ξ 3 ), and τ ∈ C, we construct the polynomial L(Ξ + τ Ξ ) in τ for the system (2.7), The Hence, the polynomial L(Ξ+τ Ξ ) in τ has three roots with positive imaginary part, namely, We begin with some computations.The (inward) normal vector to the free boundary is and the tangential vectors are Then, mod M + , (3.7) The matrix B (Ξ) corresponding to the principal part of the boundary operators in (2.10) is where I 3 is the 3 × 3 unit matrix.The adjoint matrix to l (Ξ) is where Ξ = (ξ 1 , ξ 2 , ξ 3 ) as before.Then for Ξ = m+τ n, we have m•n = 0, (Z x , Z y , −1) = 1 + Z 2 x + Z 2 y n, and where Substituting (3.5)-(3.7)into (3.12)we obtain where The system is said to satisfy the complimentary boundary condition if the rows of the matrix B (Ξ)adj l (Ξ) are linearly independent mod M + (τ ); that is, if where Qk j is the k-th row vector of the matrix Qj , then all the constants C k are zero.Suppose not all the C k are zero.Then and, in particular, det Q0 ≡ 0. However, which is a contradiction.This verifies the supplementary boundary condition.
The homogeneous problem.We claim that the homogeneous system for v, with the homogeneous boundary conditions ) has only the trivial solution v = 0. To prove it, we multiply (3.14) by v j , sum over index j, and integrate over the domain Ω t .We obtain where ν = (ν 1 , ν 2 , ν 3 ) is the outward normal at the boundary of the domain Ω t .On the boundary {z = −H 0 } ∪ { x 2 + y 2 = L}, we have v j = 0, by (3.16); on the boundary {z = 0, (x, y) ∈ A t , x 2 + y 2 < L}, we have ν = (ν 1 , ν 2 , ν 3 ) = (0, 0, 1), and we have Similarly v 2 ≡ 0, and then This completes the proof that the homogeneous system has only the trivial solution v 1 = v 2 = v 3 ≡ 0.
Then Ṽ = ( Ṽ1 , Ṽ2 , Ṽ3 ), Φ also form a solution of (4.2) -(4.5).By uniqueness, Ṽ = V, Φ = Φ.Hence v = V| z≤0 satisfies the boundary conditions (2.11) at z = 0, and thus v and φ = Φ| z≤0 form a solution of the reduced problem.The solution of the reduced problem is unique, since two different solutions of the reduced problem can be extended (by (4.1)) to two different solutions of the extended reduced problem.
Remark 1. Lemma 4.2 does not actually depend on the fact that the free boundary is of the form z = Z(t, x, y); the reflection of the free boundary across z = 0 can be done with any parametric representation.In the next section we shall, more conveniently, represent the free boundary in the form x = X(t, θ, φ).

4.1.
Existence and uniqueness of solution for the extended reduced problem.In this section, we prove that the extended reduced problem has a unique solution, for a small time interval, with the free boundary parametrized as x = X(t, λ), where λ is orthogonal to the (x, y) plane at x = X(0, 0, φ). (4.6) Then the extension of the domain by reflection across z = 0 does not have corners at the points x = X(0, 0, φ).We further assume that for some α ∈ (0, 1), and where D * = {(x, y, z), 0 ≤ x 2 + y 2 < L, −H 0 < z < H 0 }.We introduce the local orthogonal unit vectors e r (λ), e θ (λ), e φ (λ) in the directions of increasing r, θ, φ, respectively.We write surfaces X(t, λ), λ ∈ Λ in a neighborhood of X(0, λ) in the form X(t, λ) = X 0 (λ) + h(t, λ)e r (λ). (4.9) and so that, by taking scalar product with e r , It will be convenient to express the free boundary condition (4.5) as a dynamic equation for h.To do that, suppose a point X(t, λ ) moves to position X(t + ∆t, λ) after a short time ∆t, where λ = (θ, φ), λ = (θ , φ ), and introduce the point λ = (θ , φ).Writing we obtain By the continuity condition (which states that the free boundary at X(t, λ) moves with the velocity V at X(t, λ)), and Taking the scalar product with e r and using (4.10), (4.13) and (4.14), we obtain Conversely, one can deduce from (4.19), or (4.20), the free boundary condition written in the form V ν = V•ν, where V ν is the velocity of the free boundary in the normal direction ν.Indeed, by (4.15) -(4.18) we have Given any family of curves X(t, λ) as in (4.9), we denote by Ω * t the domain bounded by X(t, λ) and the boundary of the rectangle D * , and write the extended system (4.1)-(4.4) in the form of where f 1 , g 1 are linear functions.We introduce a class of functions h(t, λ) by where For any h ∈ W T,M1,M , consider the problem (4.21), (4.22).From the results of Section 3 we know that for each t, the system (4.21),(4.22) is a strongly elliptic systems satisfying the Supplementary Condition and the complimenting boundary condition.Hence we can apply the Schauder estimates [1] to obtain the following lemma.
where C(M 1 ) depends on We next extend the function V(t, x) into D * such that For simplicity, we denotes the extended function also by V. We shall need to use the same extension procedure for any V corresponding to any h in W T,M1,M .We can achieve such an extension by first extending V along each normal to the curve X(t, λ) by a polynomial of degree two in the distance along this normal (to achieve smooth C 2+α extension) and then multiply this extension by a fixed cutoff C 3 function which equals to 1 in Ω * t and vanishes outside a small neighborhood of Ω * t , independent of t, 0 ≤ t ≤ T (T small).We define a function h(t, λ) by where V is the extension defined above.We note that the assumption (4.7) will be needed in order to prove that h(t, λ) is in C 2+α .We introduce the mapping Clearly h is a fixed point of S if and only if the corresponding V and X(t, λ) form a solution of the extended reduced free boundary problem.
To prove that S has a fixed point, we use the following lemma.
Then there exists a unique solution of (4.26) satisfying where C 1 (K) depends only on K.
where C 2 (K) depends only on K.
Proof.The proof of (i) is obtained by applying the proof of Lemma 2.2 in [4], or Lemma 3.2 in [8].The only part of the lemma which requires an additional argument is the proof of (4.28).To prove this estimate we set and proceed to estimate it in the same way as . in the proof of (i).
Since the characteristic curves through (t, λ) and (t + ∆t, λ) are "as close" to each other as the characteristic curves through (t, λ) and (t, λ + ∆λ) (up to a multiplicative constant), the only difference in estimating W and W is the fact that W (0, λ) = 0 whereas W (0, λ) ≡ 0. This however causes no problem since, by integrating the derivative of W (s, λ) along the characteristic curve, which passes through (∆t, λ), we find that and then, from (4.25) we then also get the bound here the C i (M 1 ) are constants depending only on M 1 .Taking T and M such that we conclude that S maps W T,M1,M into itself.We next show that S is a contraction, and hence, it has a unique fixed point in W T,M1,M .Take h 1 , h 2 ∈ W T,M1,M and the corresponding functions , and set δ = h 1 − h 2 T .Since the equations are not in the same domain, it is difficult to make comparisons.Hence we shall make a change of variable to reduce the problems to the same domain, Ω * 1,t by the transformation (r, θ, φ) → ( r, θ, φ), where where ψ is a C ∞ function such that ψ = 1 in a neighborhood of the initial wound boundary (free boundary) and ψ = 0 near the external boundary (fixed boundary) of the domain.Clearly, this transformation pulls the domain Ω * 2,t to domain at Ω * 1,t , and, the regularity of this transformation is the same as the regularity of the function h 1 − h 2 .Furthermore, since our domain excludes the origin, this transformation will not produce any singularities for the transformed PDE system.
We shall denote by V 2 the function V 2 under the above change of variables.On the domain Ω * 1,t , the boundary conditions of V 2 under the above change of variable are the same as V 1 with error ε where The quantity C(M, P 1 ) will henceforth be used to denote any constant which depends only on M and P 1 .V 2 satisfies the same elliptic system as V 1 except for an error term (which we conveniently view as an inhomogeneous term) whose C α norm is bounded by C(M, P 1 )δ.Hence, by Lemma 4.3, Using this estimate, we can then extend the definition of V 1 and We now write h1 and h2 in integrated form along the characteristics and begin to estimate h1 − h2 , ∇ λ ( h1 − h2 ), ∇ 2 λ ( h1 − h2 ), and its Holder coefficient, for fixed t.We obtain, by calculation similar to [4,8], the estimate then we can repeat the proof of Theorem 4.5, and conclude that there exists a unique solution with norm h * T ≤ M .In this proof we need to use the estimate (4.28) in order to estimate the (α/2)-Hölder coefficient of D 2 V in t, i.e., To do that, we change the domain Ω t for V(x, t ) to Ω t for V(x, t), in the same fashion as before, and use the C α/2 Hölder condition on D 2 λ h with respect to t. Going back to the original V(t, •) − V(t , •) we then obtain Hölder estimates with respect to both x and t.We summarize this by the following theorem.Theorem 4.6.Under the conditions (4.6), (4.7), and (4.31), there exists a unique solution of the (extended) reduced problem for some time interval 0 ≤ t ≤ T (T > 0) with free boundary (4.9) satisfying where the constant depends only on the initial surface X 0 (λ) and P 2 .

5.
Existence and uniqueness for the general case.We first extend the functions ρ, w, p, e, m, f , n, b by reflection across z = 0, and denote the extended functions by the same symbols.Because of the boundary conditions for w, p, e, m, f, n, b at z = 0, the extended functions will not have any singularities at z = 0.
As in the case of the reduced problem, we solve the extended normal wound problem.Take ρ in the space For any ρ ∈ Y T,M define P by (2.4) and use Theorem 4.6 to solve the extended reduced problem, with X(t, λ), h(t, λ), V(t, x).Then (5.1) We next extend V to all of the rectangle D * = { x 2 + y 2 < L, −H 0 < z < H 0 } as in the proof of Lemma 4.3.
Next we solve the parabolic system for Y = (w, p, e, m, f, n, b) with given ρ, X(t, λ) (or h(t, λ)), V(t, x).Using parabolic estimates as in [10], we derive the estimate We extend Y into D * as we did for V in the proof of Lemma 4.2.
We define ρ as the solution of with initial data (2.12).By Lemma 4.4 we obtain ρ, ρx C α,α/2 x,t We extend ρ into the rectangle D * in the same way we extended V in the proof of Lemma 4.3, so that the estimate (5.4) holds.
Hence the mapping ρ ) where To do that we use (5.5), (5.6) to estimate the difference between Y 1 and Y 2 as defined above (using also the zero initial data for Y 1 − Y 2 ): Using (5.4) -(5.7), we can estimate the difference ρ 1 − ρ 2 from the equation (5.3) which they both satisfy, Hence if T is small enough then W is a contraction.The fixed point of W provides the unique solution of the extended normal wound problem.We have thus proved the following theorem.
Theorem 5.1.Under the condition (4.6) -(4.7) there exists a unique solution for the normal wound problem with the free boundary satisfying the estimate (4.30) for some M > 0.
6.The axially symmetric case.We claim that the PDE system is invariant under rotation about the z−axis.To prove that, let A be a constant orthonormal matrix, i.e., A • A T = I 3 .We regard v and ∇ as row vectors.Then, written in matrix from, a rotation corresponds to a change of variables x = xA and the corresponding change of the derivative operator and velocity field are given by ∇ = ∇A, v = vA.Multiplying these equations on the right by A T , we obtain We denote the free boundary Γ t by z = Z(t, r) or by φ(t, r, z) = 0, where φ(t, r, z) = Z(t, r) − z.Since the velocity of the boundary is the same as v, the equation for the free boundary is φ t + v • ∇φ = 0 on Γ t , (6.7) or Z t + v 1 Z r − v 2 = 0. (6.8)The boundary conditions for v on the free boundary is derived from the relation σ•ν = 0, where ν = (−Z r , 0, 1)/ Z 2 r + 1 is the outward normal vector (pointing into W t ).In the axially symmetric case, σ is given by The boundary conditions for v on the fixed boundaries are v 1 = v 2 = 0 on {z = −H 0 } ∪ {r = L}, (6.11) ∂v 1 ∂z = 0, v 2 = 0 on {z = 0, R(t, 0) < r < L}, (6.12) v 1 = 0, ∂v 2 ∂r = 0 on {r = 0}, (6.13) if r = R(t, z) denotes the free boundary.Equation (6.12) implies that the shear stress τ rz = 0 for z = 0.
7. Simulations of the model.In this section, we show simulation results of normal and ischemic wounds in the axially symmetric case, taking the coefficients in the PDE system the same as in [22].We solved numerically the problem by the method of Arbitrary-Lagrangian-Eulerian (ALE) using COMSOL3.5a.In order to reduce the numerical error near the boundary z = 0, we solved the problem on the extended domain Ωt .Figure 3 gives the shape and macrophage density of a normal wound (γ 1 = γ 2 = 0) at day 1, 5, and 10 solved from the model.The initial wound is given by {(r, z) : √ r 2 + z 2 < 0.9 mm}. Figure 4 gives the shape and macrophage density of a ischemic wound with γ 1 = γ 2 = 0.6, also at day 1, 5, and 10.From these figures we see that the normal wound closes on a reasonable time scale, and macrophage goes away after day 5.However, the ischemic wound does not shrink significantly, and macrophages near the wound boundary persist over the full time course of simulation.

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

5 ) 4 . 1 .Lemma 4 . 2 .
on the free boundary Γt = {Φ = 0}.(4.Definition We shall call the system (4.2) -(4.5) the extended reduced problem.If the extended reduced problem has a unique solution (in some differentiability class), then the reduced problem has a unique solution (in the same differentiability class) .

FIGURE 3 .
FIGURE 3. Normal wound healing.The boundary flux function of PDGF is given as g(z) = k pb R(t, z) with R 0 = 0.9 mm.The boundary of the white region is given by r = R(t, z), −H 0 < z < H 0 , which is the wound boundary and its reflection across z = 0.The black curve indicates the initial position of the moving boundary, which is √ r 2 + z 2 = R 0 .The color of each plot gives the macrophage density.L = H 0 = 1.8 mm.