## 1 Introduction and biological background

Among the most crucial issues in development are the evolution of spatial patterns and the mechanisms which create them. Much progress has been made in furthering our understanding of the basic principles that patterning mechanisms must possess to be able to generate specific spatial patterns. Even so, we still do not know, with any certainty, definitive details of a single pattern formation mechanism which is involved in development. Morphogenetic model mechanisms can suggest possible scenarios as to how pattern is laid down, and sometimes when (as for example, in the experimentally confirmed case of stripe patterning on the alligator Alligator mississippiensis [1]) and how the embryonic form might be created. It is now widely accepted that dramatic progress can come about through a genuine interdisciplinary approach involving experimentalists and theoreticians. The books by Murray [2,3] discuss in detail many such successful case studies.

Embryogenesis proceeds through a series of sequential processes which generate specific patterns at each stage. The network pattern observed in vasculogenesis and angiogenesis, for example, are typical examples: how they are formed is an important and question particularly in view of the work of Folkman and his colleagues (see, for example, [4,5]) on the crucial role angiogenis plays in the growth of solid tumours.

Broadly speaking, there are two prevailing views of pattern generation that have influenced the thinking of embryologists. One is the long standing and well known Turing chemical pre-pattern approach [6] and the Murray–Oster continuum mechanochemical approach (see, for example, [7–10]). The former requires cells to respond to a concentration level such as a positional information approach [11] applications of which are given in [12]. Important related mechanisms involve cells and chemotaxis, the aggregation effect whereby cells move up concentration gradients of chemoattractants. Practical examples of both experiment and theory involving bacteria are described in [13–15]. General descriptions, a full discussion of pattern formation mechanisms and real biological examples are given in [2,3]. A good view of the current state of pattern formation modelling in several different areas is given in [16,17].

## 2 Mechanochemical mechanisms for generating pattern and form

The mechanochemical modelling framework is based on the observation [18] that “relatively simple cellular forces can give rise to complex changes in form”. This concept was clearly demonstrated by experiments in which fibroblast cells are placed on an elastic substratum. The cells spread and migrate, concurrently generating strong traction forces which deform the substratum creating tension patterns on the elastic sheet [19]. The latter experiments and the (even now) illusive nature of morphogens were the impetus for the initial development of the mechanochemical modelling approach [7–9]. The quantification of these types of cellular forces has been studied, for example, by Traqui and his colleagues [20] and Moore [21].

Mesenchymal cells, on which we focus here, are spread within a three-dimensional fibrous extracellular matrix (ECM) and migrate within the matrix: they can secrete (and degrade) matrix as well as exert traction forces on the ECM environment. That is, these cells are capable of significantly altering their environment. So, within this simple framework, the cells can create spatial heterogeneity and generate form simultaneously if the parameters are in appropriate experimentally reasonable ranges. Perhaps most crucial, however, is that with such a model mechanistic framework the variables involved and the consequences of their interaction can be observed and measured experimentally.

The power of mathematical modelling lies in its ability to make predictions and suggest experiments. Mechanochemical theory is particularly well-suited to connecting theory with experiment since the model structure is firmly based on empirically measurable quantities (e.g., cell density, matrix density, cell traction forces and known chemicals) and the consequence of varying such quantities can be clearly observed.

Here we give a brief introduction to the mechanochemical theory of morphogenesis and, by way of example, describe a recent successful application to vascular network formation where the theory and experiment are sufficiently close to be able to suggest some of the essential elements in the biological mechanism which create these specific biological patterns and final structure.

The basic mechanical model hinges on two key experimentally determined properties of mesenchymal cells in vitro: (i) cells migrate within a tissue substratum made up of fibrous ECM and other cells; (ii) cells can generate large traction forces. The models try to mimic the mechanical interaction between the motile cells and the viscoelastic substratum (ECM) within which they move.

Mesenchymal cells move by exerting forces on their surrounding, consisting of the viscoelastic fibrous ECM and the surface of other cells. They use their cellular protrusions, the filopodia or lamellapodia, which stretch out from the cell in all directions, to grip whatever is available and pull. As the cells move through the ECM they deform it by virtue of their traction forces. These deformations in the ECM induce anisotropy effects which in turn affect the cell motion. Analysis of models incorporating these various effects show that coordination of such effects result in spatially organized cell aggregations. The general modelling scenario underlying the models is shown schematically in Fig. 1.

The model is a continuum one, consisting of three equations governing (i) the cell population density, (ii) the mechanical balance of the forces between the cells and the ECM and (iii) the conservation law for the ECM. Here we give only the equations in words to indicate the principles: brief mathematical details are given in the appendix with full discussions in [2,3].

### 2.1 Cell density conservation

Many factors, with varying levels of importance, affect the movement of embryonic cells within a tissue substratum which is made up of fibrous ECM. For example, some of the more obvious are:

