The characteristics of species in an evolutionary food web model

ARTICLE IN PRESS

Journal of Theoretical Biology 252 (2008) 649–661 www.elsevier.com/locate/yjtbi

The characteristics of species in an evolutionary food web model Carlos A. Lugo, Alan J. McKane Theoretical Physics Group, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK Received 20 July 2007; received in revised form 20 February 2008; accepted 20 February 2008 Available online 4 March 2008

Abstract We explore the consequences of modifying the way in which species are defined in an evolutionary food web model. In the original version of the model, the species were defined in terms of a fixed number of features, chosen from a large number of possibilities. These features represented phenotypic and behavioural characteristics of the species. Speciation consisted in occasionally replacing one of the features by another. Here we modify this scheme by firstly allowing for a richer structure and secondly by testing whether we are able to eliminate the need for an explicit choice of features altogether. In the first case we allow for changing the number of features which define a species, as well as their nature, and find that in the resulting webs the higher trophic levels typically contain species with the greatest number of features. In the second case, by a simplification of the mechanisms for inter and intra-species competition, we construct a model without any explicit features and find that we are still able to grow model food webs. We assess the quality of the food webs produced and discuss the consequences of our findings for the future modelling of food webs. r 2008 Elsevier Ltd. All rights reserved. Keywords: Co-evolutionary model; Multi-species communities

1. Introduction The subject of modelling ecological networks such as food webs has been a very active field of study in the last few years (de Ruiter et al., 2005; Pascual and Dunne, 2006), with a number of evolutionary models being proposed to describe the formation of community food webs (Drossel et al., 2001; La¨ssig et al., 2001; Yoshida, 2003b; Tokita and Yasutomi, 2003). In these models, new species are added to the food web on an evolutionary timescale, whereas extinction events, along with the increase or decrease of the population sizes, are usually governed by differential equations describing the population dynamics which operate on an underlying ecological timescale. These approaches may include several type of interactions such as mutualism, inter-specific and intra-specific competition, parasitism, etc. Several studies also include additional mechanisms such as foraging, which has been recognised to be relevant in adding stability to communities composed of 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). 0022-5193/$ - see front matter r 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2008.02.028

a large number of elements (Kondoh, 2003). In all these cases, the resultant formalism relies on the way model species are characterised. One approach favoured by several authors consists of defining model species with a number of internal degrees of freedom, which seek to describe traits at the phenotypic level or to describe some behavioural feature (Caldarelli et al., 1998; Yoshida, 2003a; Rossberg et al., 2006), and which determine the interactions between species in an explicit way. For instance, as explained later, in the Webworld model this is achieved by assigning scores to the available traits and using these to determine the interaction between species. Another common approach consists of taking a more coarse-grained characterisation of species by employing interaction matrices, where every entry represents the type and strength of the interaction between two given species (La¨ssig et al., 2001; Tokita and Yasutomi, 2003). Amongst all these approaches, the Webworld model introduced by Caldarelli et al. (1998), and subsequently refined by Drossel et al. (2001), has proved itself able to reproduce several patterns observed in real food webs (Quince et al., 2005). It models species using the first of the approaches mentioned above, and also takes foraging into

ARTICLE IN PRESS 650

C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

account. The original species definition used in the model determines in an explicit way all the quantities which play a role in the model population dynamics, such as the formation of predator–prey links, consumption rates and competition strengths, but it only plays an implicit role in the evolution rules which regulate the diet choice. It has been shown that the structure of the food webs produced by the model strongly depends on the type of population dynamics employed, for instance whether the functional response is of the Lotka–Volterra, Holling II or ratiodependent type (Drossel et al., 2004), and also on the use of an evolutionarily stable strategy (ESS) for the modelling of foraging behaviour. When the model employs those two types of dynamics it is able to produce food webs which can be compared with real webs over a wide range of parameter values (Drossel et al., 2001; Quince et al., 2005; Lugo and McKane, 2008). In this paper we study the role played by species definition in the Webworld model predictions by both making it more elaborate and also by simplifying it. In the former case, we employ model species which at every evolutionary time step, when new species are added to the system, can be created by two different speciation mechanisms. One of these is speciation by ‘‘mutation’’, where the new species is created by replacing one of the traits of the parent species—the same procedure as was adopted in the original model (Caldarelli et al., 1998; Drossel et al., 2001). The other is a new mechanism, consisting of a ‘‘complexification’’, whereby the new species is created by adding a new feature to the ‘‘phenotype’’ of its parent, instead of replacing one. The path towards a simplified species definition is investigated in the second part of the paper, and results in a model which is more closely related to those where species are defined by interaction matrix couplings. This new model variant requires a reformulation of all the aspects where the definition of species plays an important role, while at the same time keeping unchanged the dynamical aspects which do not explicitly depend on this definition. In the next section we give a detailed discussion of those elements of the model that change with the introduction of the new species definition or speciation mechanism. The introduction of the new speciation mechanism, along with the results obtained from it, is presented in Section 3, and the simplification of the model is analysed in Section 4. The paper concludes with a brief discussion of the main results obtained and the outlook for future work. For completeness, an appendix describes those aspects of the model which are not modified by the present work. 2. Model In addition to the discussion in the original papers (Caldarelli et al., 1998; Drossel et al., 2001), the model has been described in detail elsewhere (see McKane, 2004; McKane and Drossel, 2006, for example). Here we describe only those aspects which will differ from the standard

