A model for chondrogenic condensations in the developing limb: the role of extracellular matrix and cell tractions

In a previous paper (Oster, Murray & Harris, 1983; hereafter referred to as OMH) we demonstrated how cell tractions could create spatial patterns of cell aggregation. There we proposed models for the morphogenesis of feather germ patterns and the patterns of condensation of chondroblasts that presage bone formation. These models were built on earlier experimental work by Harris and his coworkers on the role of cell traction in the aggregation process of mesenchymal cells (Harris, Stopak & Wild, 1981; Harris, Stopak & Warner, 1984).

The cell traction model as presented in OMH depended on cell motion; that is, the aggregation patterns formed by the migration of motile cells which were guided by haptotactic cues set up by the cell tractions. We also briefly mentioned the possibility that condensation patterns could occur without any cell motion other than passive ‘convection’ (i.e. dragging of cells by the traction of other cells).

Furthermore, the extracellular matrix (ECM) played a passive mechanical role in OMH: it served only as an elastic substrate which could deform under the tractions of cells and transmit stresses from one point to another.

In this paper we show how the ECM may be an active participant in the formation of spatial patterns (Toole, 1972). We present a physical mechanism that can generate spatial patterns of cell aggregation in the absence of active cell motility. Of course, if cell migration is present, the pattern formation possibilities are enhanced; however, our point here is that patterning can be generated without cell locomotion.

Intuitively, the mechanism works like this. The ECM is composed of a poly electrolyte gel which is osmotically swollen. If the cells commence secreting an enzyme which digests the osmotically active component of the ECM, the gel will deswell, and bring the cells closer together. When the cells are brought close enough together by the collapsing ECM, intercellular contacts can form (e.g. via filopodia) and cell tractions between cells become more effective. By analysing the physics of this process it turns out that the aggregations of cells produced by the gel deswelling of the ECM, coupled with the cell tractions, exhibit the same spatial regularities as described in OMH for aggregations produced by cell traction and migration. We note that the model we present here is based on entirely different physical assumptions than hitherto proposed models (e.g. Ede & Law, 1969; Newman & Frisch, 1979; Wolpert, 1971).

Our model of the developing limb bud is illustrated in Fig. 1. Cells proliferate at the distal tip and emerge from the progress zone secreting extracellular matrix. It is the hydration of the ECM which keeps the limb bud ‘inflated’, and maintains the spacing between cells. There is also a ‘sleeve’ of high hyaluronate concen-tration encasing the limb bud (Feinberg & Beebe, 1983). This tends to keep the limb bud osmotically inflated independently of the internal cell-matrix consti-tution. This sleeve will provide a ‘boundary condition’ for the model. The formation of condensation patterns will emerge from the interaction of forces generated by the ECM and by the cells. In order to understand this we discuss each force separately.

The extracellular matrix is a complex of crosslinked polymer molecules whose composition and properties change during development (Stern, 1984; Hay, 1981; Alberts et al. 1983; Trinkaus, 1984). The principle components of the ECM are collagen, and various glycosaminoglycans (GAGS). The physical chemistry of the ECM is quite complex; however, for our purposes we shall focus on one particular property: the proteoglycans comprise a poly electrolyte gel which can generate a powerful osmotic swelling pressure. The mechanical effect of the swelling pressure is to distend the ECM so that it fills a much larger volume than it would otherwise occupy. The term ‘swelling pressure’ refers to the pressure that would have to be applied to the gel to prevent it from expanding; that is, it is the total force tending to dilate the gel.

The swelling pressure can be understood by referring to Fig. 2 (c.f. Oster, 1984; Tanaka, 1983; Hill, 1960). Here a piece of ECM has been placed in contact with solvent (extracellular fluid) through a porous piston. Solvent will flow down the chemical potential gradient into the gel, causing it to swell and exert pressure on the piston.

