J. theor. Biol. (1992) 158, 135-172

Continuum Model of Fibroblast-driven Wound Contraction: Inflammation-mediation ROBERT T. TRANQUILLOt§ AND J. D. MURRAY~

t Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, M N 55455 and ~ Applied Mathematics FS-20, University of Washington, Seattle, WA 98195, U.S.A. (Received on 26 February 1991, Accepted on 21 January 1992) We propose a mathematical model to aid the understanding of how events in wound healing are orchestrated to result in wound contraction. Ultimately, a validated model could provide a predictive means for enhancing or mitigating contraction as is appropriate for managing a particular wound. The complex nature of wound healing and the lack of a modeling framework which can account for both the relevant cell biology and biomechanics are major reasons for the absence of models to date. Here we adapt a model originally proposed by Murray and co-workers to show how cell traction forces can result in spatial patterns of cell aggregates since it offers a framework for understanding how traction exerted by wound fibroblasts drives wound contraction. Since it is a continuum model based on conservation laws which reflect assumed cell and tissue properties, it is readily extended to account for emerging understanding of the cell biology of wound healing and its relationship to inflammation. We consider various sets of assumed properties, based on current knowledge, within a base model of dermal wound healing and compare predictions o f the rate and extent of wound contraction to published experimental results.

I. Introduction Wound repair is a basic and vital process for homeostasis. It is also complex: in fullthickness cutaneous wounds, where loss of dermis has occurred, there are characteristic phases of immediate inflammation, short-term granulation tissue formation (simultaneous fibroplasia and angiogenesis), and long-term remodeling of the wound extracellular matrix (ECM) (Clark, 1985, 1988). Wound contraction, the centripetal movement of the wound boundary and adjacent uninjured skin towards the wound center, is another hallmark in the healing of these dermal wounds (Montandon et al., 1977). It often occurs during fibroplasia, the population of the fibrin clot by fibroblastic cells, namely, activated dermal fibroblasts and myofibroblastst. § Author to whom correspondence should be addressed. t Myofibroblasts are mesenchymal cells, perhaps deriving from dedifferentiated fibroblasts in the adjacent dermis, which share the ECM biosynthetic properties of dermal fibroblasts and the extensive actin filament motility system of smooth muscle cells. It is because of the latter property and their common observation in granulation tissue that they were originally proposed to be the cell type underlying wound contraction [Gabbiani et al. (1971); references to the original papers for all of the relevant findings stated here can be found in the review by Skalli & Gabbiani (1988)]. Subsequently, the myofibroblast concentration was reported to be relatively uniform throughout the wound granulation tissue during the time-course of wound contraction and positively correlated with the rate of contraction (McGrath & Hundahl, 1982). 135

0022-5193/92/I 80135 + 38 $08.00/0

© 1992 Academic Press Limited

136

R.T.

TRANQUILLO

AND

J. D . M U R R A Y

Compelling evidence suggests that fibroblastic cells, which secrete the collagen needed to transform the weak fibrin clot to a collagenous scar with significant mechanical strength, are also responsible for the force which drives wound contraction. However, the signals which determine the spatial and temporal distribution of fibroblastic cells in and around a wound, as well as the factors which govern the magnitude of force transmitted by the cells to local ECM fibers are poorly understood. These issues underlie the different hypotheses about where the locus of cell force exists during wound contraction (Rudolph, 1979), for example, the "picture frame" hypothesis, where fibroblastic cells pull the wound boundary in as they migrate in from the wound edge, and the "pull" hypothesis, where fibroblastic cells already within the wound generate a tension which causes the inward movement of the wound boundary. Even the nature of the exerted cell force has been debated (Ehrlich et aL, 1986; Ehrlich, 1988): "contraction force" associated with the hypothesized shortening of the entire cell body of a smooth muscle cell-like myofibroblast (Gabbiani et aL, 1979) vs. "traction force" associated with undulating pseudopods transiently adhering to or anchoring around collagen fibers, documented for fibroblasts cultured in reconstituted collagen gels (Harris et al., 1981; Harris, 1982, 1984; Stopak & Harris, 1982). The motility mechanism by which fibroblasts and/or myofibroblasts cause wound contraction is still an open question (the modeling approach described here is sufficiently abstract so as to accommodate both mechanisms, although we shall refer to the operative cell force solely as traction force, reflecting our bias, and use "fibroblast" to mean any fibroblastic cell exerting traction). Nonetheless, the traction which presumably allows a fibroblast to migrate from the surrounding dermis into the developing granulation tissue, at the same time causes some local reorganization of ECM fibers. This fiber reorganization, in turn, causes a local deformation of the ECM. The macroscopic manifestation of all the local deformations due to the traction forces exerted by all the fibroblasts involved in wound healing is wound contraction. Of course, other forces, such as viscous, elastic, and osmotic forces, also contribute to the actual rate and extent of contraction. There is strong motivation for understanding the interplay between the complex cellular, biochemical, and biomechanical phenomena which results in wound contraction, since it can be a beneficial or deleterious consequence of wound healing (Rudolph, 1980). It can be beneficial by enhancing the rate of wound closure beyond that associated with epithelialization. Further, part of the wound is then covered by functional dermis in contrast to the entire wound becoming a dense collagenous scar covered by a thin epithelium, the product of the wound healing process in the absence of contraction. However, the disfigurement caused by the adjacent skin being pulled into tension as a result of contraction (termed a contracture) can be cosmetically unacceptable and even render a proximal joint non-functional. If a predictive model of wound healing were developed, that accounted for the salient phenomena and provided a framework to accommodate the emerging understanding of the relationship between wound healing and inflammation, then a rational basis for the design and application of pharmacologic modulators of cell behavior to augment or mitigate wound contraction (as appropriate for the nature of the wound being managed) would be available. The implication of understanding the interplay between

CONTINUUM MODEL OF WOUND CONTRACTION

137

these phenomena goes far beyond improved and novel therapies for dermal wounds, however (Seemayer et al., 1981; Skalli & Gabbiani, 1988). Mechanisms related to fibroblast-driven wound contraction are believed to be operative in other pathologies involving inflammatory and wound healing components (e.g. liver cirrhosis, burn contracture, kidney fibrosis, granuloma, fibrous encapsulaton of implants, hepatic fibrosis, and pulmonary fibrosis), in fibramotoses, where a prominent inflammatory component is not evident (e.g. Dupuytren's contracture and dermatofibroma), and in the stromal response to neoplasia, where solid tumor growth is known to involve induction of a wound healing response (Dvorak, 1986). The same mechanisms are also important in the structural transformation of some artificial skin substitutes which use a dermal layer comprised of a type I collagen-based gel in which cultured fibroblasts are incorporated (Bell et aL, 1983). As indicated above, wound contraction occurs in the context of other phases of wound healing which are influentia'l in the initiation and regulation of contraction (Clark, 1985, 1988). After a fibrin clot is formed in the void from a dermal wound, there are ensuing phases of acute and immune inflammation. Inflammatory cells, perhaps most significantly activated macrophages, are believed to release a host of bioactive factors (i.e. chemotactic and growth factors) resulting in the infiltration of fibroblasts (and possibly their phenotypic transformation into myofibroblasts) and endothelial cells, thus leading to the subsequent phases of repair: granulation tissue formation characterized by the secretion of ECM macromolecules by fibroblasts and vascularization by endothelial cells, followed by an extended period of wound ECM remodeling ultimately resulting in a collagenous scar. Contraction usually parallels the phase of granulation tissue formation when the initial fibrin clot is being transformed into a network of collagen fibers and other ECM macromolecules. The close temporal relationships of contraction with inflammation and ECM biosynthesis will provide the basis for scenarios to be considered within our mathematical model of wound contraction. Wounds are not amenable to detailed quantitative study, and almost all such studies conducted to date have been performed in animal models rather than human subjects. This is somewhat problematic since most animals possess the panniculus carnosus, a skin organ beneath the dermis which is absent in humans and confers great mobility to skin in response to stress. Indeed, whereas up to 100% of animal wound area may be closed by contraction, this value is considerably less in humans, for example, 20-40% of post-excisional wound areas of the forearm (Catty, 1965; Ramirez et aL, 1969) (the balance of wound closure is attributable to epithelialization over the granulation tissue). Notwithstanding this anatomical difference, it is generally believed that the underlying mechanisms are similar, and there are several quantitative studies based on animal wounds addressing the longstanding controversy of what type of wound closes the fastest. This is typically accomplished by measuring wound area following the creation of an excisional wound, using tattoo marks to delineate the wound boundary and distinguish contraction from epithelialization. McGrath & Simon (1983) definitively showed that the contracting phase of excised dermal wounds on rats could be described with a simple exponential

138

R, T. T R A N Q U I L L O

AND

J. D . M U R R A Y

dependence on time: A ( t) = AS+ ( Ao - A I) exp ( - k d )

(1)

where A0 is the wound area when contraction begins, Af is the area remaining after contraction is completed, both areas being scaled to the excised area, and kc is the contraction rate constant. Excisional wounds undergo an immediate expansion, exhibit a "lag" phase of several days, then exhibit the exponential contracting phase described by (I), which occurs on the order of weeks (see Fig. l). AI and kc were found to be essentially independent of both the excised wound area and wound geometry for the limited cases studied (6.25cm 2 and 12.54cm 2 squares, and 12-54 cm 2 circles), as seen in Table I. The major difference in the results from these two wound geometries is associated with Ao, with the immediate expansion of square wounds being considerably greater than circular wounds of the same excised area, leading to a larger A0. Wound contraction in humans is also characterized by the lag and exponential phases (Ramirez et al., 1969) although there is no data to prove the same conclusions regarding the independence of A I and kc on initial wound area and geometry. These characteristic phases provide a basis for validation of any mathematical model of the dynamics of wound contraction. Of course, migration, division, and biosynthetic function of the fibroblast is dictated by a variety of biochemical and biophysical cues in the environment. A comprehensive review of all the known and speculated cues and the responses they

I1"2t o

1.0t0.8 0.6

o

tJ tY

0-4 0.2 ] I0

