skeletal-pattern formation. These are represented in the model by

their spatially varying (as shown in the parallelepiped de¬ned above)

and temporally varying densities. The cell types are characterized

by their expression of one of the three FGF receptors found in the

developing limb. The densities of cells expressing FGFR1 (R1 cells),

FGFR2 (R2 and R2 cells) and FGFR3 (R3 cells) are denoted respectively by

c R1 , c R2 , c R2 , and c R3 . The cells of the apical zone are R1 cells and those

of the frozen zone R3 cells (see Ornitz and Marie, 2002). The active

zone contains R2 cells and the direct products of their differentia-

tion, R2 cells. These latter cells secrete enhanced levels of ¬bronectin.

The R1 , R2 and R2 cells are mobile, while the R3 (cartilage) cells are

immobile. Thus the total density of mobile cells is c R = c R1 + c R2 + c R2 .

According to the model, the transitions and association between

the different cell types are regulated by a number of the gene prod-

ucts that we have encountered previously in early developmental pat-

tern formation: FGFs, TGF-β, a diffusible inhibitor of chondrogenesis,

and ¬bronectin. Their spatially and temporarily varying concentra-

tions are denoted c F , c A , c I , and c f respectively.

The model of Hentschel and coworkers thus comprises eight vari-

ables.

A core mechanism for chondrogenic pattern formation

The major point in formulating a model of a complex biological pro-

cess such as limb formation is to de¬ne clearly which aspects of the

8 ORGANOGENESIS 217

process one intends to investigate. The analysis of Hentschel et al.

concentrates on the mechanism that establishes the gross anatomi-

cal structure of the avian limb: the correct spatial arrangement of

its various skeletal elements, the humerus, the radius and ulna, and

the three digits (see Fig. 8.13). Once the modeling task has been de-

¬ned, one needs to identify the most relevant factors that govern

the process to be modeled. Experimental evidence suggests that the

eight variables listed above constitute a ˜˜core” set necessary to de-

scribe the development of a basic, ˜˜bare-bones” skeletal pattern (see

discussion above). These variables are intimately interconnected and

interacting: diffusing morphogen-type growth factors (i.e., FGFs and

TGF-β) and extracellular matrix molecules (i.e., ¬bronectin) secreted

by some cells represent signals for others to move, to differentiate, or

to produce other molecules or more of the same molecules. It is the

interactions among these variables (Fig. 8.14), typically described in

terms of reaction--diffusion equations (see Chapter 7) that ¬nally gen-

erate an array of possible patterns. The geometry and other physical

constraints of the limb bud select one or another of these patterns,

in the biologically relevant case the pattern of the limb skeleton. The

model of Hentschel et al. involves eight reaction--diffusion equations

(see Eq. 7.5),

‚c i D i ‚ 2c i

= + F i ({c i }),

‚t ‚ x2

where c i (i = R1 , R2 , R2 , R3 , F, A, I, f ) may denote any of the eight vari-

ables identi¬ed above, the D i are the diffusion constants (which in

the case of migrating cells may be effective parameters, as discussed

in Chapter 1), and the notation F i ({c i }) indicates that the reaction

functions may depend on any of or all the variables. (The authors

used the term ˜˜reactor--diffusion” to emphasize that the active com-

ponent is a living cell, not just a chemical reaction as in Turing™s

original formalism (Turing, 1952).)

Although the reaction functions employed in this model incor-

porate multiple cellular properties and molecular factors, much

is left out. Instead of keeping track of each of several FGFs sep-

arately, only one parameter, their global density c F , is used. Fur-

thermore, the extracellular matrix is modeled only by ¬bronectin,

which represents a marked simpli¬cation of this hydrated net-

work containing numerous types of macromolecules. Hox tran-

scription factors and nonuniformly distributed morphogens such

as Sonic hedgehog and Wnt-7A (reviewed in Tickle, 2003), which

modulate cell properties such as differential growth and the pro-

duction of cell adhesion and ECM proteins, are not considered.

These limitations do not affect the ability of this bare-bones model

to capture the essential aspect of limb formation, i.e., the proxi-

modistal generation of increasing numbers of parallel skeletal ele-

ments.

218 BIOLOGICAL PHYSICS OF THE DEVELOPING EMBRYO

Mathematical analysis of the limb bud

reactor“diffusion system

To ¬nd the general solution of a system of eight coupled partial dif-

ferential equations in a multidimensional space (spanned by all the

model parameters, such as diffusion constants and rate constants),

is a daunting task. In the present case the dif¬culty is compounded

by the fact that these equations have to describe not only the time-

dependent generation of a pattern of chemical concentrations corre-

sponding to the successive positions of prospective skeletal elements

but also the process of mesenchymal condensation, involving changes

in the local arrangement of mesenchymal cells (as discussed in Chap-

ter 6). Proving the existence of unique symmetry-breaking solutions

to a system of this complexity, a highly nontrivial endeavor, has been

accomplished by Alber et al. (2005). Their result justi¬es considering

the behavior of this model in various approximations.

