Digital Object Identifier (DOI): 10.1007/s002850000035
Mathematical Biology
K.J. Painter· P.K. Maini · H.G. Othmer
Development and applications of a model for cellular
response to multiple chemotactic cues
Received: 15 February 1999 / Revised version: 1 February 2000 / Published online: 14 September 2000 – c Springer-Verlag 2000
Abstract. The chemotactic response of a cell population to a single chemical species has been characterized experimentally for many cell types and has been extensively studied from a theoretical standpoint. However, cells frequently have multiple receptor types and can detect and respond chemotactically to more than one chemical. How these signals are integrated within the cell is not known, and we therefore adopt a macroscopic phenomenological approach to this problem. In this paper we derive and analyze chemotactic models based on partial differential (chemotaxis) equations for cell movement in response to multiple chemotactic cues. Our derivation generalizes the approach of Othmer and Stevens [29], who have recently developed a modeling framework for studying different chemotactic responses to a single chemical species. The importance of such a generalization is illustrated by the effect of multiple chemical cues on the chemotactic sensitivity and the spatial pattern of cell densities in several examples. We demonstrate that the model can generate the complex patterns observed on the skin of certain animal species and we indicate how the chemotactic response can be viewed as a form of positional indicator.
1. Introduction
A characteristic feature of a wide variety of living organisms is that they can alter their motile behavior in response to external environmental cues and, in particular, to chemical stimuli. Examples range over many orders of magnitude in spatial scale. The green turtle, Chelonia mydas travels over 1000 miles to its breeding ground, through detection of an unknown chemical source originating there [19] and species of shark have the ability to detect traces of blood miles from an injured animal. At the other extreme, fibroblast cells are thought to move into wounds in response to chemical signals emitted at the wound site, and in cellular slime mold
K.J. Painter∗, P.K. Maini: Centre for Mathematical Biology, Mathematical Institute, 24–29 St Giles’, Oxford, OX1 3LB, UK
H.G. Othmer∗∗: Department of Mathematics, 270A Vincent Hall, University of Minnesota, Minneapolis, MN 55455, USA. e-mail:othmer@math.umn.edu
∗Supported by an EPSRC earmarked studentship in mathematical biology. Present address: Department of Mathematics, University of Minnesota.
∗∗Supported in part by grant GM29123 from the National Institutes of Health. Key words: Cell movement – Chemotaxis – Multiple signals – Pattern formation
Dictyostelium discoideum, amoebae respond to the chemoattractant cAMP over
spatial scales of the order of a cell diameter [28].
The directed motile response of organisms to chemical cues is called
chemo-taxis. Frequently this individual response is manifested at the population level by
the collective movement of a population up (or down) a gradient in concentration of a chemoattractant (or chemorepellent). Organisms also rely on other forms of taxis for guidance, including light (phototaxis), electric fields (galvanotaxis), mag-netic fields (magnetotaxis), and adhesive gradients (haptotaxis). Smaller organisms such as bacteria frequently rely on some form of kinesis, which involves modula-tion of either the speed or the duramodula-tion of movement in response to environmental conditions. Both taxes and kineses can often be described by similar macroscopic equations, called chemotaxis equations (to be discussed later), even though the stochastic processes that describe movement at the single-organism level may be quite different. Thus one frequently lumps macroscopic descriptions of both taxes and kineses together under the general heading of taxes.
Chemotaxis in bacterial species such as E. coli and S. typhimurium has been extensively studied both experimentally and theoretically. E. coli cells move by ro-tating rigid helical flagella in a corkscrew-like manner [32]. Each cell contains 6–8 flagella distributed uniformly over the cell surface [37], which, when rotated coun-terclockwise, coalesce into a propulsive bundle, resulting in a relatively straight “run” [5]. When rotated clockwise they fly apart, resulting in a “tumble” which reorients the cell but causes no significant translocation. The cell thus alternates between runs and reorienting tumbles. In the absence of stimuli, the probability per unit time of a tumble is essentially independent of when the last tumble occurred [37]. In a chemoeffector gradient, the cell carries out chemotactic migration by extending runs that carry it in favorable directions. Using specific chemoreceptors to monitor its chemical environment, E. coli perceives spatial gradients as tempo-ral changes in attractant or repellent concentration. The cell in effect compares its environment during the past second with the previous 3–4 seconds and responds accordingly. Attractant increases and repellent decreases transiently raise the prob-ability of CCW rotation, or ‘bias’, and then a sensory adaptation process returns the bias to baseline, enabling the cell to detect and respond to further concentra-tion changes [36]. Attractants include sugars, amino acids and small peptides, and negative chemotaxis (repulsion) occurs in response to noxious substances.
The chemotactic response of cells to a single attractant or repellant has been characterized experimentally for many cell types and has been extensively stud-ied from a theoretical standpoint [30, 28]. However, many cell types have multiple receptor types and can detect and respond chemotactically to more than one chem-ical. For example, in wiring the nervous system of a developing embryo, tight directional control of the growth of axons is needed in order to ensure that they connect properly with nerve cells, muscles and other tissues. The four principal cues involved are short range (surface) attractants and repellents and long range (diffusive) attractants and repellents. Many diffusible attractants and repellents have been discovered, principally Netrin-1 [33] and members of the semaphorin families [12], [22]. While it functions as an attractant for certain types of axons, Netrin-1 also serves as a repellent for other types [7]. It is believed that many chemicals are
involved in guidance, and that each type of molecule attracts some axons and repels others.
Chemotaxis has also been widely-studied in the clinical sciences because it is an important component of certain aspects of the immune response. Chemat-tractants released at the site of an infection attract white blood cells in the first response [1], and it has recently been found that adaptation plays an important role in guiding leukocytes through the attractant landscape set up by spatially-separated sources [15]. In tumor-induced angiogenesis, cancerous cells secrete chemoattrac-tants into the surrounding tissue, stimulating endothelial cells in neighboring blood vessels to migrate towards the tumor, thus creating a blood vessel link to the tumor. Several chemoattractants (or tumor angiogenic factors) have been identified, in-cluding basic fibroblast growth factor, angiogenin and vascular endothelial growth factor.
Despite the evidence showing that cells respond chemotactically to multiple environmental cues, mathematical models and experimental studies typically focus on the response to a single stimulus. However, the presence of multiple chemotactic signals raises a number of important questions, viz.,
• How are the multiple signals integrated intracellularly and interpreted to initiate movement? For example, how is the occupancy of different receptor types in a bacterium integrated and transduced into a response? It is known that the chemotactic protein CheYpmodulates the rotational bias of a flagellar motor, and thus signal integration is easy for pathways using the same chemotaxis proteins, but how is it done when different receptor types generate different signal types?
• What effects do competing attractant and/or repellent signals have on macro-scopic spatial patterns? For example, bacteria form a variety of intricate spatial patterns under different nutrient levels [4, 38, 3, 6] but little is known as to to whether these patterns depend significantly on the presence of multiple sig-nals, and whether simple changes in the balance between different signals can account for observed changes in the pattern types.
• How is information filtered and processed by the organism? If signals of differ-ent magnitudes (e.g. in gradidiffer-ent strength or concdiffer-entration level), or of differdiffer-ent spatial and temporal scales are detected, how can an organism process the in-formation so as to to generate the optimal response? This question may be particularly relevant in ecology, where multiple signaling cues can vary on different spatial scales.
These questions have not been answered at the molecular or higher level in any system, and therefore we adopt a more phenomenological approach here. We begin with a master equation for a stochastic space jump process and reduce this to a partial differential equation to describe the movement of organisms in response to multiple tactic signals. We derive a number of forms for the chemotactic sensitivity (defined later), based on how information external to the organism is detected and transduced to result in movement, and we investigate in detail spatial patterning in a simple model for a cell population responding to two competing attractants.
Applications of the work are considered in bacterial cell movement and the com-plicated pigmentation patterns seen on the skins of certain animals.
How these signals are integrated within the cell is not known, and we therefore adopt a macroscopic phenomenological approach to this problem. In this paper we derive and analyze chemotactic models based on partial differential (chemotaxis) equations for cell movement in response to multiple chemotactic cues.
2. Modeling the tactic response to multiple signals
Othmer and Stevens [29] investigated continuum limits of a continuous time, dis-crete space random walk described by a master equation, in which walkers execute instantaneous jumps in space at random times. In their analysis, which was moti-vated by results of Davis [8] on reinforced random walks, the transition probabilities depend on a control substance that itself evolves in time and space. In the following we begin with the same master equation, but now under the assumption that the transition rate is dependent on multiple signals.
Suppose that the probability that a walker is atn ∈ Z at time t, given that the walk begin at 0 at t = 0, is given bypn(t). By hypothesis this satisfies the evolution equation
∂pn
∂t = T+n−1(U)pn−1+ T−n+1(U)pn+1− (T+n(U) + T−n(U))pn. (1) HereT±n are the conditional probabilities per unit time of a 1–step jump ton±1, and(T+n(U) + T−n(U))−1is the mean waiting time at thenth site. U represents the density of the control species, and is given by
U = ... ... ··· ... u1,n−1/2u2,n−1/2· · · uk,n−1/2 u1,n u2,n · · · uk,n u1,n+1/2u2,n+1/2· · · uk,n+1/2 ... ... ··· ... , (2)
fork control species. We have defined the control species on a lattice of half step size, i.e.uj,n−1/2is the density of speciesujat positionn−1/2. The three different types of models we consider here are characterized by the dependence of the transition rates onU. Following the classification of [29], we define these as strictly local,
barrier and gradient models. 2.1. Strictly Local Models
In this model the probability of a jump depends only on the level of control species at the site of the walker. We denote by unthe densities of control speciesu1, . . . , uk
at siten. We assume that the domain is homogeneous and isotropic, and therefore the transition probabilities do not depend directly on the lattice site, but only on local densities un. It follows thatT−= T+(≡ T), and consequently (1) becomes
∂pn
We definex = nh, extend the function p to R, and expand the right hand side to second order in h as a function of x. We introduce a scaling of the transi-tion rates,T(u) = λT(u), with the property that the limit limh→ 0,λ→∞λh2 = positive constant(≡ D) exists, then in this diffusion limit,
∂p ∂t = D
∂2
∂x2(T(u)p), (4)
The existence of this limit depends on the existence of a priori bounds on higher-order derivatives. In higher space dimensions the second derivative is replaced by the Laplacian. Equation (4) is subject to initial conditions and, when a finite domain is considered, boundary conditions. Furthermorep should always be nonnegative.
In one dimension the particle flux is
J = −D∂(T(u)p)∂x = −D(T(u))∂p∂x − Dp∂(T(u))∂x . (5)
We define the chemotactic sensitivity to speciesi, χi, by
χi = −DTui (6)
and then can write the flux as
J = −D(Ti(u))∂p∂x+ p k X i=1
χi(u)∂u∂xi. (7)
Ifχi > 0 (i.e.,Tui < 0), species i is an attractant and cells move up the gradient ofi at a speed that may depend on other species via the dependence of χi on u. However the global movement is determined by the sum of the responses to the individual gradients.
2.2. Barrier models
In the barrier model the transition rates for a cell at siten depend on u at n ± 1/2. One finds that without further conditions on the transition rates there is no directed movement in response to the signal; the chemotactic sensitivity vanishes identically [29]. However, if the transition rates are renormalized so that the mean waiting time is constant across the domain, then aggregation may occur. Thus we define the renormalized transition ratesN±as
N+n(U) + N−n(U) = constant, (8) and without loss of generality we take the constant as 1. The above equation is satisfied by choosing the following forms for our transition rates,
N±
n(U) = N±(un−1/2, un+1/2) = T(u T(un±1/2)
whereN± : Rk× Rk → R. We again assume that the lattice is homogeneous, we rescale the renormalized transition rates asN = λN, and we assume that
D ≡ limn→∞,h→0λh22 exists. Then in the diffusion limit we obtain
∂p ∂t = D ∂ ∂x ∂p ∂x− 2p k X j=1 N1,j− N2,j ∂uj ∂x . (10)
HereN1,jrepresents the derivative with respect to thejthcomponent of the first
argument andN2,j represents the derivative with respect to thejthcomponent
of the second argument, both evaluated at(un, un). These are determined from Equation (9), and Equation (10) becomes
∂p ∂t = D ∂2p ∂x2 − D ∂ ∂x pXk j=1 (ln T)uj · ∂uj ∂x . (11)
Thus the chemotactic sensitivity with respect to speciesi for the renormalized barrier model is given by
χi = D(ln T)ui. (12)
2.3. Non-local or gradient models
The non-local model is based on assumption that the transition rates depend on differences between lattice sites of the density of the control species. In the sim-plest case one supposes only nearest neighbor differences, and then can define T∓n±1:U → R by
T+n−1= α + β(τ(un) − τ(un−1)), (13) T−n+1= α + β(τ(un) − τ(un+1)), (14) andτ : Rk → R with α(≥ 0) and β constants. This gives the master equation
∂pn
∂t = α(pn+1− 2pn+ pn−1) − β(τ(un+1) − τ(un))(pn+1+ pn)
−β(τ(un) − τ(un−1))(pn+ pn−1). (15) We assume a scaling in the transition rates such thatD = limh→0,λ→∞λh2, and we obtain in the diffusion limit
∂p ∂t = Dα ∂2p ∂x2 − 2Dβ ∂ ∂x p k X i=1 ∂τ ∂uk · ∂uk ∂x ! . (16)
Thus the individual tactic sensitivities are given by
χi = 2Dβτui (17)
3. The chemotactic sensitivities for multiple signals
It is clear from the preceding derivations that the dependence of the chemotactic sensitivities on the control species arises from the dependence of the transition rates onU. In this section we derive various chemotactic sensitivity functions and show in particular that the presence of multiple signals carries with it the possibility of competition and nonlinear interactions in any of the detection, transduction and response steps.
3.1. Linear
The simple example of linear dependence of the transition rate on u, given by T = T0+ b · u where T0and b are such thatT > 0, illustrates key differences
in behavior between the three models. It follows from (6), (12) and (17) that the chemotactic sensitivities for the three models are
χj = −Dbj, χj = DT bj
0+ b · u, χj = 2Dβbj.
(18) In the barrier model the chemotactic sensitivity is monotone decreasing in the total chemical concentration if allbis are positive, whereas the other models have constant sensitivity with respect to all species. The sign of the sensitivities imply that the local model results in aggregation at low concentrations of speciesj if bj > 0, whereas barrier and gradient models lead to aggregation at high concentrations under the same condition [29].
3.2. Receptor-based response laws
The linear law used above is based on the assumption that the transition rate contains a basal ratea and a component proportional to the concentration. However signals are often detected via receptors or other detectors on the surface, and since their number is limited, the transition rate cannot increase indefinitely. Moreover, at the cellular level receptors frequently span the membrane and are inactive when free of ligand, but active when bound. In the absence of a detailed transduction scheme we assume that when there is a single control species the response is proportional to the fraction of receptors that is occupied by the control species [29]. Thus if the binding reaction is written
Ri+ u k +
* )
k− Ra, (19)
then the fraction of activated receptors is given by Ra = R0u/(K + u), where
K = k−/k+ andR
0 = Ra+ Ri is the total number of cell surface receptors. If
the transition rate comprises a basal rate and a component proportional to the total number of occupied receptors, then for the local and renormalized barrier models we have that
T = T0+ T1 R0u
and this leads to chemotactic sensitivities of the form
χ = −D T1R0K
(K + u)2
for strictly local sensing, and
χ = D T1KR0
(K + u)[T0(K + u) + T1R0u]
for the renormalized barrier model.
Similarly, if we assume thatτ = R0u/(K +u) for the gradient model we obtain
the sensitivity
χ = χ0
(K + u)2 (21)
whereχ0= 2DR0Kβ.
It should be noted that now the sensitivity is monotone decreasing in all cases, and does so quadratically at largeu, which reflects the fact that the detection ma-chinery saturates at large signal strengths.
3.2.1. Different receptor types
When there are two distinct signals present there are two different types of sensi-tivities, depending on whether the signals bind competitively to the same receptor or whether there is a different receptor for each signal. Experimental evidence sug-gests that both cases occur depending on the specificity of the receptors. Some receptors detect only one signal, whereas in other cases the cell surface contains several types of chemotactic receptors, each receptive to a small group of signaling chemicals (see [31] and [20]).
In the latter case, two types of chemoattractant molecules,u and v, bind to two different receptor types,RuandRv. This is a straightforward extension of the one species case, and we derive the following expressions for the number of activated receptors,Rau andRav: Rau = R0uu Ku+ u, Rav = R0vv Kv+ v. (22)
R0u andR0v are the total number of receptors of type Ru andRv, respectively.
Ku = k1−/k+1 andKv= k−2/k2+, wherek1−, k+1, k2−andk2+are the rate constants
due to rate equations similar to that given by (19). It is straightforward to extend the above ton chemical species, and we obtain the following expression for the total number of activated receptors:
X Rauj = n X j=1 R0ujuj Kj+ uj, (23)
whereRauj is the number of activated receptors of type Ruj,R0uj is the total
the transition rate for the gradient model we find that the individual chemotactic sensitivities are given by
χi =
2DβR0uiKi
(Ki+ ui)2. (24)
It is clear from this that these sensitivities have the same form as the sensitivity in the presence of a single chemical signal, which is a result of the fact that there is no interaction of the signals at the detection stage.
3.2.2. Same receptor type
However, if the signals bind competitively then we expect to see this interaction of the signals reflected in the chemotactic sensitivities. Now we suppose that an inactive receptorRiis activated by binding to either of two molecules,u or v. Thus the binding reactions are
Ri+ u k1+ * ) k1− Rau, Ri+ v k2+ * ) k2− Rav. (25)
and we assume that the total number of receptorsR0= Rau+ Rav+ Riis constant.
Assuming binding of either of the molecules occurs on the same time scale and is faster than the subsequent transduction and movement response, we have
R0− (Rau+ Rav) uk1+= k−1Rau, (26) R0− (Rau+ Rav) vk2+= k−2Rav. (27)
If the transition rate is proportional to the number of activated receptors,Rau+Rav, we find from (26) and (27) that
Rau+ Rav =
R0(k+1k2−u + k+2k−1v)
k1−k−2 + (k+1k−2u + k2+k−1v)
. (28)
More generally, for the case ofn species, u1, . . . , un, we defineRaui to be the
number of receptors activated by chemical speciesuiandR0to be the total number
of receptors. Following the above procedure we derive n X i=1 Raui = R0 Pn i=1uiKp/Ki Kp+Pni=1uiKp/Ki, (29) whereKj = k−j/kj+andKp=Qnj=1Kj.
With the above expression for the transition rate, we find that for a gradient model the chemotactic sensitivity for thejthspecies is given by
χj = 2DβKpKj Kp+Pni=1uiKp/Ki2
The only variation in the sensitivity for each individual species results from the different values ofKj.
It follows that the total chemotactic flux is given by 2Dβp
Kp+Pni=1uiKp/Ki2 X
j
KpKj∇uj, (31)
and the direction of the flux is determined by a linear function of the individual chemical gradients.
3.2.3. Nonlinear interaction sensitivity laws
Thus far the transition rates are simply assumed to be a linear combination of the total number of activated receptors, which means that are a number of independent signaling pathways whose outputs are additive. However, there may also be parallel pathways, both of which must be activated to induce movement. Thus we assume a dependence based on the product of the activated fraction of receptors, which gives
T = R0,uu
(Ku+ u) R0,vv
(Kv+ v). (32)
The chemotactic sensitivities for the gradient model are therefore given by
χu= (K χ0,uv
u+ u)2(Kv+ v), χv=
χ0,vu
(Kv+ v)2(Ku+ u). (33) whereχ0,u≡ 2DβKuR0,uR0,vandχ0,v ≡ 2DβKvR0,vR0,u. Notice that the above
implies that for lowv the chemotactic sensitivity response to species u is much weaker than the response at high concentrations ofv. Thus, by reducing the level ofv the chemotactic response of the cell may be drastically reduced despite the presence of gradients ofu. This introduces another level of control, if for example, the production of the signaling chemicals depends on the cell density.
4. Spatial response to multiple cues
The patterns of the spatial distribution of a cell population responding to multiple cues is illustrated in the following theoretical model. We consider a Turing system [39] as a mechanism for providing spatially heterogeneous chemical distributions to which a cell population chemotactically responds. More formally, a cell population density,n(x, t), is defined on the spatial domain ∈ Rn. Chemicals u and v, which are called morphogens, evolve according to a system of reaction-diffusion equations. The cell population responds chemotactically to both chemical species, leading to the equations
∂n
∂t = ∇ · {Dn∇n − nχ1(u, v)∇u − nχ2(u, v)∇v}, (34)
∂u
∂t = ∇(Du· ∇u) + f (u, v), (35)
∂v
whereχ1 and χ2 are the chemotactic sensitivity functions and f and
g
definethe chemical kinetics. Herein cell and chemical diffusion coefficients,Dn, Duand
Dv, are taken as constant. In this first model we assume that the evolution of the morphogens is independent of the cell density. We will refer to the above equations as the hybrid model.
We first consider a uniformly distributed cell population and, without loss of generality, setn(x, 0) = 1. Initial conditions for the chemical concentrations are taken as small random perturbations about the homogeneous steady state. Boundary conditions are taken to be zero flux. The application of the boundary conditions with the zero population kinetics imposes a conservation law on total cell population,
viz.
Z
n(x, t)dx = Z
n(x, 0)dx. (37)
We assume that the equations (34)–(36) have a single positive homogeneous steady state(n0, u0, v0), where without loss of generality we take n0= 1. It is also assumed
that the structure off and
g
is such that (35) and (36) have a stable nonconstant steady state solution (a Turing pattern) for suitable choices of the kinetic parameters, and that the uniform steady state is unstable to small perturbations.4.1. Patterning in one dimension
We first consider the one-dimensional domain, = [0, L], and suppose that the kinetics are described by a simplified model for the glycolysis reaction [27, 9]. The model (34)–(36) is then given by
∂n ∂t = ∂ ∂x Dn∂n∂x − nχ1(u)∂u ∂x − nχv(v) ∂v ∂x , (38) ∂u ∂t = Du ∂2u ∂x2 + δ − κu − uv 2, (39) ∂v ∂t = Dv ∂2v ∂x2+ κu + uv 2− v, (40)
with zero-flux boundary conditions:
Dn∂x∂n− nχ1(u)∂u ∂x− nχv(v) ∂v ∂x = 0 (41) ∂u ∂x = ∂v ∂x = 0 (42)
atx = 0 and x = L. Thus all three components satisfy a homogeneous Neumann boundary condition.
The absence of cellular feedback results in a decoupling of the cell move-ment equation from the reaction-diffusion system. Thus, temporarily we ignore the cell equation (38) and focus on evolution of pattern in the reaction-diffusion model. A linear analysis of Equations (39)–(40) about the homogeneous steady state,(u0, v0) = δ/(κ + δ2), δ
and the parameter region wherein Turing instabilities can develop. It does not pro-vide a specific analytical prediction of the pattern at the heterogeneous steady state because, as the solution grows, nonlinear terms become dominant and the linear analysis is no longer valid. An analytical approximation for the time independent solutions to Equations (39)–(40) can be obtained, and the results are summarized in the Appendix. The general approach can be found in a number of texts [13, 16]. To first order in the deviation from the uniform steady state, the spatial distribution of the morphogens is given by
u(x) = u0+ u1cosjx, (43)
v(x) = v0+ v1cosjx, (44)
whereu1 andv1can be determined explicitly for a particular set of parameters
through the asymptotic expansion. In Figure 1 (a) we compare analytical and nu-merical approximations for the concentrations ofu and v, using the parameters of the figure caption. The validity of solutions derived from a bifurcation analysis, given by Equations (43)–(44), is demonstrated by the very close correspondence with the numerically derived concentrations.
When the chemotactic sensitivities are non-zero it is possible that patterning of the cell population via the response to the nonuniform distribution of the mor-phogens, given (43)–(44). The steady state cell distribution that satisfies the bound-ary conditions is given by
n(x) = K exp Z x 0 1 Dn χ1du dx + χ2dv dx dx , (45)
where the constantK can be evaluated using the conservation equation (37), once theu and v distributions are known.
One of the primary characteristics of a spatial pattern of the cell population (or chemical species) is the number of maxima and minima. For chemical patterns that arise via a Turing instability (with zero flux boundary conditions and close to the bifurcation point), this is relatively straightforward: the pattern possesses a common wavelength with equally spaced peaks. Moreover, the spatial location of turning points is the same for both chemical species – they are either in phase (for pure activator-inhibitor kinetics) or in antiphase (for cross activator-inhibitor kinetics) [9]. The glycolysis reaction is an example of a cross activator-inhibitor reaction system: peaks ofu occur at troughs of v and vice versa.
4.1.1. Patterning in response to a single cue
The effect of taxis on cell patterning can be evaluated by a comparison between the critical points of the cell density and the critical points in the chemical concen-trations. The response to one chemical is considered by takingχ2= 0 in Equation
(45) and assuming thatχ1is a function ofu only. We determine
dn dx = K Dnexp 1 Dn Z χ1(u) du dxdx χ1 du dx. (46)
Thus, providingχ1is non-zero and continuous for allu(x) > 0, the turning points in
the cell density coincide with the turning points in the chemical distribution, since
K exp 1
Dn
R
χ1(u)dudxdx
is positive. All standard response laws, i.e. constant, logistic and receptor, satisfy these conditions, although one can construct a biphasic form that does not to fall into this category.
We suppose that the chemicalu binds to two different cell surface receptors R andS. The fraction of occupied receptors of each type will be given by R0u/(k1+u)
andS0u/(k2+ u) respectively. We further assume that activation of R receptors
induces an attractive response, while activation ofS receptors induces a repulsive response. Then if the transition rate is a linear combination of the activated receptors we have that T(u) = T0+ T1R0u k1+ u− T2S0u k2+ u. (47) We setα = T1R0and choosek1= 1 and T2S0= 1 for convenience. One finds that
the chemotactic sensitivity (using a gradient-based transition rule) is
χ1(u) = αk2 2− 1 + 2u(αk2− 1) + u2(α − 1) (1 + u)2(k 2+ u)2. (48)
If we setα < 1 and αk2> 1, then in addition to turning points that coincide with those of the chemical concentration, there additional critical points correspond-ing to locations at whichχ1changes sign. This biphasic law models chemotactic
attraction at lower chemical concentrations, and chemotactic repulsion at higher concentrations. Consequently, cells aggregate at an intermediate level of concen-tration. In response to the guidance cue BDNF (brain-derived neurotrophic factor), nerve growth cones have been shown to display both chemoattraction and chemore-pulsion, depending on the level of cAMP [35].
4.1.2. Cell patterning in response to both chemicals
If we specify that the chemotactic sensitivities are constantsχ1= α and χ2= β,
then the solution for cell density at a heterogeneous steady state can be obtained analytically, and one finds that
n(x) = Kcexp αu + βv Dn , (49)
whereKcis a constant. Thus the spatial distribution of the cell density is determined by the linear combinationαu + βv of the morphogens. In Figure 1 (b) – (d) we plot the functionαu + βv for constant β and increasing α, using the analytically- and numerically-computedu and v shown in Figure 1(a). If α is sufficiently different from−βv1/u1, there is little difference between the numerically- and
analytically-derived values ofαu + βv, as seen in Figure 1(b). However, as α is increased the structure of the analytically- and numerically-computed profiles diverges, as in (c), and for a small range ofα near −βv1/u1the numerical solution displays greater
structure, with two peaks in the density (cf. (d)). Although the patterns show very different structure, the maximum difference between the two profiles for allα, β is always of the order2, where a definition for can be found in the appendix. This arises from the fact that the difference between the numerical and anlytical approximations of the chemical concentrations is of order2. Asα is increased above−βv1/u1the two profiles again display the same structure.
These results suggest that the analytical approximation forn(x) is only valid outside a neighborhood of the lineα = −βv1/u1. This restriction can be
under-stood by examining the second order terms in the asymptotic expansion for the chemical concentration done in the Appendix. The second order approximations
Fig. 1. (a) Comparison between analytical and numerically derived chemical concentrations. u numerical (solid) and analytical (dashed), v numerical (dot-dash) and analytical (dotted). The discrepancy is very small, reflected by these lines overlying one another. (b) – (d) Comparison betweenαu + βv using the analytical and numerical solutions for u and v plotted in (a). This function reflects the form of the cell density at steady state under constant sensitivity rules. (b)α/Dn= β/Dn= 1.0, (c) α/Dn= 2.8 andβ/Dn= 1.0, (d) α/Dn= 2.87 andβ/Dn= 1.0. For the majority of α/Dn− β/Dnspace the two predictions match. For example, in (b) the solutions lie on each other and are difficult to discern. Asα/Dnis increased towards−v1β/Dnu1, the two forms diverge, (see (c)) until they have different structure, (see (d)). Simulations useδ = 1.2, κ = 0.06, L = π, Du= 1.0 and Dv= 0.09809. The analytical
of the chemical concentrations are given by functions of the form
u(x) = u0+ u1 cos jx + 2(u20+ u21cosjx + u22cos 2jx) (50)
v(x) = v0+ v1 cos jx + 2(v20+ v21cosjx + v22cos 2jx) . (51)
From the expansion in the Appendix, we only know the magnitudes ofu0,u1,u20,
u22,v0,v1,v20andv22. The bifurcation analysis does not determine the magnitudes
ofu21andv21, but does determine the ratio, which is given byu21/v21= u1/v1.
For the above forms, the critical points in cell density are the solutions to the equation
j (αu1+ βv1+ (αu21+ βv21)) sin jx + 2j2(αu22+ βv22) sin 2jx = 0,
(52) which, in turn, has solutions
sinjx = 0, (53)
4(αu22+ βv22) cos jx = −(αu1+ βv1+ (αu21+ βv21)). (54)
The former defines the same primary turning points as for the first order approx-imations, which correspond with those of the chemical concentrations. The latter only has solutions if αu1+ βv1 = O(), and therefore we choose chemotactic
coefficients such thatαu1+ βv1= 0. From the analysis in the appendix we have
αu21+ βv21= 0, and thus if αu22+ βv226= 0, Equation (54) predicts additional
critical points that are out of phase with the primary critical points. For the param-eters considered here, this assumption above is valid. If we fixβ > 0 and increase
α from below the line −v1β/u1, the cell density profiles change nature from a
form equivalent to the chemical profiles to one with a greater number of turning points, and the latter arise as the result of the superposition of cosjx and cos 2jx waves. Hereafter we label such cell patterns as complex to distinguish them from the simple patterns in which the critical points coincide with those of the chemical concentration.
The validity of the above analysis is strengthened by a comparison of analytical and numerical solutions for the special case,α = −v1β/u1. We use this solution
to obtain an analytical determination of the amplitude of the cell density for the cos 2jx wave type pattern, as all cos jx terms drop out (at least up to second order),
n(x) = K exp αu0+ βv0+ 2αu 20 + 2(αu 22+ βv22) cos 2jx Dn ! . (55)
When this is evaluated for the parameters of Figure 1 with α/Dn = 5.5 and
β/Dn = 1.92, we obtain an amplitude of 0.00046, which compares favorably with the amplitude of the pattern derived via numerical simulation, 0.00048. The discrepancy isO(3), where a definition for can be found in the Appendix. Thus cells may exhibit density patterns of greater spatial complexity than that of the un-derlying chemical concentration profiles in response to multiple chemotactic cues.
Next we compare these results with those obtained using sensitivities of the form χ1= α (k1+ u)2, χ 2= β (k2+ v)2, (56) whereα, β, k1andk2are constants andk1, k2> 0. In this case, the form of the cell
density at the heterogeneous steady state can be found as
n(x) = Krexp − 1 Dn α k1+ u(x)+ β k2+ v(x) . (57)
whereKr is a constant. A calculation of the locations of non-trivial turning points gives the following solutions,
αu1(v0+ k2+ v1y)2+ βv1(u0+ k1+ u1y)2= 0, (58)
wherey = cos jx. Rearranging this expression gives the following solutions for y,
y∓= ±(v0+ k2) q −αu1 βv1 − u0+ k1 u1∓ v1 q −αu1 βv1 . (59)
Of course the only real solutions lie in the range [−1, 1], and it is easy to verify that either there is no solution or one solution that satisfies this criterion, depending on the parameters. When one solution does lie in this range we predict that complex patterns develop. These predictions are confirmed by a comparison of numerical and analytical results. Examples of complex patterns are shown in Figure 2. While similar to those observed for constant laws, these patterns have greater amplitude and occur for a larger range of parameters.
Fig. 2. Patterns predicted by the model with receptor response laws given by (56). We use the analytical predictions determined in the Appendix for the parameters given in Figure 1 (u0 = 0.8, v0 = 1.2, u1 = 0.0545, v1 = −0.1562). Figures (a) through (c) demonstrate increasing complexity as we move through α-β space. Notice that the amplitude of the complex patterns is larger than corresponding patterns (Figure 1) that arise for constant sensitivity laws
4.2. Patterning in two dimensions
In two dimensions, we assume chemical concentrations of the form,
u(x, y) = u0+ u1cos nπx Lx cos mπy Ly , (60) v(x, y) = v0+ v1cos nπx Lx cos mπy Ly . (61)
For the above forms two elementary pattern types exist: Ifn and m are non-zero, we have a checkerboard type pattern composed of peaks and troughs, namely spots, whereas if only one ofn or m is non-zero then the pattern consists of stripes.
Using the receptor law sensitivities, it is straightforward to analytically demon-strate that the one-dimensional results carry over to two dimensions. The types of complex patterns formed in the cell density can be classified as (1) rings, (2) interspersing large and small spots and (3), interspersed broad and narrow stripes. Figure 3 illustrates these pattern types.
We now consider more complicated chemotactic sensitivity laws. Extending the derivation of the biphasic sensitivity law earlier to two chemical species, we assume a transitional response, T(u) = T0+T1R1u k1+ u− T2R2u k2+ u+ T3R3v k3+ v − T4R4u k4+ u. (62)
For this transitional response, the chemotactic sensitivities termsχ1 andχ2take
the form given by Equation 48. Analytical results under these laws give rise to more complicated patterns than for the “standard” laws considered earlier. In two dimensions, we observe patterns such as the spots inside rings or variant stripes, as shown in Figure 4.
Fig. 3. Typical examples of the complex patterns observed in two dimensions. (a) Rings (b) Interspersed large and small spots and (c) Alternating broad and narrow stripes. Patterns obtained using receptor type laws for underlying chemical patterns such thatu0 = 0.8, v0= 1.2, u1= 0.2 and v1 = −0.6. We take β/Dn= 1.0, k1= 0.8, k2= 1.0 and α/Dn= (a) 2.5, (b) 1.8, (c) 1.8. A threshold has been chosen for clarity of presentation
Fig. 4. Typical examples of patterns obtained under more complicated sensativity responses. (a) Spots inside rings and (b) Series of stripes. In both simulations, we use the parameters T1R1 = 1.0, T2R2 = 3.85, T3R3 = 1.0, T4R4 = 3.85, k1 = 0.5, k2 = 2.0, k3 = 0.5, k4 = 1.95, u0= v0= 1.0, u1= +0.5, v1= −0.3. Once again we have chosen a threshold in the cell density to demonstrate the pattern more clearly. For (a) we use an underlying spotted chemical pattern, whereas (b) uses an underlying striped pattern
4.3. Robustness of patterns
As we have shown, when the chemotactic sensitvity depends on two chemicals each of the different types of sensitivity laws examined leads to increased complexity of the spatial distribution of the cells. For constant chemotactic sensitivities this complexity is dependent on the fine structure of the chemical distributions (which analytically corresponds to the higher order terms of the asymptotic expansions), whereas with receptor laws it arises through nonlinear terms in the sensitivity func-tion. This fundamental difference may lead us to expect differences with respect to robustness of the patterns, and hence the applicability of different mechanisms to biological pattern formation. Robustness here refers to pattern characteristics such as amplitude of the patterns and the parameter region in which they occur. For constant sensitivity laws the amplitude of the patterns is very small, and they occur in a very restricted range of parameter space. However receptor laws produce larger amplitude pattern and these patterns exist in a larger region of parameter space.
These results are significant when the underlying chemical concentrations are less regular. For kinetic parameters that lie in a small neighborhood of the point at which the Turing instability occurs in the chemical system the pattern is very regular since the power spectrum of the solution is dominated by the amplitude of the unstable mode. However, the pattern may be much less regular far from the primary bifucation points.
In Figure 5 (a)–(c) we display the chemical concentrations derived by numerical simulation of the two dimensional reaction-diffusion system. The patterns illustrate the two pattern types commonly observed in Turing systems: A hexagonal array ((a) and (b)) and stripes (c). Increasing the domain size in (a)–(b) results in less regular
Fig. 5. Cell patterning under chemical concentrations obtained by numerical simulation of the Turing system. Left column shows heterogeneous steady state chemical patterns (v species) (a) Hexagons, δ8qwe = 2.8, κ = 0.06, Du = 1.0, Dv = 0.0125, domain
[0, π]×[0, π]. (b) For the larger domain [0, 2π]×[0, 2π], the hexagonal patterns show less regularity. (c) Stripe solutions for parameter valuesδ = 1.2, κ = 0.06, Du= 1.0, Dv= 0.08,
domain [0, 5π]×[0, 5π]. Cell patterning under constant and receptor sensitivities are shown in the middle and right columns respectively. With the regular chemical concentrations (a), both constant and receptor laws can demonstrate ring type patterns, however those obtained for receptor law have a larger diameter and greater amplitude. Under less regular patterns such as the hexagons of (b) and the stripes (c), receptor laws can give clearly defined rings or thick and thin stripes, whereas no coherent pattern emerges under the constant laws patterning. The middle and right columns show the complex patterns produced in the cell density for constant and receptor sensitivity laws, respectively. Constant laws (d) produce rings of small circumference and amplitude. Receptor laws (g), however, produce large rings of large amplitudes. Furthermore, the region ofα − β parameter space over which these patterns are observed is much greater in the latter case than in the former. This aspect is critical for the production of such patterns on a large domain. A comparison of the last two columns shows that distinct rings of large amplitude (h) or alternating broad and narrow stripes (i) are obtained with receptor laws, whereas no coherent pattern occurs for the constant sensitivity laws, (e) and (f). Thus non-constant sensitivities lead to more accurate tracking of the underlying morphogen distribution.
5. Applications
5.1. Animal skin patterns
A theoretical application of the Turing model is that of Murray ([23], [24], [25]), (see also [2] and [40]), who considered a reaction-diffusion system as the underlying mechanism for pigment patterns on mammalian coats. A reaction-diffusion sys-tem establishes an underlying nonhomogeneous morphogen field which provides the information for pigment cells to differentiate into a particular type depending on the local morphogen concentration. This theory of pigmentation is attractive due to the similarity between patterns generated by the model and coat markings, but even a cursory glance through an encyclopedia of animal patterns yields exceptions. These include alternating thick and thin stripes, which most notably occur in species of fish (e.g. the lionfish) or reptiles (for example, the tail of the gila monster), rings, which are commonly observed on species of wild cat such as the jaguar, leopard and ocelot. More complex patterns include those shown by the thirteen lined ground squirrel which display stripes interspersed with a line of spots.
5.1.1. Animal markings through the hybrid model
We can apply the complex patterns exhibited by the hybrid model in the context of animal coat markings. Essentially, we propose a mechanism whereby pigment cells organize themselves spatially via a chemotactic movement response to gradients of chemicals in the reaction-diffusion system. Further differentiation into a specific pigment type may then arise through sensing a threshold of local cell density. Bio-logical evidence suggests that the mechanism under which pigment cell precursors migrate to dermal layers of the skin may be chemotactic ([10], [18], [21], [17]) and it is not an unreasonable suppositon to assume organisation of cells in the skin also occurs through chemotaxis.
In Figure 6, we compare patterns observed with the hybrid model with patterns seen on animal coats. Ring patterns obtained by numerical simulation of the model are compared with those on the jaguar, and broad and narrow stripes predicted by the model are compared with the lionfish. Theoretically, a third type of common complex pattern was obtained in two-dimensions; the large and small spots shown in Figure 3(b). Both rings and broad and narrow stripes can be observed in large regions of parameter space (under receptor sensitivity laws), and are relatively insensitive to the regularity of the underlying chemical kinetics – see section on robustness above. The large and small spots, however show much greater sensitivity to the underlying pattern. Despite the appearance of such patterns in theory or under highly regular chemical concentrations (such as in Figure 3(b)), they are less likely to arise in biological systems, where greater variation is to be expected. Intriguingly, this may be reflected in pigmentation markings. Both rings and broad/narrow stripes have been observed on the skins of a variety of animals, yet we have been unable to find biological examples of the third type.
Fig. 6. Comparison between complicated animal skin patterns and patterns produced by the hybrid model. The right column compares thick and thin dark stripes with the pig-mentation patterns seen on the lionfish. The left column compares ring patterns with those of the young jaguar. The lionfish picture appears by courtesy of Steve Hogan (http://www.diversionoz.com/Steve/)
5.2. Bacterial chemotaxis
We now consider a simple model for pattern formation in bacteria (for exam-ple, E. coli) in an environment containing two chemotactic substances. A detailed model for signal transduction of a single attractant at the individual cell level has been analyzed [36], but here we use a simpler continuum description that uses the receptor-based schemes developed earlier. This model illustrates how differ-ences in the local movement rules with two species leads to significantly different macroscopic patterning. In one dimension, the cell density evolves according to the equation, ∂n ∂t = Dn ∂2n ∂x2 − ∂ ∂x nχ1(u, v)∂u ∂x −∂x∂ nχ2(u, v)∂v ∂x (63)
Experimentally, techniques have been developed which allow the creation of sta-tionary linear chemical gradients in chemotaxis chambers, [14]. Consequently we shall assume stationary chemical profiles of this type in the above model. The chemotactic sensitivity functions considered are derived using transitional rates of the types considered in Section 3, summarized here in Table 1. We choose
α = β = 10.0, Dn = 1.0, k1 = k2 = k3 = 1.0. We assume an initially
Table 1. Classification and functional form of the chemotactic sensitivities for the simulations
Sensitivity type χ1 χ2 Classsification
Same receptor α/(k1+ k2u + k3v)2 β/(k1+ k2u + k3v)2 Type 1 Different receptor (linear) α/(k1+ u)2 β/(k2+ v)2 Type 2 Different receptor (nonlinear) αv/(k1+ u)2(k2+ v) βu/(k
2+ v)2(k1+ u) Type 3
time-independent version of this equation can be solved analytically to obtain the heterogeneous steady state.
5.2.1. Competing attractant gradients
Chemoattractant gradients ofu and v are plotted in Figure 7 (left). Here we as-sume competing linear attractant concentrations of equal slope and maxima. The cell density at the heterogeneous steady state under the sensitivities of Table 1 in response to these attractant gradients are shown in Figure 7 (right). When both attractants are assumed to bind to the same surface receptor (Type 1), the chemo-tactic fluxes balance, resulting in a homogeneous cell density (solid line, Figure 7 (c)). For attractants binding to different cell surface receptors(Type 2 and Type 3), however, cell density peaks emerge at a point between the maxima of the two chemoattractants (Figure 7 (a) = linear, Figure 7 (b) = nonlinear). The saturation of response at higher chemical concentrations results in a weaker chemotactic effect at these regions, and cells aggregate in the centre.
By varying the maximum concentration (and hence slope) of v, we observe further differences in the macroscopic cell density for the various sensitivity types. In Figure 8, (a), we plotu (solid) and v (dashed) concentrations and now consider varying the maximum ofv, Vm. Cell densities under Type 1 sensitivities are plot-ted in Figure 8 (b). WithVm < 1.0, the gradient of u is larger and u dominates cell movement. Consequently, cell aggregations occur at maxima ofu. As Vm is increased above 1.0, however, the gradient inv is stronger, resulting in switching the cell aggregation tov maximum. A similar effect is observed for Type 2
sensi-Fig. 7. Cell density response to competing attractant signals. Left: Chemical concentrations ofu, solid line, and v, dashed line. Right: heterogeneous cell densities in response to these signals for (a) Type 2 sensitivities, (b) Type 3 sensitivities and (c) Type 1 sensitivities
Table 2. Table summarising the results for the bacterial cell model. Results classified in terms of the position of the aggregation with respect to the high concentrations ofu or v. Note that “atu” may mean the the aggregation is biased towards the higher concentrations ofu, rather than at the maximum. Where the aggregation response is particularly small, we have indicated this. Summaries in bold face demonstrate where the behaviour is particularly different for one sensitivity law, as compared with the others
Equiv. attractants Comp. attract: variousVm Gradient vs constant
Sensitivity Aggregation Agg.Vm Agg.Vmlow Agg.Vmhigh Agg.Vmlow
Type high
Type 1 None Agg. atv Agg. at u Small agg. atu Agg. at u Type 2 Centre Agg. atv Agg. at u Agg. atu Agg. atu Type 3 Centre Agg. atu Small agg. Agg. atu Small
atv agg. atu
tivities, Figure 8 (c), yet in a region aroundVm = 1.0, the cell density maximum lies between the peaks in chemical concentration. A completely different behavior occurs under Type 3 sensitivity functions, Figure 8 (d). Low values ofVm(e.g., 0.1) decrease the effective response both tou and v. Consequently, a small aggregation occurs which is biased towards the maximum concentration ofv. This may initially appear counterintuitive, however from Table 1 we see that for Type 2 sensitivities andv small, χ1is small. A reverse effect is observed at larger values ofVm, e.g. 5.0,
and the aggregation is shifted towards theu maximum. Here, a shallow gradient can be detected by the cells, despite the presence of a much stronger second signal. 5.2.2. Attractant gradient vs constant attractant level
We study the response to a single attractant gradient, yet with a uniform concen-tration of a second attractant. The forms for the concenconcen-trations ofu (solid) and v (dashed) are plotted in Figure 9 (a), and once again we examine responses to the three sensitivity laws under variation ofVm. At high concentration levels ofv and Type 1 sensitivities, cell aggregations at theu maximum is minimal, Figure 9 (b). This occurs via high levels ofv resulting in occupacy of the receptors which satu-rates the response. DecreasingVmresults in an enhanced aggregation. For Type 2 sensitivities, the aggregation tov is the same, regardless of the level of Vm, Figure 9 (c). The opposite behaviour to that described for Type 1 sensitivities occurs under Type 3 sensitivity rules. Here, the aggregation in response tou is minimised at very low levels ofVm, and maximised at high levels, Figure 9 (d).
5.2.3. Summary of the bacteria response
The aggregation in response to two gradients demonstrates highly different be-haviour, depending on the choice of sensitivity functions. These differences are summarised in Table 2. Experimental techniques which allow the creation of linear attractant gradients in chemotaxis gradients may be used to test the effect of two attractant cues, and subsequently test which sensitivity functions may apply in real
Fig. 8. Cell density response to competing attractant signals. (a) Chemical concentrations foru, solid line, and v dashed line, maximum concentration of v determined by Vm. (a) Cell densities under Type 1 sensitivities, (c) cell densities under Type 2 sensitivities, (d) cell densities under Type 3 sensitivities
systems. Type 3 sensitivities assume that the transduction from external signal to internal cell movement occurs via some nonlinear combination of the two attrac-tant cues. The response tov for these sensitivity laws demonstrates how cells can respond to very low levels of a specific signal, despite the presence of stronger signals from other sources.
6. Discussion
In this paper we developed a model for the response of a cell population to two distinct chemicals, either of which may be an attractant or repellant, and we an-alyzed some of the properties of this model. We generalized the Othmer-Stevens approach [29] by modeling the tactic response to multiple signals and considered, in turn, strictly local, barrier, and gradient models. We also considered several dif-ferent types of chemotactic sensitivity laws, and demonstrated that when organisms can respond to more than one chemical cue, the integration of signals via differ-ent combinations of receptor based-response laws leads to many new behavioral responses. We have shown that one manifestation of this more complex response is more complicated spatial patterns.
For instance, the model can exhibit spatial pattern of thick and thin stripes as was shown earlier in a model that involves a single chemotactic substance [30], and
Fig. 9. Cell density response to a single gradient, while another attractant is fixed at a uniform level. (a) Chemical concentrations foru, solid line, and v dashed line, concentration of v determined byVm. (b) Cell densities under Type 1 sensitivities, (c) cell densities under Type 2 sensitivities, (d) cell densities under Type 3 sensitivities
we suggest that this may be applicable in the context of animal coat markings. The application to bacterial chemotaxis shows how signals may be filtered such that very low levels of one stimulus can be detected in the presence of large gradients in a second stimulus. We have also shown that certain types of sensitivity laws favour robust patterning.
The patterns determined by the mechanism presented here derive from the local interpretation by cells of positional information provided by the nonuniform morphogen distribution. Because cells respond by moving, this leads to global patterning and in essence, a non-local interpretation of the positional information. This is in contrast to other models that do not involve chemotaxis, in which the complex patterns are obtained via the superposition of two linearly unstable modes [34, 26] and the response is strictly local. These models can also generate complex patterns, but the one presented here is simpler: we rely on a single pattern generator (the reaction-diffusion mechanism) to create a spatially heterogeneous field and it is the subsequent interpretation of the field that forms the complex patterning. An important question that modeling alone cannot answer is what mechanisms are used in a particular context.
7. Appendix
We obtain analytical approximations for time independent solutions to the model below using standard bifurcation analysis (see, for example, Fife (1979) or Grindrod (1993)). On the scaled one-dimensional domain [0, π] the model is given by:
∂u ∂t = ∂2u ∂x2 + γ 2δ − κu − uv2, (64) ∂v ∂t = d ∂2v ∂x2 + γ 2κu + uv2− v, (65)
with given initial conditions and zero flux at the boundaries. The parameter γ corresponds to a scaling of both spatial and temporal scales.
The homogeneous steady state solution for (64)–(65) is given by (u0, v0) =
δ κ+δ2, δ . We define f(u, v) = γ2 δ − κu − uv2 κu + uv2− v
. The dispersion relation, obtained by looking for solutions of the form exp(λt +ikx) for the system linearized around the steady state(u0, v0), is given by
(λ(k2))2+ a(k2)λ(k2) + b(k2) = 0, (66) where, a(k2) = (d + 1) k2+ γ2κ + v2 0+ 1 − 2u0v0 , (67) b(k2) = dk4+ γ2k21− 2u 0v0+ d κ + v2 0 + γ4κ + v2 0 . (68) Examination of the dispersion relation, Equation (66), demonstrates that eigenval-ues have positive real part forb < 0, where b is determined by Equation (68). Algebraic manipulation of (68) givesb = 0 when:
d = γ 2(2u 0v0− 1) k2− γ4 κ + v02 k2 k2+ γ2 κ + v2 0 , (69)
which is positive fork2 > γ 2(κ+v20)
2u0v0−1. We taked as our bifurcation parameter. We
letk = j (∈ Z) and define d∗to be the solution of Equation (69). Assume, asd is decreased, that stability of the homogeneous solution is lost to thejthmode first. Settingd = d∗− 2, for small, we consider solutions (u, v) of the form:
u v = u0 v0 +X∞ n=1 n un vn . (70)
and expand f as a Taylor series in powers of. LetL be the differential operator:
Lw = ∂ 2 ∂x2 − γ2 κ + v02 −2γ2u 0v0 γ2 κ + v2 0 d∗ ∂2 ∂x2 + γ2(2u0v0− 1) ! w (71)
for suitable functions w : R2→ R. We further define the matrix Lj by the action ofL on terms of the form cos jx.
AtO(): L u1 v1 = 0 and u1x, v1x = 0 at x = 0, π. (72)
This has a solution,
u1 v1 = A u∗ v∗ cosjx, (73)
where(u∗, v∗)T is an eigenvector parallel to u∗ v∗ = 2γ2u0v0 −j2− γ2 v2 0+ κ , (74)
and A is a constant to be determined.
The solution toO(2) is determined uniquely up to the constant B, u2 v2 =− MA2 γ2(v2 0+κ) 1 0 +B u∗ v∗ cosjx +MA 2 Q −4j2d∗− γ2 4j2 cos 2jx. (75) where M =γ22 u0v∗2+ 2u∗v∗v0 , (76) Q = 16d∗j4+ 4j2γ2 d∗ v2 0+ κ + 1 − 2u0v0 + γ4 κ + v2 0 (77) AtO(3): L u3 v3 = γ22u 0v1v2+ 2v0(u1v2+ u2v1) + u1v12 1 −1 + v1xx 0 1 = 9. (78) The Fredholm Alternative (cf. [11]) states that the above equation has a solution if and only if9 is orthogonal to solutions of the adjoint equation. This gives,
0= ν (µ − ν) γ2 Z π 0 v1xxR cos jx dx + Z π 0 2u0v1v2+ 2v0(u1v2+ u2v1) + u1v21 × R cos jx dx. (79)
where R is a constant and(µ, ν) is an eigenfunction of Lt corresponding toj2. This is parallel to γ2 v2 0+ κ j2+ γ2 κ + v2 0 .
We can evaluate the above to give: A A2− v∗νj2 2γ2C(µ − ν) = 0 (80) where C = " 2u0v∗Mj2 Q + 2v0u∗Mj2 Q − v0v∗M γ2 v2 0+κ −2v0v∗M 4Q 4j2d∗+ γ2 +3u∗v∗2 8 # .
Equation (80) implies that asd decreases through the critical value, d∗, two non-trivial equilibria bifurcate from(u0, v0) according to:
u v = u0 v0 ± v∗νj2 2γ2C(µ − ν) 1/2 u∗ v∗ cosjx + O(2). (81) References
[1] Alberts, B., Bray, D., Lewis, J., Raff, M., Roberts, K., Watson, J.D.: Molecular Biology of the Cell. Garland Publishing Inc., New York, London (1994)
[2] Bard, J.B.L.: A model for generating aspects of zebra and other mammalian coat patterns. J. Theor. Biol. 93, 363–385 (1981)
[3] Ben-Jacob, E., Cohen, I., Czirok, A., Vicsek, T., Gutnick., D.L.: Chemomodulation of cellular movement, collective formation of vortices by swarming bacteria, and colonial development. Physica A 238(1), 181– (1997)
[4] Berg, H., Budrene., E.: Complex patterns formed by motile cells in E. coli. Nature 349, 630–633 (1991)
[5] Block, S.M., Segall, J.E., Berg., H.C.: Impulse responses in bacterial chemotaxis. Cell 31, 215–226 November (1982)
[6] Michael P. Brenner, Leonid S. Levitov, Elena O. Budrene: Physical mechanisms for chemotactic pattern formation by bacteria. Biophysical Journal 74 (4), 1677–1693 (1998)
[7] Colamarino, S., Tessier-Lavigne, M.: The axonal chemoattractant netrin-1 is also a chemorepellent for trochlear motor axons. Cell 81, 621–629 (1995)
[8] Davis, B.: Reinforced random walks. Prob. Thy. Rel. Fields pages 203–229 (1990) [9] Dillon, R., Maini, P.K., Othmer., H.G.: Pattern formation in generalised turing systems
i. steady-state patterns in systems with mixed boundary conditions. J. Math. Biol. 32, 345–393 (1994)
[10] Le Douarin, N.M.: Differentiation of Normal and Neoplastic Hematopoietic Cells., chapter Ontogeny of hematopoietic organs studied in avian interspecific chimeras., pages 5–31. Cold Spring Harbor Laboratory, New York (1978)
[11] Dunford, N., Schwartz, J.T.: Linear Operators: Part I. Interscience Publishers, New York (1958)
[12] Fan, J., Raper, J.A.: Localised collapsin cues can steer growth cones without inducing their further collapse. Neuron 14, 263–274 (1995)
[13] Fife, P.C.: Mathematical Aspects of Reacting and Diffusing Systems, volume 28 of Lect. Notes in Biomath. Springer, Berlin, Heidelburg, New York (1979)
[14] Fisher, P.R., Merkl, R., Gerisch, G.: Quantitative analysis of cell motility and chemo-taxis in Dictyostelium discoideum by using an image processing system and a novel chemotaxis chamber providing stationary chemical gradients. J. Cell. Biol. 108, 973– 984 (1989)
