A user’s guide to PDE models for chemotaxis

TLDR: This paper explores in detail a number of variations of the original Keller–Segel model of chemotaxis from a biological perspective, contrast their patterning properties, summarise key results on their analytical properties and classify their solution form.

Abstract: Mathematical modelling of chemotaxis (the movement of biological cells or organisms in response to chemical gradients) has developed into a large and diverse discipline, whose aspects include its mechanistic basis, the modelling of specific systems and the mathematical behaviour of the underlying equations. The Keller-Segel model of chemotaxis (Keller and Segel in J Theor Biol 26:399–415, 1970; 30:225–234, 1971) has provided a cornerstone for much of this work, its success being a consequence of its intuitive simplicity, analytical tractability and capacity to replicate key behaviour of chemotactic populations. One such property, the ability to display “auto-aggregation”, has led to its prominence as a mechanism for self-organisation of biological systems. This phenomenon has been shown to lead to finite-time blow-up under certain formulations of the model, and a large body of work has been devoted to determining when blow-up occurs or whether globally existing solutions exist. In this paper, we explore in detail a number of variations of the original Keller–Segel model. We review their formulation from a biological perspective, contrast their patterning properties, summarise key results on their analytical properties and classify their solution form. We conclude with a brief discussion and expand on some of the outstanding issues revealed as a result of this work.

Figures

Fig. 4 a Numerical simulation of the 2D model for chemotaxis incorporating receptor binding (with α = 0.5), volume-filling (γ = 10), non-linear diffusion (n = 1.0) and saturating chemical production (φ = 1.0). b As for (a), but with the addition of logistic cell growth (r = 0.1). For both sets of numerics we let D = 0.1, χ = 5.0 on a domain of size 20 × 20 with initial data u(x, y, 0) = 1 and v(x, y, 0) = 1.0+ r(x, y), where r(x, y) is a 1% random spatial perturbation of the steady state

Fig. 5 Numerical simulation of the volume filling model (M4a) in 3D. Left hand plots show cell density distribution at various “slices” through the 3D data, right hand plots shows the contours for which u = 5. a Formation of a spherical aggregation of cells from biased initial conditions, u(x, y, z, 0) = 1 and v(x, y, z, 0) = 1+ 0.1 exp(−10((x − 1)2 + (y − 1)2+ (z − 1)2)) on the 2× 2× 2-cube. b Pattern formed from random initial data on a larger domain at t = 100.0. In both simulations we set D = 0.1, χ = 5.0 and γ = 10 and employ 41× 41× 41 spatial grid points

Fig. 1 (M1) Numerical simulation of minimal model showing evolution of cell density (solid line) and chemical concentrations (dash-dot line) to the steady state. (M2)–(M8) Numerical results for the regularised models, showing steady state cell distributions (solid line) at different values of the regularisation parameter. For comparison, the steady state distribution of the minimal model is plotted as the dash-dotted line. The cell distributions marked with an asterisk have been classified (numerically) as plateau type—all others are spikes. For all numerics, the same model set-up is considered: parameters D and χ are set at 0.1 and 5.0, respectively; initial conditions are set to be u(x, 0) = 1 and v(x, 0) = 1 + 0.1 exp(−10x2); 201 discretisation points are employed on a domain of length 1

Table 1 Summary of the submodels and their functional forms. ∗For the nonlocal model the gradient, ∇v is replaced by the non-local gradient ◦∇ρ v (see Sect. 2.4)

Table 3 Critical parameter values for the parameters D = 0.1 and χ = 5 on the interval [0, 1] with zero-flux boundary conditions

Fig. 2 Numerical simulations of models (M1)–(M8) on a larger domain from unbiased initial data. For all numerical simulations the following conditions are employed: parameters D and χ are set at 0.1 and 2.0, respectively; initial data u(x, 0) = 1 and v(x, 0) = 1.0 + r(x) where r(x) is a 1% random spatial perturbation of the steady state; domain [0, 20] with 401 grid parameters