version of the model as a consequence of changes in the model species definition. As mentioned in the Introduction, the building blocks of every community food web model involve the definition of the species. In the Webworld model species are characterised by a set of L different traits, which we call features, which are represented by integer numbers taken from a set of K available features ðKbLÞ. This definition seeks to describe species at the phenotypic or behavioural level. Along with these sets of features, every species i in the system is also characterised by a population size N i ðtÞ. Once species are defined, a key aspect of the model is to establish criteria with which to compare species in order to define the existence of predator–prey relationships, competition between species, and so on. To do this, the model compares species in two ways: it assigns a species score for each pair of species and it assigns a competition score based on how similar two species are. The species score is calculated by attributing a score to every feature against all the other features available, and then comparing two species feature by feature. To formalise this, let a and b represent two different features. Then a score mab ð¼ mba Þ is assigned by taking a random variable distributed according to some probability density function, usually taken to be Gaussian with zero mean and unit variance. This way of giving a score to features defines a set of K  K numbers which can be arranged in the form of a matrix, the feature score matrix (FSM). Then, for every two species i and j, the sum over the scores of the features characterising them: 1 XX S~ ij ¼ mab , L a2i b2j

(1)

constitutes a measure of how well adapted one species is at preying on the other. The factor L1 in front of Eq. (1) is present in order to keep the scores of the order of unity. Due the antisymmetric condition on the FSM, S~ ij 40 always implies S~ ji o0. In other words, both cases encode the same information. For this reason, the species score is defined as Sij ¼ maxf0; S~ ij g.

(2)

This definition allows for the establishment of a criterion for a possible predator–prey relationship between species i and j: i is able to obtain resources from j if Sij 40. The other way in which species are compared is in the determination of the strength of the competition between them when obtaining resources. The type of competition assumed by the model is intended to capture the fact that the more similar two species are, the stronger the competition is to gain resources from a common prey. This, of course, implies that when species are identical, the competition score takes the maximum value. To express this in the form of an equation, two species, i and j, are compared feature by feature by counting the fraction of features common to both, qij , and then evaluating the

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

competition score: aij ¼ c þ ð1  cÞqij .

(3)

Here c is a model parameter lying between 0 and 1. So aij varies from c if i and j have no features in common to 1 if i ¼ j. So far we have mentioned only those aspects of the model which are needed to discuss the modifications which we are proposing. Other aspects, such as the dynamics and foraging strategy are described in the appendix and in the references given earlier. However, it remains to discuss the speciation mechanism used in the original version of the model. This is carried out through a mutation-like event. Specifically, one of the existing species is chosen at random and one of its features, also chosen at random, is replaced by a new one taken from the set of K possible features. The population of the parent species is then decreased by 1 and the population of the new species is taken to be 1. In the next section we modify this speciation mechanism and examine the effect that this has on the web structure. In addition to the mutation-like mechanism described above, we now introduce another mechanism, which we call complexification. This is carried out by choosing one of the pre-existing species at random, and adding a new feature, without removing any of the existing ones, thus increasing the length of the feature string by 1. If the added feature is already present, the process is stopped and another attempt made. Using the above rules, as well as those given in the appendix, the model may be simulated and the average properties of the food webs generated determined. These are found for given values of the four main parameters of the model, namely the external resources R, the saturation constant b, the competition parameter c and the ecological efficiency l. In the forthcoming sections, the results of these simulations are discussed. The results are summarised in Tables 1–7 and represent the average of the outputs of the various realisations. The number of realisations (runs) performed in each case was in the range 40–70. The standard deviations obtained in each case did not tend to exhibit significant fluctuations around the mean values, except in quantities such as the fractions of top, Table 1 Results obtained for simulations employing the two speciation mechanisms described in the text, with a probability of 0.9, 0.95 and 0.99 for a mutation to occur (0.1, 0.05 and 0.01, respectively, for a complexification event) m

0.90

0.95

0.99

Number of species Fraction of basal species Fraction of intermediate species Fraction of top species Maximum level Links per species Average level Mean species length

63.0 0.12 0.85 0.04 4.0 1.6 2.4 31.0

58.0 0.10 0.83 0.07 4.0 1.6 2.4 22.0

54.0 0.12 0.79 0.08 4.0 1.6 2.3 12.0

The minimum admissible value for a species string was set to be Lmin ¼ 10.

651

Table 2 Food web properties measured for values of the competition parameter close to unity Competition parameter

0.85

0.9

0.95

0.99

Number of species Fraction of basal species Fraction of intermediate species Fraction of top species Maximum level Links per species Average level

134.0 0.15 0.81 0.04 4.0 2.8 2.1

127.0 0.16 0.81 0.04 4.0 1.9 2.2

98.0 0.19 0.77 0.05 4.0 1.9 2.2

42.0 0.22 0.70 0.09 4.0 1.4 2.2

The number of species present in the system is quite sensitive to changes in the value of c. The parameter values employed were R ¼ 105 , l ¼ 0:1 and b ¼ 5  103 .

Table 3 Web characteristics produced by the model without features for several values of c, and parameter values of R ¼ 105 , b ¼ 5  103 and l ¼ 0:1 Competition parameter

0.2

0.3

0.4

0.45

0.5

0.6

Number of species 261.0 124.0 66.0 44.0 36.0 20.0 Fraction of basal species 0.17 0.24 0.50 0.53 0.53 0.65 Fraction of intermediate species 0.67 0.57 0.31 0.25 0.24 0.11 Fraction of top species 0.16 0.19 0.19 0.22 0.23 0.24 Maximum level 3.0 3.0 3.0 3.0 3.0 2.5 Links per species 4.5 3.4 2.4 2.0 1.9 1.7 Average level 1.9 1.9 1.5 1.6 1.5 1.4

Table 4 Web characteristics produced by the model without features for several values of c, and parameter values of R ¼ 105 , b ¼ 5  103 and l ¼ 0:4 Competition parameter

0.5

0.6

0.7

Number of species Fraction of basal species Fraction of intermediate species Fraction of top species Maximum level Links per species Average level

191.0 0.11 0.79 0.09 4.0 3.2 2.5