0 i

I 20

I 30

I 40

50

II Plateau phase

Exponential phase Time (days)

FIG. I. A plot o f the wound area relative to the excised area for a full thickness rat wound (McGrath & Simon, 1983: Fig. I). A0 is the wound area when contraction begins, A[ is the area remaining after eontracti6n is completed, and A~=Ao-Ar is the contracted area. See text for details.

CONTINUUM MODEL OF WOUND CONTRACTION

139

TABLE 1 Contraction data o f McGrath & Simon (1983)for full-thickness wounds on rats Wound Large circle (12 5 4 c m 2) Large square:~ (12-54 cm 2) Small square§ (6.25 cm 2)

k~

Akc

Ao *

AI*

Air

Act

&Ait

&Act

0.234

--

1.21

0.214

0.176

0.824

--

--

0-207

-12%

1.47

0-18t

0-123

0-877

-30%

+6-5%

0.245

+18%

1.54

0.207

0" 134

0.866

+8.9%

-I.2%

* These values are relative to the pre-excision inscribed area. t These values are relative to the area post-excision after the immediate expansion. A values are relative to large circle wound. § A values are relative to large square wound.

elicit is far beyond the scope of this paper; however, a few will be listed to make a point. As mentioned above, evidence suggests that the dedifferentiation of fibroblasts [or progenitor cell(s)] into myofibroblasts is initially regulated by a diffusible factor released by inflammatory cells in the wound, most probably one or more growth factors released by macrophages (Riches, 1988). As contraction ensues, it is possible that stresses developing in the dermis around the wound also regulate the phenotype transformation. They may also modulate metabolism and proliferation as has been documented for fibroblasts cultured in reconstituted collagen gels subjected to tensile (Jain et al., 1990) and compressive (Mochitate et al., 1991) stress. The migration of fibroblasts into the wound may involve one or all of the following: chemokinesis, chemotaxis, haptotaxis, and contact guidance (McCarthy et al., 1988). There are in vitro studies which suggest or demonstrate these behaviors and provide the identities of the diffusible or immobilized factors eliciting the kinesis and taxes. Again, macrophages are believed to be primary sources of the factors in vivo, either secreted directly (i.e. cell-derived chemotactic factors) or via the action of their proteolytic function (i.e. ECM-derived chemotactic fragments). In an interesting positive feedback mechanism, the tissue contraction resulting from the traction exerted by fibroblasts in the wound can induce alignment of ECM fibers, which could provide a contact guidance cue leading to further accumulation of migrating fibroblasts in the wound (see Stopak & Harris, 1982, for an in vitro system which mimics this). Angiogenesis, occurring in parallel wth fibroplasia, causes granulation tissue to become vascularized, thus influencing the concentrations of metabolic nutrients, perhaps most importantly dissolved oxygen, available to the fibroblasts. It is known that fibroblast proliferation and collagen synthesis vary along the oxygen tension gradient in granulation tissue (Niinikoski, 1977; Bankey et al., 1989). Finally, there is the complex regulation of gene expression and consequent function of the fibroblast, embodied in the "dynamic reciprocity" principle of Bissell (Bissell & Barcellos-Hoff, 1987), as the composition of the wound is dramatically transformed during debridement and fibroplasia and the state of stress evolves with contraction. The elusive but critical point for the decrease in fibroblast numbers in the wound (including the apparent differentiation of myofibroblasts back to the progenitor cells), which is

140

R . T . TRANQUILLO AND J. D. MURRAY

correlated with a decrease in the rate of wound contraction (Rudolph, 1980; McGrath & Hundahl, 1982) may reside at this level. The collective complexity of these phenomena, given the known and as yet unknown coupling between them, is daunting. This motivates the development of a mathematical model which can be used to help elucidate the complex couplings and ultimately to exploit knowledge of them. It is appropriate that the model be abstract in the sense of omitting uncertain detail as well as flexible in the sense of being able to accommodate fundamentally different cell biology when postulated to be of significance. The mathematical model used in this work, based on the continuum model proposed by G. F. Oster, J. D. Murray, and their colleagues (Murray et al., 1983; Oster et aL, 1983; Murray & Oster, 1984; Murray et al., 1988; see Murray, 1989 for a pedagogical survey), possesses both of these features. It accounts for known cell functions, such as migration and division, and the coupling of cell traction forces to the mechanical state of the ECM, resulting in ECM deformation. In this paper, we examine predictions of contraction for an idealized dermal wound using extensions of the Oster-Murray model which incorporate salient features of wound healing suggested above and motivated further below. Our main objective here is to provide a framework for understanding how known and speculated fibroblast functions are orchestrated with other cellular, biochemical, and biomechanical events in wound healing so as to give rise to contraction. Our primary criterion of whether a particular set of fibroblast functions, ECM properties, and wound healing events is sufficient is qualitative consistency with the contraction curve typified by Fig. 1. The number of sets which could be considered is quite large given the complexity of the wound healing response, so our survey cannot be comprehensive. Our strategy has been to examine a "base model" consisting of a set on which there would seem to be general agreement (although not necessarily agreement on how they are modeled). Given that the base model does not yield a contraction curve qualitatively consistent with Fig. 1, extensions of the base model are made by incorporating additional fibroblast functions, ECM properties, and wound healing events and then checking again for consistency with Fig. 1. In addition to identifying those sets which meet that criterion, the predictive use of the model is demonstrated by suggesting how modification of fibroblast functions and ECM properties enhances or diminishes the rate and extent of contraction, whichever is desired for managing a particular dermal wound. A preliminary report of this work is found in Murray et al. (1988).

2. Base Model

q'he general features of the Oster-Murray model from which the base model for dermal wound contraction is formulated have been described at length elsewhere (e.g. Murray, 1989). Only a brief summary is given here with an emphasis on considerations peculiar to wound healing and contraction. General modeling and analysis issues which apply to all model extensions are also discussed in this section. Since we are only concerned with wound contraction here, we do not consider the subsequent phases of wound healing, and assuming wound contraction is not significantly

CONTINUUM MODEL OF WOUND CONTRACTION

141

influenced by epithelialization, we neglect the latter as well (see Sberratt & Murray, 1990, 1991). 2.1. MODEL EQUATIONS In the base model and all of its extensions, there are only three variables in the formulation: the fibroblast concentration n(_x, t), the ECM concentration p(x_, t), and the displacement vector for the fibroblast/ECM composite _u(_x,t). The current formulation reflects several significant simplifications which are consistent with the chosen degree of abstraction for the model at present: dynamics associated with migration and proliferation offibroblasts alone will be described explicitly (neglecting other cell types at this point); dynamics associated with movement and synthesis of the collective ECM will be described (rather than for its constitutive components); macroscopic stresses developed in the ECM associated with fibroblasts distributed throughout exerting traction will be described ("coarse-graining" of the point forces exerted by individual cells on individual ECM fibers). We briefly discuss these simplifications next. The first simplification may be the most severe. As mentioned above, fibroblasts appear to transform between phenotypes during the course of wound healing which may vary with respect to their traction capacity (among other properties). In this base model, however, we consider only a single population of fibroblasts with all cells exerting the same magnitude of traction. Inflammatory and immune cells and endothelial cells also exhibit transient appearances in the granulation tissue and adjacent dermis. In several extensions of the base model, an implicit account is made for the speculated role of wound macrophages in releasing soluble bioactive mediators which modulate fibroblast behavior (traction, division, and migration) by postulating that a concentration gradient of such a mediator exists during the course of wound contraction. However, no account is made of the other major cell type believed to significantly influence the wound healing process, the endothelial cell, which through the process of angiogenesis causes granulation tissue to be vascularized, thus influencing the concentrations of metabolic nutrients for fibroblasts. In a study of dermal wound healing and contraction in which a specific inhibitor of angiogenesis was employed (McGrath & Emery, 1985), there was no influence of the diminished vascularization relative to controls on the duration of the lag phase and the wounds still exhibited an exponential phase, but kc decreased in a dose-dependent fashion. Based on ultrastructural analyses of the wound fibroblasts and granulation tissue, this was believed to occur because of reduced collagen synthesis of fibroblasts in the more hypoxic granulation tissue rendering it collagen-deficient with disorganized fiber network structure, rather than because of a reduced ability of fibroblasts to exert traction on the compromised network. Instead of attempting to include angiogenesis in this first generation model, we assume that wound contraction can be described to a first approximation without its inclusion, i.e. angiogenesis is assumed to proceed such that there are no limitations of oxygen and other bloodderived nutrients which would modulate the fibroblast responses. Note that any simplifying assumption can be appropriately modified where reliable modeling can be made.

142

R.T.

TRANQUILLO

AND

J. D .

MURRAY

The second simplification means that the variable p represents the ECM concentration in both the wound and the surrounding dermis (i.e. irrespective of composition and, therefore, of location). In reality, extensive modification of the wound ECM composition occurs during wound healing, the initial fibrin clot transforming into dynamic granulation tissue with time-varying percentages of collagen, fibronectin, hyaluronic acid, and other proteoglycans and glycosaminoglycans (Clark, 1985; Weigel et al., 1986), all of which are synthesized by infiltrating cells and in turn likely to influence their behavior in a concentration-dependent manner. Consistent with the abstraction embodied in the definition of p is the assumption that the rheological properties of the wound ECM and surrounding dermis are the same. This is also a gross assumption but appropriate, given that knowledge of granulation tissue theology is virtually non-existent. Also consistent is the notion that "direct synthesis" of ECM is meaningful in the context of the model, that is, no explicit account is made for the complex transport and assembly steps of collagen and other ECM molecules. The final simplification which takes account of cell traction forces as macroscopic "active stress" is the crux of the Oster-Murray model. Cell cortical cytoplasm, an analogous "contractile fiber network", has been modeled by Dembo & Harlow (1986). They account for the forces associated with actomyosin complexes distributed throughout the submembraneous actin-rich cortex of motile cells as a macroscopic contractile stress. In contrast to our approach, they model the cytoplasm as distinct fluid and network (actin filament) phases, leading to significantly more mathematical complexity. At this stage we retain our simpler concept for the biomechanics given the other phenomena to be modeled as well. Returning to the model formulation, the three variables satisfy fundamental conservation laws: species conservation, yielding equations describing n and p, and linear momentum conservation, yielding the equation describing u. The specific terms of these equations depend, of course, on both the phenomena presumed operative as well as their mathematical modeling. The species conservation equations express the requirement that the rate of accumulation of a species in a control volume equals the net flux of the species through the boundaries of the control volume plus the net rate of generation within the control volume (e.g. Segel, 1984). For the fibroblast conservation equation, this translates into terms for active migration, division, and death. In the base model we assume only that fibroblasts migrate randomly and exhibit concentration-regulated proliferation characterized by some maximum fibroblast concentration. Thus, the speculation that gradients of mediators cause fibroblasts to accumulate in the wound by chemotaxis, haptotaxis, contact guidance, or enhanced proliferation, or that fibroblasts exhibit activation or phenotypic transformation into a traction-enhanced state is not included here. These additional properties, for which there is some suggestive if not definitive experimental verification in vitro, are considered here in extensions of this base model, except for the cases of contact guidance and haptotaxis. An account for these, which can be elicited by cues generated by the fibroblasts themselves as they exert traction and drive wound contraction, causing ECM fibers to become aligned (contact guidance) and concentrated towards the wound (haptotaxis), is made in Tranquillo et al. 1992.

CONTINUUM

MODEL

OF

WOUND

CONTRACTION

143

We take the equivalent of Fick's Law for modeling random migration, as demonstrated for fibroblasts on planar substrata (Dunn, 1983) and other cell types in isotropic collagen gels (Noble & Boyarsky, 1988), assuming sufficiently low cell concentrations such that the cells migrate independently. A reasonable choice for modeling the net rate of fibroblast proliferation is the logistic growth rate ("growth" meaning cell division), which is consistent with data for cultured dermal fibroblasts in collagen gels (Schor, 1980) and myofibroblasts on substrata (Vande Berg et al., 1989). With these assumptions, the conservation law for fibroblasts describing the time rate of change of the local fibroblast concentration is: __On+ V . n 0-U=v . DoVn + ton(No - n). Ot Ot

(2)

The cell diffusion coefficient, Do ("cell" will mean fibroblast hereafter), growth rate constant, r0, and maximum cell concentration, No, are basal values in the absence of any external fields (e.g. growth factor concentration gradient). The convective term is a passive contribution to the net cell flux due to convection of the ECM in which the cells reside with velocity Ou_/Ot, i.e. a time-dependent displacement of the celI/ECM composite. For the ECM conservation equation, terms reflecting transport and net generation of ECM must be specified. Since the ECM, a composite comprised of a fiber network embedded in ground substance, is modeled as a homogeneous compressible solid, there is only a convective component to its flux. Fibroblasts play a key, but incompletely understood, role in the regulated levels of collagen, fibronectin, and other constituents of the dermis and granulation tissue. However, to keep the base model as simple as possible, no generation term is included (see "ECM biosynthesis" extension below). With these assumptions, the ECM conservation equation describing the time rate of change of the local ECM concentration is simply: 0£ + V at

p

0u

= o.

(3)

The momentum conservation equation expresses the requirement that the rate of accumulation of momentum in the control volume equals the net momentum flux through the boundaries of the control volume associated with surface stresses plus the net rate of momentum generation within the control volume due to body forces. Inertial forces are negligible for celI/ECM mechanics (Odell et al., 1981), simplifying the momentum conservation equation to a mechanical force balance between traction forces exerted by the cells and forces associated with the physico-chemical properties of the ECM. The constitutive stress-strain relation for dermis, if one was available, would be quite complex (Fung, 1981), and this must certainly also be true for granulation tissue. Consistent with our modeling philosophy but with the loss of considerable generality, we describe the ceI1/ECM composite with an isotropic, linear viscoelastic stress tensor in the base model and in all extensions (recall that we make no rheological distinction between the dermis and granulation tissue at this stage).

144

R.T.

TRANQUILLO

AND

J. D . M U R R A Y

Even these simplifying assumptions do not uniquely define the constitutive stress tensor, and we choose here a standard form (Landau & Lifschitz, 1970). This stress tensor is modified heuristically to account for traction-exerting cells distributed in the ECM by adding a term which accounts for an isotropic active (traction) stress. We assume it is proportional to the product np (implying a rapid equilibrium of pseudopod-fiber anchorage interactions or receptor-mediated adhesions), where the proportionality constant, r, reflects the traction stress generated per unit cell and ECM concentrations. To accommodate the in vitro observation of contact inhibition of motility (Abercombie, 1970) and, therefore, of traction generation, r is taken to be a monotonic decreasing function of cell concentration, specifically r = r0/(1 + yn2). The functional form of the traction stress term is an open question, as is its validity for fibroblasts in vioo. Indeed, Gabbiani and co-workers noted that myofibroblasts possess numerous intercellular gap junctions and desmosomes at high cell concentrations which may act to synchronize and possibly enhance rather than inhibit traction. Alternative functional forms can easily be assessed within the framework of this model. Another force to be considered is the reaction force on the celI/ECM composite due to forces transmitted through attachments of the ECM to subdermai structures (e.g. muscle fasciae). This would increase beyond one the number of space dimensions required to describe the model wound (see below). To circumvent this added complexity, we model this force as a linear, elastic restoring body force, -spu_, on an element of the composite rather than as a boundary condition in a further space dimension. The elastic constant, s, can therefore account for the presence or absence of the panniculus carnosus in the simplest manner by taking on large and small values, respectively. With these assumptions, the equation governing the displacement of the celI/ECM composite becomes: V. [p, ~t +/~2 ~t0 1 + ( ( 1 - ~ v ) ) [ e + ( ~ 5 2 v ~ ) 0 , ] + r0 +1 ~ n 2 l]=sp_u

(4)

where u = x - _x0= displacement vector = ½(Vu + V_uT) = strain tensor 0 = V. _u= dilation and ~, and g2 are related to the shear, g, and bulk, K, viscosities (/1, =2/~, g2= K - 2 / 3 g ) , E is Young's modulus, v is Poisson's ratio, I is the unit tensor, and X_o denotes the initial position of the material point located at position x at time t. The linear stress-strain assumption only applies for sufficiently small strains. Further, anisotropy associated with stress-induced fiber alignment is not admissible [this to be distinguished from fiber alignment occurring local to a cell due to its tractional structuring activity (Stopak & Harris, 1982), which is not relevant for this macroscopic description and thus not inconsistent with the assumed stress tensor]. However, in wounds which exhibit significant contraction, strains are considerable

CONTINUUM MODEL OF WOUND CONTRACTION

145

and alignment of ECM fibers often occurs. Thus, a strictly valid comparison of the model predictions can only be made with wounds exhibiting minimal contraction. Similarly, significant net cell orientation associated with directed migration (e.g. chemotaxis) is not admissible given the assumed isotropic form of the traction stress. The occurrence of osmotic forces which may develop in granulation tissue, particularly in the early stages of wound healing when hyaluronic acid is present in elevated concentration (Weigel et al., 1986), are neglected here, although they can be readily incorporated into the mechanical force balance (Oster et al., 1985). In summary, the base model assumes the following cell and ECM properties: fibroblasts: ECM: celI/ECM composite:

random migration, passive convection (with ECM), logistic growth passive convection isotropic, linear viscoelastic solid with distributed cell traction stress and elastic subdermal attachments.

To compare predictions of this base model and the extensions that follow with the results of McGrath & Hundhal (1983) for full-thickness excisional wounds, the initial state of the model wound is taken to coincide with the state of the experimental wound following the immediate expansion phase. ECM is assumed to be present in the wound initially and devoid of fibroblasts, which are accurate assumptions given the rapid formation of the fibrin clot. Further, the ECM concentration in the wound is assumed to be the same initially as the dermis outside the wound boundary, defined as x = L, an assumption of unknown accuracy (x = 0 defines the wound center and L is the characteristic linear dimension of the wound). The wound is considered to be suitably isolated so that the surrounding dermis can be approximated as an infinite medium. These conditions for a one space-dimension model wound (_u is described by a scaler variable, u) are consistent with the following initial conditions: n(0, x) =0,

O0

u(0, x ) = 0 ,

x>0

(s)

and boundary conditions: On(t, O) - - - 0 ax

On(t, oo) ---~0 ax ap(t, oo) = 0 dx

u(t, 0 ) = 0

(6)

du(t, oo) ---0. 0x

The initial state of the wound is equivalent to a large perturbation of the uniform steady state (n = No, p = Po, u = 0).

146

R.T.

TRANQUILLO

AND

J, D .

MURRAY

Two wound geometries are considered here, a linear wound, which is consistent with long parallel wound boundaries, and a circular wound. Consistency with these one-space dimension descriptions also requires that the wounds are such that the granulation tissue is much greater in depth than width. Modeling the subanatomical attachments as an elastic body force relaxes the large depth requirement to some extent. We define the following dimensionless variables:

n H=--

~ =P-P-

No

1,/=-u

. ~ = x-

L

L

Po

I

?-

(7)

L 2/ Do

so that the base model equations for the linear wound (omitting tildes for notational simplicity) become:

--

Ot

ff ~

Ox

Ox 2

~-ron(1

-

(8a)

n)

=o

(8b)

O3u O2u \1 + ?'n / I~ ~-;7--C + ~-5 + = s pu 0 x Ot Ox 3x

(8c)

Ot

Ox

and for the circular wound (where x =1", the radial co-ordinate):

On

1

1 \

--'Jr

Ot

-

-

x

--

Ox

x

- -

Ox/

c3x

4- ron( 1 - n)

(9a)

0 xp

3__p_p_I 1 Ot

=0

x

3 ( x O2u]

/1

Ot Ox/

X

Ox

(9b)

Ox

f-

I

x

L °") Ox/

la Ou

u

Ox

x 2 3t

x2

\

+

\l+),n/ Ox

=spu

(9c)

CONTINUUM

MODEL

OF

WOUND

CONTRACTION

147

with the following dimensionless groupings: L2

e0 = r 0 N 0 - -

77 = r N 0 ~

Do

/~ =/-tl +/J2 ^ L2 v E-Do

_~oPoNo r0

E

~9

~ spoL2 =

E

I)