- – chemotaxis, the process by which cell motion is directed by a chemical gradient,
- – contact guidance, when the substratum defines a preferred direction of cell motion,
- – contact inhibition, inhibiting motion and reducing individual cell traction forces as a consequence of high cell densities,
- – convection, which transports cells and ECM due to local deformation in the matrix.

The forces exerted on the ECM by other cells induce passive movement of cells on the ECM. Diffusion is the usual random dispersal of cells with general motion down a cell density gradient: it can include biased or tissue-guided diffusion. Haptotaxis is the process by which ECM gradients affect cell movement. Cell tractions generate gradients of matrix density, which are in turn associated with the density of adhesive cites for cell lamellapodia to bind to. Since cells can get a stronger grip on denser matrix (up to a certain density) cells tend to move up an adhesive/matrix gradient. Galvanotaxis is movement due to the field generated by electric potentials, which are known to exist in embryos. These effects are all well documented experimentally. The model field equations encapsulate the key features which affect cell movement within its extracellular environment. For illustration we shall not include all of the effects just mentioned but it will be clear how they can be incorporated: see [2,3] for a full discussion and the mathematical representation and quantification of the terms.

Within a fixed control volume, cell density is conserved according to equation (1):

$$\mathrm{rate}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{change}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{cellular}\phantom{\rule{0.277778em}{0ex}}\mathrm{density}=\mathrm{convection}+\mathrm{diffusion}+\mathrm{haptotaxis}+\mathrm{mitosis}-\mathrm{cell}\phantom{\rule{0.277778em}{0ex}}\mathrm{death}$$ | (1) |

Contact guidance is involved in the formation of cell-matrix networks discussed below. Such guidance can be modelled in different ways. One possibility is to make the diffusion directed diffusion. In this case, the diffusion coefficient is a function of the elastic strain of the ECM (see, for example, [3]). This dependence implies that the strain of the matrix can grossly affect the direction of cell movement as we shall see in the section on network formation.

The matrix mechanical behaviour is dominated by the fibrous collagen component of the ECM. Local disturbances in the matrix density can induce convection of the matrix. Also, as (mesenchymal) cells move through the matrix environment, they can secrete and/or degrade the ECM: this occurs, for example, in wound healing.

### 2.2 Conservation of matrix density

The conservation equation for the ECM density is then given by equation (2):

$$\mathrm{rate}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{change}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{matrix}\phantom{\rule{0.277778em}{0ex}}\mathrm{density}=\mathrm{convection}+\mathrm{secretion}-\mathrm{degradation}$$ | (2) |

In many applications, the secretion term is often negligible on the time scale of the pattern formation, as is the case in the experiments associated with the matrix network related to vasculogenesis discussed below.

### 2.3 Mechanical interaction of cell and matrix

Since inertial terms are small in development, the governing equation for the mechanical interactions between the cells and matrix is simply that the sum of the various forces are in equilibrium within the tissue. We thus have (equation (3)):

$$\mathrm{forces}\phantom{\rule{0.277778em}{0ex}}\mathrm{on}\phantom{\rule{0.277778em}{0ex}}\mathrm{the}\phantom{\rule{0.277778em}{0ex}}\mathrm{body}=\mathrm{cell}\phantom{\rule{0.277778em}{0ex}}\mathrm{traction}\phantom{\rule{0.277778em}{0ex}}\mathrm{forces}+\mathrm{viscoelastic}\phantom{\rule{0.277778em}{0ex}}\mathrm{restoringforces}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{the}\phantom{\rule{0.277778em}{0ex}}\mathrm{matrix}=0$$ | (3) |

### 2.4 Viscoelastic matrix material

The stress–strain relation in such a complex tissue matrix is almost certainly nonlinear and plastic. There is still no experimental evidence which would allow us to postulate the actual form of the relation. Many patterns result from relatively small deformations and cell density differences (one example occurs in the formation of papillae, the dermal condensations of cells which presage feather and scale formation). We make the assumption that for small strains, cell-matrix material forces can be modelled by a linear viscoelastic response [22]. The matrix stress then consists of a viscous and an elastic component:

$${\mathrm{stress}}_{\mathrm{matrix}}={\mathrm{stress}}_{\mathrm{viscous}}+{\mathrm{stress}}_{\mathrm{elastic}}$$ |

### 2.5 External body forces

There are often external forces acting on the tissue body we are attempting to model. These forces may be due to attachment of the matrix to a substrate or in the case of network formation we discuss below, to the base of the experimental vessel. The body forces holding the tissue in various places can reasonably be modelled as proportional to the density of matrix and the displacement of the matrix. Alternatively in the matrix network situation discussed later it is more appropriate to relate the attachment force to the viscous drag: we come back to this below.

Although at the time of the formulation of this modelling approach [7–9] the mechanical properties of ECM and cell-matrix interactions had not been quantified, recent experimental advances have focused on these mechanical aspects of tissue interactions. As a result, many of the parameters described in the model formulation can be estimated. As well as the work mentioned above [20,21] the viscoelastic properties of mesodermal cells during early development have been quantified [22]. In the cell-matrix experiments described below parameter estimations have also been obtained.

