Hybrid modelling of individual movement and

collective behaviour

Benjamin Franz and Radek Erban

Abstract Mathematical models of dispersal in biological systems are often written

in terms of partial differential equations (PDEs) which describe the time evolution
of population-level variables (concentrations, densities). A more detailed modelling
approach is given by individual-based (agent-based) models which describe the be-
haviour of each organism. In recent years, an intermediate modelling methodology
– hybrid modelling – has been applied to a number of biological systems. These
hybrid models couple an individual-based description of cells/animals with a PDE-
model of their environment. In this chapter, we overview hybrid models in the litera-
ture with the focus on the mathematical challenges of this modelling approach. The
detailed analysis is presented using the example of chemotaxis, where cells move
according to extracellular chemicals that can be altered by the cells themselves. In
this case, individual-based models of cells are coupled with PDEs for extracellu-
lar chemical signals. Travelling waves in these hybrid models are investigated. In
particular, we show that in contrary to the PDEs, hybrid chemotaxis models only
develop a transient travelling wave.

1 Introduction

There are two fundamentally different approaches to the mathematical modelling

of systems of interacting individuals (cells, animals) in biology. If the number of
individuals is large, one often uses a continuum population-level approach, which
yields partial differential equations (PDEs) for the spatially-distributed densities of
individuals [35]. The advantage of PDE-based modelling is a well-developed mathe-
matical theory and a number of existing numerical solvers which can be used to effi-
ciently simulate the system behaviour. However, continuum approximation becomes
inaccurate if smaller groups of individuals are studied, and agent-based (individual-

Mathematical Institute, University of Oxford, 24-29 St Giles’, Oxford OX1 3LB, United Kingdom;
e-mail: [email protected], [email protected]

1
2 Benjamin Franz and Radek Erban

based) models become the method of choice [10, 44]. Examples can be found in
zoological applications, like behaviour of fish schools, bird flocks and locust groups
[9, 45]. The individual behaviour of the agents is modelled as well as the interac-
tion (e.g. attraction or repulsion) between them [11]. A number of these agents are
then simulated on the computer and their collective behaviour is analysed. This ap-
proach allows for a more detailed description of the individual behaviour and does
not discount various stochastic effects caused by a finite number of individuals. On
the other hand mathematical analysis is often hard to achieve and simulations can
be computationally intensive.
Another problem with purely agent-based models is that it is challenging to in-
corporate influences the agents might have on their environment. This is important
whenever agents interact indirectly by modifying their (evolving) environment. A
classical example is modelling chemotaxis where individual cells modify (secrete,
consume) extracellular chemical signals which diffuse in the extracellular space
[12, 18]. In this case, a hybrid modelling framework that seeks to combine the ad-
vantages of continuum and agent-based models is often used. The main idea of this
modelling approach is to describe some species as a continuum and some species as
a set of agents. For example, Dallon and Othmer [12] developed a hybrid model for
chemotaxis of slime mold Dictyostelium discoideum in which the cells are treated
as individuals in a continuum field of the chemoattractant which evolves according
to a reaction-diffusion PDE. A similar hybrid modelling framework has also been
applied to chemotaxis of bacteria [13, 51] and leukocytes [22]. The use of the hy-
brid approach allows for faster simulations than the purely agent-based model which
would treat extracellular chemicals as another set of agents. Extracellular signalling
molecules are much smaller and more abundant than cells. This property is often
used to justify that extracellular chemicals can be described as a continuum [12].
The use of hybrid models is becoming more widespread especially with the grow-
ing computational power that allows to consider more complex systems in this man-
ner, including modelling tumour growth [39] and forest dynamics [33]. In cancer bi-
ology, several hybrid cellular automaton models have been proposed in the literature
[40, 41]. For example, Smallbone et al. [41] coupled a two-dimensional cellular au-
tomaton model (describing cells) with continuum (PDE-based models) of glucose,
H+ and oxygen concentrations, building on the previous work of Patel et al. [39]
and Alarcón et al. [3]. A similar hybrid approach has been used in a number of other
studies in cancer biology [4, 21, 36]. A hybrid forest model with trees modelled as
agents and a continuum approach used for oxygen and other atmospheric gases is
presented in [33]. In economical research hybrid models are used to estimate prices
in the petrol market [25] and in general markets with a non uniform spatial demand
of products [26, 27]. In these models the demand is described as a continuous func-
tion of space whereas the retailers are considered as agents.
The term hybrid modelling is sometimes applied for models which use both
individual-based and continuum description for the same physical quantity. For ex-
ample, a “hybrid” model for the spread of an epidemic disease is presented in [8]. It
initially considers infected individuals as agents, but switches to a continuum model
when the number of infected people in an area rises above a threshold. Coupling
Hybrid modelling of individual movement and collective behaviour 3

reaction-diffusion models with a different level of detail in different parts of the

computational domain is presented in [15, 19]. “Hybrid” models of this type are
useful because they can lead to computational savings. However, in this chapter, we
will focus on hybrid models which describe some system components (e.g. cells
or animals) as individual agents and some components (e.g. external chemicals) as
continuum fields. The choice which description is used for each species is made at
the beginning and will not change during the course of the simulation. We will sum-
marise the progress in hybrid models which satisfy this definition, and clarify some
of the problems and difficulties that arise from their use.
The outline of this chapter is as follows. Section 2 will give a short overview
of the PDE-based and agent-based modelling approaches before the general math-
ematical framework for hybrid models is introduced in Section 3. Hybrid models
can be considered as extensions of (purely) agent-based models. Therefore, their
computer implementation often forms an integral part of the model. We will dis-
cuss it in detail in Section 4 where we describe the numerical simulation of hybrid
models drawing special attention to the different treatment of the continuum and the
agent-based subsystems as well as the problem of matching the two parts. In order
to give a more practical insight into the topic we will perform a case study of a hy-
brid chemotactic model in Section 5. This case study will also be used to show some
qualitative and quantitative differences that can occur when using a hybrid model
instead of the corresponding continuum model.

2 Continuum vs. agent-based models

Hybrid modelling is an intermediate approach between continuum (PDE-based)

models and agent-based models of systems of interacting individuals. In this section
we briefly review these common modelling approaches in mathematical terms. We
will make use of our notation later in Section 3 when hybrid models are considered.
Continuum (mean-field) models give rules for the evolution of the spatially de-
pendent concentration vector c ≡ c(x,t) where x ∈ Ω ⊂ Rm , m = 1, 2 or 3, and t is
the simulation time. The components of the vector c can be densities of individu-
als (cells, animals) and concentrations of extracellular signals. As the concentration
vector c can change both with position x and time t, a general continuum model
takes the form
∂c
= L (c, x,t) x∈Ω, (1)
∂t
where L is an operator on c, which in most practical cases will be a differential
or integral operator. To uniquely describe the time evolution of (1), one also has to
specify suitable initial and boundary conditions.
Example 1. Keller-Segel model
Continuum modelling is used in many areas of mathematical biology [35]. In
chemotaxis modelling (which will be the subject of Section 5), a classical exam-
ple of (1) is the Keller-Segel model of chemotaxis [31]. Here, Ω ⊂ R and the vector
4 Benjamin Franz and Radek Erban

c has two components, i.e. c = [c1 , c2 ] = [n, S] where n ≡ n(x,t) is the density of cells
and S ≡ S(x,t) is the concentration of the chemoattractant. The evolution equation
(1) is a coupled system of two PDEs for n and S:

∂n ∂ 2n ∂ ∂S
 
= Dn 2 − nχ (S) , (2)
∂t ∂x ∂x ∂x
∂S ∂ 2S
= DS 2 − k(S)n , (3)
∂t ∂x
where Dn and DS are diffusion constants of cells and chemoattractant, respectively.
The strength of chemotaxis is controlled by chemotactic sensitivity χ (S) and there-
fore by the concentration of substrate S which is consumed by cells with the rate
k(S).

The applicability of continuum modelling depends on the number of particles in

the studied system. In Example 1, the interacting “particles” are unicellular micro-
scopic organisms (n) and molecules of chemical signal (S). As there are often more
signalling molecules than cells, the validity of mean-field assumptions is dictated
by the number of cells in the system. If the system only consists of a few cells, it is
more accurate to use an individual-based approach which is introduced in the next
section.

