SIAM J. APPL. MATH. c 2000 Society for Industrial and Applied Mathematics Vol. 61, No. 2, pp. 506 527 A MATHEMATICAL MODEL FOR SPATIALLY VARYING EXTRACELLULAR MATRIX ALIGNMENT JOHN C. DALLON AND JONATHAN A. SHERRATT Abstract. Orientation of extracellular matrix fibers in the skin is a keyingredien t of tissue appearance and function, and differences in fiber alignment are one of the main distinctions between scar tissue and normal skin. In this paper, the authors develop a mathematical model for alignment of collagen fibers and the fibroblast cells that remodel them; the model extends previous work in which spatial variation was excluded. Numerical simulations of the model are presented, which show spatial variations in alignment over long transients, but with spatiallyuniform behavior in the long term. This is investigated further via asymptotic analysis, using the angular diffusion coefficient as a small parameter. This method enables calculation of the form of the steadystate orientation peaks observed numerically; by considering behavior at large times, the rate of approach to these peaks is shown to be exponential. Extension of this analysis to the spatially varying model confirms that long-time behavior will be spatiallyuniform except in one special, and biologicallyunrealistic, case. The authors conclude that the spatiallyv arying alignment patterns observed in skin are in fact slow transients, and biological implications of the modeling are discussed. Key words. Alignment, integrodifferential equations, perturbation theory, collagen, fibroblast, wound healing. AMS subject classifications. 45, 92 1. Introduction. Alignment and orientation phenomena abound in the natural world. These include ecological, biological and physical processes such as animal herds, schools of fish [18], dense fibroblast aggregates [10], actin filaments [22], extracellular matrix [25] and molecular collagen structures [2]. Understanding the factors which influence alignment allows manipulation and alteration of the system s structure, thereby influencing different factors such as survivability of a population or properties of a material. Because of its significance in applications, there have been many attempts to model and understand alignment; for example in animal herds [5], in fibroblasts [14, 15], in actin filaments [11, 3, 24] and in the extracellular matrix [1, 20]. This paper is concerned with the alignment of the extracellular matrix of skin and connective tissue which are composed primarily of the protein collagen. In these tissues it forms a fibrous network that is created and maintained by cells called fibroblasts. During wound healing many complicated processes interact to repair the wound. One of the first events is the formation of a blood clot, composed of the protein fibrin, in the wound space. Within 24 to 48 hours the fibroblasts in surrounding skin begin to invade the fibrin clot, dissolving and replacing it with new tissue as well as contracting the wound region. Thus the fibrin clot is the start of a provisional matrix that changes as the wound heals, and is eventually replaced by a fibrous matrix composed mainly of collagen. This new tissue, known as a scar, is different from the surrounding tissue. One characteristic of scars is a greater degree of alignment than normal tissue [9, 4]. By understanding the interaction between the cells and the fibrous matrix, it may be possible to manipulate the network produced in wound repair, and thereby to reduce the extent of scarring. Indeed, the first anti-scarring This work was supported byEPSR C grant GR/K71394 Department of Mathematics, Brigham Young University, Provo, UT 84602-6539 U.S.A. Department of Mathematics Heriot-Watt UniversityEdin burgh EH14 4AS U.K. 1 2 J.Copyright Dallon and J. A. Sherratt therapies are currently being used in clinical trials, and their refinement depends on a detailed understanding of collagen fiber reorientation by fibroblast cells. Because we are focusing on this one process which occurs throughout wound healing, our model does not correspond to any particular temporal phase of healing. However, a great many other processes occur in parallel with matrix production and remodelling by fibroblasts, including mechanical contraction, blood vessel ingrowth and epidermal repair. These have all been the subject of previous modeling studies, and form part of a body of work whose long-term result may be a coordinated model of the whole of wound healing. There have been earlier attempts to study fibroblast alignment using the same type of modeling we use here, but they were concerned with interactions between fibroblasts, in dense cell populations, rather than interactions between fibroblasts and the matrix on which they move [16]. In other work, the interactions between fibroblasts and the matrix are considered using very different formulations: with the matrix being modeled as a viscoelastic substance [1], through the use of reaction diffusion equations [20, 19] and incorporating a discrete formulation to model the cells [7]. All of these efforts have led to a greater understanding of alignment and the interplay between fibroblasts and the extracellular matrix which will impact wound healing and artificial tissue engineering. Previously we developed a spatially uniform model for collagen orientation [6] based on a system of integrodifferential equations. Here we extend that work by including spatial variation into the model. Although several authors have worked on alignment problems, only a few have considered spatial variation using continuous angular and spatial variation [15, 1]. Adding the spatial component involves including one to three extra dimensions to the problem, decreasing the feasibility of generating numerical solutions. In addition, the new spatial interactions also complicate any analytical treatment. Of the works cited above, [1] concentrates on tissue mechanics in a continuum dynamics framework and [15], which we discuss further in 3, consider spatial interactions which are very different from ours. The rest of the paper is organized in the following manner. First we describe the original spatially uniform model and then extend the model to include variations in space. Following that, in 4 we numerically investigate solutions of the extended model. In 5 and 7 we analyze first the spatially uniform model and then the extended model. Finally we conclude with a discussion of the implications and potential applications of the work. 2. The Spatially Uniform Model. In this model we consider collagen and fibroblasts located at the same point in space which vary in orientation. Their respective densities are represented by c(t, ) and f(t, ) where t is scaled time and is the angle of orientation with respect to an arbitrary reference direction; c and f must be periodic in . Because the collagen forms fibers, which are assumed to be smooth, we take its orientation to be in the direction of its tangent line. This direction can be represented by either or + giving c a period of , whereas f is 2 periodic. We phenomenologically model only three aspects of the system: angular diffusion of the fibroblasts, contact guidance and reorientation of collagen. All of these contribute to the angular flux of either the fibroblasts or the collagen. First, we allow the fibroblasts to reorient due to their random motion. This is modeled with a diffusive term in angle space, using a constant diffusion coefficient. Second, we model the ability of fibroblasts to align themselves in the direction of the fibers making up the substrate [12], a phenomenon called contact guidance . This is done by including the flux term f (W1 Copyright which causes fibroblasts to move up angular gradients of a weighted A spatiallyv arying model of alignment 3 average of the collagen density. The weighted average is given by convolving c with a kernel W1 in the standard way defined by: (W1 Copyright( ) = =0 (2.1) W1( Copyright ) d . Note that this integration is over one period of c. Since the fibroblasts and the collagen are at the same location in space, they will influence each other even if their orientations are different. The convolution represents this long range angular interaction in the system. Finally, as the third aspect of the system, we model the ability of the fibroblasts to alter the collagen fibers with a collagen flux term similar to that in the fibroblast equation. The use of an advective term to describe the remodeling of collagen by fibroblasts is simply a first approximation, made in the absence of detailed information about what is happening at the cellular level. In reality the process is a complex one, in which fibroblasts degrade and produce collagen as well as exerting mechanical forces on the fibers, resulting in their reorientation. For simplicity we assume no net production and consider all the remodeling as reorientation. The model is given by: f t = random orie n t at ion D f orientation by collagen f (W1 Copyright (2.2) Copyright t = orientation by fibroblasts c(W2 f) (W3 f) (2.3) where D and are positive constants. In the convolutions in (2.3), the integrals are over one period of f, that is between 0 and 2 . Note that the coefficient for the second term in equation 2.2 is taken to be one without loss of generality, by a rescaling of time (see [6] for details). The parameter is a measure of the relative strength of how effectively the fibroblasts reorient the collagen compared to how effectively the collagen aligns the fibroblasts. The orientation terms in the two equations clearly differ, with the reorientation rate of collagen depending on both fibroblast density and its derivative, whereas the reorientation rate of fibroblasts depends only on the derivative of collagen density. This is because the orientation of collagen by fibroblasts is an active process and depends on how many fibroblasts are doing the remodeling. In reality, both terms should also be multiplied by coefficients depending on local fibroblast and collagen densities, but we neglect this in the absence of any experimental data that would determine the appropriate forms. As one expects intuitively, the properties of the kernels are key in determining the behavior of the solutions of 2.2 and 2.3. IfW 1 (0) = W 2 (0) = W 3 (0) = 0, isolated peaks of collagen and fibroblasts can be steady states. These conditions are implied by the natural biological assumptions that the kernels are differentiable and even functions. Only the magnitude of the difference between the orientation of the collagen and the fibroblasts should be important, not the sign, or in other words the kernels should be even functions. In addition, the direction of the fibroblast along the fibers should not be important, and thus we assume W2 and W3 are periodic functions. Due to 4 J.Copyright Dallon and J. A. Sherratt the periodicity of c, W1 is also assumed to be periodic. For the majority of our numerical simulations we will take W1 = 2W2 = 2W3 = W where W( ) = Ce a 2 for 2 1 a ln Copyright 0 for 1 a ln C < 2 < 2 4 the periodic extension otherwise, (2.4) a = 4, = 5 10 5 and 1 Copyright = 2 2 (2.5) W(x) dx. The constant C and the factors of 2are chosen so that the area under the kernels, over the period of c for W1, and of f for W2 and W3, is 1. The affects of altering the kernels and other results are discussed in detail elsewhere [6]. We now summarize some of those results in order to give an understanding of the behavior of the solutions and thus provide the background to the new material presented in this paper. Other than the uniform steady state, which is unstable, the key type of steady state solution is one where the collagen and fibroblasts align predominantly in coinciding isolated peaks at discrete orientations. When orientation peaks are close enough, specifically within half of the support of the kernels, they will interact with each other, eventually merging to form one peak at an intermediate orientation. The two parameters in the system, D and affect the transient behavior but have little impact on the the qualitative long term behavior of the system. Of course as D increases, the fibroblasts become less localized about the orientation peak. The location and number of orientation peaks that form is heavily influenced by the initial distribution of the collagen and fibroblasts. 3. Extended Model Including Spatial Term. Because fibroblasts move as they reorganize the extracellular matrix and because the orientation of the extracellular matrix can vary in space, we extend our model to include spatial variation in two dimensions. In order to do this, we let f and c depend on x, a point in the plane, as well as orientation and time i.e., f(t, , x) and c(t, , x). Because we assume collagen is not being spatially moved but is only being reoriented, there is no spatial flux for collagen and equation 2.3 remains the same. The fibroblasts on the other hand do move actively, and equation 2.2 must now include a term for the spatial flux. The velocity of the cells is given by v = s(c) v, where s is the speed and v is a unit vector in the direction of motion. We assume the speed is a function of collagen density and can thus allow for the biologically realistic case of no cell motion when there is no collagen present. The direction of motion is determined by the cell s orientation so that v = (cos , sin ). The spatial flux is fv giving f t = angular flux D f f (W1 Copyright x spatial flux due to cell m oti on fs(c) v (3.1) in place of equation 2.2. We do not include a spatial diffusive term since biologically diffusive reorientation is likely to be the dominant cause of the spatial spread. The other work, similar to ours, considering spatial effects of orientation [15] does so in a very different manner. In that work the authors consider three models with A spatiallyv arying model of alignment 5 different spatial interactions. In the first, the spatial component includes diffusion and alignment due to spatial interactions between objects. This spatial interaction is modeled by having a multi-dimensional kernel integrated over space and time. Thus objects move in space via diffusion and in order to obtain a preferable spacing with one another. In the second model the authors include an additional spatial flux where the objects move up concentration gradients in order to aggregate, and in the final model, which is a simplification of the second, the objects instantaneously move up the gradient. They then consider the linear stability of the homogeneous solutions and compare the results with cellular automata models. As stated earlier, we incorporate the long range angular interaction in the same manner as the other authors using convolution terms, but in our model we consider interactions between two differing populations, cells and matrix. Our spatial interaction is very different and does not involve the kernels or any long range spatial attractions or motion towards dense aggregates. It is simply determined by the active motion of cells in the direction in which they are oriented. 4. Numerical Solutions. In order to numerically solve equations 3.1 and 2.3 we keep the angular flux in conservation form and use a centered difference approximation. The spatial flux is discretized, with upwinding according to the direction in which the cells are moving [13]. Other numerical methods were investigated for the spatial flux term, including a Lax Wendroff method. This resulted in less dispersion for spatial waves of fibroblasts, but was computationally less efficient. Since we are primarily interested in steady states and the long term behavior of the solutions, we opted to use the upwind method. We solve the system using Euler s method and limit the angular flux with the method described in previously [6]. Because of computational limitations, unless otherwise noted we consider a spatial domain where the solutions are homogeneous in one direction, allowing us to solve the problem in one space dimension instead of two. 4.1. No Flux Boundary Conditions. First we consider the case where there is no spatial flux at the boundaries of the spatial domain. (The angular domain has periodic boundary conditions as stated earlier.) This means that (4.1) fs(c)v d = 0 for x = 0, L. There are several pointwise boundary conditions which will satisfy equation 4.1. One interpretation of this boundary condition is a physical boundary such as the edge of a petri dish where the fibroblasts are prevented from leaving the domain. In this case, the cells will alter their direction of motion and remain in the domain. Most experiments focus on fibroblast motion away from boundaries in order to minimize their effects. Due to the lack of reported information about cellular behavior at physical obstacles, we consider several different boundary conditions which seem intuitively reasonable. 4.1.1. Reversing Cells. First we assume the that cells reverse and re-enter the domain at an angle 180 degrees different from their previous orientation. The pointwise condition satisfying equation 4.1 and representing the situation described above is (4.2) f(t, , x)s(c(t, , x))v ( ) = f(t, + , x)s(c(t, + , x))v ( + ) 6 J.Copyright Dallon and J. A. Sherratt for x = 0, L. At first look this condition does not seem sufficient to determine the solutions, but since the orientation determines the direction of cell motion for any given , only one set of boundary conditions is necessary, i.e. for [ 2 , 2] the fibroblasts are all moving to the right and only the boundary conditions at x = 0 are necessary. By considering periodic solutions of f (c is periodic and v( ) = v( + )) equation 4.2will be satisfied and it is consistent with the biological condition of having the cells reverse at the boundary. Biologically, periodic solutions of f are expected since there is no bias in the direction in which the cells line up and move along the collagen fibers. Thus we are assuming that the fibroblasts are localized in two peaks separated in orientation by . The following simulations use the boundary condition (4.1), but first we mention two other cases which satisfy equation 4.2that are not reversing cells. The first is when s = 0 at the boundary and the cells stop moving and the second consists of solutions where the fibroblasts at both boundaries are oriented at either 2 or 3 2 and move parallel to the boundary. In both of these special cases the cells end up aggregating to the boundaries and staying there. Figure 4.1 shows a solution when the collagen is initially isotropic and uniformly distributed in the domain, and the fibroblasts are localized between x = 0.5L and x = 0.8L with orientation 1.3 . As time proceeds the fibroblasts remain localized in both space and orientation moving across the domain. When they arrive at the boundary they reverse and begin moving across the domain in the opposite direction (at a new orientation a distance , in angle space, away from the previous orientation). In addition the fibroblasts are orienting the collagen into the direction of their motion. Depending on initial conditions, solutions can develop transient spatial patterns but they eventually become uniform in space with the collagen and fibroblasts forming density peaks at isolated directions. This holds true for a wide variety of initial conditions and parameter sets. Figure 4.2sho ws contour plots for the collagen density at the beginning of a simulation, an intermediate time and a long time solution that appears to be the steady state configuration for a typical simulation. The final angle at which the collagen becomes oriented, when the initial conditions are of the type shown in figure 4.2, is plotted in figure 4.3. When the collagen is given two initial orientations, 1 and 2, the final orientation is approximately the weighted average of the two initial orientations, i.e. a 1 + (1 a) 2 where a is the proportion of the spatial domain where the collagen has the orientation 1. The final orientation of the collagen does not seem to depend on the parameters and is insensitive to the choice of the speed function s. We ran several tests varying the parameters and using different functions for s which were non-decreasing with respect to c and obtained the same results. Of course the transient behavior varies with parameters, and the time taken to arrive at the steady state solution is highly dependent on the values of and s. Figure 4.4 shows the collagen density at an intermediate time for different parameters and speed. There is a pathological case, if the shorter interval between the two initial orientations, 1 and 2, includes 2 . In this case the cells will, in the one dimensional situation we consider, stop moving in space when they are oriented at 2 . Thus as the fibroblasts reorient from 1 to 2 or vice versa the majority will, as they orient in the direction 2 , stop moving, aggregate in space and dramatically slow the evolution of the system towards a uniform steady state. 4.1.2. Other no flux boundary conditions. There are other boundary conditions which would satisfy the no flux requirement given in equations 4.1. We will briefly discuss three other possibilities. The first is when the cells move parallel to A spatiallyv arying model of alignment 7 collagen density 0 o/r4ienta ti/o2n 3 /4 0 0.2 0.4 0.6 0.8 1 x/L 0 3 6 fibroblast density 0 o/r2ientat ion 3 /2 2 0 0.2 0.4 0.6 0.8 1 x/L 0 3 6 9 12 a b x/L f d as t increases Copyright Fig. 4.1. The collagen and fibroblast densities in (a) and (b) are shown at t = 4.5. The fibroblast have not yet traversed the entire spatial domain, thus the collagen is only oriented predominantly in one direction where the fibroblasts have been. The dotted line at the bottom is the contour for collagen (a) or fibroblast (b) density equal to 1. In Copyright the total fibroblast density is plotted with respect to x with two adjacent lines representing the solution at two times differing by 0.1 units. As the fibroblasts reflect offthe boundary the density at one space point increases as the fibroblasts turn around and start moving in the other direction. So near the boundary the fibroblasts are localized into two groups which spatially overlap but have different orientations. The collagen is initially uniformly distributed in the domain and the fibroblasts are localized at one orientation, 1.3 and to part of the spatial domain, x = 0.5L to x = 0.8L. The parameters used are D = .02, = 0.3, L = 1 and s = 0.2. For the discretization we use ht = 0.01, hx = L 50 and h = 50 as the time, space and angular step sizes respectively. 8 J.Copyright Dallon and J. A. Sherratt 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation a b 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation Copyright Fig. 4.2. The contour plots of the collagen density are shown for the initial conditions in (a), at t = 10 in (b) and at t = 700 in Copyright. The collagen is initially oriented in two directions 1 = 0.6 and 2 = 0.9 . In the interval x = 0 thru x = aL where a = 0.6 the collagen is oriented at angle 1 and for the rest of the interval it is oriented at angle 2. The fibroblasts are uniformly distributed in the domain. The parameters used are D = 0.01, = 2.0, L = 1 and s = 0.3. The discretization is the same as figure 4.1 except ht = 0.01. the boundary, which could result in encapsulation of regions such as occurs when connective tissue encapsulates organs [25, 8] or in tumour encapsulation [23, 17]. Another possible no flux boundary condition would have the cells randomly alter their direction when they encounter the boundary. The final no flux boundary condition we mention is one where the cells reflect off the boundary in a manner similar to the reflection of light off a mirror. This condition gives rise to the following equation: f(t, , x)s(c(t, , x)) v( ) = f(t, , (4.3) x)s(c(t, , x))v ( ). for x = 0, L. In order for solutions with localized peaks of orientation to satisfy this boundary condition and be a steady state, the peaks should be evenly spaced in the A spatiallyv arying model of alignment 9 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Final 1 Fig. 4.3. The relationship between the final orientation of the collagen and one of the two initial orientations, 1, is plotted. The black squares are data from several different simulations and the solid line is a 2+(1 a) 1 with a = 0.5 and 2 = 3.08. The initial conditions for the simulations are the same as in figure 4.2 with a = 0.5, 2 = 3.08, 1 is varied for each simulation and with the same parameters and discretization. The simulations were run up to time 1000. In all simulations the collagen density was concentrated at either one or two adjacent angular grid points for all space. The final angle was taken to be the one angle or the average of the two adjacent grid points. 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation a b Fig. 4.4. The contour plots of collagen density are shown for different simulations at time 10. In (a) = 0.6 and in Copyright s = 0.3c 0.1+c ; otherwise the initial conditions and parameters are the same as those for figure 4.2. The contour lines are for densities of 0.05, 0.5, and 5. interval 0 to . In this manner they are placed symmetrically with respect to one another and the influence of a neighboring peak aligning the densities to its orientation is balanced by another neighboring peak on the opposite side. Intuition as well as numerical simulations suggest that solutions of this type are unstable if there are more than two orientational peaks of collagen [6]. Thus the only stable configurations seem to be with the collagen oriented at 2 , 0 or a combination of the two. 10 J.Copyright Dallon and J. A. Sherratt 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation 0 0.2 0.4 0.6 0.8 1 0 /2 3 /2 2 x/L orientation a b Fig. 4.5. The contour plots of the collagen (a) and fibroblast densities (b) are shown at time 1000. The boundary conditions are given by f = 0 except at 1 for x = 0 and 2 for x = L. In (a) and (b) the solution seems to be tending towards a spatially uniform steady state except near the imposed boundary. The initial conditions are given by setting the fibroblast density to zero except at one orientation on each spatial boundary with 1 = 0.69 and 2 = 1.82. The collagen is initially uniform in space and orientation. The contour lines are for densities of 3 and 10. The parameters are the same as those given for figure 4.2 with the discretization differing only by ht = 0.001. If 1 and 2 differ by 2 then two isolated peaks of collagen will form at each of those angles. This is due to the size of the support of W (see [6]). 4.2. Dirichlet Boundary Conditions. The last boundary conditions we will mention are of Dirichlet type. Again, because the orientation determines the direction of the cells at the boundary, we need only specify half the boundary values i.e., f(t, , 0) = g1( ) for [ 2 , 2 ] f(t, ,L) = g2( ) for [ 2 , 3 2 ]. In figure 4.5, collagen and fibroblast contours are shown for simulations where g1 and g2 are zero for every point except one. This corresponds to having a constant density of fibroblasts at the boundary aligned in one direction on each boundary. As can be seen in the figure, the solution is tending to a profile with all the collagen and fibroblasts aligned in one direction except at one of the boundaries, where the alignment changes in a small spatial region to the imposed value. The alignment which dominates is not affected by the parameters but is determined by the imposed boundary conditions. If the two orientations given at the boundaries are 1 and 2 and if | cos 1| > | cos 2|, then 1 will be the predominant orientation. The reason for this is the cells with the largest horizontal component of motion out compete the other cells in reorienting the collagen. All these results indicate that the behavior of the fibroblasts at barriers could be fundamental to the type of alignment that arises. From our numerical simulations it seems that steady state solutions which vary in space are not easily formed. In the next sections we will support this with some analytical results. First we start with the spatially uniform case, developing an asymptotic solution which we will then extend to the spatially varying model. A spatiallyv arying model of alignment 11 5. Asymptotic Analysis for the Spatially Uniform Model. We know that when D = 0, peak solutions of equations 2.2 and 2.3 can be steady states [6]. Intuitively we expect the fibroblast density to spread out and be less localized as D increases, and numerically this is confirmed, with the behavior of the solutions not being sensitive to small variations in D. In order to continue with our analysis we will follow the method used in [16] where the authors consider a singular perturbation problem with a small diffusion coefficient by scaling the variable, assuming solutions of the form u(z) = exp (z) (5.1) and looking for an expansion in powers of . Because we already know that when D = 0, solutions which are localized to one orientation are possible steady states, we also look for an asymptotic expansion of the solution when D is small. This case, when the diffusion coefficient is small when compared to one ie., D = 1, means that the orienting effect of the collagen is much greater than the disorienting effect of the random diffusion. From experiments [12], an upper bound on the diffusion coefficient has been estimated to be 0.27 [6]. 5.1. Smoothform of c. We look for an asymptotic expansion of a steady state solution centered at = 0. We make the natural assumption that both f and c are smooth and differentiable; in fact, we will show that these are inappropriate and will revise the calculation in 5.2. The total mass of fibroblasts and collagen is conserved, and we let f d = Mf and c d = Mc. First we set the left hand side of equations 2.2 and 2.3 to zero and integrate giving f = f (5.2) (W1 Copyright + K( ) A( ) = c( )(W2 f)( ) (5.3) (W3 f)( ). We then rescale and the variable of integration in the convolutions, y, so that = and y = y . In addition, the supports of f and c shrink as 0, but we require f and c to remain constant, so we let F( ) = ( )f( ) and C( ) = ( Copyright ). Making these changes of variable gives F = 1 ( )F 2 2 (5.4) W1( ( y))C(y) dy + ( )K( ) ( )A( ) = C( ) ( )2 W2( ( y))F(y) dy (5.5) W3( ( y))F(y) dy . The solutions of F and C are periodic, but we are considering solutions that become more highly localized in orientation space as gets smaller. If we ignore the periodicity and extend the functions beyond their periods in such a manner that they decay like exponentials then F(y) dy = (5.6) F(y) dy + o( n) 12 J.Copyright Dallon and J. A. Sherratt for any n, with a similar equation holding for C. Expanding Wi about 0 and recalling that W i (0) = 0 due to symmetry gives Wi( ( y))F(y) dy = Wi(0) (5.7) F(y) dy + 1 2 32 W i (0) ( y)2F(y) dy + O( 2) Differentiating equation 5.7 with respect to gives Wi( ( y))F(y) dy = 32 W i (0) (5.8) F(y) dy 32 W i (0) yF(y) dy + O( 2) for i = 1, 2and 3. Similar equations hold when F is replaced by C, with the limits of integration changed to 2 . We now expand F, C and the various constants as power series in 12 , with F( ) = F0( )+ F1( )+ ,Copyright ) = C0( )+ C1( )+ , A( ) = A0+ A1+ and K( ) = K0 + K1 + . Substituting 5.7 and 5.8 into equations 5.4 and 5.5 and collecting the first order terms gives F0( ) + O( 2) = 32 ( )F0( )W 1 (0) (5.9) C0(y) dy 32 ( )F0( )W 1 (0) yC0(y) dy + 1 ( )O( 2) + ( )K0 + ( )O( ) ( )A0 + ( )O( ) = 2 2( )C0( )W2(0)W 3 (0) (5.10) F0(y) dy F0(z) dz zF0(z) dz + 1 2( )O( 2). The distinguished limit here is ( ) = and ( ) = , which by changing variables gives Mf = f( ) d = F( ) ( ) d = (5.11) F0( ) d + o( n) for any positive integer n (a similar equation holds for C) and for the first order terms we have the following two equations F0( ) = F0( )W 1 (0) Mc yC0(y) dy (5.12) + K0 A0 = C0( )W2(0)W 3 (0)Mf Mf zF0(z) dz (5.13) . A spatiallyv arying model of alignment 13 Solving these we get F0( ) = K0 exp W 1 (0) Mc 2 2 Q d + B (5.14) exp W (0) Mc 2 2 Q C0( ) = A0 W2(0)W 3 (0) [Mf R] (5.15) where Q = yC0(y) dy and R = yF0(y) dy. The variable c represents a density which is always positive, but equation 5.15 changes sign when = R/Mf unless A0 = 0. In either case our analysis has given a result which is not biologically relevant. This is due to our assumed form of c. By changing our assumptions on c, below we derive a similar form for f which is relevant. 5.2. Dirac distribution for c. By allowing D, the angular diffusion coefficient, to be nonzero we assumed the solutions of f become smooth and differentiable. Thus we did not need to treat them as distributions and consider weak solutions of the equations. The problem with the analysis in 5.1 is that we assumed the same to be true of c, that is we assumed c also became smooth when D was nonzero. In fact, neither the equations nor the numerical simulations justify this assumption and the analysis in 5.1 suggests that it is not true. Here, we again look for an asymptotic expansion of a steady state solution centered at = 0. We assume, as suggested by numerical simulations, that f is smooth but that c is a Dirac distribution weighted by Mc, the total mass of the collagen. Again the mass of fibroblasts is conserved and we let f d = Mf . Since c is the Dirac distribution we cannot, as before, simply set the left hand side of equation 2.3 to zero and integrate. This time we must consider weak solutions i.e., 0 = 0 ( )(W2 f)( ) (W3 f)( ) (5.16) ( ) d where is smooth and compactly supported, in other words a test function. Simplifying gives 0 = Mc(W2 f)(0) (W3 f)(0) d d (5.17) (0) which implies 0 = (W2 f)(0) (5.18) (W3 f)(0). We assume that f and W3 are even functions; we have chosen W3 to be even and we are assuming that f is symmetric about its peak, since the cells have no orientational bias. Taking the integral in the convolution over an interval symmetric about zero gives (5.19) (W3 f)(0) = 0 and the steady state condition given by equation 2.3 is satisfied. 14 J.Copyright Dallon and J. A. Sherratt Next we consider equation 2.2 by setting the left hand side equal to zero and integrating to get f = McfW 1 (5.20) + K( ). Now, as in 5.1, we rescale so that = and let F( ) = ( )f( ), which gives dF d = McF d d (5.21) W1( ) + ( )K( ) Expanding W1 about 0, differentiating and recalling that W 1 (0) = 0 gives d d W1( ) = W 1 (0) + O( 3 (5.22) 2 ). Again, expanding as a power series in 12 gives F( ) = F0 + F1( ) + F2( ) + and K( ) = K0 + K1 + K2 + . Substituting into equation 5.21 gives F0( ) + O( 3 2) = F0( )McW 1 (5.23) (0) + ( )K0 + ( )O( ) Again, the distinguished limit is ( ) = which gives F0( ) = K0 exp W 1 (0)Mc 2 2 d + B exp W (0)Mc 2 2 (5.24) . For a given , f is bounded so we take K0 = 0. Since f = Mf the change of variables shows that Mf = F = F0 + F1 + (5.25) F2 + implying that F0 = Mf and Fi = 0 for i = 1, 2, 3, . Thus B = Mf W (0)Mc 2 and to leading order we have f Mf W (0)Mc 2 exp W (0)Mc 2 2 (5.26) . Collecting the 32 terms gives F1 as F1 = C + B McW 1 (0) 6 3 exp(W (0)Mc 2 2 (5.27) ). However, for the kernel W1 we are considering, defined in equation 2.4, W 1 (0) = 0 and the integral condition from equation 5.25 gives C = 0. So we collect the 2 terms, giving F2 = D + B McW(iv) 1 (0) 24 4 exp(W (0)Mc 2 2 (5.28) ). A spatiallyv arying model of alignment 15 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Fibroblast density Orientation 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Fibroblast density Orientation a: = 0.5 b: = .1 Fig. 5.1. The asymptotic expansions given in equation 5.26 (solid line) and 5.29 (dotted line) are plotted and can be compared to numerical solutions of equations 2.2 and 2.3 (+) for two values of . The numerical solutions are solved with a spatial step size of 0.0157 radians and a time step of 0.0001. Here, the integral conditions 5.25 require D = BW (iv) 1 (0) 8(W (0))2Mc . Thus f Mf W (0)Mc 2 1 W(iv) 1 (0) 8(W (0))2Mc +McW(iv) 1 (0) 4 24 (5.29) exp W (0)Mc 2 2 where W(iv) 1 denotes the fourth derivative of W1. Figure 5.1 shows a comparison of the numerically calculated solution and the first two non zero terms of the asymptotic expansion shown in equation 5.29 for two values of . Provided is fairly small, the asymptotic expansions agree very well with the numerical solutions. 6. Temporal Evolution of the Steady State. When extending our asymptotic expansion to the spatially varying model (equations 3.1 and 2.3), it will be necessary to consider the way in which f and c approach their equilibria at large times. Therefore, we now consider this approach in the spatially independent model; we consider only leading order terms as 0 and t . Since c( , t) Mc ( ) as t we hypothesize that at large times, Copyright , t) = Mc 1 (t) exp 2 (t) (6.1) to leading order, where (t) 0 as t . Then Copyright W1( ) = 2 2 W1( Copyright ) d = Mc 1 (t) 2 2 W1( ) exp 2 (t) d 16 J.Copyright Dallon and J. A. Sherratt Expanding W1( ) in a power series about = 0, and applying Watson s Lemma, gives Copyright W1( ) = Mc 1 (t) W1( ) exp 2 (t) d W 1 ( ) exp 2 (t) d + 1 2W 1 ( ) 2 exp 2 (t) d + = McW1( ) +Mc 1 4 (t)W 1 ( ) + O( (t)2). Substituting this into 2.2 then gives f t = 2f 2 f McW 1( ) +Mc 1 4 (t)W 1 ( ) + O( (t)2) (6.2) In view of the leading order form given in 5.26 for f at t = , we expect a large time solution of the form f = Mf (t) 2 W 1 (0)Mc 12 exp 2 (t) 2 W1 (0)Mc (6.3) where (t) 0 as t . Then f t = (t)f 1 2 (t) 2 W 1 (0)Mc 1 + 2 (t) 2 W 1 (0)Mc 2 f = 2 f (t) 2 W 1 (0)Mc 1 2f 2 = 2f (t) 2 W 1 (0)Mc 1 + 4 2f (t) 2 W 1 (0)Mc 2 . Substituting these into 6.2giv es 1 2 (t)f (t) 2 W 1 (0)Mc 1 + 2 (t)f (t) 2 W 1 (0)Mc 2 = 2 f (t) 2 W 1 (0)Mc 1 + 4 2f (t) 2 W 1 (0)Mc 2 +2 f (t) 2 W 1 (0)Mc 1 McW 1 ( ) +Mc 1 4 (t)W 1 ( ) + O( (t)2) f McW 1 ( ) +Mc 1 4 (t)W(iv) 1 ( ) + O( (t)2) . Since f is transcendentally small away from = 0, we only require that this equation be satisfied near = 0, which gives, to leading order, (t) = 4 + 2McW 1 (0) (t) 2 W 1 (0)Mc (t) = 2McW 1 (0) (t) (t) = exp (2McW 1 (0)t)B0 A spatiallyv arying model of alignment 17 for some constant B0. Note that this decays to zero as t since W 1 (0) < 0. Therefore the equilibrium solution form for f is approached exponentially quickly as t . Turning now to the c equation 2.3, we have f Wi( ) = MfWi( ) +Mf 1 4 (t) 2 W 1 (0)Mc W i ( ) +O (t) 2 W 1 (0)Mc 2 , for i = 2, 3, applying Watson s Lemma to the convolution integral as above. Therefore f W2( ) [f W3( )] = MfW2( ) +Mf 1 4 (t) 2 W 1 (0)Mc W 2 ( ) MfW 3( ) +Mf 1 4 (t) 2 W 1 (0)Mc W 3 ( ) +O (t) 2 W 1 (0)Mc 2 . Thus to leading order at large times Copyright t = M2 fW2( )W 3( Copyright (6.4) . Substituting equation 6.1 into this gives (t)c 2 (t) + (t) 2c (t)2 = M2 f [W 2 ( )W 3( ) +W2( )W 3 ( )] c + M2 fW2( )W 3 ( )2 Copyright (t) Again, since c is transcendentally small away from = 0, this gives to leading order (t) (t) = 2 M2 f [W 2(0)W 3 (0) +W2(0)W 3 (0)] (t) = C0 exp 2 M2 f [W 2 (0)W 3 (0) +W2(0)W 3 (0)] t , for some constant C0. Thus c also approaches its final equilibrium exponentially, but at a rate that is, in general, different from that for f. 7. Asymptotic Analysis for Extended Model. In skin, the collagen matrix forms a cross weave pattern where the collagen fibers orient in primarily two directions, approximately orthogonal, with the directions varying gradually in space. Thus we are interested in finding solutions which are peak-like in orientation space, but where the orientation peak varies in space; the numerical simulations discussed in 4 did not show any solutions of this type. To consider the possibilities analytically, we look for solutions of the type found in the asymptotic expansion for the spatially homogeneous model in 5, but with the peak varying with space i.e., f = A exp W 1 (0)Mc 2 ( g(x))2 (7.1) (7.2) c = Mc [ g(x)]. 18 J.Copyright Dallon and J. A. Sherratt By changing variables so that = g(x) and y = x and setting the left hand side to zero, equation 3.1 becomes 0 = 1 f f d d W1( ) y (fs(c) v) 1 yg (7.3) (fs(c)v ) . Expanding W1 in a Taylor s series as in 5 and simplifying gives 0 = y (fs(c) v) 1 yg (fs(c) v) + O( ) = fs(c) y v 1 fs(c) yg v 1 f s(c) + f s(c) yg v + O( ) = 1 f s(c) + f s(c) yg v) + O( ), which in turn simplifies to 0 = f W 1 (0)Mc s(c) + s(c) (7.4) yg v + O( ). At least one of the three factors in the second term must be zero for the solution to be a steady state. Our assumed form of f is never zero; we want spatially varying solutions so we assume that yg is not everywhere zero; g is independent of and v depends on so yg is not orthogonal to v for all . Thus we are left with the only remaining possibility that the product s(c)f is constant. This obviously makes the spatial flux zero since in equation 3.1, the spatial flux will then become the spatial gradient of a constant multiplied by v, which only depends upon . At t = , c is a delta function, and thus the condition that s(c)f is constant does not give a meaningful relationship between s(c) and f. However, at large times, we can use the leading order solution forms calculated in 6 to extract such a relationship. This is possible when the solution forms are such that f can be written as a function of c. Comparing equations 6.1 and 6.3, this is possible only when = 0 and (t) (t) is a constant, that is when McW 1 (0) = M2 f [W 2 (0)W 3(0) +W2(0)W 3 (7.5) (0)] Then, in this special case when (t) (t) = 0, f = Mf M 0 Copyright 0 1 2 (t) 0 1 (7.6) 2 c 0 . The constant 0 will depend on initial conditions rather than parameters. If conditions are such that 0 = 1, then f = Mf Mc (7.7)Copyright and thus the requirement for a spatially varying solution becomes s(c) 1/c. For s(c) of this form, then, our analysis predicts that a range of initial conditions will tend to A spatiallyv arying model of alignment 19 spatially varying peak-like solutions of the form in equation 7.1. However, this special case is biologically unrealistic, indicating that the fibroblasts decrease in speed when there is collagen present. Although fibroblast speed probably does decrease at high collagen densities, it will also decrease at low densities (becoming zero when there is no collagen for the cells to move on), with maximum speed at an intermediate value [21]. 8. Discussion. We are interested in finding solutions with alignment that is not uniform in space. This is particularly relevant for the skin, where the orientation of the extracellular matrix varies in space. Our numerical simulations suggest that solutions tend to a spatially uniform steady state. Other than the special case mentioned above when s(c) 1/c, the analysis performed also predicts that as D 0 the solutions eventually become uniform in space. Our analysis only considers a particular solution form, with isolated peaks in orientation space and with the spatial variation only entailing a shift in the orientation peak. Although this is a very natural form to study, many other forms of spatial variation are clearly possible. However, none are found in our numerical simulations. One possibility we have not thus far considered is spatial variation of coefficients. The extracellular matrix in the skin is a very complex structure which could easily alter the effective diffusive ability of a chemical, for example via the chemical binding more readily in some regions of the matrix. In addition, , which is a measure of the ability of the fibroblasts to orient the matrix, could easily vary in space due to more cross linking of the matrix in regions or higher concentrations of the various cytokines which affect the fibroblasts. Our investigation of this found that spatially varying the diffusion coefficient D and does not seem to give rise to spatial variation in the solution. Spatial variations in D and do generate alignment patterns at intermediate times, but the long term behavior is not changed. The most interesting transients are obtained by varying D. As mentioned in 2, the isotropic or uniform steady state is unstable. For low values of the diffusion coefficient, the linear analysis suggests, and numerical simulations confirm, that two peaks of collagen orientation should grow and for higher values of D only one peak should emerge [6]. By setting D to be low in one half of the domain and large in the other half, two peaks form transiently in one region with one peak in the other, but the final orientation is still spatially uniform (figures 8.1 and 8.2). Once the peaks begin to grow, cells with an intermediate orientation move into the portion of the domain with two orientational peaks, causing these two peaks to merge and coalesce. When the skin is damaged, a host of processes are initiated which usually result in the injury being repaired so that the new tissue is different from normal tissue, i.e. scarring. There are several differences between scar tissue and normal tissue including a lack of hair follicles and sweat glands, different proportions of collagen types, a more aligned collagen matrix, and a denser collagen structure [9]. These differences can have the adverse consequences of weaker and less functional tissue, as well as cosmetic differences. We conclude that for wound healing and tissue regeneration the transient behavior of the solutions presented in this paper are most biologically relevant. This is due to the very complex and dynamic nature of skin and the wound healing process. During tissue regeneration the properties of the fibroblasts are altered and other processes such as wound contraction and inflammation interact with the system. This causes the fibroblasts to be activated and then, after a period of months, to return to a relatively inactive state [4]. Any transient behavior of the solutions in this time frame would be the final outcome for practical purposes. In addition, as 20 J.Copyright Dallon and J. A. Sherratt 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation 0 0.2 0.4 0.6 0.8 1 0 /2 3 /2 2 x/L orientation a: collagen, t = 40 b: fibroblasts, t = 40 0 0.2 0.4 0.6 0.8 1 0 /4 /2 3 /4 x/L orientation 0 0.2 0.4 0.6 0.8 1 0 /2 3 /2 2 x/L orientation Copyright collagen, t = 4000 d: fibroblasts, t = 4000 Fig. 8.1. Contour lines for the collagen (a) and Copyright and fibroblast densities (b) and (d) are shown for a simulation which has a diffusion coefficient which varies in space. In (a) and (b) t = 40, in Copyright and (d) t = 4000. The solutions are becoming almost uniform in space despite the different diffusion coefficients. The contour lines shown are for densities of 0.1 and 1. The simulation uses the reversing boundary conditions with the parameters the same as those in figure 4.2 except ht = 0.001, D = 0.01 on half the domain and D = 0.18 on the other half. The initial conditions are uniform for collagen and the fibroblasts are randomly perturbed from the uniform state with an amplitude of 0.3. shown earlier, the boundary conditions play an important role in determining how the solution will behave. For scar tissue formation this would means that the size of the wound and the properties of the adjacent tissue are important. Although much has been done in the study of alignment, there is still a need for further study. In our work an obvious next step is the further refinement of the fully two dimensional model to allow a comprehensive numerical investigation and the inclusion of a third spatial dimension. To make this modeling work more applicable to wound healing, many of the other processes which occur must be integrated into A spatiallyv arying model of alignment 21 Fig. 8.2. Isoclines for the collagen density are shown for a two-dimensional simulation similar to that in figure 8.1, at t 200. The isocline is plotted for a density of approximately 1. The parameters are the same as those in figure 8.1 except ht = 0.001, = 0.3, s = 3 and D = 0 for half the spatial domain and D = 0.3 for the other half. The y direction has periodic boundary conditions and is discretized using 30 grid points. The initial conditions are randomly perturbed about the uniform steady state with the maximum perturbation being 0.03. Two-dimensional simulations such as this are extremely time-consuming, with this case taking about 10 days on a SUN Ultra 1. the framework of the model. This means the inclusion of cell proliferation, forces generated in wound contraction, and incorporation of chemical signaling which plays a crucial role in wound healing. These are a few of the most important interactions which must be considered when extending our model towards a more realistic representation of dermal wound healing. 22 J.Copyright Dallon and J. A. Sherratt REFERENCES [1] V. H. Barocas and R. T. Tranquillo, An anisotropic biphasic theory of tissue-equivalent mechanics: the interplay among cell traction, fibrillar network deformation , fibril alignment and cell contact guidance, AMSE J. Biomech. Eng., 119 (1997), pp. 137 145. [2] L. Besseau and M. M. Giraud-Guille, Stabilization of fluid cholesteric phases of collagen to ordered gelated matrices, J. Mol. Biol., 251 (1995), pp. 197 202. [3] G. Civelekoglu and L. Edelstein-Keshet, Modelling the dynamics of f-actin in the cell, Bull. Math. Biol., 56 (1994), pp. 587 616. [4] R. A. F. Clark, Wound repair overview and general considerations, in The Molecular and Cellular Biologyof Wound Repair, R. A. F. Clark, ed., Plenum Press, 2 ed., 1996, pp. 3 50. [5] J. Cook, Waves of alignment in populations of interacting, oriented individuals, Forma, 10 (1995), pp. 171 203. [6] J.Copyright Dallon and J. A. Sherratt, A mathematical model for fibroblast and collagen orientation, Bull. Math. Biol., 60 (1998), pp. 101 129. [7] J.Copyright Dallon, J. A. Sherratt, and P. K. Maini, Mathematical modelling of extracellular matrix dynamics using discrete cells: Fiber orientation and tissue regeneration, J. Theor. Biol., 199 (1999), pp. 449-471. [8] P. J. Davies, F. J. Carter, D. G. Roxburgh, and A. Cuschieri, Mathematical modelling for keyhole surgery simulations: spleen capsule as an elastic membrane, J. Theor. Med., 1 (1999), pp. 247 262. [9] P. H. Ehrlich and T. M. Krummel, Regulation of wound healing from a connective tissue perspective, Wound Rep. Reg., 4 (1996), pp. 203 210. [10] C. A. Erickson, Analysis of the formation of parallel arrays by bhk cells in vitro, Exp. Cell Res., 115 (1978), pp. 303 315. [11] E. Geigant, K. Ladizhansky, and A. Mogilner, An integrodifferential model for orientational distributions of f-actin in cells, SIAM J. Appl. Math., 59 (1999), pp. 787 809. [12] S. Guido and R. T. Tranquillo, A methodology for the systematic and quantitative study of cell contact guidance in oriented collagen gels, J. Cell Sci., 105 (1993), pp. 317 331. [13] A. R. Mitchell and D. F. Griffiths, The Finite Difference Method in Partial Differential Equations, John Wiley& Sons, 1980. [14] A. Mogilner and L. Edelstein-Keshet, Selecting a common direction i. how orientational order can arise from simple contact responses between interacting cells, J. Math. Biol., 33 (1995), pp. 619 660. [15] , Spatio-angular order in populations of self-aligning objects: formation of oriented patches, Physica D, 89 (1996), pp. 346 367. [16] A. Mogilner, L. Edelstein-Keshet, and B. G. Ermentrout, Selecting a common direction ii. peak-like solutions representing total alignment of cell clusters, J. Math. Biol., 34 (1996), pp. 811 842. [17] I. O. L. Ng, E.Copyright S. Lai, M. T. Ng, and S. T. Fan, Tumor encapsulation in hepatocellular carcinoma, Cancer, 70 (1992), pp. 45 49. [18] A. Okubo, Dynamical aspects of animal grouping: swarms, schools flocks, and herds, Adv. Biophys., 22 (1986), pp. 1 94. [19] L. Olsen, P. K. Maini, J. A. Sherratt, and J.Copyright Dallon, Mathematical modelling of anisotropy in fibrous connective tissue, Math. Biosci., 158 (1999), pp. 145 170. [20] L. Olsen, J. A. Sherratt, P. K. Maini, and B. Marchant, Simple modelling of extracellular matrix alignment in dermal wound healing I. cell flux induced alignment, J. Theor. Med., 1 (1998), pp. 172 192. [21] S. P. Palecek, J.Copyright Loftus, M. H. Ginsberg, D. A. Lauffenburger, and A. F. Horwitz, Integrin-ligand binding properties govern cell migration speed through cell-substratum adhesiveness, Nature, 385 (1997), pp. 537 540. [22] T. D. Pollard and J. A. Cooper, Actin and actin-binding proteins. a critical evaluation of mechanisms and function, Ann. Rev. Biochem., 55 (1986), pp. 987 1035. [23] J. A. Sherratt, Travelling wave solutions of a mathematical model for tumour encapsulation, SIAM J. Appl. Math., (1999). in press. [24] J. A. Sherratt and J. Lewis, Stress-induced alignment of actin filaments and the mechanics of cytogel, Bull. Math. Biol., 55 (1993), pp. 637 654. [25] D. Stopak and A. K. Harris, Connective tissue morphogenesis by fibroblast traction, Dev. Biol., 90 (1982), pp. 383 398.