In many developmental situations (wound healing is one) chemicals can play a significant role in the pattern formation process and therefore should be included in the model. Although chemical dynamics can be incorporated in principle relatively simply, but only if we have some idea of how they interact with the cells and matrix. At this stage, these interactions are generally still uncertain.

The mathematical investigation of these model equations involves linear analysis about the steady states to obtain what is called the dispersion relation which predicts the pattern formation potential. The analysis (see, for example, [2 or 3] for a full discussion) results in a remarkably wide variety of dispersion relations which indicate a rich spectrum of nonlinear spatial patterns. The nonlinear spatial patterns almost always have to be obtained by numerical simulation of the equations. As expected the type of patterns possible depend on which effects are included in the model mechanism. Another aspect of these mechanical models is their robustness to parameter variation. It is essential to know, at least approximately, the parameter domain for specific patterns because this is one of the optimal ways to assess model robustness. In the remainder of this article we discuss an application of mechanochemical theory to vascular network formation.

## 3 In vitro vascular network formation

Mechanical and fluid mechanical forces play essential roles in the overall development of vasculature. As early as 1893 [23] the importance of (fluid) mechanical factors in the growth of blood vessels during development was pointed out. In a review [24] of this and other works on mechanical forces on angiogenesis these early observations of how vascular sprouting in the growing embryo might occur as a result of a combination of velocity of blood and pressure are described. Another review on some of the roles of mechanical forces in vascular development and remodeling is given in [25].

The development of in vitro angiogenesis systems is a controlled means of studying vessel formation. Such systems have shown the important mechanical role that the ECM plays in angiogenesis. It provides, among other things, a scaffold necessary for cell migration and morphogenesis (for example, [26,27]). As we described, cells can not only produce or degrade their ECM but also alter its structure by applying mechanical forces. Through matrix production and degradation, cells can influence the mechanical properties of their ECM and, through their mechanical forces, they can reorganise its fibrous components into matrix lines that the cells use as migratory pathways and movement cues. Such a mechanical scenario in which cell-ECM interactions can be orchestrated to form complex spatial patterns in development make it a very good candidate for the application of the Murray–Oster mechanical theory. The ability of the cell-matrix interaction to effect an alignment of the matrix thus influencing cell movement is quite common in development: some discussion of these is given in [26–29] but the key reference is the 1998 book [30] on vascular morphogenesis.

Of particular relevance to this section is in vitro work on vascular network formation [26,27]. This work helps to elucidate the mechanism of vascular network formation in vitro [26]. It is shown that networks are not specific to endothelial cells since a variety of traction exerting cells, not necessarily endothelial cells, can form networks when cultured on gelled basement membrane (Matrigel). Moreover, cells cultured on different substrates (Matrigel, type I collagen gels), also form networks provided the matrix is malleable enough. These studies suggest that the mechanism lacks cellular/matrical specificity and that a physical/mechanical mechanism could possibly explain the generality of the network forming process.

Here we describe a basic mathematical mechanical model mechanism with the aim of showing that a purely mechanical mechanism could be responsible for the observed patterns and how they are actually formed in development. Such a mechanical model does not specify the type of cells and matrix involved but rather only considers the possible mechanical interactions between the various components. We test the model against the patterns observed by particular experiments and we further use the model to quantitatively describe the interaction between the various mechanical properties and the network and their pattern forming abilities.

Experiments [26,27] show that bovine aortic endothelial cells (BAEC), cells of the murine Leydig cell line TM3, adult human dermal fibroblasts and human smooth muscle cells when cultured on Matrigel are observed to form networks. The process for all cell lines is similar: cells adhere onto the matrix and start pulling on it. The pulling results in movement of the matrix and the cells that have adhered to it, eventually forming cell aggregates with significant amounts of matrix accumulated underneath the aggregates. As a result cellular traction tension lines appear around the clusters. In time, the matrix components appear to form fibrous lines between neighboring clusters and along the tension lines. Once the lines form, the cells become actively motile and migrate onto the matrix pathways, thus forming cellular cords. Eventually the culture plane is tessellated with polygons with sides defined by the cellular cords [26]: see also Fig. 2a.

Collagen (an ingredient of tissue matrix) content influences the formation of networks, presumably by altering matrix stiffness and thereby the effect of the cell-exerted traction onto the matrix. Thickness of the matrix layer also influences the size of the networks formed: thin layers of matrix resulted in small or no network formation. On a ramp of increasing matrix thickness, larger cellular networks form on the thicker areas of matrix: refer also to Fig. 3c.

These results suggest that mechanical interactions are of primary importance for the development of pattern. Mathematical descriptions of the interactions, which we describe below, confirm the crucial role that mechanical forces play in the network pattern formation and provide a tool for assessing which mechanical properties of the ECM control the patterning process as well as how the parameters may modulate pattern size. It also provides various experimental scenarios to highlight the effect of experimentally variable parameters, such as gel thickness, density of cells and so on.

