Bulletin of Mathematical Biology (1998) 60, 126

Pattern Formation in a Generalized Chemotactic Model

M. R. MYERSCOUGH
School of Mathematics and Statistics,
F07, University of Sydney,
N.S.W. 2006, Australia
P. K. MAINI AND K. J. PAINTER
Centre for Mathematical Biology,
Mathematical Institute,
2429 St Giles,
Oxford OX1 3LB, U.K.
Many models have been proposed for spatial pattern formation in embryology
and analyzed for the standard case of zero-flux boundary conditions. However,
relatively little attention has been paid to the role of boundary conditions on
the form of the final pattern. Here we investigate, numerically, the effect of
nonstandard boundary conditions on a model pattern generator, which we choose
to be of a cell-chemotactic type. We specifically focus on the role of boundary
conditions and the effects of scale and aspect ratio, and study the spatiotemporal
dynamics of pattern formation. We illustrate the properties of the model by
application to the spatiotemporal sequence of skeletal development.
c 1998 Society for Mathematical Biology

1.

I NTRODUCTION : M ATHEMATICAL M ODELS FOR PATTERN

F ORMATION IN D EVELOPMENTAL B IOLOGY

One of the main areas of research in developmental biology seeks to understand

the key processes and mechanisms that underlie the development of structure and
form within an embryo. A number of mathematical models have been proposed
in this context. Broadly speaking, they fall into two main classes: chemical
prepattern modelsin which a uniform density of cells responds to a chemical
prepattern and differentiates accordingly, and cell movement modelsin which
the cells themselves form a spatial pattern with those cells which are in highdensity aggregates differentiating.
The two main chemical prepattern models are gradient models and reaction
diffusion models. In the former, a simple gradient in some chemical concentration
is set up and this is interpreted by cells which differentiate according to a complex mechanism involving multiple thresholds (Wolpert, 1969). In the latter, the
complex interaction of reaction and diffusion in a system of chemicals leads to
a pattern of peaks and troughs in chemical (morphogen) concentration (Turing,
.

Author to whom correspondence should be addressed.

0092-8240/98/010001 + 26

$25.00/0

bu970010

c 1998 Society for Mathematical Biology

M. R. Myerscough et al.

1952) which, it is hypothesized, is then read off by cells which differentiate

if the morphogen concentration lies above a certain threshold. Therefore, one
of the main differences between these models is that in the former the chemical
prepatterning mechanism is simple but the interpretation mechanism is complex,
while the reverse is true for the latter (Nagorcka, 1989).
The two widely known examples of cell-movement models are mechanochemical models and chemotactic models. In mechanochemical models (Oster et al.,
1983), pattern arises as the result of cell traction exerted by mesenchymal cells
on the extracellular matrix. This cell traction destabilizes the uniform steady state
which results in a pattern of peaks and troughs in cell density. In chemotactic
models, cells secrete a chemical (chemoattractant) and move up a gradient in
this chemical. It is well known that this is a mechanism for generating spatial
patterns in cell density (Keller and Segel, 1970). In cell-movement models, the
pattern is in cell density and it is assumed that cells which are in high-density
aggregates then differentiate (although this is not explicitly included in any of
these models). Murray (1989) reviewed in detail both types of model and their
application to a wide range of developmental patterning phenomena.
There are many similarities between reactiondiffusion, mechanochemical and
chemotactic pattern generators. All are based on the interaction between an activating mechanism, which is promoting growth (of chemical concentration or
cell density) and an inhibiting mechanism, which depresses growth. The uniform
steady state is destabilized when the balance between these two mechanisms is
tipped in favor of the activating mechanism due to, for example, greater diffusion of inhibitor than of activator in the case of reactiondiffusion. A spatially
homogeneous steady state then evolves into a spatially heterogeneous state. This
is a generalization of diffusion-driven instability, which applies specifically to
reactiondiffusion systems. The linear analyzes of these models are all similar
mathematically and lead to the prediction of bifurcation from spatially uniform
solutions to spatially nonuniform solutions. In the vicinity of primary bifurcation
points, these spatially nonuniform solutions are eigenfunctions of the Laplacian
with appropriate boundary conditions on the domain. This naturally leads to the
concept of developmental constraints, where certain key features of spatial patterns hold independently of the mechanism which generated them (Oster et al.,
1988).
Models of this form have been widely studied for the standard case of a homogeneous domain, that is, the parameters are constant across the whole domain,
with zero-flux boundary conditions. Recently, however, there have been a number of papers devoted to the analysis of such models on asymmetric domains,
where certain key parameters show a spatial or spatio-temporal dependence (see,
for example, Hunding and Brns, 1990; Maini et al., 1992; Kulesa et al., 1995).
The role of generalized boundary conditions on the spatial patterns exhibited by a
standard reaction diffusion model has also been analyzed by Dillon et al. (1993).
Here we focus on a cell chemotactic model for pattern formation. In Section 2
.

Pattern in a Chemotactic Model

we briefly review pervious results on general cell chemotaxis models. In Section 3

we introduce the specific model that we analyse. In subsequent sections we
investigate the role of boundary conditions and the effect of scale on the patterns
exhibited by the model using numerical methods. The one-dimensional results
provide us with crucial insight into the behavior of the model in two dimensions,
which is discussed in Section 6. A possible application to skeletal patterning in
the vertebrate limb is presented in Section 7.

2.

C HEMOTACTIC M ODELS AND PATTERN F ORMATION

Much of the recent work on chemotactic models has concentrated on waves

of migrating cells or other entities in various contexts. These include angiogenesis (Chaplain and Stuart, 1993), various cell migration assays (Ford and
Lauffenberger, 1991, Sherratt et al., 1993), propagating patterns (Myerscough
and Murray, 1992) and even herd grazing (Gueron and Liron, 1989). Rather less
attention has been given to chemotactic pattern formation. The seminal work of
Keller and Segel (1970) addressed the problem of pattern formation among aggregating Dictyostelium amoebae. More recent models for this phenomenon couple
the excitable dynamics of the signaling chemical cAMP and cell receptors to the
chemotactic cell response to cAMP to produce the streaming instability characteristic of Dictyostelium aggregation (see, for example, Ho fer et al., 1995a,b).
Chemotactic models have also been proposed for skin pigment patterns in reptiles, specifically stripe patterns in alligators (Murray, et al., 1990) and snake-skin
markings (Murray and Myerscough, 1991).
The basic general chemotaxis model takes the form:
.