1)-

(1 + v ) ( 1 - 2 v ) 1-v '

(10)

[tildes for these groupings are suppressed in (8) and (9) and hereafter].

2.2. N U M E R I C A L S O L U T I O N M E T H O D S

The full set of eqns (8) or (9), is a nonlinear coupled parabolic-hyperbolic-elliptic system and the associated conditions, (5) and (6), define an initial-boundary value problem. Even for these base model equations, analytical treatment of the simpler boundary value problem describing the steady state(s) of the model wound is not yet possible. Consequently, numerical solutions were obtained. The Numerical Algorithms Group (NAG) subroutine D03PFG was used for the initial-boundary value problem by defining the auxiliary variable v = Ou/Ot. It uses the method of lines (finite difference approximation of spatial derivatives and numerical integration in time of the resultant system of ordinary differential equations) on a fixed spatial domain using Gear's fourth-order formula. Accuracy is controlled for integration in time (the absolute value of the estimated error for each variable at each node was required to be at least 10-4 ) but not for the spatial finite differencing. Details of the code can be found in Dew & Walsh (1981). The normal number of nodes (201) was doubled to confirm accuracy by observing negligible change in the solution. The infinite domain consistent with the isolated wound assumption was approximated by requiring that the variables approach their respective "undisturbed" values (n ---' 1, p ---' 1, u--, 0) with zero slope at a "far" boundary (e.g. x~ = 10). Larger values ofxoo were used for sets of parameter values yielding solutions violating this requirement. Similar remarks apply to the use of the NAG subroutine D02GAF (finite difference approximation of spatial derivatives and Newton iteration with deferred correction for the resultant system of nonlinear algebraic equations) used for numerical solution of the associated boundary value problem. For the initial-boundary value problem, a nonlinear least squares regression fit of predicted values of the relative wound area over a series of time, A(ti), to (1) based on a time series of displacements, u(l, ti), was accomplished using the IMSL subroutine R N L I N (modified Levenberg-Marquadt method) to determine regression values for Ao, At, and k, the dimensionless contraction rate constant (k=kcL2/Do). There is variability in the goodness-of-fit among the different parameter sets, indicating that the exponential fit is more or less satisfactory depending on the parameter values. Also, the regression value for k can be quite sensitive to the number of included points u(l, b) at long times when the steady state is approached. Most of the fits underlying the parametric sensitivity plots (e.g. see Fig. 5) are at least as satisfactory