Since in the experiments [26,27] there is neither matrix production nor degradation we require only two equations, namely cell density conservation and a force balance equation. The effect of the matrix density is incorporated in the matrix thickness on the dish. As the cells pull the matrix, it moves across the dish. Experiments indicate that some matrix fibrils remain attached to the dish, while the rest are dragged across the lower parts of the matrix. The net effect on the surface of the matrix of the attachment of the fibrils on the dish, is a resistance, which in a planar situation we treat as a resistive viscous drag.

Because the thickness of the experimental gel on which the cells are embedded is very thin, and in particular, relative to the cell-matrix spatial domain, a two-dimensional model is sufficient at this stage.

### 3.1 Mathematical model

The two-dimensional mathematical model quantifies the basic experimental mechanical scenario described above. We do not include cell proliferation: on the Matrigel cultures the first cell aggregates and lines of tension appear after 4 h and the networks are complete within about 24 h. The time between two subsequent cell mitoses is about 17 hours for endothelial cells so we assume that no significant changes in the cell populations occur which could influence the pattern forming process.

Local changes in the cell density are a combination of mainly two movements, the convective flux and an anisotropic strain-dependent random motion tensor. We thus consider that cell movement can be approximated by a random walk that is biased along areas of matrix alignment. The resulting conservation equation for the cell density in words from equation (1) is then (equation (4)):

$$\mathrm{rate}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{change}\phantom{\rule{0.277778em}{0ex}}\mathrm{of}\phantom{\rule{0.277778em}{0ex}}\mathrm{cell}\phantom{\rule{0.277778em}{0ex}}\mathrm{density}=\mathrm{convection}+\mathrm{strain}-\mathrm{dependent}\phantom{\rule{0.277778em}{0ex}}\mathrm{diffusion}$$ | (4) |

In the initial analysis [28], we did not consider strain-dependent active movement so that we could test whether network formation is possible under just cell-exerted traction.

#### 3.1.1 Forces within the extracellular matrix

The width of the matrix sample in the dish is much larger than the thickness of the matrix so we approximate the matrix as a two-dimensional material. Moreover, we assume that the traction exerted by the cells remains confined in the plane parallel to the dish (or, approximately, the matrix) surface. In other words, we consider the material to be under a plane-stress assumption. The movement of the matrix is resisted by attachment of the matrix to the dish. In vivo planar angiogenesis is, by comparison, influenced by the potential attachment and continuation of the fibrous and cellular components of the adjacent tissue and is more accurately described by the three-dimensional model equivalent.

The key forces are of course those generated by the cells. The others are generated as a direct consequence. So, forces that are present in the tissue are: (i) the cell-exerted traction; (ii) the resistance due to the matrix-dish contact; (iii) the viscoelastic forces of the matrix material, which are resisting the deformation. Again since the time required for generation of the pattern is relatively long and the size of the pattern is small in absolute terms, inertia effects are negligible and the forces at any given point are again considered to be in equilibrium.

In response to a force, the matrix is displaced slowly due to its viscous properties. If strains developed in collagen gels are small (less than 10%) then the gels have been found to respond linearly to applied stresses [32,33]. We therefore again describe the matrix response as that of a linear viscoelastic body.

The equation reflecting the balance of these forces (F) is then (cf. equation (3))

$${F}_{\mathrm{cells}}\phantom{\rule{3.30002pt}{0ex}}(\mathrm{cell}-\mathrm{exerted}\phantom{\rule{0.277778em}{0ex}}\mathrm{traction})+{\mathrm{F}}_{\mathrm{anchoring}}\phantom{\rule{3.30002pt}{0ex}}\left(\mathrm{attachment}\phantom{\rule{0.277778em}{0ex}}\mathrm{on}\phantom{\rule{0.277778em}{0ex}}\mathrm{dish}\right)+{\mathrm{F}}_{\mathrm{matrix}}\phantom{\rule{3.30002pt}{0ex}}\left(\mathrm{viscoelastic}\phantom{\rule{0.277778em}{0ex}}\mathrm{restoring}\phantom{\rule{0.277778em}{0ex}}\mathrm{forces}\right)=0$$ | (5) |

_{cells}, the forces in the matrix due to the cells and F

_{matrix}the viscoelastic force in the matrix material. F

_{anchoring}is an external (body) force resisting the matrix displacement. The cell traction F

_{cells}depends on the local cell density initially increasing with the cell density eventually decreasing with large enough cell densities.

The viscoelastic stresses we consider to be made up of elastic and viscous parts as in the last section. The difference here is associated with the drag of the matrix across the dish. As the cells pull the matrix, it moves across the dish. Experiments indicate that some matrix fibrils remain attached to the dish, while the rest are dragged across the lower parts of the matrix. The net effect on the matrix surface due to the fibril attachments to the dish, is a resistant force. In this two-dimensional model we treat such resistance as a viscous drag.

Since the vertical stress is zero we can derive a relation between the strain components in the three directions [3,28].