n
= (Dn n) (n f (n, c)c) + h(n, c)
t
c
= (Dc c) + g(n, c)
t

(1)
(2)

where n is cell density and c is the concentration of chemoattractant; Dn and Dc

are diffusion coefficients, f (n, c) models the chemotactic response of the motile
cells, h(n, c) models cell division and cell death, and g(n, c) models chemoattractant production and degradation. This type of model was first proposed by
Keller and Segel (1970). Typical boundary conditions for a pattern-formation
model are traditionally periodic or zero flux. The effect of different types of
boundary conditions has not been addressed for this model.
The form of the function f (n, c) which controls the chemotactic response
of the cells to the chemoattractant ultimately depends on the sensitivity of the
cells at different concentrations of the attractant. The simplest form for f (n, c)
is f (n, c) = , a constant (see, for example, Maini et al., 1991; Chaplain
and Stuart, 1993). This assumes that the sensitivity of cells to attractant is independent of chemoattractant concentration. If cell sensitivity to chemoattractant
.

M. R. Myerscough et al.

concentration is known to decrease as concentration increases then one possible

model is f (n, c) = /c, where is constant (Keller and Segel, 1971). The action of cell-surface receptors has been modeled by setting f (n, c) = k/(k + c)2 ,
where k and are constants (Lapidus and Schiller, 1976; Ford and Lauffenberger, 1991). The above models are phenomenological. Recently, Othmer and
Stevens (1997) derived macroscopic chemotaxis equations based on microscopic
rules using a random walk approach. They showed which specific assumptions
of how a cell moves give rise to the types of phenomenological macroscopic
models above.
The coefficient Dn that governs the random component of cell motility may
depend on chemical concentration. For example, to model chemokinesis, Dn
would be an increasing function of c (Sherratt et al., 1993). For simplicity,
however, we will asssume that Dn is constant.
The function h(n, c) models cell proliferation and death. If cell numbers are
conserved, i.e. h(n, c) = 0, then in the steady state one can integrate equation (1)
and show that the final steady state is periodic or consists of a single peak or half
peak (Grindrod et al., 1989). When h(n, c) is nonzero, cell growth is usually
assumed to be logistic; i.e. h(n, c) = r n(1 n/n 0 ), where r and n 0 are usually
constants.
The production, degradation and consumption of chemoattractant are given by
the function g(n, c). The degradation or consumption term in g(n, c) is typically
linear and the production term is of the MichaelisMenten type (see, for example,
Keller and Segel, 1970; Myerscough and Murray, 1992; Maini et al., 1991;
Chaplain and Stuart, 1993; Sherratt et al., 1993). This form of g(n, c) and the
chemotactic flux of cells together form the activating processes in a chemotactic
pattern-formation model. This is in contrast to reactiondiffusion models where
the presence of one of the chemical species activates the pattern under suitable
conditions. In the next section we present the basic chemotactic model for pattern
formation that we shall analyze.
.

3.

M ODEL E QUATIONS

In order to investigate the role of boundary conditions and scale effects on

patterns in cell density we consider the model analyzed in Maini et al. (1991).
The basic two-species model is


n
n
= (Dn n) (nc) + r n 1
(3)
t
n0
.

n
c
= (Dc c) +
c
t
n+

(4)

where n is cell density and c is the concentration of chemoattractant; is the

chemotactic coefficient of the motile cells, r is the linear growth rate of the cell

Pattern in a Chemotactic Model

population, n 0 the cell density under steady-state conditions and , , and are
constants governing the rate of chemoattractant production and degradation. We
will assume that Dn and Dc are positive constants.
Equations (3) and (4) can be nondimensionalized by setting
n =

c = c

n
n0
Dn
Dc

D =

r =

Dc s

Dc

x =
=

 12

t
s

t =

x
=

.
n0

(5)

Dropping the asterisks, the nondimensional equations then become

n
= (Dn) (nc) + sr n(1 n)
t


c
n
= (c) + s
c .
t
n+

(6)
(7)

Here, s is a parameter which controls spatial and temporal scale.

We apply generalized boundary conditions

1

n
c
D
n
r
r
3


= (1 1 )(2 n s n)

(8)

c
= (1 3 )(2 cs c)
r

(9)

.
where r
is the normal derivative to the boundary, 1 , 2 and 3 are nonnegative
control parameters and n s and cs are the nontrivial values of n and c, respectively,
for which the kinetic terms vanish, i.e. n s = 1, cs = 1/(1+ ). With the boundary
conditions in this form we can impose very general boundary conditions by
varying the control parameters i , where i = 1, 2, 3.
Using homotopy methods, Dillon et al. (1993) analyzed a two-chemical reactiondiffusion model in one dimension with boundary conditions of the above
form. We shall investigate the effect of different boundary conditions on the
chemotactic pattern. In Sections 4 and 5 we present numerical solutions to the
model equations (6) and (7) with boundary conditions (8) and (9) in one spatial
dimension on the domain 0 x 1 using NAG solver D03PHF to obtain timedependent solutions and the ENTWIFE finite element package (Winters, 1985)
for steady-state solutions.
.

4.

R OLE OF B OUNDARY C ONDITIONS

The above model was extensively studied by Myerscough (1988) for the case of
no flux boundary conditions on both n and c (i.e. when 1 = 3 = 1 in equations
(8) and (9)). The model can exhibit patterns of peaks of cell density which overlie
.

M. R. Myerscough et al.

peaks of chemoattractant concentration. Depending on the parameter values and