That is, if a piece of dry gel is placed in a solvent bath it will imbibe fluid and swell. If the gel were not crosslinked, this swelling would proceed indefinitely, or until the gel expanded to fill the container. However, because they are crosslinked, the polymer strands can generate elastic forces. These forces arise from the thermal motion of the polymer fibres as they writhe about under the impact of solvent molecules. At large strains there is a small contribution from the deformation of the intermolecular bonds. Note that in equation [1] we must count the elastic pressure as negative since it tends to contract the gel.

The osmotic pressure, P_OSM, in turn has contributions from three effects: (1) the osmotic pressure arising from the mixing of the polymer with the solvent, (2) the intermolecular interactions between the polymer molecules, (3) the osmotic contribution of the counterions in the solvent (Tanaka, 1983; Hill, 1960; Flory, 1956). Of these, the dominant force is the ionic contribution; therefore, the effects of ionic strength and pH are likely to be important in regulating the hydration state of the ECM. However, for our purposes here we need only deal with the total osmotic pressure, P_OSM

The swelling pressures in the ECM derive mostly from negatively charged polymer constituents, principally hyaluronic acid (HA) and chondroitin sulphate (Grodzinsky, 1983; Alberts et al. 1983; Stern, 1985). These negative charges attract positive counterions; the effect of these counterions is to create an ‘ion pressure’ which is the major contribution to the osmotic pressure in the ECM (Tanaka, 1983). The HA macromolecule can be degraded by the enzyme hyaluronidase (HAase); the term ‘hyaluronidase’ actually refers to a class of enzymes capable of cleaving the glycosidic bonds of hyaluronic acid (c.f. Stern, 1984). The hyaluronate/hyaluronidase system admits the capacity for cells to swell and shrink their extracellular milieu: secreting HA osmotically swells the ECM, while secretion of HAase deswells it. Indeed, the hydration state of the ECM is probably determined by an equilibrium between the rates of secretion of HA and HAase.

The hydration of the ECM may exert some control over cell motility. Swelling may open sufficient intercellular space to permit cell migration. Too much hydration, however, may inhibit motion by confining the cells within a hydrated ‘coating’. Since cell migration will not play a major role in our discussion, we shall ignore this potentially major effect, although it can easily be included in the model, and will enhance its pattern-generating potential (c.f. OMH).

The major role of ECM hydration from our viewpoint is to keep the tissue inflated, and to keep the cell density low enough so that intercellular interactions are inhibited.

As discussed in OMH, cells can generate substantial traction forces on one another and on the ECM. They accomplish this by extending motile appendages such as lamellipodia and filopodia, attaching to adhesive sites and contracting (c.f. Oster, 1984). However, in this model we shall assume that cells are inhibited from protrusive activity by a superabundance of ECM. That is, when surrounded by a hydrated coat of HA-rich ECM, they are effectively isolated mechanically from their neighbours and can exert little tractions. However, if the HA coat is degraded substantially so that cells are brought into close proximity to one another, then intercellular tractions can become effective. Thus the actual cell density depends on the balance between the swelling pressure generated by the ECM and the contraction pressure produced by the cells.

It is worth noting that, as the ECM collapses and cells are brought into closer apposition, the number of intercellular contacts will not increase linearly, but will rise in a sigmoidal fashion (Perelson & Oster, 1979). That is, a decrease in tissue volume will initially increase the intercellular contacts slowly; however at some specific volume the number of contacts will suddenly jump to near its final value. This type of behaviour is characteristic of ‘phase changes’ such as the condensation of a vapour to a fluid. This phenomenon, which is purely geometrical in nature, will facilitate the sudden onset of the condensation process.

The system is in mechanical equilibrium when the condensation pressure is zero (i.e. no tendency for the volume element to change size). The formation of condensations comes about when the condensation pressure becomes negative for a time; that is, when the forces of matrix elasticity and cell tractions overpower the dilating osmotic pressure.

In the Appendix we formulate the mathematical equations that govern the condensation process. In this section we present a heuristic description of those equations so as to highlight the important physical parameters.

The model is built around the following variables which describe the condition of the limb bud in a typical volume element located at position x at time t (c.f. Fig. 3):

n(x,t) = the density of cells (chondroblasts).

m(x,t) = the density of matrix components other than hyaluronate.