The detailed quantification of the two equations for cell conservation and balance of forces together constitute the model mechanism. The equations were solved with the appropriate boundary conditions and initial conditions in cell density and matrix distribution in keeping with the experimental set-up [26,27] (see the Appendix A and [3] for a full discussion). Where these parameters specifically appear in the model is shown in the Appendix A.

### 3.2 Estimation of the parameters

Since there is a close relation between this application of the mechanical theory and in vitro experiments we were able to derive estimates for various parameters describing cell behavior and matrix properties. This is one of the most important aspects of the application of the theory. These parameters are defined in the Appendix A where it is shown where they appear in the model equations.

The stiffness of the matrix is represented by the elastic modulus of the material (E, $\mathrm{dyn}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{cm}}^{-2}$) and is determined primarily by the type, amount and organization of its fibrous components. For example, the anterior porcine vitreous body and the posterior bovine vitreous body were found to have similar collagen contents and gave elastic moduli of $26.93\phantom{\rule{3.33333pt}{0ex}}\mathrm{dyn}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{cm}}^{-2}$ and $8.01\phantom{\rule{3.33333pt}{0ex}}\mathrm{dyn}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{cm}}^{-2}$, respectively [31,34]. Since the collagen gels used in the vascular network-forming experiments have collagen content comparable to that of the vitreous bodies mentioned, we assume that these gels' elastic moduli are within a similar range.

The shear viscosity of $2.1\phantom{\rule{3.33333pt}{0ex}}\mathrm{mg}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{ml}}^{-1}$ collagen I gel has been measured using creep tests to be 7.4×10^{6} Pe [33] while that of the posterior bovine vitreous body was estimated to be 2.5×10^{2} Pe [32].

Traction per cell (τ dyn/cell) we consider to be comparable to that of human umbilical vein endothelial cells which is approximately $6.1\times {10}^{4}\phantom{\rule{3.33333pt}{0ex}}\mathrm{dyn}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{cm}}^{-2}$ [34]. From these results we estimate the traction force per cell to be in the region of τ≈0.15– 0.27 dyn/cell.

The Poisson ratio (ν, $\mathrm{cm}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{cm}}^{-1}$) can be theoretically measured by extending a material by a certain amount and measured by its compression on the perpendicular direction. The ratio of relative extension to compression on the perpendicular direction gives the magnitude of the parameter. In practice, such measurements for soft materials such as collagen gels can be very difficult, not least because the value of ν changes in time. The Poisson ratio for the fibre network of such gels is estimated at ν≈0.2 [35].

Cell motility for endothelial cells has been calculated using migration assays. The values reported for endothelial (Human Microvessel Endothelial Cells) on gelatin were $9.5\pm 1.2\times {10}^{-9}\phantom{\rule{3.33333pt}{0ex}}{\mathrm{cm}}^{2}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{s}}^{-1}$ in the presence of endothelial cell growth factor and $2.6\pm 0.6\times {10}^{-9}\phantom{\rule{3.33333pt}{0ex}}{\mathrm{cm}}^{2}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{s}}^{-1}$ in its absence and $19.3\pm 4.22\times {10}^{-9}\phantom{\rule{3.33333pt}{0ex}}{\mathrm{cm}}^{2}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{s}}^{-1}$ on fibronectin [36].

In our investigation, we were concerned with the initial stages of network pattern emergence where the strains are still small and we can assume that the parameters stay approximately constant.

## 4 Results

The model mechanism system of equations was analyzed analytically and numerically [3,28,29]. In Fig. 2 we show a selection of the results of the analysis. Analysis showed that small variations in cell density could trigger matrix displacements if the following relation between the model (and real biological) parameters obtains:

$$\frac{\tau}{E}>\frac{1}{1-{\nu}^{2}}\frac{1}{{c}_{0}}$$ | (6) |

_{0}=4×10

^{4}–$2\times {10}^{5}\phantom{\rule{3.33333pt}{0ex}}\mathrm{cells}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{cm}}^{-2}$ [26] is the initial cell seeding density.

The relation (6) suggests a basic, and intuitively reasonable, rule that could make patterns possible: patterns are more likely to form if the cell traction τ is sufficiently large, or if the matrix stiffness E is sufficiently low or if the initial seeding density is relatively large. The results are in general agreement with experimental findings. For example, when cell traction was inhibited or matrix stiffness greatly increased, no patterns formed [26]. Also, increases in the cell plating densities facilitates the formation of pattern. We can see why, since if few cells are seeded, they have difficulty in deforming the matrix and starting the patterning process.

This inequality predicts when pattern could form, but does not predict the type of pattern formed. The pattern described by the equations can be determined by solving the governing equations numerically. Numerical simulations [28] show that network formation is possible through a purely mechanical process as in Fig. 2 (see also Fig. 3).

The numerical simulations clearly showed network formation and subsequent remodelling: some polygons grew larger and others closed in, in a sphincter-like manner. There was a direct correspondence between areas of large cell density and areas of high matrix density, just as we would expect (see the experimental image in Fig. 2a).