60.0 0.31 0.58 0.11 4.1 2.1 2.1

39.0 0.37 0.52 0.11 3.8 1.9 2.0

Table 5 Food web properties obtained with the standard version of the Webworld model for several values of the competition parameter c in Eq. (3), and parameter values R ¼ 105 , b ¼ 5  103 and l ¼ 0:1 Webworld model

Number of species Fraction of basal species Fraction of intermediate species Fraction top species Max level Links per species

c ¼ 0:2

c ¼ 0:4

c ¼ 0:6

c ¼ 0:8

196 0.02 0.90 0.08 3.03 2.4

79 0.08 0.90 0.02 3.69 2.3

55 0.09 0.90 0.01 3.91 1.7

27 0.12 0.86 0.02 4.0 1.7

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

652

Table 6 Number of species, number of links and the maximum level observed by employing the Webworld model with dichotomous competition and the standard version of the model Webworld model

Number of species

Links per species

MaxLevel

Standard* D. comp 1y D. comp 2a D. comp 2b

196 134 124 191

2.4 2.8 3.4 3.1

3.7 4.0 3.0 4.0

The parameter values employed are R ¼ 105 and b ¼ 5  103 and, for the various cases: standard version (denoted by ) l ¼ 0:1 and c ¼ 0:4; dichotomous competition with features (denoted by y): l ¼ 0:1 and c ¼ 0:85; dichotomous competition without features (denoted by a): l ¼ 0:1 and c ¼ 0:4; dichotomous competition without features (denoted by b): l ¼ 0:4 and c ¼ 0:5.

Table 7 The number of species, links per species and distribution of species between trophic levels for 14 food webs (after Dunne et al., 2002), and for the model foodwebs shown in Figs. 3(b), 6(a,b), 8 and 9 Ecosystem name Number Links per Trophic level of species species 1 2

3

4

Number of species Bridge Brook Lake Scotch Broom Canton Creek Chesapeake Bay Coachella Valley El Verde Rainforest Grassland Little Rock Lake Skipwith Pond St. Marks Seagrass St. Martin Island Stony Stream Ythan Estuary 1 Ythan Estuary 2 Webworld 1 Webworld 2 Webworld 3 Webworld 4 Webworld 5

75

7.37

39

34

2

–

154 108 33 30 156

2.40 6.56 2.18 9.67 9.68

1 56 5 3 28

24 52 15 22 98

117 – 13 5 28

12 – – – 2

75 181 35 48

1.51 13.12 10.86 4.60

8 63 1 6

15 80 18 31

52 38 16 11

– – – –

44 112 134 92 69 110 29 75 32

4.95 7.43 4.46 4.58 1.79 1.8 1.3 2.2 1.78

6 63 5 5 8 18 7 36 14

29 46 44 44 32 54 12 36 10

6 3 81 42 27 36 9 3 6

3 – 4 1 2 2 1 – 2

was used in the original version of the model and the complexification mechanism. The parameters employed in the simulations were R ¼ 105 , l ¼ 0:1, c ¼ 0:6 and b ¼ 5  103 , which in the standard version of the model produce communities of a moderate size (see Quince et al., 2005, for instance). A previous study (Lugo and McKane, 2008), showed that if the length of the feature string, L, is large enough, then the evolution of the system slows down. In some cases the system remains in a state composed of only a few species. In order to prevent such a scenario, and to produce networks composed of species with different number of features, we ran simulations where speciation through mutation occurred with probability m, and speciation through complexification occurred with probability 1  m, where m is close to 1. We may characterise the system through the mean length of the feature string for species in the system. In the simulations we began with an initial string length of Lmin ¼ 10. The values chosen to perform a complexification were one, five and 10 complexifications per 100 speciation events i.e. m ¼ 0:99; 0:95 and 0.9. In order to incorporate this new speciation mechanism, it is necessary to redefine the Sij ’s and the aij ’s, because they depend explicitly on the length of the feature strings. In the former case, from Eq. (1) we see that the sum XX mab a2i b2j

intermediate and basal species. Here typical values of the standard deviation obtained were of the order of 50% in all three cases. In general, for the number of species, links per species, maximum and average level, the standard error observed were of the order of 10%, 25%, 0:04% and 0:08%, respectively. 3. Speciation–complexification results In this section, we will present results of the numerical investigations for both the mutation-like speciation which

should be normalised by the factor L1 to keep all the scores of the order of unity. If the feature strings have a different number of elements, a new normalisation factor, Qij , needs to be chosen. To do this, we define Qij ¼ ðmaxðLi ; Lj ÞÞ1 . The same choice is made for the competition matrix aij which appears in Eq. (3). This is only one of many possible definitions for Qij ; another one pffiffiffiffiffiffiffiffiffi could be Qij ¼ Li Lj , for example. The only requirement for this quantity is that it should reduce to 1=L when the species have the same number of features and that the normalisation factor should be reasonable in some sense. Since in the current model specific features are not identified with a particular biological trait, the number of possible features is a modelling choice, as is Qij . Three main types of behaviours were found, which are illustrated in Figs. 1(a) and (b), for the case of m ¼ 0:9. 1. Mature systems showing diversity and a large number of species: Large stable webs could be grown if the frequency of speciation events was such that mutationlike events were far more common than complexification events. More specifically, for m close to 1 most of the simulations were able to produce large stable webs, whereas for smaller values of m a large number of simulations had to be performed in order to obtain a reasonable number of webs with a significant number of species. The ecological properties obtained for these outputs, shown in Table 1, display a smooth downward trend in the number of species as m increases in the range

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

350

120 Mean Species String Length

A B C

100 Number of Species

653

80 60 40 20 0

A B C

300 250 200 150 100 50 0

0

2

4

6

10

0

x 104

2

4

6 Time