2.1 Agent-based modelling

In contrary to the continuum models the so-called agent-based models treat every
particle as an individual that follows an inherent set of rules. This means in par-
ticular that individual behaviour and interactions between different agents account
for the possibly complex behaviour of the system. Agent-based models are com-
monly used for systems with a small number of individuals that follow non-trivial
behavioural rules, for example in modelling of collective animal behaviour [11]
or human crowds in panic situations [23]. While continuum models have a well-
developed mathematical theory, agent-based models are sometimes written as com-
puter routines which are difficult to theoretically analyse. The literature also fails to
agree on a general definition of an agent. In this chapter, we use a definition which
is slightly adopted from [48] and used in [20].

Definition 1. An agent is a system that uses a fixed set of rules based on communi-
cation with other agents and information about the environment in order to change
its internal state and fulfil its design objective.

This definition, however, is only a formal description, which now has to be put into
a more rigorous context. Following from Definition 1, the mathematical description
of an agent has to incorporate the behavioural rules of an agent as well as the pos-
sibility of communication between them. Therefore, we assume a finite number N
Hybrid modelling of individual movement and collective behaviour 5

of agents numbered from 1 to N. In general N can depend on time, taking into ac-
count birth or death of agents. We define the current state of an agent by its internal
state variable yi (t), i = 1, . . . , N, which can describe its position, velocity and inter-
nal memory. It is this internal state and its time evolution that describes the rules
of an agent. Since these agents represent different individuals, we assume that other
agents generally have no means to access all internal state variables. In order to al-
low for communication between the agents, we define a set of external states wi (t),
which are observable by other agents. The observable states wi (t) of every agent are
in principle available to every other agent, which is ensured by creating the set of
external states X . The general agent-based model following these definitions then
takes the form

yi (t + ∆ t) = fi (yi (t),t, ∆ t, X ) , i = 1, . . . , N , (4)

wi (t) = gi (yi (t)) , i = 1, . . . , N , (5)
X = {w1 , . . . , wN } . (6)

We can see that the evolution of yi is given by the function fi , which notably de-
pends on the time step ∆ t. This general description can entail discretised versions
of ordinary differential equations (ODEs) as well as stochastic differential equa-
tions (SDEs). Additionally, agent-based systems that only change discretely can be
written in the form (4)–(6).
We understand the external states of an agent merely as an observable represen-
tation of the internal states, which is why wi (t) directly depends on yi (t) through the
function gi . The distinction between observable and non-observable states is often
used to represent internal memories that cannot be perceived by other agents.

Example 2. Animal behaviour

Agent-based models have been successfully used for the modelling of collective
animal behaviour [45]. Couzin et al. [11] showed that a relatively simplistic model
can yield complex collective behaviour and can be used to model fish schools and
bird flocks.
In this model, the internal states of an agent yi are defined to be its position
xi ∈ Rm (m = 2, 3) and its velocity vi ∈ Rm . Since both the position and velocity
of an agent potentially influence the motion of other agents, both are observable
and hence wi = yi = [xi , vi ] ∈ R2m , which means that gi = Id. The update rules
fi , i = 1, . . . , N, in this example are equivalent for each agent and incorporate the
different rules for the different zones in the model (zone of attraction, orientation
and repulsion).

Example 3. Chemotactic movement under a stationary signal

A simple agent-based model for chemotaxis in one dimension can be written as
follows [24]: the internal state yi (t) of an agent is defined as its current position in
R. Additionally, we assume that the signal S(x) is fixed and there is no interaction
between agents, hence no observable states are required. All agents start at some
initial position y0,i ∈ R and move according to the stochastic differential equation
6 Benjamin Franz and Radek Erban

∂S
dyi (t) = χ (S)
p
dt + 2Dn dW , i = 1, . . . , N , (7)
∂x
where χ (S) is the chemotactic sensitivity function introduced in Example 1, Dn is
the diffusion constant of the bacteria and dW is the Wiener-process, also known as
Brownian motion [29]. We can discretise (7) to obtain an update rule equivalent to
(4) as follows
∂S
yi (t + ∆ t) = yi (t) + χ (S) ∆ t + 2Dn ∆ t ξ ,
p
∂x
where ξ is a normally distributed random variable with zero mean and unit variance.
In the limit of infinitely many particles, this agent-based description is equivalent to
the PDE (2), which is written for the density of cells. However, if we considered a
time-evolving signal which is consumed by cells as in Example 1, a purely agent-
based model would have to simulate the trajectories of all signal molecules. This
would be computational intensive and a hybrid model which combines agent-based
simulations with PDEs can then be used to optimize computational efficiency and
accuracy.

3 Hybrid modelling – theoretical framework

Because of their hybrid nature the general framework for these models necessarily
combines the two frameworks presented in Section 2. We define a vector of contin-
uous variables c(x,t) on a domain Ω ⊂ Rm , m = 1, 2 or 3. The update rule for c is
again governed by an operator L , which now also depends on the current states of
the agents. The N agents are represented by their internal state variables yi (t) and
their set of observable states wi (t) defined in (5). To allow interactions between the
agents and the continuous variables c, the set of observable states X as defined in
(6) is used. The update rules for the system are

∂c
= L (c, x,t, X ) , x∈Ω, (8)
∂t
yi (t + ∆ t) = fi (yi (t),t, ∆ t, X , c) , i = 1, . . . , N , (9)

where X is given in (6). In (8) we see that the agents can influence the continuous
variables c through the set of observable states X . Similarly, the behaviour of the
agents can be altered by the continuous variables, as the operator fi now also depends
on c. Figure 1 shows a graphical representation of the hybrid model. It contains the
N agents represented by the internal states yi on the left. Through the function gi the
observable states wi are generated which then influence the update of the continuous
variables c as well as the agents’ behaviour themselves. We, however, encounter a
problem in this definition, as the continuous variables are defined for every time t,
while the internal agent states are only defined for discrete times. To overcome this
problem we can consider (9) in the limit ∆ t → 0, where it takes the general form of
an SDE
Hybrid modelling of individual movement and collective behaviour 7

(1) (2)
dyi = fi (yi (t),t, X , c) dt + fi (yi (t),t, X , c) dW ,
(1) (2)
where fi and fi respectively represent the stochastic and the deterministic part of
the SDE.

f1 y1 g1
w1
L c
X
wN
fN yN gN

Fig. 1 Concept of a hybrid model. Arrows symbolise direction of influence.

Example 4. Hybrid cellular automaton model for carcinogenesis

In [41] Smallbone et al. present a hybrid cellular automaton model for the formation
of cancer. This model uses reaction-diffusion equations to calculate the concentra-
tion of oxygen, glucose and hydrogen ions in the environment of the cells. The
concentrations of these chemicals therefore constitute the continuous variables c.
Each cell of the cellular automaton is represented by an agent with the internal state
yi ∈ N defining which of the finite number of possible phenotypes the cell at this
position has (including the “phenotype” empty). As these phenotypes are observ-
able by neighbouring cells, we have wi = yi . This cellular automaton model has a
generation-based update rule, which means that the states yi are only updated once
every time step. The rules of the model then represent the probabilistic functions fi
in equation (4), where the change depends on the current phenotype, the neighbour-
ing cells and the concentrations of the considered chemicals at the cell position.

Example 5. Hybrid model for chemotaxis of Dictyostelium discoideum

Dallon and Othmer developed a hybrid model for the chemotaxis of Dictyostelium
discoideum [12] that combines individual cell movement with a continuous extra-
cellular concentration of cAMP modelled by a PDE. The internal states of the agents
are the position of the individual xi , as well as the variables representing the intra-
cellular processes. Only the position and one of the intracellular variables influence
the external field and therefore form the observable states wi . The update rules fi are
given through ODEs for the internal dynamics and rules of motion for the position.
8 Benjamin Franz and Radek Erban

3.1 A position-based hybrid model