The choice of parameters influenced the formation of pattern and rate at which the network pattern formed. If traction was too low, the plating density too low, or the stiffness too high, no patterns formed. If, on the other hand, stiffness was too low or traction too high, cells and the matrix formed clusters but no matrix cords were observed.

Viscosity of the matrix, as to be expected intuitively, influenced the time it took for the pattern to appear, provided traction and stiffness were within the appropriate range, but did not determine whether pattern would form or not.

The matrix density played a significant role in the size of the network formed. Since the matrix density in these experiments is directly related to the thickness of the membrane we simulated a situation in which there was a gradient in matrix density. We predicted that a network with a gradient in network size would be formed. Fig. 3 shows the results [28] of such a numerical simulation of the model equations. Subsequent experiments confirmed this prediction.

We then investigated how essential the directed diffusion was by solving the equations without any diffusion. One of the initially surprising and important findings was that random motility of cells was not necessary for the formation of pattern. Networks would still form, provided the seeding density was sufficiently large.

### 4.1 Results from the mechanical theory

Mathematical analysis and simulation of the model mechanism give insight as to which elements are crucial to the development of pattern, help assess the relative importance of such elements, and predict how changes in these factors affect pattern growth and shape. In this way, the model is useful in assessing the biological conditions under which a pattern may or may not form.

We have presented a simple, purely mechanical model for network pattern formation in vitro, in which cellular traction forces play a crucial role in forming the planar pattern of cellular aggregates, cords and matrix deformations observed in vitro. In spite of the simplifying assumptions in the model the results show that in the absence of traction, as in the in vitro system, no pattern forms, supporting our hypothesis that the pattern forms as a result of cellular traction forces. The evolution of the pattern from the model proceeds in the same manner as the in vitro pattern, with irregular polygons both increasing in size and decreasing in number, with the smallest polygons pinching off and vanishing.

As in the in vitro system, pattern forms in our system only when the ratio of traction forces and seeding cell density to gel stiffness is above a certain threshold. Our model results confirm that on thicker gel larger polygons will form. It also indicates that isotropic cellular traction is sufficient for pattern formation. We further found that if cellular traction is sufficiently high or the plating densities large, then traction forces alone can lead to formation of networks. Biased cellular migration may be a component of pattern formation in vivo and in vitro, but we have shown it is not a necessary feature.

We also found that matrix thickness is an important factor influencing the pattern. Cells spread on a matrix whose thickness increased from one end of the dish to the other, formed polygon networks whose size increased with increasing thickness of the gel. The numerical simulation results are in agreement with the experimental results [26].

## 5 Discussion and some general concluding remarks

The mechanochemical theory of pattern formation has been proven to be important in the study of a spectrum of biological processes such as wound healing, angiogenesis, and limb development, to name just a few. One of the significant features of the applicability of this theoretical framework to a wide variety of problems is that the definable parameters are in principle all measurable experimentally. Although a great deal of experimental work has been devoted to the molecular biology of development, few studies have focused on mechanical aspects: they are technically difficult.

Although certain morphogenetic events can be described by chemical and/or mechanical mechanisms, the vast majority of pattern formation events are the result of the interaction of chemical and mechanical processes. For this reason, mechanochemical theory is a broadly applicable, powerful tool in defining mechanisms for numerous developmental events.