Full text

J. Math. Biol. (2009) 58:183–217
DOI 10.1007/s00285-008-0201-3
Mathematical Biology
A user’s guide to PDE models for chemotaxis
T. Hillen · K. J. Painter
Received: 6 June 2007 / Revised: 18 February 2008 / Published online: 15 July 2008
© Springer-Verlag 2008
Abstract Mathematical modelling of chemotaxis (the movement of biological cells
or organisms in response to chemical gradients) has developed into a large and diverse
discipline, whose aspects include its mechanistic basis, the modelling of specific sys-
tems and the mathematical behaviour of the underlying equations. The Keller-Segel
model of chemotaxis (Keller and Segel in J Theor Biol 26:399–415, 1970; 30:225–
234, 1971) has provided a cornerstone for much of this work, its success being a
consequence of its intuitive simplicity, analytical tractability and capacity to replicate
key behaviour of chemotactic populations. One such property, the ability to display
“auto-aggregation”, has led to its prominence as a mechanism for self-organisation of
biological systems. This phenomenon has been shown to lead to finite-time blow-up
under certain formulations of the model, and a large body of work has been devoted
to determining when blow-up occurs or whether globally existing solutions exist. In
this paper, we explore in detail a number of variations of the original Keller–Segel
model. We review their formulation from a biological perspective, contrast their pat-
terning properties, summarise key results on their analytical properties and classify
their solution form. We conclude with a brief discussion and expand on some of the
outstanding issues revealed as a result of this work.
Mathematics Subject Classification (2000) 92C17
T. Hillen was supported by NSERC and K. J. Painter was partially supported by NIH Grant CA113004.
T. Hillen (
B
)
Department of Mathematical and Statistical Sciences,
University of Alberta, Edmonton T6G2G1, Canada
e-mail: thillen@ualberta.ca
K. J. Painter
Department of Mathematics and Maxwell Institute for Mathematical Sciences,
Heriot-Watt University, Edinburgh EH14 4AS, UK
e-mail: painter@ma.hw.ac.uk
123

184 T. Hillen, K. J. Painter
1 Introduction
From microscopic bacteria through to the l argest mammals, the survival of many
organisms is dependent on their ability to navigate within a complex environment
through the detection, integration and processing of a variety of internal and exter-
nal signals. This movement is crucial for many aspects of behaviour, including the
location of food sources, avoidance of predators and attracting mates. The ability to
migrate in response to external signals is shared by many cell populations, playing
a fundamental role co-ordinating cell migration during organogenesis in embryonic
development and tissue homeostasis in the adult. An acquired ability of cancer cells to
migrate is believed to be a critical transitional step in the path to tumour malignancy.
The directed movement of cells and organisms in response to chemical gradients,
chemotaxis, has attracted significant interest due to its critical role in a wide range
of biological phenomena (the book of Eisenbach [28] provides a detailed biological
comparison between chemotactic mechanisms across different cells and organisms).
In multicellular organisms, chemotaxis of cell populations plays a crucial role through-
out the life cycle: sperm cells are attracted to chemical substances released from the
outer coating of the egg [33]; during embryonic development it plays a role in organ-
ising cell positioning, for example during gastrulation (see [26]) and patterning of the
nervous system [85]; in the adult, it directs immune cell migration to sites of inflam-
mation [111] and fibroblasts into wounded regions to initiate healing. These same
mechanisms are utilised during cancer growth, allowing tumour cells to invade the
surrounding environment [19] or stimulate new blood vessel growth (angiogenesis)
[56]. Extensive research has been conducted into the mechanistic and signalling pro-
cesses regulating chemotaxis in bacteria, particularly in E. coli [4], and in the life cycle
of cell slime molds such as Dictyostelium discoideum [26]. While the biochemical and
physiological bases are less well understood, chemotaxis also plays a crucial role in
the navigation of multicellular organisms. The nematode worm C. elegans undergoes
chemotaxis in response to a variety of external signals [67] while in insects, the fruit
fly Drosophila melanogaster navigates up gradients of attractive odours during food
location [10] and male moths follow pheromone gradients released by the female
during mate location [47].
Theoretical and mathematical modelling of chemotaxis dates to the pioneering
works of Patlak in the 1950s [86] and Keller and Segel in the 1970s [44,45]. The
review article by Horstmann [40] provides a detailed introduction into the mathe-
matics of the Keller–Segel (KS) model for chemotaxis. In its original form this model
consists of four coupled reaction-advection-diffusion equations. These can be reduced
under quasi-steady-state assumptions to a model for t wo unknown functions u and v
which will form the basis for our study in this article. The general form of the model
is
u
t
=∇(k
1
(u,v)∇u − k
2
(u,v)u∇v) + k
3
(u,v),
v
t
= D
v
v + k
4
(u,v)− k
5
(u,v)v,
(1)
where u denotes the cell (or organism) density on a given domain  ⊂ R
n
and v
describes the concentration of the chemical signal. The cell dynamics derive from
123