So far, we have defined a general framework for a hybrid model that allows for
a great freedom in the choice of internal and external states of the agents. In the
next step we want to refine this framework for the more specific models used in
chemotaxis modelling [12, 17, 43]. In order to be able to interpret the agents as part
of a species situated inside the domain Ω , we need to introduce the notion of an
agent’s position in Ω . Moreover, we assume that all agents are equal for an external
observer except for their position, or in other words the set of observable states of
the agents wi (t) is the position xi (t) of the agents inside Ω , i.e.

wi (t) ≡ xi (t) .

This definition excludes Couzin et al. models for animal behaviour [11] as well
as cellular automaton models [41], but it is sufficient for the chemotaxis example
studied in Section 5.
Because of the agents’ similarity, we no longer need to define an abstract set X ,
but can instead define a density function ρδ on Ω through
N
ρδ (x,t) = ∑ δ (x − xi (t)) , x∈Ω. (10)
i=1

When discussing numerical simulations of hybrid models, we will see that this def-
inition of ρδ is already a first step towards obtaining a continuous density function
for the agents. With this definition we can redefine the operator L , which governs
the behaviour of the continuous variables c and (1) reads as follows

∂c
= L (c, x,t, ρδ ) .
∂t
For the evolution of the internal agent states yi we assume now that every agent can
only perceive information about the continuous variables c at its current position.
Hence, the operator fi no longer depends on c on the whole domain, but only on
c(xi ) and the first spatial derivative in this point, i.e. fi , i = 1, . . . , N, are functions
for all further considerations. Equation (9) therefore becomes

yi (t + ∆ t) = fi (yi ,t, ∆ t, ρδ , c(xi ,t), ∇c(xi ,t)) . (11)

This special type of hybrid systems still allows for a wide range of flexibility and
can therefore be used to model a variety of different processes. In Section 5 we study
position-based models for chemotaxis in more depth.
Hybrid modelling of individual movement and collective behaviour 9

3.2 Initial and boundary conditions

An important aspect of modelling is the incorporation of initial and boundary condi-

tions. Hybrid models necessarily combine the conditions from the two different ap-
proaches. For the continuous variables one usually has an initial value c0 (x), while
for the agents an initial distribution of their position and internal states is given,
which is then used to generate each agents’ position at the beginning of the simula-
tion. In some applications the agents can be born during the course of the simulation.
In this case, we have to ensure the appropriate initialisation of its internal variables.
A similar idea of independent conditions for the continuum and the agent-based
parts of the hybrid model is used for the boundary conditions. The values of the
continuous variables on the boundary usually have to satisfy an equation of the type

G (c, x,t) = 0 , x ∈ ∂Ω , (12)

where G is a general operator. In the most commonly used cases (12) enforces
certain values on c or its gradient on the boundary. For the agents the boundary
conditions are often given in a more descriptive manner. For example, agents can
leave the domain through one end and automatically reappear on the other end. This
periodic boundary condition implies that the number of agents in the system is con-
served. Periodic boundaries are widely used because of their simplicity and because
they effectively shape an infinite domain. Reactive boundaries absorb agents with a
probability p, while reflecting them with probability 1 − p [14]. If p = 0, one often
speaks of a reflecting boundary, while for p = 1 the condition is called an absorbing
boundary.

4 Hybrid modelling – numerical implementation

For similar reasons as in purely agent-based models it is often very hard to obtain
analytic results for hybrid models. This increases the importance of numerical simu-
lations for gaining insight into the behaviour of the system. The mixture of different
modelling frameworks, however, renders the process of setting up a numerical simu-
lation non-trivial. Each part of the model has to be considered differently and a way
of matching the two parts has to be developed. In this section we discuss a numeri-
cal framework and evaluate difficulties one has to overcome when implementing a
hybrid model.
The general task for the numerical simulation of a hybrid model is to calculate
approximations for both c and yi at times t j = j∆ t, j = 1, 2, . . . given initial data for
each of these variables according to Section 3.2. We additionally assume that the do-
main Ω can (for the continuous part of the hybrid model) be adequately represented
by the points r1 , . . . , rL ∈ Ω , which means that we seek to compute approximate val-
ues for c(t j , rl ), j = 1, 2, . . ., l = 1, . . . , L and yi (t j ), i = 1, . . . , N. In order to simplify
the notation, we introduce
10 Benjamin Franz and Radek Erban

C j = [c(t j , r1 ), . . . , c(t j , rL )] j = 0, 1, . . . .

Due to the different characters of the continuous and the agent-based subsystems,
different approaches have to be used for their numerical solutions. For each of the
subsystems one tries to answer the question of how to get from t j to t j+1 still guar-
anteeing an accurate approximation of the system. For the continuous variables this
means, we seek a solver that generates the values of C j+1 using the values C0 , . . . ,C j
and the current distribution of the agents ρδ (·,t j ) given by (10), which can be sym-
bolised as
Ld
C0 , . . . , C j , ρδ (·,t j ) 7−→

C j+1 . (13)
In (13) we introduced the operator Ld , which is a discretised version of the contin-
uous operator L used in (8). In the most common case, where L is a differential
operator, Ld could be a finite element or finite difference approximation of L . Note
that in (13) we have made the implicit assumption that the solver used for (8) only
takes the positions of the agents at time t j into account. For the agents equation (11)
is already given in a time-discrete way and can therefore be used directly to update
the internal states.
The introduction of this general scheme raises some immediate problems, which
we will discuss in the remainder of this section. The first difficulty are the differing
spatial resolutions for the two subsystems, which we address in Section 4.1. Other
problems like time stepping, choices of solvers and the influence of stochastic ef-
fects are presented in Section 4.2.

4.1 Spatial matching in numerical simulations

A spatial matching between the continuous variables and the agents is required dur-
ing a numerical simulation of a hybrid system, because different spatial resolutions
are applied. The agents can be positioned at an arbitrary point inside the domain
Ω , while the data for c is only calculated at the points rl . This triggers a two-way
matching problem, as one has to generate estimates for the agent distribution at the
points rl as well as for the continuous variables c everywhere inside Ω .
First, let us consider estimating the agent density distribution throughout Ω and
especially at the points rl , which is necessary for the update relation (13). So the
general mapping we are trying to achieve is
N
Ψ
ρδ (x) = ∑ δ (x − xi ) 7−→ ρ (x) ∈ C0 (Ω ) .
i=1

The requirements for the estimated density function ρ (x) can alter for different ap-
plications, but here we require it to be at least a continuous function. One way to
achieve such a mapping is the so-called kernel density estimation [46]. In general
the kernel density estimation can be used to estimate the probability density function
of a random process, if one has been given a number of realisations of this process.
Hybrid modelling of individual movement and collective behaviour 11

The name stems from the use of a kernel K(x), which is typically a continuous,
symmetric and normalised function. Let us for simplicity assume a one dimensional
random process, in which case these conditions take the form
Z
K(x) ∈ C0 (R), K(−x) = K(x), K(x)dx = 1 . (14)
R

Additionally, K(x) is often required to be non-negative in order to generate a non-

negative estimate. Most commonly used kernels include a Gaussian kernel and a
piecewise linear kernel with compact support. In practice a scaled version of K is
used, which leads to the introduction of a bandwidth parameter h. We define
1 x
Kh (x) = K ,
h h
which still satisfies the conditions (14). With given positions x1 , . . . , xN , an estimate
of the probability density function is then given by
N
ρ (x) = ∑ Kh (x − xi ) = Kh (x) ∗ ρδ (x) . (15)
i=1

Figure 2 shows an example of a kernel density estimation for 100 normally dis-
tributed random variables using a Gaussian kernel with different bandwidths h. In
Figure 2(a) we can see that the choice of a very small h leads to a highly oscillating
estimate, while a very big h can lead to the estimate being too wide as shown in Fig-
ure 2(c). An optimal choice for the parameter h and the kernel itself always depends
on the nature of the problem and the number of samples N.

0.6 0.6 0.6

0.5 0.5 0.5

0.4 0.4 0.4

ρ

0.3 0.3 0.3

0.2 0.2 0.2

0.1 0.1 0.1

0 0 0
−5 0 5 −5 0 5 −5 0 5
x x x

(a) h = 0.1 (b) h = 0.25 (c) h = 0.9

