<<

. 43
( 66 .)



>>

Four main types of mesenchymal cell are involved in chick-limb
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

<<

. 43
( 66 .)



>>