8

10 x 104

Relative Frequency

Relative Frequency

Time

8

0

0

0.6

0.8 Competition Score

1

0.8 Competition Score

1

Fig. 1. (a) The three regimes found for the two competing evolution processes described in the text, with m ¼ 0:9. (b) Evolution of the mean lengths of the feature strings for m ¼ 0:9. A corresponds to a process driven by complexification of a very well-settled ancestor (item 3 in Section 3). B corresponds to systems where a well-adapted species settles into the system leading to a diverse final state of species with a significant memory of its common ancestor (item 2 in Section 3). C corresponds to a system that evolves without any preference for a particular kind of species, so that the stationary state consists of a diverse ecosystem (item 1 in Section 3). (c) Histogram of the competition scores in a case where the system evolved into a diverse stationary state. (d) Histogram of the competition scores in a case where the system was unable to grow.

½0:90; 0:99, however, all the other properties remain fairly robust. The percentage of systems that give rise to large diverse webs increases with the value of m, and for the set of values explored here (0.9, 0.95 and 0.99) only the last two produced rich systems in the majority of the runs (90% and 94%, respectively, for 40 runs). On the other hand, for m ¼ 0:9 only 60% of 70 runs grew this type of web. Fig. 2 shows the averaged mean length for m ¼ 0:95 and 0:99. We can see that the averaged mean length in this regime grows very slowly (each of them is fitted to logarithmic functions of the form y ¼ A þ B ln x), contrary to what is observed in the cases discussed below. Typical behaviour of the number of species and mean length is illustrated in Figs. 1(a) and (b) (label C). One can ask if long string lengths would have behavioural consequences or if they would affect the survival rates of the species. These are complex questions, but in a previous study (Lugo and McKane, 2008) we have shown that species with large L have an advantage, which leads to better survival rates. Fig. 1(c) is a histogram obtained in the mature configuration for the competition scores, which shows that the distribution of scores take values in all of the bins in the interval

½c; 1. This indicates that a diverse community exists in this mature state. The other types of behaviour observed are discussed in the next two paragraphs. 2. Systems where a predominant species does not allow new species to settle into the network: This case is illustrated in Fig. 1(a) (label B), where a typical time-series of the number of species is shown. There we can see that the system only contains a small number of species, which are primary producers. Contrary to the previous case, this behaviour is found less frequently as m approaches 1 (40%, 9% and 8% for m ¼ 0:9; 0:95 and 0.99, respectively). It occurs when, at an early stage in the evolution, a species evolves into one that is better adapted to the system, displacing many of the other species in the system. This drives the system into a regime where only basal species survive, but as soon as new species are introduced, they take over the system or are eliminated. From Fig. 1(b) (label B), it is possible to conclude that this continuous overthrow of species is carried out by continuous ‘‘complexification’’ of a ‘‘superspecies’’, rather than by mutation-like events. In the figure, the time-series of the average species length is shown: it grows almost linearly at the beginning of the

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

654

stage in the evolution of the system, a given species with a feature string of a certain length finds itself very well adapted to the system, and all the mutants descended from it populate the system, competing between themselves for resources. Figs. 1(a) and (b) (label A) show typical outputs for the number of species and the mean feature-string length as functions of time in this case. This behaviour occurs with a very low frequency (1% and 2% for m ¼ 0:95 and 0.99, respectively), and was mainly found in the range of m40:9, but more like an exception than a trend.

Averaged Mean Length µ=0.95 30

Features String Mean Length

6.0+1.4ln(t)

25

20

15

10 0

2

4

6 Time Averaged Mean Length µ=0.99

8

10 x 104

14

Features String Mean Length

9.35+0.24ln(t)

13

12

11

10

9 0

2

4

6 Time

8

10 x 104

Fig. 2. The evolution of the mean length of the feature string for systems which create stable, diverse stationary states (a) for m ¼ 0:95 and (b) for m ¼ 0:99. The mean length will keep increasing due to the way speciation events are defined, however, for this case this growth is very slow. Here we compare both cases with functions of the form y ¼ A þ B ln x.

evolutionary process. This linear growth is eventually attenuated because, at some stage, the creation of a species by increasing the parent species features set by one, or by replacing one its features, ceases to make a significant difference in determining the score of the new species. However, since now the system only has very few species and species with large feature sets, the creation of new species which are better than those present prevails over the possibility of creating species with very different scores which are strong enough to coexist with the new species. 3. Middle-sized systems with low diversity: This is an intermediate situation between the two previous cases; it may be considered to be transitional behaviour between these two. It occurs when, at a relatively early

The introduction of the new speciation mechanism, as an alternative to the original mutational one, is in one way a more realistic assumption for the evolving nature of ecosystems. This is, however, a very simple and heuristic way of proceeding and, as we have noted, the model is sensitive to this novel speciation device. In order to produce a mature system which is able to exhibit diversity and an interesting topology, one must let this mechanism act in a slower way than the mutational one, i.e. the probability of adding a feature to a given species has to be much smaller than the probability of replacing one. It is interesting to note that in this case, for large enough times, it is possible to have systems of coexisting species with feature strings that possess very different lengths. As shown in Fig. 3 another observed pattern concerning the level species occupy in the resultant food webs produced by the model, consists of a tendency for species with fewer features to settle in the lowest level of the web, whereas those with more features tend to occupy higher levels. Since, as we have already stressed, features are not—at least in the current version of the model—mapped on to specific phenotypic or behavioural traits in real systems, this phenomenon cannot be interpreted in terms of real species at present. Others have been more ambitious, particularly with regard to the evolution of complexity in digital organisms, for instance models such as Tierra (Ray, 1991) and Avida (Adami and Brown, 1994). Furthermore, a similar approach to ours is discussed by Maron (2004). This latter study took the form of an exploration of the basic idea, rather than a systematic investigation. Having studied in some detail the response of the model to a more complicated definition of species, we now investigate the model in a direction which is in some sense opposite: by employing a simplified definition of species. 4. Dichotomous competition and species without features In this section we begin by introducing an intermediate version of the model, which involves a re-definition of the competition matrix which does not make use of species features, but where the model species are still defined by sets of features. The competition matrix assumed by Drossel et al. (2001) is given by Eq. (3) which, as discussed in Section 2, works under the premise that the strength of the competition to obtain resources from a common prey

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