148

R.T.

TRANQUILLO

AND

J. D . M U R R A Y

as the one presented later in Fig. 3(e); however, variations in k of up to ±0. l can be associated with different time cut-offs. Since the displacement, u, defined in the conservation equations above is in reference to spatial (Eulerian) and not material (Langrangian) co-ordinates, u(l, t) is the displacement of the material point at x = 1 at time t, which except for t = 0 is generally not the material point located on the wound boundary. The position of the wound boundary, xb(t), which is the measured variable (Fig. l), is defined implicitly by xb(t) = 1 +u(xb(t), t).

(11)

However, since only small strain solutions were considered (see below), the use of

u(1, t) in lieu ofu(xb(t), t) is a good approximation. This was confirmed by obtaining comparable values for k and ,,If using xb(t) obtained by solving (11) and by using the approximation based on u(l, t) [e.g. for the case illustrated below in Figs 3 and 4, k=0.944 vs. 0.910 and/ff=0-912 vs. 0-906 using xb(t) vs. u(1, t), respectively]. 2.3. P A R A M E T E R V A L U E S

There is little data currently available from which to define values for the model parameters, leaving a large parameter space for characterization. However, the desirable feature of the model that unwounded skin should be stable to small perturbations serves to constrain the permissible parameter space. Given the morphogenetic properties of the Oster-Murray model equations, this is not a trivial consideration (e.g. Murray & Oster, 1984). A linear stability analysis of the uniform state (i.e. unwounded dermis) for the base model results in the following necessary and sufficient conditions for instability: r0 ) < 0

