Mathematical Modelling of Extracellular Matrix Dynamics using Discrete Cells: Fiber Orientation and Tissue Regeneration
Matrix orientation plays a crucial role in determining the severity of scar tissue after dermal wounding. We present a model framework which allows us to examine the interaction of many of the factors involved in orientation and alignment. Within this framework, cells are considered as discrete objects, while the matrix is modeled as a continuum. Using numerical simulations, we investigate the effect on alignment of changing cell properties and of varying cell interactions with collagen and fibrin.
Published in Journal of Theoretical Biology (1999) 199, 449 471 Mathematical Modelling of Extracellular Matrix Dynamics using Discrete Cells: Fiber Orientation and Tissue Regeneration John C. Dallon1, Jonathan A. Sherratt1 and Philip K. Maini2 1Department of Mathematics Heriot-Watt University Edinburgh, EH14 4AS U.K 2Centre for Mathematical Biology Mathematical Institute University of Oxford 24-29 St Giles Oxford OX1 3LB, U.K. July 12, 2000 1 RUNNING TITLE: Models of Extracellular Matrix Dynamics Abstract Matrix orientation plays a crucial role in determining the severity of scar tissue after dermal wounding. We present a model framework which allows us to examine the interaction of many of the factors involved in orientation and alignment. Within this framework, cells are considered as discrete objects, while the matrix is modeled as a continuum. Using numerical simulations, we investigate the effect on alignment of changing cell properties and of varying cell interactions with collagen and fibrin. keywords: cell flux, cell speed, cell polarization, tissue regeneration, collagen production, fibrin degradation. 2 1 Introduction Alignment phenomena are found throughout our environment, in contexts as varied as liquid crystals (Priestly et al., 1975), intracellular actin networks (Small et al., 1995), fish and insect swarming (Okubo, 1986), cellular biology (N ubler-Jung, 1987) and wound healing (Ehrlich & Krummel, 1996). In this paper, we focus on wound healing, where a higher degree of alignment is a key characteristic of scar tissue. An understanding of how processes interact to produce alignment is a crucial step in the development of anti-scarring therapies. Dermal skin tissue is composed of a relatively sparse cell population surrounded by a fibrous network of proteins called the extracellular matrix. This matrix is primarily composed of collagen that interacts with cells in a dynamic way which is not fully understood (Hay, 1991). It is clear that the extracellular matrix acts as a scaffolding, providing directional cues via a phenomenon known as contact guidance (Guido & Tranquillo, 1993; Clark et al., 1990; Hsieh & Chen, 1983). Fibroblasts, the cells which create and maintain the extracellular matrix, help to orient and give structure to the fibrous network (Clark, 1996). This complex fibroblast-extracellular matrix interaction, known as dynamic reciprocity, is the focus of this paper. We are specifically concerned with the orientation of the collagen matrix and learning what factors are important in alignment. The fibroblasts influence the collagen matrix both biochemically and mechanically. Evidence suggests that as the cells move and produce collagen the fibrils are oriented in the direction of motion (Markwald et al., 1979; Birk & Trelstad, 1986; Trelstad & Hayashi, 1979). Whether or not this is due to tractional forces is unclear, although as fibroblasts move they do exert tractional forces which can alter the structure of the collagen matrix (Stopak & Harris, 1982; Harris et al., 1981). In this paper, for simplicity, we only consider alignment processes which are associated with moving cells and effect the matrix near the cell, leaving an investigation of more long range alignment and contraction for future work. Mathematical modelling of alignment has been prevalent in the past few years, focusing on a wide variety of applications. Of particular relevance to this study is work on fibroblast orientation (Edelstein- Keshet & Ermentrout, 1990; Mogilner & Edelstein-Keshet, 1995; Mogilner et al., 1996; Mogilner & Edelstein-Keshet, 1996); this focuses on orientation due to cell-to-cell contact inhibition rather than directional cues taken from their substrate, on which we concentrate. In contrast to this previous work, we are not primarily concerned with fibroblast orientation as such, but rather on how this affects the orientation of the fibrous network. The work which most closely relates to the modelling in this paper is that which examines alignment of the extracellular matrix (Olsen et al., 1999; Olsen et al., 3 1998; Dallon & Sherratt, 1998; Barocas & Tranquillo, 1997). Each of these represent very different approaches to modelling similar systems. Olsen et al. (1998) model the system with reaction diffusion equations where the species can take one of two orthogonal directions, while Olsen et al. (1999) use a related model including the long range mechanical forces within the fibrous network. Dallon & Sherratt (1998) formulate the long range angular interactions as integro-differential equations, ignoring spatial variation, and in Barocas & Tranquillo (1997) the extracellular matrix is treated as a biphasic medium with viscoelastic properties which both orients, and is oriented, by the cells. All these models use continuum descriptions for variables, whereas the models described in this paper represent the cells as discrete entities and the extracellular matrix as a continuum vector field. In the next section we start our study by describing a simple model involving only one fibrous matrix component (collagen). In Section 3, results from numerical simulations of the model are presented. Section 4 introduces a more complicated model which represents the extracellular matrix as two fibrous networks, one of collagen and one of fibrin. By examining a simple model initially, and then a more complex extension, more is learned about the underlying mechanisms for orientation. The results of the more complicated model are explained in Section 5 and we conclude with a discussion of applications to wound healing. 2 The Matrix Orientation Model Initially we only model simple interactions between fibrous collagen matrix and the fibroblasts. Although we have a collagen network in mind, the model could be applied to any cellular reorganization of a fibrous matrix. The cells are treated as discrete objects whose paths are given by f i(t) = (fi1 , fi2 ), where the superscript denotes cell i. The collagen is represented as a continuous vector field denoted by c(x, t) where x represents the cartesian coordinates of a point in the plane. The vectors c are in 2, have unit length and their direction represents the predominant direction of the collagen fibers at that point in space and time. As mentioned previously, the cells receive directional cues from the extracellular matrix. These cues are modeled by df i(t) dt f i(t) = s v(t) v(t) (1) with v(t) = (1 )c(f i(t), t) + f i(t ) f i(t ) (2) 4 where and s are positive constants with s representing the cell speed and a time lag. The first term in equation 2 represents the influence of the collagen matrix on the direction of the cell ( contact guidance ) and if = 0 the cell moves exactly in the direction of the collagen. Fibroblasts on fibrous gels become very elongated and maintain a leading edge (Friedl et al., 1998) which we model by the second term of equation 2 which gives the cell a bias in the direction it was moving. The parameter determines the strength of this bias, so if = 1, the cell does not change direction and moves in a straight line determined by the initial conditions, independent of its environment. In equation 1 the linear combination of directions is normalized and scaled to give the appropriate speed. The fibroblasts alter the collagen matrix by changing its orientation. As mentioned previously we model only the local flux-induced alignment. Since each fibroblast can contribute to the reordering of the collagen, the overall effect of the cell population on the matrix is represented by the vector f f (x, t) = N i=0 wi(x, t) f i(t ) f i(t ) (3) where is again a time lag and N is the total number of cells. This gives the cumulative effect of all the fibroblasts on the collagen direction. We assume that the effect of the cells acts in this additive fashion; this simplifying assumption is reasonable in view of the low cell density in our simulations. The weight function, wi, interpolates the influence of the discrete cells to the continuum collagen variable and keeps the influence of the fibroblasts local. The issue of how discrete variables interact with continuum variables is discussed in Dallon (2000) which deals with numerical issues arising from discrete and continuum hybrid model formulations. The weight function used here is graphed in Fig. 1 and defi ned by wi(x, t) = a1a2 with aj = max 1 |fi j (t) xj | L , 0 , (4) where x = (x1, x2) and L = 10 m. Thus the support of the function is a 20 by 20 m square. Ideally the support of the weight function should correspond to the shape of the cell, but since it is constantly changing and unknown we use the above simplification. The time lag represents the time taken for the cell to change direction after obtaining the directional cues from its environment. For example, the head of an elongated fibroblast may be traveling in a different direction than its tail. The orientation of c(x, t) evolves according to d (x, t) dt = f sin ( ) (5) Here is the angle of c, the vector representing the collagen direction and (x, t) is the angle of f (x, t). Thus when the difference between the angles is small the derivative is small and when the directions are 5 -10 0 m 10 -10 0 10 m 0 0.5 1 Figure 1. A graph of the weight function wi. The cell location f i(t) = (0, 0). orthogonal the rate of change reaches a maximum. In addition it is periodic so that when a fiber and cell are oriented in either the same direction or 180 degrees apart the fiber does not change directions. A vector field representation of the fibers is not unique because the fibers are filaments which extend in both the positive and negative direction of its associated vector. The representation of the fibers which we use is cell dependent and satisfies Copyright f i(t) > 0 . (6) Here , denotes the standard Euclidean inner product in the plane and some small positive constant. We choose this representation to ensure that a cell moves along the fiber in the manner which would require the cell to make the smallest change in its direction. In other words, if the vector in the direction of the cell and the vector representing the fiber are placed head to tail, the inner angle, , is obtuse (see Fig. 2). If the cell direction and the fiber direction are at approximately right angles (approximately defined by the choice of ), then the fiber orientation is randomly chosen to be the positive or negative direction. The numerical algorithm consists of the following steps: 1. Interpolate the cell s direction to the grid. 2. Orient the fibers with respect to the cell s direction 3. Interpolate fibers to the cell s location 6 Fiber orientation Cell direction a b Figure 2. A schematic depicting a cell direction and two possible fiber orientations for the same fiber. Of the two possible representations for the fiber orientation (a) and (b), only (b) satisfies equation 6 for the shown cell direction. For this representation the inner angle is greater than 90 degrees which means that when following the fiber, the cell alters its course the least rather than turning back on itself. 4. Normalize the interpolated direction 5. Change the fibers direction and normalize. 6. Move the cell in the direction found in step 3. Further details of the numerical scheme are given in the Appendix. 3 Results of Realignment In this section we describe results from the model of collagen realignment by considering the effect of important parameters, the initial fiber orientation, cell density, how cells enter the domain and finally different boundary conditions. A typical simulation result is shown in Fig. 3. There are 600 cells which start uniformly distributed throughout the domain of 0.5 mm by 1 mm. The boundary conditions are reflective with the component of motion of a cell in the direction perpendicular to the boundary encountered being reflected. This boundary condition was chosen because it is the most realistic for the epidermal boundary in the wound, can be interpreted as a no flux boundary for the other wound edges with an equal number of cells entering and leaving the domain (for further details see the Appendix). The collagen orientation is represented in the figures by drawing curves whose tangent corresponds to the direction of the vector field at the point, i.e. streamlines. 7 a b Figure 3. The collagen orientations and cell positions for a typical simulation. In (a) the initial random collagen orientation is shown and in (b) the collagen orientation is shown after 100 hours of remodelling by the fibroblasts on a domain of 0.5 mm by 1.0 mm. The cells have a speed of 15 microns per hour, = 5, = 0, k = 0.15 hours and the numerical grid for the vector field is 51 by 101. 3.1 Rate at which the fibroblasts alter the collagen: The parameter represents a cell s ability to alter the direction of the collagen fibers and is determined by properties of both the fibroblasts and the collagen matrix. For instance, the stiffness of the collagen matrix, the local tractional forces the cells create, the manner in which the cells produce collagen fibers and the degree the cells chemically alter the matrix will all influence how easily the cells reorient the collagen. Of these four properties mentioned, only the tractional forces of fibroblasts have been measured, giving insufficient data for determining . These forces have been measured by placing cells on floating silicon rubber films and measuring the distortment (Harris et al., 1980) and by measuring the overall contractile forces on collagen gels with moving fibroblasts (Eastwood et al., 1996). Although this information is important, it is not sufficient to determine . Even though quantitative values for are not currently possible, we can check the model s sensitivity to and later we give suggestions of how to roughly determine experimentally. An important result of our numerical simulations is that the maximal smoothness of the extracellular matrix occurs at a finite value of (Fig. 4). That is, the degree by which a randomly oriented vector field is smoothed, due to the cells, increases with to a critical level and then decreases. Large means that the cells quickly reorganize the collagen into the direction they are traveling. The reduced alignment 8 a: = 2.5 b: = 5 c: = 20 Figure 4. The effect of altering the rate at which the fibroblasts change the fiber direction. It is seen that as the influence of the cells on the collagen orientation increases the pattern has more structure in (a) where = 2.5, the most structure in (b) where = 5 (same as Fig. 3[b]) and then less structure in Copyright where = 20. The simulations shown here have the same parameters and setup as that shown in Fig. 3. at large occurs because the cells re-randomize the field and do not create as much alignment. The orientation of the fibrous gel generated by our model can, when started with an initially random configuration, remain random at a global level, yet have a lot of structure (see Figs. 3 and 4). Although globally there is no predominant orientation, on a finer scale different orientations will dominate resulting in local structure and pattern. In order to compare different simulations we want to quantify the amount of structure or alignment, but global averages are of little use. They can not even distinguish between the initial and final configurations. Thus we use the following simplistic method to measure the structure which was inspired by the technique described in Rand & Wilson (1995). We subdivide the domain into uniform, non-intersecting regions. In each region a fiber direction is chosen and the absolute value of the differences with all the other fiber directions within the region are averaged. This gives some measure for the alignment of the fibers within the subregion. The values are then averaged for all the subregions giving one value for the overall domain. A priori, we do not know the scale of the structure, so this procedure is carried out several times with different sized subregions. We consider 2, 4, 8, 16, 32 and finally 64 subregions. The results comparing the simulations with different values of are shown in Fig. 5. Of course, as the regions become too small the measure becomes meaningless since for smaller regions there are fewer points to average and in the extreme case of a region with one grid point, the fiber is truly aligned with itself. One can see by the results that this crude measure confirms what 9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 5 10 15 20 25 30 35 40 45 50 55 Measure of alignment Window size random =2.5 =5 =20 Figure 5. The measure of alignment applied to the results for simulations with different values of . The alignment is greatest for = 5 and the corresponding alignment patterns are shown in Figs. 3 and 4. The x-axis denotes the number of grid points in the square subregions used. The symbols represent the values for in the following manner: (+) the random initial conditions, ( ) = 2.5, (*) = 5 and ( ) = 20. is visually obvious. Throughout the rest of the paper this measure of alignment is used to confirm comparisons which are visually apparent and determine which simulations have more structure when a visual inspection is not sufficient. 3.2 Cell speed: s Increasing the speed of the cell has an obvious effect on the cell paths: they move farther and in longer straight line runs. Because the time step in the simulations is kept constant, the cell samples the environment with the same frequency to determine its direction at higher speeds, but moves more quickly between samplings. Similar results can be obtained by fixing the speed and altering the time step. A more subtle effect of variations in cell speed is on the collagen alignment, and this depends on the support of the weight functions wi. If a cell moves sufficiently fast that the distance traveled in one time step is greater than the support of wi, then at its new location the cell does not continue to alter the fibers at its previous location; we denote this critical value of the speed above by sw. Thus increases in speed do not alter the extent of fiber alignment but simply align the fibers in a different manner. At very low speeds the cell will alter its course many times in a small area causing little overall alignment in the fibrous gel. As the speed increases the cell will alter its course fewer times in a given area but 10 a: s = 5 b: s = 15 c: s = 60 Figure 6. Collagen alignment patterns for different cell speeds. As the cell speed increases the alignment also increases. The speed varies from 5 microns per hour in (a), 15 microns per hour in (b) to 60 microns per hour in Copyright. Other than the speed the setup and parameters are the same as Fig. 3. modify the region within the support of wi causing more alignment. Increasing the speed further will decrease the effect since the cell modifies the fiber at its previous positions to a lesser extent. When it moves faster than sw it no longer modifies the previous position. Thus we expect there to be an optimal cell speed that aligns the gel to the greatest extent. This behavior at low speeds is easily seen in Fig. 6 where an increase from 5 microns per hour to 15 microns per hour clearly increases the alignment of the gel. Yet as the speed increases the boundary conditions begin to significantly affect the patterns making it difficult to confirm whether there is a maximum speed which gives optimal alignment. For our choice of wi and time step, the critical speed sw 132 microns per hour which is out of the range of biologically realistic speeds. Fibroblasts move anywhere from 12 to 60 microns per hour in three dimensional collagen gels (Friedl et al., 1998). For some speed ranges, the simulations are not sensitive to small changes in speed. With the parameter sets tested, we find that for speeds ranging from 15 to 60 microns per hour, small changes in speed have little effect on the overall amount of alignment. 3.3 Cell polarization: Fibroblasts moving on a fibrous substratum tend to become highly elongated with a clear leading edge in what is called a bipolar form (Grinnell & Minter, 1978), as opposed to having a more spread morphology with a large, nonlocalized, leading lamella when placed on a smooth surface (Grinnell & 11 0.1 0.3 0.5 0.7 0.9 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 y position (microns) x position (microns) =0 =0.1 =0.3 =0.7 =0.9 =1.0 Figure 7. A cell path is shown for several simulations where only the cell polarization parameter is varied. The cell path straightens considerably once is non-zero. The various symbols show the cell s location at various times during a simulation. The data shown represents a run of about 33 minutes (k = 0.045 hours). Bennett, 1981). As stated previously the model cells affect the matrix in the support of wi, which realistically should more closely match the area covered by the cell membrane. Yet having a greatly elongated support for wi not only complicates the model but is also an oversimplification of the collagen production and orientation processes. Production occur in involutions within the trailing edge of the cell and not throughout its entire length (Birk & Trelstad, 1986). Thus we take a square support for wi and account for the elongation of fibroblasts by allowing the model cells to be polarized, which means they have a preferential direction; the degree of polarization is indicated in our model by the parameter (see equation 2). As is evident by Fig. 7, this parameter greatly influences the path and final position of the cells. When = 0 the cell turns back on itself and its final position is close to the initial position. As is increased the cell path becomes smoother and the final location of the cell is farther from its initial position, but the most significant effect is during the transition in from zero to non-zero. Even for relatively high values of , such as 0.9, the cell deviates significantly from a straight course, since at each time step there is a slight alteration in direction, and the cumulative effect of this rapidly becomes significant. It is clear from the cell paths that including a non-zero is more realistic since in practice cell types which tend to polarize do not take such circuitous paths (Friedl et al., 1998). Interestingly, although the cell paths are very sensitive to changes from zero to non-zero values of , the collagen alignment pattern is not. In fact, as increases from zero the degree of alignment changes slightly and only as approaches 1 does the effect on the alignment become significant (see Fig. 8). 12 a: = 0 b: = 0.1 c: = 0.7 d: = 0.9 Figure 8. The alignment patterns for simulations which differ only in the amount of cell polarization assumed. Collagen alignment is not sensitive to the polarization parameter until it approaches 1: (a) = 0, (b) = 0.1, Copyright = 0.7, and (d) = 0.9. The parameters and setup are the same as those for Fig. 3. 3.4 Initial fiber orientation An initial aligned region can determine the overall alignment of the domain. We tested the influence of aligned regions by using initial matrix orientations which included such regions but otherwise have the same setup as the previous simulations. In the standard setup although the cells are uniformly distributed in the domain, they are polarized in randomly chosen directions. The initial collagen gel is randomly oriented with a strip running vertically through the domain which has a horizontal orientation. If the width of this strip is sufficiently large it will determine the overall alignment of the region regardless of cell polarization. As the width is decreased the effect on the alignment also decreases until for sufficiently small strips of alignment there is no effect. Cell polarization changes the rate at which the overall alignment decreases as the width of the strip decreases, with the rate being smaller with higher cell polarization. An example of this is shown in Fig. 9. In a particular set of simulations a strip of alignment 20 microns wide did not produce any overall alignment effect with either = 0 or = 0.7. As the width of the strip was increased to 50 microns the case with no cell polarization showed only slightly more horizontal alignment than the normal case, whereas when = 0.7 the horizontal alignment dominated the final fiber configuration. A strip of 100 microns in width was sufficiently large to align the overall final pattern. 13 a: = 0 b: = 0.7 Figure 9. The effect of having an initial aligned region in simulations with and without cell polarization. The alignment patterns are shown for simulations with an initial 30 micron strip located approximately halfway in the horizontal direction and running throughout the domain vertically where the collagen is aligned horizontally. The amount of polarization is the only difference between (a) where = 0 and (b) where = 0.7. Other than the initial orientation and the setup is the same as that for Fig. 3. These results are also influenced by the ability of cells to reorient collagen fibers, as represented by . In some cases which result in global alignment, an increase in the value of alters the result so there is less or no alignment. This is due to the same effect shown in Fig. 4, but here the cells more effectively realign and thus randomize the initially aligned strip. One possible implication of the influence of aligned strips is that very thin scars will diminish with time while the fibroblasts remain active, whereas large scars will become more severe. Thus if the fibroblasts initially setup a small region of alignment, in the months after wounding as they continue to remodel the tissue the alignment should diminish causing less scarring, whereas if the initial aligned region is large enough the remodelling of the fibroblasts will enhance the alignment and cause more severe scarring. This may be a factor in wound healing abnormalities such as keloid scars. In addition, due to the influence of on this process it may be possible to devise experiments which would determine this parameter by seeing if initially aligned regions grow or diminish. 3.5 Cell density As expected, an increase in cell density increases the alignment of the collagen matrix. Biologically this effect should be enhanced since fibroblasts exhibit contact inhibition. In other words, when two cell 14 a b Figure 10. The alignment patterns when cells enter the domain from the boundaries. In (a) = 0.9 and the cells enter uniformly along all the boundaries and in (b) = 0 and the cells enter only at the base. In both cases the cells previous direction is randomly set to be a direction in the half plane pointing into the domain. All other parameters as in Fig. 3. touch each other they retract and go in altered directions. At low densities this is will not change our results, but it has been shown that for high cell densities contact inhibition alone can cause cells to align in parallel arrays due to direct interaction with neighboring cells (Mogilner & Edelstein-Keshet, 1996). This would certainly enhance the alignment of the collagen gel, but is an aspect which we do not include in our model since the density of fibroblasts in the dermis is relatively sparse when compared to those used in the in vitro alignment studies (Erickson, 1978). 3.6 Cell flux The manner in which cells are introduced into the domain plays an important part in determining the alignment of the collagen gel. When the cells enter from the boundary they all have a component of motion in the same direction and tend to reinforce that component leading to more alignment. In the simulations presented thus far, the cells have been initially placed uniformly in the gel. This eliminates the aforementioned alignment effect and also gives a lower local density for the cells. The simulations in Fig. 10 illustrate cases where the cells enter the domain from the boundaries. When the cells are not polarized, = 0, the effect on alignment is small, but with polarization included, the manner in which the cells enter the domain strongly influences the overall alignment patterns. 15 3.7 Boundary conditions With the current domain size and cell speed, the boundary conditions do not significantly affect the results during the time of interest for our applications. However, results of the simulations after long times can be strongly influenced by the boundary conditions because the number of times a cell interacts with the boundary increases with cell speed. We have limited our runs to 100 hours since this is the time of greatest biological significance; moreover experiments dealing with fibroblasts moving in gels typically run on this time scale. In wound healing, experimental data is typically collected up to 21 days post wounding, but it is believed that the important factors determining alignment occur soon after wounding and so we only consider the first 6 days post wounding. Fibroblasts typically begin entering the wound region within 24 to 48 hours post wounding leaving the remaining four to five days as the relevant time. 4 The Tissue Regeneration Model Having acquired a basic understanding of the matrix orientation model, we now develop a more complicated model by adding features such as matrix production. In order to motivate the additional features, we begin by discussing some details of the wound healing process. When the skin is wounded, a sequence of complex and overlapping processes are initiated which result in the repair of the wound. One of the first events is the polymerization of the protein fibrin into a fibrous mesh forming the blood clot (Clark, 1996). This temporary extracellular matrix has several functions including the provision of a framework or scaffolding which the fibroblasts and other cells can use to move into the region. Fibroblasts typically enter the wound region between 24 to 48 hours after wounding. At this time they begin to replace the fibrin clot with an extracellular matrix composed of different proteins, until eventually they construct a matrix primarily of collagen which restores the integrity of the skin. In order to incorporate tissue regeneration the matrix orientation model is modified in three fundamental ways. First the extracellular matrix is considered to be composed of two fibrous networks, namely collagen and fibrin. In this way the domain can simulate normal tissue composed of collagen, a blood clot composed of fibrin or some combination of the two. Second, we extend the model to include the ability of fibroblasts to alter the composition of the matrix by producing and degrading the proteins. Finally the speed of the fibroblasts is determined by the composition of the extracellular matrix. It 16 is well known that fibroblasts move at different speeds on different fibrous substrates. We retain the assumptions from the previous model that the fibroblasts alter the matrix alignment and obtain contact guidance cues from the matrix. In this amended model the fibroblasts are again discrete objects whose paths are given by f i(t). The extracellular matrix is represented by two vector fields, c(x, t) representing the collagen network and b(x, t) representing the fibrin network or blood clot. As before the direction of the vector denotes the predominant direction of the fibers at the specified point in the plane, and the vector representation is cell dependent, satisfying equation 6 for each vector field c and b. Unlike the previous model the vectors do not have unit length; rather, the length of the vectors represents the density of the protein which is altered by the fibroblasts. The fibroblasts receive directional cues and speed information from both components of the extracellular matrix. Thus equations 1 and 2 are modified in the following manner: f i(t) = s(Copyright, b ) v(t) v(t) (7) with v(t) = (1 ) u(f i(t), t) u(f i(t), t) + f i(t ) f i(t ) (8) and u(x, t) = (1 )c(x, t) + b(x, t). (9) As before f represents the time derivative of f , is a positive constant representing the degree to which the cells are polarized, is a time lag, s represents the speed (no longer constant) and is a positive parameter which determines how much influence the fibrin network has on the cell s direction. Blood clots have a high concentration of fibronectin, an adhesive protein on which fibroblasts migrate more readily than on collagen (Grinnell & Bennett, 1981). We assume that the fibronectin is uniformly distributed throughout the blood clot and is proportional to the density of fibrin. Thus the speed of the fibroblasts depends directly on the collagen and fibrin density, although it is really the fibronectin and not the fibrin which is important. The speed function is assumed to be the product of two functions, one which depends on the collagen density and is decreasing and another which depends on the fibrin density and is increasing. The functional form of the collagen dependence is shown in Fig. 11 with the important feature of being relatively insensitive to low collagen densities. The fibrin dependence is taken to be linear for simplicity. Overall the speed can range from a maximum value, denoted by , to a minimum of about 9 (for more details see the Appendix). 17 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 speed factor collagen density Figure 11. The graph of the function showing the dependency of the cells speed on the collagen density. The collagen density is scaled so that it remains between 0 and 1.5. When fibroblasts migrate into the wound region they remove the blood clot by degrading the fibrin and eventually replacing it with a collagen network, via both production and degradation of the extracellular matrix (Jeffrey, 1992). The rate at which they produce collagen is highly dependent on both the composition of the extracellular matrix surrounding them and their chemical environment (Clark et al., 1995; Tuan et al., 1996). For simplicity we model the changes in collagen and fibrin density with the following ordinary differential equations: d c(x, t) dt = (pc dc c(x, t) ) N i=1 wi(x, t) (10) d b(x, t) dt = db b(x, t) N i=1 wi(x, t), (11) where wi is defined in equation 4, N is the total number of cells, pc, dc, and db are all positive constants. The first term on the right-hand side of equation 10 models collagen production at a constant rate by each cell, and the second term models degradation at a rate proportional to the density of collagen already present. The contribution of a cell to the collagen modification is given by multiplying the two terms by wi, then each cell s contribution is added up to give the total production rate. As before wi interpolates the influence of the discrete cells to the nearby continuum fiber variables. Equation 11 for the evolution of the fibrin density is similar to equation 10 for the collagen density evolution but the cells do not produce fibrin so there is no production term. The extracellular matrix is reoriented in the same manner as before and is described by equations 5 and 3. This completes the description of the tissue regeneration model. 18 5 Results of the Tissue Regeneration Model In the following simulations our standard setup has initial conditions where at any point there is only one matrix type. Thus the domain is divided into non overlapping regions with either collagen or fibrin networks. The fibrin is randomly oriented, representing a blood clot, and is devoid of fibroblasts which enter from the domain edges. There are two basic cases: one which can be used to compare with the results of the previous model and one which is the standard setup for this amended model. The first basic case has an initial randomly oriented fibrin matrix of uniform density in the domain. The cells are initially uniformly placed within the blood clot (not relevant to wound healing) and the results are shown in Fig. 12. One can see that the collagen density is fairly uniform in its distribution and the fibrin has been effectively removed. By comparing the results to Fig. 8(d) it is obvious that there is much less alignment with the new model. This is primarily due to the changes in the cell speed, which in this case gives an average cell speed of 3.37, and can be demonstrated by fixing the cell speed at 15 microns per hour and comparing the results (see Fig. 13). The second base case is the same as the first except the cells enter the domain from the boundaries with the results shown in Fig. 14. In this more biologically realistic case, the effects of the cells entering from the boundary on the extracellular matrix can be clearly seen at the earlier time. 5.1 Alterations in collagen production and fibrin degradation: pc, dc and db As the collagen production is increased the average collagen density increases and the fibrin density is unaffected. This is demonstrated in Figs. 15(a) and Copyright, which show plots of the average collagen and fibrin densities with respect to the parameter pc. Figs. 15(b) and (d) show time plots of the average collagen and fibrin protein densities respectively. Each line corresponds to the time plot for one simulation. Thus the time plots for several simulations with different values of pc are plotted on the same graph. In this way the range of densities obtained as the parameters are changed is more easily seen and surprisingly there are significant fluctuations. Likewise when db is increased the fibrin density decreases and the collagen density remains roughly the same (see Fig. 16) but again with fluctuations. These fluctuations are explained by the discrete nature of the cells. When parameter values are changed, the cells traverse different paths giving solutions which pointwise in space are very different. Fig. 17 shows the average collagen density for several simulations which differ only in the random initial conditions used for the matrix orientation. The fluctuations in collagen densities in Figs. 16 and 17 as well as in the fibrin density in Fig. 15 are of the same magnitude confirming that they are independent of the 19 a b Copyright Figure 12. Results from a typical simulation where the cells are initially uniformly placed in the domain. In (a) the collagen orientation and density are shown. The lines representing the collagen orientation are streamlines for the vector field and are grey scaled to represent the density with white representing no collagen and black representing a collagen density of 1.5. In (b) the fibrin density is shown with black representing fibrin density of 1.0 and white representing no fibrin. In Copyright the cell positions are shown with the area in black corresponding to the support of the weight functions wi. The average cell speed is 3.37 microns per hour and the average collagen density at t = 96 is 1.30. The time shown is t = 100 hours, k = 0.15 hours, N = 600, = 15 microns per hour, = 0.5, = 5, pc = 0.64, dc = 0.44, df = 0.6, and = 0.9. 20 Figure 13. The collagen alignment when the cell speed is fixed. This simulation is the same as that shown in Fig. 12 except the cell speed is fixed at 15 microns per hour and the average collagen density at t = 96 is 1.44. parameter being altered. The other effects of altering these parameters are predominantly due to changes in the cell speed and can also be obtained by altering the functional dependence of the speed on the densities. This includes the degree the cells are able to invade the blood clot which is discussed in the next section. 5.2 Invasiveness The different setup for the tissue regeneration model introduces the additional feature of invasiveness which needs to be considered. Since the blood clot is initially devoid of fibroblasts, the degree to which the fibroblasts in the collagen matrix are able to move into the fibrin becomes an important issue. There are three factors which play a role in determining the invasiveness of the cells: cell polarization, cell speed and protein production. The protein production plays a part since it alters the cell s speed by altering the densities. By changing either the speed function or the production and degradation parameters any degree of invasiveness can be obtained. Cell polarization plays a role by determining the frequency with which a cell turns around rather than penetrating into a region. 21 a b c d e f Figure 14. Results from a typical simulation where the cells enter the fibrin matrix from the domain edges. In (a), (b), and Copyright the collagen orientation and density, the fibrin density and the cell positions, respectively, are shown at t = 20 hours. In (d), (e) and (f) the collagen orientation and density, the fibrin density and the cell positions, respectively, are shown at t = 100 hours. Compare with Fig. 10(a). The parameters are the same as those shown in Fig. 12. The average cell speed is 3.24 microns per hour and the average collagen density at t = 96 is 1.25. 22 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Average collagen density Rate of collagen production t=20 t=50 t=80 t=140 t=200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 50 100 150 200 Collagen density Time a b 0 0.1 0.2 0.3 0.4 0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Average fibrin density Rate of collagen production t=20 t=50 t=80 t=140 t=200 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 Fibrin density Time c d Figure 15. Graphs of the average protein densities as the collagen production is varied. It is seen that the collagen density increases as the collagen production increases and the fibrin density remains relatively unchanged. In (a) and Copyright the average collagen and fibrin densities are plotted with respect to the collagen production rate respectively at several different times. The collagen density is scaled to range from 0 to 1.5 and fibrin density of 1 represents a typical clot density. The time plots (b) and (d) show the evolution of the collagen and fibrin densities respectively with each line representing a different simulation as the collagen production rate is changed. Other than varying pc and dc while keeping their ratio a constant 1.5, the parameters and setup is the same as that for Fig. 14. The values of pc used in the simulations start at 0.1 and are incremented by about .05 up to 0.85. 23 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Average collagen density Rate of fibrin degradation t=20 t=50 t=80 t=140 t=200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 50 100 150 200 Collagen density Time a b 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Average fibrin density Rate of fibrin degradation t=20 t=50 t=80 t=200 t=140 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 50 100 150 200 Fibrin density Time c d Figure 16. Graphs of the average protein densities as the fibrin degradation is varied. The fibrin density decreases as the degradation parameter increases and the collagen density remains relatively unchanged. In (a) and Copyright the average collagen and fibrin densities are plotted with respect to the fibrin degradation rate respectively at several different times. The collagen density is scale to range from 0 to 1.5 and fibrin density of 1 represents a typical clot density. The time plots (b) and (d) show the evolution of the collagen and fibrin densities respectively with each line representing a different simulation as the fibrin degradation rate is changed. As before other than varying db the parameters and setup are the same as Fig. 14. The values of db used in the simulations start at 0.1 and are incremented by about .05 up to 0.85. 24 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 -300 -250 -200 -150 -100 -50 0 Average collagen density Seed for random initial orientation t=20 t=50 t=80 t=140 t=200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 50 100 150 200 Collagen density Time a b Figure 17. The effect of changing the random initial orientations on the average collagen density. In (a) the collagen density is plotted with respect to the seed for the random number generator at several different times. The collagen density is scale to range from 0 to 1.5. In (b) time plots of the collagen density in simulations using different seeds for the random number generator determining the initial fibrin orientation are superimposed. The parameters and setup are the same as Fig. 14 other than the particular random initial conditions. Fifteen different seeds were used. 5.3 Influence of the clot on cell direction: The fibrin clot has a random structure created by the polymerization of the fibrin as the blood plasma fills the wound. The directional cues of the fibrin clot have a randomizing effect on the direction of the cells. In equation 9 the influence of the fibrin on fibroblast direction is determined by the parameter . When is small, so that the fibrin does not contribute much to the cell direction, there is more alignment of the collagen fibers and as increases the collagen becomes more random. An example of this can be seen in Fig. 18 where the value of is varied for two simulations. Even when is small the randomizing influence will, with time, have a significant impact on the cells directions and will consequently alter the overall collagen alignment. This indicates that either the mechanisms creating alignment in scar tissue are strong or the fibrin clot plays a minor role in directing the fibroblasts. 5.4 Interface The different way in which the cells interact with the two types of protein fibers causes different alignment properties. Most significantly, the collagen fibers direct the cells and are reoriented by them, whereas the fibrin fibers are not reoriented. This causes more alignment when the cells are in collageneous regions of the domain and the randomizing influence of the fibrin is not reduced as the cells flow through those regions. This is demonstrated by simulations in which the domain is initially divided into a random collagen gel and a random fibrin gel, with a distinct interface between the two in the 25 a: = 0.1 b: = 0.9 Figure 18. Changes in how the blood clot influences cell direction are shown in these collagen plots. When the fibrin matrix has less influence on the cell direction (a) the collagen is more aligned and when the fibrin has more influence (b) the collagen is less aligned. Compare (a) and (b) with Fig. 14 where = 0.1, = 0.9 and = 0.5 respectively. The average cell speed in (a) is 3.24 and the average collagen density is 1.23. In (b) the average cell speed is 3.25 and the average collagen density 1.13. vertical direction. Fig. 19 shows the collagen alignment patterns for such a simulation with 300 cells initially positioned at the collagen-fibrin interface. There is more alignment in the collagen region with the cells being guided into aligned streams, whereas on the fibrin side a less oriented collagen matrix is produced. By considering only the speed effect, one would expect opposite results since the cells are moving more slowly in the collagen region than in the fibrin region. By eliminating the speed dependence and making the cells travel at a constant speed, which is close to the average speed in the previous simulation, similar results are obtained (Fig. 19(b)). This demonstrates that speed is not the controlling factor. However, if the speed is a constant 15 microns per hour strong alignment is attained on both sides showing how different effects combine to influence the outcome. The higher speed tends to increase the alignment and, for this example, the effect of the high speed overcame the randomizing influence of the fibrin clot. At the lower speed the randomizing influence of the fibrin clot determines the collagen alignment on the fibrin side, but on the other region the alignment is due to the mechanism relating to how the cells interact with the fibers and not the speed. 26 a b Figure 19. Alignment patterns showing the effect of an initial transition from collagen to fibrin. The alignment is stronger on the left side of the domain where the initial protein was collagen and weaker on the right side where the initial protein was fibrin. The position of the initial interface between collagen and fibrin is indicated by the arrow. In (a) the cell speed depends on the protein densities with an average speed of 3.11 microns per hour, = 15 and the total collagen density is 1.14. In (b) the cell speed is constant and set at 5 microns per hour. The remaining parameters are the same as those given in Fig. 12 with 300 cells initially placed along the collagen fibrin interface and polarized in randomly chosen directions. 27 6 Discussion We have developed and studied a matrix orientation model which includes the interactions of cells being directed by their substrate and simultaneously reorienting it. Furthermore we have extended this model to allow the cells to alter the composition as well as the orientation of the matrix. Simulations of these models have important implications concerning the mechanisms which cause alignment. The simple model clearly shows that cell speed, flux, polarization, density, initial matrix orientation and the influence of cells on the matrix can all have a significant impact on the overall alignment of the collagen. As was indicated by both models, depending on the circumstances, some of these factors will contribute more than others in producing alignment. Most, if not all, of these can be tested experimentally, but the complicated interactions mean that the results require careful interpretation; also the effects of gel contraction should be minimized as far as possible in experimental tests. We suggest four types of possible experiments: Altering the speed of fibroblasts. Our model predicts that increasing the speed of the cells should cause greater alignment. This could be done either with a chemoattractant (Knapp et al., 1999) or by altering the integrin expression levels of the fibroblasts (Palecek et al., 1997). By choosing a chemoattractant such as epidermal growth factor which can increase the speed of the fibroblasts up to 3-fold depending on the substratum (Ware et al., 1998) the collagen alignment should be increased. Of course this chemoattractant also decreases the directional persistence of cells modeled here as polarization. These two properties of the cells have opposite effects on the collagen alignment, but our model predicts that the effect of the increased speed outweights the effect of the decreased cell polarization. For example, increasing the cell speed by a factor of three and halving the polarization parameter results in a significantly more aligned collagen matrix. Reducing the contact guidance of the cell. This can be accomplished by inhibiting the formation of microtubules with colcemid (Oakley et al., 1997). When this is done the cells are still mobile but do not align along very narrow grooves in the substratum. Alternately, treatment with colchicine causes a rounder morphology of fibroblasts in a fibrous gel but does not significantly alter the collagen production (Unemori & Werb, 1986). A final possibility would be to alter the substrate by using monomeric collagen (Mercier et al., 1996). On this type of substrate the fibroblasts take a more rounded morphology and presumably have less contact guidance. All of these interventions alter not only the contact guidance but also the degree of cell polarization. Our model suggests both of these act to reduce alignment. Reducing cell polarization is modeled by reducing , and 28 although in our model we assume that the cells align exactly in the predominant direction of the collagen, reducing should effectively model a reduction in contact guidance. Effects of initial collagen orientation. This could be tested using the approach of Matsumoto et al. (1998), who transplanted pieces of tendon where the fibroblasts had been killed into normal tendon. The manner in which the structure of the transplanted tendon was altered by the invading cells was examined in two cases: one where the transplant was the same size as the defect into which it was transplanted and another where the transplant was larger and therefore lax. Altering this procedure by rotating the grafted tendon, one can setup initial conditions where the collagen has regions with differing orientations. The results of two simulations with this type of initial conditions are shown in Fig. 20. The initial collagen gel is oriented vertically with a square patch in the center oriented horizontally. The different values of show very different results. A similar idea would be to place pieces of oriented gel (Guido & Tranquillo, 1993) in randomly oriented gel. By changing the size of the differently oriented regions the magnitude of could be estimated, helping to determine how effectively the cells can reorient collagen. Determining the effects of flux. There are many experiments which could be devised to determine the importance of the way in which the cells enter the region. Gels can be prepared void of cells and with cell distributed throughout. Our simulations suggest that gels which have a continuous flux of cells from the edges will result in more aligned collagen near the gel edges than gels which have cells distributed throughout. The tissue regeneration model showed that the rate at which the cells alter the matrix composition has little direct effect on the matrix alignment. In the more complicated setting of tissue regeneration there was less alignment due to both the reduced cell speed and the randomizing influence of the blood clot. Thus our modelling shows that matrix alignment is a complicated process with many factors contributing to the overall structure including both the inherent properties of the fibroblasts as well as the conditions in which they find themselves. Even in our simple model the complexity of alignment was illustrated. In the context of wound healing we conclude that cell flux is the most significant alignment mechanism, especially when coupled with cell polarization. Yet, the initial orientation of the collagen should not be overlooked. As the cells move into a wound they may establish an initial region of alignment when regenerating the tissue which is then propagated or stabilized during the months of remodelling that occur. Thus inherent properties of the fibroblasts are of secondary importance to the conditions in 29 a: = 0.5 b: = 5.0 Figure 20. Simulations designed to show the effect of an initial region of collagen with different orientation. Initially all the collagen is oriented vertically except a 200 by 200 micron square in the center, which is oriented horizontally. In (a) = 0.5 and the region of horizontally oriented collagen remains largely unchanged after 100 hours of simulation; whereas, in (b) = 5.0 and the region with horizontal orientation has changed shape dramatically. In these simulations = 0.7 and the remaining parameters are the same as those given in Fig. 3. which the cells are placed in determining matrix alignment for wound healing. In addition the alignment mechanism must be strong to overcome randomizing influences such as the blood clot. In a complex biological process such as wound repair, mathematical modelling plays a valuable role because of its ability to study specific parts of the process in isolation. This approach has been used to study epithelial repair (Dale et al., 1994), cytokine activity in scar tissue (Dale et al., 1996), wound angiogenesis (Pettet et al., 1996; Chaplain & Byrne, 1996), and the tissue mechanics underlying wound contraction (Olsen et al., 1995; Tracqui et al., 1995; Tranquillo & Murray, 1992). Almost all of this modelling work excludes the effects of tissue anisotropy, which will vary dynamically during healing as a result of cell matrix interactions. However, experimental evidence now suggests that such anisotropy plays an important role in scar formation (Ehrlich & Krummel, 1996; Whitby & Ferguson, 1991). The models we have developed here provide a framework that enables realistic modelling of these phenomena, including dynamically varying anisotropy. This is particularly important because of recent progress in the development of potential anti-scarring therapies (Shah et al., 1995), which act by altering the growth factor profile within the wound during healing. A realistic mathematical model of the scar formation process would provide a powerful tool in this development, enabling the identification 30 of optimal therapeutic regimes, and our work provides one component of such a realistic model. Acknowledgements. This research (PKM) was supported in part by the Institute for Mathematics and its Applications (University of Minnesota) with funds provided by the National Science Foundation. Part of this work was supported by the London Mathematical Society, under scheme 3. JCD acknowledges support from EPSRC grant GR/K71394 7 Appendix Here we describe the numerical implementation of the models and their details. Before doing so we note that for simple configurations of the extracellular matrix such as oriented in one direction, local realignment only causes local changes which do not propagate into the entire domain. That is if the initial orientation is smooth on a scale of sk where s is the cell speed and k is the time step, it remains smooth in a stable manner such that small perturbations are smoothed. 7.1 Protein Densities The vector fields c and b are discretized in space so that Copyright,m(t) = c( hx,mhy, t) and similarly for b ,m where hx and hy are the space steps in the x and y direction. Equations 10 and 11 describing the evolution of the protein densities are solved at the grid locations using Euler s method. The densities are constrained to be non-negative and the grid used has 51 by 101 grid points. 7.2 Protein Fiber Directions Equation 3 is implemented with the time lag = 0.15 hours. This represents a biologically reasonable lag of about 9 minutes. It also corresponds to the time step k for the calculations making the time lagged direction readily available. Multiplying equation 5 by cos , applying the change rule and trigonometric identities give d dt (sin ) = f sin sin sin2 cos sin cos , (12) and similarly multiplying equation 5 by sin gives d dt (cos ) = f cos sin sin cos cos cos2 . (13) 31 By changing to Cartesian coordinates where c = c Copyright , c = (cos , sin ), f = f (cos , sin ) and f is defined in equation 3, equations 12 and 13 become d dt ( Copyright + f Copyright f c = f . (14) Discretizing equation 14 using a first order approximation to the derivative and adding terms of O(k) and higher results in Copyright,m((n + 1)k) = d(nk) + kf (nk) d(nk) + kf (nk) + O(k) (15) where Copyright,m is the discretized version of c and d(t) = Copyright,m(t) (16) with the sign being chosen so that d, f > 0 . (17) Recall that the vector representation is cell dependent (see equation 6), but the cumulative effect of all the cells is being used in equation 15. Thus the representation needs to be defined with respect to the cumulative direction of all the contributing cells. This formulation is valid for both the matrix orientation and the tissue regeneration models and is used for pragmatic reasons in the numerical algorithm. When the collagen density is allowed to vary the cells will remodel the existing collagen to the same degree regardless of the density. If there is no collagen then d is taken to be zero. 7.3 Interpolation The interpolation from the discrete cells to the vector field is a simple evaluation of the weight function defined in equation 4. Other weight functions have been tested including a step function which is constant over the same support of the one defined earlier. The same qualitative results are obtained. The interpolator I from the discretized vector fields to the cell location is a tensor product interpolant using quartic Lagrangian interpolation in each direction (Ralston & Rabinowitz, 1978). It is defined by I( Copyright x, y) = 2 n= 2 2 m= 2 cj m,k n j m(x) k n(y) (18) j(x) = (x xj 2)(x xj 1)(x xj+1)(x xj+2) (xj xj 2)(xj xj 1)(xj xj+1)(xj xj+2) (19) where xj = (j 1)hx and the analogous yk = (k 1)hy are grid points. 32 7.4 Cell Paths The paths of the fibroblasts, f i(t) are computed as continuous piecewise linear curves. We define the direction of the linear segment starting at f i(nk) to be Ti n where n = 0, 1, 2 . . . K and k is the time step. This direction is determined by implementing the equations 7, 8 and 9 as follows: Ti n+1 = vn vn (20) with vn = (1 )un(f i[(n + 1)k]) + Ti n Ti n (21) and un(x) = (1 )I( c(nk), x) + I( b(nk), x) (22) with similar equations for equations 1 and 2 of the matrix orientation model. The cell takes its directional cues from the matrix before remodelling it, although it is the altered density which the cell uses to determine its speed. The speed is constant in the matrix orientation model and for the tissue regeneration model is taken to be sn = 1 8 + 7 8 + 3 I[ c(nk, f i(nk))] 6 1 + 2 I[ b(nk, f i(nk))] 3 (23) where is a positive constant denoting the maximum speed of the cell and sn is the speed corresponding to the line segment starting at f i(nk) . The speed and k determine the length of the line segments which define f i. The speed function is the product of an increasing linear function depending on the fibrin density and a decreasing function of collagen density with the form illustrated in Fig. 11. It is relatively insensitive to low collagen densities and decreases the speed by about 1/4 when the collagen density is 1 and up to about 2/3 when the collagen density is 1.5 (the maximum allowed). The dependency of the speed on the fibrin is a linear function which goes from 1/3 when the fibrin density is 0 to 1 when the fibrin density is 1. Thus the speed can range from a maximum when b = 1 and Copyright = 0 to a minimum of about 9 when b = 0 and Copyright = 1.5. 7.5 Boundary Conditions In wounds the cells potentially can enter or leave the region from all of the wound boundaries with the exception of the epidermal boundary. This is clearly not the case for controlled experiments of cells in collagen gels. It is not understood how the cells respond to the boundaries in the wound, and in 33 controlled experiments the focus is on the interior of the gel and not the edges. We have tried various boundary conditions including using a toroidal surface, absorbing boundaries and two types of no flux boundaries. The two no flux boundary conditions are first one where the cells re-enter the domain at randomly chosen sites after leaving it. The other one is used in the simulations for this paper and has the cell reflect off the boundary by changing the sign of the component of motion perpendicular to the boundary. This boundary condition was chosen because it is the most realistic for the epidermal boundary in the wound, can be interpreted as a no flux boundary for the other wound edges with an equal number of cells entering and leaving the domain. Furthermore, it simplifies the problem by having the same boundary condition for all the edges. In addition the simulation results are very similar for either no flux boundary condition and on the time scale of 100 hours, they are similar for all the boundary conditions examined. 7.6 Parameter Values Although there is a significant amount of experimental data dealing with fibroblasts and collagen, establishing parameter values is very difficult due to the nature of the experiments. Here we justify, to the extent possible, the parameter values used in the model. In the matrix orientation model the speed of the fibroblasts is generally taken to be 15 microns per hour. Fibroblasts moving in three dimensional collagen gels move from 12 to 60 microns per hour in vitro (Friedl et al., 1998), but fibroblasts migrating through a fibrin clot while remodelling the composition of the extracellular matrix certainly move with speeds in the lower range if not slower. In the tissue regeneration model the maximum cell speed, is taken to be 15 microns per hour with typical average cell speeds over 3 microns per hour. The effects of varying this parameter are discussed in Section 3.2. Total cell number N is taken to be 600. This number gives the maximum fibroblast density in a wound region (Adams, 1997). The cell polarization parameter is discussed in Section 3.3. The parameter representing the influence of the fibroblasts on orienting the collagen is . For = 1, a cell would reorient the collagen at its location to its direction within 4-5 hours, as can be demonstrated by solving equation 5 with f = 1. Detailed estimates of this parameter are unavailable, and effects of this parameter are discussed in Section 3.1. 34 The extent to which the fibrin clot guides the fibroblasts is governed by the parameter . We take = 0.5 assuming the contact guidance influence on the fibroblasts to be equal for the collagen and fibrin, although numerical simulations indicate that only near the extreme values of = 0 and = 1 are the results of the simulation significantly affected (see Section 5.3). The time lag can be thought of as the time difference from when the cell senses the contact guidance cues to the time it take the cell doing the remodelling to reorient. It is taken to be 0.15 hours (9 minutes) which means that the cell moves at most 2.25 microns in this time interval. The collagen production and degradation rates, pc and dc as well as the fibrin degradation rate, db are not easy quantities to determine. The experiments which measure the collagen production of fibroblasts in wounds collect data several days apart. Thus with the range of errors for the measurements and only one data point other than zero collagen density at the time of wounding, there was little use in trying to fit the data. The effect of altering the parameters is shown in Section 5.1. There the values range from 0.1 to 0.9 which in the spatially homogeneous case give transition times from the initial conditions to steady state ranging between 5 hours and 60 hours. References Adams, J. J. (1997). The cell kinetics of murine incisional wound healing. Ph.D. thesis, University of Manchester. Barocas, V. H., & Tranquillo, R. T. (1997). An anisotropic biphasic theory of tissue-equivalent mechanics: the interplay among cell traction, fibrillar network deformation , fibril alignment and cell contact guidance. AMSE J. Biomech. Eng., 119(May), 137 145. Birk, D. E., & Trelstad, R. L. (1986). Extracellular compartments in tendon morphogenesis: Collagen, fibril, bundle, and macroaggregate formation. J. Cell Biol., 103, 231 240. Chaplain, M. A. J., & Byrne, H. M. (1996). Mathematical modelling of wound healing and tumour growth - 2 sides of the same coin. Wounds: A Compendium of Clinical Research and Practice, 8, 587 616. Clark, P., Connolly, P., Curtis, A. S. G., & Wilkinson,CopyrightCopyright W. (1990). Topographical control of cell behaviour. ii. multiple grooved substrata. Dev., 108, 635 644. Clark, R. A. F. (1996). Wound repair overview and general considerations. Pages 3 50 of: Clark, R. A. F. (ed), The molecular and cellular biology of wound repair, 2 edn. Plenum Press. 35 Clark, R. A. F., Nielsen, L. D., Welch, M. P., & McPherson, J. M. (1995). Collagen matrices attenuate the collagen-synthetic response of cultured fibroblasts to TGF- . J. Cell Sci., 108, 1251 1261. Dale, P., Maini, P. K., & Sherratt, J. A. (1994). Mathematical modeling of corneal epithelial wound healing. Math. Biosci., 124, 127 147. Dale, P. D., Sherratt, J. A., & Maini, P. K. (1996). A mathematical model for collagen fibre formation during foetal and adult dermal wound healing. Proc. R . Soc. Lond. B, 263, 653 660. Dallon, J.Copyright (2000). Numerical aspects of discrete and continuum hybrid models in cell biology. Applied Numerical Mathematics, 32, 137 159. Dallon, J.Copyright, & Sherratt, J. A. (1998). A mathematical model for fibroblast and collagen orientation. Bull. Math. Biol., 60(1), 101 129. Eastwood, M., Porter, R., Khan, U., McGrouther, G., & Brown, R. (1996). Quantitative analysis of collagen gel contractile forces generated by dermal fibroblasts and the relationship to cell morphology. J. Cell. Physiol., 166(1), 33 42. Edelstein-Keshet, L., & Ermentrout, B. G. (1990). Models for contact-mediated pattern formation: cells that form parallel arrays. J. Math. Biol., 29, 33 58. Ehrlich, P. H., & Krummel, T. M. (1996). Regulation of wound healing from a connective tissue perspective. Wound Rep. Reg., 4(April-June), 203 210. Erickson,Copyright A. (1978). Analysis of the formation of parallel arrays by bhk cells in vitro. Exp. Cell Res., 115, 303 315. Friedl, P., Z anker, K. S., & Br ocker, E. B. (1998). Cell migration strategies in 3-d extracellular matrix: Differences in morphology, cell matrix interactions, and integrin function. Microsc. Res. Tech., 43, 369 378. Grinnell, F., & Bennett, M. H. (1981). Fibroblast adhesion on collagen substrata in the presence and absence of plasma fibronectin. J. Cell Sci., 48, 19 34. Grinnell, F., & Minter, D. (1978). Attachment and spreading of baby hamster kidney cells to collagen substrata: Effects of cold-insoluble globulin. Proc. Nat. Acad. Sci. USA, 75(9), 4408 4412. Guido, S., & Tranquillo, R. T. (1993). A methodology for the systematic and quantitative study of cell contact guidance in oriented collagen gels. J. Cell Sci., 105, 317 331. Harris, A. K., Wild, P., & Stopak, D. (1980). Silicone rubber substrata: A new wrinkle in the study of cell locomotion. Science, 208, 177 179. 36 Harris, A. K., Stopak, D., & Wild, P. (1981). Fibroblast traction as a mechanism for collagen morphogenesis. Nature, 290, 249 251. Hay, E. D. (ed). (1991). Cell biology of extacellular matrix. 2 edn. Plenum Press. Hsieh, P., & Chen, L. B. (1983). Behavior of cells seeded on isolated fibronectin matrices. J. Cell Biol., 96, 1208 1217. Jeffrey, J. J. (1992). Collagen degradation. Chap. 10, pages 177 194 of: Cohen, I. K., Diegelmann, R. F., & Lindblad, W. J. (eds), Wound healing biochemical and clinical aspects. Saunders. Knapp, D. M., Helou, E. F., & Tranquillo, R. T. (1999). A fibrin or collagen gel assay for tissue cell chemotaxis: Assessment of fibroblast chemotaxis to grgdsp. Exp. Cell Res., 247, 543 553. Markwald, R., Fitzharris, T., Bolender, D., & Bernanke, D. (1979). Structural analysis of cell:matrix association during the morphogenesis of atrioventricular cushion tissue. Dev. Biol., 69, 634 654. Matsumoto, N., Horibe, S., Nakamura, N., Senda, T., Shino, K., & Ochi, T. (1998). Effect of alignment of the transplanted graft extracellular matrix on cellular repopulation and newly synthesized collagen. Arch. Orhtop. Trauma Surg., 117, 215 221. Mercier, I., Lechaire, J.-P., Desmouliere, A., Gaill, F., & Aumailley, M. (1996). Interactions of human skin fibroblasts with monomeric or fibrillar collagens induce different organization of the cytoskeleton. Exp. Cell Res., 225, 245 256. Mogilner, A., & Edelstein-Keshet, L. (1995). Selecting a common direction i. how orientational order can arise from simple contact responses between interacting cells. J. Math. Biol., 33, 619 660. Mogilner, A., & Edelstein-Keshet, L. (1996). Spatio-angular order in populations of self-aligning objects: formation of oriented patches. Physica D, 89, 346 367. Mogilner, A., Edelstein-Keshet, L., & Ermentrout, B. G. (1996). Selecting a common direction ii. peak-like solutions representing total alignment of cell clusters. J. Math. Biol., 34, 811 842. N ubler-Jung, K. (1987). Tissue polarity in an insect segment: denticle patterns resemble spontaneously forming fibroblast patterns. Dev., 100, 171 177. Oakley,Copyright, Jaeger, N. A. F., & Brunette, D. M. (1997). Sensitivity of fibroblasts and their cytoskeletons to substratum topographies: Topographic guidance and topographic compensation by micromachined grooves of different dimensions. Exp. Cell Res., 234, 413 424. Okubo, A. (1986). Dynamical aspects of animal grouping: swarms, schools flocks, and herds. Adv. Biophys., 22, 1 94. 37 Olsen, L., Sherratt, J. A., & Maini, P. K. (1995). A mechanochemical model for adult dermal wound contraction and the permanence of the contracted tissue displacement profile. J. Theor. Biol., 177, 113 128. Olsen, L., Sherratt, J. A., Maini, P. K., & Marchant, B. (1998). Simple modelling of extracellular matrix alignment in dermal wound healing I. cell flux induced alignment. J. Theor. Med., 1, 172 192. Olsen, L., Maini, P. K., Sherratt, J. A., & Dallon, J.Copyright (1999). Mathematical modelling of anisotropy in fibrous connective tissue. Math. Biosci., 158, 145 170. Palecek, S. P., Loftus, J.Copyright, Ginsberg, Mark, H., Lauffenburger, D. A., & Horwitz, A. F. (1997). Integrinligand binding properties govern cell migration speed through cell-substratum adhesiveness. Nature, 385, 537 540. Pettet, G. J., Chaplain, M. A. J., McElwain, D. L. S., , & Byrne, H. M. (1996). On the role of angiogenesis in wound healing. Proc. R . Soc. Lond. B, 263, 1487 1493. Priestly, E. B., Wojtowiez, P. J., & Sheng, P. (eds). (1975). Introduction to liquid crystals. New York: Plenum Press. Ralston, A., & Rabinowitz, P. (1978). A first course in numerical analysis. New York: McGraw Hill Book Company. Rand, D., & Wilson, H. (1995). Using spatio-temporal chaos and intermediate-scale determinism to quantify spatially-extended ecosystems. Proc. R . Soc. Lond. B, 259, 111 117. Shah, M., Foreman, D. M., & Ferguson, M. W. J. (1995). Neutralisation of TGF- 1 and TGF- 2 or exogenous addition of TGF- 3 to cutaneous rat wounds reduce scarring. J. Cell Sci., 108, 985 1002. Small, J. V., Herzog, M., & Anderson, K. (1995). Actin filament organization in the fish keratocyte lamellipodium. J. Cell Biol., 129, 1275 1286. Stopak, D., & Harris, A. K. (1982). Connective tissue morphogenesis by fibroblast traction. Dev. Biol., 90, 383 398. Tracqui, P., Woodward, D. E., Cruywagen, G.Copyright, Cook, J., & Murray, J. D. (1995). A mechanical model for fibroblast-driven wound healing. J. Biol. Systems, 3, 1075 1084. Tranquillo, R. T., & Murray, J. D. (1992). Continuum model of fibroblast-driven wound contraction: Inflammation-mediation. J. Theor. Biol., 158, 135 172. 38 Trelstad, R. L., & Hayashi, K. (1979). Tendon collagen fibrillogenesis: Intracellular subassemblies and cell surface changes associated with fibril growth. Dev. Biol., 71, 228 242. Tuan, T.-L., Song, A., Chang, S., Younai, S., & Nimni, M. E. (1996). In vitro fibroplasia: Matrix contraction, cell growth, and collagen production of fibroblasts cultured in fibrin gels. Exp. Cell Res., 223, 127 134. Unemori, E. N., &Werb, Z. (1986). Reorganization of polymerized actin: A possible trigger for induction of procellagenase in fibroblasts cultured in and on collagen gels. J. Cell Biol., 103, 1021 1031. Ware, M. F., Wells, A., & Lauffenburger, D. A. (1998). Epidermal growth factor alters fibroblast migration speed and directional persistence reciprocally and in a matrix-dependent manner. J. Cell Sci., 111, 2423 2432. Whitby, D. J., & Ferguson, M. W. J. (1991). Immunohistochemical localization of growth factors in fetal wound healing. Dev. Biol., 147, 207 215. 39