A user’s guide to PDE models for chemotaxis 185
population kinetics and movement, the latter comprising a diffusive flux modelling
undirected (random) cell migration and an advective flux with velocity dependent on
the gradient of the signal, modelling the contribution of chemotaxis. k
1
(u,v)describes
the diffusivity of the cells (sometimes also called motility) while k
2
(u,v)is the che-
motactic sensitivity; both functions may depend on the levels of u and v. k
3
describe
cell growth and death while the functions k
4
and k
5
are kinetic functions that describes
production and degradation of the chemical signal. A key property of the above equa-
tions is their ability to give rise to spatial pattern formation when the chemical signal
acts as an auto-attractant, that is, when cells both produce and migrate up gradients
of the chemical signal.
1.1 Derivation and applications of chemotaxis models
Whilst a number of further approaches have been developed (e.g. stochastic and
discrete approaches such as those in [22,42,64,66,84,99]), it is the deterministic
Keller–Segel continuum model that has become the prevailing method for representing
chemotactic behaviour in biological systems on the population level. A large amount
of effort has been expended on explaining their origin from a mechanistic/microscopic
description of motion. The review by Horstmann [40] considers five methods in detail
and we refer to this work for further details of this significant area. Briefly, these are
(i) arguments based on Fourier’s law and Fick’s law [45], (ii) biased random walk
approaches [78], (iii) interacting particle systems [97], (iv) transport equations [2]or
[35], and (v) stochastic processes [86]. A more recent derivation from multi-phase
flow modelling has been proposed by Byrne and Owen [15].
As mentioned above, Keller–Segel type equations have become widely utilised in
models for chemotaxis, a r esult of their ability to capture key phenomena, intuitive
nature and relative tractability (analytically and numerically) as compared to dis-
crete/individual based approaches. To illustrate the breadth of this field, we describe
some of those areas that have benefited from the use of KS equations, apologising to
those whose works have been omitted for succinctness.
In response to starvation, the slime mold Dictyostelium discoideum initiates an
aggregation process conducted by relay of and migration to the chemical cAMP. A
number of models have been developed based on systems of equations similar to
(1) that successfully capture many key features of the lifecycle [38,44]. Understand-
ing bacterial pattern formation has also benefited from modelling: certain bacteria,
including E. coli and S. typhimurium, can be induced to form a variety of spatial
patterns when provided a suitable environment [11,12,108]. Mathematical models
indicate a chemotactic process in which cells produce an auto-attractant may underlie
this patterning ([101,108], see also [68]). Models based on the Keller–Segel equations
have also been developed to understand whether chemotaxis may underpin embryonic
pattern forming processes, such as the formation and dynamics of the primitive streak
(an early embryonic s tructure that co-ordinates tissue movements) [81], pigmentation
patterning in snakes [69] and fish [83] and cell colonisation and neural crest migration
[54].
123