Fig. 2 Kernel density estimate for N = 100 agents, which are placed according to a normal distri-
bution with different bandwidths h. Crosses along the x-axis represent the agents, the dashed line
is the underlying Gaussian probability density function and the solid line is the generated estimate
by (15).

The second spatial matching problem that occurs when simulating a combined
continuous and agent-based system is the need to estimate the values of the contin-
uous variables (and possibly their derivatives) at an arbitrary position inside Ω . The
operator we are looking for can be symbolised through
12 Benjamin Franz and Radek Erban

Θ
(r1 , c(r1 )) , . . . , (rL , c(rL )) 7−→ ĉ(x) ∈ C0 (Ω ) .

Though similar to the operator Ψ , we here have the advantage that we know the
positions of the points r1 , . . . , rL beforehand and that we know they give an adequate
representation of the domain Ω . With this additional information, one can argue that
the problem at hand represents an interpolation problem from the grid points rl onto
the whole domain Ω . This result allows for the use of approaches from the well-
studied fields of interpolation and approximation theory [47]. In some cases the
interpolation regime is already implicitly incorporated in the numerical solution of
the update equation for the continuous variables, for example if one chooses to use
a finite element approach.

Example 6. Numerical realisation of Example 5

In Example 5 we presented a hybrid model for chemotaxis of slime mold Dic-
tyostelium discoideum developed by Dallon and Othmer [12]. To generate a discre-
tised operator Ld they used the particle-in-cell method [38]. For the kernel density
estimation Ψ they use a piecewise linear kernel and for the interpolation operator Θ
a fifth order spline interpolation was employed.

4.2 Other aspects of numerical simulations

The spatial matching between the two parts is the biggest additional challenge posed
by the use of a hybrid model. Here, we discuss some other problems that occur dur-
ing this process. The first problem is the choice of a solver both for the continuous
variables and for the internal states of the agents. One can choose from a wide range
of standard approaches for both problems. The way the two parts are interwoven,
however, sets some restrictions. It is, for example, almost always impossible to use a
fully implicit solver for both parts, especially if the functions fi for the internal agent
states contain random variables. Additionally, one has to consider the accuracy of
the different solvers and should ideally try to match these to prevent unnecessary
computational effort that does not lead to more accurate results.
The discrete nature of the agent-based parts automatically introduces stochastic
effects into the system. Various examples of these effects will be discussed in Sec-
tion 5. It is important to consider these effects when choosing the time stepping and
the spatial resolution for the simulation. In particular, these choices will depend on
the number of agents in the system. It is generally possible to allow different time
steps for different parts of the system, for example the agents could be simulated
with a finer time stepping than the continuous variables or vice-versa. For each part
of the system the time steps have to be chosen in a way that ensures an accurate
solution depending on the spatial resolution and the solver that is used. In Section 5
we study one application area of hybrid systems in more detail and analyse the effect
of some of these choices on the system.
Hybrid modelling of individual movement and collective behaviour 13

5 Case study: Hybrid modelling of chemotaxis

In Section 3 we introduced a general framework for hybrid models that combine

agent-based models with mean-field equations and we now concentrate on one ap-
plication area for hybrid modelling: cell migration. In particular we focus on the
movement of cells induced by gradients in the concentration of extracellular chemi-
cals, a process that is known as chemotaxis. Chemotaxis is one of the main forms of
cell migration and is used in a variety of cells, including bacteria cells [7]. Hybrid
models of chemotaxis have been successfully used in the literature [12, 13, 22, 51].
The first notion of chemotaxis goes back to the late 19th century, when Engel-
mann and Pfeffer detected the process. In the late 1960s it was Adler [1, 2] who per-
formed experiments with the bacteria E. coli that helped understanding and quan-
tifying the process and was later used as comparison for the early mathematical
models. Adler placed a colony of E. coli at one end of a long thin pipe that was
filled with oxygen and an additional energy source. Through the process of chemo-
taxis the colony started to move with a constant speed away from the closed end
forming a narrow band of bacteria. The band was visible to the naked eye and Adler
was able to measure the speed with which it moved forward.
In the 1970s the first mathematical descriptions of chemotaxis were formulated,
with the Keller-Segel model, which we will discuss in Section 5.1, as one of the
early breakthroughs. A review of the impact this first model had on the modelling
of chemotaxis is given in [28]. Section 5.2 will introduce a hybrid version of this
model, which we will further investigate and analyse in Section 5.3.

5.1 The Keller-Segel model

As mentioned above, Keller and Segel developed the first mathematical model to de-
scribe the process of chemotaxis in 1971 [31]. The original model considers both the
bacteria and the chemotactic substrate in a continuum limit, which therefore results
in a coupled system of two PDEs. The original form of the system only considers
one spatial dimension and gives a way to compute the concentration of bacteria de-
noted by n(x,t) and the concentration of substrate S(x,t) through the PDEs (2)–(3),
introduced in Example 1. In equation (2) we can see that the behaviour of the bacte-
ria is governed by two independent effects and therefore takes the form of a general
advection-diffusion equation. The diffusion of the bacteria occurs with the diffusion
constant Dn , while the advection is governed by the chemotactic sensitivity χ (S).
The substrate, as seen in (3), diffuses with the diffusion constant DS and is consumed
by the bacteria with a consumption rate k(S) that depends on the concentration of
substrate itself.
In a follow-up to the paper [31], Keller and Segel showed that under certain
conditions the developed system of partial differential equations yields travelling
wave solutions [32]. In particular they were able to proof that travelling wave solu-
tions can only exist if χ (S) has a singularity at some critical value Scrit . For reasons
14 Benjamin Franz and Radek Erban

of simplicity they concentrated on the simplest such functions χ (S) = κS with the
critical concentration at Scrit = 0. In their analysis Keller and Segel made some ad-
ditional assumptions for the various parameter values and simplified (2)–(3) to the
nondimensionalised PDEs
∂n ∂ ∂n κ ∂S
 
=µ −n , (16)
∂t ∂x ∂x S ∂x
∂S
= −n . (17)
∂t
The nondimensionalised system is set up for x ∈ [0, 1] with an initial value of
S(x, 0) = 1 and no-flow boundary conditions. As initial distribution of the agents
we choose n(x, 0) = δ (x), which corresponds to the initial state of Adler’s experi-
ments where all bacteria were inserted at one end of the tube. We consider reflective
boundary conditions for both bacteria and extracellular signal at x = 0 and x = 1.
In order to investigate the influence of the two dimensionless parameters µ and κ
on the travelling wave, Figures 3–4 show the concentration of n and S at t = 0.5 for
various values of µ and κ . In Figure 3 we can see that the parameter µ influences
the width of the wave while leaving its general shape untouched. Increasing µ leads
to a wider wave and a decrease in the maximum of n. Accordingly, the gradient in
S is higher for the narrower bands caused by smaller values of µ . As can be seen
in Figure 4, the parameter κ influences the general shape of the wave. In the case
κ = 2 the travelling band of bacteria is symmetric, while a κ bigger than two leads
to a wave that is steeper in the front (right) and falls slowly in the back (left) of the
wave. Choosing κ smaller than two causes an opposite effect with the wave being
bent backwards.

15 1.5
µ = 1/60 µ = 1/60
µ = 1/30 µ = 1/30
µ = 1/15 µ = 1/15
10 1
n

5 0.5

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
(a) n(x, 0.5) (b) S(x, 0.5)

Fig. 3 Travelling wave solution of the Keller-Segel model (16)–(17) for different values of the
parameter µ , where κ = 2.
Hybrid modelling of individual movement and collective behaviour 15

10 1.5
κ = 4/3 κ = 4/3
κ=2 κ=2
8
κ=4 κ=4
1
6

S
4
0.5
2

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
(a) n(x, 0.5) (b) S(x, 0.5)

Fig. 4 Travelling wave solution of the Keller-Segel model for different values of the parameter κ ,
where µ = 1/30.

5.2 Hybrid models of chemotaxis

One of the assumptions made by Keller and Segel in their original model is to con-
sider the bacteria as a continuum rather than explicitly describe their individual
behaviour. For systems that do not satisfy this assumption hybrid chemotaxis mod-
els have been developed in the literature [12, 13, 22, 51]. In this section we present
three of them. The bacteria are modelled as agents with varying numbers of internal
states and their position xi ∈ Ω , as the only observable state. All three models con-
sider the substrate in a continuum limit and the PDE (17) takes the role of equation
(8) in our description of the hybrid modelling framework.