h(x,t) = the density of hyaluronate.

a(x,t) = the density of hyaluronidase.

σ (x,t) = the stress (force per unit area) due to the expansive osmotic forces and the compressive elastic and traction forces.

where ϵ is the mechanical strain (i.e. the fractional deformation) and δ is a ‘diffusion coefficient’ which contains the elastic moduli of the extracellular matrix material. G(ϵ,a) is the ‘reaction’ which accounts for production and degradation of hyaluronidase, but F(ϵ,a) is a ‘reaction’ in formal terms only: it contains the osmotic, elastic and cell traction forces. Thus it has a form that has no obvious analog in chemical kinetics.

It should be emphasized that there is only a formal mathematical similarity between equations [4] and those of reaction-diffusion models based on mor-phogens. The motivation and mechanisms underlying both models are quite different. Furthermore, in our model the variables (cells, matrix, HA, HAase, and strain) are all quantities whose properties are readily measurable with conventional assays.

In the Appendix we show how these equations conspire to produce spatial patterns from an initially homogeneous distribution of cells and matrix. That is, the equations mimic the scenario described in Fig. 1:

• Cells emerge from the progress zone secreting hyaluronate. The hydration of the hyaluronate inflates the tissue, and inhibits intercellular contacts.

• After leaving the progress zone the balance between the cells’ production of hyaluronate and hyaluronidase shifts toward hyaluronidase. This initiates an osmotic deswelling of the tissue which brings the cells closer together.

• When the collapsing matrix increases the cell density so that intercellular contacts commence to increase significantly, the intercellular traction forces come into play.

• At this stage, the deswelling is sufficient so that the interplay between the osmotic swelling, tissue elasticity, and cell traction conspire to cause a condensation of cells into the dense aggregation which presages chondrogenesis.

What is not apparent from the verbal description of the model equations is why any particular spatial pattern of condensations should emerge. In particular, if condensation does occur, why should more than one appear, and why in any particular spatial arrangement?

The reason for this phenomenon is not easy to explain qualitatively. Math-ematically, the evolution of the homogeneous state (i.e. uniform cell and matrix density) to a non-homogeneous state occurs for very much the same reasons familiar from chemical systems: at some set of conditions the uniform state becomes unstable and ‘bifurcates’ into a spatially non-uniform state (c.f. OMH; Segel, 1984). Roughly, the reason is as follows.

At the onset of condensation, each volume element is in a state of mechanical balance between the expansive osmotic pressure and the contractive cell tractions. As the cell tractions come to dominate, and condensation commences, it does so about certain foci. The nature of the contraction focus is autocatalytic: once a focus of contraction nucleates, it recruits to itself cells and matrix from the surrounding tissue. However, the ‘range of influence’, or mechanical domain, of a contraction focus is limited by the elastic nature of the cell-matrix medium. Therefore, some distance away from a contraction focus, another focus can form. The distance between foci depends on the properties of the tissue, such as its elastic properties and the scale and geometry of the system. Foci need not be points; line condensations will form in a cylindrical domain. Subsequently, such line condensations can break up into smaller condensations. Indeed, the nature of the bifurcation process dictates that only certain patterns of chondrogenic condensations are admissible (Alberch et al. 1985).

According to the model, in a domain of cylindrical or elliptical cross section, there are three basic types of bifurcation, as illustrated in Fig. 4 (c.f. fig. 10 of OMH). We shall call these axial (A), transverse (T), and longitudinal (L) bifurcations. Type A bifurcations produce the first axial condensations (femur, humerus) from an initially homogeneous tissue. The splitting of the femur into the tibia and fibula is an example of a type T bifurcation, and the formation of the phalanges exemplifies a type L bifurcation.

Successive condensations can start in one of three ways. A type A condensation may appear de novo in a previously homogeneous region. A type T condensation may branch off from an existing condensation, forming a Y-shaped juncture joining the branches. This junctional region may later disappear as the cells are recruited into the main condensations. Alternatively, the junctional region may be devoid of cells (c.f. OMH). A type L condensation may also appear in two ways: as a splitting of an existing columnar condensation, or a new condensation appearing as an extension of an existing one. The model generally precludes ‘trifurcations’, i.e. triple splittings of an existing condensation. However, because the model equation [4a] is a tensor equation it is possible -albeit unlikely -that trifurcations could exist. More likely, an L bifurcation that follows a T which appears as a trifurcation may be resolved into binary processes.