initial conditions used, these peaks can be internal or can lie on the boundary.
Although such patterns are appropriate for applications such as skin pigmentation
patterns, they are not appropriate for cases in which patterns are internal to the
domain, for example, in skeletal development in the limb. Here we investigate
the role of alternative boundary conditions.
4.1. The role of symmetric mixed boundary conditions. We now impose mixed
boundary conditions on the system of equations (6) and (7). Specifically we take
1 = 1, 2 = 3 = 0 in equations (8) and (9) so that c = 0 on both boundaries,
while n still has no flux boundary conditions. These mixed boundary conditions
eliminate the boundary peaks. This can be seen in Fig. 1. The stability of these
solutions can be determined either by solving the time-evolution problem or by
analyzing the steady-state problem using the ENTWIFE numerical bifurcation
package.
If cells are present, i.e. n 6= 0, then c will not be zero except at the boundary. Hence, near each boundary there is a gradient in chemoattractant c which
increases away from the boundary. Cells move up this gradient and away from
the boundary. This forces pattern to occur within the domain and also eliminates
boundary peaks. Dillon et al. (1993) found that mixed boundary conditions can
also internalize peaks in the standard Turing model equations.
Symmetric boundary conditions always result in a symmetric pattern in the
sense that all peaks have the same amplitude. We next explore how asymmetric
boundary conditions give rise to asymmetric patterns.
.

4.2. The role of asymmetric mixed boundary conditions. To impose asymmetric mixed boundary conditions the values of the parameter set {1 , 2 , 3 }
in equations (8) and (9) must be different at x = 0 from those at x = 1.
We shall investigate asymmetric boundary conditions on c. Therefore, we fix
1 = 1 to impose zero-flux conditions on n. For c we choose 2 = 3 = 0
at x = 0. At x = 1 we set 3 = 0 and 2 > 0. With these parameter values we have c(0) = 0, c(1) = 2 cs . The solutions to equations (6)
and (7) with these boundary conditions are spatially asymmetric. Furthermore,
the solutions evolve in a specific temporal fashion. This temporal sequence is
shown in Fig. 2. The asymmetry in the boundary conditions for c imposes
an initially large gradient in c at the boundary x = 0. This initiates a peak
in cell density near this boundary. If further peaks form, they do so in a
temporal sequence from x = 0 to x = 1; that is, the first peak to form is
nearer the x = 0 boundary provided that 2 is close to 1. Their amplitudes
are different and their maxima do not occur in regular intervals along the domain. Thus, we have a genuine spatial asymmetry. If c(1) > cs , i.e. 2 > 1,
there is an outward gradient in c at x = 1 initially, as c reaches cs internally. This results in a peak in cell density on the boundary. If c(1) < cs ,
( < 1) then there is an inward gradient in c and peaks in n do not form on the

Pattern in a Chemotactic Model

A
1.5

1.0

0.5

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

x
B
1.5

1.0

0.5

0.2

0.4
x

Figure 1. Numerical simulations of the cell-chemotactic model equations (6) and (7)
in one dimension for different boundary conditions, with parameters D = 0.25, r =
0.04, = 2, = 1, s = 600. For the initial conditions we imposed random perturbations
of maximum size 0.01 about the uniform state n = 1 for cell density, and set c to the
uniform state of 0.5. The time evolution equations were solved using NAG library
routine D03PHF until a steady state was reached: (a) Zero flux boundary conditions,
(i.e. in equations (8) and (9) 1 = 3 = 1). (b) Mixed boundary conditions with 1 =
1, 2 = 3 = 0 in equations (8) and (9) so that c = 0 on the boundary. In each case
the final steady state is shown. Note that there is a boundary peak in (a), but in (b) the
pattern is internalized. Note also that the pattern in (b) is simpler (i.e. it has fewer
peaks), than that in (a). n, ; c, .

M. R. Myerscough et al.

A
1.5

1.0

0.5

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

x
0.6

0.4

0.2

0
0

0.2

0.4
x

Figure 2. Time evolution for the model equations (6) and (7) in one dimension with zero
flux boundary conditions on n [1 = 1 in equation (8)], but with asymmetric boundary
conditions on c. Specifically, for the left-hand boundary x = 0, 2 = 0, 3 = 0 in
equations (8) and (9), while for the right-hand boundary x = 1, 2 = 0.9, 3 = 0.
With the parameters D = 0.25, r = 0.04, = 2, = 1, s = 1000, this corresponds to
c(0) = 0, c(1) = 0.45. Solutions are shown for (a) cell density n, and (b) chemoattractant
density c. Initial conditions as in Fig. 1. t = 200, ; t = 2000 ; steady
state, .

Pattern in a Chemotactic Model

boundary (Fig. 3). As c(1) is decreased to zero the spatial asymmetry in n is

lost.

5.

T HE E FFECTS OF S CALE

For both symmetric and asymmetric boundary conditions, the nature of the
pattern, in particular the number of peaks within the pattern, will depend on the
domain size, or in the non-dimensional equations (6) and (7) on the parameter
s. Increasing s effectively increases the dimensional size of the domain and the
dimensional evolution time.
The time-evolution problem was solved for various values of s with both symmetric and asymmetric mixed boundary conditions; a selection of these solutions
are presented in Section 4. For a more comprehensive picture of how domain
size influences pattern, bifurcation diagrams were generated for steady-state solutions of equations (6) and (7) in one dimension with different sets of boundary
conditions. These are shown in Figs 47. When mixed boundary conditions are
used for equations (6) and (7), no nontrivial uniform steady state exists and a
nonuniform steady state always exists even for very small s. Where both n and
c have no flux boundary conditions, and hence a nontrivial, uniform steady state
exists, the bifurcation diagram is likely to have many branches which bifurcate
from the uniform steady state (cf. Maini et al., 1991; Dillon et al., 1993). For
mixed boundary conditions there is one unbranched curve for most values of s.
In Figs 4 and 6 the measure of the steady-state solution is taken as the average
of n over the domain, i.e.
.

knk =

n(x) d x.
0