A number of useful simpli¬cations are possible in formulating a

model for this process without signi¬cantly departing from the un-

derlying physics and biology. Speci¬cally, we will illustrate how the

separation of time scales and linear stability analysis was used to re-

duce the number of equations in the model of Hentschel and cowork-

ers. We will then indicate how the outlined mathematical formalism

leads to the desired pattern.

The equation that describes the spatiotemporal changes in FGF

concentration in the model has the form

‚c F ‚ 2c F

= D F 2 ’ kF c F + J F (x, t). (8.8)

‚t ‚x

To simplify the notation, we will indicate the spatial variable in the

proximodistal direction only, noting however that the equations must

be solved in at least two dimensions. (Variations in the dorsoventral

direction, i.e., the z direction, are ignored in this model; this is equiv-

alent to assuming that the cross-sections of the bones are symmet-

rical along that axis, not a serious departure from biological real-

ity.) In Eq. 8.8, kF and J F (x, t) are respectively the decay constant and

the ¬‚ux of FGFs. The latter quantity represents the overall (cellular)

source of the growth factors. Provided that kF and J F (x, t) are given

and that the boundary conditions on the surface of the developing

limb parallelepiped are speci¬ed for the function c F (x, t), the above

equation, in principle, uniquely determines the spatial distribution of

FGFs everywhere in the growing limb at any instant of time. (In their

two-dimensional analysis Hentschel et al. assumed zero-¬‚ux bound-

ary conditions on the anteroposterior borders of the limb, i.e., that

morphogens cannot leak out of those surfaces.) As discussed above,

it is the AER that is the major source of FGFs. Therefore, as far as

FGF production is concerned, this structure can be represented by a

steady boundary ¬‚ux of these molecules into the apical and active

zones. The diffusion of FGFs is assumed to be a ˜˜fast process” in com-

parison to the growth of the limb; that is, the motion and eventual

condensation of cells in the active zone is assumed to take place in

8 ORGANOGENESIS 219

the presence of an already equilibrated FGF distribution determined

by the simpli¬ed equation

d2 c F kF

= c F. (8.9)

2

dx DF

This equation, which describes the (time-independent) FGF density

pro¬le in the active zone, is supplemented by the boundary condi-

tion at the distal tip of the active zone, D F dc F /dt = j, giving the

in¬‚owing FGF diffusion current (see Chapter 1); j has a ¬xed value

that is related to the equilibrium value of the ¬‚ux J F at the bound-

ary. The approximation used to arrive at Eq. 8.9 from Eq. 8.8 is called

the separation of time scales. Since Eq. 8.9 is an ordinary differential

equation, its solution is relatively simple. Once c F (x) is determined

from Eq. 8.9, it can be inserted into the other equations in which it

appears, thus leading to fewer coupled reaction--diffusion equations.

To illustrate how the dynamics of cells is incorporated into the

model of Hentschel et al., we consider the R2 cells. The dynamics of

these cells includes diffusion, haptotaxis, mitosis, and differentiation.

It is governed by the following ˜˜reactor--diffusion” equation:

‚c R2 ‚ 2c R ‚c R2 ‚c f

= D R2 2 2 ’ χ + r c R2 (c eq ’ c R2 )

‚t ‚x ‚x ‚x

+ k12 c R1 ’ k22 c R2 . (8.10)

The second term on the right-hand side describes condensation result-

ing from haptotaxis up the ¬bronectin concentration gradient and

thus depends on ‚c f /‚ x; it also depends on ‚c R2 /‚ x, because this fac-

tor characterizes the developing pattern of R2 cell condensation. The

factor χ is a constant. Mitosis is modeled by the term r c R2 (c eq ’ c R2 ).

This speci¬c form implies that an equilibrium density of the R2 cells

is stable in the absence of mesenchymal condensation. Such an equi-

librium would come about as a result of the suppression of mitosis at

high cell densities, c R2 > c eq . The last two terms describe, respectively,

the differentiation of R1 cells into R2 cells, and the differentiation of

R2 cells into R2 cells. Since cells differentiate in response to molec-

ular signals, the rate constants k 12 and k 22 are functions of FGF and

TGF-β (for further details see Hentschel et al., 2004). The equations

governing the dynamics of the other molecules (TGF-β, inhibitor, and

¬bronectin) and cell types are constructed in a similar manner, by

taking into account their roles and biological speci¬city.

The separation of time scales and the use of conveniently scaled

density functions ¬nally reduces the number of variables in the model

of Hentschel and coworkers to four, the concentrations of the three

morphogens FGF, TGF-β, and the inhibitor and the density of the

total mobile cell population R . These were ¬rst analyzed by using the

method of linear stability analysis (for details of this method see the

Appendix to Chapter 5 and for a speci¬c application see the model

of juxtacrine signaling of Sherratt et al. in Chapter 7). As discussed in

Chapter 7, linear stability analysis establishes which spatial patterns

might develop as possible solutions to the system of equations. These

220 BIOLOGICAL PHYSICS OF THE DEVELOPING EMBRYO