As the limb bud grows, it is reasonable to suppose that the model parameters vary smoothly from the progress zone proximally. If so, one expects that the condensations will form sequentially in proximal to distal order. Therefore, the progression of condensations will appear as a sequence of T and L bifurcations. In general, we expect that at a level where two or more condensations coexist (e.g. radius and ulna) when one of the elements bifurcates (either A or T), the bifurcation of the other element will be delayed distally. This is because a condensation has a ‘domain of influence’ wherein it recruits cells into itself. An initiating condensation will tend to coopt recruitment by creating a focus of compressive stress. This will tend to inhibit nearby condensations until they grow past the stress focus and can commence their own bifurcation centre. That is, condensation distorts the stress field near the bifurcation site that inhibits the surrounding foci from bifurcating.

In a broad, flat tissue expanse, such as the distal ‘palm’ area of the forelimb, it is possible for several isolated condensations to arise independently. That is, there is enough space and tissue for separate foci to appear, and there need not be a strict proximal to distal order (Alberch, Murray & Oster, 1986). If this happens, then each focus can initiate its own progression of T and L bifurcation structures. These structures will grow outward from each centre, and there is some ambiguity as to what happens when the two islands eventually merge. We will discuss how the physics of the condensation process may constrain the possible limb morphologies in a subsequent publication. Here we shall only point out that the model puts severe constraints on the developmental process.

Each of these bifurcations can be triggered by variations in material parameters, geometry, and/or tissue size. As discussed in OMH, the quantities that control the bifurcations are dimensionless ratios of the physical parameters. Thus a bifurcation may not be attributable to a unitary cause, but may be brought about by the interaction of several effects. One of the contributions of the model is to delineate the balance of effects which may precipitate a condensation or a splitting of an existing condensation.

For example, in most species, as the limb grows distally it generally flattens into a paddle shape. In OMH we attributed this phenomenon to the tractions generated by the condensing chondroblasts. In the present model the deswelling of the ECM is a force tending to flatten the limb from its cylindrical proximal geometry into its elliptical distal cross section. At some degree of ellipticity a bifurcation can be triggered (e.g. femur to tibia-fibula). However, the bifurcation may be deferred or promoted by variations in other parameters, according to how they enter into the dimensionless ratios.

Following the formation of a condensation, the perichondrial membrane forms around the cell aggregate. If this constraining envelope is prevented from forming, the condensation frequently disappears. A clue to the possible origin of this structure lies in the orientation of the chondrocytes within the condensation.

Cells at the centre of the condensation tend to become rounded, while per-ipheral cells are flattened circumferentially (Rooney, Archer & Wolpert, 1984). The model suggests a mechanism for this. If, after condensation, the cells recom-mence secreting hyaluronate (and/or other hydrophilic matrix components) then the centre of the condensation will tend to reswell. When this happens, a pattern of strains will be set up in the ECM. This pattern of strains can be computed from elementary considerations, as described in the Appendix.

It turns out that the stress pattern in a swelling cylinder is such that the circumferential (‘hoop’) stress is twice the longitudinal (axial) stress (this is why a boiled hot dog always splits lengthwise, not transversely). Thus the reexpanding matrix in the centre of the condensation will cause fibres to align circumferentially at the periphery of the condensation. Since cells embedded in a fibrous matrix are known to align along strain directions, the reexpansion will create the observed pattern of cell geometries. If there is some mechanism that causes the flattened peripheral cells to differentiate into perichondrial tissue, then the condensation will be consolidated and will not disperse. While it is tempting to speculate on the role of cell shape in triggering such differentiations (c.f. Zanetti & Solursh, 1984), the scope of this model is limited to generating the appropriate geometrical configurations of cell aggregations.

