• The robustness of the Webworld model to changes in its structure

  • FileName: lm08a.pdf [read-online]
    • Abstract: 1 distributed entries and C ¼ 0:2; 0:5; 1:0. In both cases the ... For large C values (e.g. C > 0:25), the model succeeds in producing web structures, for values less than this bound, the model is unable to ...

Download the ebook

Alan J. McKane Theoretical Physics Group, School of Physics and Astronomy University of Manchester Carlos A. Lugo population dynamics Manchester M13 9PL, UK article info species ecological complexity 5 (2008) 106–120
available at www.sciencedirect.com
journal homepage: http://www.elsevier.com/locate/ecocom
The robustness of the Webworld model to changes in its
Carlos A. Lugo, Alan J. McKane *
Theoretical Physics Group, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK
article info abstract
Article history: We investigate the robustness of a co-evolutionary model of multi-species communities to
Received 10 January 2007 changes in its structure. The model has both a conventional population dynamics and a
Received in revised form network dynamics through which new species are introduced and unsuccessful species
21 June 2007 removed. The effect of choosing different forms of population dynamics has already been
Accepted 21 June 2007 investigated for this model. Here we study the effects of changing the nature of species
Published on line 3 April 2008 interactions and thresholds such as the number of individuals of a species required to be
present before that species is said to be extinct. We find that the model is remarkably robust
PACS: to these changes, with only a very few of the modifications leading to food webs which differ
87.23.Kg appreciably from those found in the original model. This indicates that many aspects of the
87.23.Cc structure of the model are not, in general, a factor in determining the nature of the food webs
89.75.Fb constructed.
# 2008 Elsevier B.V. All rights reserved.
Food webs
Evolutionary dynamics
Coarse-grained speciation
Extinction thresholds
1. Introduction Tokita and Yasutomi, 2003; Yoshida, 2003a). In these models,
food webs are recognised as networks that are the result of
In the last decade models of food-web structure have moved some dynamical process, rather than static entities, which
from simple, and mainly static (Polis and Winemiller, 1996), to was the assumption made in previous approaches (Williams
more realistic, and dynamic, representations of predator–prey and Martinez, 2000; Cattin et al., 2004).
networks (de Ruiter et al., 2005; Pascual and Dunne, 2006). The This increased realism has invariably meant that more
latter include models which describe food-web evolution choices have to be made in the construction of models, and
based on the aggregation of species by some evolutionary part of the modelling procedure should include an investiga-
mechanism, and operate at a time-scale we denote as tion into the robustness of the results produced by the model
evolutionary (e.g. Ito and Ikegami, 2003; Rossberg et al., when these underlying choices and assumptions are altered.
2006a,b and references therein), as well as models which seek In this paper we will carry out this procedure for one of the
to describe the co-evolution of ecological communities based most studied models of this type, the Webworld model of
on predator–prey relationships, by coupling the evolutionary food-web dynamics introduced by Caldarelli et al. (1998) and
dynamics with population dynamics (Bastolla et al., 2002; further improved by Drossel et al. (2001). This later version is
* Corresponding author. Tel.: +44 161 275 4192; fax: +44 161 306 4303.
E-mail addresses: [email protected] (C.A. Lugo), [email protected] (A.J. McKane).
1476-945X/$ – see front matter # 2008 Elsevier B.V. All rights reserved.
ecological complexity 5 (2008) 106–120 107
Box 1. The Webworld model dynamics is closely coupled to the population dynamics by
virtue of a relation between the fraction of time that a given
species predates on a given prey and the corresponding
functional response. As a consequence, this adaptive dynamics
is determined after a choice for the population dynamics has
been made. This relationship is an evolutionarily stable strategy
(ESS); a strategy which cannot be invaded by another, and so
this aspect of the model seems to us to be on a reasonably sound
footing. It is modelled by Eqs. (1) and (2) in Box 1. There gi j ðtÞ and
f i j ðtÞ represent respectively the consumption rate and the
amount of effort that a predator i puts in searching and hunting
its prey j. These quantities depend on each other, so, in practice,
both equations are iterated until a fixed point of Eqs. (2) is
reached, and therefore until the evolutionary stable strategy
condition is fulfilled (for further details see Drossel et al. (2001)).
For the population dynamics the model employs Eq. (3) in
Box 1 to describe the change of population sizes Ni with time.
The motivation for this choice is discussed by Drossel and
McKane (2003). The consequences of changing the form of this
equation and the functional response Eq. (1) has already been
investigated in some detail by Drossel et al. (2004). This study
showed that not all types of population dynamics – within an
evolutionary context – lead to non-trivial food-web structures.
In essence it was found that an artificial mechanism had to be
introduced which restricted predator choice to prey to which
they were best suited to exploit, or a particular type of
functional response which achieved this naturally had to be
utilised. The main question which we wish to address here is:
are there corresponding restrictions for the other types of
dynamics in the Webworld model?
This question leads us to investigate the more novel, and
less well-studied, evolutionary dynamics of the model. The
dynamic modelling of food webs has had to come to terms
with describing species in a more fine-grained way than is
common in textbooks on theoretical ecology, where species
may be referred to only as ‘‘predators’’ or ‘‘prey’’ (Maynard
Smith, 1974; Pielou, 1977) or in models of food-web assembly
where they may be referred to as ‘‘plants’’, ‘‘herbivores’’ or
‘‘carnivores’’ (Law, 1999). If the intention is to construct food
webs where new species are introduced into the network as
variants of existing species or as immigrants, we have to
characterise the species in some way, so that we specify these
variants. It is the nature of this characterisation, the way that
the one we refer to as the standard version of the model. the effectiveness of predators is determined and other choices
Webworld has essentially the two types of dynamics men- of the mapping from the definition of species into the
tioned above (Box 1 summarises the algorithm and the parameters of the predator–prey dynamics, that will be of
corresponding mathematical expressions). Along with these, interest to us here. We will vary the number of features that
the model includes an extra dynamical process at a shorter characterise species, the distribution of the scores of these
timescale than the population dynamics. This is adaptive features against each other, and generally investigate the
foraging, whereby predators constantly change their choice of robustness of the model webs to variation of the underlying
prey according to a strategy which adapts to the changing definition of species, their traits and the effectiveness of these
nature and population sizes of the prey. traits against each other. All the results we present were
Adaptive foraging has been studied in a similar way and at obtained by performing over 30 runs of the model and
this timescale in other models (Kondoh, 2003). It has proved to averaging the results. In each of the tables the mean and
be a valuable aid to the understanding of the stability of the standard deviations are shown.
ecosystems, despite their high degree of complexity. A static The outline of the paper is as follows. In Section 2 we
model employing foraging theory (Matsuda and Namba, 1991) modify the actual definition of species themselves, by
has also been put forward in order to understand some changing the number of features required to define a species,
topological properties of food webs, giving results which are as well as the overall number of species which are available. In
in good agreement with real data. In the Webworld model, this Section 3, we present results obtained by modifying the way
108 ecological complexity 5 (2008) 106–120
scores of features against each other are chosen, and examine
the case where only a fraction of the features are correlated
with each other. The final modification, described in Section 4,
involves changing the threshold for the number of new
species introduced after a speciation, and also the extinction
threshold. We end with a comparison with other models, a
discussion of the results obtained and their implications in
Section 5. A straightforward calculation, leading to an
approximate analytic expression for the distribution of
species’ scores in a particular case, is described in Appendix A.
2. Changes in the basic elements of the model
2.1. Modification of species’ definition
Fig. 1 – Averaged number of species for different values of
In this section we will begin by discussing the results of the K. A smooth variation in system size occurs as K
modification of the structure of the original Webworld model. increases. This constitutes the main trend observed. For
In this model species are made up of features: integers each of larger values of K the changes are less marked. This is also
which represents a distinct phenotypic or behavioural trait. true for the other web properties.
We start with a rather straightforward modification of the
model, consisting of exploring its sensitivity to different
values of the number of possible features, K, and the number
of distinct features that each species has, L. In early studies, other web properties remain essentially unchanged. Fig. 1
the values taken for these numbers were K ¼ 500 and L ¼ 10. shows time-series of the number of species for several values
Now, if the number of distinct species which are possible is of K. The role played by NRS seems not to be important; it lies in
denoted by NRS , then NRS ¼ K!=L!ðK À LÞ! For K ¼ 500 and L ¼ 10, the interval ð1013 ; 1020 Þ when K varies from 100 to 500 for L
NRS is of the order of 1020 , and so the ‘‘choice’’ of possible new fixed and equal to 10. So, as long as we keep the number of
species is to all intents and purposes infinite. Therefore one realisable species large enough, it seems that the model will
would expect that the choice of K and L would not be too reproduce webs with the same ecological and topological
important as long as NRS is sufficiently large, but clearly this properties.
needs to be investigated.
2.1.2. Variation of L
2.1.1. Variation of K Results obtained keeping the number of possible features K
The first results we discuss correspond to the case of varying fixed and allowing the length of a species string to increase its
the size, K, of the urn from which features are chosen, keeping value are shown in Table 2. In contrast to the case of K,
the species string length L fixed and equal to 10. Table 1 variations in L may have a larger role in the dynamics of the
summarises these results. Note that Webworld has four main system, since it appears explicitly when evaluating quantities
parameters: the rate of input of external resources, R, the such as the score matrix, Si j , and the competition matrix, ai j ,
saturation constant, b (see Eq. (1)), the average ecological which are discussed in more detail below. As we shall see, in
efficiency, l (see Eq. (3)) and the competition constant, c (see simple cases it is possible to obtain approximate analytic
Eq. (5)) (Drossel et al., 2001). expressions for the distribution of the Si j . This allows us to
We can see that the number of species present in the gain insight into the role played by L and K in determining the
system is slightly sensitive to this variation, however, all the strengths of the scores.
Table 1 – Web characteristics and standard deviations obtained varying the urn size of available features, keeping the
species string length fixed and equal to L ¼ 10 (the parameter values employed in this case were R ¼ 105 , l ¼ 0:1, c ¼ 0:6
and b ¼ 5 Â 10À3 )
500 400 300 200 100
20 19 18 16 13
Approximate value of NRS 10 10 10 10 10
Number of species 48:0 Æ 5:0 46:0 Æ 11:0 49:0 Æ 7:0 42:0 Æ 6:0 38:0 Æ 6:0
Fraction of basal species 0:12 Æ 0:02 0:13 Æ 0:03 0:12 Æ 0:03 0:14 Æ 0:02 0:14 Æ 0:02
Fraction of intermediate species 0:79 Æ 0:05 0:78 Æ 0:05 0:80 Æ 0:05 0:80 Æ 0:05 0:77 Æ 0:05
Fraction of top species 0:09 Æ 0:04 0:09 Æ 0:05 0:07 Æ 0:05 0:06 Æ 0:04 0:09 Æ 0:04
Max level 4:0 Æ 0:1 4:0 Æ 0:1 4:0 Æ 0:1 4:1 Æ 0:1 4:0 Æ 0:1
Links per species 1:6 Æ 0:2 1:6 Æ 0:2 1:7 Æ 0:2 1:7 Æ 0:2 1:7 Æ 0:2
Average level 2:4 Æ 0:1 2:3 Æ 0:1 2:3 Æ 0:1 2:4 Æ 0:1 2:4 Æ 0:1
ecological complexity 5 (2008) 106–120 109
Table 2 – Web features and standard deviations obtained by varying the length of the species feature array
5 10 20 50 70 100
11 20 35 69 86 107
Approximate value of NRS 10 10 10 10 10 10
Number of species 51:0 Æ 9:0 48:0 Æ 5:0 54:0 Æ 7:0 58:0 Æ 7:0 57:0 Æ 8:0 58:0 Æ 7:0
Fraction of basal species 0:13 Æ 0:03 0:12 Æ 0:02 0:12 Æ 0:04 0:09 Æ 0:04 0:09 Æ 0:04 0:09 Æ 0:03
Fraction of intermediate species 0:77 Æ 0:06 0:79 Æ 0:05 0:81 Æ 0:05 0:87 Æ 0:06 0:87 Æ 0:05 0:87 Æ 0:07
Fraction of top species 0:10 Æ 0:05 0:09 Æ 0:04 0:07 Æ 0:03 0:05 Æ 0:02 0:02 Æ 0:02 0:02 Æ 0:03
Max level 3:9 Æ 0:2 4:0 Æ 0:1 4:0 Æ 0:1 4:0 Æ 0:2 4:0 Æ 0:2 4:0 Æ 0:1
Links per species 1:8 Æ 0:2 1:6 Æ 0:2 1:6 Æ 0:2 1:6 Æ 0:2 1:6 Æ 0:2 1:6 Æ 0:2
Average level 2:4 Æ 0:1 2:4 Æ 0:1 2:4 Æ 0:1 2:4 Æ 0:1 2:4 Æ 0:2 2:4 Æ 0:2
We begin by recalling (Drossel et al., 2001) that the non-zero entries of the feature matrix, not to the connectance
competition scores Eq. (5) in Box 1 involve the fraction of of the food web itself.
common features qi j ¼ ni j =L, where ni j is the number of To analyse the role of L in Eq. (4), we consider the simple
common features between species i and j, and 0 c 1, the case of a feature matrix given by:
competition parameter. To define the species scores, Si j , a real
number mab is assigned to each pair of features a and b. This 8
defines a K Â K matrix whose entries are constrained by maa ¼ > À1 with probability C;
< 2
0 and mab ¼ Àmba . This antisymmetric matrix, called the mgd ¼ C g >d (7)
feature score matrix, or feature matrix, has off-diagonal >1
> with probability ;
> 2
elements which follow some probability distribution. Then 0 with probability 1 À C;
two different species, say i and j, are given a score according to
˜ ˜
Si j ¼ max ð0; Si j Þ, where Si j is given by Eq. (4) in Box 1. The along with the antisymmetric condition mgd ¼ Àmdg . To obtain
prefactor, LÀ1 , is included to keep the score of the order of ˜
an idea of how the distribution of the Si j depends on L (and K),
unity. Using this convention, we specify that a predator–prey we relabel the indices in the sum, obtaining a single variable
link between i and j may exist if Si j > 0, with i being the which runs from 1 to L2 . We now define the random variable
predator and j being the prey. Along with this, a subsidiary Y L2 ¼ L xa , where xa is given by mab =L. For simplicity, we
condition is required to establish a link between two species. will make the assumption that the xa ’s are independent vari-
This requirement is imposed on the fractional effort that every ables. This is clearly only approximately true: the species may
species assigns to each of its prey and is done to avoid the have features in common and we are ignoring the evolution-
proliferation of weak links across the network. So the ary nature of the original process. We show in Appendix A
condition Si j > 0, is a necessary, but not a sufficient, condition that, if we make this assumption, it is straightforward to
to create a link between species. If this is fulfilled, the criteria obtain a probability distribution for Y L2 , simply by evaluating
for link establishment is chosen to be f i j > 0:01. the expression:
In the standard version of the model, the feature score
matrix has Gaussian distributed random entries. This means X L2 
1X l l À 2j

that every two different phenotypical features are related by a PðyÞ ¼ pL
plÆ1 d yÀ (8)
l 2l j L
l¼0 j¼0
number taken from such a statistical distribution. This
assumption can be easily criticised. In the first place, it might
be that any two arbitrary features are not comparable. They with p0 ¼ ð1=KÞ þ ð1 À ð1=KÞÞð1 À CÞ and pÆ1 ¼ Cð1 À ð1=KÞÞ,
could describe a completely different aspect of a species respectively. Although Eq. (8) applies to the case of a very
phenotype. For example, it does not seem to be a correct way simple mab (Eq. (7)), and it depends on the assumption of
to compare some kind of aspect of the vision of an individual independence of the xa ’s, it predicts qualitatively the beha-
to some other, unrelated, anatomical trait. Secondly, if they viour observed during the actual process. Fig. 2 shows com-
are related, why would they follow this specific distribution? parisons between the observed distributions and the ones
Dealing with the second point first, we can modify the model obtained with Eq. (8) for several connectances. We can see
by testing distributions which are different from Gaussian. To that in this approximation the role played by L is more com-
deal with the first point, we can define a variant of the model, plex than the one played by K. The latter mainly narrows or
by introducing a new parameter: the connectance, C, of the widens the distribution through the result s 2 ¼ Cð1 À ð1=KÞÞ for
feature matrix, which specifies the fraction of entries that are the variance, meanwhile increasing L gives more possible
different from zero: ways to obtain a given result. At least for this simple choice
of the mab matrix entries distribution, K does not have any

x with probability C; influence whatsoever on determining the mean of the distri-
mab ¼ (6)
0 with probability ð1 À CÞ; bution of the scores.
For the competition matrix ai j , we can measure to some
where x is random variable following some given probability extent the way the competition evolves during an evolu-
distribution, not necessarily Gaussian. It is important to note tionary time-step, and the role played by L. Suppose we have a
that this definition of connectance refers to the fraction of parent species p, with competition matrix a p j , when a new
110 ecological complexity 5 (2008) 106–120
Fig. 2 – Left column shows histograms obtained from sampling the random variable Si j , for different values of the
connectance C. The right column shows histograms obtained by plotting the approximate analytic expression
given by Eq. (8).
child species c is added to the system. We can only have three given by
possibilities for the number of common features, with the
existing species j:
> 1 þ ð1 À kÞ;
< L
8 ð1 þ kÞ
< n p j þ 1; qch ¼ 1À ; (10)
> L
nc j ¼ n p j À 1; (9) >
> k
:n : >
:1 À :
If we assume that for a particular prey h, the parent From Eq. (10) it is obvious that, the larger the value of L, the
species differs in k features, then n ph ¼ L À k, and qch will be smaller the change in the fraction of common features. The
ecological complexity 5 (2008) 106–120 111
Fig. 4 – Average number of species for several values of L.
The number of species in the mature state and the number
of trophic species increase as L increases. This leads to
communities with a diverse configuration and many
equivalent species for L > 20.
2.1.3. Results for the variation of L
Table 2 summarises the results obtained using the same set of
parameter values employed studying variations in K (see
Table 1), which was fixed at a value K ¼ 500. We should stress
that these properties, as well as the properties presented in
Table 1 and the remaining tables, refer to trophic webs, which
Fig. 3 – (a) Typical time evolution of the number of species are webs composed of trophic species (species having the
for L ¼ 5 and 20. For these values the system reaches same set of predators and prey are identified). As can be seen
maturity within an interval of 105 evolutionary time-steps. from Table 2, there are variations in the fractions of top,
(b) Same as (a), but for L ¼ 70: Here the system evolves intermediate and basal species as L changes. Along with these,
more slowly. For large L, it may take a considerable time to as shown in Fig. 3 (a) and (b), not only does the number of time-
produce a set of species able to populate the system in a steps required to reach a mature state depend on L, but also
stable and diverse way. The transient states may consist the ‘‘total’’ number of species seems to increase as L does.
of systems composed of very few species. Fig. 4 shows that the difference in the number of species and
trophic species becomes more noticeable as L increases. This
is just a manifestation of what was stated before: for large
values of L, the system evolution is slower; mature states with
L large, consist of communities composed by species which
smaller the latter, the smaller the change in competition are equivalent, in terms of trophic species. In the cases L ¼ 5
scores and the scores Si j . So for a new species entering the and 10 this difference is in general no more than a few species
system, the effect of changes in the value of L is to slow down and the food webs and trophic food webs possess the same
the creation of a diverse community, because the larger L is, statistical properties (and in most cases are the same). All the
the smaller the change in the quantities that regulate the other values of L explored showed this proliferation of
interaction strengths among species. In other words, at every equivalent species, which may be the origin of the differences
evolutionary time-step, the scores and overlaps of a new in the fractions of top, intermediate and basal species shown
species entering the system are calculated, and they will not in the table. Fig. 5 illustrates this. There model food webs, both
be substantially different to those of the parent species if the total and trophic for the same system, are shown for L ¼ 20. As
value of L is large enough. This means that in the early stages can be observed, the number of species and trophic species
of the evolution, the system will be composed of very similar differ by about 20 species.
species, and then it may be very hard to grow a diverse Fig. 6 (a)–(d) shows results obtained for K ¼ 500 and L ¼ 100
ecosystem. However, the continuous iterative process of for the competition matrix distribution and the number of
speciation will, eventually, create species which are distant species present in the system for two different sets of random
enough in both overlaps and scores from those already numbers. The first case (Fig. 6(a) and (b)) illustrates the
existing. Then diversification and system growth will be situation where the system began to grow appreciably at a
observed. certain time during the evolution and eventually reached a
112 ecological complexity 5 (2008) 106–120
Fig. 5 – Model food web obtained for L ¼ 20. (a) Food-web structure showing all the species present in the system and (b) food
web obtained by including only trophic species. It can be seen that the number of nodes in the network is considerably
reduced, since for large L species are produced whose descendants are similar in terms of scores and fraction of common
stationary state (Fig. 6(b)). We can observe that the competi- variations which in all cases have an error interval with a non-
tion matrix distribution (Fig. 6(a)) has several broad peaks, empty intersection.
which should belong to very similar species, but it is more or
less well distributed over the whole interval ½c; 1Š. In contrast to
this, Fig. 6(c) and (d), show the case where, for the same 3. Connectance and distribution of the feature
parameters, the system did not grow in the same time interval. matrix
The competition matrix distribution shows a single broad
peak around a certain value, whilst the other regions of the In this section we examine the results obtained by employing
interval are empty, apart from near 1. This suggests that early different distributions for the mab matrix. In the original
in the history of the system a well adapted species appeared, Webworld model the entries for a > b (the matrix is antisym-
and that the only subsequently successful species were metric) were taken to be Gaussian random variables with zero
variants of this one which took over the community in a kind mean and unit variance. How does the structure of the webs
of red queen situation, since all the other ai j ’s lie near 1. In any change if we choose other distributions?
case, even for the situation where the system reaches a diverse For this purpose we kept all the model parameters fixed
mature configuration in reasonably short times (3 Â 105 for and equal to R ¼ 105 , l ¼ 0:1, c ¼ 0:7 and b ¼ 5 Â 10À3 . This
L ¼ 70 and 5 Â 105 for L ¼ 100), the measured properties show choice was made because it produces webs which are
ecological complexity 5 (2008) 106–120 113
Fig. 6 – (a) Histogram of the competition scores for a typical realisation with L ¼ 100 and K ¼ 500, where the system was able
to form a web after a long-time interval. (b) Time-series corresponding the histogram in (a). (c and d) Same as (a) and (b), but
for a case when the system was unable to grow within the same time period (105 time steps).
moderately large in the standard Gaussian case, so we can We see that all the average properties of the mature
compare the results to those obtained previously. Table 3 ecosystems do not vary greatly from one case to another. This
shows the different web properties obtained for Laplace, is reasonable, and expected, because in each case we are
Cauchy, Gaussian and uniformly distributed entries for the employing smooth, symmetric and well behaved continuously
mab matrix. In all the cases the mean of these distributions was distributed random entries, which are being added to form the
zero and the scale parameter of the Laplace and Cauchy scores Si j , and we would not expect any dynamical aspect to be
distributions were unity. very sensitive to this. This led us to ask if we could choose the
Table 3 – Results obtained employing three different continuous random entries with the same set of model parameter
values: R ¼ 105 , l ¼ 0:1, c ¼ 0:7 and b ¼ 5 Â 10À3 (the standard deviations corresponding to each quantity are also shown)
Feature score distribution
Laplace Cauchy Gaussian Flat (À1,1)
Number of species 42:0 Æ 7:0 38:0 Æ 9:0

Use: 0.3233