patterns are identi¬ed as modes that are exponentially growing in

time; they are spatially and temporally varying linear perturbations

to the homogeneous solutions. Such modes can be selected by varying

the ratio of the length of the active zone and the dimension along

the anteroposterior axis, l(t) = l x (t)/l y . Biology thus dictates how l(t)

should change with time. If a skeletal pattern with one skeletal el-

ement (i.e., the humerus) is to be followed by two skeletal elements

(i.e., the radius and ulna), and then by three (i.e., the digits), then the

choice of model parameters should allow l(t) to vary in such a way as

to make the modes representing the various skeletal elements grow

exponentially.

The function l(t) was constrained, but not uniquely determined,

by realistic values for the model parameters. Numerical integration

of the three nonlinear reactor--diffusion equations with plausible

choices for the model parameters and the time variation of l(t) led

to a density pro¬le of mobile cells that varied inversely, and in a

˜˜quantal” fashion, with l x , the proximodistal length of the active

zone. Speci¬cally, one, two, and three stripe-like distributions of con-

densed cells (representing successive modes of the reactor--diffusion

system under the chosen parameter and boundary values) formed as

the active zone length diminished continuously with time. (Note that

the nonlinear dynamics changes the prediction of the linear stability

analysis in that the modes saturate in concentration instead of con-

tinuing to grow at an exponential rate. For details of the numerical

analysis, see Hentschel et al., 2004).

To visualize the actual pattern of cartilage elements produced by

the model over the duration of limb outgrowth, simulations of this

system were performed with different initial conditions of activator

concentration. The simulations were carried out under the biologi-

cally justi¬ed assumption that the elements in the frozen zone keep

growing while successive elements in more distal regions are being

generated. The results of three of these two-dimensional simulations

are shown in Fig. 8.15A with, for comparison, a diagram of a section of

the actual skeletal pattern in the plane spanned by the proximodistal

and anteroposterior axes at the completion of chondrogenesis (see Fig.

8.13). Although all the simulations look limb-like, the dependence of

the pattern™s details on initial conditions point to the incompleteness

of the model. The ˜˜robustness” (developmental stability; see Chapter

10) of embryonic limb development depends on a multiplicity of bio-

chemical interactions (see Tickle, 2003) that extend beyond the bare-

bones mechanism modeled by Hentschel and coworkers.

Simulations of the biologically more realistic three-dimensional

version of this model are computationally much more demanding, re-

quiring a ˜˜multimodel” computational framework (Chaturvedi et al.,

2005; Cickovski et al., 2005). Using a simpli¬ed version of the sys-

tem of Hentschel and coworkers to model the activator and inhibitor

chemical ¬elds, Chaturvedi et al. (2005) obtained the time series

of pro¬les of the activator (i.e., TGF-β) concentration across succes-

sively distal positions in the developing limb shown in Fig. 8.15B.

Fig. 8.15 Simulations of limb skeletal development in the model of Hentschel et al.

(2004). (A) On the left, a longitudinal section of the skeleton of the chicken limb at seven

days of development; anterior is at the top of the diagram and posterior is at the

bottom. On the right, typical examples of skeletal structures comparable to the

seven-day pattern generated by the model, using different initial conditions. The cell

density in the active zone was calculated from Eq. 8.10 and then ¬xed for cells leaving the

active zone and entering the frozen zone (see Fig. 8.14). Growth in all zones was

assumed to occur at a constant rate. Earlier-formed cartilage elements were thus subject

to more growth than later ones. The distribution of cartilage is shown in a continuous

grayscale in the simulation panels, black representing the highest cartilage density. The

skeletal form in the model is dependent on the parameter values and time-dependent

changes in the active zone (which were the same for all three simulations shown) and on

the initial conditions (which differed, yielding slightly different patterns). (B) Time series

of the concentration of the diffusible morphogen TGF-β displayed in cross-sections of

the active zone at successive stages of development, with time increasing in the upward

direction, in a three-dimensional model based on the analysis of Hentschel et al.

(Chaturvedi et al., 2005). (C) A three-dimensional simulation of cell distribution during

chick-limb skeletal patterning, based on a modi¬ed version of the model of Hentschel

et al. (2004), using the CompuCell3D multimodel simulation framework (Cickovski et al.,

2005). Whereas the original model represented the cell density as a continuous variable,

this simulation used the cellular Potts model (see Chapter 6) to represent cell position

and motion. Three successive developmental stages are shown, the ¬nal stage

corresponds to those in panel A. Cells that have undergone condensation are shown

here in gray. (The simulations shown on the right in panel A are from Hentschel et al.,

2004, by permission. Panel B is from Chaturvedi (2005) by permission. Panel C is given

courtesy of T. Cickovski and J. A. Izaguirre, University of Notre Dame).

222 BIOLOGICAL PHYSICS OF THE DEVELOPING EMBRYO

Cickovski et al. (2005), using a slightly different version of the model

of Hentschel et al. (2004) in conjunction with an energy-minimizing

cellular automaton model for cell condensation like that discussed in