Why use mathematics to study something as intrinsically complicated and ill-understood as development, angiogenesis, wound healing and so on? We suggest that mathematics, rather mathematical modelling, must be used if we ever hope to genuinely and realistically convert an understanding of the underlying mechanisms into a predictive science. Mathematics is required to bridge the gap between the level on which most of our knowledge is accumulating (cellular and below) and the macroscopic level of the patterns we see. In wound healing, scar formation itself is a major concern of the surgeon and the patient. A mathematical approach lets us explore the logic of pattern formation. Even if the mechanisms were well understood – and they certainly are not at this stage – mathematics would be required to explore the consequences of manipulating the various parameters associated with any particular scenario. In the case of wound healing, and it will be increasingly so in angiogenesis with the related cancer component, the number of options that are fast becoming available to wound and cancer managers will become overwhelming unless we can find a way to simulate particular treatment protocols before applying them in practice. The latter has been already of use in understanding the efficacy of various treatment scenarios with brain tumours (glioblastomas) (see [3] for a full discussion and new two-step treatment regimes for skin cancer [37].

There is no doubt that we are a long way from being able to reliably simulate actual developmental scenarios. Not only is the active cellular control of the key processes poorly understood, in the case of wounds, for example, they are difficult to reproduce, with similar wounds on different parts of the same body healing at different rates and with different results. Despite these limitations, we argue that exploring the logic of biological processes is worthwhile even in our present state of knowledge. It allows one to take a hypothetical mechanism and examine its consequences in the form of a mathematical model, making predictions and suggesting experiments that would verify or invalidate the model; even the latter casts light on the biology. Indeed, the very process of constructing a mathematical model can be useful in its own right. Not only must one commit to a particular mechanism, but one is also forced to consider what is truly essential to the process, the key players (variables) and the key processes by which they evolve. We are thus involved in constructing frameworks on which we can hang our understanding. The equations, the mathematical analysis and the numerical simulations that follow serve to reveal quantitatively as well as qualitatively the consequences of that logical structure.

Theoretical studies have served to highlight where our knowledge is deficient and to suggest directions in which fruitful experimentation might lead us. A crucial aspect of this research is the interdisciplinary content. A critical test of all of these theoretical constructs is in their impact on the experimental community. We believe that genuine collaboration between experimentalists and applied mathematicians will lead us more rapidly towards a fuller understanding, if not a complete one, of the biological processes involved in pattern formation and a vast array of medical problems.

## Appendix A

Here we briefly give the mathematical description of the mechanical theory as it applies to vascular network theory.

The two-dimensional mathematical model quantifies the basic experimental mechanical scenario described above. We do not include cell proliferation: on the Matrigel cultures the first cell aggregates and lines of tension appear after 4 h and the networks are complete within about 24 h. The time between two subsequent cell mitoses however is about 17 h for endothelial cells so we assume that no significant changes in the cell populations occur which could influence the pattern forming process.

Local changes in the cell density are then a combination of mainly two movements, the convective flux and an anisotropic strain-dependent random motion tensor. We thus consider that cell movement can be approximated by a random walk that is biased along areas of matrix alignment. The conservation equation for the cell density, denoted by c ($\mathrm{cells}\phantom{\rule{1.69998pt}{0ex}}{\mathrm{mm}}^{-2}$), is

$$\frac{\partial \mathrm{c}}{\partial \mathrm{t}}=-\nabla \xb7\left(c\frac{\partial \mathrm{u}}{\partial \mathrm{t}}\right)+\nabla \xb7\nabla \xb7\left(\mathrm{D}\left(\u03f5\right)\mathrm{c}\right)$$ |

$$\mathrm{D}\left(\u03f5\right)={\mathrm{D}}_{0}\left(\begin{array}{cc}\multicolumn{1}{c}{1+\frac{{\u03f5}_{\mathrm{xx}}-{\u03f5}_{\mathrm{yy}}}{2}}& \multicolumn{1}{c}{\frac{{\u03f5}_{\mathrm{xy}}+{\u03f5}_{\mathrm{yx}}}{2}}\\ \multicolumn{1}{c}{\frac{{\u03f5}_{\mathrm{xy}}+{\u03f5}_{\mathrm{yx}}}{2}}& \multicolumn{1}{c}{1-\frac{{\u03f5}_{\mathrm{xx}}-{\u03f5}_{\mathrm{yy}}}{2}}\end{array}\right)$$ |

_{xx}, ε

_{yy}, ε

_{xy}, ε

_{yx}represent the two-dimensional components of the strain tensor, ε, and D

_{0}is the motility coefficient when no strain is present.

The key forces are those generated by the cells. The others are generated as a direct consequence. So, forces that are present in the tissue are: (i) the cell-exerted traction; (ii) the resistance due to the matrix-dish contact; (iii) the viscoelastic forces of the matrix material which are resisting the deformation. Since inertia effects are negligible the forces at any given point are considered to be in equilibrium.

In response to a force, the matrix is displaced slowly due to its viscous properties. We describe the matrix response as that of a linear viscoelastic body. The equation reflecting the balance of these forces is then

$${F}_{\mathrm{cells}}\phantom{\rule{3.30002pt}{0ex}}(\mathrm{cell}-\mathrm{exerted}\phantom{\rule{0.277778em}{0ex}}\mathrm{traction})+{\mathrm{F}}_{\mathrm{anchoring}}\phantom{\rule{3.30002pt}{0ex}}\left(\mathrm{attachment}\phantom{\rule{0.277778em}{0ex}}\mathrm{on}\phantom{\rule{0.277778em}{0ex}}\mathrm{dish}\right)+{\mathrm{F}}_{\mathrm{matrix}}\phantom{\rule{3.30002pt}{0ex}}\left(\mathrm{viscoelastic}\phantom{\rule{0.277778em}{0ex}}\mathrm{restoring}\phantom{\rule{0.277778em}{0ex}}\mathrm{forces}\right)=0$$ |

_{cells}=∇·σ

_{cells}, with σ

_{cells}the stress tensor in the matrix due to the cells and F

_{matrix}=∇·σ

_{matrix}with σ

_{matrix}the viscoelastic stress tensor in the matrix material. F

_{anchoring}is an external (body) force resisting the matrix displacement.

The cell traction F_{cells} depends on the local cell density as described earlier and it decreases with large enough cell densities. We take the cell-exerted stress on the matrix due to the cells as

$${\sigma}_{\mathrm{cells}}=\frac{\tau \phantom{\rule{1.69998pt}{0ex}}c}{1+\lambda \phantom{\rule{1.69998pt}{0ex}}{c}^{2}}\phantom{\rule{1.69998pt}{0ex}}\mathrm{I}$$ |

The viscoelastic stresses we consider to be made up of elastic and viscous parts as described generally in the section on the mechanical theory. The stress–strain relation in such a complex tissue matrix is almost certainly nonlinear and plastic. There is still no experimental evidence which would allow us to postulate the actual form of the relation. Since many patterns result from relatively small deformations and cell density differences we make the assumption that for small strains, cell-matrix material forces can be modelled by a linear viscoelastic response [22]. The matrix stress tensor then consists of a viscous and an elastic component:

$${\sigma}_{\mathrm{matrix}}={\sigma}_{\mathrm{viscous}}+{\sigma}_{\mathrm{elastic}}$$ |

$${\sigma}_{\mathrm{viscous}}={\mu}_{1}\phantom{\rule{1.69998pt}{0ex}}{\u03f5}_{t}+{\mu}_{2}\phantom{\rule{1.69998pt}{0ex}}{\theta}_{t}\phantom{\rule{1.69998pt}{0ex}}\mathrm{I}$$ |

_{1}and μ

_{2}represent the bulk and shear viscosities of the matrix, respectively.

The elastic component of the stress is given by

$${\sigma}_{\mathrm{elastic}}=\frac{E}{1+\nu}\left(\u03f5+\frac{\nu}{1-2\nu}\phantom{\rule{1.69998pt}{0ex}}\theta \phantom{\rule{1.69998pt}{0ex}}\mathrm{I}\right)$$ |

In this application the difference is associated with the drag of the matrix across the dish. As the cells pull the matrix, it moves across the dish. Experiments indicate that some matrix fibrils remain attached to the dish, while the rest are dragged across the lower parts of the matrix. The net effect on the matrix surface due to the fibril attachments to the dish, is a resistant force. In this two-dimensional model we treat such resistance as a viscous drag:

$${F}_{\mathrm{anchoring}}=-\mathrm{s}\phantom{\rule{1.69998pt}{0ex}}\frac{{\mathrm{u}}_{t}}{h}$$ |

### Matrix thickness

Since the vertical stress σ_{zz}=0, we get a relation between the strain components in the three directions, ε_{xx}, ε_{yy}, ε_{zz}. From Hooke's law in three dimensions we have

$${\sigma}_{\mathrm{zz}}=\frac{E}{1+\nu}\left[(1-\nu ){\u03f5}_{\mathrm{zz}}+\nu ({\u03f5}_{\mathrm{xx}}+{\u03f5}_{\mathrm{yy}})\right]=0.$$ |

$${\u03f5}_{\mathrm{zz}}=\frac{-\nu}{1-\nu}({\u03f5}_{\mathrm{xx}}+{\u03f5}_{\mathrm{yy}})$$ | (A.1) |

When there is a strain ε_{zz} in the z-direction, the thickness is calculated using h(x,y,t)=h_{0}(x,y) (1+ε_{zz}) from which we derive the relation giving the thickness h as a function of space and time:

$$\mathrm{h}(\mathrm{x},\mathrm{y},\mathrm{t})={\mathrm{h}}_{0}(\mathrm{x},\mathrm{y})(1+{\u03f5}_{\mathrm{zz}})={\mathrm{h}}_{0}\left[1-\frac{\nu}{1-\nu}({\u03f5}_{\mathrm{xx}}+{\u03f5}_{\mathrm{yy}})\right]\Rightarrow \mathrm{h}(\mathrm{x},\mathrm{y},\mathrm{t})={\mathrm{h}}_{0}\left(1-\frac{\nu}{1-\nu}\theta \right)$$ |

_{xx}+ε

_{yy}.

### Boundary conditions for the model equations

Again, being guided by the experiments, the matrix barely moves along the edge of the dish. We incorporate this into the model by assuming zero displacement as boundary conditions. Hence, for a square dish (square rather than circular simply for ease of numerical programming and reduction of computer time) of dimensions [ab]×[ab] the condition becomes:

$$\mathrm{u}(\mathrm{x},\mathrm{y}=\mathrm{a},\mathrm{t})=\mathrm{u}(\mathrm{x},\mathrm{y}=\mathrm{b},\mathrm{t})=\mathrm{u}(\mathrm{x}=\mathrm{a},\mathrm{y},\mathrm{t})=\mathrm{u}(\mathrm{x}=\mathrm{b},\mathrm{y},\mathrm{t})=0$$ | (A.2) |

$${J}_{\mathrm{cells}}=\frac{\partial \mathrm{u}}{\partial \mathrm{t}}\mathrm{c}-\nabla \xb7\left(\mathrm{D}\left(\u03f5\right)\mathrm{c}\right)=0$$ |

_{cells}is the flux of cells.

The equations for cell conservation and balance of forces together constitute the model mechanism. The equations were solved with the above boundary conditions and initial conditions in cell density and matrix distribution in keeping with the experimental set-up and the results shown in Figs. 2 and 3 above.