Model I

The first approach to design a hybrid version of the Keller-Segel model, is to in-
terpret the evolution equation for n as a Fokker-Planck equation for a number of
randomly moving particles similarly to the idea presented in Example 3. A chemo-
taxis model of this form was formulated by Stevens [43]. The movement of each of
the agents is described by the stochastic differential equation

κ ∂ S(xi )
dxi = µ dt + 2µ dW .
p
(18)
S(xi ) ∂ x

The parameters used in (18) correspond to the ones in the dimensionless Keller-
Segel equations (16)–(17). This particle-based description of equation (16) shows
one of the weaknesses in the original Keller-Segel model. According to (18) an agent
can theoretically jump any given distance in one time step, implying that some of
them can move with a speed that is not achievable for bacteria.
16 Benjamin Franz and Radek Erban

Model II

Driven by weaknesses of the first model, a different type of random walk, known
as velocity-jump process, seems a more realistic choice for the bacterial behaviour.
The motion of bacteria E. coli consists of two phases [7]. During a run-phase the
bacterium moves with a constant speed straight into a chosen direction. This run
lasts for a randomly distributed time before the bacterium enters the tumble-phase
in which it chooses a new direction randomly [6]. As we are considering a one-
dimensional model, there are only two possible directions of motion: to the left and
to the right. A right-moving agent continues to the right for a time that is given by an
exponentially distributed random variable before it switches its direction. In order
to incorporate the bias of bacteria towards higher concentrations of chemoattrac-
tants, Othmer et al. [37] introduced a biased velocity-jump process. In this biased
random walk the duration for the run phase depends on information gathered at the
current position of the individual. In particular, the model in [37] allows the agents
to directly measure the gradient of the substrate concentration at their current po-
sition. The run-phase then tends to be longer, if the concentration increases in the
current direction of motion, while for a decreasing signal, the turning probability is
increased.
The turning frequency λ is therefore adjusted according to the current movement
direction, the value and the gradient of S. To represent the direction of motion, the
velocity vi (t) = ±s is introduced, where s denotes the constant speed. In terms of
the hybrid modelling framework introduced in Section 3, the internal variable is
yi = [xi , vi ]. The agent-based description of the bacteria can be written in the form

xi (t + ∆ t) = xi (t) + vi (t)dt ,
− vi (t) with probability λ ± ∆ t

vi (t + ∆ t) = ,
vi (t) otherwise

where
κs ∂ S
λ ± = λ0 ∓ .
2S ∂ x
In a continuum limit this velocity-jump process is equivalent to the hyperbolic
chemotaxis equation [16]:

1 ∂ 2n ∂ n s2 ∂ ∂ n κ ∂S
 
+ = −n , (19)
2λ0 ∂ t ∂t 2λ0 ∂ x ∂ x S ∂x

where n is the concentration of bacteria. This shows that changing the type of
random-walk used for the agents can influence the corresponding continuum equa-
tion. Nevertheless (19) can be used to adjust the parameters of the agent-based
model to match the parameters of the Keller-Segel model, as the large time be-
haviour of (19) is given by the classical chemotaxis equation (16), where we have
µ = s2 /(2λ0 ) [30]. Lui et al. [34] showed that coupling the hyperbolic chemotaxis
equation (19) with (3) for the substrate also yields travelling wave solutions similar
Hybrid modelling of individual movement and collective behaviour 17

to the original Keller-Segel system. An investigation of this case for a more general
dependence of the turning frequency is given in [50].

Model III

More accurate descriptions of the individual behaviour of bacteria incorporate the

sensing and processing of extracellular signals [5, 42]. Hybrid models with descrip-
tions of these intracellular processes have been used by Dallon and Othmer [12] as
well as Xue et al. [49]. Erban and Othmer [17, 18] used an agent with a toy version of
the internal dynamics that includes two main features of the sensing process: a fast
excitation and a slower adaptation. We will use a simple model with one additional
internal variable zi that acts as a memory and allows the agent to identify increasing
or decreasing signal concentrations [17]. The model is based on a velocity-jump
process with a turning frequency λ , which depends on zi . This internal variable is
chosen to follow the value of a sensing function g(S) with the adaptation time ta .
Thus, the model can be written in the hybrid form presented in Section 3, using
yi = [xi , vi , zi ] as follows:

xi (t + ∆ t) = xi (t) + vi (t)dt ,
− vi (t) with probability λ ∆ t ,

vi (t + ∆ t) =
vi (t) otherwise ,
g(S(xi (t))) − zi (t)
zi (t + ∆ t) = zi (t) + ∆t ,
ta
where
λ = λ0 + zi − S(xi ) .
In the limit ∆ t → 0 and N → ∞ this process can be described by the chemotaxis
equation
∂n s2 ∂ ∂ n dg ∂ S
 
2ta
= − , (20)
∂t 2λ0 ∂ x ∂ x 1 + 2λ0ta dS ∂ x
provided that t is large (t ≫ 1/λ0 ) and the gradient of S is shallow [17]. Choices
for the parameters of this model can be made by matching (20) with the classi-
cal chemotaxis equation (16), which especially indicates that g is given through
dg/dS ∼ χ (S).
In Figure 5(a) a simulation of the hybrid model of type III is shown. Simula-
tions of the other two models were also performed, with results almost identical to
the one seen in Figure 5(a). We simulate N = 104 agents with the dimensionless
model parameters 1/ta = λ0 = 1.5 × 10−3 , s = 10−2 , g(S) = 4.5 × 10−3 log(S) and
∆ t = 10−4 . These parameters were chosen in such a way that they match the global
parameters µ = 1/30 and κ = 2 used for the classical Keller-Segel model. On a first
impression, it looks as though the resulting agent distribution at t = 0.5 matches the
predicted concentration of the Keller-Segel system well except for some stochastic
effects. In Figure 5(b), however, we show the agent distribution in the region behind
18 Benjamin Franz and Radek Erban

the travelling band. Further analysis of this region showed that here the extracel-
lular signal is completely exploited. Some agents are left in this zone and undergo
an unbiased random walk without a chemotactic signal to guide them. This means
that these agents do not necessarily manage to catch up with the travelling wave
again but instead stay in the exploited region. In the remainder of this section, we
study this effect, which we refer to as dropout in more detail. We will show that it
significantly influences the system dynamics for large times.

10 wave direction 0.1

8
0.08

6 0.06
n

n
4 0.04

2 0.02

0 0
0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2
x x
(a) (b)

Fig. 5 Numerical Simulation of the hybrid Keller-Segel model with internal dynamics (Model III).
Parameters are N = 104 , 1/ta = λ0 = 1.5 × 10−3 , s = 10−2 , g(S) = 4.5 × 10−3 log(S), ∆ t = 10−4 .
(a) Distribution of agents at time 0.5 (solid line) and the results given by the Keller-Segel model
(16)–(17) (dashed line, which is almost indistinguishable from the solid line).
(b) Histogram of agent positions in subinterval [0, 0.2].

5.3 Analysis of the dropout

In Figure 5(b) we saw that the hybrid model, in contrast to the original Keller-Segel
model, creates a region behind the wave where the substrate is completely exploited.
The main assumptions for a mean-field approach are violated in this region, namely
the number of bacteria and the concentration of extracellular material are very small,
which renders a continuum approach here not applicable. Stochastic effects due to
the small number of bacteria then lead to the complete exploitation of S, which
causes the dropout of some of the agents. These agents can no longer sense any
gradient in extracellular substrate and are therefore moving completely randomly,
which makes it very unlikely for them to become part of the travelling band again.
Due to the constant loss of agents, the velocity and the height of the wave will de-
crease as the wave moves along. Note that a complete exploitation in these models
is only possible under the assumption that S does not diffuse, which was made by
Keller and Segel and is incorporated in the PDE (17). The dropout effect is interest-
ing for us, because it shows a qualitative difference between the hybrid model and
the original Keller-Segel model, as the hybrid model only yields transient travelling
Hybrid modelling of individual movement and collective behaviour 19