186 T. Hillen, K. J. Painter
Keller–Segel type models have been developed by Lauffenburger and others [3,57]
to describe the inflammatory response of leukocytes to bacterial infection and by var-
ious authors [14,96] to model their migration in a Boyden-chamber. Modelling the
role of chemotaxis in pathological processes is a large field: Luca et al. [61] consid-
ered whether the chemotactic aggregation of microglia may provide a mechanistic
basis for senile plaques during progression of Alzheimer’s disease, while chemotaxis
has been incorporated into the modelling of a number of distinct stages of tumour
growth, including the migration of invasive cancer cells [88], tumour-induced angio-
genesis (see the reviews [16,63]) and macrophage invasion into tumours [79]. Finally,
we should not neglect the modelling of taxis in the context of spatial ecological pro-
cesses, including “prey-taxis” [43,58], herd grazing [31] and the spatial dynamics of
mountain pine beetle attacks [60].
1.2 The model variations
In addition to their utilisation within models for biological systems, a large body of
work has emerged on the mathematical properties of the Keller–Segel equations (1)
and, in particular, on t he conditions under which specialisation or variations of (1)
either form finite-time blow-up or have globally existing solutions. The majority of
this work has been devoted to a special case of (1), in which the functions k
j
are
assumed to have linear form [see model (M1) below], a model we shall refer to as the
minimal model following the nomenclature of Childress and Percus [18].
The minimal model has rich and interesting properties including globally exist-
ing solutions, finite time blow-up and spatial pattern formation. Detailed reviews can
be found in the survey of Horstmann [40], and in the textbooks of Suzuki [98] and
Perthame [87]. We shall discuss further details of these aspects in Sect. 2.1.
The minimal model is derived according to a limited set of conjectures and a num-
ber of variations have been described based on additional biological realism. In this
paper, we systematically consider some of these variations. For obvious reasons, it is
impossible to cover all variations and, as an example, the limiting case of zero diffu-
sion in the chemical signal studied by Levine and Sleeman [59] and coworkers will not
be considered here (although we will discuss this case further within the concluding
discussion). A non-diffusing signal represents an immovable entity and is perhaps
more appropriately an example of haptotaxis.
The variations are each introduced in a form that includes a single additional param-
eter that, under an appropriate limit, reduces the system to the minimal f orm. In many
cases this modification regularises the problem such that solutions exist globally in
time. Hence we call the corresponding parameter for each of the extended models the
regularisation parameter. The regularisation parameter allows us to study in detail
bifurcation conditions, pattern formation and properties of the nonuniform solutions.
We will not discuss certain questions such as the convergence of solutions of the vari-
ations to the minimal model in the corresponding limit case and leave this for future
studies.
Below we list the ten models studied in this paper. We give explanations, motiva-
tions and literature references in Sect. 2.
123

A user’s guide to PDE models for chemotaxis 187
The minimal model
u
t
=∇
(
D∇u − χu∇v
)
,
v
t
=∇
2
v + u − v.
(M1)
Signal-dependent sensitivity models
We study two versions of signal-dependent sensitivity, the “receptor” model,
u
t
=∇

D∇u −
χu
(1 + αv)
2
∇v

,
v
t
=∇
2
v + u − v,
(M2a)
where for α → 0 the minimal model is obtained, and the “logistic” model
u
t
=∇

D∇u − χu
1 + β
v + β
∇v

,
v
t
=∇
2
v + u − v,
(M2b)
where for β →∞the minimal model follows and for β → 0 we obtain the classical
form of χ(v) = 1/v.
Density-dependent sensitivity models
We study two models with density-dependent sensitivity, the “volume-filling” model,
u
t
=∇

D∇u − χu

1 −
u
γ

∇v

v
t
=∇
2
v + u − v,
(M3a)
where the limit of γ →∞leads to the minimal model, and
u
t
=∇

D∇u − χ
u
1 + u
∇v