655

55 50

Trophic level 1 Trophic level 2 Trophic level 3 Trophic level 4

45

Species Length

40 35 30 25 20 15 10 5 0

0

20

40 Species Label

60

80

Fig. 3. The length of the species feature set (a) and model food web (b) for typical outputs at the mature state of the community (105 time steps), where a new species enters the system as a mutant of its parent species with probability m ¼ 0:95 or as a more complex version of it, with probability of 1  m ¼ 0:05. In (a) the horizontal axis represents the species label, shown in (b) and it is assigned arbitrarily. The symbols represent the level on which each species are located in the food web. Species with more features tend to occupy higher levels than those with fewer features.

between individuals is larger if the degree of similarity between the competing species is greater. This requires a way of measuring similarities between species, which is given by the number of common features qij in Eq. (3). Here, in order to preserve the essence of the competition mechanism in the original, we propose a simpler method of modelling this, but which can be thought of as a limiting case of Eq. (3). The idea is to reduce the level of description of the competition mechanism to aij ¼ c þ ð1  cÞdij ,

(4)

which replaces Eq. (3). Here dij is defined to take the value 1 if i ¼ j and 0 if iaj. Therefore for inter-specific competition, the score takes the value c ð0oco1Þ, whereas the intra-specific score equals 1. We call this form of competition ‘‘dichotomous’’. Apart from this, all the other aspects of the model are retained. In particular, speciation takes place only by mutation-like events.

4.1. Dichotomous competition results The above modification of the model was studied mainly by investigating the response to changes in the value of c. As reported in Quince et al. (2005), the effect of varying the value of c employing Eq. (3) mainly consists of variations in the number of species allowed in the system; it decreases as c increases (Table 5). Therefore, for this version of the model, it should be expected that this effect will also occur and be more noticeable. However, now competition scores are no longer distributed throughout the interval ½c; 1, but only at the extreme values c and 1. Fig. 4 shows time-series for the number of species for several values of c, and parameter values equal to R ¼ 105 , l ¼ 0:1 and b ¼ 5  103 , which in the standard formulation of the model produce moderately large communities. In the case being considered here, the number of species present is considerably higher than in the original version of the

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

656

Competition parameter c=0.99

900

c=0.2

800

c=0.5

700

c=0.6

600

c=0.7

500 400 300

60 50 Number of species

Number of species

1000

40 30 20 10

200 100

0 0

0 500 1000 1500 Time (evolutionary time−steps)

Fig. 4. Response of the model to changes in c employing Eq. (4) to calculate the competition scores with parameter values R ¼ 105 , l ¼ 0:1 and b ¼ 5  103 .

model for almost every value of c; it is only for values in the range 0:8oco1:0, that the model produces food webs with sizes comparable to those reported for the original version. Table 2 summarises food web properties obtained for values of c in this range; the results represent the average values taken from around 50 independent realisations of the model. It should be mentioned that for values of c close to one, the statistics collected all belong to the case where the model was able to grow diverse communities. However, this was not the only observed scenario: there were also runs where the system did not grow in a reasonable number of time steps and cases where a mature system crashed into a configuration with very few species (see Fig. 5). The results obtained also show a trend in the number of links per species, which increases as the value employed for c decreases, which is in general larger than that obtained with the standard version for the same value of c. Table 6 shows results obtained with different versions of the model where the number of links and the average number of levels are similar, but when this occurs the number of species in the standard version of the model is greater than in the other cases. Recalling that the model uses a generalised version of the ratio-dependent functional response (Eq. (A.3)) to calculate the consumption rates, which are calculated in such a way as to fulfil the foraging strategy condition (Eq. (A.2)), it seems reasonable to believe that the criterion for link establishment, f ij 40:01, is satisfied by more pairs of species in the new version. This is because all the contributions in the second term of the sum in the denominator of Eq. (A.3) contain aki ¼ c which, as mentioned earlier, is the minimum allowed value for the competition scores in the standard version of the model. Fig. 6 shows typical model food webs obtained for values of c ¼ 0:95 and 0:99. They illustrate in a clear way both trends: the increase in the number of species and in the

0.5

1

1.5

Time (evolutionary time−steps)

2000

2 x 105

Competition parameter c=0.99 50 Average number of species

0

40 30 20 10 0 0

0.5 1 1.5 Time (evolutionary time−steps)

2 x 105

Fig. 5. The effect of employing a large value of the competition parameter is a decrease in the number of species present. However, it is not always possible to reach stable mature states. (a) Typical behaviours found employing c ¼ 0:99: formation of a diverse community, crash and no growth. (b) Average of the number of species for cases where the communities reach a mature state with several species in a time interval of 2  105 evolutionary time steps.

number of links per species as c decreases. Both web structures possess the same number of levels, but with a very different degree of complexity. As discussed in the next section this effect was also observed in the version of the model where species were not defined in terms of features.

4.2. Species without features Having defined a version of the model where the definition of competition is freed from the need for features to exist, we now introduce a version in which species features are not explicitly present at all. As described earlier, when a standard speciation event occurs, one feature of a randomly chosen species (denoted p, for parent) is replaced, creating a new species (denoted c, for child), which is related to all the other species by Eqs. (1)

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

657