wave solutions. In this section we create measures for this dropout in order to get
an estimate of the number of lost agents from the simulations. We will then move
on to analyse the effect of some system parameters on the dropout. Finally, some
theoretical results about the loss of agents are presented and compared to numerical
results.

Dropout measures

In order to be able to quantify the dropout of agents from the travelling wave, we
need to investigate certain conditions that render an agent as dropped out. A condi-
tion of this form allows us to define an index set Γ (t) that contains the agents who
are currently part of the wave.
However, before defining and comparing different conditions for the dropout,
we investigate some global statistical values of the agent set. The first measure to
indicate the fact that agents have dropped out is the position of the centre of the wave
c(t). From [32] we know that the theoretical wave speed of the nondimensionalised
Keller-Segel system is 1 and therefore the predicted position of the centre of the
wave is cm f (t) = t. In comparison to that the actual position of the wave can be
measured from the agents’ positions via

1 N
c1 (t) = ∑ xi (t) .
N i=1
(21)

The problem with this option is that it includes dropped out agents for the calculation
of the wave centre, which can bias the calculation. To overcome this problem, a
second option for finding the centre of the wave is given through
1
c2 (t) = ∑ xi (t) ,
|Γ | i∈
(22)
Γ

which implies that the found centre position depends on the choice for the index set
Γ . For short times c1 (t) and c2 (t) give similar results, but will differ for large times.
Using this wave centre c1 (t), we can calculate the variance of the agent positions as
an indicator for the width of the wave and therefore for the dropout. In Figure 6(a)
this variance is compared to the variance of the travelling wave solution found by
Keller and Segel, which is σm f = (π µ )2 /3. Initially the measured variance increases
linearly towards the theoretical value, which is caused by the start of the agents on
the boundary x = 0. After the wave is fully developed, the variance starts to rise
over the theoretical value, which indicates a significantly wider wave and therefore
dropout of agents.
With these statistical values for the agent set we have now different options to
define an agent as dropped out from the wave and therefore to define the index set
Γ . The first option is to allow an agent to have a certain distance r from the centre
of the wave. Agents with a distance bigger than r are therefore considered to be
20 Benjamin Franz and Radek Erban

dropped out. Hence,

Γ1 ≡ Γ1 (t; r) = {i ∈ {1, . . . , N} | xi (t) ≥ c1 (t) − r} . (23)

Because of the non-finite support of the travelling wave solution for the original
Keller-Segel system, the measure defined in (23) is strongly dependent on r, which
makes the choice of r important. One should choose r in a way that the solution
of the original Keller-Segel model only predicts a very small number of dropout
agents. One way to pick r is to use a multiple of the theoretical standard deviation
of the wave.
A second option of defining an agent as dropped out is to use the observation that
S is exploited behind the wave. An agent is then considered to be dropped out of the
wave if the value of S at its current position is 0. Thus,

Γ2 ≡ Γ2 (t) = {i ∈ {1, . . . , N} | S(xi (t)) = 0} . (24)

Using the sets Γ1 and Γ2 we can now define 2 dropout measures d1 (t; r) and d2 (t) by

1 1
d1 (t; r) = 1 − |Γ1 (t; r)| , and d2 (t) = 1 − |Γ2 (t)| . (25)
N N
Figure 6(b) shows plots of the behaviour of d1 (t; r) and d2 (t). We can see that after
the initial period of adjustment due to the start on the boundary x = 0, all measures
have an increasing trend with some fluctuations around it. The measure d1 (t; 0.15)
matches well with d2 (t), but has less fluctuations.

−3
x 10 0.015
2.5

2
0.01
1.5
d(t)
σ(t)

1
0.005
0.5

0 0
0 0.1 0.2 0.3 0.4 0.5 0.2 0.3 0.4 0.5
t t
(a) (b)

Fig. 6 Simulation results of the variance and dropout for short times, where the same parameter
values as in Figure 5 are used.
(a) Variance σ (t) estimated from the simulation (solid line) and variance of the stationary wave
given by the mean-field model σm f (dashed line).
(b) Dropout given by (25): d1 (t; 0.1) (dash-dotted line), d1 (t; 0.15) (solid line), d1 (t; 0.2) (dotted
line) and d2 (t) (dashed line).
Hybrid modelling of individual movement and collective behaviour 21

Large time behaviour

In this section we investigate the large time behaviour of the travelling wave in the
hybrid chemotaxis Model III. We study the behaviour of the bacteria and the signal
in the half-line [0, ∞]. For large times the definitions c1 (t) and c2 (t) given by (21)–
(22) differ significantly because many agents drop behind the wave. Therefore, c2 (t)
is more meaningful to describe the centre of the wave in this case. However, as c2 (t)
depends on Γ , we can no longer use Γ ≡ Γ1 to find the agents that have dropped
out, because Γ1 depends on the definition of the centre of the wave. We therefore
use d2 (t) given by (25) as measure for the dropout in the analysis of large time
behaviour, where we are particularly interested in the slowing down of the wave.
Hence, we define the velocity of the wave v(t) through

c2 (t + ∆ T ) − c2 (t)
v(t) = , (26)
∆T
where ∆ T is chosen to be much larger than ∆ t in order to minimise the fluctuations
in v(t). We simulate N = 104 agents with the same parameters as before. The results
of one simulation are shown in Figure 7. We see that after t = 50 about 40% of the
agents have dropped out from the wave. The predicted slowing down of the wave
is demonstrated in Figure 7(b), where we plot v(t) as a function of time. We use
∆ T = 0.1 in the definition (26). As the velocity shrinks with the number of agents
in the wave, we have v(t) ≈ 1 − d2 (t), which is also demonstrated in Figure 7(b).

0.5
1 v(t)
1−d (t)
2
0.4
0.9
v(t), 1−d (t)

0.3
2

0.8
d (t)
2

0.2
0.7

0.1 0.6

0 0.5
0 10 20 30 40 50 0 10 20 30 40 50
t t
(a) (b)

Fig. 7 Dropout and velocity of the travelling wave for large time, where the same parameter values
as in Figure 5 are used.
(a) Dropout d2 (t) given by (25).
(b) Velocity of the wave v(t) given by (26) (solid line) compared with 1 − d2 (t) (circles).

Dropout in dependence on N and ∆ x

In the next step we use the derived measure (25) in order to analyse the influence
of certain system parameters on the dropout. In particular, we are interested in the
22 Benjamin Franz and Radek Erban

dependence of the dropout on the number of agents N and the gridsize ∆ x. The
variation of the number of agents N in the system is a way of comparing the hybrid
with the continuum model. One would expect that the dropout goes to 0 as N goes
to infinity. On the other hand the ∆ x dependence is a problem of the hybrid model,
as one would ideally want the dropout to be independent of the chosen grid. We
performed a number of simulations for various values of N (200 simulations for
each value) and ∆ x (100 simulations for each value) and in each case measured the
value of d1 (0.5; 0.15). The results are plotted in Figure 8. In Figure 8(a) we plot

−1
10
d (0.5; 0.15)

d (0.5; 0.15)
−2
10
1

1
−2
10

3 4 −3 −2
10 10 10 10
N ∆x

(a) (b)

Fig. 8 Dropout d1 (0.5; 0.15) given by (25) as a function of (a) N and (b) ∆ x. In each figure we
show results given by individual simulations (crosses), average values of d1 (0.5; 0.15) estimated
from simulations (circles) and linear fits explained in the text (dashed-line).

the average values of d1 (0.5; 0.15) estimated from simulations as circles. The best
linear fit in the double logarithmic plot,
√ shown as the dashed line, has a gradient of
−0.53, which indicates that d1 ∼ 1/ N. This relationship can be explained through
the central
√ limit theorem, which predicts that the noise in the system should shrink
with N.
The plot in Figure 8(b) shows a more complicated dependence. For larger values
of ∆ x the dashed line with gradient −1 can be fitted indicating that a finer grid
leads to an increase in dropout, which seems slightly surprising at first glance, as
one expects a finer grid to allow for a more accurate representation of the original
PDE. This effect can, however, be explained by the lower number of agents per
gridpoint and therefore the higher noise expected at each gridpoint. As ∆ x decreases
the dropout seems to level off, meaning that the choice of a finer grid at this point
does not influence the dropout drastically. Bearing in mind that we ideally wanted
the dropout to be independent of ∆ x, this levelling off effect seems to indicate the
region of choice for ∆ x in order to get an accurate solution.