,
v
t
=∇
2
v + u − v,
(M3b)
where ε → 0 leads to the minimal model.
The non-local model
u
t
=∇

D∇u − χu
◦
∇
ρ
v

,
v
t
=∇
2
v + u − v,
(M4)
The non-local gradient
◦
∇
ρ
v is defined in Section (2.4) and chosen such that the
minimal model follows for ρ → 0.
123

Frequently asked questions

Q1. What are the contributions mentioned in the paper "A user’s guide to pde models for chemotaxis" ?

The Keller-Segel model of chemotaxis ( Keller and Segel in J Theor Biol 26:399–415, 1970 ; 30:225– 234, 1971 ) has provided a cornerstone for much of this work, its success being a consequence of its intuitive simplicity, analytical tractability and capacity to replicate key behaviour of chemotactic populations. This phenomenon has been shown to lead to finite-time blow-up under certain formulations of the model, and a large body of work has been devoted to determining when blow-up occurs or whether globally existing solutions exist. In this paper, the authors explore in detail a number of variations of the original Keller–Segel model. The authors review their formulation from a biological perspective, contrast their patterning properties, summarise key results on their analytical properties and classify their solution form. The authors conclude with a brief discussion and expand on some of the outstanding issues revealed as a result of this work. 

Q2. Why do the authors restrict their numerical exploration to just 2 cases?

Due to the computationally exhaustive nature of 3D simulations the authors restrict their numerical exploration to just 2 cases: the minimal model (M1), for which finite time blow-up occurs in 3D, and the volume-filling model (M3a), which is known to have globally existing solutions [36,110]. 

Q3. What is the key to the success of the Keller–Segel model?

Fundamental to the success of the Keller–Segel model is its ability to demonstrate aggregation in certain parameter regions, a phenomenon that can be intuitively understood through the process of chemotactic migration up gradients of a self-produced chemical. 

Q4. What is the role of chemotaxis in multicellular organisms?

In multicellular organisms, chemotaxis of cell populations plays a crucial role throughout the life cycle: sperm cells are attracted to chemical substances released from the outer coating of the egg [33]; during embryonic development it plays a role in organising cell positioning, for example during gastrulation (see [26]) and patterning of the nervous system [85]; in the adult, it directs immune cell migration to sites of inflammation [111] and fibroblasts into wounded regions to initiate healing. 

Q5. What is the main method for demonstrating global existence of solutions?

Exploration of the literature reveals two principal methods for demonstrating the global existence of solutions; (i) finding an L∞ a-priori estimate for the chemotaxis term in the population flux, i.e. the term −k2(u, v)u∇v in (1), and (ii) to find a Lyapunov function. 

Q6. What is the case for a kinetic term describing an Allee effect?

The case for a cell kinetics term describing an Allee effect ( f (u) = u(1−u)(u−α)) has been studied in a series of papers [23,32,48,65], where it is assumed that the kinetics term acts on a faster time scale compared to that of cell movement. 

Q7. What is the qualitative behaviour of solutions on bounded domains?

For (M1) on bounded domains it has been shown that the qualitative behaviour of solutions strongly depends on the space dimension. 

Q8. What is the simplest way to model the volume filling model?

the volume-filling model (M3a) reduces to the minimal model as γ → ∞, corresponding to allowing an unlimited number of cells to accumulate at each location. 

Q9. What is the simplest way to model volume filling?

Following the derivation of the continuum limit and nondimensionalising as before, the authors obtain the general volume filling modelut = ∇ (D (q − uqu)∇u − χuq(u)∇v) , vt = ∇2v + u − v.(5)Their prototype volume filling model (M3a) appears for the simple and plausible choice of q(u) = 1 − u/γ , for 0 ≤ u ≤ γ , where γ ≥ 1 denotes the maximum cell density. 

Q10. What is the minimum model of chemotactic sensitivity?

With this version, the authors obtain the minimal model as β → ∞, while for β = 0 the authors obtain the “logistic” chemotactic sensitivity mentioned briefly above.