Fig. 6. Food webs produced by the model employing the dichotomous competition score Eq. (4), for values of c ¼ 0:95 and 0:99. The radii of the spheres are proportional to the logarithm of species population sizes.

and (2), i.e. S cj ¼ maxf0; S~ cj g,

(5)

with S~ cj given by XX mab a2c b2j

(similarly for S~ jc ). Now S~ cj and S~ pj are related by ! X 1 X ~ ~ S cj ¼ S pj þ mnf b  m pf b , L b2j b2j

(6)

where nf and pf represent the new feature and the replaced feature of the parent species, respectively. The second term on the right-hand side of Eq. (6) is a random variable, and represents the change between the score of the parent species and that of the child species. We may therefore write Eq. (6) as S~ cj ¼ S~ pj þ Zj ,

(7)

but now choose Zj without making any use of species features. That is, we choose a parent species as in the original model, but now take the Zj ’s from some

probability distribution. The actual scores S ij are related to S~ ij through Eq. (2) as before. This choice of the scores, along with the competition score assignation discussed earlier, constitutes a formulation of the model which preserves its evolutionary nature, without making any use of internal degrees of freedom to characterise species. For this study, we employed Zj ’s taken from a Gaussian distribution with zero mean and unit variance. This choice for the evolutionary dynamics follows from studying in some detail the statistical properties of the scores in the original model (Lugo and McKane, 2008). There it was found that, in general, the scores exhibit distributions which are approximately Gaussian, but with ‘‘fat tails’’. The distribution has zero mean and variance very close to C, where C is connectance or fraction of nonzero elements in the FSM. In this paper we have taken C ¼ 1, which motivates the expectation that the distribution of scores will be approximately Gaussian with zero mean and unit variance. With the method of generating scores settled, the assignment of all the other quantities (efforts, initial population, etc.) for the new species is performed in the same way as the original model, where the initial values are inherited from the parent species where

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

658

possible, or they take an initial value which is randomly assigned in the case of the efforts, and equal to 1 in the case of the initial populations. Finally, the initial configuration assumed for this version of the model is almost the same as in the standard version (see Appendix A) with the only difference that the initial score S10 4b=l is a random number taken from the same distribution as the Zj ’s. 4.3. Results Table 3 shows the results obtained from the model using the rules given by Eqs. (4) and (7) to calculate the competition and species scores, respectively, and averaging the results over an ensemble of 50 webs. As expected, there is a trend in the number of species as the value of c changes (see Fig. 7). The community sizes are now more moderate in comparison to the ones shown in Table 2, and the number of time steps required to reach the mature configuration appears to be less than in the previous case. It should be noticed that the fraction of basal species is now larger than the fraction of intermediate species for values of c greater than or equal to 0.4 in the table. For c ¼ 0:2 (see Table 3), the model once again has a majority of species in the intermediate category. The number of links per species observed in this case does not differ too much compared to those obtained using the standard version of the model, but it is still slightly higher (see Table 6). Again this increases as the value of c decreases. Figs. 8 and 9 show a typical model food web obtained for a value of c ¼ 0:4. In Fig. 8(a) only intra-level links are shown, whereas Fig. 8(b) shows only inter-level links. Another relevant result of this exploration of the model is that it does not seem to produce food webs with more than three trophic levels, whereas in the standard version of the model the creation of webs with four levels is not unusual for the same set of parameter values. To produce food webs with more than three trophic levels, keeping the 350

Average number of species.

300 250 C=0.5

200

C=0.4 150

C=0.7 C=0.2

100 50 0 0

2000

4000 6000 8000 Time (evolutionary time−steps)

10000

Fig. 7. Time series of the number of species for several values of c and parameter values R ¼ 105 , b ¼ 5  103 and l ¼ 0:1.

statistics of the Zj ’s fixed along with the simple model of competition, some of the other parameters need to be modified. A natural choice is the ecological efficiency, l, since a more efficient flow of resources would be expected to allow species to be further away from the external resources. To illustrate this in the clearest possible way, a rather large value of l was chosen. It should be stressed that this value was taken to illustrate this effect, and we are not suggesting that l should in reality take on this value. Fig. 9(b) shows a typical food web obtained setting parameters R ¼ 105 , l ¼ 0:4, b ¼ 5  103 and c ¼ 0:7. This food web possess four trophic levels and features summarised in Table 4, along with results for some other values of c. These represent the average values obtained from a sample of around 40 webs. In all the cases studied for the latter parameter values, the model food webs typically produced contained four trophic levels and the number of links per species is larger than those in the standard version of the model (see Table 6). 5. Discussion and outlook We have explored the role played by species definition in one of the most studied evolutionary models of food web dynamics. Firstly, we allowed species to evolve by introducing additional features, along with the standard speciation by mutation. The results (Table 1) for this variant of the model showed that it is possible to produce stable mature communities in which species possess differing number of features. This occurs when the mutation-like speciation mechanism is allowed to act faster than the one which increases the number of features (Figs. 1 and 2). The structure of the resultant food webs remains very close to those produced by the standard version of the model, with the added result that species with a more complex structure tend to occupy the higher levels in their respective food chains (Fig. 3). Within the current formulation of the model it is not possible to compare this directly with patterns observed in real systems. A version of the model employing a dichotomous competition score and species definition with and without features was then introduced and studied in some detail. The main results (Tables 2–4) showed an increase in the size of the community in both cases, due to a global reduction in the inter-specific competition, and the production of webs with a slightly larger number of links per species than those produced in the original model (Table 5). In the case of the model employing features to calculate the species scores, the results relating to the food web structural properties are not very different to those produced by the standard version of the model. In the model without species features, it was also found that the number of time steps required to produce systems which reach a mature state was less than that required in the dichotomous competition version with features. The community sizes in this latter version are also more moderate in comparison to the dichotomous version with

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

