arXiv:0704.0271v2  [cond-mat.stat-mech]  27 Feb 2008
An individual based model with global competition interaction: fluctuation effects in
pattern formation
E. Brigatti ⋆±, V. Schwämmle† and Minos A. Neto⋆
†Centro Brasileiro de Pesquisas F́ısicas, Rua Dr. Xavier Sigaud 150, 22290-180, Rio de Janeiro, RJ, Brazil
⋆Instituto de F́ısica, Universidade Federal Fluminense, Campus da Praia Vermelha, 24210-340, Niterói, RJ, Brazil
±e-mail address: edgardo@if.uff.br
(October 26, 2018)
We present some numerical results obtained from a simple individual based model that describes
clustering of organisms caused by competition. Our aim is to show that, even when a deterministic
description developed for continuum models predicts no pattern formation, an individual based
model displays well defined patterns, as a consequence of fluctuation effects caused by the discrete
nature of the interacting agents.
87.17.Aa,87.23.Kg, 87.23.-n, 05.10.Ln
I. INTRODUCTION
Birth and death processes are two of the most relevant
characteristics of the dynamics of biological populations
and can be responsible for the emergence of stable spa-
tial patterns [1]. In fact, the intrinsic asymmetry in the
nature of birth and death processes can enhance small
initial differences in the spatial population density and
lead to the formation of structures [2–4]. These clusters
are resistant to some levels of diffusion and emerge as
soon as the birth of new individuals outcompetes their
movements. For this reason, simplified models combin-
ing birth and death processes with Brownian movement
are able to describe aggregation of individuals [5].
Another central ingredient, present in ecological sys-
tems, that can cause the generation of spatial structures,
is the competition for resources [6–8]. Different individ-
uals struggle for nutrients with a competition strength
directly dependent on the individuals’ spatial density
within the competing range. Reproduction and/or death
rates depending on the number of individuals in the sur-
roundings can stand for this kind of interaction. This fea-
ture has attracted the interest of experts from a variety of
fields, ranging from pure mathematics [9] and non-linear
physics [10–14,27] to population biology [16,17] and the-
oretical studies in evolutionary theory [20–25]. In ad-
dition, similar behaviors can be found in physical sys-
tems, such as, for example, in mode interaction in crys-
tallization fronts [18] and in spin-wave patterns [19]. It
is remarkable that the structured state generated by this
kind of frequency-dependent interaction exists only for
some specific form of the interaction [24] and is reached
through a transition in the parameter space. This transi-
tion, (segregation transition [10]), drives the steady state
of the system from a spatially homogeneous distribution
to one marked by some clearly distinguishable inhomo-
geneities.
All these models, characterized by diffusion effects and
an implementation of frequency dependent birth and
death processes, permit multiple interpretations.
In a common interpretation the system space directly
represents the physical space where the organisms live
and the diffusion represents their spatial movement.
Competition between individuals corresponds to a mech-
anism of growth control caused by limited common re-
sources. In this case, pattern formation can reproduce
the evolution of bacterial colonies [6], plankton concen-
tration [5], development of vegetation [8] or spatial dis-
tribution of predators [7].
On the other hand, a different interpretation enables
us to describe the speciation process: the generation of
two different species starting from one single continuous
population of interbreeding organisms. To be specific,
we can describe the speciation process by representing
all the phenotypic characteristics that determine the bi-
ological success of an individual by a number, the strat-
egy value, that labels each individual. By reproduction,
that includes a mutation process, an offspring inherits a
strategy that slightly differs from that of its parent. In
order to model natural selection, a frequency-dependent
mechanism that mimics competition, completes the in-
gredients necessary for the emergence of population clus-
tering. In this scenario, the generation of a new cluster
is interpreted, in a broad sense, as a speciation event.
Now, if the model space represents the mentioned strat-
egy space and the diffusion models the mutation process
during reproduction, we can identify the mechanism of
growth control with natural selection and the branching
events with the speciation process. This different inter-
pretation justifies the analogy between a model that de-
scribes the speciation process and the ones that describe
spatial pattern formation in the evolution of bacterial
colonies, vegetation or predation.
We must remember that, since this model does not in-
clude sexual reproduction, we are describing trait diver-
gence in an asexual population, rather than speciation.
Anyway, apart from effects strictly related to sexual re-
production, the dynamics characterized by the individu-
http://arxiv.org/abs/0704.0271v2
als’ diffusion from regions of low viability in favor of more
viable ones, is the essential core of these two phenomena.
A detailed and motivated discussion of these processes
can be found in the following references [20–25].
We can describe such processes starting from an indi-
vidual based model, that yields information on the be-
havior of a finite system (finite population) and accounts
for fluctuation effects caused by the discrete nature of the
interacting agents. Another approach, that neglects fluc-
tuations, describes individuals just with the use of a field
that represents the population density at each position in
space over time. This method, usually denominated con-
tinuous mean-field description [15,27], becomes exact in
the infinite-size limit if fluctuations are small compared
to averages [4]. Note that we are not making a mean-
field approximation in the nature of the interaction (see,
for instance, ref. [26]), that, in our model, is local and
different for each individual.
Choosing this second strategy, the generalization
of a well-investigated equation (Fisher-Kolmogoroff-
Petrovsky-Piscounoff [28]) is quite common at present.
In addition to a diffusion process with coefficient D and
a population growth mechanism (rate a), this equation
incorporates a growth-limiting process controlled by the
parameter b [11,12,14]:
∂ρ(x, t)
= D∇2ρ(x, t) + aρ(x, t) − (1.1)
b ρ(x, t)
ρ(y, t)F (x, y)dy ,
where ρ is the density of individuals at position x and
time t. Competition is obtained by varying the death
probability for each individual, and is controlled through
the influence function F (x, y). Let us focus on the shape
of the influence function: it can range from a simple box-
like function to a globally uniform interaction. However,
the Gaussian function should be considered particularly
relevant. If, for instance, we need to represent the ac-
tivity of a sedentary animal the interaction represented
in the influence function should take into account the in-
dividual’s daily excursion around the fixed breeding site
that can be represented by Brownian motion and, for this
reason, by a Gaussian distribution. In the same way, if
we want to represent the habitat degeneration induced by
the growth of a colony of plants [29] we can think that the
colony is originated by a single individual that disperses
its seeds in a way also well described by Brownian motion.
More generally, for a biological interaction that does not
stop at some defined length (presence of a cutoff), and
that is nonlocal and controlled by a purely stochastic
process, the Gaussian function should be the most natu-
ral choice. On the other hand, this choice is a source of
complications. Deterministic descriptions, in the case of
a Gaussian influence function, predict no pattern forma-
tion [12,30]. However, such descriptions do not take into
account fluctuations arising from the discrete character
of individuals. The importance of these fluctuations has
been recently pointed out in a quite paradigmatic exam-
ple, where random walking organisms that reproduce and
die at a constant rate spontaneously aggregate [2–5,31].
The deterministic approximation is not able to show
this behavior, incapable of capturing the essential asym-
metry between birth, a multiplicative process that incre-
ments the density in the regions adjacent to the parent,
and death events, that occur anywhere. Even when
the patterns can be obtained within the deterministic
description, a recent work outlines the importance of
fluctuations by showing their impact on affecting transi-
tion points and amplitudes [14,32,33].
In this work we present some numerical results ob-
tained by means of a simple individual based model. Our
aim is to show the appearance of a segregation transition
in a model where the deterministic instability, produced
by the non-local interaction, is not sufficient for gener-
ating inhomogeneities [12], but the superimposed micro-
scopic stochastic fluctuations permit the emergence of
patterns. Moreover, we compare the model implemen-
tation in the strategy space with the implementation in
the physical space. In the first, used to characterize the
speciation process, diffusion corresponds to a mutation
phenomenon, operating just one time in each individ-
ual’s life. In the second, directly related to the reaction-
diffusion equation (eq. (1.1)), diffusion describes a typical
Brownian motion.
The paper is organized as follows. The next section
describes the model used in our simulations. Section III
shows, for specific values of the parameters, the emer-
gence of spatially inhomogeneous steady states. In sec-
tion IV we prove that these patterns are not caused by
some finite size effect. Section V is devoted to illustrate
the segregation transition and general conditions that al-
low spatial segregation of arbitrary wavelengths. In sec-
tion VI we describe, in the light of the existing literature,
the cluster size dependence on diffusion rate and popu-
lation size and give some hints related to the behavior of
fluctuations. Conclusion are reported in section VII.
II. THE MODEL
The simulations start with an initial population of N0
individuals located along a ring of length L, i.e. we take
periodic boundary conditions. At each time step, our
model is controlled by the following microscopic rules:
1) Each individual, characterized by its position x, dies
with probability P ,
P = K ·
exp(−
(x− yj)2
) , (2.1)
whereN(τ) is the total number of individuals at the ac-
tual step τ and yj is the respective individuals’ positions.
The distance between two individuals is obtained by tak-
ing the shorter distance on the ring. The strength of
competition declines with increasing distance according
to a Gaussian function with deviation C. The parameter
K depicts the carrying capacity.
2) If the individual survives this death selection, it
reproduces. The newborn, starting from the parent’s lo-
cation, moves in a random direction a distance obtained
from a Gaussian distribution of standard deviation σ.
This change represents the effect of mutations in the
offspring phenotype.
As soon as all the individuals have passed the death
selection and eventually reproduced, the next time step
begins. This model implementation is analogous to the
diffusion process described by eq. (1.1).
To establish a more direct comparison between that
mean field description and the individual based simula-
tion, we have also implemented the model with an exact
microscopic representation of the diffusion term. In this
case, at any given time step, we perform first a loop over
all particles where individuals move some distance, in a
random direction, chosen from a Gaussian distribution of
standard deviation σ. At the end of this loop, a second
one starts, where:
1) Each individual with strategy x dies with a proba-
bility obtained from eq. (2.1).
2) If the individual survives the death selection pro-
cess, it reproduces and the newborn maintains the same
location of the parent.
If, in eq. (1.1), we measure time in units of the simula-
tion time step, the coefficient D is related to our simula-
tion parameter through D ∝ σ2. The influence function
is given by a Gaussian of standard deviation C and the
effective growth rate is 1− P , with P given by eq. (2.1).
This non constant growth rate can be represented by
eq. (1.1) for a = 1 and the frequency dependent part
included in the integral term (see ref. [14]).
Essentially, the only difference between these two ver-
sions of our model is that, in the first one, individuals
move only at birth, while in the second version they can
move at every time step throughout their life. Since the
death probability, at equilibrium, is approximately 1/2,
usually, in the second version, an individual will move
between one and two times during its entire life. For
this reason, there should be no relevant differences in the
qualitative behavior of the two model implementations.
As shown by the measures reported in section IV, the
only significant effect is the appearance of slightly wider
distributions.
III. MODULATION
In the following we present some typical examples of
steady states generated by the dynamics of the model
that clearly show the emergence of patterns for some spe-
cific values of the parameters.
For a global competition that results to be extremely
long-ranged (large C values, in relation to the values of
parameters σ and K), the steady state is characterized
by a spatially homogeneous occupancy. If the σ value is
sufficiently large or/and the K value sufficiently small,
totally homogeneous distributions are obtained (Figure
1), otherwise the solution is smooth but with the pop-
ulation concentrated in one region of the ring, with its
width controlled by σ.
As stated above, a simple heuristic analysis of eq. (1.1)
in the Fourier space shows that there is a necessary con-
dition for the emergence of inhomogeneity: the Fourier
transform of the influence function must have negative
values and large enough magnitude [12,13]. A Gaus-
sian in an infinite domain has a positive counterpart in
Fourier space and so does not match such requirements.
In contrast with these results, the fluctuations present in
our individual based model arrive to excite one specific
mode and modulations of this wavelength appear. The
tuning of the parameters allows modulations of arbitrary
wavelengths (Figure 1).
When C is decreased, the competition between modes
becomes stronger and no single mode dominates. In this
situation, some small regions of the ring are occupied
forcing all the remaining areas, up to some range, to
be nearly empty. The landscape results populated by
several living colonies divided by dead regions. There is
almost no competition between individuals of different
colonies and the space separating them can be identified
with an effective interaction length. This steady state
(spiky state [10]) corresponds to a sequence of isolated
colonies (spikes) and seen in the Fourier space, many
active wavelengths contribute to it (Figure 2).
Finally, for extremely short-ranged competition, in re-
lation to the σ value, no collective cooperation between
different excited modes emerges and a noisy spatially
homogeneous distribution appears (Figure 2).
We describe these paradigmatic steady states of the
system by characterizing the related spatial structure
with the help of a structure function S(q) [14] defined
as follows:
S(q) =
exp[i2πq · xj(τ)]
(3.1)
where the sum is performed over all individuals j with
their positions determined by xj(τ). Note, that, for con-
venience, the structure function is calculated only over
the closed interval [0, L]. The function is averaged over
some time interval T in order to avoid noisy data. S(0)
corresponds to the square of the mean number of individ-
uals in the system. The maxima of S(q) identify the rel-
evant periodicity present in the steady state. We will see
that the position of the global maximum (qM ) provides
an appropriate order parameter for the identification of
the segregation transition.
In our study we explored two different initial condi-
tions. In the first (local i.c.), the colony is located in
a finite and compact region of the space. In the second
(global i.c.), the individuals are spread all over the space.
The final distribution is independent of this choice and,
generally, local initial conditions make the system reach
the steady state more slowly. For this reason, if not dif-
ferently specified, our results are obtained from global
initial conditions.
IV. FINITE SIZE EFFECTS
Our analysis starts by exploring the model dependence
on the space size L. The reason for such interest is given
by the necessity of testing if the pattern formation is not
merely a product of some finite size effect. This is im-
portant, in the light of what was reported by Fuentes
et al. in ref. [12]. In their work, a numerical solution
of eq. (1.1) with a Gaussian influence function, showed a
segregation transition. But such a transition was just the
effect of the finite domain size that acted like a cut-off
for the Gaussian. Evidence of this interpretation came
from the observation that the amplitude of the pattern
depended on the ratio of the standard deviation of the in-
fluence function to the domain size - the critical values of
the standard deviation corresponding to the segregation
transition depended linearly on the domain size - and the
same patterns appeared for a modified Gaussian, which
vanishes abruptly beyond a cutoff.
The study of our individual based model gave differ-
ent results. By running some simulations with exactly
the same parameters but changing the ring extension, we
were able to show that the system is not influenced by the
domain size. If we choose data from spiky steady states,
that permit clear quantitative measures, it is possible to
remark that the general morphology of the patterns does
not change increasing the L value. In fact, both the pop-
ulation density and the mean number of peaks per space
interval remains constant. Moreover, in order to provide
a more precise test of possible small variations in the
distribution, we measured the cluster size. This quan-
tity was calculated by evaluating the standard deviation
(< x2 >i −< x >i2)1/2 of the position of the i individu-
als confined in each peak, then averaged over the different
peaks present at step τ and, finally, averaged over many
time steps after the system has reached the steady state.
Varying the system size caused no changes in the clus-
ters size (see Figure 3). From this result, we concluded
that the general aspect of the steady state does not
change with L. In particular, in contrast with what hap-
pens when the mean field equation is solved numerically,
the patterns do not depend on the ratio C
. For example,
for C = 0.2 and L = 50 we obtained a spiky steady state,
for C = 0.004 and L = 1 (same ratio) we obtained an
homogeneous steady state. Taking into account these re-
sults, from now on, all our simulations are implemented
on a ring of size 1.
We have just shown how the average of the population
size < N > scales with L in the steady state. Now, we
will give, through a simple heuristic argument, an esti-
mation of < N > as a function of the parameters K and
C, that will be useful also in the following of our analysis.
We can assume that, locally and in the steady state,
the number of deaths must be, on average, compensated
by the number of newborns, in order to comprise a sta-
ble population. For this reason, the death probability P
must equal 1/2. Assuming that the number of neighbors
that compete with a single individual are the ones living
up to a distance C and that, in these surroundings, the
average density N/L can be considered to be uniform, P
reduces to K ·2C ·N/L. Thus, N ∝ L/(CK). Looking at
Figure 4 we can see that this crude evaluation, that ne-
glects diffusion, inhomogeneity and reduces the influence
function to a box, describes well the general behavior of
the data obtained from our simulations.
V. SEGREGATION TRANSITION
In the previous paragraphs we supported the fact that
steady states, depending on the parameter values, can
assume inhomogeneous spatial distributions. Now, we
will try to describe the transition towards these states
(segregation transition). The structure function intro-
duced in eq. (3.1) provides a proper order parameter
to describe this transition. Different regions in the pa-
rameter space, coinciding with different steady states,
correspond to different positions of the global maximum
(obviously we are not taking into account S(q) at q = 0)
of the structure function. The transition from a homo-
geneous to an inhomogeneous distribution (see Figure 2)
matches the jump of the position of the global maxi-
mum (qM ) to a clear integer value, corresponding to the
number of clusters present in the space. For this rea-
son we can characterize the transition by looking at the
shape assumed by S(q), or looking at the value of qM .
If the space is homogeneously occupied, the structure
function does not present an integer maximum. On the
contrary, the maximum is located at qM ≃ 1.4. This
value corresponds to an uniform distribution of individ-
uals in the interval [0,1], approximated by the expression
exp(i2πq ·x)dx
. The segregation transition is char-
acterized by the passage of qM from 1.4 to an integer
value as soon as a modulation becomes dominant. In
Figure 5 we show qM as a function of C, varying K and
σ. First of all, from the analysis of these data, we can
observe that the number of clusters scales as C−1 (or,
equivalently, the periodicity of the inhomogeneous phase
has wave lengths proportional to C). Moreover, a criti-
cal value of C exists for which the transition takes place.
This Ccritic grows with 1/K and with σ. An analysis
of the available data suggests the possible dependence:
Ccritic ∝ σ2/3K−1/3 . Finally, for larger values of the
parameter C, in this range of the parameters σ and K,
the distributions are characterized by just one peak.
For any value of the competition strength, as can be
seen in Figure 6, there exists a critical value σcritic, de-
pendent on C and K, above which no spatial structures
emerge. Another measurement, that permits us to state
this relation in a different and clearer way, is presented
in the next paragraph.
VI. CLUSTER SIZE AND FLUCTUATIONS
In the following we describe, in the light of the existing
literature, the cluster size dependence on diffusion rate
and population size for the two different implementations
of the model.
We start by analyzing the typical size S of the clus-
ters that appear in the spiky phase. The data exposed in
Figure 7 show a dependence of the cluster size on the dif-
fusion coefficient: S ∝ σ, equivalent to S ∝
D. These
results are in accordance with the data presented in ref.
[33] obtained from an individual based model. In addi-
tion, this work pointed out how this behavior deviates
from the conclusions obtained from the deterministic ap-
proximation, where the cluster size was only weakly de-
pendent on the diffusion coefficient; another fact support-
ing the relevance of fluctuation effects in these systems.
We can easily interpret the dependence of the clus-
ter size on the diffusion coefficient by assuming that the
individuals confined in a cluster diffuse a distance pro-
portional to
DJ where J is the number of jumps the
individual performs in its life. In the case of the first
implementation of the model, where the diffusion is due
to the mutation process, J is obviously 1. Similar results
are obtained with the second implementation (see Fig-
ure 7), apart from a slightly wider cluster size (in this
case, individuals move, on average, more than just once
during their lifetime). Even so, the data show the same
dependence on the diffusion coefficient.
Finally, we present S as a function of K: S ∝
The reasons for this behavior are already explained in
ref. [34]: the cluster size is not controlled just by the
single individual’s number of jumps; in fact the diffusive
process continues with its descendants. For this reason,
it is proportional to the mean lifetime of a family, esti-
mated to be proportional toK−1 (see ref. [34] for details).
We introduce a new quantity that is useful for describ-
ing the existence of a critical diffusion value, giving an
estimation of its dependence on others parameters and
a confirmation of previous results. This quantity, which
we call mobility M(τ), estimates the mean mobility of
individuals. At a given time step τ , we choose an indi-
vidual i. Then we look for the closest agent, among all
the population, at time step τ − 1. We identify as di the
distance between these two individuals. Averaging over
the entire population N(τ) we obtain:
M(τ) =
di. (6.1)
The values assumed by M on varying the parameters σ
and C are shown in Figure 8. It is easy to distinguish two
clearly different behaviors. If the system is organized in a
spiky state (when σ ≤ σcritic), M(τ) ∝ σ. M is another
way of measuring the mean distance that an individual
moves during its lifetime inside the region defined by the
cluster. For this reason, this measure is coherent with the
data obtained from the direct evaluation of the standard
deviation of the clusters. In contrast, when the system is
organized in the homogeneous phase (when σ > σcritic)
M becomes independent of σ and is proportional to the
inverse of the occupation density M(τ) ∝ KC. The val-
ues of the mobility obtained from simulations with dif-
ferent values of C and K can be easily collapsed into
one function (see the inset of Figure 8). The collapse is
performed using the scaling σ → σ/C
K. This indicates
that the characteristic value of the crossover, σcritic, that
separates the two different behaviors of M(τ), scales as:
σcritic ∝ C
K. (6.2)
We conclude our study with some measurements try-
ing to catch some properties of the system fluctuations.
First, we estimated the fluctuations of the total popula-
tion, averaging over different simulations. The variance
turned out to be constant throughout the time evolution
and of the order of the square root of the total population.
The mechanism of auto-regulation of the population di-
mension does not allow the growth of big differences in
the total number of individuals.
For this reason, we focused our attention on the spa-
tial distribution of the population and tried to measure
some properties of these fluctuations. We studied the
variation of the local number of individuals in the same
simulation, for different times. We analyzed the evo-
lution of the system starting from local initial condi-
tions, with the population concentrated in the interval
[0.49, 0.51]. In this situation, the system evolves in time
with a small cluster fluctuating around the initial space
interval. This situation changes when a branching event
occurs that generates two well-defined clusters. Our in-
terest is in showing the behavior of local space fluctua-
tions and capturing possible variations in correspondence
with the branching event. First of all, we looked at the
mean value of the spatial local fluctuations Fs(τ), defined
as Fs(τ) = (
f2j )
1/2, where fj is the occupancy varia-
tion of the bin j from time step τ − 1 to time step τ . We
performed the average over all the b bins the ring was
divided in and obtained Fs(τ) =
N(τ), with no rel-
evant variations throughout the time evolution, even in
the time interval corresponding to the branching event.
More interesting is the shape of the frequency distribu-
tion of the size of fj . In fact, a simple Gaussian does
not fit this distribution, that presents extended tails (see
Figure 9). Throughout the system time evolution, the
shape of the normalized distribution is conserved. For
global initial conditions the same frequency distribution,
with extended tails, is recovered at the steady state. It
is identical to the one obtained with local initial condi-
tions and measured at the steady state. We think that
the deviation of the distribution from a Gaussian can be
considered as a hint that fluctuations play a relevant role
in the dynamics of the systems.
VII. CONCLUSIONS
We presented some results regarding clustering of or-
ganisms caused by a frequency-dependent interaction
that represents competition. We showed how this way
of modeling competition can be used not only to de-
scribe spatial phenomena in population biology, but also,
through a more abstract interpretation, to test ideas of
evolutionary theory (for example, studying the speciation
process).
From this unifying perspective, our study, obtained
from an extensive collection of data coming from sim-
ulations of an individual based model with global com-
petition, pointed out the relevance of fluctuation effects
in pattern formation. For the influence function adopted,
the mean field description predicts the absence of spatial
structures. On the contrary, fluctuations are able to ex-
cite the emergence of well defined patterns, which can
not be generated from a deterministic instability.
Furthermore, we discussed other fundamental proper-
ties of our model in the light of the existing literature,
unfolding a comparison with other models that describe
spatial segregation originated by some deterministic in-
stability. We showed that the observed patterns are not
due to a finite size effect, we characterized the behavior
of the segregation transition in various regions of the pa-
rameter space and we studied the existence of a critical
diffusion value. We analyzed the dependence of the clus-
ter size on the diffusion coefficient and pointed out some
characteristics of the fluctuations of the system.
ACKNOWLEDGMENTS
We are grateful to J.S. Sá Martins for a critical reading
of the manuscript and thank the Brazilian agency CNPq
for financial support.
[1] S. A. Levin and L. A. Segel, SIAM Rev. 27, 45 (1985); R.
Durrett and S. A. Levin, Philos. Trans. R. Soc. London,
Ser. B 343, 329 (1994).
[2] Y.-C. Zhang, M. Serva and M. Polikarpov, J. Stat. Phys.
58, 849 (1990).
[3] N. M. Shnerb, Y. Louzoun, E. Bettelheim and S.
Solomon, Proc. Natl. Acad. Sci. U.S.A. 97, 10322 (2000).
[4] B. Houchmandzadeh, Phys. Rev. E 66, 052902 (2002).
[5] W. R. Young, A. J. Roberts and G. Stuhne, Nature 412,
328 (2001).
[6] E. Ben-Jacob, I. Cohen, and H. Levine, Adv. Phys. 49,
395 (2000); J.Wakita, K. Komatsu, A. Nakahara, T. Mat-
suyama, and M. Matsushita, J. Phys. Soc. Jpn. 63, 1205
(1994).
[7] K. S. McCann, J. B. Rasmussen, and J. Umbanhowar,
Ecol. Lett. 8, 513 (2005).
[8] J. B. Wilson and A. D. Q. Agnew, Adv. Ecol. Res. 23,
263 (1992); R. Lefever and O. Lejeune, Bull. Math. Biol.
59, 263 (1997); N. M. Shnerb, P. Sarah, H. Lavee, and
S. Solomon, Phys. Rev. Lett. 90, 038101 (2003).
[9] N. F. Britton, SIAM J. Appl. Math. 50, 1663 (1990), S.
A. Gourley and N. F. Britton, J. Math. Biol. 34, 297
(1996).
[10] Y. E. Maruvka and N. M. Shnerb, Phys. Rev. E 73,
011903 (2006).
[11] A. Sasaki, J. Theor. Biol. 186, 415 (1997).
[12] M. A. Fuentes, M. N. Kuperman and V. M. Kenkre,
Phys. Rev. Lett. 91, 158104 (2003).
[13] M.A. Fuentes, M.N. Kuperman and V.M. Kenkre, J.
Phys. Chem. B 108, 10505 (2004).
[14] E. Hernandez-Garcia and C. Lopez, Phys. Rev. E 70,
016216 (2004).
[15] J. Billingham, Nonlinearity 17, 313 (2004).
[16] G. Flierl, D. Grunbaum, S. Levin and D. Olson, J. Theor.
Biol. 196, 397 (1999).
[17] B. M. Bolker, Theor. Popul. Biol. 64, 255 (2003).
[18] D. A. Kurtze, Phys. Rev. B 40, 11104 (1989).
[19] F.-J. Elmer, Phys. Rev. Lett. 70, 2028 (1993).
[20] J. Roughgarden, Am. Nat. 106, 683 (1972).
[21] F. Bagnoli and M. Bezzi, Phys. Rev. Lett. 79 3302
(1997).
[22] U. Dieckmann and M. Doebeli, Nature 400, 354 (1999);
M. Doebeli and U. Dieckmann, Nature 421, 259 (2003).
[23] E. Brigatti, J. S. Sá Martins and I. Roditi, Physica A
376, 378 (2007); V. Schwämmle and E. Brigatti, Euro-
phys. Lett. 75, 342 (2006).
[24] S. Pigolotti, C. Lopez and E. Hernandez-Garcia, Phys.
Rev. Lett. 98 258101 (2007).
[25] H. Sayama, M.A.M. de Aguiar, Y. Bar-Yam and M.
Baranger, Phys. Rev. E 65, 051919 (2002)
[26] H. Sayama, L. Kaufman and Y. Bar-Yam, Phys. Rev. E
62, 7065 (2000).
[27] D. R. Nelson and N.M. Shnerb, Phys. Rev. E 58, 1383
(1998).
[28] R. A. Fisher, Ann. Eugenics 7, 353 (1937); A. Kolmogo-
roff, I. Petrovsky and N. Piscounoff, Moscow Univ. Bull.
Math. 1, 1 (1937).
[29] F. J. Weissing and J. Huisman, J. Theor. Biol. 168, 323
(1994).
[30] J. Polechova and N. H. Barton, Evolution, 59, 1194
(2005).
[31] M. Meyer, S. Havlin and A. Bunde, Phys. Rev. E 54,
5567 (1996).
[32] C. Lopez and E. Hernandez-Garcia, Physica D 199, 223
(2004).
[33] E. Hernandez-Garcia and C. Lopez, Physica A 356, 95
(2005).
[34] E. Hernandez-Garcia and C. Lopez, J. Phys.: Condens.
Matter 17, S4263 (2005).
FIG. 1. Homogeneous steady state distribution ( top,
C = 4.0, 1/K = 80000, σ = 0.01) and modulated steady state
distribution (bottom, C = 0.9, 1/K = 18000, σ = 0.01). The
insets show the structure functions, S(q), of the correspond-
ing simulations. We show the distributions at time step 1000,
whereas the structure functions are averaged over 500 time
steps.
FIG. 2. Spiky (top, C = 0.059, 1/K = 500, σ = 0.001) and
homogeneous (bottom, C = 0.005, 1/K = 500, σ = 0.001)
steady states. The insets show the structure functions S(q).
We show the distributions at time step 2000, whereas the
structure functions are averaged over 1000 time steps. The
transition between these two states, in this typical range of
parameters, has been extensively studied.
0 10 20 30 40 50
0.005
FIG. 3. Dependence of the cluster size on the ring size
L. We present data from the model with mutation (trian-
gles) and from the one that implements diffusion (circles),
C = 0.2, K = 0.0029, σ = 0.001. The average is carried out
over all the clusters present at a given time step and over
different time steps.
FIG. 4. The number of individuals N present in the steady
state is proportional to (CK)−1. This result is in accordance
with the one obtained for a box-type influence function of
length C (see ref. [32]). We present data for different simu-
lations with 1/K ∈ [50, 500] and σ ∈ [0.0001, 0.01]. This last
parameter does not influence the final number of individuals.
The dashed line has slope -1.
FIG. 5. Segregation transition at Ccritical. Upper figure:
variation in dependence of K, where 1/K = 50, 150, 400, 500
and σ = 0.001. Lower figure: variation in dependence of σ,
where σ = 0.0005, 0.001, 0.002, 0.005, 0.01 and 1/K = 200.
FIG. 6. As shown in these figures, a critical σ value exists,
above which no spatial structures emerge. Upper figure: vari-
ation in dependence on C; we set 1/K = 100. Lower figure:
variation in dependence on K; we set C = 0.01.
FIG. 7. Top: Cluster size as a function of σ; 1/K = 300.
The solid line has slope 1. Bottom: cluster size as a function
of 1/K; σ = 0.0001. The solid line has slope 1/2. Trian-
gles represent data from the simulations where diffusion is
implemented through mutations, circles for the direct imple-
mentation of the diffusive process; we set C = 0.09.
FIG. 8. Mobility dependence on the diffusion parameter
σ for different C values; K = 0.01. In the inset, the data
collapse for arbitrary values of the parameters C and K.
FIG. 9. Spatial local fluctuation distribution at the steady
state: deviation from a Gaussian (C = 0.1, K = 0.00005,
σ = 0.0017). Data are averaged over 5 time steps. The con-
tinuous line is the best fit Gaussian.

  We present some numerical results obtained from a simple individual based