For the symmetric mixed boundary conditions of Section 4, with parameter

values D = 0.25, = 2, r = 0.04, and = 1, increases in the number of peaks
take place via a series of limit points. This is very different from the case in
Dillon et al., (1993), where transitions from one type of pattern to another can
occur in the absence of a bifurcation. For example, for s < 450 the solution
has only one peak. As s increases, this peak becomes indented at the top. This
indentation deepens on rounding the limit point at s = 544.11 to give an unstable
solution with two distinct peaks. At the next limit point, s = 505.52, this twopeak solution becomes stable. Thus, this fold in the bifurcation curve is associated
with the formation of a minimum in n which separates the split peak into two
distinct peaks. The time-evolved solution tends to the split-peak solution which
corresponds to the upper part of the fold. This is usually, but not always, the
case (Figs 4 and 5).
The second fold in the bifurcation curve in Fig 4 is associated with the twoto three-peak transition. Here, the third peak begins to form in the minimum
.

10

M. R. Myerscough et al.

A
1.5

1.0

0.5

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

x
B
1.5

1.0

0.5

0.2

0.4
x

Figure 3. Steady-state solutions of equations (6) and (7) in one dimension with zeroflux boundary conditions on n [1 = 1 in equation (8)], but with asymmetric boundary
conditions on c: c(0) is kept fixed at 0 but different boundary conditions are imposed
at x = 1: (a) 2 = 0.9, 3 = 0 in equation (9), (b) 2 = 1.1, 3 = 0 in equation (9).
Parameters and initial conditions are as in Fig. 1. n, ; c, .

Pattern in a Chemotactic Model

0.8

11

E
G

0.7
||n||

B
A

0.6

0.5
0

1000

2000
s

3000

Figure 4. Bifurcation diagram for steady-state solutions of equations (6) and (7) with
no flux boundary conditions on n and homogeneous Dirichlet boundary conditions on c.
Parameter values used are D = 0.25, = 2, r = 0.04 and = 1.

between the existing two peaks. The formation of this peak is associated with a
Hopf bifurcation in which the primary branch becomes unstable. (We will not
investigate Hopf bifurcations further in this paper as our main concern is with
stationary pattern.) Moving around the two limit points of this fold, the minima
in n deepen to give three distinct peaks of roughly equal height.
As s increases toward the next fold, the peak in the center of the domain
becomes split, and as for the first fold, this split deepens as the limit points at
s = 2028.99 and s = 1672.95 are traversed to give four distinct peaks.
The bifurcation diagram becomes more complex for s > 2000. The next limit
point, following the curve from the third fold, is close to another Hopf bifurcation
which is associated with the formation of a fifth peak at the center of the domain.
As before, the minima on either side of this peak deepen on rounding the limit
point to give a five-peak pattern. These straightforward four- and five-peak
patterns continue to be the attractors for the time-evolved solutions despite the
existence of less regular patterns on other parts of the bifurcation curve between
s = 2000 and s = 3100.
The bifurcation diagram for mixed asymmetric boundary conditions is shown in
Fig. 6 with the corresponding steady-state solutions for n being shown in Fig. 7.
There are no Hopf or secondary bifurcations for s < 3000 in this case. The first
three folds are associated with the growth of minima in n to give new peaks.
Here, because of the boundary conditions, the new peak always forms on the
right-hand side of the domain where c = 0.45 at x = 1. (On the left-hand side
c = 0 at x = 0.) The fourth and fifth sets of limit points, with their associated
unstable regions at approximately s = 1900 and 2500, are associated with a slight
leftward movement of the four peaks rather than the generation of extra peaks.
In this particular asymmetric case, when three solutions exist, the time-evolved
solutions are attracted to the steady-state solutions on the upper part of the curve.

12

M. R. Myerscough et al.
A

Figure 5. Steady-state solutions in n corresponding to marked points on Fig. 4 (the

corresponding chemical concentration profiles are qualitatively similar to those of n).

This is not the case for every parameter set. For example, when D = 0.35,
= 2, r = 0.01 and = 1, the first three folds in the bifurcation curve appear
similar to those in Fig. 6. The time-evolved solutions, however, tend towards
solutions on the lower part of the fold in two out of the three cases.
In this example, the case with asymmetric boundary conditions differs from the
symmetric case in that the new peaks form on one side of the domain rather than
at its center and there are no Hopf or secondary bifurcations. This lateral (as
opposed to central) formation of new peaks in cell density, as dimensional domain
size increases, distinguishes pattern formed in asymmetrical environments from
pattern formed in strictly symmetrical domains.
Application of this type of model to a developmental phenomenon implicitly
assumes that cells at high density differentiate so that the pattern in n exhibited
by the model is the pattern observed in development. (There is evidence that in
limb development, for example, cells in high-density aggregates will differentiate
(see Maini and Solursh, 1991, for review).) Note that there is the possibility that
transient solutions that are different from the final steady state may have an effect
on the pattern of differentiated cells. However, under the biologically realistic
assumption that differentiation occurs on a time scale slower than that of pattern
formation, one can show that transient solutions in n have only a negligible effect
.

Pattern in a Chemotactic Model

C

13

0.8
G
H

||n||

D
B

0.7

0.6
0

1000

2000

3000

s
Figure 6. Bifurcation diagram for steady-state solutions of equations (6) and (7) with no
flux boundary conditions on n and asymmetric boundary conditions on c; that is, c = 0
at x = 0 and c = 0.45 at x = 1. Parameter values used are D = 0.25, = 2, r = 0.04
and = 1.

on the final solution profile of differentiated cells. Differentiation is essentially a

slave process and does not influence the qualitative form of the n and c solution
profiles. Therefore, we are justified in analyzing only the steady-state problem
for the n c model and then inferring the differentiated state from the n profiles.

6.

T WO -D IMENSIONAL S OLUTIONS

We now consider two-dimensional domains with different types of boundary