Theoretical analysis

More theoretical insight into the dropout effect can be obtained by considering a
simplified system, where the concentration of extracellular material S is a given
Hybrid modelling of individual movement and collective behaviour 23

function that does not change over time. A natural choice for the function S(x) is
the travelling wave solution found by Keller and Segel [32]. Using the knowledge
of the exploited region behind the wave, we can adjust this function slightly to allow
for the analysis of the dropout effect. We therefore define S to be equal to 0 for x
smaller than some critical position xc and to take the form of the travelling wave
solution everywhere else. In this section we will use κ = 2, thus we put

(1 + exp(−x/µ ))−1 ,

x ≥ xc ,
S(x) = (27)
0, x < xc .

To be able to use a time-independent function for S we have to make adjustments to

the movement of the agents, as they would otherwise follow the increasing gradient
towards the right of the real axis. Therefore, we subtract the expected wave speed of
1 from the movement velocity of the agents in order to keep them in a position that
is realistic for the travelling wave. In other words, we use a coordinate system that
moves with the travelling wave solution. For example, for an agent of Model I the
evolution equation becomes
 
2 dS(xi )
dxi = µ − 1 dt + 2µ dW .
p
S(xi ) dx

With the help of this simplified system we can now make further analytic and sim-
ulative investigations into the effect of different xc on the quantity of the dropout. If
an agent enters the exploited region x < xc , two behaviours are considered. In the
first case, the agent would be considered dropped out and is absorbed by the bound-
ary, so that it has no chance of becoming part of the wave again. The second case
allows the agent to randomly move around in the exploited area and therefore al-
lows the agent to enter the non-exploited region again. For both cases we performed
100 simulations for each of the considered values of xc and measured the value of
d1 (0.5; 0.15) as defined before, this time using 0 as the mean position. The average
values of d1 (0.5; 0.15) estimated from the simulations are shown in Figure 9 as cir-
cles. To analyse the case of an absorbing boundary at x = xc we consider the system
in the limit N → ∞, which is described by the following equation (compare to (16))

∂n ∂n ∂ ∂n
 
2 dS
− =µ −n . (28)
∂t ∂x ∂x ∂x S dx

The boundary condition on the left-hand boundary can be written in the form
n(xc ) = 0. Further conditions for x → ∞ can be introduced. We look for a separable
solution of the form
n(x,t) = exp(−λ t)M(x) ,
where λ is a positive constant. Plugging this ansatz into (28) leads to
 ′ ′
S
µ M + M − 2µ M
′′ ′
+λM = 0, (29)
S
24 Benjamin Franz and Radek Erban
0 0
10 10

d (0.5; 0.15)

d (0.5; 0.15)
−1 −1
10 10

1
−2 −2
10 10

−0.12 −0.1 −0.08 −0.06 −0.12 −0.1 −0.08 −0.06

xc xc

(a) (b)

Fig. 9 Dropout d1 (0.5; 0.15) defined by (25) as a function of xc for static signal given by (27)
where the same parameter values as in Figure 5 are used. In each figure we show average values of
d1 (0.5; 0.15) estimated from 100 simulations (circles).
(a) Simulations where no comeback from x = xc is allowed. The dashed line is a result of the
theoretic analysis given by (30).
(b) Simulations where dropout agents can return. The dashed line is 50% of the dropout predicted
by (30).

where primes denote derivatives with respect to x. For the ODE (29) a non-negative
solution is sought that satisfies M(xc ) = 0 and M(x) → 0 as x → ∞. The general
solution for (29) is
γ 1+γ
   
2λ µ exp x 3+
2µ + (1 + γ + 2λ µ ) exp x 2µ
M(x) = C1    2
exp µx + 1
γ 1−γ
   
2λ µ exp x 3−
2µ + (1 − γ − 2λ µ ) exp x 2µ
− C2    2 ,
exp µx + 1

where γ = 4λ µ + 1. The integration constants C1 and C2 have to be chosen to

p

satisfy the boundary conditions. Because of the nature of (29) as an eigenvalue

problem, only the quotient C1 /C2 can be determined uniquely, which also means
that the condition M(x) → 0 as x → ∞ is satisfied for all values C1 ,C2 ∈ R. Taking a
closer look at the limit x → ∞, we can see that the direction of the approach changes
in dependence of λ , in particular, a non-negative solution can only be obtained for
λ smaller than a critical value λc (xc ). This critical value λc (xc ) is achieved for the
case where C1 /C2 turns out to be 0. Applying the left-hand boundary condition
M(xc ) = 0 for this case yields to the unique value λc (xc ) given through
   −2
1 xc xc
λc (xc ) = − exp − 1 + exp − = −S′ (xc ) .
µ µ µ
Hybrid modelling of individual movement and collective behaviour 25

This value λc (xc ) can now be used to get a predicted value of the dropout d pred (xc ,t)
via
d pred (xc ,t) = 1 − exp(λ (xc )t) . (30)
The function d pred (xc , 0.5) is plotted as the dashed line in Figure 9(a). We can see
that it matches well with the simulation results. The slight overestimation given by
d pred (xc , 0.5) can be explained through the time it takes before the first agents start
reaching the critical position xc from the starting position at x = 0. For the situation
with comeback, we choose a value λ = αλc (xc ) to predict the dropout, where α
is a constant. Matching this with the simulation results as shown in Figure 9(b) we
found that α ≈ 0.5, which indicates that about 50% of agents come back into the
wave after they have dropped out. This effect could be modelled by using a reactive
boundary [14] instead of the free diffusion zone behind the wave.

6 Conclusion

In this chapter we reviewed the advances that have been made in the field of hy-
brid modelling of collective behaviour. Hybrid models combine agent-based models
with mean-field concentration models and allow a more accurate description of cer-
tain systems than the general mean-field approach. Compared to purely agent-based
models hybrid models have the advantage of a reduced computational complexity
and a wider range of applicability. As hybrid models explicitly consider individual
behaviour as well as interactions between individuals, stochastic effects are incor-
porated which can alter the behaviour from that of the corresponding continuum
model. This became especially clear during the studies of hybrid chemotaxis mod-
els in Section 5. We showed that the hybrid models do not produce a travelling wave
in the classical sense, as agents are dropping out behind the wave. This effect leads
to a decrease in the number of agents in the wave, which also slows down the wave,
as demonstrated in Figure 7. We also discussed some of the problems and difficul-
ties related to the use of hybrid models. In particular the spatial matching between
the discrete agents and the continuous variables has to be considered. We showed in
Figure 8 that the choice of the gridsize can have a significant effect on the behaviour
of hybrid models and has to be handled with care.

Acknowledgements

The research leading to these results has received funding from the European Re-
search Council under the European Community’s Seventh Framework Programme
(FP7/2007-2013) / ERC grant agreement No. 239870. This publication was based
on work supported in part by Award No KUK-C1-013-04, made by King Abdul-
26 Benjamin Franz and Radek Erban

lah University of Science and Technology (KAUST). RE would also like to thank
Somerville College, University of Oxford, for a Fulford Junior Research Fellowship.

References

1. Adler, J.: Chemotaxis in Bacteria. Science 153, 708–716 (1966)