model that describes clustering of organisms caused by competition. Our aim is
to show how, even when a deterministic description developed for continuum
models predicts no pattern formation, an individual based model displays well
defined patterns, as a consequence of fluctuations effects caused by the
discrete nature of the interacting agents.

arXiv:0704.0271v2  [cond-mat.stat-mech]  27 Feb 2008
An individual based model with global competition interaction: fluctuation effects in
pattern formation
E. Brigatti ⋆±, V. Schwämmle† and Minos A. Neto⋆
†Centro Brasileiro de Pesquisas F́ısicas, Rua Dr. Xavier Sigaud 150, 22290-180, Rio de Janeiro, RJ, Brazil
⋆Instituto de F́ısica, Universidade Federal Fluminense, Campus da Praia Vermelha, 24210-340, Niterói, RJ, Brazil
±e-mail address: edgardo@if.uff.br
(October 26, 2018)
We present some numerical results obtained from a simple individual based model that describes
clustering of organisms caused by competition. Our aim is to show that, even when a deterministic
description developed for continuum models predicts no pattern formation, an individual based
model displays well defined patterns, as a consequence of fluctuation effects caused by the discrete
nature of the interacting agents.
87.17.Aa,87.23.Kg, 87.23.-n, 05.10.Ln
I. INTRODUCTION
Birth and death processes are two of the most relevant
characteristics of the dynamics of biological populations
and can be responsible for the emergence of stable spa-
tial patterns [1]. In fact, the intrinsic asymmetry in the
nature of birth and death processes can enhance small
initial differences in the spatial population density and
lead to the formation of structures [2–4]. These clusters
are resistant to some levels of diffusion and emerge as
soon as the birth of new individuals outcompetes their
movements. For this reason, simplified models combin-
ing birth and death processes with Brownian movement
are able to describe aggregation of individuals [5].
Another central ingredient, present in ecological sys-
tems, that can cause the generation of spatial structures,
is the competition for resources [6–8]. Different individ-
uals struggle for nutrients with a competition strength
directly dependent on the individuals’ spatial density
within the competing range. Reproduction and/or death
rates depending on the number of individuals in the sur-
roundings can stand for this kind of interaction. This fea-
ture has attracted the interest of experts from a variety of
fields, ranging from pure mathematics [9] and non-linear
physics [10–14,27] to population biology [16,17] and the-
oretical studies in evolutionary theory [20–25]. In ad-
dition, similar behaviors can be found in physical sys-
tems, such as, for example, in mode interaction in crys-
tallization fronts [18] and in spin-wave patterns [19]. It
is remarkable that the structured state generated by this
kind of frequency-dependent interaction exists only for
some specific form of the interaction [24] and is reached
through a transition in the parameter space. This transi-
tion, (segregation transition [10]), drives the steady state
of the system from a spatially homogeneous distribution
to one marked by some clearly distinguishable inhomo-
geneities.
All these models, characterized by diffusion effects and
an implementation of frequency dependent birth and
death processes, permit multiple interpretations.
In a common interpretation the system space directly
represents the physical space where the organisms live
and the diffusion represents their spatial movement.
Competition between individuals corresponds to a mech-
anism of growth control caused by limited common re-
sources. In this case, pattern formation can reproduce
the evolution of bacterial colonies [6], plankton concen-
tration [5], development of vegetation [8] or spatial dis-
tribution of predators [7].
On the other hand, a different interpretation enables
us to describe the speciation process: the generation of
two different species starting from one single continuous
population of interbreeding organisms. To be specific,
we can describe the speciation process by representing
all the phenotypic characteristics that determine the bi-
ological success of an individual by a number, the strat-
egy value, that labels each individual. By reproduction,
that includes a mutation process, an offspring inherits a
strategy that slightly differs from that of its parent. In
order to model natural selection, a frequency-dependent
mechanism that mimics competition, completes the in-
gredients necessary for the emergence of population clus-
tering. In this scenario, the generation of a new cluster
is interpreted, in a broad sense, as a speciation event.
Now, if the model space represents the mentioned strat-
egy space and the diffusion models the mutation process
during reproduction, we can identify the mechanism of
growth control with natural selection and the branching
events with the speciation process. This different inter-
pretation justifies the analogy between a model that de-
scribes the speciation process and the ones that describe
spatial pattern formation in the evolution of bacterial
colonies, vegetation or predation.
We must remember that, since this model does not in-
clude sexual reproduction, we are describing trait diver-
gence in an asexual population, rather than speciation.
Anyway, apart from effects strictly related to sexual re-
production, the dynamics characterized by the individu-
http://arxiv.org/abs/0704.0271v2
als’ diffusion from regions of low viability in favor of more
viable ones, is the essential core of these two phenomena.
A detailed and motivated discussion of these processes
can be found in the following references [20–25].
We can describe such processes starting from an indi-
vidual based model, that yields information on the be-
havior of a finite system (finite population) and accounts
for fluctuation effects caused by the discrete nature of the
interacting agents. Another approach, that neglects fluc-
tuations, describes individuals just with the use of a field
that represents the population density at each position in
space over time. This method, usually denominated con-
tinuous mean-field description [15,27], becomes exact in
the infinite-size limit if fluctuations are small compared
to averages [4]. Note that we are not making a mean-
field approximation in the nature of the interaction (see,
for instance, ref. [26]), that, in our model, is local and
different for each individual.
Choosing this second strategy, the generalization
of a well-investigated equation (Fisher-Kolmogoroff-
Petrovsky-Piscounoff [28]) is quite common at present.
In addition to a diffusion process with coefficient D and
a population growth mechanism (rate a), this equation
incorporates a growth-limiting process controlled by the
parameter b [11,12,14]:
∂ρ(x, t)
= D∇2ρ(x, t) + aρ(x, t) − (1.1)
b ρ(x, t)
ρ(y, t)F (x, y)dy ,
where ρ is the density of individuals at position x and
time t. Competition is obtained by varying the death
probability for each individual, and is controlled through
the influence function F (x, y). Let us focus on the shape
of the influence function: it can range from a simple box-
like function to a globally uniform interaction. However,
the Gaussian function should be considered particularly
relevant. If, for instance, we need to represent the ac-
tivity of a sedentary animal the interaction represented
in the influence function should take into account the in-
dividual’s daily excursion around the fixed breeding site
that can be represented by Brownian motion and, for this
reason, by a Gaussian distribution. In the same way, if
we want to represent the habitat degeneration induced by
the growth of a colony of plants [29] we can think that the
colony is originated by a single individual that disperses
its seeds in a way also well described by Brownian motion.
More generally, for a biological interaction that does not
stop at some defined length (presence of a cutoff), and
that is nonlocal and controlled by a purely stochastic
process, the Gaussian function should be the most natu-
ral choice. On the other hand, this choice is a source of
complications. Deterministic descriptions, in the case of
a Gaussian influence function, predict no pattern forma-
tion [12,30]. However, such descriptions do not take into
account fluctuations arising from the discrete character
of individuals. The importance of these fluctuations has
been recently pointed out in a quite paradigmatic exam-
ple, where random walking organisms that reproduce and
die at a constant rate spontaneously aggregate [2–5,31].
The deterministic approximation is not able to show
this behavior, incapable of capturing the essential asym-
metry between birth, a multiplicative process that incre-
ments the density in the regions adjacent to the parent,
and death events, that occur anywhere. Even when
the patterns can be obtained within the deterministic
description, a recent work outlines the importance of
fluctuations by showing their impact on affecting transi-
tion points and amplitudes [14,32,33].
In this work we present some numerical results ob-
tained by means of a simple individual based model. Our
aim is to show the appearance of a segregation transition
in a model where the deterministic instability, produced
by the non-local interaction, is not sufficient for gener-
ating inhomogeneities [12], but the superimposed micro-
scopic stochastic fluctuations permit the emergence of
patterns. Moreover, we compare the model implemen-
tation in the strategy space with the implementation in
the physical space. In the first, used to characterize the
speciation process, diffusion corresponds to a mutation
phenomenon, operating just one time in each individ-
ual’s life. In the second, directly related to the reaction-
diffusion equation (eq. (1.1)), diffusion describes a typical
Brownian motion.
The paper is organized as follows. The next section
describes the model used in our simulations. Section III
shows, for specific values of the parameters, the emer-
gence of spatially inhomogeneous steady states. In sec-
tion IV we prove that these patterns are not caused by
some finite size effect. Section V is devoted to illustrate
the segregation transition and general conditions that al-
low spatial segregation of arbitrary wavelengths. In sec-
tion VI we describe, in the light of the existing literature,
the cluster size dependence on diffusion rate and popu-
lation size and give some hints related to the behavior of
fluctuations. Conclusion are reported in section VII.
II. THE MODEL
The simulations start with an initial population of N0
individuals located along a ring of length L, i.e. we take
periodic boundary conditions. At each time step, our
model is controlled by the following microscopic rules:
1) Each individual, characterized by its position x, dies
with probability P ,
P = K ·
exp(−
(x− yj)2
) , (2.1)
whereN(τ) is the total number of individuals at the ac-
tual step τ and yj is the respective individuals’ positions.
The distance between two individuals is obtained by tak-
ing the shorter distance on the ring. The strength of
competition declines with increasing distance according
to a Gaussian function with deviation C. The parameter
K depicts the carrying capacity.
2) If the individual survives this death selection, it
reproduces. The newborn, starting from the parent’s lo-
cation, moves in a random direction a distance obtained
from a Gaussian distribution of standard deviation σ.
This change represents the effect of mutations in the
offspring phenotype.
As soon as all the individuals have passed the death
selection and eventually reproduced, the next time step
begins. This model implementation is analogous to the
diffusion process described by eq. (1.1).
To establish a more direct comparison between that
mean field description and the individual based simula-
tion, we have also implemented the model with an exact
microscopic representation of the diffusion term. In this
case, at any given time step, we perform first a loop over
all particles where individuals move some distance, in a
random direction, chosen from a Gaussian distribution of
standard deviation σ. At the end of this loop, a second
one starts, where:
1) Each individual with strategy x dies with a proba-
bility obtained from eq. (2.1).
2) If the individual survives the death selection pro-
cess, it reproduces and the newborn maintains the same
location of the parent.
If, in eq. (1.1), we measure time in units of the simula-
tion time step, the coefficient D is related to our simula-
tion parameter through D ∝ σ2. The influence function
is given by a Gaussian of standard deviation C and the
effective growth rate is 1− P , with P given by eq. (2.1).
This non constant growth rate can be represented by
eq. (1.1) for a = 1 and the frequency dependent part
included in the integral term (see ref. [14]).
Essentially, the only difference between these two ver-
sions of our model is that, in the first one, individuals
move only at birth, while in the second version they can
move at every time step throughout their life. Since the
death probability, at equilibrium, is approximately 1/2,
usually, in the second version, an individual will move
between one and two times during its entire life. For
this reason, there should be no relevant differences in the
qualitative behavior of the two model implementations.
As shown by the measures reported in section IV, the
only significant effect is the appearance of slightly wider
distributions.
III. MODULATION
In the following we present some typical examples of
steady states generated by the dynamics of the model
that clearly show the emergence of patterns for some spe-
cific values of the parameters.
For a global competition that results to be extremely
long-ranged (large C values, in relation to the values of
parameters σ and K), the steady state is characterized
by a spatially homogeneous occupancy. If the σ value is
sufficiently large or/and the K value sufficiently small,
totally homogeneous distributions are obtained (Figure
1), otherwise the solution is smooth but with the pop-
ulation concentrated in one region of the ring, with its
width controlled by σ.
As stated above, a simple heuristic analysis of eq. (1.1)
in the Fourier space shows that there is a necessary con-
dition for the emergence of inhomogeneity: the Fourier
transform of the influence function must have negative
values and large enough magnitude [12,13]. A Gaus-
sian in an infinite domain has a positive counterpart in
Fourier space and so does not match such requirements.
In contrast with these results, the fluctuations present in
our individual based model arrive to excite one specific
mode and modulations of this wavelength appear. The
tuning of the parameters allows modulations of arbitrary
wavelengths (Figure 1).
When C is decreased, the competition between modes
becomes stronger and no single mode dominates. In this
situation, some small regions of the ring are occupied
forcing all the remaining areas, up to some range, to
be nearly empty. The landscape results populated by
several living colonies divided by dead regions. There is
almost no competition between individuals of different
colonies and the space separating them can be identified
with an effective interaction length. This steady state
(spiky state [10]) corresponds to a sequence of isolated
colonies (spikes) and seen in the Fourier space, many
active wavelengths contribute to it (Figure 2).
Finally, for extremely short-ranged competition, in re-
lation to the σ value, no collective cooperation between
different excited modes emerges and a noisy spatially
homogeneous distribution appears (Figure 2).
We describe these paradigmatic steady states of the
system by characterizing the related spatial structure
with the help of a structure function S(q) [14] defined
as follows:
S(q) =
exp[i2πq · xj(τ)]
(3.1)
where the sum is performed over all individuals j with
their positions determined by xj(τ). Note, that, for con-
venience, the structure function is calculated only over
the closed interval [0, L]. The function is averaged over
some time interval T in order to avoid noisy data. S(0)
corresponds to the square of the mean number of individ-
uals in the system. The maxima of S(q) identify the rel-
evant periodicity present in the steady state. We will see
that the position of the global maximum (qM ) provides
an appropriate order parameter for the identification of
the segregation transition.
In our study we explored two different initial condi-
tions. In the first (local i.c.), the colony is located in
a finite and compact region of the space. In the second
(global i.c.), the individuals are spread all over the space.
The final distribution is independent of this choice and,
generally, local initial conditions make the system reach
the steady state more slowly. For this reason, if not dif-
ferently specified, our results are obtained from global
initial conditions.
IV. FINITE SIZE EFFECTS
Our analysis starts by exploring the model dependence
on the space size L. The reason for such interest is given
by the necessity of testing if the pattern formation is not
merely a product of some finite size effect. This is im-
portant, in the light of what was reported by Fuentes
et al. in ref. [12]. In their work, a numerical solution
of eq. (1.1) with a Gaussian influence function, showed a
segregation transition. But such a transition was just the
effect of the finite domain size that acted like a cut-off
for the Gaussian. Evidence of this interpretation came
from the observation that the amplitude of the pattern
depended on the ratio of the standard deviation of the in-
fluence function to the domain size - the critical values of
the standard deviation corresponding to the segregation
transition depended linearly on the domain size - and the
same patterns appeared for a modified Gaussian, which
vanishes abruptly beyond a cutoff.
The study of our individual based model gave differ-
ent results. By running some simulations with exactly
the same parameters but changing the ring extension, we
were able to show that the system is not influenced by the
domain size. If we choose data from spiky steady states,
that permit clear quantitative measures, it is possible to
remark that the general morphology of the patterns does
not change increasing the L value. In fact, both the pop-
ulation density and the mean number of peaks per space
interval remains constant. Moreover, in order to provide
a more precise test of possible small variations in the
distribution, we measured the cluster size. This quan-
tity was calculated by evaluating the standard deviation
(< x2 >i −< x >i2)1/2 of the position of the i individu-
als confined in each peak, then averaged over the different
peaks present at step τ and, finally, averaged over many
time steps after the system has reached the steady state.
Varying the system size caused no changes in the clus-
ters size (see Figure 3). From this result, we concluded
that the general aspect of the steady state does not
change with L. In particular, in contrast with what hap-
pens when the mean field equation is solved numerically,
the patterns do not depend on the ratio C
. For example,
for C = 0.2 and L = 50 we obtained a spiky steady state,
for C = 0.004 and L = 1 (same ratio) we obtained an
homogeneous steady state. Taking into account these re-
sults, from now on, all our simulations are implemented
on a ring of size 1.
We have just shown how the average of the population
size < N > scales with L in the steady state. Now, we
will give, through a simple heuristic argument, an esti-
mation of < N > as a function of the parameters K and
C, that will be useful also in the following of our analysis.
We can assume that, locally and in the steady state,
the number of deaths must be, on average, compensated
by the number of newborns, in order to comprise a sta-
ble population. For this reason, the death probability P
must equal 1/2. Assuming that the number of neighbors
that compete with a single individual are the ones living
up to a distance C and that, in these surroundings, the
average density N/L can be considered to be uniform, P
reduces to K ·2C ·N/L. Thus, N ∝ L/(CK). Looking at
Figure 4 we can see that this crude evaluation, that ne-
glects diffusion, inhomogeneity and reduces the influence
function to a box, describes well the general behavior of
the data obtained from our simulations.
V. SEGREGATION TRANSITION
In the previous paragraphs we supported the fact that
steady states, depending on the parameter values, can
assume inhomogeneous spatial distributions. Now, we
will try to describe the transition towards these states
(segregation transition). The structure function intro-
duced in eq. (3.1) provides a proper order parameter
to describe this transition. Different regions in the pa-
rameter space, coinciding with different steady states,
correspond to different positions of the global maximum
(obviously we are not taking into account S(q) at q = 0)
of the structure function. The transition from a homo-
geneous to an inhomogeneous distribution (see Figure 2)
matches the jump of the position of the global maxi-
mum (qM ) to a clear integer value, corresponding to the
number of clusters present in the space. For this rea-
son we can characterize the transition by looking at the
shape assumed by S(q), or looking at the value of qM .
If the space is homogeneously occupied, the structure
function does not present an integer maximum. On the
contrary, the maximum is located at qM ≃ 1.4. This
value corresponds to an uniform distribution of individ-
uals in the interval [0,1], approximated by the expression
exp(i2πq ·x)dx
. The segregation transition is char-
acterized by the passage of qM from 1.4 to an integer
value as soon as a modulation becomes dominant. In
Figure 5 we show qM as a function of C, varying K and
σ. First of all, from the analysis of these data, we can
observe that the number of clusters scales as C−1 (or,
equivalently, the periodicity of the inhomogeneous phase
has wave lengths proportional to C). Moreover, a criti-
cal value of C exists for which the transition takes place.
This Ccritic grows with 1/K and with σ. An analysis
of the available data suggests the possible dependence:
Ccritic ∝ σ2/3K−1/3 . Finally, for larger values of the
parameter C, in this range of the parameters σ and K,
the distributions are characterized by just one peak.
For any value of the competition strength, as can be
seen in Figure 6, there exists a critical value σcritic, de-
pendent on C and K, above which no spatial structures
emerge. Another measurement, that permits us to state
this relation in a different and clearer way, is presented
in the next paragraph.
VI. CLUSTER SIZE AND FLUCTUATIONS
In the following we describe, in the light of the existing
literature, the cluster size dependence on diffusion rate
and population size for the two different implementations
of the model.
We start by analyzing the typical size S of the clus-
ters that appear in the spiky phase. The data exposed in
Figure 7 show a dependence of the cluster size on the dif-
fusion coefficient: S ∝ σ, equivalent to S ∝
D. These
results are in accordance with the data presented in ref.
[33] obtained from an individual based model. In addi-
tion, this work pointed out how this behavior deviates
from the conclusions obtained from the deterministic ap-
proximation, where the cluster size was only weakly de-
pendent on the diffusion coefficient; another fact support-
ing the relevance of fluctuation effects in these systems.
We can easily interpret the dependence of the clus-
ter size on the diffusion coefficient by assuming that the
individuals confined in a cluster diffuse a distance pro-
portional to
DJ where J is the number of jumps the
individual performs in its life. In the case of the first
implementation of the model, where the diffusion is due
to the mutation process, J is obviously 1. Similar results
are obtained with the second implementation (see Fig-
ure 7), apart from a slightly wider cluster size (in this
case, individuals move, on average, more than just once
during their lifetime). Even so, the data show the same
dependence on the diffusion coefficient.
Finally, we present S as a function of K: S ∝
The reasons for this behavior are already explained in
ref. [34]: the cluster size is not controlled just by the
single individual’s number of jumps; in fact the diffusive
process continues with its descendants. For this reason,
it is proportional to the mean lifetime of a family, esti-
mated to be proportional toK−1 (see ref. [34] for details).
We introduce a new quantity that is useful for describ-
ing the existence of a critical diffusion value, giving an
estimation of its dependence on others parameters and
a confirmation of previous results. This quantity, which
we call mobility M(τ), estimates the mean mobility of
individuals. At a given time step τ , we choose an indi-
vidual i. Then we look for the closest agent, among all
the population, at time step τ − 1. We identify as di the
distance between these two individuals. Averaging over
the entire population N(τ) we obtain:
M(τ) =
di. (6.1)
The values assumed by M on varying the parameters σ
and C are shown in Figure 8. It is easy to distinguish two
clearly different behaviors. If the system is organized in a
spiky state (when σ ≤ σcritic), M(τ) ∝ σ. M is another
way of measuring the mean distance that an individual
moves during its lifetime inside the region defined by the
cluster. For this reason, this measure is coherent with the
data obtained from the direct evaluation of the standard
deviation of the clusters. In contrast, when the system is
organized in the homogeneous phase (when σ > σcritic)
M becomes independent of σ and is proportional to the
inverse of the occupation density M(τ) ∝ KC. The val-
ues of the mobility obtained from simulations with dif-
ferent values of C and K can be easily collapsed into
one function (see the inset of Figure 8). The collapse is
performed using the scaling σ → σ/C
K. This indicates
that the characteristic value of the crossover, σcritic, that
separates the two different behaviors of M(τ), scales as:
σcritic ∝ C
K. (6.2)
We conclude our study with some measurements try-
ing to catch some properties of the system fluctuations.
First, we estimated the fluctuations of the total popula-
tion, averaging over different simulations. The variance
turned out to be constant throughout the time evolution
and of the order of the square root of the total population.
The mechanism of auto-regulation of the population di-
mension does not allow the growth of big differences in
the total number of individuals.
For this reason, we focused our attention on the spa-
tial distribution of the population and tried to measure
some properties of these fluctuations. We studied the
variation of the local number of individuals in the same
simulation, for different times. We analyzed the evo-
lution of the system starting from local initial condi-
tions, with the population concentrated in the interval
[0.49, 0.51]. In this situation, the system evolves in time
with a small cluster fluctuating around the initial space
interval. This situation changes when a branching event
occurs that generates two well-defined clusters. Our in-
terest is in showing the behavior of local space fluctua-
tions and capturing possible variations in correspondence
with the branching event. First of all, we looked at the
mean value of the spatial local fluctuations Fs(τ), defined
as Fs(τ) = (
f2j )
1/2, where fj is the occupancy varia-
tion of the bin j from time step τ − 1 to time step τ . We
performed the average over all the b bins the ring was
divided in and obtained Fs(τ) =
N(τ), with no rel-
evant variations throughout the time evolution, even in
the time interval corresponding to the branching event.
More interesting is the shape of the frequency distribu-
tion of the size of fj . In fact, a simple Gaussian does
not fit this distribution, that presents extended tails (see
Figure 9). Throughout the system time evolution, the
shape of the normalized distribution is conserved. For
global initial conditions the same frequency distribution,
with extended tails, is recovered at the steady state. It
is identical to the one obtained with local initial condi-
tions and measured at the steady state. We think that
the deviation of the distribution from a Gaussian can be
considered as a hint that fluctuations play a relevant role
in the dynamics of the systems.
VII. CONCLUSIONS
We presented some results regarding clustering of or-
ganisms caused by a frequency-dependent interaction
that represents competition. We showed how this way
of modeling competition can be used not only to de-
scribe spatial phenomena in population biology, but also,
through a more abstract interpretation, to test ideas of
evolutionary theory (for example, studying the speciation
process).
From this unifying perspective, our study, obtained
from an extensive collection of data coming from sim-
ulations of an individual based model with global com-
petition, pointed out the relevance of fluctuation effects
in pattern formation. For the influence function adopted,
the mean field description predicts the absence of spatial
structures. On the contrary, fluctuations are able to ex-
cite the emergence of well defined patterns, which can
not be generated from a deterministic instability.
Furthermore, we discussed other fundamental proper-
ties of our model in the light of the existing literature,
unfolding a comparison with other models that describe
spatial segregation originated by some deterministic in-
stability. We showed that the observed patterns are not
due to a finite size effect, we characterized the behavior
of the segregation transition in various regions of the pa-
rameter space and we studied the existence of a critical
diffusion value. We analyzed the dependence of the clus-
ter size on the diffusion coefficient and pointed out some
characteristics of the fluctuations of the system.
ACKNOWLEDGMENTS
We are grateful to J.S. Sá Martins for a critical reading
of the manuscript and thank the Brazilian agency CNPq
for financial support.
[1] S. A. Levin and L. A. Segel, SIAM Rev. 27, 45 (1985); R.
Durrett and S. A. Levin, Philos. Trans. R. Soc. London,
Ser. B 343, 329 (1994).
[2] Y.-C. Zhang, M. Serva and M. Polikarpov, J. Stat. Phys.
58, 849 (1990).
[3] N. M. Shnerb, Y. Louzoun, E. Bettelheim and S.
Solomon, Proc. Natl. Acad. Sci. U.S.A. 97, 10322 (2000).
[4] B. Houchmandzadeh, Phys. Rev. E 66, 052902 (2002).
[5] W. R. Young, A. J. Roberts and G. Stuhne, Nature 412,
328 (2001).
[6] E. Ben-Jacob, I. Cohen, and H. Levine, Adv. Phys. 49,
395 (2000); J.Wakita, K. Komatsu, A. Nakahara, T. Mat-
suyama, and M. Matsushita, J. Phys. Soc. Jpn. 63, 1205
(1994).
[7] K. S. McCann, J. B. Rasmussen, and J. Umbanhowar,
Ecol. Lett. 8, 513 (2005).
[8] J. B. Wilson and A. D. Q. Agnew, Adv. Ecol. Res. 23,
263 (1992); R. Lefever and O. Lejeune, Bull. Math. Biol.
59, 263 (1997); N. M. Shnerb, P. Sarah, H. Lavee, and
S. Solomon, Phys. Rev. Lett. 90, 038101 (2003).
[9] N. F. Britton, SIAM J. Appl. Math. 50, 1663 (1990), S.
A. Gourley and N. F. Britton, J. Math. Biol. 34, 297
(1996).
[10] Y. E. Maruvka and N. M. Shnerb, Phys. Rev. E 73,
011903 (2006).
[11] A. Sasaki, J. Theor. Biol. 186, 415 (1997).
[12] M. A. Fuentes, M. N. Kuperman and V. M. Kenkre,
Phys. Rev. Lett. 91, 158104 (2003).
[13] M.A. Fuentes, M.N. Kuperman and V.M. Kenkre, J.
Phys. Chem. B 108, 10505 (2004).
[14] E. Hernandez-Garcia and C. Lopez, Phys. Rev. E 70,
016216 (2004).
[15] J. Billingham, Nonlinearity 17, 313 (2004).
[16] G. Flierl, D. Grunbaum, S. Levin and D. Olson, J. Theor.
Biol. 196, 397 (1999).
[17] B. M. Bolker, Theor. Popul. Biol. 64, 255 (2003).
[18] D. A. Kurtze, Phys. Rev. B 40, 11104 (1989).
[19] F.-J. Elmer, Phys. Rev. Lett. 70, 2028 (1993).
[20] J. Roughgarden, Am. Nat. 106, 683 (1972).
[21] F. Bagnoli and M. Bezzi, Phys. Rev. Lett. 79 3302
(1997).
[22] U. Dieckmann and M. Doebeli, Nature 400, 354 (1999);
M. Doebeli and U. Dieckmann, Nature 421, 259 (2003).
[23] E. Brigatti, J. S. Sá Martins and I. Roditi, Physica A
376, 378 (2007); V. Schwämmle and E. Brigatti, Euro-
phys. Lett. 75, 342 (2006).
[24] S. Pigolotti, C. Lopez and E. Hernandez-Garcia, Phys.
Rev. Lett. 98 258101 (2007).
[25] H. Sayama, M.A.M. de Aguiar, Y. Bar-Yam and M.
Baranger, Phys. Rev. E 65, 051919 (2002)
[26] H. Sayama, L. Kaufman and Y. Bar-Yam, Phys. Rev. E
62, 7065 (2000).
[27] D. R. Nelson and N.M. Shnerb, Phys. Rev. E 58, 1383
(1998).
[28] R. A. Fisher, Ann. Eugenics 7, 353 (1937); A. Kolmogo-
roff, I. Petrovsky and N. Piscounoff, Moscow Univ. Bull.
Math. 1, 1 (1937).
[29] F. J. Weissing and J. Huisman, J. Theor. Biol. 168, 323
(1994).
[30] J. Polechova and N. H. Barton, Evolution, 59, 1194
(2005).
[31] M. Meyer, S. Havlin and A. Bunde, Phys. Rev. E 54,
5567 (1996).
[32] C. Lopez and E. Hernandez-Garcia, Physica D 199, 223
(2004).
[33] E. Hernandez-Garcia and C. Lopez, Physica A 356, 95
(2005).
[34] E. Hernandez-Garcia and C. Lopez, J. Phys.: Condens.
Matter 17, S4263 (2005).
FIG. 1. Homogeneous steady state distribution ( top,
C = 4.0, 1/K = 80000, σ = 0.01) and modulated steady state
distribution (bottom, C = 0.9, 1/K = 18000, σ = 0.01). The
insets show the structure functions, S(q), of the correspond-
ing simulations. We show the distributions at time step 1000,
whereas the structure functions are averaged over 500 time
steps.
FIG. 2. Spiky (top, C = 0.059, 1/K = 500, σ = 0.001) and
homogeneous (bottom, C = 0.005, 1/K = 500, σ = 0.001)
steady states. The insets show the structure functions S(q).
We show the distributions at time step 2000, whereas the
structure functions are averaged over 1000 time steps. The
transition between these two states, in this typical range of
parameters, has been extensively studied.
0 10 20 30 40 50
0.005
FIG. 3. Dependence of the cluster size on the ring size
L. We present data from the model with mutation (trian-
gles) and from the one that implements diffusion (circles),
C = 0.2, K = 0.0029, σ = 0.001. The average is carried out
over all the clusters present at a given time step and over
different time steps.
FIG. 4. The number of individuals N present in the steady
state is proportional to (CK)−1. This result is in accordance
with the one obtained for a box-type influence function of
length C (see ref. [32]). We present data for different simu-
lations with 1/K ∈ [50, 500] and σ ∈ [0.0001, 0.01]. This last
parameter does not influence the final number of individuals.
The dashed line has slope -1.
FIG. 5. Segregation transition at Ccritical. Upper figure:
variation in dependence of K, where 1/K = 50, 150, 400, 500
and σ = 0.001. Lower figure: variation in dependence of σ,
where σ = 0.0005, 0.001, 0.002, 0.005, 0.01 and 1/K = 200.
FIG. 6. As shown in these figures, a critical σ value exists,
above which no spatial structures emerge. Upper figure: vari-
ation in dependence on C; we set 1/K = 100. Lower figure:
variation in dependence on K; we set C = 0.01.
FIG. 7. Top: Cluster size as a function of σ; 1/K = 300.
The solid line has slope 1. Bottom: cluster size as a function
of 1/K; σ = 0.0001. The solid line has slope 1/2. Trian-
gles represent data from the simulations where diffusion is
implemented through mutations, circles for the direct imple-
mentation of the diffusive process; we set C = 0.09.
FIG. 8. Mobility dependence on the diffusion parameter
σ for different C values; K = 0.01. In the inset, the data
collapse for arbitrary values of the parameters C and K.
FIG. 9. Spatial local fluctuation distribution at the steady
state: deviation from a Gaussian (C = 0.1, K = 0.00005,
σ = 0.0017). Data are averaged over 5 time steps. The con-
tinuous line is the best fit Gaussian.