conditions. Of course, there is now much more variety in the different types of
conditions that may be imposed. To investigate how the one-dimensional results
generalize, we choose to focus on two types of boundary condition, which we
denote by Type A and Type B (see Fig. 8). In both boundary conditions, we
impose zero flux for cell density but vary the conditions for the chemical. In
Type A, we have zero-flux boundary conditions along two of the boundaries and
a source and sink at the other boundaries. More specifically, at the left-hand
boundary we set c = 0 and at the right-hand boundary c = 0.45. For Type B,
the sink extends to three of the boundaries and there is a chemical source at one
end.
We define the aspect ratio a as L/W (Fig. 8), and in all the simulations
presented below we fix the parameter values of D, r, and as for the onedimensional case and vary a and s. The initial conditions for the chemical
concentration and cell density are set at the respective steady-state levels with a
small random perturbation (1%) everywhere for the cell density.
We first consider Type A boundary conditions and investigate the extent to
which stripes are maintained as we move from one to two dimensions. We find
that as the parameter value s increases, stripes are maintained only if a lies above

14

M. R. Myerscough et al.
A

Figure 7. Steady-state solutions in n corresponding to marked points on Fig. 6.

a critical value, which can be determined from linear analysis (see Appendix).
We illustrate the results in Fig. 9. Intuitively, this is as we would expect from
analogy with the standard zero-flux case. In the latter, linear analysis predicts the
formation of spots or stripes, depending on the parameters and initial conditions.
Spots are effectively the superposition of stripes of different orientations. As a
increases, the space to form multioriented stripes decreases, and at a critical value
of a, stripes of only one orientation will form.
Figure 10 shows the time evolution of the pattern for the parameters s =
900, a = 1. Initially, the pattern is essentially one dimensional with three parallel
stripes. This transient solution eventually breaks up into a mixed stripe/spot
pattern.
The solutions for Type B boundary conditions are shown in Fig. 11 for increasing s. We increase a with s so that the domain width is kept constant. As s
increases, a central spot of high-cell density splits into two spots and then a third
spot emerges between these pre-existing spots.

Pattern in a Chemotactic Model

15

Chemical source

Chemical sink
Zero flux

W
Type A

Type B

Figure 8. Schematic representation of the boundary conditions imposed on c for the

two-dimensional simulations (see text for details). Boundary conditions for n are zero
flux everywhere.

1.5

1
0.5
0
A

0.5

0.5

0.5

1
B
0

1
C
0

Figure 9. Steady-state solutions for n under Type A boundary conditions for varying s
and a. (a) s = 300, a = 1; (b) s = 400, a = 5.0; (c) s = 900, a = 7.41. Corresponding
one-dimensional cross sections illustrating cell density are shown in (d)(f). [White
denotes high cell density (n > 1.1), black denotes low cell density (n < 1.0)].

16

M. R. Myerscough et al.

Figure 10. Time evolution for n under Type A boundary conditions for s = 900,
a = 1.0. (a) t = 0.6, (b) t = 1.0, (c) t = 2.0, (d) t = 4.0. [White denotes high cell
density (n > 1.1), black denotes low cell density (n < 0.8)]. Here the sink for c is at
the top of each square and the source is at the bottom (that is, the figure on the left-hand
side of Fig. 8 has been rotated through 90 clockwise).

7.

A PPLICATION TO L IMB D EVELOPMENT

Although a great deal of experimental and theoretical work has been carried
out on skeletal patterning in the chick limb, several key questions remain to
be answered (see Maini and Solursh, 1991, for review). Specifically, there is
no theoretical model, to our knowledge, that can account for the sequence of
transitions from one to two to three skeletal elements along the anteriorposterior
(AP) axis of the limb, their spatial asymmetry, and the temporal sequence of their
development across the AP axis. The anterior half of the limb gives rise to most
of the humerus, the radius and digit II. The posterior half gives rise to the ulna and
digits III and IV. The pattern is laid down in a proximodistal (PD) and posterior
anterior sequence. For example, digit IV appears before digit III, which appears
before digit II.
Bard and Lauder (1974) showed that a Turing model was unrealistic as a mechanism for digit development because of its sensitivity to initial conditions and
parameter values. Newman and Frisch (1979) proposed a Turing model for skeletal development along the limb, but based their analysis on a linearized version of
the model and did not fully address the role of nonlinear effects (Othmer 1986).
.

Pattern in a Chemotactic Model

A

17

Figure 11. Steady-state solutions for n under Type B boundary conditions for varying
s and a. (a) s = 150, a = 1; (b) s = 225, a = 1.25; (c) s = 300, a = 1.43; (d)
s = 400, a = 1.67; (e) s = 550, a = 1.96; (f) s = 700, a = 2.22. [White denotes high
cell density (n > 1.6), black denotes low cell density (n < 1.4)].

More recently, Dillon et al. (1993) showed that different types of boundary
conditions can greatly reduce the sensitivity of patterns produced by a reaction
diffusion system, and, under appropriate boundary conditions, that the model can
generate in a robust and controlled fashion the sequence of transitions from one
to two to three peaks. Their model produces spatially symmetric patterns. Maini
et al. (1992) showed that the introduction of spatial asymmetry in the diffusion
coefficients of the morphogens, through a control chemical, could lead to spatial
asymmetry in the solution. Such spatially asymmetric diffusion is biologically
realistic (Brummer et al, 1991, Maini, 1994). Some of the models discussed
above can account for the spatiotemporal sequence of pattern formation along
the PD axis. At each AP cross section, however, the patterns formed by the
above models arise simultaneously. This is inconsistent with the skeletal pattern
in the limb which occurs in a temporal sequence across the AP axis. In fact,
none of the above models addresses this issue.
In this article we have shown that spatially asymmetric boundary conditions
can give rise to genuine spatio-temporal patterning. If we consider the model
equations (6) and (7) in one dimension with Neumann boundary conditions on n
and Dirichlet boundary conditions on c, then Figs 47 show that, as the domain
length increases, the steady-state solution profile of cell density changes from one
to two to three peaks. The nature of this change depends crucially on the form of
.