We have described the phenomenon of matrix deswelling coupled with cell tractions as a mechanism for generating the patterns of cell aggregation that accompany chondrogenesis. The model does not depend on active cell crawling, as does the previous model discussed in OMH. However, there is no difficulty in adding cell motility to the model, and as expected, the pattern-forming capabilities are enlarged. The point of the present model is to demonstrate that cell crawling is not strictly necessary for producing regular aggregations of mesenchymal cells.

There are other instances of morphogenetic processes which commence as aggregations of mesenchymal cells. In OMH we discussed the formation of feather germs by the cell traction mechanism which depended on active cell crawling. However, the same deswelling/traction mechanism discussed here can accomplish the formation of the feather germ papillae. Indeed, it is easy to see that a combination of the two processes would be an even more effective pattern generator. Theory can only point out the physical possibilities; experiments are required to distinguish between them. If the osmotic swelling of the ECM is a crucial factor in creating cell aggregations, then experiments which modulate this force should have a profound morphogenetic effect. In this context the polyelectrolyte nature of the hyaluronate may permit experiments which modify the extracellular ionic strength and/or pH without unduly disrupting normal cell function. Antibodies to hyaluronidase are also a possible mode of intervention. We note that cytoskeletal agents such as cytochalasin which disrupt the contractile abilities of cells should be distinguishable from agents which affect the swelling pressure of the ECM. Indeed, recent experiments suggest that the profound morphogenetic aberrations produced by retinoic acid may derive from its influence on cellular production of hyaluronate (Kochhar, Penner & Hickey, 1984). From the viewpoint of the present model it is easy to see how disrupting the hyaluronate/hyaluronidase system can alter morphogenetic patterns.

Finally, the model does not address the issue of what cues the switch in the hyaluronate/hyaluronidase system which initiates condensation near the progress zone. This may well be a proximodistal gradient of some sort, or a cell lineage and/or ageing effect.

P. Alberch first called our attention to the relevant experiments on the role of ECM in chondrogenesis, and suggested modifying the original cell traction model to include the effects of matrix deswelling. Conversations with Julian Lewis were crucial to the development of the model, as were conversations with Nigel Holder, Claudio Stern and Lewis Wolpert. Albert Harris provided valuable advice and criticism. Support for this work was provided by grants from N.S.F. (MCS-8110557) to GFO and from the Science and Engineering Research Council of Great Britain (GR/c/63595). PKM would like to acknowledge the research studentship support from the Department of Education of Northern Ireland. This work was performed at the Centre for Mathematical Biology, at the University of Oxford.

The model equations

A.l The stress equation

A.2 Conservation equations for cells and matrix

The equation for cell density

The non-osmostic component of the ECM

The hyaluronate component of the ECM

Hyaluronidase (HAase)

B. l One-dimensional equations

A simplified version of the model

Dimensionless equations

B.2 Spatial patterns

B.3 Stress patterns during formation of the perichondrion

The model is built on the following field variables:

n(x,t) = the density of chondrocytes at position x and time t.

m(x,t) = the density of the non-osmotic components of the extracellular matrix at position x and time t.

h(x,t) = the density of hyaluronate at position x and time t.

a(x,t) = the density of hyaluronidase at position x and time t.

u(x,t) = the displacement of a material point initially at×at time t.

where e= Δ u+ Δu^T is the (linear) strain, E the Young’s modulus, μ is the viscosity, I is the unit tensor, and the coefficients L₁ and L₂ govern the magnitude of the second order strains, which arise from long range interactions (c.f. Oster et al. 1983). The qualitative shapes of the functions τ (·) and II (·, ·) are shown in Fig. 5. For computational purposes we choose the following convenient parametrizations

where τ, K and Π are constant parameters -physical properties of the medium which are measureable by mechanical means.

Here we assume that mitosis and matrix secretion are not major effects during the aggregation stage.

where S_h is the secretion rate of HA per cell and D_h is the rate of degradation of HA by HAase. The qualitative shapes of these functions are shown in Fig. 6A. For computational purposes we have parametrized these functions as follows:

