Abstract
We present a generalized Keller–Segel model where an arbitrary number of chemical compounds react, some of which are produced by a species, and one of which is a chemoattractant for the species. To investigate the stability of homogeneous stationary states of this generalized model, we consider the eigenvalues of a linearized system. We are able to reduce this infinite dimensional eigenproblem to a parametrized finite dimensional eigenproblem. By matrix theoretic tools, we then provide easily verifiable sufficient conditions for destabilizing the homogeneous stationary states. In particular, one of the sufficient conditions is that the chemotactic feedback is sufficiently strong. Although this mechanism was already known to exist in the original Keller–Segel model, here we show that it is more generally applicable by significantly enlarging the class of models exhibiting this instability phenomenon which may lead to pattern formation.
1. Introduction
Pattern formation is a fascinating area of biology with many unanswered questions. An early example is furnished by morphogenesis, which focuses on the following question: How can an initially spherically symmetric system, such as an embryo in early stages of development, lead to a creature as spherically asymmetric as a person, or an animal with a spotted coat? One of the first attempts to theoretically explain morphogenesis (by means of a mathematical model) was by Alan Turing. In Citation18, he proposed that an instability mechanism in specific systems of chemical reaction–diffusion systems is the underlying cause of morphogenesis by showing that under certain conditions on both the reaction and diffusion terms, some such systems which are stable without diffusion are destabilized with diffusion Citation13. Over time, pattern formation in various biological systems has been attributed to such Turing instabilities in systems of reaction–diffusion equations.
The main goal of this article is to present an alternate avenue potentially leading to pattern formation via chemotaxis. The key idea of this mechanism is already present in the well-known Keller–Segel model Citation11 (reviewed below), as elucidated by Schaaf Citation17. It relies on the destabilization via a chemotactic term of an otherwise stable uniform steady state. This mechanism of destabilization based on the chemotactic sensitivity is also analysed for the case of the two equation system in one spatial dimension in Citation19. In this article, we will identify a larger class of systems in higher spatial dimensions where a similar chemotactic destabilization mechanism may lead to pattern formation.
Chemotaxis occurs when the movement of a species is influenced by chemicals in the environment Citation9. It naturally arises in a wide variety of interesting biological settings, including cancer metastasis, angiogenesis, immune system function and egg fertilization Citation8. Considering reaction–diffusion systems with chemotactic terms, we show that there are quite general classes of chemical reaction networks (CRNs) which can be destabilized by means of an inherent feedback mechanism. The feedback arises as follows. A population of cells produces a chemical, which in turn participates in a CRN. The CRN leads to one or more chemical products. One of these products serves as the chemotactic signal for the cells, thereby completing the feedback loop. Under certain conditions this feedback loop has unstable homogeneous steady states which, in the absence of the feedback, would be stable. The destabilization of homogeneous stationary states, caused by the feedback mechanism built into the system, may give rise to pattern formation. As mentioned earlier, such a mechanism is already present in a simple CRN appearing in the well-known Keller–Segel model Citation11. In this model, the signal emitted by the cells is the same as the chemotactic signal, and it is subject only to decay. Since this CRN consists of a single decay reaction, it is perhaps one of the simplest possible CRNs. In this article, we find a larger class of CRNs with destabilizable homogeneous steady states.
To further explain the difference with the traditional destabilization mechanism, let us first review Turing's diffusion-driven instability for a two-species system (generalizations to multi-species systems are also well known), as explained in Citation14, using a reaction–diffusion system of the form
The question Turing asked is whether (u*, v*), considered as a steady state of EquationEquation (1), can be linearly unstable, even if it is a stable steady state of EquationEquation (2)
. In other words, is it possible that the eigenvalue problem arising from linearizing EquationEquation (1)
at (u*, v*), namely
The model that we propose for pattern formation on the other hand is a fully nonlinear reaction–diffusion system which also contains a nonlinear chemotaxis term. Destabilization of homogeneous steady states in our model does not require any restrictions on the diffusion coefficients, although we do impose some conditions on the sign pattern of the Jacobian corresponding to some of the reaction terms. These may, however, be less crucial than the role played by the chemotactic term.
The Keller–Segel model for chemotaxis is one of the most widely studied models. It was developed focusing on the cellular slime mold Dictyostelium discoideum. As explained in Citation9 Citation11, the derivation of the Keller–Segel model is based on the aggregation stage of D. discoideum’s life cycle. The unicellular slime mould grows by cell division until it depletes its food source, and begins to enter a starvation mode. After some time, one cell will emit a signal of cyclic adenosine monophosphate (cAMP). The other cells are chemotactically attracted to this signal, and in turn begin to emit cAMP themselves. The cells also produce an enzyme which degrades the cAMP by first binding to it and forming a complex which in turn breaks up into the enzyme plus what we call a ‘degraded product’. To model this, the following assumptions are made in Citation9 Citation11:
• | Let u(x, t) denote the cell density of the slime mold, v(x, t) the density of the chemoattractant cAMP, η(x, t) the density of the enzyme, c(x, t) the density of the complex formed by the binding of cAMP and the enzyme, and d(x, t) denote the density of the degraded product. | ||||
• | The chemoattractant is produced by the amoeba at a rate f(v) per amoeba and the enzyme is produced at a rate g(v, η) per amoeba, which, according to [11 is allowed to vary with both the chemoattractant and enzyme concentrations. | ||||
• | The chemoattracant, enzyme, and complex react as shown below. | ||||
• | The cAMP, the enzyme and the complex all diffuse according to standard Fick's law. | ||||
• | The cell concentration changes by random diffusion, as well by chemotaxis, due to the cells ‘climbing the gradient’ of the chemoattractant cAMP. |
The last assumption is again modelled using Fick's law of diffusion, but with the amoebic flux J (u) consisting of a purely diffusive part and a chemotactic contribution, namely
Keller and Segel Citation11 use a steady-state assumption to reduce the above set of four equations to a two-equation system which models only the density of the chemoattractant and the cell density. Further simplifications result in the so-called Citation8 Citation9 ‘minimal system’,
Our work fits in the context of the more recent attempts to investigate multi-species, multi-chemical systems with chemotactic terms Citation3–6,Citation10 Citation12 Citation20. Just as in the case of the Keller–Segel model, the blow-up phenomenon in these generalized models has received considerable attention Citation3 Citation5 Citation10 Citation20. Of particular interest to the issue we address here is the recent work by Fasano et al. Citation6, Horstmann Citation10, and Liu and Wang Citation12. All articles consider the existence of non-homogeneous or periodic steady states, a signature of pattern formation for multi species models. In Citation6, a two-species two-chemical system is proposed where each species produces one of the two chemicals, which in turn serve as attractive chemotactic signals for the species that produce them. In addition, however, these chemicals are also chemotactic repellants for the species that do not produce them. The model is analysed for scalar spatial domains, and focuses on the existence of non-uniform and periodic steady-state patterns. Motivated in part by Fasano et al. Citation6, Horstmann Citation10 proposed an n species, m chemicals model where all chemicals are allowed to serve as chemotactic signals for the species, and the species themselves can chemotactically attract or repel each other. In addition, the chemicals may be produced and/or consumed by the species and, in principle, could react with each other in fairly general ways. Several natural questions that have attracted the attention of several researchers in the context of the Keller–Segel model are addressed in Citation10 for the multi-species, multi-chemical models. They include blow-up, global existence of solutions, existence and non-existence of non-uniform steady states, destabilization of uniform steady states as a way to achieve pattern formation, and the existence of Lyapunov functionals. Although Citation10 begins with this very general model, the focus of the article quickly turns to specific cases, motivated by the literature, where only a limited number of non-reacting chemicals (often just one or two) appear in the model. One notable exception is in [Section 3.1]Citation10, where a one-species, two chemical system is proposed, where the two chemicals react in a non-trivial way. The model is inspired by work of Citation2, where an E. coli population is studied in the presence of two chemotactic attractants, namely glucose and oxygen. In Citation12 a system with two chemical agents is also considered, but in this case one is a chemoattractant while the other acts as a chemorepellent. In contrast, the emphasis of the work presented here is on possible pattern formation in models with several chemicals which react according to fairly general reaction networks, but only a single species that chemotactically responds to one of these chemicals.
In Section 2 of this article, we explain our generalized model and the system of partial differential equations which arises from it. In Section 3, we discuss conditions for the existence of steady-state solutions to the system in a special case. In Section 4 we state our main results, which give sufficient conditions under which steady states are unstable. Section 5 gives examples of theoretical biological systems which meet the conditions of our theorems and also provides a motivating example for generalizing our results.
2. The generalized model
In this section we present a generalization of the Keller–Segel model Equation(5) where the modelled species interacts with several chemical compounds.
As before, consider a population of some species occupying Ω, an open bounded connected n-dimensional domain. Let the population density at a point x in Ω and time t be u(x, t). The population climbs the gradient of a chemical of concentration v
N
. To model the more general and the more likely biological scenario, where this chemical is in reaction with several other compounds in the environment, we introduce the concentrations of N−1 other chemicals, denoted by (). These chemicals can interact with each other as well as with v
N
creating a CRN. Writing
, we model the rate of change of
due to the CRN by
for some function
. In applications, this function
is often determined by mass action kinetics. Note that in general
includes decay of the chemicals.
Some of these compounds can also be produced by the species. To model this scenario, consider a general subset and assume that the species produce the chemical v
l
at a rate
for all l∈S. The remaining compounds are not generated by the species, so we set α
l
=0 for l¬∈S. Let
be the vector whose ith component is the α
i
introduced above. Thus, the rate of change of
would be governed by
, in the absence of any diffusive effects.
However, we also want to account for spatial diffusion. In general, both the species population and the chemical concentrations diffuse. Assume that the diffusion coefficient for u is given by the constant D>0, and the diffusion coefficients for each v i are given by the constants [Dtilde] i >0. Let [Dtilde] be the diagonal matrix whose ith diagonal entry is [Dtilde] i . Assuming that the chemotactic sensitivity function is χ u, for some constant χ>0, we arrive at our generalized model:
As an immediate example of an application of our generalized model, consider the unreduced Keller–Segel model Equation(4) with the chemotactic sensitivity function
and the rates f(v) and g(v, η) assumed to be constant. We change the notation to let v
1 denote the density of the emitted enzyme, v
3 the density of the chemoattractant cAMP and v
2 the density of complex formed from the cAMP and the enzyme. Then setting
In the following sections, we will examine the constant stationary states of the model Equation(6). We will be especially concerned with conditions for the destabilization of such constant states.
3. Existence of homogeneous steady states
A steady state of our generalized model Equation(6)
satisfies
Before we do so, we need to introduce some terminology from matrix theory that we will use throughout. For any square matrix A, let ρ(A) denote the spectral radius of A. We write A≥0 for a real matrix (or vector) A if each entry of A is non-negative. Similarly, A>0 signifies that every entry of A is positive.
Definition 3.1
A Metzler matrix
is a matrix having non-negative off-diagonal entries, i.e.
We will also need M-matrices, which can be defined in many equivalent ways. Following Citation1, we introduce the definition below.
Definition 3.2
A matrix
is an M
-matrix if it satisfies
(1) | |||||
(2) A ii ≥0. | |||||
(3) There exists a non-negative matrix B≥0 and a number
|
Clearly, if A is an M-matrix, then −A is a Metzler matrix. An M-matrix A is invertible if and only if there exists a B≥0 and s>ρ(B) such that A=sI−B. This is because by the Perron–Frobenius Theorem Citation1, the spectral radius of a matrix B≥0 is also an eigenvalue of the matrix. We will also need the following well-known result [p. 137]Citation1:
Lemma 3.3
Suppose
is non-singular. Then A is an M
-matrix if and only if
A
−1≥0.
We are now ready to address the above-mentioned question of the existence of homogeneous steady states, under the assumption that each component of is a linear function.
Proposition 3.4
Suppose
for some A in
and
. If A is a non-singular M
-matrix, then for every positive constant u
*, there exists a unique non-negative homogeneous stationary solution
of Equation
Equation(6)
. Conversely, A is a non-singular M
-matrix if Equation
Equation(6)
has a non-negative homogeneous stationary solution
for every u
*>0 and every
.
Proof
First, suppose that A is a non-singular M-matrix. Let u * be any positive constant. Then setting
Conversely, suppose is a non-negative homogeneous stationary solution obtained using u
*=1 and
, the unit vector in the ith coordinate direction. Then, by EquationEquation (8)
,
. Letting C denote the matrix whose ith column is
, this implies that AC=I, the identity. Hence A is invertible and, moreover, since
, we have A
−1=C≥0. Hence by Lemma 3.3, A is a non-singular M-matrix. ▪
Having given a sufficient condition under which plenty of non-negative homogeneous steady states exist, we now proceed to study the stability of such states, assuming they exist. Note that we no longer assume that is linear in the ensuing analysis.
4. Instability of homogeneous steady states
In order to analyse the stability of a given homogeneous steady state , we begin by linearizing the spatial partial differential operator in EquationEquation (6)
about
. Let λ be an eigenvalue of the linearized operator, i.e. λ satisfies
We will use the weak formulation of EquationEquation (10), where u and the components v
i
are all sought in the (complex) Sobolev space H
1(Ω) consisting of square integrable functions whose first-order distributional derivatives are also square integrable. By definition of the weak eigenproblem, the number λ is an eigenvalue of the linearized operator if there exists a non-trivial
in
satisfying
4.1. Reduction to finite dimensions
The first step in our study of EquationEquation (10) proceeds by generalizing a method of Schaaf Citation17. This allows us to relate the eigenvalues of the infinite dimensional eigenproblem Equation(12)
to a finite dimensional eigenproblem. The latter is easier to analyse. The finite dimensional eigenproblem involves an
matrix, written in block form as
Theorem 4.1
The number λ solves the eigenproblem
Equation(12)
if and only if it is an eigenvalue of M(μ
i
) for some integer
i≥0.
Proof
Suppose λ satisfies EquationEquation (12) with a non-trivial
in
. Then, since
form a complete orthonormal set in L
2(Ω), it is possible to find an i such that not all of the N+1 numbers
To prove the converse, suppose λ satisfies
Remark 1 Recalling that the mean of u in Equation (6a) does not vary with t, we may restrict ourselves while linearizing to perturbations U with . Suppose we also restrict the perturbations in V
k
to be spatially non-homogeneous by imposing
. Then, following the proof of Theorem 4.1, it is easy to see that λ is an eigenvalue of M(μ
i
) for some i≥1 if and only if it solves the eigenproblem Equation(12)
with the additional conditions
for all k=1, …, N.
4.2. A sufficient condition for destabilization
We now focus on conditions under which the linearized operator has at least one eigenvalue λ in the right half of the complex plane. In such a case, the corresponding stationary state is called linearly unstable. In view of Theorem 4.1, we only need to study the eigenvalues of the parametrized matrix M(μ). Let us begin by reviewing a few useful results from matrix theory that we shall need. The following definitions and the next two theorems are well known – see e.g. [Theorems 1.4 and 2.35 in Chapter 2]Citation1.
Definition 4.2
The directed graph G(A) associated to an N×N matrix A, consists of N vertices
where an edge leads from V
i
to V
j
if and only if a
ij
≠0.
Definition 4.3
A directed graph G is strongly connected if for any ordered pair (V i , V j ) of vertices of G (with i≠j), there exists a sequence of edges (a path ) which leads from V i to V j .
Definition 4.4
A matrix A is called irreducible if G(A) is strongly connected.
Theorem 4.5 (Perron–Frobenius)
Let A≥0 be a square matrix.
(1) If A>0, then ρ(A) is a simple eigenvalue of A, greater than the magnitude of any other eigenvalue. | |||||
(2) If A is irreducible, then ρ(A) is a simple eigenvalue, any eigenvalue of A of the same modulus is also simple, A has a positive eigenvector
|
Theorem 4.6 (Spectral radius bounds)
Let A≥0 be an irreducible N×N matrix. Letting s i denote the sum of the elements of the i th row of A, define
With the help of these well-known results, we can now give a sufficient condition for the destabilization.
Theorem 4.7
Assume that the system
Equation(6)
has a positive homogeneous steady-state solution
and let J be as in Equation
Equation(11)
. Suppose J and
satisfy the following conditions:
(1) There is an
i
*
with
| |||||
(2) J is irreducible, and | |||||
(3) J is Metzler. |
Remark 2 The biological interpretation of the first condition is that the organism must produce at least one of the chemicals in the reaction network. The second condition is a technical condition which we will relax in the next theorem. The third condition is the most restrictive of the conditions, but we give an example of a large class of CRNs for which it holds in Corollary 5.2. We also give examples of CRNs which do not meet condition three, but for which a similar result still holds (see Section 5).
Remark 3 For any fixed i
* we may rescale u
* by to rewrite Equation (7a) as
Proof
We use Theorem 4.1. Accordingly, it is enough to prove that M(μ i ) has a positive eigenvalue λ for some μ i . We will now show that for any μ<0, the matrix M(μ) has a positive eigenvalue under the given assumptions.
To this end, fix some μ<0, and set . Let
In fact, B is an irreducible non-negative matrix. To see this, consider the directed graph of J made of vertices . The graph G(B) is obtained after augmenting G(J) by a vertex corresponding to the first row and column, say V
0, and adding vertex loops as needed to account for the positive diagonal of B. Since J is irreducible, G(J) is strongly connected. Moreover, since
, we conclude that there is a path from V
j
to V
0 for any j. Similarly, since K>0, there is a path from V
0 to V
j
for any j. Hence, G(B) is strongly connected and B is irreducible.
We now claim that
In the remainder of this proof, we establish the claim Equation(14). To this end, let us denote by B
t the transpose of the matrix B and note that for any integer m≥1,
To prove EquationEquation (15), we start by observing that the (i, j)th entry of (B
t)
m
is
This proves that the first limit in EquationEquation (14) is infinity. The proof that the second limit is also infinity follows in the same manner by finding, for each 0≤i≤N, a path of length m from V
i
to V
i*. ▪
4.3. Relaxing the irreducibility condition
As in the previous proof, let be the vertices of the G(J). One way to check the second condition (irreducibility) of Theorem 4.7 is to verify that there is a path from any V
i
to V
j
for all i≠j. In this subsection, we will give another condition involving only one path that is often easier to check.
Let us first recall some further terminology. A strongly connected component of a graph is a maximal strongly connected subgraph. Note that a strongly connected graph consists of exactly one strongly connected component. We consider the following equivalence relation on vertices of a graph: Two vertices V i and V j in a graph are equivalent if there is a directed path from V i to V j and a path from V j to V i . Then, the equivalence classes created by this equivalence relation consist exactly of the vertices of the strongly connected components of the graph. Keeping these in mind, recall the following definition of Citation1.
Definition 4.8
The classes of an N×N non-negative matrix A are the disjoint subsets
corresponding to the vertices of the equivalence classes of
G(A). We also identify a class with its strongly connected component.
If A is reducible, then by reordering the vertices, a permutation matrix P can be found so that
Theorem 4.9
Assume that the system
Equation(6)
has a positive homogeneous steady-state solution
and suppose that J and
satisfy the following conditions:
(1) There is an i
*
with
| |||||||||||||||||
(2) J is Metzler. |
Then, if either the product χ u
*
or α
i
*
is sufficiently large,
is linearly unstable.
Proof
As in the proof of Theorem 4.7, fix μ<0, choose r large enough so that B=M(μ)+r I has positive diagonal entries, and consider the graph G(B
t) with vertices .
We proceed by identifying a cycle in the graph. Since , there is an edge from V
0 to V
i
*
. Also, by our assumption, there is a path from V
i
*
to V
N
. Moreover, since
there is an edge from V
N
to V
0, thus completing a cycle. In particular, we have shown that V
0, V
i
*
, and V
N
are part of the same strongly connected component in G(B
t).
Therefore, finding a permutation matrix P as in EquationEquation (17), we conclude that there is a triangular matrix
with a diagonal block T
j
corresponding to the strongly connected component containing V
0, V
i
*
, and V
N
. Due to the strong connectivity, T
j
is irreducible. Now, we can prove that
Thus, for sufficiently large K or α i * , the triangular matrix T has an eigenvalue larger than r, and so does B. Hence, M(μ)=B−rI has a positive eigenvalue. ▪
5. Examples
In this section, we present several examples of how to use the theory developed in the previous sections. We consider the system Equation(6) with a
that describes the kinetics of some CRN linking the chemicals
. In all the examples,
is obtained by the laws of mass action kinetics.
Example
[(The linear case)]5.1 Suppose is linear, i.e.
Let us first observe a consequence of mass action kinetics. The reaction with a rate constant k
ji
gives rise to an off-diagonal entry k
ji
in the matrix A. Since the rate constants are non-negative, this implies that A must be a Metzler matrix.
If, in addition, we know that −A is a non-singular M-matrix, then Proposition 3.4 can be applied to conclude the existence of non-negative steady states.
To investigate stability, we proceed by studying the paths in the CRN and applying our previous results, as exemplified next.
Corollary 5.2
Assume that Equation
Equation(6)
has a positive homogeneous steady-state solution
. Assume that the kinetics of a CRN on the chemicals
is described by a linear function
which is obtained by the law of mass action kinetics. If, for some chemical
v
i
*
, there is a path in the CRN from v
i
*
to v
N
, and α
i
*
>0, then
is linearly unstable whenever χ u
*
or α
i
*
is sufficiently large.
Proof
We only need to verify the conditions of Theorem 4.9. View the CRN as a directed graph on vertices with an edge from V
i
to V
j
if and only if there is a reaction producing the chemical v
j
from the chemical v
i
. This graph, except possibly for diagonal loops, is exactly G(A
t), where A is as in EquationEquation (18)
. The linearity of
also implies that the Jacobian J in EquationEquation (11)
coincides with A. Hence, the given path condition verifies Condition (4.9) of Theorem 4.9. We have already seen above that the assumptions of mass action kinetics imply that J=A is Metzler, thus verifying Condition (4.9) of Theorem 4.9. ▪
Example
[(Dimerization and decay)]5.3 One of the simplest chemical reactions is dimerization or the reaction
Recall from Section 3 that homogeneous steady states are solutions of EquationEquation (8),
Eliminating v 2 from Equation (21a) using Equation (21b), we find that v 1 satisfies
We next turn to determining stability for this example. Clearly is not linear, so we cannot apply Corollary 5.2. However, we can apply Theorem 4.9. The Jacobian of
evaluated at a positive homogeneous steady state
is
Notice that the Metzler matrix J above has two negative eigenvalues (since the trace of J is negative, and the determinant is positive). Therefore, the CRN by itself (without the introduction of chemotaxis) is stable. However, once we introduce chemotaxis, the arguments above show that with a high enough chemotactic sensitivity (χ u *), or a high enough chemical production by the species (α1), the stable system can be destabilized.
Example 5.4 Now we give an example where a hypothesis of Theorems 4.7 and 4.9 fails. Nonetheless, we can analyse the system by applying Theorem 4.1.
Consider a reaction of the form
Then the function describing the reaction network is
To verify that positive steady states exist, we look for u>0 and satisfying
Notice that the third equation of is omitted because it is the same as the second up to the sign. By subtracting Equation (23a) from Equation (23b), we obtain the equivalent system
From these equations it is clear that
Having verified that positive steady states exist, we turn to determine their stability. The Jacobian of the function evaluated at a positive homogeneous steady state
is
However, by virtue of Theorem 4.1, the linearized stability of steady states is determined by the eigenvalues of the matrix
First consider the case K=0, when there is no chemotactic term. Then, M(μ) is block triangular and its eigenvalues are μ D<0 and those of its next diagonal block, which can be expressed as
We next turn to the case where K>0. Using cofactor expansion along the first row, we can compute det M(μ) to be
To conclude, we claim that M(μ) must have a negative eigenvalue when K≥K
0. If all eigenvalues of M(μ) are real, the claim is obvious. If not, complex eigenvalues occur in conjugate pairs, so in order that their product is negative, at least two of the four eigenvalues must be real. Furthermore, these two real eigenvalues must have opposite signs to satisfy .
Therefore, for all sufficiently large values of K, the steady state is linearly unstable.
Acknowledgements
This work was partially supported by the NSF under grants DMS-0818050, DMS-1014817 and DMS-1211635. PDL wishes to thank the VLAC, the Flemish Academic Centre for Science and the Arts, for hosting and supporting him during a sabbatical leave from the University of Florida, the Université Catholique de Louvain-la-Neuve for awarding him a visiting professorship in the fall of 2011, and the University of Florida for a Faculty Enhancement Opportunity (FEO) fund.
References
- Berman , A. and Plemmons , R. 1994 . Nonnegative Matrices in the Mathematical Sciences , Philadelphia : SIAM .
- Boon , J.-P. and Herpigny , B. 1986 . Model for chemotactic bacterial bands . Bull. Math. Biol. , 48 : 1 – 19 .
- Calvez , B. and Perthame , B. 2006 . A Lyapunov function for a two-chemical species version of the chemotaxis model . BIT , 46 : S85 – S97 .
- Dung , L. and Smith , H. L. 1999 . Steady states of models of microbial growth and competition with chemotaxis . J. Math. Anal. Appl. , 229 : 295 – 318 .
- Espejo , E. E. , Stevens , A. and Velázquez , J. J.L. 2010 . A note on non-simultaneous blow-up for a drift-diffusion model . Differ. Integral Equ. , 23 : 451 – 462 .
- Fasano , A. , Mancini , A. and Primicerio , M. 2004 . Equilibrium of two populations subject to chemotaxis . Math. Models Methods Appl. Sci. , 14 : 503 – 533 .
- Gantmacher , F. R. 1959 . Applications of the Theory of Matrices , New York : Interscience Publishers .
- Hillen , T. and Painter , K. J. 2009 . A user's guide to PDE models for chemotaxis . J. Math. Biol. , 58 : 183 – 217 .
- Horstmann , D. 2003 . From 1970 until present: The Keller–Segel model in chemotaxis and its consequences. I . Jahresber. Deutsch. Math.-Verein. , 105 : 103 – 165 .
- Horstmann , D. 2011 . Generalizing the Keller–Segel model: Lyapunov functionals, steady state analysis, and blow-up for multi-species chemotaxis models in the presence and repulsion between competitive interacting species . J. Nonlinear Sci. , 21 : 231 – 270 .
- Keller , E. and Segel , L. 1970 . Initiation of slime mold aggregation viewed as an instability . J. Theoret. Biol. , 26 : 399 – 415 .
- Liu , J. and Wang , Z. 2012 . Classical solutions and steady states of an attraction-repulsion chemotaxis in one dimension . J. Biol. Dyn. , 6 : 31 – 41 .
- Maini , P. K. 2004 . The impact of Turing's work on pattern formation in biology . Math. Today , 40 : 140 – 141 .
- Murray , J. D. 2003 . Mathematical Biology II: Spatial Models and Biomedical Applications , 3 , New York : Springer .
- Perthame , B. 2007 . Transport Equations in Biology , Basel : Birkhauser-Verlag .
- Satnoianu , R. A. , Menzinger , M. and Maini , P. K. 2000 . Turing instabilities in general systems . J. Math. Biol. , 41 : 493 – 512 .
- Schaaf , R. 1985 . Stationary solutions of chemotaxis systems . Trans. Amer. Math. Soc. , 292 : 531 – 556 .
- Turing , A. M. 1952 . The chemical basis of morphogenesis . Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. , 237 : 37 – 72 .
- Wang , X. and Xu , Q. 2012 . Spiky and transition layer steady states of chemotaxis systems via global bifurcation and Helly's compactness theorem . J. Math. Biol. , Available at http://dx.doi.org/10.1007/s00285-012-0533-x
- Wolansky , G. 2002 . Multi-components chemotactic system in absence of conflicts . European J. Appl. Math. , 13 : 641 – 661 .