18

M. R. Myerscough et al.

the boundary conditions for c. For symmetric Dirichlet boundary conditions, new
peaks arise from the splitting of other peaks, whereas for asymmetric boundary
conditions peaks arise near the boundary where x = 1. The latter is more
consistent with the spatiotemporal patterning observed in the skeletal development
in the vertebrate limb.
It should be noted that the patterns arising from model equations (6) and (7),
with asymmetric boundary conditions appear in a temporal sequence across the
AP axis which is not entirely consistent with that observed in the limb. However,
we can modify the model by explicitly including cell differentiation and assuming
that differentiation not only depends on the density of undifferentiated cells,
but also on an external factor which exists in a gradient across the AP axis.
Specifically, the model then becomes
n
= (Dn) (nc) + s (r n(1 n) M(x)(n))
t
N
= s(n)M(x)
t


n
c
2
= c+s
c
t
n+

(10)
(11)
(12)

where N is the nondimensional density of differentiated cells, (n) is the function

controlling cell differentiation and M(x) is the function which modulates the rate
of differentiation across the AP axis. The function (n) must have the property
that it takes very low values for small n and saturates to a high value for large
n. Specifically, we consider
(n) =

(n n c )m
+ (n n c )m

(13)

where , n c , , and m are suitably chosen parameters. The modulating function

M(x), which models the effect of an external factor which diffuses across the
domain from one boundary, is given by
M(x) =

cosh(mx)
cosh(20m)

(14)

where m is a constant. This form for the function M(x) has been chosen because
it represents the steady-state solution of a reactiondiffusion equation where a
modulating substance has a fixed concentration (taken to be 1) at x = 1 with zero
flux at x = 0. The only kinetic term is a linear degradation term. This modified
model of equations (10)(12) gives a temporal sequence which is consistent with
those observed experimentally in the chick limb (Fig. 12).
The inclusion of cell differentiation alone results in the model having no nontrivial steady state. However, the inclusion of a small amount of dedifferentiation

Pattern in a Chemotactic Model

19

A
1.5

1.0

0.5

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

x
B
1.5

1.0

0.5

0.2

0.4
x

Figure 12. Simulation of the model system (10)(12) in one spatial dimension for
Neumann boundary conditions on n and asymmetric Dirichlet conditions on c, with
c(0) = 0, c(1) = 0.45. Parameter values are D = 0.25, = 2, r = 0.04, = 1, Md =
0.3, M0 = 1, = 0.1, n c = 0.3, m = 8; in (a) s = 400 and in (b) s = 1000. Transient
solutions are shown in (a) for t = 12 000 and in (b) for t = 6000. In each case the
rightmost aggregate in undifferentiated cells, N , dominates the initial transient solution.
The final steady-state patterns (not shown) in (a) have two peaks in N and in (b) have
three peaks in N . In both cases the amplitude decreases from left to right. n, ;
N [the profiles for c in (a) and (b) are not shown but have two and three peaks,
respectively].

20

M. R. Myerscough et al.

leads to steady-state solutions in n, c and N which are qualitatively similar to

those for the case of no dedifferentiation. In these cases, the number of peaks,
their position and their temporal sequence remain unaltered, but away from the
peaks, N eventually tends to zero exponentially.
While pattern is being laid down the limb bud is growing. This can be modeled
by taking the scale factor s to be a function of s(t) in equations (10)(12).
Solving the model in this case leads to results that are qualitatively similar to
those described.
The above arguments are based on one-dimensional simulations, while the limb
is, of course, a three-dimensional structure. To make our model more realistic,
we now consider the two-dimensional simulations presented in Section 6. We
may think of Type B conditions as corresponding to a two-dimensional cross
section of a limb perpendicular to the PD axis, with the source of chemical being
at one end of the AP axis. Type A conditions correspond to a cross section
perpendicular to the dorsoventral axis, with the source once more located on
one side of the AP axis. Figure 11 shows that, as the domain size increases,
there is a transition sequence from one to two to three spots, consistent with the
sequence observed along the PD axis, whereas Fig. 9 shows a sequence similar
to that observed along the plane of the limb parallel to the dorsoventral axis.
The correct sequence in this case will form only if the aspect ratio is sufficiently
small, suggesting that pattern forms in a small strip and not in the full domain.
This is consistent with the idea that pattern is laid down in the progress zone.

8.

D ISCUSSION

In this article we have considered the effect of different boundary conditions

on a cell-chemotaxis pattern generator. The introduction of more general boundary conditions leads to significant differences in the spatiotemporal properties of
the patterns produced compared with those exhibited under standard Neumann
boundary conditions. These are the internalization of the peaks in cell density
(that is, no peaks form on the boundary) and the formation of stable patterns
with asymmetric peaks. We have chosen to consider the particular case in which
the cell density obeys zero flux on the boundary, which is a biologically realistic
assumption. Hence, mixed symmetric or asymmetric boundary conditions are
introduced by imposing Dirichlet conditions on c.
Imposing homogeneous Dirichlet boundary conditions on chemoattractant concentration c while continuing to use Neumann conditions for cell density n results
in peaks occuring within the domain rather than on the boundary. Such boundary
conditions eliminate solutions with peaks on the boundary and therefore enhance
the robustness of the patterns which do occur. Dillon et al. (1993) showed a
similar phenomenon in reactiondiffusion models.
Asymmetric fixed boundary conditions on c produce a spatially asymmetric
pattern. Such boundary conditions can actually enhance the pattern formation
.

Pattern in a Chemotactic Model

21

potential of the cell-chemotactic model. With these boundary conditions peaks