659

Fig. 8. Model food web produced by the model without species features. (a) Food web structure omitting the intra-level links. (b) Subgraphs showing only links between intra-level species. The corresponding parameter values are R ¼ 105 , b ¼ 5  103 , c ¼ 0:4 and l ¼ 0:1.

Fig. 9. Multi-levelled food web produced by the model without features with parameters R ¼ 105 , l ¼ 0:4, b ¼ 5  103 and c ¼ 0:7.

features. However, some of the food web structural properties are only comparable to those obtained with the original model by increasing the value of some of the parameters, such as the ecological efficiency (Table 4). These results indicate that how species are defined have a moderate effect on the formation of communities. How-

ever, they have relatively less effect on the predictions of the model than the employment of adaptive foraging or the population and evolutionary dynamics do, which depend only in an implicit way on the species definition. So these novel versions of the model produce food webs which compare well with real food webs (see Table 7, for instance)

ARTICLE IN PRESS 660

C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

if we tune some of the parameter values, but keep the dynamical aspects intact. The two versions of the model using dichotomous competition with and without features, showed a slightly higher number of links per species than the original model. Whether this is more in line with data on real food webs is difficult to say. The agreement of the original model with data was most recently discussed in Quince et al. (2005), but the whole issue of such comparisons is fraught with difficulty (Drossel and McKane, 2003). The results shown in Table 7 for the number of links per species is taken from the 14 ecosystems studied by Dunne et al. (2002), and were the ones used for the comparison in Quince et al. (2005). The number of links per species spans a range between 1.51 and 13.12. Moreover, the definition of a ‘‘link’’ in models is just as problematical in models as in data, since typically there are very many weak links. This means that a change in the definition of the cut-off where flow of resources are large enough to be classified as a link makes a considerable difference in the number of links in the web (Quince et al., 2005). In the modelling of complex systems, such as food webs, there is always a dichotomy between including more aspects of ‘‘reality’’ into the model and simplifying it to make it easier to analyse. A prerequisite for the latter is that the essential structures that the model produces are not significantly altered because of the simplification. This is true of the modifications which we have discussed in this paper: eliminating features from the model still leads to the formation of non-pathological food webs. Having a version of the model which does not make any use of the species features, opens up the possibility of investigating and making a more direct contact between the Webworld model and other approaches. One of these is the model discussed by Tokita and Yasutomi (2003). There, a similar way of introducing species to that described by Eq. (7) was employed and generalised to mimic the processes of invasion (local and global), along with mutations. Therefore a version of the Webworld model where species could not just evolve from a species currently present in the community, but invade from an external community, could be investigated in a straightforward manner by following a similar procedure to that described by Tokita and Yasutomi (2003). There are even simpler stochastic models which the model can be compared to: where scores are chosen to be Gaussian random variables, but where there is no population dynamics (McKane et al., 2000) or neutral theories where all species are equivalent (Hubbell, 2001). These models can be systematically analysed, and although they are considerably simpler than the simplest model we consider here, it still is valuable to keep comparisons in mind. These models and the ones we have introduced in this paper form a sequence multispecies coevolutionary models which range from the complicated to the simple. We have argued that there is a continuity in the structure of the food webs they produce, which have the potential to allow any insights we might gain from the analysis of simple models

to still be potentially applicable to more complicated models. We hope to try to make these links more concrete in the future and so gain a deeper understanding of the modelling of predator–prey networks. Acknowledgements We thank C. Powell for useful discussions. C.A.L. acknowledges the award of a studentship from CONACYT (Mexico). A.J.M. wishes to thank the EPSRC (UK) for partial support under Grant no. GR/T11784/0. Appendix A. Webworld model: dynamical aspects While the aspects of the original version of the model (Drossel et al., 2001) which are modified in this paper are discussed in the Introduction and in Section 2, other aspects which are left unchanged will be discussed in this appendix. Further details are also to be found in McKane (2004) and McKane and Drossel (2005, 2006). A.1. Population dynamics Species obtain resources from their prey and give up resources to their predators. If we let gij ðtÞ represent the rate of consumption of individuals of species j by those of species i, then the updating of the species’ populations can be modelled in a general way through the balance equation: X X N_ i ðtÞ ¼ N i ðtÞ þ l gij ðtÞN j ðtÞ  gji ðtÞN i ðtÞ. (A.1) j

j

Here the first term represents loss of population from natural death, the second represents the amount of population gain per unit time due to predation, with the parameter l playing the role of average ecological efficiency. The last term represents the loss of resource due to predation by other species. The choice of the functional form for the gij leads to a particular population dynamics model. Eq. (A.1) represents a seasonal update of the population, so it is a population change that takes place at an intermediate timescale. However, in the Webworld model there is also an underlying process governing the dynamics of the gij : a foraging strategy, which occurs on a still shorter timescale. A.2. Foraging strategy We denote the fraction of the searching time that species i devotes to preying on species j to be f ij . We call this the effort that species i puts into preying on species j. One foraging strategy consists of assuming that the searching continues until the gain per unit effort gij =f ij that a predator puts into preying on any one of its prey are equal. It can be shown that this choice constitutes an evolutionarily stable strategy (ESS) (Drossel et al., 2001). Now gij has some functional form which will depend on the efforts

ARTICLE IN PRESS C.A. Lugo, A.J. McKane / Journal of Theoretical Biology 252 (2008) 649–661

f ij , and the efforts may themselves may be expressed as , X f ij ¼ gij gik , k

in order that X f ik ¼ 1.

661

The foraging strategy and the population updates are repeated until Eq. (A.1) reaches a stationary state. A speciation event is then carried out with the new individual being introduced into the community with an initial population of N initial , which here is taken to be 1.0. References

