A Mathematical Model for Fibroblast and Collagen Orientation

Published in Bulletin of Mathematical Biology (1998)60, 101 129 A Mathematical Model for A Mathematical Model for Fibroblast and Collagen Orientation
Fibroblast and Collagen Orientation John C.
Dallon and Jonathan A. Sherratt Mathematics Institute University of Warwick Coventry CV4 7AL, U.K. February 15, 2000 RUNNING TITLE: An orientation model 1 Abstract Due to the increasing importance of the extracellular matrix in many biological problems, in this paper we develop a model for fibroblast and collagen orientation with the ultimate objective of understanding howfibroblasts form and remodel the extracellular matrix, in particular its collagen component. The model uses integro-differential equations to describe the interaction between the cells and fibers at a point in space with various orientations. The equations are studied both analytically and numerically to discover different types of solutions and their behavior. In particular we examine solutions where all the fibroblasts and collagen have discrete orientations, a localized continuum of orientations and a continuous distribution of orientations with several maxima. The effect of altering the parameters in the system is explored, including the angular diffusion coefficient for the fibroblasts, as well as the strength and range of the interaction between fibroblasts and collagen. We find the initial conditions and the range of influence between the collagen and the fibroblasts are the two factors which determine the behavior of the solutions. The implications of this for wound healing and cancer are discussed, including the conclusion that the major factor in determining the degree of scarring is the initial deposition of collagen. 1 Introduction Alignment within biological systems has been the subject of considerable recent interest. Applications vary over a wide variety of systems including herd movement, flocks of birds, schools of fish, insect swarms (Grunbaum, 1994; Okuba, 1986), cellular movement, actin networks and collagen networks (Elsdale, 1973; Pollard & Cooper, 1986; Besseau & Giraud-Guille, 1995). This paper focuses on cellular alignment with respect to collagen fibers. There are numerous applications of this specific system including tumor growth, angiogenesis, scar tissue formation, connective tissue formation and embryonic morphogenesis. Extracellular matrix is increasingly being identified as playing a complex and important role in many biological processes. The collagen proteins are a major component of the extracellular matrix in all mammalian connective tissues and contribute significantly to its structure by forming collagen fibers. Collagen is produced by fibroblasts in the form of procollagen precursors and polymerizes into fibrils, which combine to form a fibrous network or matrix (Alberts et al. , 1994). The procollagen molecules are released via secretory vesicles, which fuse with the cell membrane to create deep, narrow recesses in the fibroblast cell surface. It is in these recesses that the collagen fibrils are formed. Birk & Trelstad, (1986) theorize that these deep recesses give the fibroblast control over the micro-environment in which 2 the collagen fibrils are forming, and thus control over the structure of the collagen matrix. This provides a link between the fibroblast and collagen orientations. Conversely, the collagen matrix is an essential framework which the fibroblasts use as scaffolding to crawl along. Thus the collagen orientation also influences the orientation of fibroblasts and their ability to move. Among the various biological alignment systems, one that has been extensively modeled is the intracellular actin filament network, which shows pronounced alignment patterns in response to the local stress field; stress can be either self-generated or externally applied. This was originally modeled by Sherratt & Lewis, (1993) using a phenomenological spatiotemporal model, enabling in particular the pronounced alignment localized at epithelial wound edges to be studied. In this model, actin alignment is taken to be a function of the ratio of the local, instantaneous, principle components of stress. A more detailed but spatially homogeneous framework formulated as integro-differential equations for the densities of bound and free actin filaments as a function of orientation and time is developed by Civelekoglu & Edelstein-Keshet, (1994). This work has recently been extended by Geigant et al. , (1997), focusing on bifurcations from disorder to alignment. In ecological swarming of macro-organisms and bacteria, orientation plays an important role; here cellular automata is the most prevalent modeling tool (Stevens, 1995; Deutsch, 1995). Cook, (1995) and Grunbaum, (1997) have recently proposed frameworks for reducing integro-differential equations for alignment phenomena to reaction-diffusion-advection equations, which are widely applicable within both ecological and fibroblast culture contexts. In work more directly related to ours, Edelstein-Keshet & Ermentrout, (1990) model the orientation between fibroblasts mediated by cell-to-cell contact. Variables are defined representing densities at one spatial location which depend upon time and an angle of orientation. They use convolutions with different kernels to model the long range angular interaction of the cells, i.e, cells at the same spatial location with different orientations. The variables in those models represent either cells and lamellipodia or bound cells and free cells. One of the major results of their work is to find that patterns can form as a result of only cellular contact responses. This work is generalized and extended by Mogilner & Edelstein- Keshet, (1995) and Mogilner et al. , (1996). In these papers the authors develop two modifications of the model for bound cells and free cells using different assumptions. They find that all three models behave similarly, and they study this behavior by looking at peak like solutions which represent almost complete alignment. In Mogilner & Edelstein-Keshet, (1996) the model is modified to include a spatial component. This work provides a good foundation from which we develop a model specifically suited for 3 fibroblast collagen interactions. Two key differences between our model and this previous work are the focus on the cell-to-cell interactions, which we ignore, and more importantly, the fact that in our system there is no conversion between the variables, i.e., fibroblasts do not become collagen and vice-versa. The rest of the paper is organized in the following manner. In section 2 we present the mathematical model with section 3 giving some analysis of the equations. Section 4 discusses the numerical implementation of the model, and in section 5 some simple numerical tests are given which confirm both the analysis and the numerical scheme. The remaining sections deal with different numerical simulations of the model. First, in section 6 we examine numerical solutions which evolve from initial conditions in which collagen and fibroblasts have identical alignment. Finally, in section 7 we examine the interactions of fibroblasts with collagen. Our findings are discussed in section 8. 2 The Model In our model, f( , ) and c( , ) denote the densities of fibroblasts and collagen fibers respectively at time , oriented at an angle with respect to an arbitrary reference direction. For simplicity we restrict our attention to spatially homogeneous situations. The fibroblasts can be oriented in any direction, so that f is defined for [0, 2 ]. However, collagen fibers are inherently non-directional making fiber orientations of and + equivalent. Thus the model can be formulated either by taking this equivalence into account when formulating the effect of the collagen on fibroblasts, or by defining c only for values of [0, ]; we choose the latter approach. This feature of the system, namely that collagen is nondirectional, has some important consequences for the model predictions, as well as causing technical difficulties. In the formulation which we have chosen, difficulties which recur throughout the paper are associated with the different domains of the variables. Both variables satisfy periodic boundary conditions for themselves and their angular derivatives. The equation for the fibroblasts has a diffusion term, modeling random reorientation, and also a flux term, modeling directed orientation. Biologically, it is known that fibroblasts move up ridges in the substratum (Bray, 1992) and more specifically, experiments by Guido & Tranquillo, (1993) showthat within oriented collagen gels, fibroblasts move preferentially in the direction of collagen orientation by pulling themselves along the fibers. Thus the flux term is due to this tendency of fibroblasts to move in the direction of collagen fibers, so that if there is a gradient of collagen (with orientation, ), the fibroblasts tend to reorient so as to move up that gradient. The flux term can be better understood by 4 comparing it to the standard chemotaxis flux term (Murray, 1993) J = f Copyright Copyright (1) Since we are considering a fixed spatial location, even collagen which is far away in angle space influences the fibroblasts. This is modeled in the attractant term by replacing c with a convolution (defined formally below) which represents a weighted averaging over other orientations. This long range interaction is a fundamental feature of alignment models and we have followed the same approach as Edelstein-Keshet & Ermentrout, (1990) and Mogilner & Edelstein-Keshet, (1995) where they also use convolutions. Thus our model has a flux term involving the density of fibroblasts and the gradient of a convolution term involving the collagen. For simplicity, we take the factor Copyright to be constant. Within the extracellular matrix, collagen takes the form of a fibrous network with an elaborate structure including cross links, and consequently there is essentially no random reorientation of the collagen. Because the fibroblasts degrade and produce collagen, thus reforming the network with collagen oriented in the direction of the fibroblasts, the equation for c has an angular flux term. For simplicity, in our system we do not allow any net change in the amount of collagen; the fibroblasts simply remodel the existing network. As a result there is only one term in the c equation which has three factors: the density of collagen, a convolution with fibroblasts and the gradient of a convolution with fibroblasts. This term is similar to the flux term in the equation for the fibroblasts, but it contains an additional convolution which arises because the rate of collagen remodeling depends on the density of fibroblasts doing that remodeling. The factor is comparable to the Copyright factor in equation 1, but since this is in the evolution equation for collagen, it depends on f, which is the attractant for the c variable. The model consists of the following evolution equations: f = random orie n t at ion D f orientation by collagen 1f (W1 Copyright for [0, 2 ] (2) Copyright = 2 orientation by fibroblasts c(W2 f) (W3 f) for [0, ]. (3) For an alternate motivation of the terms in the above equations the reader is refered to Mogilner & Edelstein-Keshet, (1995) where the authors interpret similar terms using force and angular velocities. 5 The boundary conditions are periodic in , namely f( , 0) = f( , 2 ), f ( , 0) = f ( , 2 ), Copyright , 0) = c( , ),Copyright ( , 0) = c ( , ), (4) and the normalization conditions are given by 2 0 f(0, ) d = 1, 0 c(0, ) d = 1. (5) The normalization conditions simply avoid the necessity to keep track of the initial total densities. The convolution is defined by (W u)( ) = W( x)u(x) dx where the integral is over the domain of u; thus if the convolution involves f the limits of integration are 0 and 2 , while for c the limits are 0 and . In these weighted averages, the kernel W1 is determined by the way the orientation of the fibroblasts is changed due to collagen and W2 and W3 are determined by howthe fibroblasts reorient the collagen. Rescaling time by setting t= 1 simplifies the system, giving the final form of the equations as f t = D f f (W1 Copyright for [0, 2 ] (6) Copyright t = Copyright )(W2 f)( ) (W3 f)( ) for [0, ] (7) where D = D/ 1 and = 2/ 1. The only conditions on W1, W2 and W3 are that W1 is periodic and the others are 2 periodic. Again for convenience, and without loss of generality, we add a normalization requirement so that 2 2 W1( ) d = 1 and W2( ) d = W3( ) d = 1. (8) The detailed behavior of the model of course depends on the forms of W1, W2 and W3, which are determined by biological features of the system that are currently unknown. However some basic properties of the kernels can be deduced from intuitive expectations. In particular, the influence of collagen fibers on fibroblasts and vice-versa should depend only on the magnitude of the angular separation and not on its sign; thus the kernels should be even, Wi(x) = Wi( x) for i = 1 to 3. If the kernels are differentiable at zero then being even implies that W i (0) = 0 for i = 1 to 3. This means that due to the symmetry there is a turning point at zero, which we expect intuitively to be a local maximum. The 2 periodicity of f constrains the kernels W2 and W3 to be 2 periodic but whether the angular distance between a collagen fiber and a fibroblast is or + , the fibroblast has the same influence 6 on the collagen. Taking this into account requires that W2 and W3 be periodic. For simplicity we let W1( ) = 2W2( ) = 2W3( ) and unless stated otherwise we use W1( ) = Ce a 2 for 2 2 the periodic extension otherwise (9) if a is sufficiently small that C exp( a 2/4) > , and otherwise W1( ) = Ce a 2 for 2 1 a ln Copyright 0 for 1 a ln C < 2 < 2 4 the periodic extension otherwise. (10) Thus we use kernels in the form of a Gaussian, but truncated so that the kernels are set to zero whenever the Gaussian would be less than ; in all our simulations we take = 5 10 5. This truncation is significant in a number of our numerical simulations; biologically, it allows the possibility of a limited range of orientations within which fibroblasts and collagen fibers affect one another. The constant C is chosen so that the kernels satisfy (8). Notice that the parameter a determines the support, or range of influence, for these kernels. As was stated earlier, the details of the appropriate shape for these kernels are unknown, and the above form is chosen simply as being intuitively reasonable. In ecological models integro-differential equations are widespread and the kernels are called redistribution kernels, representing dispersal in physical space. In that setting, where more is known experimentally about the shape of redistribution kernels, different kernel shapes have been derived from first principles (Neubert et al. , 1995). A similar approach may be possible in cell biology, and one of the aims of this paper is to suggest potential experimental approaches for this (see section 8). In addition to parameters affecting the kernels, the model contains two dimensionless parameters, namely D, which reflects the angular diffusion coefficient of the cells, and , which reflects the rate at which collagen is re-aligned by the cells. Experimental data is already available in the work of Guido and Tranquillo (1993) which enables the value of D to be estimated as 0.27; we derive this estimate in section 8 of the paper. However, the experimental procedure prevents estimation of the underlying timescales, so that the corresponding dimensional value cannot be estimated. In addition, we are unaware of existing data from which can be estimated, although we suggest appropriate experimental approaches in section 8. 7 3 Analysis of the Model 3.1 General Properties Solutions of the evolution equations (6) and (7) with the associated boundary conditions (4) and initial conditions which satisfy equations (5) have the following properties: 1. If f and c, the fibroblast and collagen densities, are initially non negative, they remain so. This is true provided that the solutions are continuous and differentiable. 2. The total density of collagen and fibroblasts remains constant in time, i.e., 2 0 f(t, ) d = 1 and 0 c(t, ) d = 1 (11) provided f and c are integrable with respect to for each t and ct and ft are continuous. 3. The constant solution f = 1 2 and c = 1 is a steady state. This corresponds to a fully isotropic representation of fibroblasts and collagen fibers. 4. When the fibroblast s diffusion coefficient D = 0, solutions in which all the collagen and fibroblasts are oriented in parallel are also steady states; of course the fibroblast orientations can be in either of the directions parallel to the collagen. More precisely, by considering weak solutions when D = 0 showthat a steady state solution is (f( ),Copyright )) = ( ( a), ( a)) where is the Dirac distribution. This requires that W 1 (0) = 0 and W2(0)W 3 (0) = 0, which are expected intuitively (see section 2) and are satisfied by the kernels in (9, 10). If in addition, we assume that W2 and W3 are periodic, then ( ( + a), ( a)) and ( ( a) + ( + a), ( a)) are steady state solutions also. In some cases, series of delta like functions, corresponding to several isolated orientations, are also steady states; this is discussed in more detail later in the paper (see section 6). 5. Solutions where the delta functions are symmetric with respect to one another are steady states when there is no diffusion. In other words, f = k i=1 ai (x bi) where bi [0, 2 ] can be a steady state solution if for each i, f(bi x) = f(bi + x), with similar conditions for c. These types of solutions are encountered in section 6.3 where they are referred to as type II solutions. 8 3.2 Linear Stability As a first step to understanding the behavior of the system we examine the stability of the constant steady state solution u0 = (fo, co) = ( 1 2 , 1 ) with respect to angularly inhomogeneous perturbations. In order to do this we first linearize equations (6) and (7) about u0. Note that any angularly homogeneous perturbations would be steady state solutions except that those other than uo do not satisfy the normalization condition. By observing that Wi k = k, with all its derivatives zero, for i = 1,2 or 3 and k a constant, one can see that the linearized equations are: f t = D f f0 (W1 Copyright for [0, 2 ] (12) Copyright t = c0(W2 f0)( ) (W3 f)( ) for [0, ]. (13) We assume perturbations of the form f = k=0 e kt (ak cos k + bk sin k ) (14) c = k=0 e kt ( k cos k + k sin k ) (15) where f is 2 periodic and c is periodic, i.e., k = k = 0 for k odd. The periodicity of c implies that ct is also periodic which in turn implies that ak = bk = 0 for k odd. In order to meet the boundary conditions, only even wave numbers need to be considered. First we observe that if W and u are periodic with period T, then when evaluating W u it does not matter over which interval of length T the integral is taken. We chose the interval which is symmetric about zero. Simplification using trigonometric identities then gives W3 f(t, ) = k=0 k even e kt ak cos k W3(k) + bk sin k W3(k) (16) W1 c(t, ) = k=0 k even e kt k cos k W1(k) + k sin k W1(k) (17) where W (k) = a a W(x) cos kx dx. (18) When equations (14) and (15) are substituted into the linearized equations (12) and (13), the coefficients must satisfy the following equations: kak = k2Dak + k2f0 k W1(k) (19) k k = k2 f0c0ak W3(k) (20) 9 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 9 10 Growth rate Wave number 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Growth rate Diffusion a b Figure 1. The growth rate, k, predicted by the linear analysis is plotted for different wave numbers. In graph (a) the growth rate is plotted as a function of the wave number k when D = 1. In graph (b) the maximal growth rates, k = 4 (the dotted line) and k = 2 (the solid line), are plotted as a function of D, the diffusion coefficient. In both cases = 1 and the kernels are defined by a = 4 in equations 9. where ak = (ak, bk), k = ( k, k) and ak = k = 0 for k odd. For k even, ak = 0 and k = 0, k must satisfy the condition k = k2 2 D D2 + 4 c0f2 0 W 1(k) W3(k) . (21) Our assumption that W1 = 2W3 implies 2 W3(k) = W1(k). Thus for k even there is always one positive and one negative k. The mode with the largest positive k is expected to be dominant. Analytical comparison of the growth rates is not possible because of the complex form of W1(k), but by numerically computing k the wave number with the maximum growth rate can easily be found; typical results are shown in figure 1a. Since the odd wave numbers do not satisfy the boundary conditions, the even wave number with the maximum growth rate is 2. This is more clearly seen in figure 1b where 2 and 4 are plotted against D. The maximum dominant mode switches from 4 to 2 as D is increased and once they switch the growth rates remain very similar. Thus the homogeneous steady state is always unstable, and our linear analysis suggests that depending upon the diffusion coefficient, either one or two peaks should begin to evolve corresponding to either one or two orientations for the collagen. This is verified numerically in section 5, and the dependence on the diffusion coefficient is explored more fully later in the paper. By decreasing the range of influence that collagen and fibroblasts have on each other, which means decreasing the support of the kernels, the graph in figure 1a shifts to the right. In the case where D = 0, the mode with the maximum growth rate seems to be the maximum number of independent peaks 10 possible, subject to peak separation being greater than half the support of the kernel (see section 6). 3.3 Nonlinear Analysis Having determined the stability of the uniform steady state using linear stability analysis, we now consider the full nonlinear problem. By so doing we try to understand how the nonlinearities affect the system. Being motivated by previous work where patterns with two characteristic wave lengths are seen (Maini, 1990), we proceed by considering solutions of the form f = a0 + a2(t) cos 2 + a4(t) cos 4 + . . . (22) c = b0 + b2(t) cos 2 + b4(t) cos 4 + . . . . (23) Substituting these into the nonlinear equations (6) and (7), using equations (16) and (17) and integrating the equations after multiplying by 1, cos 2 and cos 4 gives a 0(t) = 0 (24) a 2(t) = 4Da2(t) + 4a0b2(t) W1(2) + 4a2(t)b4(t) W1(4) + . . . (25) a 4(t) = 16Da4(t) + 16a0b4(t) W1(4) + 4a2(t)b2(t) W1(2) + . . . (26) b 0(t) = 0 (27) b 2(t) = 2 b0a2(t) W3(2) a0 4a4(t) W3(4) 4 b2(t)a4(t) W3(4) a0 3a4(t) W3(4) + . . . (28) b 4(t) = 2 b0 a22 (t) W 2 3 (2) + 2a0a4(t) W3(4) b2(t)a2(t) W3(2) 2a0 3a4(t) W3(4) + . . . . (29) The same procedure could be carried out for a more general form of the solution but this suffices for our purposes. The interesting feature of these equations can be seen in equation (26) where the last term depends only on the mode 2 terms and similarly in equation (29) where the second to last term depends only on terms from lower modes. This means that in the full nonlinear model, initial conditions with a purely mode 2 term induce a small mode 4 contribution to the solution. In principle, this in turn induces a mode 8 term in the solution, etc. but because of constraints imposed by the conservation of mass and the lowgro wth rates of the very high modes, these high mode contributions are negligible. We confirm the growth of a mode 4 solution out of mode 2 initial data numerically in section 5. 4 Numerical Implementation Now that we have a basic analytical understanding of the system, we proceed to investigate its behavior in more detail using numerical methods. Before looking at simulations we describe the numerical method 11 used to solve the equations by first dealing with the convolutions and then with the equations. The convolutions are calculated by discretizing the domain with a uniform grid of mesh length h = N and using the left-hand rectangular rule. We denote half the length of the support of the kernel by se. Henceforth, any reference to W1, W2 or W3 means the discretized version, in which the equations are only evaluated at points on the spatial grid. In order to solve the model numerically, the equations are discretized using a centered difference formula while keeping them in conservation form. The time derivative is calculated implicitly using the trapezoidal rule. Thus the diffusive part of the system is simply a Crank-Nicolson discretization. The domain [0, ] is discretized into N points and the domain [0, 2 ] is discretized into 2N points. The resulting discretized nonlinear equations are solved using the software package nksol (Brown & Saad, 1987). This software uses an inexact Newton method to solve the nonlinear system with linear Krylov iterations used to approximate the Newton equations. We refer to solutions of the discretized equations as fn and cm where n [0, 2N 1] and m [0,N 1] denote the grid points (unless otherwise noted N = 100 in our simulations). One property of the true solutions f and c is that if they are initially non negative then they remain so. However, this is not the case for the numerical solutions fn and cm: w hen the densities become sufficiently aggregated and thus there is an abrupt transition from lowdensit y to high density, the numerical solutions become negative and begin to oscillate. This is a common problem in discretized systems with abrupt transitions (Osher & Chakravarthy, 1984). In our case the oscillations are reinforced due to the nature of the equations and quickly grow. Thus, the discretized system does not maintain an important feature of the system - non negative solutions. In the discretized system, the case corresponding to total alignment in one direction is a discrete version of the delta function, in defined to be in = 1 for i n = 0 0 otherwise , (30) where i defines the grid point at which the function is located, and n [0, 2N 1] denotes the grid point. These are the types of solutions we expect to develop in the model, based on our preliminary analysis and intuitive expectation. Yet these are precisely the type of solutions which are difficult to reproduce in a discretized system due to the discontinuities. To get around this problem, a common practice in advective schemes is the use of flux limiters which smooth the oscillations that occur near abrupt changes in the solutions (Thuburn, 1996; Sweby, 1984). We also adopt this approach albeit in a very crude manner. When the solution is about to become negative, the algorithm limits the flux to 12 keep the solution non negative. In practice this is implemented in the following manner: if fn (or cm) is negative, it is added to the larger of its neighbors, fn+1 or fn 1, and then it is set to zero. If its modified neighbor is negative, say fn+1, then fn+1 is also set to zero. The numerical solution nowremain non negative and still conserves mass. This not only acts as a crude flux limiter, but also makes sense in the context of the behavior of the underlying equations. Those equations have the property that when f = 0 then f is also zero, and our ad hoc flux limiter has the effect of reducing the numerical derivative at the point where the solution is set to zero. As always, but particularly in light of the unsophisticated flux limiter added to our algorithm, we interpret and accept our numerical results only as far as they are consistent with the analytical predictions of the model. 5 Numerical Confirmation of the Basic Properties of the Model In order to verify both the numerical scheme and the analysis of section 3 we look at perturbations about the homogeneous steady state ( 1 2 , 1 ). When solving the nonlinear system given by equations (6) and (7) with initial conditions corresponding to a mode 2-type perturbation about the steady state, one sees both the growth of this term, and the appearance of a mode 4 term as predicted in section 3.3. To ensure that the mode 4 term is due to the nonlinearity and is not being introduced through numerical errors, we solve the linear system (12) and (13), with the initial conditions again having purely mode 2 terms. In this case the mode 2 terms growand no other modes appear. As verification of our numerical scheme, we confirmed that the growth rate in the linear and non linear problem are consistent with the analysis. We nowconsider the effect of a single point perturbation to the steady state in both f and c at the same point. The modes with = 4 and = 2 initially grow, and depending on the diffusion coefficient D either two peaks form (small diffusion) or one peak forms (high diffusion). This again agrees with the analysis, although it is hard to determine if the resulting peak formation is due to the different growth rates of the modes or to the constraints imposed by the support of the kernels; one problem is that the two growth rates are very similar for most parameter values. When small random perturbations to the steady state are used as the initial conditions which maintain c(0, ) f(0, ) = f(0, + ), similar results are obtained either one or two peaks form depending on the value of the diffusion coefficient (see figure 2). The results from this section help us to understand the model in two ways. Depending on the diffu- 13 0 /2 Angle 0 5 10 15 20 25 30 Time 0 /2 Angle 0 50 100 150 200 250 Time a b Figure 2. Collagen density is plotted for two simulations with the same initial conditions and different fibroblast diffusion coefficients, illustrating that the collagen evolves to two distinct isolated orientations when the diffusion coefficient is small and to a single orientation when the diffusion coefficient is larger. In (a) D = 0 and the collagen evolves from an almost isotropic state to two distinct orientations; whereas in (b) where D = 0.5, the collagen evolves to a single orientation. This agrees with the results of the linear analysis; in fact in (b) the mode 4 terms (corresponding to 2 peaks on [0, ]) can be seen growing as predicted. In both simulations = 1 and the kernels are defined by a = 4 in (9). The initial conditions are random perturbations of magnitude 0.03 about the homogeneous steady state keeping c(0, ) f(0, ) = f(0, + ). The darker shading between the contour lines indicates higher density. 14 sion coefficient, the long term behavior of the collagen is to have either one or two isolated orientations. The other important result is that even with the ad hoc modification, the numerical scheme captures the basic features of the system predicted by the analysis. This gives us confidence in the numerical simulations, enabling further numerical experiments. 6 Patterns Generated fromSim ple Initial Conditions In order to understand more fully the types of behavior which the system can exhibit, we examine in this section the evolution of numerical solutions starting from simple initial conditions, with f c on the interval [0, ] in all cases. With this condition the fibroblasts and collagen initially have the same orientation, so that one intuitively expects the behavior to be less complicated than when they have different initial orientations; the latter case is considered in section 7. First we shall examine what happens when the collagen and fibroblast orientations are limited to one, two, or three isolated initial directions. We then address the case in which the collagen and fibroblasts have a continuous interval of orientations, and finally we look at solutions which arise from initial conditions with several local maxima. 6.1 Isolated Initial Orientations Of particular interest is the case when all the collagen and fibroblasts are ordered in a few directions. This is modeled by weighted delta functions located at each orientation, with the weights reflecting the fraction of density oriented at each of the angles. Starting with the simplest case in which everything is oriented in one direction, = ih, we have initial conditions of the form fn = a in and cm = d im. In the case D = 0, corresponding to no random reorientation in the fibroblast population, these solutions are steady states, as predicted by the analysis in section 3. As D is increased above zero, the long term solution for the fibroblast density becomes more uniform, eventually approaching a constant f( ) 1 2 ; the collagen density remains localized at a single grid point. Because the domain of f is [0, 2 ], we also consider the initial conditions for f of the type fn = a1 in + a2 i+N n (recall that Nh = ). This corresponds to having fibroblasts oriented both in the direction of the collagen and opposite to that direction. Biologically, one expects this scenario to behave exactly as that above, in which the fibroblasts are all oriented with the collagen, and this is confirmed in simulations. Moving to the next simplest case in which the collagen and fibroblasts have two orientations, = ih 15 and = jh, the initial conditions are fn = a in + a jn and cm = d im + d jm . In this case the two peaks move together until they merge and form a single peak in each variable (see figure 3a). This behavior is determined by the kernels W1, W2 and W3; in the simulations shown in figure 3, these are given by equation (9) with a = 4. This makes the support, se 1.28 radians for W1. Biologically this means that collagen with orientation = influences any fibroblast with orientation [ se, +se]. This influence brings the fibroblasts to an orientation closer to that of the collagen, and the influence of the fibroblasts on the collagen has a similar effect. If se is less than the separation of the two initially imposed peaks, the long range interaction does not extend from one peak to the other, making them independent of each other and one would expect the initial conditions to be a steady state. In figure 3a, the peak separation is about 0.94, which is less than se 1.28, and thus the peaks should influence each other resulting in their merging. Finally we consider the case of three initial localized orientations. If two of these orientations are symmetric about the third then this middle peak remains stationary. Each neighboring peak pulls on the middle peak with the same intensity resulting in no net change in position (see figure 3b). However, the two outer peaks are drawn into the middle, resulting in the eventual merging of all three peaks. There is of course a special case when the three peaks are symmetric about each other which results in an unstable steady state of type II which is more fully explained in section 6.3. This description depends critically on the initial outer peaks having the same intensities, and if the neighboring peaks have differing heights, their influence is no longer symmetric. In order to understand this we mention here two points. The first is that there may be up to three factors which influence a fibroblast peak - the random reorienting of fibroblasts due to a non zero D, the pull of the coinciding collagen peak due to the special initial conditions f Copyright and the pull of neighboring collagen peaks. The second point is that for most parameters, the collagen changes more slowly than the fibroblasts. This is because the rate of collagen remodeling depends on fibroblasts density, reflected mathematically by the flux term having an additional convolution in equation (7) compared to the flux term in equation (6). When neighboring peaks have different heights they pull each other with differing strengths. The smaller peaks of fibroblasts tend to reorient to the direction of the larger neighboring peaks of collagen, yet the collagen peak at the same orientation tends to keep the fibroblasts at their original orientation. If the coinciding collagen peak did not exist, there would not be any reason for the fibroblasts to maintain the original orientation and they would not remain localized, but would rather almost immediately reorient to the direction of the other collagen peaks. The result is that the smaller fibroblast peaks remain localized but they spread out in angle space as they slowly move towards the orientation of the larger 16 0 /2 Angle 0 2 4 6 Time 0 10 20 30 0 0 a 0 /2 Angle 0 5 10 15 Time 0 10 20 30 0 0 b Figure 3. The evolution of collagen densities when collagen and fibroblasts initially have (a) two and (b) three distinct, isolated orientations. In both cases the peaks in orientation merge to give a single isolated collagen orientation, with the fibroblast orientations localized around this. The collagen densities are plotted for simulations where the parameters are the same: D = 0, = 1 and in equation (9) a = 4 making se 1.28 radians. The initial conditions are f = c = 0 except for f30 = c30 = f60 = c60 = 50 in (a) and f30 = c30 = f60 = c60 = f90 = c90 = 100 3 in (b) giving a peak separation of about 0.94 radians. The fibroblast densities look similar. 17 peaks until all the peaks merge at a neworien tation. 6.2 Intervals of Initial Orientation Having confirmed and clarified the behavior of the model for discrete initial orientations, we now consider more complex scenarios. Specifically, we consider initial conditions representing collagen and fibroblasts which have a localized range of orientations. This is described by a hat function which we define on the grid as Hi,k n = 1 for k n i k 0 otherwise . (31) Thus, i denotes the grid point at the center of the hat, 2k+1 denotes the width of the hat, in grid points, and n denotes the grid point where the function is being evaluated. Starting from initial conditions fn = d1Hi,k n and cm = d2Hi,k m with kh< 2 , these initial orientations evolve into densities with all the collagen and fibroblasts oriented in one direction, = ih, described by in (see figure 4a). As with the previous simulations, this result is determined by the type of kernels used in the convolutions, with the support of the kernels being crucial. Recall that the support of the kernel, W1, determines the range of orientations over which collagen influences the fibroblasts, and vice-versa for W2 and W3. When this support is less than the support of the hat function (see figure 4b, c and d), two peaks form initially. This is a continuum version of the case discussed above, in which two peaks are placed symmetrically about a third. At the edge of the hat function the pull is only in one direction - towards the center of the hat function - while far enough into the hat function the pull is equal in both directions, causing no net change. By looking at the convolution of the initial conditions, it is easy to predict where these peaks will form. The convective term draws the densities toward the center of the hat function until the convolution plateaus. Peaks continue to form in the hat until it is filled with non interacting peaks. Peaks are non interacting depending on their separation and the support of the kernels. If the separation is too small, the kernel causes the peaks to interact and merge, otherwise they persist. This link can be easily seen in figure 5 where the support of the hat function, the length of the support of W1 or 2se, and the last time plot of the collagen density are shown. Knowing this, it seems reasonable to expect the wave number with the maximum growth rate to increase as the support of the kernel decreases as was found in section 3.2. Increasing the diffusion coefficient, D, causes the fibroblast density to spread out, rendering a broad range of fibroblast orientations. The main consequence is that peaks may merge (see figure 6b). As the diffusion is increased, the shape of the kernel becomes more important. If the kernel is steep enough 18 0 /2 Angle 0 0.5 1 1.5 2 Time 0 10 20 0 0 /2 Angle 0 0.2 0.4 Time 0 5 0 a b 0 /2 Angle 0 0.05 0.1 0.15 0.2 0.25 Time 0 2 4 6 0.2 0 /2 Angle 0 0.05 0.1 Time 0 1 2 3 T c d Figure 4. Collagen densities in simulations where the effect of varying the angular distance over which fibroblasts and the collagen influence one another is studied. As the distance is decreased, reflected by increasing a in equation (9), more peaks which are independent of one another are able to form. Halving the support almost doubles the number of peaks which are formed. In these four simulations, in order to alter this distance, only a is changed: in (a) a = 4 making se 1.28 radians, in (b) a = 20 making se 0.6 radians, in Copyright a = 75 making se 0.38 radians and in (d) a = 250 making se 0.16 radians. The initial conditions are f = v = H59,22 n with D = 0 and = 1. The fibroblast densities look very similar. These results are plotted in a different way in figure 5 19 0 /2 Angle 0 5 10 15 20 25 30 Density 0 /2 Angle 0 2.5 5 7.5 10 12.5 15 Density a b 0 /2 Angle -2 0 2 4 6 8 10 Density 0 /2 Angle -1 0 1 2 3 4 5 6 7 Density c d Figure 5. Collagen density plotted at one instant in time to showthe relationship between the peak separation and the range of influence between the collagen and the fibroblasts. The dark bar on top shows the length of the support of W1 or 2se and the grey bar at the bottom shows the support of the initial hat function. These simulations are the same as those shown in figure 4. 20 around zero then it keeps the fibroblast peaks very narrow, counteracting the effect of diffusion, but if the kernel is not steep near zero then the diffusion spreads the fibroblast density more readily. As the regions of high fibroblast density become close, they merge into one peak. This occurs to a greater extent as D is increased, and helps to explain why the wave number with the maximum growth rate depends on D, i.e., as D is increased the dominant wave number decreases (see section 3.2). The more spread out the density, the more slowly the peaks form. Finally, at very high diffusion coefficients, peaks in f which are a distance away from those imposed initially are formed (see figure 6c). This is due to the periodic nature of the kernels and of c. If the diffusion causes enough of the density to spread into the interval [ , 2 ], it starts to form peaks corresponding to those in the interval [0, ], provided that the diffusion is not strong enough to level them out. Biologically this corresponds to fibroblasts traveling in either direction with the same orientation as the fibers. 6.3 A Continuumof Orientations with Several Maxima Having established the evolution from a fewdiscrete initial orientations as well as the behavior of an initial localized continuum of orientations, we now look at initial conditions which are a mixture of the two - continuous, with several local maxima. Specifically, we consider initial conditions where fn and cm are proportional to 1 + cos knh and 1 + coskmh respectively; recall that n and m denote mesh points, whereas k and h represent the wave number and the angular step size respectively. Here k must be even to satisfy the periodic boundary conditions for c. Our simulations suggest that these initial conditions evolve into two different types of solutions: delta functions which are independent (see figure 7b) and delta functions which are symmetric with respect to each other (see figure 7a). We refer to these as type I and type II, respectively. When D = 0, the long term behavior is either type I or type II. In the previous simulations all the solutions have been of type I, where the peaks are separated by a distance larger than se making them independent of each other. Type II solutions were mentioned briefly in section 3. The solutions of type II are unstable steady states and as such do not persist biologically. However, studying the evolution of type II solutions from initial data gives valuable insight, which is particularly useful in section 7 where transients similar to type II solutions persist for a long time. Type II solutions are formed with k 2 peaks when k is fairly small: in the case a = 4, we observe such solutions when k < 10. The same factors as before determine the solution type, namely the range of the collagen and fibroblast interactions and the separations of the peaks in their initial densities. If this separation 21 0 Angle 2 0 0.05 0.1 0.15 0.2 Time 0 2 4 6 8 2 0 T 0 Angle 2 0 0.1 0.2 0.3 0.4 Time 0 2 4 2 0 T a b 0 Angle 2 0 1 2 3 Time 0 0.2 0.4 0.6 .8 2 0 T .2 .4 6 Copyright Figure 6. An illustration of the effect of increasing the angular diffusion coefficient on the distribution of fibroblast orientations. As the diffusion coefficient is increased, the orientations become more spread out around the localized collagen orientations. In (a) D = 0 and the fibroblasts end up having four isolated orientations (this is the same simulation as figure 4c where the collagen density is shown). As the diffusion coefficient is increased in (b) to D = 0.1 the peaks no longer represent isolated orientations since they are spread over several grid points and the middle two have merged leaving only three peaks. Increasing D further in Copyright to D = 1.0 shows the formation of three additional peaks, at a distance away from the first set representing fibroblasts orienting in the direction opposite the collagen fibers. The initial conditions and parameters (other than D) are the same as in figure 4c. 22 0 /2 Angle 0 10 20 30 40 50 60 70 Time 0 /2 Angle 0 20 40 60 80 100 Time a b 0 /2 Angle 0 2 4 6 8 10 12 Time 0 /2 Angle 0 5 10 15 20 Time c d Figure 7. An illustration of type I (a,Copyright and type II (b, d) solutions evolving from initial conditions with several orientation maxima. We plot contour lines of the collagen density with darker shading indicating higher densities. As this shows, the solutions change from type II in (a) when there are four local maximum in the initial conditions, or k = 8, to type I in (b) when there are 5 local maximum, or k = 10, and a = 4 in equation (9). By increasing a to a = 10, the range of influence of the collagen and the fibroblasts is decreased and the transition from type II solutions to type I solutions occurs at higher k values, between k = 12 and k = 14. In Copyright k = 10 nowforms a type II solution and in (d) k = 14 forms a type I solution with three independent peaks. The initial conditions are for both f and c are proportional to sin kx + 1. In all cases D = 0 and = 1. The fibroblast densities look similar. 23 is sufficiently large, it allows a quick consolidation to delta functions separated by equal distances. Since the peaks are distributed about each other symmetrically and are of equal heights, they form an unstable steady state (see section 6.1). When the diffusion coefficient is increased from zero, the value of k where the transition from simulations which develop type II solutions to simulations which develop type I solutions will decrease. This is due to the fact that the random reorientation keeps the maximum from coalescing as quickly and compactly, which allows the noise (introduced via numerical errors) to destroy the symmetry. Solutions of type I are stable steady states and are therefore biologically attainable and more significant. Such solutions are formed when k is sufficiently large: for a = 4 we observe these solutions when k > 8, with the solutions always having two peaks. This latter observation is because when a = 4, se 1.28 for W1, which determines that at most two peaks can form on the interval [0, ]. However, because k is large in these cases, the initial densities are sufficiently spread out that it takes a long time for the localization to occur. Mathematically, the explanation for this is that the convolutions of the initial conditions with the kernels are sufficiently smooth that they have very small derivatives, which thus concentrate the density only very slowly. This is so slow that the small asymmetries, arising from the discretization, perturb the system, leading to the formation of type I solutions. In applications, of course, any symmetric initial conditions would be perturbed by natural fluctuations. By decreasing the range of influence of collagen and fibroblasts, more peaks in type I solutions form and type II solutions should form from initial conditions where the density is more evenly distributed. Changing a from 4 to 10 verifies this by causing k = 10 to form five peaked solution of type II (see figure 7c) and k = 14 to form a three peaked solution of type I (see figure 7d). 6.4 Conclusions All of the types of initial conditions considered so far have resulted in solutions where the collagen densities are concentrated at discrete, isolated orientations and the fibroblast densities are localized around these discrete orientations, with the degree of aggregation increasing as the angular diffusion coefficient goes down. There are two types of solutions which have been observed: type I solutions where the peaks are independent of one another, and type II solutions where the solution is symmetric about each peak. Only type I solutions are physically relevant since type II solutions are unstable, to asymmetric perturbations, developing into type I solutions. The independence of the peaks in type I solutions is due to the half-range of influence of the collagen and fibroblasts being less than the 24 separation of the peaks. Thus the support of the kernel, which corresponds biologically to the range of directions over which the collagen and fibroblasts are able to reorient one another, determines the maximum number of peaks that can occur in type I solutions. 7 Interaction between Fibroblasts and Collagen Previously the initial conditions have been chosen such that the collagen and fibroblast densities were proportional to each other, in order to reduce the extent of their interaction, and thus simplify the results. Now we change this strategy and investigate the interaction between the fibroblasts and the collagen. When their initial configurations are different, the parameter becomes very important. Recall that determines howstrongly the collagen is reordered. Figure 8 illustrates the results from a simulation in which the collagen is initially all set at one angle and the fibroblasts at another. When is small the collagen s initial conditions are more important in determining the final solution, and the fibroblasts reorient in the direction of the collagen (figure 8c). When is increased sufficiently, the situation is reversed and the fibroblasts initial conditions become more important, causing the collagen to reorient in the direction of the fibroblasts (figure 8d). At intermediate values of both the fibroblast and the collagen orientations alter, stabilizing at some intermediate orientation (figure 8a). Increasing the diffusion coefficient simply makes the fibroblast density more spread out, as seen in figure 8b. If the initial conditions are changed to hat functions, similar results are obtained, but with some additional complications. The hat functions can become split depending on their width and the separation. Figure 8e illustrates the case when the initial densities for fibroblasts and collagen are both hat functions, but with different weights. As expected, the hat function with more density draws the collagen to a greater extent than the function with less density. They all merge to form single coinciding peaks. The final type of simulation that we have used to understand the interaction of fibroblasts and collagen have initial conditions with fn proportional to 1+sin jnh and cm proportional to 1+sinkmh. Here we are using the same notation as above, but the wave numbers j for the variable f and k for the variable c can be different. In these simulations the collagen has several orientation maxima which are different from the orientation maxima of the fibroblasts. The short term behavior of the solution depends on : if is sufficiently small, then the fibroblasts initially try to reorient to the form of the initial collagen density (see figure 9a), and vice-versa if is large (see figure 9b). However, these initial reorientations do not persist long enough to significantly alter the long term behavior of the collagen and 25 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time a ( = 3, D = 0) b ( = 3, D = 0.1) 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time Copyright( = 0.1, D = 0) d ( = 30, D = 0) e ( = 20, D = 0) Figure 8. Fibroblast and collagen densities in simulations which demonstrate how they influence each other. Initially, the fibroblasts and collagen have separate orientations. The evolution depends on the parameters , which determines how strongly the fibroblasts reorient the collagen, and D, the angular diffusion coefficient of the fibroblasts. The simulations shown in (a), Copyright and (d) differ only in the value of which is used. In (a) where = 3, the fibroblasts and collagen roughly have the same influence on each other causing the final orientation to be at an intermediate value from the initial ones. In Copyright where = 0.1, the fibroblasts are drawn to the orientation of the collagen, whereas in (d) where = 30, the collagen is drawn to the orientation of the fibroblasts. By comparing (b) with (a) the effect of the diffusion coefficient is shown. The only difference between the simulation shown in (a) and (b) is that the diffusion coefficient is changed from D = 0 to D = 0.1 respectively. In (b) the fibroblasts have a greater range of orientations. In (e) the collagen is drawn closer to the orientation of the greater mass of fibroblasts. Here = 20 and the initial conditions are hat functions with different heights, that is f 2H5,2 n + H40,2 n and c H20.5,5 n . (The non integer value in the superscript for H indicates that the hat function is centered between grid points.) The region for the collagen density when c > 1 is shaded, and the contour for f = 1 is a bold shaded line. Unless otherwise stated, D = 0 and the initial conditions are c 35 n and f 15 n . 26 only moderately influence the fibroblasts. The variables reorient in a manner largely determined by the initial conditions; this can be seen in figure 9. In the context of scar tissue formation, this suggests that the initial deposition of collagen is the most important factor in determining its long term orientation. More specifically, the model shows that the remodeling of the collagen by the fibroblasts alters the collagen orientation in a manner which is less dependent on the properties and initial conditions of the fibroblasts than on the initial conditions of the collagen. Yet in scar tissue formation, where the fibroblasts produce and degrade collagen (processes which are purposely not included in this model) the initial conditions of the fibroblasts may be crucial to the initial deposition of collagen. Despite this, once the collagen alignment is established and possibly just initiated, the conclusion of our model holds and remains valid for the period of collagen remodeling which continues for months post wounding (Mast, 1992). These results are consistent with the observation that anti-scarring therapies such as application of TGF s are only effective in the very early phase of wound repair (Shah et al. , 1994), despite the fact that scar remodeling occurs on a relatively long time scale. Restating, this model indicates that the parameter is important for the initial behavior of the solutions, but the initial density distributions are more significant for the final form of the solution. It is interesting to note that in these simulations, the fibroblasts and collagen can form isolated orientations which persist for long transient periods. This is contrary to behavior described in section 6.3, where the initial maxima did not become compact before coalescing. The fact that the peaks here do become compact explains why they persist for long periods. Their influence on each other is in the tail of the kernel, making it small until they move closer. This type of solution could be biologically significant, for example if the fibroblasts become inactive during this long transient. A final situation to consider is when both fibroblast and collagen densities have random initial conditions, set by randomly choosing a value between zero and one for each grid point. Those values are then rescaled so that fn and cm satisfy the normalization condition. As expected from previous considerations, we have found in a large number of simulations that the long time behavior involves the maximum number of peaks which can form in a type I solution (not illustrated for brevity). Changes in the parameter simply change the time it takes for the peaks to form; if is small then f forms peaks very quickly, while if is large then c forms peaks very quickly. If the diffusion coefficient is increased, the f solution smooths out and the c solution takes longer to form peaks. 27 0 /2 Angle 0 10 20 30 40 50 60 Time 81% of the fibroblast density divided between this peak and the corresponding peak in the interval to 2 . Similarly 19% of the fibroblast density is divided between this peak and that on to 2 . The two collagen peaks shown below eventually merge into one. 0 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 96% of the collagen density at this peak 4% of the collagen density at this peak a b Figure 9. Collagen and fibroblast densities from simulations which illustrate that initial conditions have a greater influence than the parameter on the final solutions. In (a) although = 0.1 causing the fibroblasts to try and quickly reorient in the direction of the collagen, the final solution has most of the fibroblast density at two peaks corresponding to the fibroblasts initial condition with two maxima, or k = 2 (only half the fibroblast domain is shown). The long term behavior of the collagen is two peaks with roughly half the collagen density in each, again a consequence of the initial conditions for collagen with four maxima, or k = 8. In (b) = 150 causing the collagen to try and take the orientation of the fibroblasts, but again the initial conditions prevail and most of the collagen density is oriented at one peak corresponding to the initial conditions with k = 2. The long term behavior of the fibroblast density is four peaks, two have about 75 percent of the total density divided equally between them and the rest is divided equally between the remaining two peaks. The initial conditions for the fibroblast has four maxima, or k = 8. In (a) the shaded region shows where the collagen density is greater than 0.4 and the dark line is the contour for fibroblast density 0.1 and in (b) the shaded region shows where the collagen density is greater than 0.5 and the dark line is the contour for fibroblast density 0.27. Both simulations have D = 0 and initial conditions where f and c are proportional to sin k + 1. 28 8 Discussion In this paper we have developed a simple model for fibroblast and collagen alignment interactions and investigated its behavior using a combination of analytical and numerical techniques. The results have yielded a number of insights into the alignment process, and we begin the discussion by considering the application of these insights to wound healing and cancer. Dermal wound repair is currently a very active research topic, due to recent advancements which promise to lead to newclinical techniques for reducing scarring (McCallion & Ferguson, 1996). Despite a large volume of experimental research (see Clark, 1996a for review) and some mathematical modeling (Olsen et al. , 1995; Dale et al. , 1996), many details of the process remain poorly understood. It is known that collagen alignment plays a key role in the healing process; in fact, collagen alignment is one method for characterizing scar quality. In humans and other tight skinned animals, collagen has a cross weave structure in normal tissue, whereas in scar tissue it is aligned parallel to the plane of the skin (Harmon et al. , 1995; Welch et al. , 1990). As a dermal wound is repaired, fibroblasts replace the provisional matrix of fibrin with a collagen matrix (Clark, 1996b). The collagen is then reorganized for months by the fibroblasts until at some time they become quiescent and the matrix remains relatively unchanged (Mast, 1992), corresponding mathematically to a steady state. In our model, the configuration observed in normal skin can be represented by the collagen density having two peaks of orientation, with roughly 90o separation, and scar tissue can be represented by a collagen density profile with one alignment peak. Both of these solutions are indeed steady state solutions of our model. Moreover the model predicts that both of these steady states can be stable for the same parameters. Thus the properties of the fibroblasts, characterized by the parameters D, and the kernels W1, W2 and W3, need not be changed in order to obtain either type of solution. Rather, our model predicts that it is the initial conditions which determine which of the steady states form. This is consistent with biological observations that transient application of growth factors can permanently alter the quality of repair (Shah et al. , 1994). Cancer invasion is another area of application in which the key process is the movement of a cell population through a collagen-dominated extracellular matrix. For most cancers, the relevant cell type is transformed epithelial cells, but these share many phenotypic similarities with fibroblasts, and our model is quite applicable in this context. In tumor invasion, the role of cell and matrix orientation has received relatively little attention, with recent work focusing instead on the details and interaction of protease production, directed cell movement and altered cell adhesion; a reviewof recent experimental 29 work in this area is given in Jiang & Mansel, (1996), and modeling approaches are described in Byrne & Chaplain, (1996) and Perumpanai et al. , (1996). Determination of the role of collagen reorientation during invasion is an important modeling challenge for which our work lays the foundations. A number of investigators have studied the interaction between fibroblasts and collagen using in vitro experiments involving fibroblasts in collagen gels (Clark et al. , 1995; Guido & Tranquillo, 1993; Stopak & Harris, 1982). Our model relates directly to this type of experiment, and our results suggest a series of experiments that could be performed in order to estimate the model parameters: 1. The first and most fundamental experiment is to introduce fibroblasts into a collagen gel in which all the collagen is oriented in a single direction. This procedure has in fact been performed by Guido & Tranquillo, (1993), and their data can be used to estimate the dimensionless diffusion coefficient D, as mentioned in section 2. Specifically, Guido & Tranquillo present a histogram of fibroblast orientations, which can be conveniently summarized by the standard deviation away from the mean, which is approximately 0.39 radians (the mean is of course the predominant collagen direction). Simulations of our model imply that in a uni-directional collagen network, the standard deviation of the fibroblast orientations is an increasing function of D, as expected intuitively, and the standard deviation of 0.39 radians implies that D 0.27. Strictly, this is only an upper bound on D, because in the experiments, the collagen deviates to some extent from being uni-directional, to an extent that cannot be determined quantitatively. However, model simulations showthat fibroblast distribution is in fact relatively insensitive to this deviation in collagen distribution, suggesting that D 0.27 is a reasonable estimate. Unfortunately, the experimental procedure in Guido & Tranquillo, (1993) means that only the steady state is considered, so that no dimensional information is available. However, if it were possible to modify their procedure and introduce the cells after aligning the gel, this type of information would be accessible. 2. Further experiments would require extensions of the procedure of Guido & Tranquillo, (1993), so that the gel has two (or more) isolated collagen orientations. The natural approach, in keeping with the development of this paper, is to begin by constructing experimental initial conditions in which there are two different isolated collagen orientations, with fibroblast directions localized around these as in the simulations for figure 3a. Measurement of the time evolution of the distribution of fibroblast directions could be compared with model simulations, to delineate D and W1. We suggest that this may be achievable by obtaining two separate uni-directional collagen gels with corresponding cell populations, as in (1) above, and then placing one above the other, at an angle. 30 3. As in the theoretical development in the paper, the next experimental step would be to construct experimental initial conditions in which fibroblasts and collagen had different orientations. One possible approach to achieving this would be to juxtapose two collagen gels, with the collagen unidirectional in one, and oriented in two different isolated directions in the other. The fibroblasts could be placed in the uni-directional gel, from where they would enter the other gel with one predominant direction of motion. Measurement of the long-term distribution of collagen and fibroblast orientations near the interface, compared to model simulations, would then enable , W1, W2 and W3 to be determined. For instance, if is very large the collagen in the dualdirectional gel should become uni-directional, whereas if is small the fibroblasts should become oriented in both directions of the gel. Although the spatial aspect, which is ignored in our model, will certainly play a role, in this region the local dynamics should dominate. The work in this paper provides a theoretical framework on which more complex models could be built. There are three important modifications that could naturally be made to the model in order to make it more widely applicable. The most important is to add a spatial component, and study the way in which orientations develop as the fibroblasts move spatially. We are particularly interested in the transitions between two regions with different characteristic patterns of orientation, such as from scar tissue to normal tissue. Such an extension would greatly increase the complexity of the model, requiring the collagen and fibroblast densities to be a function of time, orientation, and two spatial coordinates. The second important modification would be to add a term to the model representing production of collagen by fibroblasts. This is relevant in a number of applications: for example, scar tissue has a greater density of collagen than does normal tissue, due to production by fibroblasts entering the wound (Shah et al. , 1992). Finally, we would like to extend our model to three space dimensions by using two angular variables. Conceptually this is straightforward, but the resulting model would be computationally much more intensive, and experimental verification of the two-dimensional framework is an essential precursor to this extension. Acknowledgements. We thank Philip Maini (Oxford), Mark Ferguson (Manchester) and Markus Owen (Warwick) for helpful discussions. This work was supported by grant GR/K71394 from the EPSRC, and by a grant from the London Mathematical Society (scheme 3). 31 References Alberts, B., Dennis, B., Lewis, J., Raff, M., Roberts, K., & Watson, J. D. 1994.Molecular Biology of the Cell. 3 edn.Garland Publishing.Chap. Cell Junctions, Cell Adhesion, and the Extracellular Matrix, pages 978 986. Besseau, L., & Giraud-Guille, M. M. 1995. Stabilization of Fluid Cholesteric Phases of Collagen to Ordered Gelated Matrices. J. Mol. Biol., 251(2), 197 202. Birk, D. E., & Trelstad, R. L. 1986. Extracellular Comparments in Tendon Morphogenesis: Collagen Fibril, Bundle, and Macroaggregate Formation.J. Cell Biol., 103, 231 240. Bray, D. 1992.Cell Movements.Garland. Brown, P. N., & Saad, Y. 1987.Hybrid Krylov Methods for Nonlinear Systems of Equations.Tech. rept. Lawrence Livermore National Laboratory. Byrne, H. M., & Chaplain, M. A. J. 1996. Modelling the role of cell-cell adhesion in the growth and development of a carcinoma.Math. Comp. Modelling, 24, 1 17. Civelekoglu, G., & Edelstein-Keshet, L. 1994.Modelling the dynamics of F-actin in the cell.Bull. Math. Biol., 56, 587 616. Clark, R. A. F. (ed). 1996a.The Molecular and Cellular Biology of Wound Repair. 2 edn.Plenum Press. Clark, R. A. F. 1996b.Wound Repair Overviewand General Considerations. Pages 3 50 of: Clark, R. A. F. (ed), The Molecular and Cellular Biology of Wound Repair, 2 edn.Plenum Press. 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. Cook, J. 1995.Waves of alignment in populations of interacting, oriented individuals. Forma, 10, 171 203. 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. Deutsch, A. 1995.Towards analysing complex swarming patterns in biological systems with the help of lattic-gas cellular automata.J. Biol. Systems, 3, 947 956. Edelstein-Keshet, L., & Ermentrout, B. G. 1990.Models for contact-mediated pattern formation: cells that form parallel arrays. J. Math. Biol., 29, 33 58. 32 Elsdale, T. 1973. The generation and maintenance of parallel arrays in cultures of diploid fibroblasts. In: Kulonen, E., & Pikkarainen, J. (eds), Biology of Fibroblast.NewY ork: Academic Press. Geigant, E., Ladizhansky, K., & Mogilner, A. 1997. An integro-differential model for orientational distributions of F-actin in cells.SIAM J. Appl. Math.To appear. Grunbaum, D. 1994. Swarming behaviour as an Aide to Chemotaxis. In: Parrish, J., & Hammner, W. (eds), 3D Animal Aggregations.Cambridge University Press. Grunbaum, D. 1997. Advection-diffusion equations for generalised tactic searching behaviors. J. Math. Biol. in press. 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. Harmon,Copyright B., Zelickson, B. D., Roenigk, R. K., Wayner, E. A., Hoffstrom, B., Pittelkow, M. R., & Brdoland, D. G. 1995.Dermabrasive scar revision - immunochemical and ultrastructural evaluation. Dermatol. Surg., 21, 503 508. Jiang, W. G., & Mansel, R. E. 1996.Progress in anti-invasion and anti-metastasis research and treatment. Int. J. Oncology, 9, 1013 1028. Maini, P. K. 1990.Superposition of modes in a caricature of a model for morphogenesis.J. Math. Biol., 28, 307 315. Mast, B. A. 1992.The Skin.Chap. 22, pages 344 355 of: Cohen, I. K., Diegelmann, R. F., & Lindblad, W. J. (eds), Wound Healing Biochemical and Clinical Aspects.Saunders. McCallion, R. L., & Ferguson, M.W. J. 1996.FetalWound Healing and the Development of Antiscarring Therapies for Adult Wound Healing. Pages 561 600 of: Clark, R. A. F. (ed), The Molecular and Cellular Biology of Wound Repair, 2 edn.Plenum Press. Mogilner, A., & Edelstein-Keshet, L. 1995.Selecting a common direction I. Howorien tational 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. Murray, J. D. 1993.Mathematical Biology. 2 edn.Springer.Chap. 9, pages 241 244. 33 Neubert, M. G., Kot, M., & Lewis, M. A. 1995. Dipersal and Pattern Formation in a Discrete-Time Predator-Prey Model.Theoretical Population Biology, 48(1), 7 43. Okuba, A. 1986.Dynamical aspects of animal grouping: swarms, schools flocks, and herds.Adv. Biophys., 22, 1 94. Olsen, L., Sherratt, J. A., & Maini, P. K. 1995. A Mechanochemical Model fo Adult Dermal Wound Contraction and the Permanence of the Contrated Tissue Displacement Profile. J. Theor. Biol., 177, 113 128. Osher, S., & Chakravarthy, S. R. 1984. High resolution schemes and the entropy condition. SIAM J. Numer. Anal., 21(5), 955 984. Perumpanai, A. J., Sherratt, J. A., Norbury, J., & Byrne, H. M. 1996. Biological inferences from a mathematical model for malignant invasion.Invasion and Metastasis, 16, 209 221. Pollard, T. D., & Cooper, J. A. 1986.Actin and actin-binding proteins. A critical evaluation of mechanisms and function.Ann. Rev. Biochem., 55, 987 1035. Shah, M., Foreman, D. M., & Ferguson, M. W. J. 1992.Control of scarring in adult wounds by neutralising antibody to transforming growth factor .Lancet, 339, 213 214. Shah, M., Foreman, D. M., & Ferguson, M. W. J. 1994. Neutralising antibody to TGF- 1,2 reduces cutaneous scarring in adult rodents. J. Cell Sci., 107, 1137 1157. Sherratt, J. A., & Lewis, J. 1993. Stress-induced alignment of actin filaments and the mechanics of cytogel.Bull. Math. Biol., 55, 637 654. Stevens, A. 1995.Trail following and aggregation of myxobacteria.J. Biol. Systems, 3, 1059 1068. Stopak, D., & Harris, A. K. 1982.Connective Tissue Morphogenesis by Fibroblast Traction.Devel. Biol., 90, 383 398. Sweby, P. 1984.High Resolution Schemes using Flux Limiters for Hyperbolic Conservation Laws.SIAM J. Numer. Anal., 21(5), 995 1011. Thuburn, J. 1996.Multidimensional Flux-Limited Advection Schemes.J. Comp. Phys., 123, 74 83. Welch, M. P., Odland, G. F., & Clark, R. A. F. 1990. Temporal Relationships of F-actin bundle formation, collagen and fibronectin matrix assembly, and fibronectin receptor expression to wound contraction. J. Cell Biol., 110, 133 145. 34

Click tabs to swap between content that is broken into logical sections.

Due to the increasing importance of the extracellular matrix in many biological problems, in this paper we develop a model for fibroblast and collagen orientation with the ultimate objective of understanding how fibroblasts form and remodel the extracellular matrix, in particular its collagen component. The model uses integro-differential equations to describe the interaction between the cells and fibers at a point in space with various orientations. The equations are studied both analytically and numerically to discover different types of solutions and their behavior. In particular we examine solutions where all the fibroblasts and collagen have discrete orientations, a localized continuum of orientations and a continuous distribution of orientations with several maxima. The effect of altering the parameters in the system is explored, including the angular diffusion coefficient for the fibroblasts, as well as the strength and range of the interaction between fibroblasts and collagen. We find the initial conditions and the range of influence between the collagen and the fibroblasts are the two factors which determine the behavior of the solutions. The implications of this for wound healing and cancer are discussed, including the conclusion that the major factor in determining the degree of scarring is the initial deposition of collagen.

Published in Bulletin of Mathematical Biology (1998)60, 101 129 A Mathematical Model for A Mathematical Model for Fibroblast and Collagen Orientation
Fibroblast and Collagen Orientation John C.
Dallon and Jonathan A. Sherratt Mathematics Institute University of Warwick Coventry CV4 7AL, U.K. February 15, 2000 RUNNING TITLE: An orientation model 1 Abstract Due to the increasing importance of the extracellular matrix in many biological problems, in this paper we develop a model for fibroblast and collagen orientation with the ultimate objective of understanding howfibroblasts form and remodel the extracellular matrix, in particular its collagen component. The model uses integro-differential equations to describe the interaction between the cells and fibers at a point in space with various orientations. The equations are studied both analytically and numerically to discover different types of solutions and their behavior. In particular we examine solutions where all the fibroblasts and collagen have discrete orientations, a localized continuum of orientations and a continuous distribution of orientations with several maxima. The effect of altering the parameters in the system is explored, including the angular diffusion coefficient for the fibroblasts, as well as the strength and range of the interaction between fibroblasts and collagen. We find the initial conditions and the range of influence between the collagen and the fibroblasts are the two factors which determine the behavior of the solutions. The implications of this for wound healing and cancer are discussed, including the conclusion that the major factor in determining the degree of scarring is the initial deposition of collagen. 1 Introduction Alignment within biological systems has been the subject of considerable recent interest. Applications vary over a wide variety of systems including herd movement, flocks of birds, schools of fish, insect swarms (Grunbaum, 1994; Okuba, 1986), cellular movement, actin networks and collagen networks (Elsdale, 1973; Pollard & Cooper, 1986; Besseau & Giraud-Guille, 1995). This paper focuses on cellular alignment with respect to collagen fibers. There are numerous applications of this specific system including tumor growth, angiogenesis, scar tissue formation, connective tissue formation and embryonic morphogenesis. Extracellular matrix is increasingly being identified as playing a complex and important role in many biological processes. The collagen proteins are a major component of the extracellular matrix in all mammalian connective tissues and contribute significantly to its structure by forming collagen fibers. Collagen is produced by fibroblasts in the form of procollagen precursors and polymerizes into fibrils, which combine to form a fibrous network or matrix (Alberts et al. , 1994). The procollagen molecules are released via secretory vesicles, which fuse with the cell membrane to create deep, narrow recesses in the fibroblast cell surface. It is in these recesses that the collagen fibrils are formed. Birk & Trelstad, (1986) theorize that these deep recesses give the fibroblast control over the micro-environment in which 2 the collagen fibrils are forming, and thus control over the structure of the collagen matrix. This provides a link between the fibroblast and collagen orientations. Conversely, the collagen matrix is an essential framework which the fibroblasts use as scaffolding to crawl along. Thus the collagen orientation also influences the orientation of fibroblasts and their ability to move. Among the various biological alignment systems, one that has been extensively modeled is the intracellular actin filament network, which shows pronounced alignment patterns in response to the local stress field; stress can be either self-generated or externally applied. This was originally modeled by Sherratt & Lewis, (1993) using a phenomenological spatiotemporal model, enabling in particular the pronounced alignment localized at epithelial wound edges to be studied. In this model, actin alignment is taken to be a function of the ratio of the local, instantaneous, principle components of stress. A more detailed but spatially homogeneous framework formulated as integro-differential equations for the densities of bound and free actin filaments as a function of orientation and time is developed by Civelekoglu & Edelstein-Keshet, (1994). This work has recently been extended by Geigant et al. , (1997), focusing on bifurcations from disorder to alignment. In ecological swarming of macro-organisms and bacteria, orientation plays an important role; here cellular automata is the most prevalent modeling tool (Stevens, 1995; Deutsch, 1995). Cook, (1995) and Grunbaum, (1997) have recently proposed frameworks for reducing integro-differential equations for alignment phenomena to reaction-diffusion-advection equations, which are widely applicable within both ecological and fibroblast culture contexts. In work more directly related to ours, Edelstein-Keshet & Ermentrout, (1990) model the orientation between fibroblasts mediated by cell-to-cell contact. Variables are defined representing densities at one spatial location which depend upon time and an angle of orientation. They use convolutions with different kernels to model the long range angular interaction of the cells, i.e, cells at the same spatial location with different orientations. The variables in those models represent either cells and lamellipodia or bound cells and free cells. One of the major results of their work is to find that patterns can form as a result of only cellular contact responses. This work is generalized and extended by Mogilner & Edelstein- Keshet, (1995) and Mogilner et al. , (1996). In these papers the authors develop two modifications of the model for bound cells and free cells using different assumptions. They find that all three models behave similarly, and they study this behavior by looking at peak like solutions which represent almost complete alignment. In Mogilner & Edelstein-Keshet, (1996) the model is modified to include a spatial component. This work provides a good foundation from which we develop a model specifically suited for 3 fibroblast collagen interactions. Two key differences between our model and this previous work are the focus on the cell-to-cell interactions, which we ignore, and more importantly, the fact that in our system there is no conversion between the variables, i.e., fibroblasts do not become collagen and vice-versa. The rest of the paper is organized in the following manner. In section 2 we present the mathematical model with section 3 giving some analysis of the equations. Section 4 discusses the numerical implementation of the model, and in section 5 some simple numerical tests are given which confirm both the analysis and the numerical scheme. The remaining sections deal with different numerical simulations of the model. First, in section 6 we examine numerical solutions which evolve from initial conditions in which collagen and fibroblasts have identical alignment. Finally, in section 7 we examine the interactions of fibroblasts with collagen. Our findings are discussed in section 8. 2 The Model In our model, f( , ) and c( , ) denote the densities of fibroblasts and collagen fibers respectively at time , oriented at an angle with respect to an arbitrary reference direction. For simplicity we restrict our attention to spatially homogeneous situations. The fibroblasts can be oriented in any direction, so that f is defined for [0, 2 ]. However, collagen fibers are inherently non-directional making fiber orientations of and + equivalent. Thus the model can be formulated either by taking this equivalence into account when formulating the effect of the collagen on fibroblasts, or by defining c only for values of [0, ]; we choose the latter approach. This feature of the system, namely that collagen is nondirectional, has some important consequences for the model predictions, as well as causing technical difficulties. In the formulation which we have chosen, difficulties which recur throughout the paper are associated with the different domains of the variables. Both variables satisfy periodic boundary conditions for themselves and their angular derivatives. The equation for the fibroblasts has a diffusion term, modeling random reorientation, and also a flux term, modeling directed orientation. Biologically, it is known that fibroblasts move up ridges in the substratum (Bray, 1992) and more specifically, experiments by Guido & Tranquillo, (1993) showthat within oriented collagen gels, fibroblasts move preferentially in the direction of collagen orientation by pulling themselves along the fibers. Thus the flux term is due to this tendency of fibroblasts to move in the direction of collagen fibers, so that if there is a gradient of collagen (with orientation, ), the fibroblasts tend to reorient so as to move up that gradient. The flux term can be better understood by 4 comparing it to the standard chemotaxis flux term (Murray, 1993) J = f Copyright Copyright (1) Since we are considering a fixed spatial location, even collagen which is far away in angle space influences the fibroblasts. This is modeled in the attractant term by replacing c with a convolution (defined formally below) which represents a weighted averaging over other orientations. This long range interaction is a fundamental feature of alignment models and we have followed the same approach as Edelstein-Keshet & Ermentrout, (1990) and Mogilner & Edelstein-Keshet, (1995) where they also use convolutions. Thus our model has a flux term involving the density of fibroblasts and the gradient of a convolution term involving the collagen. For simplicity, we take the factor Copyright to be constant. Within the extracellular matrix, collagen takes the form of a fibrous network with an elaborate structure including cross links, and consequently there is essentially no random reorientation of the collagen. Because the fibroblasts degrade and produce collagen, thus reforming the network with collagen oriented in the direction of the fibroblasts, the equation for c has an angular flux term. For simplicity, in our system we do not allow any net change in the amount of collagen; the fibroblasts simply remodel the existing network. As a result there is only one term in the c equation which has three factors: the density of collagen, a convolution with fibroblasts and the gradient of a convolution with fibroblasts. This term is similar to the flux term in the equation for the fibroblasts, but it contains an additional convolution which arises because the rate of collagen remodeling depends on the density of fibroblasts doing that remodeling. The factor is comparable to the Copyright factor in equation 1, but since this is in the evolution equation for collagen, it depends on f, which is the attractant for the c variable. The model consists of the following evolution equations: f = random orie n t at ion D f orientation by collagen 1f (W1 Copyright for [0, 2 ] (2) Copyright = 2 orientation by fibroblasts c(W2 f) (W3 f) for [0, ]. (3) For an alternate motivation of the terms in the above equations the reader is refered to Mogilner & Edelstein-Keshet, (1995) where the authors interpret similar terms using force and angular velocities. 5 The boundary conditions are periodic in , namely f( , 0) = f( , 2 ), f ( , 0) = f ( , 2 ), Copyright , 0) = c( , ),Copyright ( , 0) = c ( , ), (4) and the normalization conditions are given by 2 0 f(0, ) d = 1, 0 c(0, ) d = 1. (5) The normalization conditions simply avoid the necessity to keep track of the initial total densities. The convolution is defined by (W u)( ) = W( x)u(x) dx where the integral is over the domain of u; thus if the convolution involves f the limits of integration are 0 and 2 , while for c the limits are 0 and . In these weighted averages, the kernel W1 is determined by the way the orientation of the fibroblasts is changed due to collagen and W2 and W3 are determined by howthe fibroblasts reorient the collagen. Rescaling time by setting t= 1 simplifies the system, giving the final form of the equations as f t = D f f (W1 Copyright for [0, 2 ] (6) Copyright t = Copyright )(W2 f)( ) (W3 f)( ) for [0, ] (7) where D = D/ 1 and = 2/ 1. The only conditions on W1, W2 and W3 are that W1 is periodic and the others are 2 periodic. Again for convenience, and without loss of generality, we add a normalization requirement so that 2 2 W1( ) d = 1 and W2( ) d = W3( ) d = 1. (8) The detailed behavior of the model of course depends on the forms of W1, W2 and W3, which are determined by biological features of the system that are currently unknown. However some basic properties of the kernels can be deduced from intuitive expectations. In particular, the influence of collagen fibers on fibroblasts and vice-versa should depend only on the magnitude of the angular separation and not on its sign; thus the kernels should be even, Wi(x) = Wi( x) for i = 1 to 3. If the kernels are differentiable at zero then being even implies that W i (0) = 0 for i = 1 to 3. This means that due to the symmetry there is a turning point at zero, which we expect intuitively to be a local maximum. The 2 periodicity of f constrains the kernels W2 and W3 to be 2 periodic but whether the angular distance between a collagen fiber and a fibroblast is or + , the fibroblast has the same influence 6 on the collagen. Taking this into account requires that W2 and W3 be periodic. For simplicity we let W1( ) = 2W2( ) = 2W3( ) and unless stated otherwise we use W1( ) = Ce a 2 for 2 2 the periodic extension otherwise (9) if a is sufficiently small that C exp( a 2/4) > , and otherwise W1( ) = Ce a 2 for 2 1 a ln Copyright 0 for 1 a ln C < 2 < 2 4 the periodic extension otherwise. (10) Thus we use kernels in the form of a Gaussian, but truncated so that the kernels are set to zero whenever the Gaussian would be less than ; in all our simulations we take = 5 10 5. This truncation is significant in a number of our numerical simulations; biologically, it allows the possibility of a limited range of orientations within which fibroblasts and collagen fibers affect one another. The constant C is chosen so that the kernels satisfy (8). Notice that the parameter a determines the support, or range of influence, for these kernels. As was stated earlier, the details of the appropriate shape for these kernels are unknown, and the above form is chosen simply as being intuitively reasonable. In ecological models integro-differential equations are widespread and the kernels are called redistribution kernels, representing dispersal in physical space. In that setting, where more is known experimentally about the shape of redistribution kernels, different kernel shapes have been derived from first principles (Neubert et al. , 1995). A similar approach may be possible in cell biology, and one of the aims of this paper is to suggest potential experimental approaches for this (see section 8). In addition to parameters affecting the kernels, the model contains two dimensionless parameters, namely D, which reflects the angular diffusion coefficient of the cells, and , which reflects the rate at which collagen is re-aligned by the cells. Experimental data is already available in the work of Guido and Tranquillo (1993) which enables the value of D to be estimated as 0.27; we derive this estimate in section 8 of the paper. However, the experimental procedure prevents estimation of the underlying timescales, so that the corresponding dimensional value cannot be estimated. In addition, we are unaware of existing data from which can be estimated, although we suggest appropriate experimental approaches in section 8. 7 3 Analysis of the Model 3.1 General Properties Solutions of the evolution equations (6) and (7) with the associated boundary conditions (4) and initial conditions which satisfy equations (5) have the following properties: 1. If f and c, the fibroblast and collagen densities, are initially non negative, they remain so. This is true provided that the solutions are continuous and differentiable. 2. The total density of collagen and fibroblasts remains constant in time, i.e., 2 0 f(t, ) d = 1 and 0 c(t, ) d = 1 (11) provided f and c are integrable with respect to for each t and ct and ft are continuous. 3. The constant solution f = 1 2 and c = 1 is a steady state. This corresponds to a fully isotropic representation of fibroblasts and collagen fibers. 4. When the fibroblast s diffusion coefficient D = 0, solutions in which all the collagen and fibroblasts are oriented in parallel are also steady states; of course the fibroblast orientations can be in either of the directions parallel to the collagen. More precisely, by considering weak solutions when D = 0 showthat a steady state solution is (f( ),Copyright )) = ( ( a), ( a)) where is the Dirac distribution. This requires that W 1 (0) = 0 and W2(0)W 3 (0) = 0, which are expected intuitively (see section 2) and are satisfied by the kernels in (9, 10). If in addition, we assume that W2 and W3 are periodic, then ( ( + a), ( a)) and ( ( a) + ( + a), ( a)) are steady state solutions also. In some cases, series of delta like functions, corresponding to several isolated orientations, are also steady states; this is discussed in more detail later in the paper (see section 6). 5. Solutions where the delta functions are symmetric with respect to one another are steady states when there is no diffusion. In other words, f = k i=1 ai (x bi) where bi [0, 2 ] can be a steady state solution if for each i, f(bi x) = f(bi + x), with similar conditions for c. These types of solutions are encountered in section 6.3 where they are referred to as type II solutions. 8 3.2 Linear Stability As a first step to understanding the behavior of the system we examine the stability of the constant steady state solution u0 = (fo, co) = ( 1 2 , 1 ) with respect to angularly inhomogeneous perturbations. In order to do this we first linearize equations (6) and (7) about u0. Note that any angularly homogeneous perturbations would be steady state solutions except that those other than uo do not satisfy the normalization condition. By observing that Wi k = k, with all its derivatives zero, for i = 1,2 or 3 and k a constant, one can see that the linearized equations are: f t = D f f0 (W1 Copyright for [0, 2 ] (12) Copyright t = c0(W2 f0)( ) (W3 f)( ) for [0, ]. (13) We assume perturbations of the form f = k=0 e kt (ak cos k + bk sin k ) (14) c = k=0 e kt ( k cos k + k sin k ) (15) where f is 2 periodic and c is periodic, i.e., k = k = 0 for k odd. The periodicity of c implies that ct is also periodic which in turn implies that ak = bk = 0 for k odd. In order to meet the boundary conditions, only even wave numbers need to be considered. First we observe that if W and u are periodic with period T, then when evaluating W u it does not matter over which interval of length T the integral is taken. We chose the interval which is symmetric about zero. Simplification using trigonometric identities then gives W3 f(t, ) = k=0 k even e kt ak cos k W3(k) + bk sin k W3(k) (16) W1 c(t, ) = k=0 k even e kt k cos k W1(k) + k sin k W1(k) (17) where W (k) = a a W(x) cos kx dx. (18) When equations (14) and (15) are substituted into the linearized equations (12) and (13), the coefficients must satisfy the following equations: kak = k2Dak + k2f0 k W1(k) (19) k k = k2 f0c0ak W3(k) (20) 9 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 9 10 Growth rate Wave number 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Growth rate Diffusion a b Figure 1. The growth rate, k, predicted by the linear analysis is plotted for different wave numbers. In graph (a) the growth rate is plotted as a function of the wave number k when D = 1. In graph (b) the maximal growth rates, k = 4 (the dotted line) and k = 2 (the solid line), are plotted as a function of D, the diffusion coefficient. In both cases = 1 and the kernels are defined by a = 4 in equations 9. where ak = (ak, bk), k = ( k, k) and ak = k = 0 for k odd. For k even, ak = 0 and k = 0, k must satisfy the condition k = k2 2 D D2 + 4 c0f2 0 W 1(k) W3(k) . (21) Our assumption that W1 = 2W3 implies 2 W3(k) = W1(k). Thus for k even there is always one positive and one negative k. The mode with the largest positive k is expected to be dominant. Analytical comparison of the growth rates is not possible because of the complex form of W1(k), but by numerically computing k the wave number with the maximum growth rate can easily be found; typical results are shown in figure 1a. Since the odd wave numbers do not satisfy the boundary conditions, the even wave number with the maximum growth rate is 2. This is more clearly seen in figure 1b where 2 and 4 are plotted against D. The maximum dominant mode switches from 4 to 2 as D is increased and once they switch the growth rates remain very similar. Thus the homogeneous steady state is always unstable, and our linear analysis suggests that depending upon the diffusion coefficient, either one or two peaks should begin to evolve corresponding to either one or two orientations for the collagen. This is verified numerically in section 5, and the dependence on the diffusion coefficient is explored more fully later in the paper. By decreasing the range of influence that collagen and fibroblasts have on each other, which means decreasing the support of the kernels, the graph in figure 1a shifts to the right. In the case where D = 0, the mode with the maximum growth rate seems to be the maximum number of independent peaks 10 possible, subject to peak separation being greater than half the support of the kernel (see section 6). 3.3 Nonlinear Analysis Having determined the stability of the uniform steady state using linear stability analysis, we now consider the full nonlinear problem. By so doing we try to understand how the nonlinearities affect the system. Being motivated by previous work where patterns with two characteristic wave lengths are seen (Maini, 1990), we proceed by considering solutions of the form f = a0 + a2(t) cos 2 + a4(t) cos 4 + . . . (22) c = b0 + b2(t) cos 2 + b4(t) cos 4 + . . . . (23) Substituting these into the nonlinear equations (6) and (7), using equations (16) and (17) and integrating the equations after multiplying by 1, cos 2 and cos 4 gives a 0(t) = 0 (24) a 2(t) = 4Da2(t) + 4a0b2(t) W1(2) + 4a2(t)b4(t) W1(4) + . . . (25) a 4(t) = 16Da4(t) + 16a0b4(t) W1(4) + 4a2(t)b2(t) W1(2) + . . . (26) b 0(t) = 0 (27) b 2(t) = 2 b0a2(t) W3(2) a0 4a4(t) W3(4) 4 b2(t)a4(t) W3(4) a0 3a4(t) W3(4) + . . . (28) b 4(t) = 2 b0 a22 (t) W 2 3 (2) + 2a0a4(t) W3(4) b2(t)a2(t) W3(2) 2a0 3a4(t) W3(4) + . . . . (29) The same procedure could be carried out for a more general form of the solution but this suffices for our purposes. The interesting feature of these equations can be seen in equation (26) where the last term depends only on the mode 2 terms and similarly in equation (29) where the second to last term depends only on terms from lower modes. This means that in the full nonlinear model, initial conditions with a purely mode 2 term induce a small mode 4 contribution to the solution. In principle, this in turn induces a mode 8 term in the solution, etc. but because of constraints imposed by the conservation of mass and the lowgro wth rates of the very high modes, these high mode contributions are negligible. We confirm the growth of a mode 4 solution out of mode 2 initial data numerically in section 5. 4 Numerical Implementation Now that we have a basic analytical understanding of the system, we proceed to investigate its behavior in more detail using numerical methods. Before looking at simulations we describe the numerical method 11 used to solve the equations by first dealing with the convolutions and then with the equations. The convolutions are calculated by discretizing the domain with a uniform grid of mesh length h = N and using the left-hand rectangular rule. We denote half the length of the support of the kernel by se. Henceforth, any reference to W1, W2 or W3 means the discretized version, in which the equations are only evaluated at points on the spatial grid. In order to solve the model numerically, the equations are discretized using a centered difference formula while keeping them in conservation form. The time derivative is calculated implicitly using the trapezoidal rule. Thus the diffusive part of the system is simply a Crank-Nicolson discretization. The domain [0, ] is discretized into N points and the domain [0, 2 ] is discretized into 2N points. The resulting discretized nonlinear equations are solved using the software package nksol (Brown & Saad, 1987). This software uses an inexact Newton method to solve the nonlinear system with linear Krylov iterations used to approximate the Newton equations. We refer to solutions of the discretized equations as fn and cm where n [0, 2N 1] and m [0,N 1] denote the grid points (unless otherwise noted N = 100 in our simulations). One property of the true solutions f and c is that if they are initially non negative then they remain so. However, this is not the case for the numerical solutions fn and cm: w hen the densities become sufficiently aggregated and thus there is an abrupt transition from lowdensit y to high density, the numerical solutions become negative and begin to oscillate. This is a common problem in discretized systems with abrupt transitions (Osher & Chakravarthy, 1984). In our case the oscillations are reinforced due to the nature of the equations and quickly grow. Thus, the discretized system does not maintain an important feature of the system - non negative solutions. In the discretized system, the case corresponding to total alignment in one direction is a discrete version of the delta function, in defined to be in = 1 for i n = 0 0 otherwise , (30) where i defines the grid point at which the function is located, and n [0, 2N 1] denotes the grid point. These are the types of solutions we expect to develop in the model, based on our preliminary analysis and intuitive expectation. Yet these are precisely the type of solutions which are difficult to reproduce in a discretized system due to the discontinuities. To get around this problem, a common practice in advective schemes is the use of flux limiters which smooth the oscillations that occur near abrupt changes in the solutions (Thuburn, 1996; Sweby, 1984). We also adopt this approach albeit in a very crude manner. When the solution is about to become negative, the algorithm limits the flux to 12 keep the solution non negative. In practice this is implemented in the following manner: if fn (or cm) is negative, it is added to the larger of its neighbors, fn+1 or fn 1, and then it is set to zero. If its modified neighbor is negative, say fn+1, then fn+1 is also set to zero. The numerical solution nowremain non negative and still conserves mass. This not only acts as a crude flux limiter, but also makes sense in the context of the behavior of the underlying equations. Those equations have the property that when f = 0 then f is also zero, and our ad hoc flux limiter has the effect of reducing the numerical derivative at the point where the solution is set to zero. As always, but particularly in light of the unsophisticated flux limiter added to our algorithm, we interpret and accept our numerical results only as far as they are consistent with the analytical predictions of the model. 5 Numerical Confirmation of the Basic Properties of the Model In order to verify both the numerical scheme and the analysis of section 3 we look at perturbations about the homogeneous steady state ( 1 2 , 1 ). When solving the nonlinear system given by equations (6) and (7) with initial conditions corresponding to a mode 2-type perturbation about the steady state, one sees both the growth of this term, and the appearance of a mode 4 term as predicted in section 3.3. To ensure that the mode 4 term is due to the nonlinearity and is not being introduced through numerical errors, we solve the linear system (12) and (13), with the initial conditions again having purely mode 2 terms. In this case the mode 2 terms growand no other modes appear. As verification of our numerical scheme, we confirmed that the growth rate in the linear and non linear problem are consistent with the analysis. We nowconsider the effect of a single point perturbation to the steady state in both f and c at the same point. The modes with = 4 and = 2 initially grow, and depending on the diffusion coefficient D either two peaks form (small diffusion) or one peak forms (high diffusion). This again agrees with the analysis, although it is hard to determine if the resulting peak formation is due to the different growth rates of the modes or to the constraints imposed by the support of the kernels; one problem is that the two growth rates are very similar for most parameter values. When small random perturbations to the steady state are used as the initial conditions which maintain c(0, ) f(0, ) = f(0, + ), similar results are obtained either one or two peaks form depending on the value of the diffusion coefficient (see figure 2). The results from this section help us to understand the model in two ways. Depending on the diffu- 13 0 /2 Angle 0 5 10 15 20 25 30 Time 0 /2 Angle 0 50 100 150 200 250 Time a b Figure 2. Collagen density is plotted for two simulations with the same initial conditions and different fibroblast diffusion coefficients, illustrating that the collagen evolves to two distinct isolated orientations when the diffusion coefficient is small and to a single orientation when the diffusion coefficient is larger. In (a) D = 0 and the collagen evolves from an almost isotropic state to two distinct orientations; whereas in (b) where D = 0.5, the collagen evolves to a single orientation. This agrees with the results of the linear analysis; in fact in (b) the mode 4 terms (corresponding to 2 peaks on [0, ]) can be seen growing as predicted. In both simulations = 1 and the kernels are defined by a = 4 in (9). The initial conditions are random perturbations of magnitude 0.03 about the homogeneous steady state keeping c(0, ) f(0, ) = f(0, + ). The darker shading between the contour lines indicates higher density. 14 sion coefficient, the long term behavior of the collagen is to have either one or two isolated orientations. The other important result is that even with the ad hoc modification, the numerical scheme captures the basic features of the system predicted by the analysis. This gives us confidence in the numerical simulations, enabling further numerical experiments. 6 Patterns Generated fromSim ple Initial Conditions In order to understand more fully the types of behavior which the system can exhibit, we examine in this section the evolution of numerical solutions starting from simple initial conditions, with f c on the interval [0, ] in all cases. With this condition the fibroblasts and collagen initially have the same orientation, so that one intuitively expects the behavior to be less complicated than when they have different initial orientations; the latter case is considered in section 7. First we shall examine what happens when the collagen and fibroblast orientations are limited to one, two, or three isolated initial directions. We then address the case in which the collagen and fibroblasts have a continuous interval of orientations, and finally we look at solutions which arise from initial conditions with several local maxima. 6.1 Isolated Initial Orientations Of particular interest is the case when all the collagen and fibroblasts are ordered in a few directions. This is modeled by weighted delta functions located at each orientation, with the weights reflecting the fraction of density oriented at each of the angles. Starting with the simplest case in which everything is oriented in one direction, = ih, we have initial conditions of the form fn = a in and cm = d im. In the case D = 0, corresponding to no random reorientation in the fibroblast population, these solutions are steady states, as predicted by the analysis in section 3. As D is increased above zero, the long term solution for the fibroblast density becomes more uniform, eventually approaching a constant f( ) 1 2 ; the collagen density remains localized at a single grid point. Because the domain of f is [0, 2 ], we also consider the initial conditions for f of the type fn = a1 in + a2 i+N n (recall that Nh = ). This corresponds to having fibroblasts oriented both in the direction of the collagen and opposite to that direction. Biologically, one expects this scenario to behave exactly as that above, in which the fibroblasts are all oriented with the collagen, and this is confirmed in simulations. Moving to the next simplest case in which the collagen and fibroblasts have two orientations, = ih 15 and = jh, the initial conditions are fn = a in + a jn and cm = d im + d jm . In this case the two peaks move together until they merge and form a single peak in each variable (see figure 3a). This behavior is determined by the kernels W1, W2 and W3; in the simulations shown in figure 3, these are given by equation (9) with a = 4. This makes the support, se 1.28 radians for W1. Biologically this means that collagen with orientation = influences any fibroblast with orientation [ se, +se]. This influence brings the fibroblasts to an orientation closer to that of the collagen, and the influence of the fibroblasts on the collagen has a similar effect. If se is less than the separation of the two initially imposed peaks, the long range interaction does not extend from one peak to the other, making them independent of each other and one would expect the initial conditions to be a steady state. In figure 3a, the peak separation is about 0.94, which is less than se 1.28, and thus the peaks should influence each other resulting in their merging. Finally we consider the case of three initial localized orientations. If two of these orientations are symmetric about the third then this middle peak remains stationary. Each neighboring peak pulls on the middle peak with the same intensity resulting in no net change in position (see figure 3b). However, the two outer peaks are drawn into the middle, resulting in the eventual merging of all three peaks. There is of course a special case when the three peaks are symmetric about each other which results in an unstable steady state of type II which is more fully explained in section 6.3. This description depends critically on the initial outer peaks having the same intensities, and if the neighboring peaks have differing heights, their influence is no longer symmetric. In order to understand this we mention here two points. The first is that there may be up to three factors which influence a fibroblast peak - the random reorienting of fibroblasts due to a non zero D, the pull of the coinciding collagen peak due to the special initial conditions f Copyright and the pull of neighboring collagen peaks. The second point is that for most parameters, the collagen changes more slowly than the fibroblasts. This is because the rate of collagen remodeling depends on fibroblasts density, reflected mathematically by the flux term having an additional convolution in equation (7) compared to the flux term in equation (6). When neighboring peaks have different heights they pull each other with differing strengths. The smaller peaks of fibroblasts tend to reorient to the direction of the larger neighboring peaks of collagen, yet the collagen peak at the same orientation tends to keep the fibroblasts at their original orientation. If the coinciding collagen peak did not exist, there would not be any reason for the fibroblasts to maintain the original orientation and they would not remain localized, but would rather almost immediately reorient to the direction of the other collagen peaks. The result is that the smaller fibroblast peaks remain localized but they spread out in angle space as they slowly move towards the orientation of the larger 16 0 /2 Angle 0 2 4 6 Time 0 10 20 30 0 0 a 0 /2 Angle 0 5 10 15 Time 0 10 20 30 0 0 b Figure 3. The evolution of collagen densities when collagen and fibroblasts initially have (a) two and (b) three distinct, isolated orientations. In both cases the peaks in orientation merge to give a single isolated collagen orientation, with the fibroblast orientations localized around this. The collagen densities are plotted for simulations where the parameters are the same: D = 0, = 1 and in equation (9) a = 4 making se 1.28 radians. The initial conditions are f = c = 0 except for f30 = c30 = f60 = c60 = 50 in (a) and f30 = c30 = f60 = c60 = f90 = c90 = 100 3 in (b) giving a peak separation of about 0.94 radians. The fibroblast densities look similar. 17 peaks until all the peaks merge at a neworien tation. 6.2 Intervals of Initial Orientation Having confirmed and clarified the behavior of the model for discrete initial orientations, we now consider more complex scenarios. Specifically, we consider initial conditions representing collagen and fibroblasts which have a localized range of orientations. This is described by a hat function which we define on the grid as Hi,k n = 1 for k n i k 0 otherwise . (31) Thus, i denotes the grid point at the center of the hat, 2k+1 denotes the width of the hat, in grid points, and n denotes the grid point where the function is being evaluated. Starting from initial conditions fn = d1Hi,k n and cm = d2Hi,k m with kh< 2 , these initial orientations evolve into densities with all the collagen and fibroblasts oriented in one direction, = ih, described by in (see figure 4a). As with the previous simulations, this result is determined by the type of kernels used in the convolutions, with the support of the kernels being crucial. Recall that the support of the kernel, W1, determines the range of orientations over which collagen influences the fibroblasts, and vice-versa for W2 and W3. When this support is less than the support of the hat function (see figure 4b, c and d), two peaks form initially. This is a continuum version of the case discussed above, in which two peaks are placed symmetrically about a third. At the edge of the hat function the pull is only in one direction - towards the center of the hat function - while far enough into the hat function the pull is equal in both directions, causing no net change. By looking at the convolution of the initial conditions, it is easy to predict where these peaks will form. The convective term draws the densities toward the center of the hat function until the convolution plateaus. Peaks continue to form in the hat until it is filled with non interacting peaks. Peaks are non interacting depending on their separation and the support of the kernels. If the separation is too small, the kernel causes the peaks to interact and merge, otherwise they persist. This link can be easily seen in figure 5 where the support of the hat function, the length of the support of W1 or 2se, and the last time plot of the collagen density are shown. Knowing this, it seems reasonable to expect the wave number with the maximum growth rate to increase as the support of the kernel decreases as was found in section 3.2. Increasing the diffusion coefficient, D, causes the fibroblast density to spread out, rendering a broad range of fibroblast orientations. The main consequence is that peaks may merge (see figure 6b). As the diffusion is increased, the shape of the kernel becomes more important. If the kernel is steep enough 18 0 /2 Angle 0 0.5 1 1.5 2 Time 0 10 20 0 0 /2 Angle 0 0.2 0.4 Time 0 5 0 a b 0 /2 Angle 0 0.05 0.1 0.15 0.2 0.25 Time 0 2 4 6 0.2 0 /2 Angle 0 0.05 0.1 Time 0 1 2 3 T c d Figure 4. Collagen densities in simulations where the effect of varying the angular distance over which fibroblasts and the collagen influence one another is studied. As the distance is decreased, reflected by increasing a in equation (9), more peaks which are independent of one another are able to form. Halving the support almost doubles the number of peaks which are formed. In these four simulations, in order to alter this distance, only a is changed: in (a) a = 4 making se 1.28 radians, in (b) a = 20 making se 0.6 radians, in Copyright a = 75 making se 0.38 radians and in (d) a = 250 making se 0.16 radians. The initial conditions are f = v = H59,22 n with D = 0 and = 1. The fibroblast densities look very similar. These results are plotted in a different way in figure 5 19 0 /2 Angle 0 5 10 15 20 25 30 Density 0 /2 Angle 0 2.5 5 7.5 10 12.5 15 Density a b 0 /2 Angle -2 0 2 4 6 8 10 Density 0 /2 Angle -1 0 1 2 3 4 5 6 7 Density c d Figure 5. Collagen density plotted at one instant in time to showthe relationship between the peak separation and the range of influence between the collagen and the fibroblasts. The dark bar on top shows the length of the support of W1 or 2se and the grey bar at the bottom shows the support of the initial hat function. These simulations are the same as those shown in figure 4. 20 around zero then it keeps the fibroblast peaks very narrow, counteracting the effect of diffusion, but if the kernel is not steep near zero then the diffusion spreads the fibroblast density more readily. As the regions of high fibroblast density become close, they merge into one peak. This occurs to a greater extent as D is increased, and helps to explain why the wave number with the maximum growth rate depends on D, i.e., as D is increased the dominant wave number decreases (see section 3.2). The more spread out the density, the more slowly the peaks form. Finally, at very high diffusion coefficients, peaks in f which are a distance away from those imposed initially are formed (see figure 6c). This is due to the periodic nature of the kernels and of c. If the diffusion causes enough of the density to spread into the interval [ , 2 ], it starts to form peaks corresponding to those in the interval [0, ], provided that the diffusion is not strong enough to level them out. Biologically this corresponds to fibroblasts traveling in either direction with the same orientation as the fibers. 6.3 A Continuumof Orientations with Several Maxima Having established the evolution from a fewdiscrete initial orientations as well as the behavior of an initial localized continuum of orientations, we now look at initial conditions which are a mixture of the two - continuous, with several local maxima. Specifically, we consider initial conditions where fn and cm are proportional to 1 + cos knh and 1 + coskmh respectively; recall that n and m denote mesh points, whereas k and h represent the wave number and the angular step size respectively. Here k must be even to satisfy the periodic boundary conditions for c. Our simulations suggest that these initial conditions evolve into two different types of solutions: delta functions which are independent (see figure 7b) and delta functions which are symmetric with respect to each other (see figure 7a). We refer to these as type I and type II, respectively. When D = 0, the long term behavior is either type I or type II. In the previous simulations all the solutions have been of type I, where the peaks are separated by a distance larger than se making them independent of each other. Type II solutions were mentioned briefly in section 3. The solutions of type II are unstable steady states and as such do not persist biologically. However, studying the evolution of type II solutions from initial data gives valuable insight, which is particularly useful in section 7 where transients similar to type II solutions persist for a long time. Type II solutions are formed with k 2 peaks when k is fairly small: in the case a = 4, we observe such solutions when k < 10. The same factors as before determine the solution type, namely the range of the collagen and fibroblast interactions and the separations of the peaks in their initial densities. If this separation 21 0 Angle 2 0 0.05 0.1 0.15 0.2 Time 0 2 4 6 8 2 0 T 0 Angle 2 0 0.1 0.2 0.3 0.4 Time 0 2 4 2 0 T a b 0 Angle 2 0 1 2 3 Time 0 0.2 0.4 0.6 .8 2 0 T .2 .4 6 Copyright Figure 6. An illustration of the effect of increasing the angular diffusion coefficient on the distribution of fibroblast orientations. As the diffusion coefficient is increased, the orientations become more spread out around the localized collagen orientations. In (a) D = 0 and the fibroblasts end up having four isolated orientations (this is the same simulation as figure 4c where the collagen density is shown). As the diffusion coefficient is increased in (b) to D = 0.1 the peaks no longer represent isolated orientations since they are spread over several grid points and the middle two have merged leaving only three peaks. Increasing D further in Copyright to D = 1.0 shows the formation of three additional peaks, at a distance away from the first set representing fibroblasts orienting in the direction opposite the collagen fibers. The initial conditions and parameters (other than D) are the same as in figure 4c. 22 0 /2 Angle 0 10 20 30 40 50 60 70 Time 0 /2 Angle 0 20 40 60 80 100 Time a b 0 /2 Angle 0 2 4 6 8 10 12 Time 0 /2 Angle 0 5 10 15 20 Time c d Figure 7. An illustration of type I (a,Copyright and type II (b, d) solutions evolving from initial conditions with several orientation maxima. We plot contour lines of the collagen density with darker shading indicating higher densities. As this shows, the solutions change from type II in (a) when there are four local maximum in the initial conditions, or k = 8, to type I in (b) when there are 5 local maximum, or k = 10, and a = 4 in equation (9). By increasing a to a = 10, the range of influence of the collagen and the fibroblasts is decreased and the transition from type II solutions to type I solutions occurs at higher k values, between k = 12 and k = 14. In Copyright k = 10 nowforms a type II solution and in (d) k = 14 forms a type I solution with three independent peaks. The initial conditions are for both f and c are proportional to sin kx + 1. In all cases D = 0 and = 1. The fibroblast densities look similar. 23 is sufficiently large, it allows a quick consolidation to delta functions separated by equal distances. Since the peaks are distributed about each other symmetrically and are of equal heights, they form an unstable steady state (see section 6.1). When the diffusion coefficient is increased from zero, the value of k where the transition from simulations which develop type II solutions to simulations which develop type I solutions will decrease. This is due to the fact that the random reorientation keeps the maximum from coalescing as quickly and compactly, which allows the noise (introduced via numerical errors) to destroy the symmetry. Solutions of type I are stable steady states and are therefore biologically attainable and more significant. Such solutions are formed when k is sufficiently large: for a = 4 we observe these solutions when k > 8, with the solutions always having two peaks. This latter observation is because when a = 4, se 1.28 for W1, which determines that at most two peaks can form on the interval [0, ]. However, because k is large in these cases, the initial densities are sufficiently spread out that it takes a long time for the localization to occur. Mathematically, the explanation for this is that the convolutions of the initial conditions with the kernels are sufficiently smooth that they have very small derivatives, which thus concentrate the density only very slowly. This is so slow that the small asymmetries, arising from the discretization, perturb the system, leading to the formation of type I solutions. In applications, of course, any symmetric initial conditions would be perturbed by natural fluctuations. By decreasing the range of influence of collagen and fibroblasts, more peaks in type I solutions form and type II solutions should form from initial conditions where the density is more evenly distributed. Changing a from 4 to 10 verifies this by causing k = 10 to form five peaked solution of type II (see figure 7c) and k = 14 to form a three peaked solution of type I (see figure 7d). 6.4 Conclusions All of the types of initial conditions considered so far have resulted in solutions where the collagen densities are concentrated at discrete, isolated orientations and the fibroblast densities are localized around these discrete orientations, with the degree of aggregation increasing as the angular diffusion coefficient goes down. There are two types of solutions which have been observed: type I solutions where the peaks are independent of one another, and type II solutions where the solution is symmetric about each peak. Only type I solutions are physically relevant since type II solutions are unstable, to asymmetric perturbations, developing into type I solutions. The independence of the peaks in type I solutions is due to the half-range of influence of the collagen and fibroblasts being less than the 24 separation of the peaks. Thus the support of the kernel, which corresponds biologically to the range of directions over which the collagen and fibroblasts are able to reorient one another, determines the maximum number of peaks that can occur in type I solutions. 7 Interaction between Fibroblasts and Collagen Previously the initial conditions have been chosen such that the collagen and fibroblast densities were proportional to each other, in order to reduce the extent of their interaction, and thus simplify the results. Now we change this strategy and investigate the interaction between the fibroblasts and the collagen. When their initial configurations are different, the parameter becomes very important. Recall that determines howstrongly the collagen is reordered. Figure 8 illustrates the results from a simulation in which the collagen is initially all set at one angle and the fibroblasts at another. When is small the collagen s initial conditions are more important in determining the final solution, and the fibroblasts reorient in the direction of the collagen (figure 8c). When is increased sufficiently, the situation is reversed and the fibroblasts initial conditions become more important, causing the collagen to reorient in the direction of the fibroblasts (figure 8d). At intermediate values of both the fibroblast and the collagen orientations alter, stabilizing at some intermediate orientation (figure 8a). Increasing the diffusion coefficient simply makes the fibroblast density more spread out, as seen in figure 8b. If the initial conditions are changed to hat functions, similar results are obtained, but with some additional complications. The hat functions can become split depending on their width and the separation. Figure 8e illustrates the case when the initial densities for fibroblasts and collagen are both hat functions, but with different weights. As expected, the hat function with more density draws the collagen to a greater extent than the function with less density. They all merge to form single coinciding peaks. The final type of simulation that we have used to understand the interaction of fibroblasts and collagen have initial conditions with fn proportional to 1+sin jnh and cm proportional to 1+sinkmh. Here we are using the same notation as above, but the wave numbers j for the variable f and k for the variable c can be different. In these simulations the collagen has several orientation maxima which are different from the orientation maxima of the fibroblasts. The short term behavior of the solution depends on : if is sufficiently small, then the fibroblasts initially try to reorient to the form of the initial collagen density (see figure 9a), and vice-versa if is large (see figure 9b). However, these initial reorientations do not persist long enough to significantly alter the long term behavior of the collagen and 25 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time a ( = 3, D = 0) b ( = 3, D = 0.1) 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 0 /4 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time Copyright( = 0.1, D = 0) d ( = 30, D = 0) e ( = 20, D = 0) Figure 8. Fibroblast and collagen densities in simulations which demonstrate how they influence each other. Initially, the fibroblasts and collagen have separate orientations. The evolution depends on the parameters , which determines how strongly the fibroblasts reorient the collagen, and D, the angular diffusion coefficient of the fibroblasts. The simulations shown in (a), Copyright and (d) differ only in the value of which is used. In (a) where = 3, the fibroblasts and collagen roughly have the same influence on each other causing the final orientation to be at an intermediate value from the initial ones. In Copyright where = 0.1, the fibroblasts are drawn to the orientation of the collagen, whereas in (d) where = 30, the collagen is drawn to the orientation of the fibroblasts. By comparing (b) with (a) the effect of the diffusion coefficient is shown. The only difference between the simulation shown in (a) and (b) is that the diffusion coefficient is changed from D = 0 to D = 0.1 respectively. In (b) the fibroblasts have a greater range of orientations. In (e) the collagen is drawn closer to the orientation of the greater mass of fibroblasts. Here = 20 and the initial conditions are hat functions with different heights, that is f 2H5,2 n + H40,2 n and c H20.5,5 n . (The non integer value in the superscript for H indicates that the hat function is centered between grid points.) The region for the collagen density when c > 1 is shaded, and the contour for f = 1 is a bold shaded line. Unless otherwise stated, D = 0 and the initial conditions are c 35 n and f 15 n . 26 only moderately influence the fibroblasts. The variables reorient in a manner largely determined by the initial conditions; this can be seen in figure 9. In the context of scar tissue formation, this suggests that the initial deposition of collagen is the most important factor in determining its long term orientation. More specifically, the model shows that the remodeling of the collagen by the fibroblasts alters the collagen orientation in a manner which is less dependent on the properties and initial conditions of the fibroblasts than on the initial conditions of the collagen. Yet in scar tissue formation, where the fibroblasts produce and degrade collagen (processes which are purposely not included in this model) the initial conditions of the fibroblasts may be crucial to the initial deposition of collagen. Despite this, once the collagen alignment is established and possibly just initiated, the conclusion of our model holds and remains valid for the period of collagen remodeling which continues for months post wounding (Mast, 1992). These results are consistent with the observation that anti-scarring therapies such as application of TGF s are only effective in the very early phase of wound repair (Shah et al. , 1994), despite the fact that scar remodeling occurs on a relatively long time scale. Restating, this model indicates that the parameter is important for the initial behavior of the solutions, but the initial density distributions are more significant for the final form of the solution. It is interesting to note that in these simulations, the fibroblasts and collagen can form isolated orientations which persist for long transient periods. This is contrary to behavior described in section 6.3, where the initial maxima did not become compact before coalescing. The fact that the peaks here do become compact explains why they persist for long periods. Their influence on each other is in the tail of the kernel, making it small until they move closer. This type of solution could be biologically significant, for example if the fibroblasts become inactive during this long transient. A final situation to consider is when both fibroblast and collagen densities have random initial conditions, set by randomly choosing a value between zero and one for each grid point. Those values are then rescaled so that fn and cm satisfy the normalization condition. As expected from previous considerations, we have found in a large number of simulations that the long time behavior involves the maximum number of peaks which can form in a type I solution (not illustrated for brevity). Changes in the parameter simply change the time it takes for the peaks to form; if is small then f forms peaks very quickly, while if is large then c forms peaks very quickly. If the diffusion coefficient is increased, the f solution smooths out and the c solution takes longer to form peaks. 27 0 /2 Angle 0 10 20 30 40 50 60 Time 81% of the fibroblast density divided between this peak and the corresponding peak in the interval to 2 . Similarly 19% of the fibroblast density is divided between this peak and that on to 2 . The two collagen peaks shown below eventually merge into one. 0 /2 Angle 0 0.2 0.4 0.6 0.8 1 Time 96% of the collagen density at this peak 4% of the collagen density at this peak a b Figure 9. Collagen and fibroblast densities from simulations which illustrate that initial conditions have a greater influence than the parameter on the final solutions. In (a) although = 0.1 causing the fibroblasts to try and quickly reorient in the direction of the collagen, the final solution has most of the fibroblast density at two peaks corresponding to the fibroblasts initial condition with two maxima, or k = 2 (only half the fibroblast domain is shown). The long term behavior of the collagen is two peaks with roughly half the collagen density in each, again a consequence of the initial conditions for collagen with four maxima, or k = 8. In (b) = 150 causing the collagen to try and take the orientation of the fibroblasts, but again the initial conditions prevail and most of the collagen density is oriented at one peak corresponding to the initial conditions with k = 2. The long term behavior of the fibroblast density is four peaks, two have about 75 percent of the total density divided equally between them and the rest is divided equally between the remaining two peaks. The initial conditions for the fibroblast has four maxima, or k = 8. In (a) the shaded region shows where the collagen density is greater than 0.4 and the dark line is the contour for fibroblast density 0.1 and in (b) the shaded region shows where the collagen density is greater than 0.5 and the dark line is the contour for fibroblast density 0.27. Both simulations have D = 0 and initial conditions where f and c are proportional to sin k + 1. 28 8 Discussion In this paper we have developed a simple model for fibroblast and collagen alignment interactions and investigated its behavior using a combination of analytical and numerical techniques. The results have yielded a number of insights into the alignment process, and we begin the discussion by considering the application of these insights to wound healing and cancer. Dermal wound repair is currently a very active research topic, due to recent advancements which promise to lead to newclinical techniques for reducing scarring (McCallion & Ferguson, 1996). Despite a large volume of experimental research (see Clark, 1996a for review) and some mathematical modeling (Olsen et al. , 1995; Dale et al. , 1996), many details of the process remain poorly understood. It is known that collagen alignment plays a key role in the healing process; in fact, collagen alignment is one method for characterizing scar quality. In humans and other tight skinned animals, collagen has a cross weave structure in normal tissue, whereas in scar tissue it is aligned parallel to the plane of the skin (Harmon et al. , 1995; Welch et al. , 1990). As a dermal wound is repaired, fibroblasts replace the provisional matrix of fibrin with a collagen matrix (Clark, 1996b). The collagen is then reorganized for months by the fibroblasts until at some time they become quiescent and the matrix remains relatively unchanged (Mast, 1992), corresponding mathematically to a steady state. In our model, the configuration observed in normal skin can be represented by the collagen density having two peaks of orientation, with roughly 90o separation, and scar tissue can be represented by a collagen density profile with one alignment peak. Both of these solutions are indeed steady state solutions of our model. Moreover the model predicts that both of these steady states can be stable for the same parameters. Thus the properties of the fibroblasts, characterized by the parameters D, and the kernels W1, W2 and W3, need not be changed in order to obtain either type of solution. Rather, our model predicts that it is the initial conditions which determine which of the steady states form. This is consistent with biological observations that transient application of growth factors can permanently alter the quality of repair (Shah et al. , 1994). Cancer invasion is another area of application in which the key process is the movement of a cell population through a collagen-dominated extracellular matrix. For most cancers, the relevant cell type is transformed epithelial cells, but these share many phenotypic similarities with fibroblasts, and our model is quite applicable in this context. In tumor invasion, the role of cell and matrix orientation has received relatively little attention, with recent work focusing instead on the details and interaction of protease production, directed cell movement and altered cell adhesion; a reviewof recent experimental 29 work in this area is given in Jiang & Mansel, (1996), and modeling approaches are described in Byrne & Chaplain, (1996) and Perumpanai et al. , (1996). Determination of the role of collagen reorientation during invasion is an important modeling challenge for which our work lays the foundations. A number of investigators have studied the interaction between fibroblasts and collagen using in vitro experiments involving fibroblasts in collagen gels (Clark et al. , 1995; Guido & Tranquillo, 1993; Stopak & Harris, 1982). Our model relates directly to this type of experiment, and our results suggest a series of experiments that could be performed in order to estimate the model parameters: 1. The first and most fundamental experiment is to introduce fibroblasts into a collagen gel in which all the collagen is oriented in a single direction. This procedure has in fact been performed by Guido & Tranquillo, (1993), and their data can be used to estimate the dimensionless diffusion coefficient D, as mentioned in section 2. Specifically, Guido & Tranquillo present a histogram of fibroblast orientations, which can be conveniently summarized by the standard deviation away from the mean, which is approximately 0.39 radians (the mean is of course the predominant collagen direction). Simulations of our model imply that in a uni-directional collagen network, the standard deviation of the fibroblast orientations is an increasing function of D, as expected intuitively, and the standard deviation of 0.39 radians implies that D 0.27. Strictly, this is only an upper bound on D, because in the experiments, the collagen deviates to some extent from being uni-directional, to an extent that cannot be determined quantitatively. However, model simulations showthat fibroblast distribution is in fact relatively insensitive to this deviation in collagen distribution, suggesting that D 0.27 is a reasonable estimate. Unfortunately, the experimental procedure in Guido & Tranquillo, (1993) means that only the steady state is considered, so that no dimensional information is available. However, if it were possible to modify their procedure and introduce the cells after aligning the gel, this type of information would be accessible. 2. Further experiments would require extensions of the procedure of Guido & Tranquillo, (1993), so that the gel has two (or more) isolated collagen orientations. The natural approach, in keeping with the development of this paper, is to begin by constructing experimental initial conditions in which there are two different isolated collagen orientations, with fibroblast directions localized around these as in the simulations for figure 3a. Measurement of the time evolution of the distribution of fibroblast directions could be compared with model simulations, to delineate D and W1. We suggest that this may be achievable by obtaining two separate uni-directional collagen gels with corresponding cell populations, as in (1) above, and then placing one above the other, at an angle. 30 3. As in the theoretical development in the paper, the next experimental step would be to construct experimental initial conditions in which fibroblasts and collagen had different orientations. One possible approach to achieving this would be to juxtapose two collagen gels, with the collagen unidirectional in one, and oriented in two different isolated directions in the other. The fibroblasts could be placed in the uni-directional gel, from where they would enter the other gel with one predominant direction of motion. Measurement of the long-term distribution of collagen and fibroblast orientations near the interface, compared to model simulations, would then enable , W1, W2 and W3 to be determined. For instance, if is very large the collagen in the dualdirectional gel should become uni-directional, whereas if is small the fibroblasts should become oriented in both directions of the gel. Although the spatial aspect, which is ignored in our model, will certainly play a role, in this region the local dynamics should dominate. The work in this paper provides a theoretical framework on which more complex models could be built. There are three important modifications that could naturally be made to the model in order to make it more widely applicable. The most important is to add a spatial component, and study the way in which orientations develop as the fibroblasts move spatially. We are particularly interested in the transitions between two regions with different characteristic patterns of orientation, such as from scar tissue to normal tissue. Such an extension would greatly increase the complexity of the model, requiring the collagen and fibroblast densities to be a function of time, orientation, and two spatial coordinates. The second important modification would be to add a term to the model representing production of collagen by fibroblasts. This is relevant in a number of applications: for example, scar tissue has a greater density of collagen than does normal tissue, due to production by fibroblasts entering the wound (Shah et al. , 1992). Finally, we would like to extend our model to three space dimensions by using two angular variables. Conceptually this is straightforward, but the resulting model would be computationally much more intensive, and experimental verification of the two-dimensional framework is an essential precursor to this extension. Acknowledgements. We thank Philip Maini (Oxford), Mark Ferguson (Manchester) and Markus Owen (Warwick) for helpful discussions. This work was supported by grant GR/K71394 from the EPSRC, and by a grant from the London Mathematical Society (scheme 3). 31 References Alberts, B., Dennis, B., Lewis, J., Raff, M., Roberts, K., & Watson, J. D. 1994.Molecular Biology of the Cell. 3 edn.Garland Publishing.Chap. Cell Junctions, Cell Adhesion, and the Extracellular Matrix, pages 978 986. Besseau, L., & Giraud-Guille, M. M. 1995. Stabilization of Fluid Cholesteric Phases of Collagen to Ordered Gelated Matrices. J. Mol. Biol., 251(2), 197 202. Birk, D. E., & Trelstad, R. L. 1986. Extracellular Comparments in Tendon Morphogenesis: Collagen Fibril, Bundle, and Macroaggregate Formation.J. Cell Biol., 103, 231 240. Bray, D. 1992.Cell Movements.Garland. Brown, P. N., & Saad, Y. 1987.Hybrid Krylov Methods for Nonlinear Systems of Equations.Tech. rept. Lawrence Livermore National Laboratory. Byrne, H. M., & Chaplain, M. A. J. 1996. Modelling the role of cell-cell adhesion in the growth and development of a carcinoma.Math. Comp. Modelling, 24, 1 17. Civelekoglu, G., & Edelstein-Keshet, L. 1994.Modelling the dynamics of F-actin in the cell.Bull. Math. Biol., 56, 587 616. Clark, R. A. F. (ed). 1996a.The Molecular and Cellular Biology of Wound Repair. 2 edn.Plenum Press. Clark, R. A. F. 1996b.Wound Repair Overviewand General Considerations. Pages 3 50 of: Clark, R. A. F. (ed), The Molecular and Cellular Biology of Wound Repair, 2 edn.Plenum Press. 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. Cook, J. 1995.Waves of alignment in populations of interacting, oriented individuals. Forma, 10, 171 203. 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. Deutsch, A. 1995.Towards analysing complex swarming patterns in biological systems with the help of lattic-gas cellular automata.J. Biol. Systems, 3, 947 956. Edelstein-Keshet, L., & Ermentrout, B. G. 1990.Models for contact-mediated pattern formation: cells that form parallel arrays. J. Math. Biol., 29, 33 58. 32 Elsdale, T. 1973. The generation and maintenance of parallel arrays in cultures of diploid fibroblasts. In: Kulonen, E., & Pikkarainen, J. (eds), Biology of Fibroblast.NewY ork: Academic Press. Geigant, E., Ladizhansky, K., & Mogilner, A. 1997. An integro-differential model for orientational distributions of F-actin in cells.SIAM J. Appl. Math.To appear. Grunbaum, D. 1994. Swarming behaviour as an Aide to Chemotaxis. In: Parrish, J., & Hammner, W. (eds), 3D Animal Aggregations.Cambridge University Press. Grunbaum, D. 1997. Advection-diffusion equations for generalised tactic searching behaviors. J. Math. Biol. in press. 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. Harmon,Copyright B., Zelickson, B. D., Roenigk, R. K., Wayner, E. A., Hoffstrom, B., Pittelkow, M. R., & Brdoland, D. G. 1995.Dermabrasive scar revision - immunochemical and ultrastructural evaluation. Dermatol. Surg., 21, 503 508. Jiang, W. G., & Mansel, R. E. 1996.Progress in anti-invasion and anti-metastasis research and treatment. Int. J. Oncology, 9, 1013 1028. Maini, P. K. 1990.Superposition of modes in a caricature of a model for morphogenesis.J. Math. Biol., 28, 307 315. Mast, B. A. 1992.The Skin.Chap. 22, pages 344 355 of: Cohen, I. K., Diegelmann, R. F., & Lindblad, W. J. (eds), Wound Healing Biochemical and Clinical Aspects.Saunders. McCallion, R. L., & Ferguson, M.W. J. 1996.FetalWound Healing and the Development of Antiscarring Therapies for Adult Wound Healing. Pages 561 600 of: Clark, R. A. F. (ed), The Molecular and Cellular Biology of Wound Repair, 2 edn.Plenum Press. Mogilner, A., & Edelstein-Keshet, L. 1995.Selecting a common direction I. Howorien tational 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. Murray, J. D. 1993.Mathematical Biology. 2 edn.Springer.Chap. 9, pages 241 244. 33 Neubert, M. G., Kot, M., & Lewis, M. A. 1995. Dipersal and Pattern Formation in a Discrete-Time Predator-Prey Model.Theoretical Population Biology, 48(1), 7 43. Okuba, A. 1986.Dynamical aspects of animal grouping: swarms, schools flocks, and herds.Adv. Biophys., 22, 1 94. Olsen, L., Sherratt, J. A., & Maini, P. K. 1995. A Mechanochemical Model fo Adult Dermal Wound Contraction and the Permanence of the Contrated Tissue Displacement Profile. J. Theor. Biol., 177, 113 128. Osher, S., & Chakravarthy, S. R. 1984. High resolution schemes and the entropy condition. SIAM J. Numer. Anal., 21(5), 955 984. Perumpanai, A. J., Sherratt, J. A., Norbury, J., & Byrne, H. M. 1996. Biological inferences from a mathematical model for malignant invasion.Invasion and Metastasis, 16, 209 221. Pollard, T. D., & Cooper, J. A. 1986.Actin and actin-binding proteins. A critical evaluation of mechanisms and function.Ann. Rev. Biochem., 55, 987 1035. Shah, M., Foreman, D. M., & Ferguson, M. W. J. 1992.Control of scarring in adult wounds by neutralising antibody to transforming growth factor .Lancet, 339, 213 214. Shah, M., Foreman, D. M., & Ferguson, M. W. J. 1994. Neutralising antibody to TGF- 1,2 reduces cutaneous scarring in adult rodents. J. Cell Sci., 107, 1137 1157. Sherratt, J. A., & Lewis, J. 1993. Stress-induced alignment of actin filaments and the mechanics of cytogel.Bull. Math. Biol., 55, 637 654. Stevens, A. 1995.Trail following and aggregation of myxobacteria.J. Biol. Systems, 3, 1059 1068. Stopak, D., & Harris, A. K. 1982.Connective Tissue Morphogenesis by Fibroblast Traction.Devel. Biol., 90, 383 398. Sweby, P. 1984.High Resolution Schemes using Flux Limiters for Hyperbolic Conservation Laws.SIAM J. Numer. Anal., 21(5), 995 1011. Thuburn, J. 1996.Multidimensional Flux-Limited Advection Schemes.J. Comp. Phys., 123, 74 83. Welch, M. P., Odland, G. F., & Clark, R. A. F. 1990. Temporal Relationships of F-actin bundle formation, collagen and fibronectin matrix assembly, and fibronectin receptor expression to wound contraction. J. Cell Biol., 110, 133 145. 34