do not form simultaneously, but rather form in a sequence which usually starts
with the first large peaks forming nearer the boundary where c is lower.
Maini et al. (1992) investigated the spatially asymmetric patterns produced
by a modified reactiondiffusion system. In their model, two chemicals (morphogens) interacted to give diffusion-driven instability in the standard fashion.
The introduction of a third chemical which modulated the diffusion coefficients
of the morphogens and had a spatially asymmetric profile resulted in spatially
asymmetric steady-state solutions in morphogen concentrations. In their case
the introduction of spatially varying diffusion coefficients splits the domain into
pattern-forming and nonpattern-forming subdomains. A similar result cannot be
obtained for the chemotactic system because a spatial variation in the diffusion
or chemotactic coefficients results in one part of the domain being more strongly
chemotactic than the other part. The former actually induces gradients in the latter leading to patterns in both parts of the domain. This is an inherent difference
between these models.
In the cell-chemotactic model peaks are set up in a specific temporal sequence. In contrast, in the three chemical reactiondiffusion system discussed in
Maini et al. (1992), the peaks in morphogen concentration form simultaneously.
Therefore, the mechanism by which the asymmetry is set up in these two different systems is intrinsically different. In the reactiondiffusion system, a spatially
varying diffusion coefficient essentially modifies the underlying pattern formed
due to diffusion-driven instability. In the chemotactic system, heterogeneous
fixed-boundary conditions in c create gradients in c which implicitly influence
the patterning mechanism. Therefore, although the final steady states of the models may look qualitatively similar, the processes that lead to the asymmetry are
essentially very different.
.

A CKNOWLEDGEMENTS
This work was supported by an Australian Research Council Small Grant (to
MRM). PKM would like to thank the School of Mathematics and Statistics,
University of Sydney, for their hospitality. KJP was supported by an EPSRC
Earmarked Studentship in Mathematical Biology.

A PPENDIX
A.1. Analytical Determination of Aspect Ratio. We begin our determination of
the aspect ratio through consideration of the more standard case, in which the
boundary conditions for both the cell density and chemical concentration are zero
flux everywhere. A standard linearization of the model system (6)(7) about the

22

M. R. Myerscough et al.

homogeneous steady state (1, cs ), yields the following linearized set of equations:
u
= D 2 u 2 v sr u
t


v

2
uv
= v+s
t
( + 1)2

(15)
(16)

where u = n 1, v = c cs . We seek solutions to the linearized system of the

form:
w exp(ik x + t),

(17)

subject to the boundary conditions, where w = (u, v)T and x = (x, y)T .
The dispersion relation is given by:
2 + a(k 2 ) + b(k 2 ) = 0

(18)

a(k 2 ) = (D + 1)k 2 + s(r + 1)



s
2
4
k 2 + s 2r
b(k ) = Dk + Ds + r s
(1 + )2

(19)

where

(20)

Heterogeneous solutions will develop only when the real part of becomes
positive for admissible values of k 2 . For zero-flux boundary conditions on the
domain [0, 1] [0, W ], such k 2 are of the form

k =
2

q2
p + 2
W

(21)

where p, q N.
We know that a(k 2 ) > 0 for all k 2 > 0, and so heterogeneous patterns are
possible only if there exists a range of k 2 for which b(k 2 ) < 0.
The conditions are thus given by:

s
(1 + )2

s
> (D + r )s
(1 + )2
2
(D + r )s > 4Ds 2r

(22)
(23)

When these conditions are satisfied, we obtain a range of k 2 R+ for which

the steady state is unstable in the linearized system. When a finite domain is

Pattern in a Chemotactic Model

23

considered, however, the only heterogeneous solutions that grow are those that
satisfy the boundary conditions. This equates to requiring the set given by:



( p, q) N :
2

q2
p + 2
W
2


(kl , ku )
2

(24)

to be nonempty, where kl and ku are, respectively, the lower and upper intercepts
of the dispersion curve with the k 2 axis.
The linearized solutions give us an idea to the type of solutions we may expect
to the system for a given set of parameters. We therefore fix all model parameters
apart from W and consider the range of unstable modes when W is small.
2
The function given by Wq 2 tends to infinity as W tends to zero and thus for
sufficiently small W (and q > 0),
2

q2
> ku 2
W2

(25)

In this case, the only unstable modes are those for which q = 0. Obviously,
2

2
q2
2 2
2 q
<

p
+

,
W2
W2

(26)

and so the first mode that becomes unstable as W is increased, for which q 6= 0
will be the (0, 1) mode.
Thus to find the minimum aspect ratio, we determine the size of W for which
the (0, 1) mode becomes unstable. The problem effectively becomes a onedimensional problem. The critical length for stripe solutions to the zero-flux
problem is given by

Wc =
(27)
ku
where ku is given by
ku2



s
1
=
(D + r )s
2D
(1 + )2
v
!!
u 
2
u
s
t
+
(D + r )s 4Ds 2r
(1 + )2

(28)

Using the parameters in the simulations presented in this paper, the critical
length is determined by

5 2
Wc = q
(29)

s(21 + 41)

24

M. R. Myerscough et al.
Table A1. Numerical and analytical comparison of maximal domain width for stripe
solutions. Numerical results are to two decimal place accuracy and given are the smallest
Wc for which stripes break down in two dimensions.

s
300
400
900
3000

Analytical Wc
0.245
0.212
0.141
0.078

Numerical Wc

0.21
0.14
0.08

We have thus shown how the linear analysis can predict the onset of patterning
for the model with zero flux at the boundary.
In our system, however, only the cell population has zero-flux conditions everywhere. The chemical concentration is bounded on two sides of the domain with
asymmetric Dirichlet conditions, and zero flux at the other boundaries. The effect
of increasing s, however, is to increase the domain size and, effectively, increase
the length of domain between the Dirichlet boundaries. We can thus argue that, as
s is increased, the influence of the boundaries on the patterning at the center of the
domain weakens, and so the system can be approximated by the one-dimensional
problem of searching for the point of instability of the (0, 1) mode.
This hypothesis is strongly supported by the results of numerical simulations.
A comparison between critical aspect ratios predicted by the above analysis and
the numerically determined values is shown in Table A1.
The analytical prediction apparently breaks down near s = 300. We have
analyzed aspect ratios down to 0.5, (corresponding to Wc = 2), with no apparent
breakdown of the stripe solutions. This suggests that the boundary conditions are
strong enough to influence the entire domain.
R EFERENCES