2. Adler, J.: Chemoreceptors in Bacteria. Science 166, 1588–1597 (1969)
3. Alarcón, T., Byrne, H., Maini, P.: A cellular automaton model for tumour growth in inhomo-
geneous environment. Journal of Theoretical Biology 225, 257–274 (2003)
4. Anderson, A.: A hybrid mathematical model of solid tumour invasion: the importance of cell
adhesion. Mathematical Medicine and Biology : a Journal of the IMA 22, 163–186 (2005)
5. Barkai, N., Leibler, S.: Bacterial chemotaxis - united we sense ... Nature 393, 18–21 (1998)
6. Berg, H.: How bacteria swim. Scientific American 233, 36–44 (1975)
7. Berg, H., Brown, D.: Chemotaxis in Esterichia coli analysed by three-dimensional tracking.
Nature 239, 500–504 (1972)
8. Bobashev, G., Goedecke, M., Yu, F., Epstein, J.: A Hybrid Epidemic Model: Combining the
Advantages of Agent-Based and Equation-Based Approaches. In: S. Henderson (ed.) Pro-
ceedings of the 2007 Winter Simulation Conference, pp. 1532–1537. IEEE (2007)
9. Buhl, J., Sumpter, D., Couzin, I., Hale, J., Despland, E., Miller, E., Simpson, S.: From disorder
to order in marching locusts. Science 312, 1402–1406 (2006)
10. Couzin, I., Krause, J., Franks, N., Levin, S.: Effective leadership and decision-making in ani-
mal groups on the move. Nature 433, 513–516 (2005)
11. Couzin, I., Krause, J., James, R., Ruxton, G., Franks, N.: Collective memory and spatial sort-
ing in animal groups. Journal of Theoretical Biology 218(1), 1–11 (2002)
12. Dallon, J., Othmer, H.: A discrete cell model with adaptive signalling for aggregation of dic-
tyostelium discoideum. Philosophical Transactions of the Royal Society B: Biological Sci-
ences 352(1351), 391–417 (1997)
13. Erban, R.: From individual to collective behaviour in biological systems. Ph.D. thesis, Uni-
versity of Minnesota (2005)
14. Erban, R., Chapman, S.J.: Reactive boundary conditions for stochastic simulations of reaction-
diffusion processes. Physical Biology 4(1), 16–28 (2007)
15. Erban, R., Chapman, S.J.: Time scale of random sequential adsorption. Physical Review E
75(4), 041,116 (2007)
16. Erban, R., Kevrekidis, I., Othmer, H.: An equation-free computational approach for extracting
population-level behavior from individual-based models of biological dispersal. Physica D
215(1), 1–24 (2006)
17. Erban, R., Othmer, H.: From individual to collective behaviour in bacterial chemotaxis. SIAM
Journal on Applied Mathematics 65(2), 361–391 (2004)
18. Erban, R., Othmer, H.: From signal transduction to spatial pattern formation in E. coli: A
paradigm for multi-scale modeling in biology. Multiscale Modeling and Simulation 3(2),
362–394 (2005)
19. Flegg, M., Chapman, J., Erban, R.: The Two Regime Method for optimizing stochastic
reaction-diffusion simulations (2011). Submitted
20. Franz, B.: Synchronisation Properties of an Agent-Based Animal Behaviour Model. Master’s
thesis, University of Oxford (2009)
21. Gerlee, P., Anderson, A.: An evolutionary hybrid cellular automaton model of solid tumour
growth. Journal of Theoretical Biology 246, 583–603 (2007)
22. Guo, Z., Sloot, P., Tay, J.C.: A hybrid agent-based approach for modeling microbiological
systems. Journal of Theoretical Biology 255(2), 163–175 (2008)
23. Helbing, D., Farkas, I., Vicsek, T.: Simulating dynamical features of escape panic. Nature
407, 487–490 (2000)
Hybrid modelling of individual movement and collective behaviour 27

24. Hellweger, F., Bucci, V.: A bunch of tiny individuals: Individual-based modeling for microbes.
Ecological Modelling 220, 8–22 (2009)
25. Heppenstall, A., Evans, A., Birkin, M.: A Hybrid Multi-Agent / Spatial Interaction Model
System for Petrol Price Setting. Transactions in GIS 9, 35–51 (2005)
26. Heppenstall, A., Evans, A., Birkin, M.: Using Hybrid Agent-Based Systems to Model
Spatially-Influenced Retail Markets. Journal of Artificial Societies and Social Simulation
9, 2 (2006)
27. Heppenstall, A., Evans, A., Birkin, M.: Genetic algorithm optimisation of an agent-based
model for simulating a retail market. Environment and Planning B: Planning and Design
34, 1051–1070 (2007)
28. Horstmann, D.: From 1970 until present: the Keller-Segel model in chemotaxis and its conse-
quences i. Jahresbericht der DMV 105(3), 103–165 (2003)
29. van Kampen, N.: Stochastic Processes in Physics and Chemistry, 3rd edn. North-Holland,
Amsterdam (2007)
30. Karch, G.: Selfsimilar profiles in large time asymptotics of solutions to damped wave equa-
tions. Studia Mathematica 143, 175–197 (2000)
31. Keller, E., Segel, L.: Model for chemotaxis. Journal of Theoretical Biology 30, 225–234
(1971)
32. Keller, E., Segel, L.: Traveling bands of chemotactic bacteria: A theoretical analysis. Journal
of Theoretical Biology 30, 235–248 (1971)
33. Landsberg, J., Waring, R.: A generalised model of forest productivity using simplified con-
cepts of radiation-use efficiency, carbon balance and partitioning. Forest Ecology and Man-
agement 95, 209–228 (1997)
34. Lui, R., Wang, Z.A.: Traveling wave solutions from microscopic to macroscopic chemotaxis
models. Journal of Mathematical Biology 61(5), 739–761 (2010)
35. Murray, J.: Mathematical Biology. Springer-Verlag (2002)
36. Osborne, J., Walter, A., Kershaw, S., Mirams, G., Fletcher, A., Pathmanathan, P., Gavaghan,
D., Jensen, O., Maini, P., Byrne, H.: A hybrid approach to multi-scale modelling of cancer.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering
Sciences 368, 5013–5028 (2010)
37. Othmer, H., Dunbar, S., Alt, W.: Models of dispersal in biological systems. Journal of Math-
ematical Biology 26, 263–298 (1988)
38. P., O., Brackbill, J.: On particle-grid interpolation and calculating chemistry in particle-in-cell
methods. Journal of Computational Physics 103, 37–52 (1993)
39. Patel, A., Gawlinski, E., Lemieux, S., Gatenby, R.: A cellular automaton model of early tumor
growth and invasion. Journal of Theoretical Biology 213, 315–331 (2001)
40. Ribba, B., Alarcón, T., Marron, K., Maini, P., Agur, Z.: The use of hybrid cellular automaton
models for improving cancer therapy. In: P. Sloot, B. Chopard, A. Hoekstra (eds.) ACRI 2004,
Lecture Notes in Computer Science, vol. 3305, pp. 444–453. Springer-Verlag (2004)
41. Smallbone, K., Gatenby, R., Gillies, R., Maini, P., Gavaghan, D.: Metabolic changes during
carcinogenesis: Potential impact on invasiveness. Journal of Theoretical Biology 244, 703–
713 (2007)
42. Spiro, P., Parkinson, J., Othmer, H.: A model of excitation and adaptation in bacterial chemo-
taxis. Proceedings of the National Academy of Sciences USA 94, 7263–7268 (1997)
43. Stevens, A.: The derivation of chemotaxis equations as limit dynamics of moderately inter-
acting stochastic many-particle systems. SIAM Journal on Applied Mathematics 61, 183–212
(2000)
44. Sumpter, D.: The principles of collective animal behaviour. Philosophical Transactions of The
Royal Society B: Biological Sciences 361, 5–22 (2006)
45. Sumpter, D.: Collective Animal Behavior. Princeton University Press (2010)
46. Wand, M., Jones, M.: Kernel smoothing. Chapman & Hall/CRC (1994)
47. Wiener, N.: Extrapolation, Interpolation, and Smoothing of Stationary Time Series. MIT Press
(1964)
48. Wooldridge, M.: An introduction to multi-agent systems. Wiley (2002)
28 Benjamin Franz and Radek Erban

49. Xue, C., Budrana, E., Othmer, H.: A hybrid model that explains radial and spiral stream for-
mation in Proteus mirabilis colonies (2011). Submitted
50. Xue, C., Hwang, H.J., Painter, K., Erban, R.: Travelling Waves in Hyperbolic Chemotaxis
Equations. Bulletin of Mathematical Biology 1, 1–29 (2010)
51. Xue, C., Othmer, H.: Multiscale models of taxis-driven patterning in bacterial populations.
SIAM Journal on Applied Mathematics 70(1), 133–167 (2009)