where B_o, K_0,K₁ B₁ and K₂ are constants.

where C₀, C₁₅ and K₃ are constants. The quantities τ (cell traction) and C₀ (HAase production) are the principle parameters controlling the emergence of spatial patterns. Therefore, we shall investigate the model’s behaviour as these parameters are varied.

In [B2] we have set L₂ = 0; this will not qualitatively affect our results.

where N and M are positive constants.

Where P₀ is a constant. The rationale is as follows: (a) the presence of h invokes a production of a which in turn degrades h; this is equivalent in effect to the term — aa in [B6c]. (b) The inverse relationship between strain and osmotic pressure is contained in the term P₀/(H-ϵ).

where σ₀ is the constant stress at the boundary; we discuss this below.

where v is the production rate constant for a and co is just Ci from [A5b].

Thus we have reduced the model to a pair of equations for the strain, e, and the hyaluronidase, a. Note that these equations have the structure of a diffusion reaction system; however, the ‘reaction’ term in the mechanical equation [B7] is not one that would arise obviously from chemical kinetics.

where L, T and a₀ set the scales for size, time and HAase concentration.

We can reduce the number of parameters still further by selecting appropriate time (T) and length (L) scales. For example, if we are interested in the process on a time scale associated with the elastic response time, we choose T = γ/E; then E* = 1 in [BIO],

The possibility for generating spatial patterns by the model equations [Bll—12] can be appreciated by examining the nullclines of the spatially independent system, as shown in Fig. 7A. The mathematics which demonstrates the existence of spatially structured solutions to equations [Bl 1-12] is quite standard (see, for example, Murray, 1981). These equations are capable of generating a wide variety of spatial patterns; indeed, all of the results obtained by the previous model (c.f. the Appendix of OMH) emerge from the present model, although they are based on different physical assumptions.

Note that the spatial patterns generated by [Bl 1,12] are described by the strain (i.e. the spatial displacement of material points) and the concentrations of HAase. From [B6] the cell density may be computed directly. The scenario for chon-drogenesis which emerges from this model is identical to that described in figures 8– 10 in OMH.

Following the condensation of chondroblasts into the dense aggregations that presage cartilage secretion, the perichondrium forms. This is accompanied by a striking sequence of cell shape changes. The cells in the centre of the aggregation round up, and those around the periphery flatten circumferentially. These peripheral cells differentiate into the perichondrium -a membranous container for the cartilage capsule. Indeed, if the perichondrium is prevented from forming, the aggregation frequently disperses, and chondrogenesis fails.

The shapes of the cells at this stage may be a clue as to the physical forces the cells experience. For example, it is known that cells in a fibrous matrix will tend to align themselves along strain directions. A hypothesis that is consistent with the observed distribution of cell shapes following condensation is that following aggregation the central cells recommence secreting hyaluronate (or, what amounts to the same thing, the balance between HA and HAase secretion swings back to HA).

If this happens, the aggregation will swell osmotically, and set up a particular stress distribution in the cell-matrix medium. This stress distribution will be such that the longitudinal stresses will be about half of the circumferential stresses. That is, the ‘hoop’ stresses are twice the axial stresses (see, for example, Wainwright, Biggs, Currey & Gosline, 1976, pp. 293–294). This is true in any cylinder under internal pressure, which is why a boiled hot-dog always bursts lengthwise, rather than circumferentially.

Thus σ_θ = 2σ_L.

The consequence of this stress distribution is that cells and ECM at the periphery of the condensation will be under larger hoop stresses than longitudinal stresses. Thus alignment in the circumferential direction is expected for the same reason that cells embedded in a fibrous material will align along the direction the material is stretched. Conversely, cells in the middle of the condensation will experience a more or less isotropic stress environment, and so rehydration of the ECM will lead them to simply round up.

That the perichondrium is a tension-induced structure has been suggested previously (e.g. Wolpert, pers. comm.); we simply point out that, in the context of the present model, the cell orientations accompanying perichondrium formation can be understood in terms of the stress distribution accompanying rehydration of the chondrogenic condensation. This should be a testable hypothesis.