Bard, J. and I. Lauder (1974). How well does Turings theory of morphogenesis work?
J. Theor. Biol. 45, 501531.

Brummer, F., G. Zempel, P. Buhle, J. C. Stein and D. F. Hulser (1991). Retinoic acid
modulates gap junction permeability: A comparative study of dye spreading and ionic
coupling in cultured cells. Exp. Cell. Res. 196, 158163.

Chaplain, M. A. J. and A. M. Stuart (1993). A model mechanism for the chemotactic

response of endothelial cells to tumour angiogenesis factor. IMA J. Math. Appl. Med.
Biol. 10, 149168.

Dillon, R., P. K. Maini, and H. G. Othmer (1993). Pattern formation in generalised

Turing systems: I. Steady-state patterns in systems with mixed boundary conditions.
J. Math. Biol. 32, 345-393.

Ford, R. M. and D. A. Lauffenberger (1991). Analysis of chemotactic bacterial distributions in population migration assays using a mathematical model applicable to steep
or shallow attractant gradients. Bull. Math. Biol. 53, 721749.

Grindrod, P., S. Sinha and J. D. Murray (1989). Steady state spatial patterns in a cellchemotaxis model. IMA J. Math. Appl. Med. & Biol. 6, 6979.

Gueron, S. and N. Liron (1989). A model of herd grazing as a travelling wave, chemotaxis and stability. J. Math. Biol. 27, 595606.
.

Pattern in a Chemotactic Model

25

Hofer, T., J. A. Sherratt and P. K. Maini (1995a). Dictyostelium discoideum: Cellular

self-organisation in an excitable medium . Proc. Roy. Soc. Lond. B259, 249-257.

Hofer, T., J. A. Sherratt and P. K. Maini (1995b). Cellular pattern formation during
Dictyostelium aggregation. Physica D85, 425-444.

Hunding, A. and M. Brns (1990). Bifurcation in a spherical reaction diffusion system

with imposed gradient. Physica D44, 285302.

Keller, E. F. and L. A. Segel (1970). Initiation of slime mold aggregation viewed as an

instability. J. Theor. Biol. 26, 399415.

Keller, E. F. and L. A. Segel (1971). Travelling bands of bacteria: a theoretical analysis.

J. Theor. Biol. 30, 235248.
Kulesa, P., G. C. Cruywagen, S. R. Lubkin, P. K. Maini, J. S. Sneyd and J. D. Murray

(1995). Modelling spatial patterning of the primordia in the lower jaw of Alligator
mississippiensis. J. Biol. Syst. 3, 975985.

Lapidus, I. R. and R. Schiller (1976). A model for the chemotactic response of a bacterial
population. Biophys. J. 16, 779789.

Maini, P. K., M. R. Myerscough, K. H. Winters and J. D. Murray (1991). Bifurcating spatially heterogeneous solutions in a chemotaxis model for biological pattern
formation. Bull. Math. Biol. 53, 701719.

Maini, P. K. (1994). Coupled models for spatial organization in development. Ber.

Bunsenges. Phys. Chem. 98, 1172-1175.

Maini, P. K. and M. Solursh (1991). Cellular mechanisms of pattern formation in the

developing limb. Int. Rev. Cytol. 129, 91-133.

Maini, P. K., D. Benson and J. A. Sherratt (1992). Pattern formation in reaction diffusion
models with spatially inhomogeneous diffusion coefficients. IMA J. Math. Appl. Med.
& Biol. 9, 197213.

Murray, J. D. (1989). Mathematical Biology, New York: Springer-Verlag

Murray, J. D., D. C. Deeming and M. J. W. Ferguson (1990). Size dependent pigmentation pattern formation in embryos of Alligator mississipiensis: Time of initiation of
pattern formation mechanism. Proc. Roy. Soc. (London) B239, 279293.

Murray, J. D. and M. R. Myerscough (1991). Pigmentation pattern formation on snakes.

J. Theor. Biol. 149, 339360.

Myerscough, M. R. (1988). A chemotactic model for biological pattern formation. DPhil

thesis, Corpus Christi College, Oxford, U.K.

Myerscough, M. R. and J. D. Murray (1992). Analysis of propagating pattern in a

chemotaxis system. Bull. Math. Biol. 54, 7794.

Nagorcka, B. N. (1989). Wavelike isomorphic prepatterning in development. J. Theor.

Biol. 137, 127162.

Newman, S. A. and H. L. Frisch (1979). Dynamics of skeletal pattern formation in

developing chick limb. Science 205, 662668.

Oster, G. F., J. D. Murray and A. K. Harris (1983). Mechanical aspects of mesenchymal

morphogenesis. J. Embryol. Exp. Morphol. 78, 83-125.

Oster, G. F., N. Shubin, J. D. Murray and P. Alberch (1988). Evolution and morphogenetic rules. The shape of the vertebrate limb in ontogeny and phylogeny. Evolution
45, 862884.

Othmer, H. G. (1986). On the Newman-Frisch model of limb chondrogenesis. J. Theor.

Biol. 121, 505508.

Othmer, H. G. and A. Stevens (1997). Aggregation, blow-up and collapse: the ABCs
of taxis in reinforced random walks. SIAM J. Appl. Math 57, 10441081.
.

26

M. R. Myerscough et al.

Sherratt, J. A., E. H. Sage and J. D. Murray (1993). Chemical control of eucaryotic cell
movement: a new model. J. Theor. Biol. 162, 2340.

Turing, A. (1952). The chemical basis of morphogenesis. Phil. Trans. Roy. Soc. Lond.
B237, 3772.

Winters, K. H. (1985). Entwife User Manual (Release 1), Theoretical Physics Division
AERE, Harwell, Oxfordshire, U.K.

Wolpert, L. (1969). Positional information and the spatial pattern of cellular differentiation. J. Theor. Biol. 25, 147.
.

Received 13 August 1995 and accepted 12 June 1997