s+ro 1 - 1 + 7 , or

(12) 2T0

l+pr0

-(1- 0-85 for the cases investigated. 2.4. R E S U L T S

Although the base model has a stable, spatially non-uniform steady state evolving from the prescribed initial wound state, (5), its qualitative features are not consistent

CONTINUUM

MODEL

OF

WOUND

CONTRACTION

149

with those of a contracting wound for the regions of parameter space examined. Typically, there is a transient retraction of the wound boundary outward (u > 0) and relaxation inward back to the undisturbed state (u = 0) rather than any contraction of the wound boundary inward (u < 0). An example is given in Fig. 2. The population of the wound by fibroblasts resulting from random migration and proliferation is seen in Fig. 2(a), restoring the cell concentration to that characteristic of the dermis (n --~ 1). The absence of any contraction is seen from the displacement profile in Fig. 2(b) [i.e. u(1, t ) > 0 for all t]. It thus seems that the base model set of fibroblast and ECM properties, which we consider to be a minimal set of properties applicable for dermal wound healing, are inadequate for simulating wound contraction. 3. Model Extensions

All extensions of the base model add to the fibroblast functions, ECM properties, and wound healing events assumed for the base model. Except for the "ECM biosynthesis" extension, they are motivated by the widely recognized influence of inflammation-derived biochemical mediators on fibroblast functions. Mathematical models of inflammation provide a basis for predicting the concentration field of an inflammatory mediator around the wound (e.g. Lauffenburger & Alt, 1987). However, we have chosen here to circumvent the added complexity of a variable concentration field of the mediator and simply assume that a steady-state concentration exponentially decreasing from the wound center exists over the time-course of contraction. A future paper will explore the causal factors in wound healing and contraction by accounting for the dynamics of such mediators. Here the mediator concentration, c, is taken to have the following dependence on distance from the wound center: c=c0 exp ( - -~)

(13)

co is the concentration at the wound center, c and Co are scaled to the concentration for half-maximal enhancement of the cell behavioral property [e.g. (14)]. tr is a parameter specifying the spatial domain of influence of the mediator from the wound center. In all results presented here, we take co = o- = 1.

3.1. T R A C T I O N V A R I A T I O N

Among the various ways an inflammatory mediator could influence the behavioral properties of fibroblasts, it has been proposed that fibroblasts transform into traction-enhanced myofibroblasts around the wound due to some inflammationderived growth factor (Skalli & Gabbiani, 1988; see also Montesano & Orci, 1988). A simple alternative to modeling the population dynamics of the myofibroblasts explicitly is to make the traction parameter, r0, of the fibroblasts a function of the mediator concentration, c. So we assume the following saturable dependence:

0,015

0,020

0-025

0.0 0.0;30

0-1

0.2

\

1

J.

-0.005

0.000

0-005

0.0

I

2.0

I

1.0

//,,:-\,, \\ \ .

// \\

/ \ \ /r\ \

/..-.,

I

3,0

1

4.0

1

5.0

I

6,0

(h)

(o)

/

o 0,r:'-\\,\ \

0.04

(x)

"-0-02 0.O

-O-OI -

Di•tonce

7-0

A 4,%

/

O'051N

0.955

0,965

0.975

0-985

0.995

1,005

1.0i5

1.0

I

I

I

2.0

I

I

3-0

I

I

4-0

l

5-0

I

...__..=

I

_

I

6-0

I

(d)

(c)

7.0

FIG. 2. Representative prediclions of the "base model" for a linear wound: (a) fibroblast concentration, , ; (b) displacement of the celI/ECM composite, u; (c) ECM concentration, p; (d) and strain of the cell/ECM composite, ~u/ax, are plotted vs. distance from the wound center for a series of times after wounding, x = 0 is the w o u n d center and x = l is the initial location o f the wound boundary. Since the profiles are symmetrical about x = 0 , only the results for x>__0 are shown. Parameter values are r0 = I , / l = I, s = I, t o = 0 . 5 , 7 = I.

el 0

.........._._~ ~

4"0

5.5

2-5 3.0

0-3

Z-O

I-5

I-0

0.5

0.0

Time ( t )

0.4

0.5

0-6

0.7

0.8

o" O-OIC Q

E

g

o

o

.0

Im

A

0.9

1.0

I.I

CONTINUUM

MODEL

OF

WOUND

CONTRACTION

151

where r/is a traction enhancement coefficient. Another justification for making r0 dependent on the concentration of an inflammatory mediator is the observation that fibroblasts in vitro express an increased number of integrins upon exposure to transforming growth factor-fl (TGF~) (Heino et al., 1989). Integrins are cell receptors for cell adhesion molecules such as fibronectin, which is known to be abundant around the wound during the early phases of wound healing (Grinnell et al., 1981). T G F n is a growth factor which elicits many of the same cell behavioral effects as platelet-derived growth factor (PDGF) and macrophage-derived growth factor (MDGF), P D G F being present in largest amounts during the platelet rupture associated with initial clot formation and M D G F presumably being present in largest amounts over the latter stages of inflammation which are dominated by wound macrophages (Huang et al., 1988). An example of the predicted profiles for this case is given in Fig. 3 for the linear wound. The population of the wound by fibroblasts resulting from random migration and proliferation is again seen in Fig. 3(a). The evolution of the celI/ECM composite to a contracted steady state [i.e. u(1, t) < 0] is apparent from the displacement profile in Fig. 3(b) and the simulated time-course of displacing tattoo marks in Fig. 3(f) [derived from Fig. 3(b)]. Associated with the contraction is an increase in the ECM concentration [Fig. 3(c)] and development of negative (compressive) strain [Fig. 3(d)] in the wound. The plot of the wound area (based on the location of the wound margin) over time in Fig. 3(e) is qualitatively consistent with the data in Fig. 1, suggesting that the fibroblast-ECM-inflammatory mediator mechanisms of this scenario can interact to yield the empirical exponential fit which McGrath & Simon (1983) characterized. The various stress and force profiles during the time-course of contraction are presented in Fig. 4, showing the decreasing viscous stress as the rate of contraction slows, and the increasing traction stress (associated with increasing cell concentration in the wound) balanced by increasing elastic stress in the ceil/ ECM composite and increasing elastic subdermal attachment restoring force. Given that these results are in accord with the simple criterion of contraction, it is of interest to determine how other parameters of the model affect the qualitative and quantitative features of the predictions for this case. The dependence of the dimensionless contraction rate constant, k, and the wound area remaining at steady state, As, on the traction enhancement factor, ry, are presented in Fig. 5. A i decreases for increasing r i [Fig. 5(a)], as expected since a greater traction exerted per cell suggests an increased final extent of contraction. Figure 5(b) indicates only a weak dependence of k on Z~rover the range investigated. A stronger dependence might occur for larger rr since that would generate greater convective velocities which directly affect the dynamics reflected in k, but larger rr yields involves strains too large to be consistent with the current constitutive stress tensor. A useful result can be obtained under the assumption that strains are sufficiently small so that the ECM conservation eqns, (8b) and (9b), can be linearized: p = l ---

OU Ox

(8b')

0-." 0-02

o.~

o.I

0,! O~

0.1

o.!

-0-I0

-0-08

-O.OG

'o

~

"'r-._

i

J

i

,'.o

/-"

~!o

,'.o

,'.o

,.o

3"5 4,0

\

z.o

~.0

o

"y*

.--Z

J

(t,)

[a )

o,ltonce

(x)

-o.ro

O.~4~/ o,10.

/

I ' ~l ~ ) l

I

\,,

I

I

I

I

I (dl

(c)

.~

-

0.0"

5-¢

6
1 is associated with wound expansion, A < 0 is associated with wound contraction). The dashed line is the nonlinear least squares regression fit of the numerical data u(I, tj) to the exponential form of (1) investigated by McGrath & Simon (1983). In ( f ) the movement of tattoo marks located across the wound and adjacent dermis is simulated for a series of times over the course of contraction. Parameter values are ro = !, p = 1, s = I, ro = 0- 5, 7"= 1.

i

-0"04

?'7"',.

.:-':',,-'

0,0

O-!

z -o~

~

i g

=

-~

I-('

> ,
Z

C t-

> Z

bO

CONTINUUM

MODEL

OF

WOUND

153

CONTRACTION

1.2 (a) T: I.C

(c)

I. 0

r=].O

Tractlo~ slress Viscous stress E l a s t i c stress .....................

0.8

Restoring force O'E

0-4

0.2

= .

O.C

I

-0.2 I-2

.

I

.

.

.

.

.

.

-

.

I

"~'--_-"S- . . . . I

I

I

(b)

(d)

r=2,0

r,4.0

I-0

0-8

0.6

0.4

0.2

0,0

,...,~ .-r:f--""

......--*

-0-2 0.0

1.0

2 0

3.0

4-0

I bO

I 2.0

I 3.0

4.0

Distance ( x )

FIG. 4. Profiles for the various stresses and the elastic restoring force for the case of the "variable traction" extension presented in F i g . 3.

in the linear case, and p = 1--

1 Ou

-x Ox

(9b')

in the axisymmetric case (a lucid account of the self-consistent use of the small strain approximation can be found in Lin & Segel, 1974). Note that a uniform steady-state solution n,(x)= 1 satisfies the cell conservation eqns, (8a) and (9a). Substituting this along with the above linearized approximations for p into the mechanical force balances, (8c) and (9c), yields boundary value problems for the steady-state displacement, u,(x), i.e. the final contracted state o f the linear and circular wounds given the

154

R

T.

TRANQU1LLO

AND

J.

D.

MURRAY

(a

Co=l /z=l 5 =1 ro -0.5 y=l

0.98 A

(-0.027,0.038) ~ 0.96

"o..

0,94

g ~ c

"e.

0-92

t,

(-Oq 3 6 , 0 . 0 5 5 7 0 \ \

0.9

"ID\ \

0.88 4

l

l

[

I

l

l

l

l

l

l

[

l

l

I

l

t

I

l

I

1

I

I

J

I\

(b) 3.5

~

2.5

iiJ

P, ~

1.5

(.) 0 - ~ - - 0

0 - - - - 0 - - -

•

-0------

0 . .~ I

0

I

1

I

I

i

0-5

1

i

i

I

1

i

I

i

I

I

I

I

I

i

1,5

Traction force enhancement

)

2

I

I

I

I

2.5

coefficient ( T t )

FIG. 5. Dependence for the "variable traction" extension of (a) the final wound area at steady state ("uncontracted" area), As, and (b) the dimensionless contraction rate constant, k, on the traction enhancement coefficient, rs, Other parameter values are the same as Fig. 3. The range of strain is indicated below selected points.

sustained concentration gradient o f inflammatory mediator. The equation describing u,(x) for the linear w o u n d is given by:

I + ~, a--Zx~- su'~ ~

,Ix/

1 + ~,

dx

The displacement profiles obtained from numerically solving (15) agree very well with those obtained from numerical solution o f the full system o f equations, (8a-c),

CONTINUUM MODEL OF WOUND CONTRACTION

155

as a confirmation of the latter, for the values of r0 and rfinvestigated in Fig. 5 where traction-induced strains are sufficiently small. The same is true for the corresponding equations for a circular wound, (9a-c) and the axisymmetric equivalent of (15). Choosing r f = 2 , we consider the predictions for two potentially-manipulatable parameters in dermal wound therapy: r0, representing the ratio of characteristic times for fibroblast population of the wound by random migration across the wound space to proliferation, and/1, representing the ratio of viscous force to elastic force associated with displacement of the wound margin during the characteristic time L2/ Do [ro and p are defined in (10)]. We would expect both ro and p to influence the results for k, but not very significantly for Af given their absence from the linearized boundary value problem for u~(x), (15). The results of Figs 6 and 7 are in accord with these expectations. Only a weak dependence of Af on 1"0[Fig. 6(a)] and p [Fig. 7(a)] is observed. The dependence of k on ro [Fig. 6(b)] is consistent with the definitions of k and r0. For r0 decreasing from 1 (i.e. Do increasing), where the rate of populating the wound becomes increasingly determined by fast random migration, k becomes independent of r0. Recall k=kcff/Do, so the dimensional contraction rate constant, kc, is increasing as r0 decreases. Thus, the rate of wound contraction is not limited by the relaxation rate of the cell/ECM composite in this case. For r0 increasing from l, the population of the wound is not mainly determined by random migration due to the contribution from proliferation, and k increases. For large r0 the rate of populating the wound becomes increasingly limited by slow random migration, and k again approaches an asymptotic value. Thus, k~ is now decreasing as r0 increases. Likewise, the dependence of k on p [Fig. 7(b)] is consistent with the definition of p : for decreasing values of p, the rate of wound contraction becomes less dependent on the relaxation rate of the cell/ECM composite and increasingly determined by the rate of fibroblast population of the wound, reflected in the increase of k to an apparent asymptotic value (in this inviscid limit, which is computationally more difficult, the mechanics of wound contraction become governed by a balance of only elastic, restoring, and traction forces at each instant, with only the last being explicitly time-dependent via its coupling to the time-varying cell concentration); for increasing values of p, the converse is true, with k decreasing asymptotically to 0. It is of interest to compare predictions from this case with the findings of McGrath & Simon (1983) that kc and ,4i are essentially independent of initial wound size and geometry. If parameters for the illustrative case of r l - - 2 in Fig. 3 (k=0.910, AI = 0-906) are accordingly adjusted to consider a linear wound of one-half the initial area, Ao/2 (i.e. L/2), the model predictions are k = 0 • 196 and Ar=0-826 (Table 2). Accounting for the scaling dependence of k on L [since (1) in terms of dimensionless time is used to fit the numerical data), the expected value of k, if k~ is indeed independent of geometry, is k = 0.222. Thus, in the case of a linear wound, the model predicts a slight dependence of k,. on initial wound area (the relative change in k~, Ak,. = - 12%), at least for this region of parameter space (a change of this magnitude may not be significant given the variability in k associated with an imperfect exponential fit which was noted earlier). This is consistent with the lack of dependence found by McGrath & Simon for a square wound with one-half the initial area (Table l; Akc = + 18% was determined not statistically significant, nor were any of the values

156

R.T. I

TRANQUILLO

AND

J.

D.

MURRAY

/.L=I 5=1

(a)

0.98

X=I

=0.5 rt=2

ro

0.96 0.94 § E

0,92 I

iT

0"9I 0-88

•

i

, ,,

ill

I

I

I

I

I IIII

I

I

I

I

I

I

IIIII

J

,

I

I II

4 5.5 5 E 0

2-5 0 V 0

1 2~ 1.5 O~/0--

E 0 U

.i

-- • . . . .

---O--~ •

I

I-

0"2I O-Ol

.... Illl

I

o.I

I

I

I

Illll

I

I

I

I I IIII

I

I

I0

I

I

I I III

I00

Cell growth rate constant (ro)

F I G . 6. Dependence for the "variable traction" extension of (a) the final wound area at steady state ("uncontracted" area), A/, and (b) the dimensionless contraction rate constant, k, on the dimensionless cell growth rate constant, to. Other parameter values are the same as F i g . 3.

for Akc or AAI). Comparisons of the wound area are more difficult to interpret since the values of AI observed by McGrath & Simon are much smaller than those predicted by the model under the small strain constraint: comparing values of AAI, the relative change in final wound area, the model predictions agree with McGrath & Simon's conclusion that the extent of contraction is independent of initial wound size; agreement does not exist for AAc values, the relative change in contracted area, the model predictions showing a strong dependence on initial wound size. For a circular wound with the same linear dimension as the linear wound (i.e. same diameter as width) the model predicts k=0.952, quite close to the value predicted

CONTINUUM

MODEL

OF WOUND

157

CONTRACTION

(o) $=1 To.o.5 ),=1

0.98 E b

rf-2

0.96 0.94 0-92

u_

0.9 1

•

0"88

I

I

i

I I Iltll

41 (b) 3.5

I

I

I

iiiii

1

I

I

I

ii

Ill

i

i

i

1

•

A

"" ~ 0

v

3

~

~

c o

§

\

\

2.5

\

\ \

o

g

\

1"51

\ I

g

\

0

0.5 0

).01

i

I

i

i

illil

O-I

I|

I

CelI/ECM composite

IO viscosity

iO0

(#)

FIG. 7. Dependence for the "variable traction" extension of (a) the final wound area at steady state ("uncontracted" area), As, and (b) the dimensionless contraction rate constant k, on the dimensionless celI/ECM composite viscosity, #. Other parameter values are the same as Fig. 3.

for the linear wound (Akc=-4.6%), in accord with McGrath & Simon's finding concerning the independence of kc from wound geometry ( A k , = - 12%). Again, it is difficult to draw conclusions from comparing the AA/and AA~ values for this case. Although comparison data are not available, model predictions were obtained for a circular wound with initial area Ao/2 (i.e. L/2~/2), where k = 0.422 (vs. prediction of 0.455, Ak, = - 7.2%), and for initial area 2A0 (i.e. 2½L), where k = 2.12 (vs. prediction of !- 82, Akc = + 16%). Thus, the model is consistent with their conclusion concerning independence of k, from A0 for circular as well as linear geometry. As for the linear

158

R. T. TRANQUILLO AND J. D. MURRAY TABLE 2 "Traction Variation" model predictions

Wound Large circular (A =Ao) Large linear§ (A =Ao t) Small linear 'l (A=AÜ/2) Small circular§ ( A = A t / 2 ) Very large circular§ (A =2Ao)

k*

kpr~ *

hk~

At t

A, *

AA r *

AAt

0.952 0.910 0.196 0.422 2.12

---0.222 0.455 1.82

--4.6% -12% -7.2% +16%

0.858 0.906 0-826 0-814 0.900

0.142 0.904 0.174 0.186 0.100

-+5.6% -8-8% -5-1% +4-9%

--34% +85°/, +31% -30%

* k is a dimensionless form of kc ; kpred is the predicted value of k assuming k, is constant. t These values are relative to the maximum area pre-contraction. Based on width equal to diameter of large circular wound. § A values are relative to large circular wound. II A values are relative to large linear wound.

geometry, the predicted values of AAr for the circular geometry are consistent with their conclusion that AI is independent of Ao, whereas the predicted values of AA,. are inconsistent with their data indicating that A,. is independent of A0. It is reasonable to question whether their finding that k,. and Arare independent of Ao, based on square and circular wounds should, in fact, apply to linear wounds as well, since now the geometries are fundamentally different. Our model predictions suggest so, but it remains to be verified.

3.2. G R O W T H V A R I A T I O N

Given the feasibility of topical application of patient-derived growth factors to dermal wounds (Knighton et al., 1989). it is of interest to understand the consequences of a concentration gradient of growth factor(s) around the wound which affect the growth characteristics of infiltrating fibroblasts (as noted in "traction variation", a growth factor may affect cell behavior in ways besides modulating proliferation, e.g. traction capacity). Unfortunately, there are no in vivo data from which to deduce how ro and No in the logistic growth rate term of (2) may depend on growth factor concentration, nor to our knowledge in vitro data from cultured fibroblasts in tissueequivalent collagen gels. Here we assume for illustration that fibroblasts will proliferate to a higher maximum cell concentration in the presence of a growth factor, and take the functional dependence of No on c again to be a simple saturable dependence as was done for r in (14):

oE,+

(16)

To meet the criterion of contraction, it was necessary to reduce y (the parameter measuring the extent of contact inhibition of traction) from the value 1 assumed until now to be O. 1 (see below). To facilitate the computations, p was increased from 1 to 10. Otherwise, the parameter values are as for the previous cases (n.b. rf--O

CONTINUUM

MODEL

OF

WOUND

CONTRACTION

159

applies here). As seen in Figs 8 ( N I = 5 ) and 9 ( N I = 10), the profiles are distinctly different from the "traction variation" extension, notably the increased cell concentration in the wound. Also, at least for these parameter values, the point of maximum displacement occurs in the dermis (at x~2-5) rather than at the initial location of the wound margin ( x = 1) [cf. Fig. 3(b)]. Unlike the "traction variation" results, the dependence of AI on the enhancement factor is not monotonic, an optimal Ns = N ~ 5 existing for minimum As, [Fig. 10(a)]. This requires further explanation since contraction (i.e. displacement) with respect to the point x = 2 . 5 is almost the same for N I = 5 and N r = 10. The reason for the biphasic AI dependence and non-intuitive k dependence [Fig. 10(b)] is that the displacement profile around the wound boundary is quite sensitive to Ny [cf. Figs 8(b) and 9(b)]. This is due to the assumed occurrence of contact inhibition of traction. For a constant ECM concentration, the traction stress term is easily seen to have a maximum value of (4X) -I/2 per unit ECM concentration at a value n* = X-t/2, or n* = 3"2 for the value X =0.1 illustrated here. Thus, traction inhibition should be significant for Ny= 10 since n ~ 5 > n * in the wound, leading to smaller displacement of the wound boundary relative to Ny= 5 where n ~n* in the wound. Additional computations show that the value of N~ is closely tied to the value of X, as expected. As before, the final contracted state obtained from the long time solution of the initial-boundary value problem agrees well with that determined from the corresponding linearized boundary value problem, here a system of coupled ordinary differential equations for ns(x) and us(x), but only when Ns< 3. Convergence to a meaningful steady-state solution was not obtained for NI> 3. In fact, there is evidence of singular behavior for these cases in the numerical solution of the time-dependent equations, with a growing maximum in ECM concentration (coinciding with a growing minimum in strain) occurring near the wound boundary [Figs 8(c) and 9(c)], suggestive of a pathophysiological occurrence. Unfortunately, the model cannot address the long-term behavior in these cases because of the small strain constraint. Note that the present form of the traction stress term does not allow for traction inhibition at high ECM concentration, which makes this singularity admissible. However, we suspect that in reality, the pseudopodal activity of a fibroblast underlying traction becomes inhibited when surrounded by sumciently high ECM (fiber) concentrations.

3.3. C H E M O T A X I S

Fibroblasts are considered to be chemotactic to a number of factors which are generated in the wound [i.e. ceil-derived chemotactic factors or proteolytic ECM fragments (Aibini et al., 1985; Grotendorst et al., 1985)], and it is typically assumed that chemotaxis is an important mechanism by which fibroblasts, as well as most other cells involved in wound healing, localize in the wound. Although the modeling of chemotaxis is far more developed than for other tactic mechanisms, it is still not a completely solved problem (Tranquillo & Lauffenburger, 1990). Here we assume the simplest description of chemotaxis, based on adding a chemotactic flux to the

-O-O9

-0-O4

0-0

0"°1 l

o.ol

~-0 0-5

"

"

N

1.0

2-0

~

%,

,

~

,

~

3.0

4.0

:¢...]

~0

S.O

7.O

,

I

8.0

___~_

_.Z._..e__.

z.~

_oL_

~o ~.~,

I

t )

---*---

I

0

Time( I

{hi

~.o

i

(o

I

,

0.0

- 0-OP!

I

t-O

~

I

,

2.0

I

---

I

I

I (*ll

1-0

i

4.0

i

5-0

i

a-O

i

7.0

i

g-O

,

t.O

,

. . . . . . . . . . . . . . . . . . . . . . . . . . .

~.,~____~_

*

*

J

--- _ ---~=

' , '.,,',.-.i~:-Z~.-,.~"

-Q.OSI-~'~/~/~ k . .

-0-02

o*~tQnc* ( x )

tO.O

"-

.

t'7-

~o,

0"04 1

O.H[

o.*~[

•

I0-0

1.0;

O.t$

0.94

0,gS

O.g6

(>s7

O.Se

O.Ss

,0 0

O.0

, 1.0

2.0

i

!gO

,

T~m* { r )

4.0

,

5.O

,

60

,

7.O

,

FIG. 8. Representative predictions of the "growth variation" extension for N / = 5 and a linear wound: see Fig. 3 legend. The profile of the maximum fibroblast concentration. No, associated with the inflammatry mediator gradient is indicated [not to scale; actual range is 1.0-3.5, sec (13) and (16)]. Other parameter values arc ro -~ 1, # = 10, s = I, t o = 0 . 5 , ~, = 0 . !.

i

~-_~

~

~

4.0

4.S

5.0

>

/:

O

Z O

t-" O

N

JO

>

,H

,H

o~

CONTINUUM

MODEL

OF WOUND

161

CONTRACTION

E

e~

o dl

m

~r °

~

i

"

" v~ o

0

o I

o m c,

Ii

II

'k , 4 ;:",,'I~ ~ I

.2

":(._'.~ ~. ( i

~. . . . I"" " " ' "

/--",

,~.~-.~

,

?

i

,..,

; ~ ,

0

2 6,-k

~4 .'d,

/ I ~ / I i

\

','k.~a.-. I '~ ,\'S,~,:

~

-~.o

--

0

t~

162

R. T. T R A N Q U I L L O

AND

J. D. M U R R A Y

(a) 0.9~

//e /

\ o

o

\

0.96

\

• //

5\ -0.044,0.015) ~

0-94

/

(-0.075,0-050)

/

/

/

g v c D

-6

(-0.065,0-025) •

~•~

~"~

(-0-075, 0-036)

0.92

ro:l /a=lO S=I

c h

ro:O. 5 7" :O.I

0.9

i

0.88

i

I

I

i

I

I

l

l

n

l

l

t

l

l

l

l

l

l

l

n

t

l

(b) :3.5 4c

2.5 o

2

•ru

1.5

c

o

t.)

I

/O / J

0-5 n

L

0

a

[

n

t

2

I

I

t

4

t

n

[

6

I

t

I

I

t

i

l

8

I

I

rO

t

I

12

Moximum cell concentration enhoncement coefficient (/Vr)

FIG. 10. Dependence for the "growth variation" extension of (a) the final wound area at steady state ("uncontracted" area), At, and (b) the dimensionless contraction rate constant, k, on the maximum cell concentration enhancement coefficient, NI. Other parameter values are the same as Figs 8 and 9.

random migration flux yielding [cf. (2)]:

On+V • n __zOu= Ot

Ot

V- (DoVn - y nVc) + ron(No - n)

(17)

where c is the chemotactic factor concentration, described as before, (13), with a steady-state exponentially decaying concentration, and Z is the chemotaxis

CONTINUUM MODEL OF WOUND CONTRACTION

163

coefficient, reflecting the directional bias exhibited by a cell in a chemotactic gradient [generally, Z will be a function of c, as documented for several leukocytes (e.g. Tranquillo et al., 1988), but here we take it to be constant for simplicity]. As for the "growth variation" extension, X=0 • 1 is used (n.b. r i = 0 and N,~=0 apply here). Recognizing that chemotaxis serves as a mechanism of fibroblast accumulation, it is expected that results for this extension, presented in Figs l l and 12, would share many similarities with the "growth variation" extension, for the same reasons. Over the range of Z examined (Z = 0-1.2), the dependencies of k and Af o n 27 are qualitatively similar to those for r I (cf. Fig. 5). Although dependencies similar to those for NI (cf. Fig. 10) appear for larger values of 27 (greater chemotaxis), the small strain assumption becomes more violated and the results less valid. There is also considerably more strain at the wound center, where a singular increase in ECM concentration is indicated for larger values of ZIn passing, note that there is certainly an orthokinesis aspect offibroblast migration in wound healing (i.e. dependence of cell speed on the local magnitude of a stimulus), as fibroblasts resident in the dermis are believed to be essentially non-motile (zero speed) in the absence of any stimulation. An account for this requires yet another cell flux term (Tranquillo & Lauffenburger, 1990), but the theoretical consequences are well known : if the cell speed only increases monotonically in the direction of the wound, orthokinesis will counteract the accumulation of fibroblasts in the wound associated with chemotaxis; if it increases then decreases in the wound, it will enhance the accumulation.

3.4. E C M B I O S Y N T H E S I S

Here, the base model is extended to include an ECM generation term to account for the extensive deposition of collagen and other ECM macromolecules by fibroblasts in granulation tissue. Since the regulation of this process is not yet well understood, a particular functional dependence of the net ECM generation rate on n, u, p and c cannot be independently justified. We assume, solely for the purpose of illustration here, that this rate is proportional to that of fibroblast proliferation (i.e. a logistic growth-type rate), so that there is ECM accumulation until n ~ 1 . In this case, the ECM conservation equation becomes [cf. (3)]: oO,.+V " - = b n ( l - n ) . Ot Ot

(18)

Although this rate form yields a contracted steady state (Figs 13 and 14) without having to invoke the influence of an inflammatory mediator, the fact that, by construction, p attains a steady state before the completion of wound contraction is not consistent with the observation that collagen continues to accumulate for extended periods after completion. A greater understanding of ECM regulation in wound healing is needed to guide any future proposition of realistic rate forms for ECM generation.

/\.o

~

~

.... ,~o

I

....

DJatl,¢! ( • |

-°'o.o

,'o

(c)

I

I (d )

.~o

:.o

,,~. . . .

""........................................

1

............................

_~,.~. ~%'...:~

°'tO k

0.10|

--F -"/"~"~" I

1.20 ~ ' ' ,

" )~

|

.'o .......... ;': ......... ~-o-- ..... ~o

~

(b )

:j._-::

___z:o__

(o)

0-,2 0.0

0.I14

0-H

0.90

1,00

1.02

1,04

1,06

2.5

I!O

l;.5

2!0 Tim

(r)

2!S

3;.0

$!~

4.O

FIG. 1 I. Representative predictions of the "chemotaxis" extension for Z = 1 and a linear wound: see Fig. 3 legend. The profile of the chemotactic factor concentration, c, defined by (13) is indicated (not to scale). Other parameter values are ro = 1, p = 1, s = I, t o = 0 . 5 , ~,=0.1 (same as Fig. 3 except r / = 0 ) .

--'~

-

-°'"~""'I~:-'}~"', ".,.' ;~',---:,.'' --"'.-r-.

-~ "v's---::~

/

" . ] ~ ~

.~ ,~,

o.0([...

O.3

0.6

0.g

I'5

>

O > Z

,O

> Z

CONTINUUM

MODEL

OF

WOUND

165

CONTRACTION

!

(a)

ro=l

#:1 s=l To =0.5

0.98 IL

~'~'~.

,q

y =0.I

0.96 (-0.083,0.075) " ~

o

0.94 {-0d4,0-056) "~ ~.....

§

i--..

0.92

(-0.18,0-070) " " ' ~ . .

t~ 0.9

0.88

I

l

l

l

l

t

l

i

l

l

l

l

l

l

z

l

l

l

l

l

l

l

l

4

(b) 3.5 3

§

2.5 2

--~

1.5

g o

U 0"5 I

0

I

I

I

0.2

I

I

e

I

,

0-4 Chemotoxis

I

I

[

I

0.6 coefficient

I

I

I

0.8

!

I

I

I

I

I

I

I

1.2

(X)

FIG. 12. Dependence for the "chemotaxis" extension of (a) the final wound area at steady state ('uncontracted" area), A/, and (b) the dimensionless contraction rate constant, k, on the dimensionless chemotaxis coefficient, X(Xco/Do). Other parameter values are the same as Fig. 11.

4. Discussion This work marks an initial effort aimed at elucidating how the cellular, biochemical, and biomechanical phenomena known or speculated in the dermal wound healing process are Orchestrated to give rise to wound contraction. Given the complexity of the phenomena, compounded by the uncertainty o f many regulatory points, the strategy we have adopted is to establish a continuum model framework based on a set o f cell functions and ECM properties wh/ch seem intrinsic to any scenario o f

hi

-o.0~

.:.

0"0

,."

x\ \ ~ ' \

~

;.o

,'.o

../-~:

Lo

..~_~ ~'~:~"

Lo

-

,o

-~-

o~stoR¢*

,t. . . . . (• )

- 0 , 1 7 / I

....

-0-08

I ..... i-/,\

....

\ ?~

1,0

V

7

1

F ~,~,~

I 2.0

I 3.0

I 4-0

I $.0

I S-o

7.0

/

0.0

0,5

I-0

I-5

2.0 rim*

( r I

2.5

3.0

3.5

4.0

FIG. 13. Representative predictions for the