k

Then the equations for the ESS read gij f ij ¼ P , k gik gij ¼

S ij f N j P ij , bN j þ k aki Skj f kj N k

(A.2)

(A.3)

where the sum in the denominator runs over all predators k which prey on species j. The functional response, gij , used above is the generalised ratio-dependent functional response frequently used in the Webworld model. Here b and alm are, respectively, the saturation constant and the competition between species l and m, the latter being defined by Eq. (3). Eqs. (A.2) and (A.3) are iterated until a fixed point of the efforts is reached and the ESS condition fulfilled. Only then is an update of the population sizes performed using Eq. (A.1) with the revised rates. A species which has a population less than some threshold, N death , is removed from the system. Normally this value is taken to be N death ¼ 1:0. Finally, in reference to the way we define a link between two given species, we have already mentioned that the condition Sij 40 must be fulfilled for i to be a predator and j to be a prey. However, in addition, a second requirement relating to the efforts has to be imposed. We take this to be f ij 40:01. This choice means that very weak links are not counted as links, specifically, species that 1 devote less than 100 of their available search time hunting a particular prey, are not classified as predators of that species. This is in line with many other studies, both empirical and modelling, however, as discussed in Section 5, and in more detail in Quince et al. (2005), such choices are arbitrary and have important consequences when making comparisons with data. A.3. Species aggregation The last step in the evolution process is the introduction of new species into the system. The evolution starts from an original state which consists of two species, say, species zero and species one. Species zero plays the role of the environmental resources, so that species one has to able to feed from this environment. The feature set of species zero remains fixed in time, as well as its population size, which is given by R=l, where R is the rate of input of resources into the system. Species one is a randomly created species which should satisfy the condition S 10 4b=l, in order to be able to survive (Drossel et al., 2001).

Adami, C., Brown, C.T., 1994. Evolutionary learning in the 2D artificial life system Avida. In: Brooks, R.A., Maes, P. (Eds.), Artificial Life IV. MIT Press, Cambridge, MA, pp. 377–381. Caldarelli, G., Higgs, P.G., McKane, A.J., 1998. Modelling coevolution in multispecies communities. J. Theor. Biol. 193, 345–358. de Ruiter, P.C., Wolters, V., Moore, J.C., 2005. Dynamical Food Webs. Elsevier, Amsterdam. Drossel, B., McKane, A.J., 2003. Modelling food webs. In: Bornholdt, S., Schuster, H.G. (Eds.), Handbook of Graphs and Networks. WileyVCH, pp. 218–247. Drossel, B., Higgs, P.G., McKane, A.J., 2001. The influence of predator–prey population dynamics on the long term evolution of food web structure. J. Theor. Biol. 208, 91–107. Drossel, B., McKane, A.J., Quince, C., 2004. The impact of nonlinear functional responses on the long-term evolution of food web structure. J. Theor. Biol. 229, 538–548. Dunne, J.A., Williams, R.J., Martinez, N.D., 2002. Network structure and biodiversity loss in food webs: robustness increases with connectance. Ecol. Lett. 5, 558–567. Hubbell, S.P., 2001. The Unified Neutral Theory of Biodiversity and Biogeography. Princeton University Press, New Jersey. Kondoh, M., 2003. Foraging adaptation and the relationship between food web complexity and stability. Science 299, 1388–1391. La¨ssig, M., Bastolla, U., Manrubia, S.C., Valleriani, A., 2001. Shape of ecological networks. Phys. Rev. Lett. 86, 4418–4421. Lugo, C.A., McKane, A.J., 2008. The robustness of the Webworld model to changes in its structure. Ecol. Complexity, in press, doi:10.1016/ j.ecocom.2007.06.012. Maron, M., 2004. An ecological approach to the evolution of complexity. Master’s Thesis, University of Sussex. McKane, A.J., 2004. Evolving complex food webs. Eur. Phys. J. B 38, 287–295. McKane, A.J., Drossel, B., 2005. Modelling evolving food webs. In: de Ruiter, P.C., Wolters, V., Moore, J.C. (Eds.), Dynamical Food Webs:. Elsevier, Singapore, pp. 74–88. McKane, A.J., Drossel, B., 2006. Models of food web evolution. In: Pascual, M., Dunne, J.A. (Eds.), Ecological Networks: Linking Structure to Dynamics in Food Webs. Oxford University Press, Oxford, pp. 223–243. McKane, A.J., Alonso, D., Sole´, R.V., 2000. A mean field stochastic theory for species rich assembled communities. Phys. Rev. E 62, 8466–8484. Pascual, M., Dunne, J.A., 2006. Ecological Networks: Linking Structure to Dynamics in Food Webs. Oxford University Press, Oxford. Quince, C., Higgs, P.G., McKane, A.J., 2005. Topological structure and interaction strengths in model food webs. Ecol. Model. 187, 389–412. Ray, T.S., 1991. Evolution and optimization of digital organisms. In: Billingsley, K.R., Derohanes, E., Brown, I.H. (Eds.), Scientific Excellence in Supercomputing: The IBM 1990 Contest Prize Papers. The Baldwin Press, The University of Georgia, pp. 489–531. Rossberg, A.G., Matsuda, H., Amemiya, T., 2006. Some properties of the speciation model for food-web structure: mechanisms for degree distributions and intervality. J. Theor. Biol. 238, 401–415. Tokita, K., Yasutomi, A., 2003. Emergence of a complex network in a model ecosystem with extinction and mutation. Theor. Popul. Biol. 63, 131–146. Yoshida, K., 2003a. Dynamics of evolutionary patterns of clades in a foodweb system model. Ecol. Res. 18, 625–637. Yoshida, K., 2003b. Evolutionary dynamics of species diversity in an interaction web system. Ecol. Model. 163, 131–143.