Tellus (1999). 51B. 815-829 Printed in UK - ail rights reserved Copyright <8 Munksgaard, 1999 TELLUS ISSN 0280-6495 Food web complexity enhances community stability and climate regulation in a geophysiological model By STEPHAN P. HARDING, Schumacher College, The Old Postern, Darlington, Devon TQ9 6EA, England, UK (Manuscript received 7 August 1998; in final form 25 February 1999) ABSTRACT A central debate in community ecology concerns the relationship between the complexity of communities and their stability. How does the richness of food web structures affect their resistance and resilience to perturbation? Most mathematical models of communities have shown that stability declines as complexity increases but so far, modellers have not included the material environment in their calculations. Here an otherwise conventional community ecology model is described, which includes feedback between the biota and their climate. This "geophysiological" mode! is stable in the sense that it resists perturbation. The more complex the community included in the model, the greater its stability in terms of both resistance to perturbation and rate of response to perturbation. This is a realistic way to model the natural world because organisms cannot avoid feedback to arid from their material environment. 1. Introduction Alfred Lotka (1925), one of the founders of population biology, suggested that the dynamics of populations cannot be considered in isolation from the material environment in which they are embedded and with which they interact. He saw organisms and environment as a single evolving system in which each influences the other in a relationship of mutual feedback. He suggested that the mathematics of the evolution of this whole system would be more tractable than that of either population or environment taken separately. Thus, so far most theoretical ecologists have not heeded Lotka's advice. Instead, they have modelled communities without considering the feedbacks between organisms and their chemical and physical environment. Such work has given fruitful insights into the dynamics of populations in which coupling with the material environment can be considered so weak as to be effectively non-existent (May, 1981; Pimm, 1979a; Pimm, 1979b). Complexity is often defined in terms of the presence of more species, stronger inter-species interactions and greater connectance between species (Begon et al., 1996). Using these definitions of complexity, modellers have generally found that simple model communities are more likely to be stable than complex ones, apart from the exceptions by (De Angelis, 1975). Pimm's (1984) measures of stability, resilience and resistance are used in the model described here. Resilience is taken to mean the speed with which a community returns to a former state having been displaced from it by perturbation, whilst resistance is thought of as the ability to avoid such displacement. Stable communities are those with high resilience and resistance. Resilience and resistance are not necessarily positively correlated. For example, a highly resistant system, when eventually displaced, may not be resilient because it requires a long time to return to steady state. In this paper resilience is measured by the rate of return of a community to steady state after a perturbation. The familiar result, of complexity leading to instability in ecosystem models, is not inevitable Tellus 51B (1999), 4 816 S. P. HARDING (Johnson et.al., 1996). For example, complexity can lead to stability in food webs where the food supply is not afTected by consumers (De Angelis, 1975), or in food webs perturbed by the removal of species in lower trophic levels (Pimm, 1979a), or where non-linear per capita growth rates, non-equilibrium dynamics and realistic feeding interactions are incorporated (Polis, 1998; McCann et al„ 1998). A qualitatively different class of strongly stabilised models are those that include environmental feedback (Watson and Lovelock, 1983). In these, instability is the rare condition' In this paper, experiments with a community ecology model are described which incorporate explicit mathematical feedbacks between biota and environment. In this class of model greater complexity gives rise to greater stability. 2, Daisyworld The numerical model, daisyworld (Lovelock, 1983; Watson and Lovelock, 1983; Lovelock,' 1992) was used as the basis of the experiments' (Fig. 1). The original daisyworld is a mean field model of an abstract flat surface "sown" with seeds of 2 daisy types, one dark and the other light, A model sun warms the planet. Daisies germinate and grow between 5°C-40°C with peak rates at 22.5°C, following a Gaussian curve. Once daisies begin to grow, they compete for space and as the available space decreases, so does their growth rate. The interaction between solar luminosity and the planetary albedo sets the temperature of the planet. Plants with low albedos (darker types) absorb more solar radiation than would the bare surface and thus warm both themselves and the planet. Those with high albedos (i.e., whiter types) reflect solar radiation back to space and therefore cool both themselves and their environment. This simple model system responds to changes in solar luminosity and keeps planetary temperature well within the limits tolerable by life even when solar luminosity varies over a wide range. The model is unusually stable to changes in initial conditions and to perturbation. Temperature regulation via the biota's albedo in daisyworld is a particular example of the general biota-environment feedbacks of the real world. The mathematical -foundations of daisyworld are increasingly being used in other Earth-system models, such as Elapsed Time toiy£^£S^j££ rtht0n,y 2 da'Sy tyPeS' °ne dafk and °"e "S11'- ^ the "-rsity, of the s(cadjly iaTZnS 7ot1Z:J^ emPeratUrC Cl°Se l° thC °PtimUm f°r dais* growth over a wide range Tdlus 51B (1999), 4 FOOD WEB COMPLEXITY IN A GEOPH YSJOLOGICAL MODEL 817 those of Lovelock and Kump (1994) and Lenton (1998). A frequent criticism of daisyworld is that it is a toy model and does not apply to the real world. However, a similar vegetation albedo feedback on climate in the boreal forests has been modelled in a manner inspired by daisyworld (Betts and Cox, 1997). Daisyworld is useful and practical because its stability and self-regulating properties do not depend on a particular mathematical formulation. Daisyworld's growth equation was taken from a deterministic plant population dynamics model (Carter and Prince, 1981). Furthermore, Maddock (1991) has shown that self-regulation around a steady state arises in a variety of well known population ecology models where there is tight coupling between biota and their environment. 3. Methods 23 different daisy types were added to the original daisyworld model, each with a unique albedo ranging between 0.2 (dark) and 0.775 (light). Daisy-eating herbivores were also included, each type preying on daisies according to one of three classic functional responses, depicted in Fig. 2. (Rolling, 1959). Previous studies have looked at the impact of these three herbivore feeding styles acting separately (Harding and Lovelock, 1996), but in this model they were introduced together and competed for food and space. A carnivore, which consumed the herbivores according to a type 3 functional response could also be introduced. Neither herbivores nor carnivores could change their colours, i.e., their albedos, 0.5 0.4 0.3 0.2 0.1 / / 0.2 0.4 0.6 0.8 1.0 M to O c o '€ o Q. o 0.5 0.4 0.3 0.2 0.1 0 0.2 0.4 0.6 0.8 1.0 Daisy abundance 0.2 0.4 0.6 0.8 Fig. 2. The functional response types which could be allocated to the three herbivores. The four figures show how uncnonal response configuration was systematically varied in order to explore the model's response to this manipula-/iil' A , lolIowmg F values were substituted into the type 1, type 2 and type 3 functional response eqs.(il) (12) (13)::m) type 1: F = 1; type 2: F = 3; type 3: F = -20. Iffc^c1^"U^TT= 3;7ypeS^-^1 = F - 0.5, type 2: F = 1.5; type 3: F = -3. (2d) type 1: F ^03T5F^'T^l3:lyFJ:'T^ -15.--------- Tellus 51B (1999), 4 818 S. P. HARDING and so they had no direct feedback effects on climate. However, like the daisies, their growth rate was temperature dependent. The equations used in this model were based on those of the original model (Lovelock, 1992): re = (SJ(l-A)/a)0 25 - 273, (1) A =xAt + (aiAi + a2A2 + ... + akAt), (2) 7J = 9{J4-J4i) + T„ /J,=(l -0.00326(22.5- 7J)2), daJdt^ciiixPi-yi), x = 1 — {«! + a2 + ... + at), (3) (4) (5} (6) where % = the effective temperature of the planet and A is its average albedo, and J represents the solar luminosity relative to the actual luminosity, taken as 1.0. The stellar constant is S, and a is the Stcfan-Boltzmann constant. The total area of the planet is taken as unity, the fraction of the planet which bare ground is x,. and _the_albede^. of this bare ground area is ^g.(Thfi.fraction of the pianef covered by each of the k daisylype is af, and the constant albedo of each daisy type is At. The temperature of each daisy type is 7;, and q is a constant, 20, representing the re-distribution of solar energy. Finally, each daisy type has a constant death rate, y, (set at 0.15) and a temperature dependent, helmet shaped growth rate, /?,, optimum at 22.5°C. The quadratic growth equation (4) was replaced with a more biologically realistic Gaussian temperature dependent growth curve according to the following equation: ft = exp(H22.5-7i)20.01). I (7) The equations for including higher trophic levels followed eqs.{8) to (10) in Lovelock (1992), but with two modifications. Firstly, herbivore and carnivore growth rates (jih and /?c, respectively) depended only on, the Tc overall temperature of the system, according to the following equations: )3b = exp(-(22.5- rj20.01), &=exp(-(22.5-Te)20.01). 1 (8) 3 (9) The second modification was to introduce the three classical functional responses (Holling, 1959) to model density-dependent impacts of herbivores on daisy populations and of the carnivore on herbivore populations, where A (plants), B (herbi- vores) and C (carnivores) ^efer to populations of each trophic level, and (x, y] z) to the space unoccupied by these populations. The modified growth equation for daisy types was: dajdt = a,l0tx - y, - (B^puX,) + B2(PllX,) + B3(P3i^H. (10) ^ where X, is determined by one of the eqs. (11), (12) or (13), and p is a couling matrix defined as: 5 where ph 6 (0, 1). If p,j = 0, daisy j was not eaten by herbivore i, whereas if p;j = 1, the daisy was eaten. Here X-t represents the functional response for a given herbivore, where X; modified the impact of Bir the herbivore population, according to the \- following equations, in which F determined the slope of the functional' response, up to a plateau of 0.5: For a type 1 response: X{ = a,F, up to a plateau of X( = 0.5. (11) ( For a type 2 response: i^d-expf-fl.-F))^, (12) -? For a type 3 response: ^, = (I-exp(flf-f))/2. (13)J: The values of F used in the experiments are given in the legend of Fig. 2. The equation used to model the growth of the three herbivores was a version of Lovelock's (1992) equation (9), modified to incorporate kt, the density dependent predation from the carnivore, when present, and Git a growth increment, itself a function of Ej, the total amount of daisy biomass consumed by the herbivore. The new equation was as follows, where i varies from 1 to 3, corresponding to the three herbivores, and where the herbivores' death rate was held constant at 7, = 0.3: dVd( = W,-y + Gi-y,.-CA,.), (14) y where: G,= 1 -exp(-1.02£(), and: Xi — I — exp(—20B2). (15) /, a (i6) /; Te!lus51B(1999),4 FOOD WEB COMPLEXITY IN A. GEOFHYSIOLOGICAL MODEL 819 d£(/dt = B,. Y PiA Mi The total amount of daisy biomass for each herbi- due to natural selection amongst daisies for anti-vore was calculated as follows: herbivore responses such as distasteful secondary compounds which are known to influence plant-(17)/ ^ herbivore dynamics in the field (Huntly, 1991). The stabilising effect of environmental feedback where i refers to herbivores and n refers to daisy makes possible large models involving many evol-lypes, and: ving non-linear differential equations. The price for this advantage is the need to explore the Ej = fe + B; (£pi„/t„)d(, (lg) extended space occupied by the model variables, J Q and it was decided to explore the effects of varying where it was a constant set at k = 0.0001, used to connectance patterns. There were 300 possible initialise the herbivore population. Similarly, Lovelock's (1992) eq. (10) for describing the growth of the carnivore population was modified to incorporate G, the carnivore's growth increment, a function of E, the total herbivore biomass eaten by the carnivore, with its death rate held constant at y = 0.05: dc/dt = c(B(tz+G-y), where: G= 1 -exp(-1.02£), where E is the total amount of herbivore biomass eaten by the carnivore, calculated as follows: ways to link the 23 daisy types with the three herbivores, depending on the number of daisy types eaten by each herbivore. This was in the absence of carnivores and for all food webs where each herbivore was allocated the same type 2 functional response, with slope F = 2.2, and a plateau of 0.5. The connectance c was set as either c = 1 or c = 2. Each of these 300 "food web (19) ^configurations" specified the numbers of daisy types eaten by each herbivore in a given food web. The triangular diagram (Fig. 3) illustrates the part (20) , of the model space occupied by each food web d£/d£ = C 3 I 4 dr, where: A; = (l — exp(-Finally: E = k + C Q»d( (22) j (23) where k was a constant, set at k = 0.001, used to initialise the carnivore population. A critical parameter in the experiments reported here was food web connectance, c, which was varied from the maximum of c = 3, where each daisy was eaten by each of the three herbivores, to the minimum of c= I, where each daisy was eaten by only one randomly selected herbivore. The connectance parameter, c, specified the values (0 or 1) in the coupling matrix, p. Intermediate values of c, that is, "c = 1 or 2" and "c = 2" denoted that each daisy was eaten by one or two randomly selected herbivores, or by two randomly selected herbivores respectively. These low (c < 3) food web connectances can be thought as having arisen configuration. This space was explored systematically by a computer program that stepped through a self-similar portion of the total space shown in Fig 3. Linking the three herbivores with randomly chosen daisy types made connectance patterns. The program produced three of these patterns for each food web configuration. For all 3 patterns all variables were held constant whilst solar luminosity, J, was increased in 20 equal increments, from J = 0.61 to J — 0.93. The time taken for the diversity index to return to its steady state was Aised to measure the food web's resilience. This was observed for each luminosity step after the introduction of the three herbivores, and where appropriate, after the subsequent introduction of - the carnivore. Resistance, the ability of the system to resist change, was measured by the amplitude of the temperature deviation from the optimum growth temperature of 22.5"C after the experimental perturbations. To calculate the stability of the model the means of the 20 resilience and resistance measurements of each pattern were taken and then the grand means of the three replicates were calculated as the final measures of resilience and resistance. The effects of variations in herbivore functional response type and slope on the behaviour of the model were explored by assigning the three herbi- Teilus51B(1999), 4 820 S. P. HARDING Number of daisy types eaten by Herbivore 2 or Herbivores I & 3 Fig. 3. (a) The parameter space of all possible food web configurations for c = 1 and c = 2 food webs, in the absence of carnivores, and with each herbivore allocated the same functional response, To obtain data for the entire parameter space it was sufficient to sample the top left sub-triangle (in heavy line}, and then to extrapolate the results to the remaining space. This is because the food web configurations within each sub-triangle are self-similar in terms of herbivore impacts on daisy types. vores in a given food web either identical or different functional responses, of high or low slope, and then measuring resilience. As before, three replicate connectance patterns were generated for c = 1 and c = 1 or 2 webs, with food web configuration held constant (see legend of Fig. 7 for details). In addition, the effects of perturbing the model with a 5% step increase in solar luminosity were explored. As before, the functional response curves of Fig. 2b were used to generate these experiments. Five c = 1 webs and five c = 2 webs, each with a different randomised connectance pattern, were constructed, each with constant food web configuration (see legend to Fig, 9 for details). Each web was allowed to reach steady state under a low initial solar luminosity and was then perturbed. The time taken to return to a steady state daisy diversity index was then rneasured; a rapid return meant high resilience. These resilience measurements were repeated for each of 20 equal increments in solar luminosity. For each web, the mean Fig. 4. Population dynamics and temperature regulation in food webs with different values of c, the connectance parameter, amongst 23 daisy types, 3 herbivores and one carnivore. Herbivores were introduced after 40 time steps and a carnivore after 120 time steps. The herbivores' functional response configuration was that of Fig. 2b, in which herbivore 1 was allocated the type 1 response, herbivore 2 the type 2 response and herbivore 3 the type 3 response, (a) Dynamics of a fully connected food web (c = 3). Solar luminosity was held constant at J = 0.744. (b) Dynamics of a more loosely connected food web (c= 1 or 2), with the food web configuration: herbivores 1, 2 and 3 ate n = 14, n = 12 and n = 14 randomly assigned daisy types respectively. Solar luminosity was held constant at / = 0.744. (c) Dynamics of a food web with the lowest connectance (c = 1), herbivores 1, 2 and 3 ate n = 8, n = 8 and n = 7 randomly assigned daisy types respectively. Solar luminosity was held constant at J = 0.826. Tellus51B(1999),4 FOOD WEB COMPLEXITY IN A GEOPHYSIOLOGICAL MODEL 821 1_ 0.6 Temperalure Sasped Time 6|-- c l/'-v W\ rv'V aA\ Temperature Daisy Type 1 herbivore Type 2 hertuVofe Daisy i\^M\[\i\i\i\f\ f i 1 I 1' 11 li 1 \ > [I l' > l; i »' 4 1! * II ! Type 3 herbivore 3 herbivore types Daisies 100 150 200 Elapsed Time Tellus SIB (1999), 4 822 S. P. HARDING of these 20 measurements was taken as an index of its resilience. 4. Results Fig. 4a shows a typical example of the population dynamics and temperature regulation for a fully connected (c = 3) food web, with solar luminosity held constant at a level comfortable for the biota. The figure illustrates the stability of models incorporating environmental feedback. Less stable behaviour was found with models where the food webs were more loosely connected over a range of solar luminosities, from 0.7 to 1.2 times the present solar output. Model stability was not much altered by changes in herbivore functional responses but was sensitive to different food web configurations and connectance patterns, although temperature regulation was still well within the bounds tolerable by life. At the extremes of solar luminosity just before the system failed it became unstable. There were rapid oscillations in daisy and herbivore populations and in temperature. This was particularly noticeable with steep type 1 and 2 functional responses, but no attempt was made to analyse these odd effects. As the system in Fig. 4a evolved, it was dominated by 2 daisy types whose albedos best matched the given solar luminosity. The population dynamics were stable and temperature regulation was at the optimum. Upon adding the three herbivores their grazing reduced the abundance of the few dominant daisies, and allowed the less well matched rarer types to flourish. Daisy diversity increased with a smooth gradation of daisy type abundances, with no single type overly dominant. Competition between the herbivores led to the dominance of a single type. Upon introducing the carnivore the dominant herbivore became its prey and so the other two herbivores increased until all three reached equal abundance. These changes in herbivore abundance had only a small effect on the daisies although a few types became extinct. The model system kept its temperature close to optimal for the biota throughout these perturbations. Herbivores and the carnivore had little effect on daisyworld's capacity to regulate climate. When the food web connectance was decreased from c = 3 to c = 1 or 2, the system was less stable {Fig. 4b). Although the herbivore dynamics were similar to those in the fully connected web, the daisy dynamics were markedly different. After the introduction of the herbivores, 4 daisy types reached high abundance, with a roughly even gradation in the abundance of less common types. The appearance of the herbivores led to a temperature increase to a steady state about 1.0°C higher than the model optimum for plant growth of 22.5°C. In this experiment, one of the herbivores again became dominant. This in turn decreased the grazing of daisy types eaten by its competitors and these increased. Among the now more abundant daisy types were those less well matched to the solar luminosity and therefore, temperature regulation was impaired. Adding the carnivore improved the model performance. It reduced the abundance of the dominant herbivore, allowing the others to increase until the 3 types were equally abundant. The abundant daisy types that had previously flourished due to competitive suppression of their herbivores now declined. There was again an even gradation of daisy type abundances, and those with best matched albedos predominated, thereby improving temperature regulation. Models with the most loosely connected (c = 1) webs were markedly less stable, At certain solar luminosities and in the absence of the carnivore they sustained high amplitude periodic dynamics. When this happened (Fig. 4c), the populations of two herbivores regularly oscillated out of phase with each other, while the third herbivore's numbers gently oscillated in step with one of the other herbivores. The oscillations of the model system also affected the temperature, which was lightly modulated at the same frequency and over a range of 0.8°C. The mean temperature remained constant. These periodic fluctuations emerged unexpectedly in c = 1 webs. There appeared to be no relationship between the emergence of periodicity and solar luminosity (Fig. 5). This form of instability seemed to reflect the randomness in the underlying connectance patterns. Sustained periodic dynamics were not seen in the more fully connected food webs, where c> 1. Adding the carnivore to c = 1 food webs obliterated periodic dynamics where they had been present (Fig. 4c). The presence of the carnivore brought the system to a steady state. The three herbivore types then became equally abundant, there was an even Tellus51B(1999),4 FOOD WEB COMPLEXITY IN A GEOPHYSIOLOGICAL MODEL 823 high amplitude cycles low ampltide cycles no cycles • «^« • * ». ^ »♦♦^ ••• v«v * ♦ < • n / \ w A •-• •-•-•-•-•-•-♦-•-•V V- • ♦ ♦ * • * • ^ • * »-W #-«-• V • « 0.6 0.96 Solar Luminosity Fig. 5. The frequency of occurrence of high amplitude cycles (as in Fig. 3c), low amplitude cycles and no cycles in five randomly assembled c = I food webs. There was no relationship between solar luminosity and the presence high amplitude cycles, which appear in roughly 25% of solar luminosity steps gradation in the abundances of daisy types and smooth temperature regulation. The model displayed stability once more, as it did with the more connected food webs. Even with the least stable, loosely connected daisyworlds the departure of the temperature from optimum was seldom more than 2.0°C over a wide range of solar luminosities. Fig. 6 shows that despite some scatter in c = l webs, resilience and resistance showed a statistically significant, positive correlation throughout the parameter space of c = l and c = 2 food webs. Thus, in what follows, changes in resilience can be understood to imply similar changes in resistance, and hence in the overall stability of the model. Despite wide variations in critical initial conditions of the model (i.e., in food web configuration, connectance pattern and slope and type of the herbivore functional responses), low connectance generally gave rise to low resilience (Figs.r>, 8). However there were minor exceptions. Firstly, in food webs where the three herbivores followed type 3 functional responses (Fig. 7b), c= l webs were approximately equally resilient to c ~ 1 or 2 webs, irrespective of the slope of the functional response. This was due to the fact that type 3 functional responses had very little impact on rare daisy types, amongst which were types whose albedos matched a given luminosity well enough to engender good temperature regulation, thereby improving resilience. Secondly, food webs varied in their resilience depending on their "herbivore connectance diversity". This quantity was highest in food webs in which each herbivore was linked to the same number of daisy types. For example, high herbivore connectance diversity would occur in a c = 1 food web if two of the herbivores ate 8 daisy types each and the third herbivore ate the remaining 7 daisy types. Low herbivore con- Tel!us5IB(1999), 4 824 S. P. HARDING 4 160 120 food webs, r =0.303, p<0.05 □ c=2 food webs, r2=0.531, p<0.001 the sct of samP>ed < = ' ^ ^ebs, i , ,h °f;amPledfc=u2 food webs- The two variables show a s.gnificant positive correlation in both cases although the relat.onsh.p for the c = 1 food webs shows a relaüvely high degree of scatter ' nectance diversity would occur in c = 1 webs if herbivore 1 ate all 23 daisy types, with the other two herbivores eating no daisies. Fig. 8 reveals that c = 1 webs were markedly less resilient than c = 2 webs over a wide range of herbivore con-nectance diversities. Furthermore, in c = 1 webs, despite some scatter in the data, it was clear that resilience continually increased as herbivore con-nectance diversity declined. Thus, for c = 1 food webs, return time to steady state decreased as herbivore 1 became more prominent in a food web. By contrast, c = 2 webs displayed an extensive region of roughly constant high resilience, covering almost all of the parameter space, with far less variability in the data. In the real world climate, change is a frequent perturbation. Increasing the solar luminosity by 5% simulated a change of roughly equivalent magnitude to that between glacial and interglacial climates. Fully connected food webs showed high restlicnce, reaching steady state almost instantly after the perturbation (Fig. 9a-c). The presence of carnivores made no difference to resilience, but at some solar luminosities did slightly improve temperature regulation (Fig. 9d-f). In the absence of the carnivore and in food webs with the lowest connectance (c = 1), resilience fell to low levels. In these webs, at some luminosities, the climate perturbation set off a long term periodic oscillation in population and in temperature. The oscillations did decay but only after many time steps {Fig. 9g i). Adding a carnivore improved resilience, although the return time to equilibrium was longer than in more highly connected food webs with a carnivore present (Fig. 9j-l). At other luminosities, the perturbation did not give rise to periodic behaviour in c = 1 webs, although both population and temperature changed and then took a relatively long time to return to steady state. Adding the carnivore again improved resilience. Resilience was least for c = 1 food webs with carnivores absent, and increased when the carnivore was present, or when connectance was increased (Fig. 10). 5. Discussion So far, almost all model experiments in community ecology have been made with communities where the only selection pressures were internal, that is where the environment was defined by the TeI!us51B(1999), 4 FOOD WEB COMPLEXITY IN A OF.OPHYSIOLOGICAL MODEL 825 Type 1 only Type 3 only 80 60 40 20 80 H) E £