- Email: [email protected]
The effect of surface curvature on wound healing in bone: II. The critical size defect.
MATHEMATICAL AND COMPUTER MODELLING
PERGAMON
Mathematical and Computer Modelling 35 (2002) 1085-1094 www.elsevier.com/locate/mcm
T h e Effect of Surface C u r v a t u r e W o u n d H e a l i n g in B o n e : II. T h e Critical Size D e f e c t .
on
J . A . ADAM D e p a x t m e n t of M a t h e m a t i c s a n d S t a t i s t i c s O l d D o m i n i o n University, Norfolk, VA 23529, U.S.A.
(Received July 2001; accepted August 2001) Abstract--Two related models are studied for wound healing on a spherical surface. The first, and simpler model is based on a reaction-diffusion equation describing the process of healing as the inward growth of new cells from a circular wound on a spherical surface. Several aspects of the wound healing problem are considered: by interpreting the healing process as a "pseudowave" propagating across the spherical surface a heuristic account of the "speed" of healing is possible, and the corresponding healing time is characterized. Of paxticular importance in relation to animal models is the existeuce (or not) of a critical size defect (a CSD, defined below); this is discussed as a consequence of the stability of the steady states of the system to nonuniform spatial (or angular) perturbations. Explicit criteria are derived under which a CSD exists (within the model) in terms of the skull radius and wound radius. Although neither model explicitly distinguishes between wound healing in bone or tissue, it is in the former t h a t CSDs are known to occur. The second model invokes a weighted spatial average cell density which permits the presence of both a short-range activation t e r m (as in the first model) and a long-range inhibition term. Under these circumstances, within a suitable parameter range, the phenomenon of aggregation may occur in addition to the behavior predicted by the first model. It is suggested t h a t such aggregation is manifested in the case of keloid scarring, which can occur as a result of wound healing in tissue. (~ 2002 Elsevier Science Ltd. All rights reserved.
Keywords--Reaction-diffusion Keloid scar.
equation, Surface curvature, Critical size defect, Wound healing,
1. I N T R O D U C T I O N In an earlier paper [1], the time-independent distribution of wound-induced growth factor (GF) over a spherical cap was considered. This and related papers (see [2] and references therein) are associated with the existence of the so-called critical size defect, which has been defined as the smallest interosseous wound that does not heal by bone formation during the lifetime of the animal [3] (see also [2] for further clinical references). For practical purposes, this timescale can usually be taken as one year. In [4], the definition was further extended to a defect which has less than ten percent bony regeneration during the lifetime of the animal. CSDs can "heal" by fibrous c o n n e c t i v e t i s s u e f o r m a t i o n , b u t s i n c e t h i s is n o t b o n e , i t d o e s n o t h a v e t h e p r o p e r t i e s ( s t r e n g t h , etc.) that
a completely healed defect would.
Some typical CSD diameters
0895-7177/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S0895-7177(02)00073-0
a r e , for r a t , r a b b i t ,
Typeset by ¢4A/rS-TEX
1086
J.A. ADAM
dog, and monkey calvaria (skullcap), respectively: 8mm, 15 mm, 20mm, and 15 mm (details can be found in [3]). A phenomenological time-dependent problem is now examined with intent to determine, somewhat heuristically, the speed of wound healing, the time for wound healing to occur, and the existence of a critical size defect (CSD). The starting point is the logistic reaction-diffusion equation for the spread of new bone cells inward from a circular wound edge on the surface of a sphere of radius R. The distribution of new cells is considered to be azimuthally symmetric at all times. The initial wound angle subtended at the center of the spherical skull is 2~0 where 80 E (0, 7r/2); if total healing occurs, this angle decreases to zero (which corresponds to a polar angle of 7r/2 in standard notation, so 8 -- 0 is the "north" pole of the skull). Because of the azimuthal symmetry, the angular domain is 0 < 8 < 80 < zr/2. In terms of the bone cell density u(0, t), the governing partial differential equation is
Ou D ('O~u Ou) (i T Ot - - ~ \082 +cotOÕ-~ + a u _ u),
(1)
where D is the coefficient of diffusion for the GF, a is the linear growth rate, and K is the saturation limit for the bone cell density. The quantities D, a, and K are all constants in this model. However, it will be useful (in connection with the steady states ~ below) to allow these constants (and also 0) to be piecewise constant for the purposes of studying the evolution of the healing process. Eliminating the first derivative term via the change of dependent variable
u(8, t) = f(8)y(8, t) - (csc 0)1/2y(8, t),
(2)
equation (1) becomes
õ~=ZõB+~
l + ~ c o t 28 ~+«y
1---
,
(3)
where/~ = D / R 2. Homogeneous steady states ~ exist when ~ = 0 and
~{~(1 )}~.
=ä-]
2. P S E U D O W A V E
l + ~ c o t 28
SPEEDS;
+a
-
HEALING
(4) TIMES
By analogy with travelling wave solutions for infinite domains, we seek "pseudowave" solutions of the form y = y(~), ~ = - R S - ct. (5) This corresponds to the pseudowave moving in the direction of decreasing 8 with positive "speed" c. The qualification "pseudowave" is of course necessary because (i) the domain cannot be (-c~, oo), and also because of (ii) the presence of O-dependence in the new reactive term in equation (3). However, as noted by Murray [5] in connection with a mathematical model for calcium waves on amphibian eggs, some insight into mechanisms may still be obtained by using similar arguments for piecewise constant values of ~. The governing ordinary differential equation is
d2y dy a f , , D - ~ + C-d-~+ ---~y~y - y) = O.
(6)
A standard phase plane analysis for the (y(~), y~(~))-plane reveals that the nature of the equilibrium points (0, 0) and (y*, 0) is determined by the characteristic equation
a f . = 0, DA2 + cA =t: -~-y
(7)
T h e Effect of Surface C u r v a t u r e
1087
where the + and - signs refer to (0,0) and (y*,0), respectively. Thus, (y*,0) is a saddle point, and (0, 0) is a stable node provided that the pseudowave "speed"
{( ° ( ' D
Cw~amin:2
a+~-~-
1+~cot28
)))'~~
(s)
where ¢min is a constant for each fixed 8, where 0 < 8 < 80 < ~r/2 as noted above. This condition guarantees that the heteroclinic trajectory connecting these two steady states is the only bounded positive trajectory. Were the domain ( - c o , co), we would be completely justified in claiming that a "healing" wave moves from the stare y = y* to y = 0 as 8 decreases (and ~ increases). Thus, by analogy we consider the healed state to spread as a moving front throughout the (now) restricted domain. Note that Cmin increases for given R as 8 decreases, i.e., the wound heals laster the smaller is the angular extent of the wound. This also implies that if 8 is allowed to decrease in equation (8) in a piecewise constant manner, the healing process speeds up as the wound radius contracts. In addition, Cmin always exceeds its planar value 2x/-Dä for finite R, and increases away from this value as D increases. An upper bound for the time taken for a wound to heal (if indeed it does so) may be reasonably defined from this model as
Bth
L°
th :
dt-----
R
d8
o Cmin(8)
RSo = Th.
(9)
< Cmin(8)
The lower limit 80 is the initial (half-) angular size of the wound. Strictly, of course, this pseudowave analysis is based on c being constant, which it will be for a given value of 8; it would be of interest, therefore, to compare observed healing times in an experimental model with th and Th as calculated here. Based on studies of c in planar models, it appears likely that Cmin is the speed that will be attained for a given value of 8, though it is not clear whether the dependence on 8 will complicate the evolution of the pseudowave as it moves across the spherical cap. The integral can be evaluated to yield
th(R,8(0))
R
-
x/D(4a
L ° sin 8 + J) vv(0) X / # - 2 - - c o s 2 0
R
x/D( 4a + e)
{arcsin #
d8
- arcsin[# cos 80] },
(lOa) (10b)
where # =
\4äa ~ 2~,]
'
0 < # < 1.
(11)
For numerical purposes, the upper bound th can be expressed in dimensionless form as
4ath
x2
= ~(arcsin
and for comparison, the supreme
Th
#
arcsin[# cos 00])
(10c)
can be similarly expressed as
4aTh =
x2
~/x 2 + 2 + c o t 2 00
where x = 2Rv/-~~.
0o,
(10d)
1088
J.A.
ADAM
0.06 /
0.05 0.04
/
0.03• /
0.02.
./
/
II////
///
//
/ //i/
/
//
/
/
tl ¸
/
/
/
/
/
/
/
/
/
/'
/
~ " j"
/
0.01. tl ~
.I
_~'S~":
.lf" ........
. [ .........
..........-.-
0.2
,
__-.---.----
, i
0.4
0.8
0.6
Figure 1. The dimensionless upper bound on the healing time, 4ath(x), where x ---2Rv/äß, for fixed values of semiwound angle 00 (00 = 0.1, 0.3, and 0.5; these are the lower three curves in ascending order). The supremum Th for 00 ---- 0.5 is also shown for comparison (topmost curve). // 0.7" /
/
/
/
/
0.6
/ / /
/
0.5
/
/
/
0.4 0.3 /
0.2 0.1 0
//
//
_ y---
/
.j-
i
0.2
0.4
. . . .
0.6
I
,
~
0.8
I
1
Figure 2. The dimensionless upper bound on the healing time, 4ath(0O), for fixed values of x (x ---- 0.5, 1, 2, in ascending order; the supremum Th is not shown).
Graphs of 4ath(x) are presented in Figure 1 for fixed values of 00 (and the supreme Th in orte case is also shown for comparison). In Figure 2, graphs of 4ath(OO) are presented for fixed values of R (but Th is not shown). The excess Th -- th is greater for larger values of x (or R) and 0o than for smaller values, as would be expected for the case of larger wounds. In each case, the value of Th is about a factor of two larger than the corresponding value of th.
3. L I N E A R
STABILITY:
EXISTENCE
OF THE
CSD
We return to the governing differential equation (3) and explore the linear stability of the steady states ~ -- 0 and ~ -- y*. While it is customary to make modal expansions in terms of surface spherical harmonics (or in the case of azimuthal symmetry, as hefe, Legendre polynomials) [6-8], note that equation (3) is isomorphic to a one-dimensional diffusion equation, and a simple Fourier component of the form below will suffice for present purposes. The governing equation is, from
T h e Effect of Surface Curvature
1089
(3) and (4),
Oy
02y af , . = Zõõ~ + ~ y t y - y).
(12)
It is clear that spatially (or more accurately, angularly) uniform solutions satisfy a logistic equation, with the steady state ~ = 0 being linearly unstable, and ~ = y* being linearly stable to spatially uniform perturbations. Regarding nonuniform perturbations, let us consider stability to perturbations of wavenumber k by writing
Y = fl + eYle «t+iÆR° + 0 @2) . Substituting this expression into equation (12) for 9 = 0 and ignoring terms of O(e 2) and higher, the following dispersion relation arises, again valid for eaeh fixed value of 0 and R:
a=-Dk
«,
2 +-~y
=-Dk 2 +~
°(~
l + ~ c o t 20
)
+a,
(13)
so ~ = 0 is linearly stable, and no growth occurs, i.e., a CSD exists if
-~
> ~ l + ~ e o t 20 ,
(14)
i.e., for a sufficiently large wound radius r (= ROo) for a given R (assuming k 2 > a/D), or for a suffieiently large "skull" radius R for a given wound "angle" 00. One implication of this result, for a given value of k is that there may be partial healing only if there is a loeal deerease in the ratio a R 2 / D or (for constant a/D) an increase in the loeal radius of curvature in the wound environment. Note also that this result does uot imply that a larger wound radius will neeessarily result in a CSD as eompared with a smaller radius; given two spherical surfaeës of radii R1 and R 2 > Æ1, a wound of radius Sl = R101 may constitute a CSD, whereas a larger wound of radius s2 = R202 > sl may not, if 02 < 01. The criterion for wound healing to be initiated is instability of the equilibrium state ~ = 0, i.e., a > 0 or the reverse of inequality (14), but sinee we are primarily interested in the eonditions under whieh the existenee of a CSD is predieted in this model, the eonsequences of (14) as written are more useful. We state these in the form of a theorem. THEOREM.
(i) For a given value of O, and hence, wound radius, a CSD will exist provided
{ D ( l + l/2cot20) } 1/2, (b)
k2 > ~ .
ff Condition (b) is not satisfied, a > 0 and so healing occurs regardless of the values of R and 0. (ii) For a given value of R, and hence, "sku11" size, a CSD will exist provided the wound radius r = RO satisfies the criteria
(e)
a
r _> Rarecot {4R 2 @ 2 _ ~ )
_ 2} 1/2
(16)
and (d)
k2 > ~
1
a
+ ~.
(17)
If Criterion (d) is not satisfied, aga/n, as above, a > 0 and so healing occurs regardless of the values of R and Ô. The proofs follow from the dispersion relation (13) for a < O. Note that in each case, the subsidiary Conditions (b) and (d) require the perturbation wavenumbers to be sufficiently /arge, or equivalently, the/r wavelengths to be sufBciently smM1 for the existence of a CSD. A similar analysis about fl = y* shows that
af. , cr = - D k 2 - --~y so « < 0 and the "heMed" state y* is stable.
1090
J.A. ADAM
4.
BONE
CELL
AGGREGATION
It is known that when bone heals, as in the case of fractures, a coarsening or thickening of the bone frequently occurs in the vicinity of the original injury. Sometimes the resulting bone growth is referred to as a "collar', and may be irregular in its distribution in any given direction. We shall examine the aggregation or "clumpiness" of such new bone growth in the 8 direction, maintaining for simplicity the original assumption of azimuthal symmetry, though clearly this will be violated in general. In order to represent this aggregation mathematically within the existing model, the approach developed for an entirely different context by Britton [9,10] will be adopted (see also [11]). In those papers, the aggregation term was used to describe social advantages of grouping together (e.g., to reduce the risk of predation) and to incorporate effects on the specific growth rate of the population in the neighborhood of a point and not merely at the point of interest. Further ecological implications of this and related rnodels may be found in the above-mentioned references. Returning for the moment to the original equation (1), we focus attention on the logistic term which we write as aus(u), where s(u) = 1 - u / K is the specific growth rate for the logistic equation. Clearly, this depends pointwise on 8; the modification introduced in [9] is to extend this to include some type of weighted spatial average of the cell density near the point. Accordingly, in place of the above logistic term the expression s(u, v) is used, where
s(u, v)
1 + ~ K - (1 + ~)K'
(is)
being a positive parameter, and v being the spatial average defined in equation (19) below. The average is more or less local depending on the nature of the weighting; if v = u, then (18) reduces to the logistic case, the most localized form possible. At the other extreme, if u = K everywhere, then by definition v = K , and s(u, v) = 0; there is no net "reaction" term because resources are consumed and renewed at identical rates. In the present model, these resources are implicit, being the wound-induced growth factors that initiate and maintain bone cell proliferation until healing occurs. Subsequent work will focus on incorporating these resources explicitly. An implication of s being an increasing function of u is that there is a local advantage afforded by aggregation (the more cells in a neighborhood, the stronger the new bone is likely to be). On the other hand, the growth rate of cells might be expected to decrease as the average population increases, and this is reflected in s being a decreasing function of v. Note also (following [9]) from equation (18) that e ~ u / K is a short-range activation term while - ( 1 + o~)v/K is a long-range inhibition term, and this is a common feature of systems supporting pattern formation [12]. Following [9], we define v as a convolution u(8, t)
g * u(8, t) - f g(I 8 - ¢1)u(¢, t) d¢, Je
(19)
where the condition o g ( l e [ ) d¢ = 1
(20)
is appropriate for g to be considered an average. The population of new cells arising in the healing process is in the domain O = (0, 8o) C [0, 7r/2), 8o being the angle subtended by the wound radius prior to healing (if this occurs). The choice of g being a decreasing function of its argument reflects the likelihood that the more cells there are in a region around 8, the more competition there will be for the available resources (nutrients, growth factors, etc.) than if they are further away. A simple choice is, therefore, that of an exponentially decaying function of the form g(8) = 5e -te, (21) where 5=
e 1 -- e -cOo
The Effect of Surface Curvature
in accordance
with constraint
1091
(20). Thus,
v(e,t) =
s
s 80
e-+%((b,
t) dq5.
(22)
0
By expanding
this expression
as
s I3
v(e,t)
= 6
s 00
e-+@h(~,
t) df$ + b
0
e-+%(c#), t) &#I,
e
it follows that d%
s= and the reaction
term in equation
(1) is now.written
au
As before,
a change
in the dependent
(234
c2v - 2c6u as
,+Epyg*u]. [ variable
u is in order; using
(2) the modified
equation
(3 )
is (24) which,
together
with
(23b) forms the basis for the analysis below. A fuller study of the problem would require supplemental information concerning initial conditions on u (or y) and v, but that will not be necessary here. As above, spatially uniform solutions (for given 0) satisfy
2Sf
v=-y=ny
??J
and
C
afy
dt=K
(y* - m-ly)
)
(250)
where m-1
=
2[1 + Iy]S C
5. EXISTENCE
- cr.
OF A
CSD
As above, we will investigate stability to nonuniform perturbations about steady states, but now for a two-dimensional system (y,v) about (0,O) and (my*,nmy*). Note that if 26/c = 1, equation (25b) (with the angular derivative now included) reduces to (12), but for the present system, the definition of b precludes this limiting case, so the aggregation model is similar to but independent of the one-dimensional system. For the limiting case to exist, v E u and this occurs only when g(8,4) = &(0 - 4), the Dirac delta function. Linearizing
equations
(24) and (23b) about a+Dlc2-
(O,O), we obtain
a.f
_i7y*
Yl
=
the system
0,
>
-2cdfy1
+ (c” + kW)
W&b)
211= 0.
Nontrivial positive solutions exist if cr = -Dk2 + (af/K)y*, so (0,O) is linearly stable, and no growth occurs, i.e., a CSD exists, as in the single-variable problem, if inequality (14) holds as before.
1092
J. A. ADAM
A similar linearization about (my*, nmy*) yields the system o+Dk2-
y*2I1 = 0,
Wa,b)
Nontrivial solutions exist if aafm K ‘*
a+Dk2--
a(1 + o)m K ‘*
-2cbf
zz 0,
c2 + k2R2
i.e., if aafm 0 = a(k2) = - K ‘* -Dk2 EF-Dk2-
2cSfa( 1 + a)m -
K(C2
+k2R2)
y*
G
(284 (28b)
c2 + k2R2’
Note that a(O) = -afK-‘y* < 0, so the system is stable for k2 = 0, and also for k2 sufficiently large. There is a region of instability if a(k2) = 0 has two distinct positive roots k2 which are solutions of DR2k4 + (Dc2 - FR2) k2 + H = 0,
(29)
where H = G - c2F = c2afmy*Ke1. This will occur provided o > Q, = M2c2 + 2Mc,
(30)
where M2 =
DK = D{R2m R2af my*
[i
(I+
icot28)
+a]}-‘_
(31)
This condition assures the existence of a range of wavenumbers kl < k < kp for which the healing steady state (my*, nmy*) is unstable, i.e., aggregation (clumpiness) occurs. In terms of M2, the biquadratic equation (29) reduces to the tidier form M2R4k4 + (M2c2 - a) R2k2 + c2 = 0
(32)
and the stability diagram in the R2k2(a) plane is shown in the figure below for several values of M2 (0.6,1,2) with c2 = 1; similar types of graph are obtained for different values of c2. The regions of instability are to the right of each particular curve; when the tangent line to the curve is vertical, Q = a,. To interpret these results, notice that the determining quantity in inequality (30) is the parameter MC; c is a measure of how localized the average v of u is [9], so if c is large, only points close to 8 are strongly weighted, whereas if c is small, values of B further away become relatively more important. The parameter M is a measure of the geometric features of the wound and its environment; if in particular the skull radius R is sufficiently large or the wound “angle” 190is sufficiently small (or both), then M is small and this tends to counteract the consequences of c being large. Conversely, sufficiently small values of R and large values of #a would counteract the effects of c being small. In other words, aggregation will occur in this model if the product MC is sufficiently small.
The Effect of Surface Curvature
1093
3.5 3 2.5 21.510.5 0
2
4
6
8
Figure 3. The stability diagram in the R2k2(a) plane for M 2 = 0.6, 1 and 2 (increasing to the right) and c2 = 1; c~ = c~c at the vertex of each curve. The system is stable to aggregation on the left of each particular curve, and unstable on the right. 6.
DISCUSSION
B y describing the healing of a circular w o u n d on a spherical surface as t y p e of w a r e of advance, various criteria have been established: (i) the m i n i m u m wave speed as a function of skull radius and w o u n d angle, (ii) u p p e r b o u n d s on the time for healing of a given w o u n d to be complete (if the w o u n d does in fact heal), (iii) conditions on skull radius and w o u n d size for a CSD (critical size defect) to exist, and (iv) conditions under which cell "aggregation" occurs during the healing process. B u t w h a t does aggregation mean in the present context? T h e steady state (my*, n m y * ) , referred to above as the "healed" state is p r e s u m a b l y in general stable, i.e., c~ < sc, because one does not expect, in biological terms at least, such situations to be unstable. However, a possible (but speculative) interpretation is the following: the present model does not distinguish between w o u n d healing of bone or tissue on a spherical surface; the consequences of surface c u r v a t u r e are the m a j o r contributions of the model. However, it is known t h a t certain types of w o u n d healing in tissue can give rise to keloid scars. These are heavy scars in which much more scar tissue is p r o d u c e d t h a n is generally required for total healing, or more accurately, '% sharply e!evated , irregularly-shaped, progressively enlarging scar due to the formation of excessive a m o u n t s of collagen in the corium during connective tissue repair" [13] (the eorium is often called the dermis; it is "the layer of skin deep to the epidermis, consisting of a dense bed of vascular connective tissue" [13]). If the instability of the steady state can be identified with keloid scarring, t h e n the condition c~ < C~c would be expected to be the governing condition as rar as healing in bone is concerned. For ~ # 0, this will be the case for sufficiently large values of M e , i.e., for small e n o u g h values of R and large enough values of ~ or large enough values of e, or some suitable c o m b i n a t i o n of these quantities. If c~ > c~c, t h e n according to this interpretation, keloid scarring might be manifested in healed tissue. T h e present model does not include sufficient biological information to distinguish between tissue and bone; it would be of interest to a t t e m p t to explore such differences in future studies.
REFERENCES 1. J.A. Adam, The effect of surface curvature on wound healing in bone, Appl. Math. Lett. 15 (1), 59-62, (2002). 2. J.A. Adam, Healing times for circular wounds on plane and spherical bone surfaces, Appl. Math. Lett. 15 (1), 55-5S, (2002).
1094
J . A . ADAM
3. J.P. Schmitz and J.O. HoUinger, The critical size defect as an experimental model for craniomandibulofacial nonunions, Clinical Orthopaedics and Related Research 205, 299-308, (1986). 4. J.O. Hollinger and J.C. Kleinschmidt, The critical size defect as an experimental model to test bone repair materials, J. Craniofacial Surg. 1, 60-68, (1990). 5. J.D. Murray, Mathematical Biology, Chapter 11, Springer-Verlag, Heidelberg, (1989). 6. H.P. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theor. Biol. 56, 229-242, (1976). 7. M.A.J. Chaplain, The development of a spatial pattern in a model for cancer growth, In Experimental and Theoretical Advances in Biological Pattern Formation, (Edited by H.G. Othmer, P.K. Maini and J.D. Murray), pp. 45-60, Plenum Press, (1993). 8. J.A. Adam, General aspects of modeling tumor growth and immune response, In A Survey of Models for Tumor-Immune System Dynamics, (Edited by J.A. Adam and N. Bellomo), Chapter 2, Birkhaüser, Boston, MA, (1997). 9. N.F. Britton, Aggregation and the competitive exclusion principle, J. Theor. Biol. 136, 57-66, (1989). 10. N.F. Britton, Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model, Siam J. Appl. Math. 50, 1663-1688, (1990). 11. S.A. Gourley and N.F. Britton, A predator-prey reaction-diffusion system with nonlocal effects, J. Math. Biol. 34, 297-333, (1996). 12. A. Gierer and H. Meinhardt, A theory of biological pattern formation, Kybernetik 12, 20-39, (1972). 13. Dorland's Illustrated Medical Dictionary, (27 th Edition), W.B. Saunders Co., Philadelphia, PA, (1988).
PERGAMON
Mathematical and Computer Modelling 35 (2002) 1085-1094 www.elsevier.com/locate/mcm
T h e Effect of Surface C u r v a t u r e W o u n d H e a l i n g in B o n e : II. T h e Critical Size D e f e c t .
on
J . A . ADAM D e p a x t m e n t of M a t h e m a t i c s a n d S t a t i s t i c s O l d D o m i n i o n University, Norfolk, VA 23529, U.S.A.
(Received July 2001; accepted August 2001) Abstract--Two related models are studied for wound healing on a spherical surface. The first, and simpler model is based on a reaction-diffusion equation describing the process of healing as the inward growth of new cells from a circular wound on a spherical surface. Several aspects of the wound healing problem are considered: by interpreting the healing process as a "pseudowave" propagating across the spherical surface a heuristic account of the "speed" of healing is possible, and the corresponding healing time is characterized. Of paxticular importance in relation to animal models is the existeuce (or not) of a critical size defect (a CSD, defined below); this is discussed as a consequence of the stability of the steady states of the system to nonuniform spatial (or angular) perturbations. Explicit criteria are derived under which a CSD exists (within the model) in terms of the skull radius and wound radius. Although neither model explicitly distinguishes between wound healing in bone or tissue, it is in the former t h a t CSDs are known to occur. The second model invokes a weighted spatial average cell density which permits the presence of both a short-range activation t e r m (as in the first model) and a long-range inhibition term. Under these circumstances, within a suitable parameter range, the phenomenon of aggregation may occur in addition to the behavior predicted by the first model. It is suggested t h a t such aggregation is manifested in the case of keloid scarring, which can occur as a result of wound healing in tissue. (~ 2002 Elsevier Science Ltd. All rights reserved.
Keywords--Reaction-diffusion Keloid scar.
equation, Surface curvature, Critical size defect, Wound healing,
1. I N T R O D U C T I O N In an earlier paper [1], the time-independent distribution of wound-induced growth factor (GF) over a spherical cap was considered. This and related papers (see [2] and references therein) are associated with the existence of the so-called critical size defect, which has been defined as the smallest interosseous wound that does not heal by bone formation during the lifetime of the animal [3] (see also [2] for further clinical references). For practical purposes, this timescale can usually be taken as one year. In [4], the definition was further extended to a defect which has less than ten percent bony regeneration during the lifetime of the animal. CSDs can "heal" by fibrous c o n n e c t i v e t i s s u e f o r m a t i o n , b u t s i n c e t h i s is n o t b o n e , i t d o e s n o t h a v e t h e p r o p e r t i e s ( s t r e n g t h , etc.) that
a completely healed defect would.
Some typical CSD diameters
0895-7177/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S0895-7177(02)00073-0
a r e , for r a t , r a b b i t ,
Typeset by ¢4A/rS-TEX
1086
J.A. ADAM
dog, and monkey calvaria (skullcap), respectively: 8mm, 15 mm, 20mm, and 15 mm (details can be found in [3]). A phenomenological time-dependent problem is now examined with intent to determine, somewhat heuristically, the speed of wound healing, the time for wound healing to occur, and the existence of a critical size defect (CSD). The starting point is the logistic reaction-diffusion equation for the spread of new bone cells inward from a circular wound edge on the surface of a sphere of radius R. The distribution of new cells is considered to be azimuthally symmetric at all times. The initial wound angle subtended at the center of the spherical skull is 2~0 where 80 E (0, 7r/2); if total healing occurs, this angle decreases to zero (which corresponds to a polar angle of 7r/2 in standard notation, so 8 -- 0 is the "north" pole of the skull). Because of the azimuthal symmetry, the angular domain is 0 < 8 < 80 < zr/2. In terms of the bone cell density u(0, t), the governing partial differential equation is
Ou D ('O~u Ou) (i T Ot - - ~ \082 +cotOÕ-~ + a u _ u),
(1)
where D is the coefficient of diffusion for the GF, a is the linear growth rate, and K is the saturation limit for the bone cell density. The quantities D, a, and K are all constants in this model. However, it will be useful (in connection with the steady states ~ below) to allow these constants (and also 0) to be piecewise constant for the purposes of studying the evolution of the healing process. Eliminating the first derivative term via the change of dependent variable
u(8, t) = f(8)y(8, t) - (csc 0)1/2y(8, t),
(2)
equation (1) becomes
õ~=ZõB+~
l + ~ c o t 28 ~+«y
1---
,
(3)
where/~ = D / R 2. Homogeneous steady states ~ exist when ~ = 0 and
~{~(1 )}~.
=ä-]
2. P S E U D O W A V E
l + ~ c o t 28
SPEEDS;
+a
-
HEALING
(4) TIMES
By analogy with travelling wave solutions for infinite domains, we seek "pseudowave" solutions of the form y = y(~), ~ = - R S - ct. (5) This corresponds to the pseudowave moving in the direction of decreasing 8 with positive "speed" c. The qualification "pseudowave" is of course necessary because (i) the domain cannot be (-c~, oo), and also because of (ii) the presence of O-dependence in the new reactive term in equation (3). However, as noted by Murray [5] in connection with a mathematical model for calcium waves on amphibian eggs, some insight into mechanisms may still be obtained by using similar arguments for piecewise constant values of ~. The governing ordinary differential equation is
d2y dy a f , , D - ~ + C-d-~+ ---~y~y - y) = O.
(6)
A standard phase plane analysis for the (y(~), y~(~))-plane reveals that the nature of the equilibrium points (0, 0) and (y*, 0) is determined by the characteristic equation
a f . = 0, DA2 + cA =t: -~-y
(7)
T h e Effect of Surface C u r v a t u r e
1087
where the + and - signs refer to (0,0) and (y*,0), respectively. Thus, (y*,0) is a saddle point, and (0, 0) is a stable node provided that the pseudowave "speed"
{( ° ( ' D
Cw~amin:2
a+~-~-
1+~cot28
)))'~~
(s)
where ¢min is a constant for each fixed 8, where 0 < 8 < 80 < ~r/2 as noted above. This condition guarantees that the heteroclinic trajectory connecting these two steady states is the only bounded positive trajectory. Were the domain ( - c o , co), we would be completely justified in claiming that a "healing" wave moves from the stare y = y* to y = 0 as 8 decreases (and ~ increases). Thus, by analogy we consider the healed state to spread as a moving front throughout the (now) restricted domain. Note that Cmin increases for given R as 8 decreases, i.e., the wound heals laster the smaller is the angular extent of the wound. This also implies that if 8 is allowed to decrease in equation (8) in a piecewise constant manner, the healing process speeds up as the wound radius contracts. In addition, Cmin always exceeds its planar value 2x/-Dä for finite R, and increases away from this value as D increases. An upper bound for the time taken for a wound to heal (if indeed it does so) may be reasonably defined from this model as
Bth
L°
th :
dt-----
R
d8
o Cmin(8)
RSo = Th.
(9)
< Cmin(8)
The lower limit 80 is the initial (half-) angular size of the wound. Strictly, of course, this pseudowave analysis is based on c being constant, which it will be for a given value of 8; it would be of interest, therefore, to compare observed healing times in an experimental model with th and Th as calculated here. Based on studies of c in planar models, it appears likely that Cmin is the speed that will be attained for a given value of 8, though it is not clear whether the dependence on 8 will complicate the evolution of the pseudowave as it moves across the spherical cap. The integral can be evaluated to yield
th(R,8(0))
R
-
x/D(4a
L ° sin 8 + J) vv(0) X / # - 2 - - c o s 2 0
R
x/D( 4a + e)
{arcsin #
d8
- arcsin[# cos 80] },
(lOa) (10b)
where # =
\4äa ~ 2~,]
'
0 < # < 1.
(11)
For numerical purposes, the upper bound th can be expressed in dimensionless form as
4ath
x2
= ~(arcsin
and for comparison, the supreme
Th
#
arcsin[# cos 00])
(10c)
can be similarly expressed as
4aTh =
x2
~/x 2 + 2 + c o t 2 00
where x = 2Rv/-~~.
0o,
(10d)
1088
J.A.
ADAM
0.06 /
0.05 0.04
/
0.03• /
0.02.
./
/
II////
///
//
/ //i/
/
//
/
/
tl ¸
/
/
/
/
/
/
/
/
/
/'
/
~ " j"
/
0.01. tl ~
.I
_~'S~":
.lf" ........
. [ .........
..........-.-
0.2
,
__-.---.----
, i
0.4
0.8
0.6
Figure 1. The dimensionless upper bound on the healing time, 4ath(x), where x ---2Rv/äß, for fixed values of semiwound angle 00 (00 = 0.1, 0.3, and 0.5; these are the lower three curves in ascending order). The supremum Th for 00 ---- 0.5 is also shown for comparison (topmost curve). // 0.7" /
/
/
/
/
0.6
/ / /
/
0.5
/
/
/
0.4 0.3 /
0.2 0.1 0
//
//
_ y---
/
.j-
i
0.2
0.4
. . . .
0.6
I
,
~
0.8
I
1
Figure 2. The dimensionless upper bound on the healing time, 4ath(0O), for fixed values of x (x ---- 0.5, 1, 2, in ascending order; the supremum Th is not shown).
Graphs of 4ath(x) are presented in Figure 1 for fixed values of 00 (and the supreme Th in orte case is also shown for comparison). In Figure 2, graphs of 4ath(OO) are presented for fixed values of R (but Th is not shown). The excess Th -- th is greater for larger values of x (or R) and 0o than for smaller values, as would be expected for the case of larger wounds. In each case, the value of Th is about a factor of two larger than the corresponding value of th.
3. L I N E A R
STABILITY:
EXISTENCE
OF THE
CSD
We return to the governing differential equation (3) and explore the linear stability of the steady states ~ -- 0 and ~ -- y*. While it is customary to make modal expansions in terms of surface spherical harmonics (or in the case of azimuthal symmetry, as hefe, Legendre polynomials) [6-8], note that equation (3) is isomorphic to a one-dimensional diffusion equation, and a simple Fourier component of the form below will suffice for present purposes. The governing equation is, from
T h e Effect of Surface Curvature
1089
(3) and (4),
Oy
02y af , . = Zõõ~ + ~ y t y - y).
(12)
It is clear that spatially (or more accurately, angularly) uniform solutions satisfy a logistic equation, with the steady state ~ = 0 being linearly unstable, and ~ = y* being linearly stable to spatially uniform perturbations. Regarding nonuniform perturbations, let us consider stability to perturbations of wavenumber k by writing
Y = fl + eYle «t+iÆR° + 0 @2) . Substituting this expression into equation (12) for 9 = 0 and ignoring terms of O(e 2) and higher, the following dispersion relation arises, again valid for eaeh fixed value of 0 and R:
a=-Dk
«,
2 +-~y
=-Dk 2 +~
°(~
l + ~ c o t 20
)
+a,
(13)
so ~ = 0 is linearly stable, and no growth occurs, i.e., a CSD exists if
-~
> ~ l + ~ e o t 20 ,
(14)
i.e., for a sufficiently large wound radius r (= ROo) for a given R (assuming k 2 > a/D), or for a suffieiently large "skull" radius R for a given wound "angle" 00. One implication of this result, for a given value of k is that there may be partial healing only if there is a loeal deerease in the ratio a R 2 / D or (for constant a/D) an increase in the loeal radius of curvature in the wound environment. Note also that this result does uot imply that a larger wound radius will neeessarily result in a CSD as eompared with a smaller radius; given two spherical surfaeës of radii R1 and R 2 > Æ1, a wound of radius Sl = R101 may constitute a CSD, whereas a larger wound of radius s2 = R202 > sl may not, if 02 < 01. The criterion for wound healing to be initiated is instability of the equilibrium state ~ = 0, i.e., a > 0 or the reverse of inequality (14), but sinee we are primarily interested in the eonditions under whieh the existenee of a CSD is predieted in this model, the eonsequences of (14) as written are more useful. We state these in the form of a theorem. THEOREM.
(i) For a given value of O, and hence, wound radius, a CSD will exist provided
{ D ( l + l/2cot20) } 1/2, (b)
k2 > ~ .
ff Condition (b) is not satisfied, a > 0 and so healing occurs regardless of the values of R and 0. (ii) For a given value of R, and hence, "sku11" size, a CSD will exist provided the wound radius r = RO satisfies the criteria
(e)
a
r _> Rarecot {4R 2 @ 2 _ ~ )
_ 2} 1/2
(16)
and (d)
k2 > ~
1
a
+ ~.
(17)
If Criterion (d) is not satisfied, aga/n, as above, a > 0 and so healing occurs regardless of the values of R and Ô. The proofs follow from the dispersion relation (13) for a < O. Note that in each case, the subsidiary Conditions (b) and (d) require the perturbation wavenumbers to be sufficiently /arge, or equivalently, the/r wavelengths to be sufBciently smM1 for the existence of a CSD. A similar analysis about fl = y* shows that
af. , cr = - D k 2 - --~y so « < 0 and the "heMed" state y* is stable.
1090
J.A. ADAM
4.
BONE
CELL
AGGREGATION
It is known that when bone heals, as in the case of fractures, a coarsening or thickening of the bone frequently occurs in the vicinity of the original injury. Sometimes the resulting bone growth is referred to as a "collar', and may be irregular in its distribution in any given direction. We shall examine the aggregation or "clumpiness" of such new bone growth in the 8 direction, maintaining for simplicity the original assumption of azimuthal symmetry, though clearly this will be violated in general. In order to represent this aggregation mathematically within the existing model, the approach developed for an entirely different context by Britton [9,10] will be adopted (see also [11]). In those papers, the aggregation term was used to describe social advantages of grouping together (e.g., to reduce the risk of predation) and to incorporate effects on the specific growth rate of the population in the neighborhood of a point and not merely at the point of interest. Further ecological implications of this and related rnodels may be found in the above-mentioned references. Returning for the moment to the original equation (1), we focus attention on the logistic term which we write as aus(u), where s(u) = 1 - u / K is the specific growth rate for the logistic equation. Clearly, this depends pointwise on 8; the modification introduced in [9] is to extend this to include some type of weighted spatial average of the cell density near the point. Accordingly, in place of the above logistic term the expression s(u, v) is used, where
s(u, v)
1 + ~ K - (1 + ~)K'
(is)
being a positive parameter, and v being the spatial average defined in equation (19) below. The average is more or less local depending on the nature of the weighting; if v = u, then (18) reduces to the logistic case, the most localized form possible. At the other extreme, if u = K everywhere, then by definition v = K , and s(u, v) = 0; there is no net "reaction" term because resources are consumed and renewed at identical rates. In the present model, these resources are implicit, being the wound-induced growth factors that initiate and maintain bone cell proliferation until healing occurs. Subsequent work will focus on incorporating these resources explicitly. An implication of s being an increasing function of u is that there is a local advantage afforded by aggregation (the more cells in a neighborhood, the stronger the new bone is likely to be). On the other hand, the growth rate of cells might be expected to decrease as the average population increases, and this is reflected in s being a decreasing function of v. Note also (following [9]) from equation (18) that e ~ u / K is a short-range activation term while - ( 1 + o~)v/K is a long-range inhibition term, and this is a common feature of systems supporting pattern formation [12]. Following [9], we define v as a convolution u(8, t)
g * u(8, t) - f g(I 8 - ¢1)u(¢, t) d¢, Je
(19)
where the condition o g ( l e [ ) d¢ = 1
(20)
is appropriate for g to be considered an average. The population of new cells arising in the healing process is in the domain O = (0, 8o) C [0, 7r/2), 8o being the angle subtended by the wound radius prior to healing (if this occurs). The choice of g being a decreasing function of its argument reflects the likelihood that the more cells there are in a region around 8, the more competition there will be for the available resources (nutrients, growth factors, etc.) than if they are further away. A simple choice is, therefore, that of an exponentially decaying function of the form g(8) = 5e -te, (21) where 5=
e 1 -- e -cOo
The Effect of Surface Curvature
in accordance
with constraint
1091
(20). Thus,
v(e,t) =
s
s 80
e-+%((b,
t) dq5.
(22)
0
By expanding
this expression
as
s I3
v(e,t)
= 6
s 00
e-+@h(~,
t) df$ + b
0
e-+%(c#), t) &#I,
e
it follows that d%
s= and the reaction
term in equation
(1) is now.written
au
As before,
a change
in the dependent
(234
c2v - 2c6u as
,+Epyg*u]. [ variable
u is in order; using
(2) the modified
equation
(3 )
is (24) which,
together
with
(23b) forms the basis for the analysis below. A fuller study of the problem would require supplemental information concerning initial conditions on u (or y) and v, but that will not be necessary here. As above, spatially uniform solutions (for given 0) satisfy
2Sf
v=-y=ny
??J
and
C
afy
dt=K
(y* - m-ly)
)
(250)
where m-1
=
2[1 + Iy]S C
5. EXISTENCE
- cr.
OF A
CSD
As above, we will investigate stability to nonuniform perturbations about steady states, but now for a two-dimensional system (y,v) about (0,O) and (my*,nmy*). Note that if 26/c = 1, equation (25b) (with the angular derivative now included) reduces to (12), but for the present system, the definition of b precludes this limiting case, so the aggregation model is similar to but independent of the one-dimensional system. For the limiting case to exist, v E u and this occurs only when g(8,4) = &(0 - 4), the Dirac delta function. Linearizing
equations
(24) and (23b) about a+Dlc2-
(O,O), we obtain
a.f
_i7y*
Yl
=
the system
0,
>
-2cdfy1
+ (c” + kW)
W&b)
211= 0.
Nontrivial positive solutions exist if cr = -Dk2 + (af/K)y*, so (0,O) is linearly stable, and no growth occurs, i.e., a CSD exists, as in the single-variable problem, if inequality (14) holds as before.
1092
J. A. ADAM
A similar linearization about (my*, nmy*) yields the system o+Dk2-
y*2I1 = 0,
Wa,b)
Nontrivial solutions exist if aafm K ‘*
a+Dk2--
a(1 + o)m K ‘*
-2cbf
zz 0,
c2 + k2R2
i.e., if aafm 0 = a(k2) = - K ‘* -Dk2 EF-Dk2-
2cSfa( 1 + a)m -
K(C2
+k2R2)
y*
G
(284 (28b)
c2 + k2R2’
Note that a(O) = -afK-‘y* < 0, so the system is stable for k2 = 0, and also for k2 sufficiently large. There is a region of instability if a(k2) = 0 has two distinct positive roots k2 which are solutions of DR2k4 + (Dc2 - FR2) k2 + H = 0,
(29)
where H = G - c2F = c2afmy*Ke1. This will occur provided o > Q, = M2c2 + 2Mc,
(30)
where M2 =
DK = D{R2m R2af my*
[i
(I+
icot28)
+a]}-‘_
(31)
This condition assures the existence of a range of wavenumbers kl < k < kp for which the healing steady state (my*, nmy*) is unstable, i.e., aggregation (clumpiness) occurs. In terms of M2, the biquadratic equation (29) reduces to the tidier form M2R4k4 + (M2c2 - a) R2k2 + c2 = 0
(32)
and the stability diagram in the R2k2(a) plane is shown in the figure below for several values of M2 (0.6,1,2) with c2 = 1; similar types of graph are obtained for different values of c2. The regions of instability are to the right of each particular curve; when the tangent line to the curve is vertical, Q = a,. To interpret these results, notice that the determining quantity in inequality (30) is the parameter MC; c is a measure of how localized the average v of u is [9], so if c is large, only points close to 8 are strongly weighted, whereas if c is small, values of B further away become relatively more important. The parameter M is a measure of the geometric features of the wound and its environment; if in particular the skull radius R is sufficiently large or the wound “angle” 190is sufficiently small (or both), then M is small and this tends to counteract the consequences of c being large. Conversely, sufficiently small values of R and large values of #a would counteract the effects of c being small. In other words, aggregation will occur in this model if the product MC is sufficiently small.
The Effect of Surface Curvature
1093
3.5 3 2.5 21.510.5 0
2
4
6
8
Figure 3. The stability diagram in the R2k2(a) plane for M 2 = 0.6, 1 and 2 (increasing to the right) and c2 = 1; c~ = c~c at the vertex of each curve. The system is stable to aggregation on the left of each particular curve, and unstable on the right. 6.
DISCUSSION
B y describing the healing of a circular w o u n d on a spherical surface as t y p e of w a r e of advance, various criteria have been established: (i) the m i n i m u m wave speed as a function of skull radius and w o u n d angle, (ii) u p p e r b o u n d s on the time for healing of a given w o u n d to be complete (if the w o u n d does in fact heal), (iii) conditions on skull radius and w o u n d size for a CSD (critical size defect) to exist, and (iv) conditions under which cell "aggregation" occurs during the healing process. B u t w h a t does aggregation mean in the present context? T h e steady state (my*, n m y * ) , referred to above as the "healed" state is p r e s u m a b l y in general stable, i.e., c~ < sc, because one does not expect, in biological terms at least, such situations to be unstable. However, a possible (but speculative) interpretation is the following: the present model does not distinguish between w o u n d healing of bone or tissue on a spherical surface; the consequences of surface c u r v a t u r e are the m a j o r contributions of the model. However, it is known t h a t certain types of w o u n d healing in tissue can give rise to keloid scars. These are heavy scars in which much more scar tissue is p r o d u c e d t h a n is generally required for total healing, or more accurately, '% sharply e!evated , irregularly-shaped, progressively enlarging scar due to the formation of excessive a m o u n t s of collagen in the corium during connective tissue repair" [13] (the eorium is often called the dermis; it is "the layer of skin deep to the epidermis, consisting of a dense bed of vascular connective tissue" [13]). If the instability of the steady state can be identified with keloid scarring, t h e n the condition c~ < C~c would be expected to be the governing condition as rar as healing in bone is concerned. For ~ # 0, this will be the case for sufficiently large values of M e , i.e., for small e n o u g h values of R and large enough values of ~ or large enough values of e, or some suitable c o m b i n a t i o n of these quantities. If c~ > c~c, t h e n according to this interpretation, keloid scarring might be manifested in healed tissue. T h e present model does not include sufficient biological information to distinguish between tissue and bone; it would be of interest to a t t e m p t to explore such differences in future studies.
REFERENCES 1. J.A. Adam, The effect of surface curvature on wound healing in bone, Appl. Math. Lett. 15 (1), 59-62, (2002). 2. J.A. Adam, Healing times for circular wounds on plane and spherical bone surfaces, Appl. Math. Lett. 15 (1), 55-5S, (2002).
1094
J . A . ADAM
3. J.P. Schmitz and J.O. HoUinger, The critical size defect as an experimental model for craniomandibulofacial nonunions, Clinical Orthopaedics and Related Research 205, 299-308, (1986). 4. J.O. Hollinger and J.C. Kleinschmidt, The critical size defect as an experimental model to test bone repair materials, J. Craniofacial Surg. 1, 60-68, (1990). 5. J.D. Murray, Mathematical Biology, Chapter 11, Springer-Verlag, Heidelberg, (1989). 6. H.P. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theor. Biol. 56, 229-242, (1976). 7. M.A.J. Chaplain, The development of a spatial pattern in a model for cancer growth, In Experimental and Theoretical Advances in Biological Pattern Formation, (Edited by H.G. Othmer, P.K. Maini and J.D. Murray), pp. 45-60, Plenum Press, (1993). 8. J.A. Adam, General aspects of modeling tumor growth and immune response, In A Survey of Models for Tumor-Immune System Dynamics, (Edited by J.A. Adam and N. Bellomo), Chapter 2, Birkhaüser, Boston, MA, (1997). 9. N.F. Britton, Aggregation and the competitive exclusion principle, J. Theor. Biol. 136, 57-66, (1989). 10. N.F. Britton, Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model, Siam J. Appl. Math. 50, 1663-1688, (1990). 11. S.A. Gourley and N.F. Britton, A predator-prey reaction-diffusion system with nonlocal effects, J. Math. Biol. 34, 297-333, (1996). 12. A. Gierer and H. Meinhardt, A theory of biological pattern formation, Kybernetik 12, 20-39, (1972). 13. Dorland's Illustrated Medical Dictionary, (27 th Edition), W.B. Saunders Co., Philadelphia, PA, (1988).