Abstract

The interest of environmental management is in the long-term health of populations and ecosystems. However, toxicity is usually assessed in short-term experiments with individuals. Modelling based on dynamic energy budget (DEB) theory aids the extraction of mechanistic information from the data, which in turn supports educated extrapolation to the population level. To illustrate the use of DEB models in this extrapolation, we analyse a dataset for life cycle toxicity of copper in the earthworm Dendrobaena octaedra. We compare four approaches for the analysis of the toxicity data: no model, a simple DEB model without reserves and maturation (the Kooijman–Metz formulation), a more complex one with static reserves and simplified maturation (as used in the DEBtox software) and a full-scale DEB model (DEB3) with explicit calculation of reserves and maturation. For the population prediction, we compare two simple demographic approaches (discrete time matrix model and continuous time Euler–Lotka equation). In our case, the difference between DEB approaches and population models turned out to be small. However, differences between DEB models increased when extrapolating to more field-relevant conditions. The DEB3 model allows for a completely consistent assessment of toxic effects and therefore greater confidence in extrapolating, but poses greater demands on the available data.

1. Introduction

The interest of environmental management is in the long-term health of populations and ecosystems. However, toxicity tests at these higher levels of organization are costly, impractical, difficult to interpret and for some species ethically unjustifiable or even prohibited. Toxicity testing therefore usually comprises experiments with a (more or less) homogeneous cohort of a selected test species. These experiments are (relatively) short term and conducted under controlled laboratory conditions, often following standard test protocols (e.g. OECD 1984). To extrapolate the toxic effects at the individual level to meaningful consequences at the population level, model approaches are essential.

Populations in the field may show complex dynamic behaviour in response to toxic stressors, depending on (time-varying) environmental conditions, interactions with conspecifics and presence of other species (competitors, predators and parasites). However, environmental management requires simple and general summary statistics at the population level. Under constant environmental conditions (and excluding intra- and interspecific interactions and density effects), all populations will eventually grow exponentially, with rate constant r (the asymptotic population growth rate). This statistic is related to another popular statistic, the population growth factor λ (technically not a rate), as r = (log λ)/Δt (where Δt is the time interval implied for λ). The influence of toxicants on life-history traits (survival and reproduction) can directly be expressed as a decrease in r. Even though it is clearly unrealistic to expect prolonged exponential growth in real populations, r reflects the inherent capacity of a population to bounce back from a drastic decline in numbers. Therefore, r can be considered a general fitness measure for the population, and its use in environmental risk assessment has been advocated (Forbes & Calow 2002).

For single-celled organisms such as algae and bacteria, population models are often unstructured, leading to simple expressions for the toxic response at the population level (e.g. Kooijman et al. 1996). However, population-level questions for most multicellular organisms require structured populations, in which life-history characteristics differ with age or size. Two popular approaches for structured population predictions in response to toxicants are simple matrix models (individuals divided in discrete classes and discrete time) and the Euler–Lotka equation (traits as continuous functions of age and continuous time), discussed in more detail later in this paper.

Toxicity tests can yield data on survival and reproduction over (a part of) the life cycle, which can be used directly to parameterize a population model (e.g. Bindesbøl et al. 2007). However, such a descriptive analysis only provides a projection of the observed data to population growth rates. It does not provide information on the underlying processes of the toxic response, and therefore no extrapolation to field conditions can be made. Far more insight can be obtained by combining population modelling with a biology-based analysis of the toxic effect at the individual level. Dynamic energy budget (DEB) theory provides a solid framework to deal with toxic effects, treating them as a disruption of regular metabolic processes such as maintenance and assimilation (Kooijman & Bedaux 1996; Jager et al. 2006). A DEB-based treatment of toxicity data is specifically relevant for population projections as it provides a causal link between body size and reproduction, which helps to estimate vital rates from observations. Furthermore, a DEB approach allows for educated extrapolations of effects for untested concentrations and other environmental conditions. Moreover, a DEB approach also makes it possible to estimate the interaction of multiple stress factors (Jager et al. 2004; Klok 2008). As shown by Kooijman & Metz (1984), the population response to toxicants under food limitation depends critically on which specific metabolic process is affected by the compound.

To shed some light on the pros and cons of various DEB models and population approaches, we take an existing dataset for life cycle toxicity of copper to the earthworm Dendrobaena octaedra. This dataset originates from Bindesbøl et al. (2007), who also made population projections using a simple two-stage model. Klok et al. (2007) used the same dataset with a four-stage matrix model after analysis with a simple DEB model. In this contribution, we will analyse the same data with more complex DEB models and calculate r using a matrix model as well as the Euler–Lotka equation. We focus our discussion on the estimation of r and thus on constant environmental conditions.

2. Deb approaches for analysing toxicity data

An overview of the three DEB approaches used in this study is given in table 1; the corresponding equations are provided in the electronic supplementary material.

(a) Kooijman–Metz formulation

Kooijman & Metz (1984) published one of the first versions of a DEB model (here referred to as KM) and discussed the impact of toxicants on individuals and the consequences for the population level. This approach incorporates the idea of a ‘mode of action’; i.e. toxicants affect a particular metabolic process that leads to a concerted change in life-history traits. The KM formulation does not consider reserves or maturity maintenance, does not include toxicokinetics for sublethal effects (DEB parameters under stress are constant in time) and assumes no particular relationship between exposure concentration and the underlying DEB parameters (the data at each exposure concentration are treated separately). The latter implies that inter- or extrapolation to other exposure concentrations is not possible within this approach. Furthermore, the paper by Kooijman and Metz is a theoretical one; no fits to experimental toxicity data are shown. Klok & De Roos (1996) tested the sublethal part of the KM model for copper effects on earthworms. They simplified the KM formulation for unlimited food conditions (as is standard practice in toxicity tests), and only used five compound DEB parameters (such as the Von Bertalanffy growth rate constant and the maximum reproduction rate). Depending on the physiological mode of action, these compound parameters change in a systematic manner under stress. For example, if stress influences assimilation, this decreases maximum length and maximum reproduction rate by the same factor.

Overview of three DEB approaches for ecotoxicity, as discussed in this paper.

For our comparison with the other DEB models, we slightly extend the approach of Klok and De Roos by assuming a threshold for effects (no-effect concentration, NEC) and a linear relationship between the external concentration above the threshold and the effect on the DEB parameters (see electronic supplementary material).

(b) DEBtox

The DEBtox approach (Kooijman & Bedaux 1996; Jager et al. 2006) follows from two basic assumptions: a chemical needs to be taken up to cause effects and the internal concentration affects one (or more) primary DEB parameter(s) in a linear fashion above an NEC. Five basic modes of action were originally proposed (Kooijman & Bedaux 1996). The toxicants are assumed to act on the primary DEB parameters, but the model is formulated in terms of compound parameters such as maximum body length and maximum reproduction rate. For toxicokinetics, a first-order one-compartment model is used, accounting for effects of growth (dilution and change in surface–volume ratio). DEBtox applies a simplified DEB model to accommodate the limited datasets that result from standard toxicity tests; reserve density is constant, reproduction starts at a fixed body size and egg costs are treated as a primary parameter (not logically linked to the DEB parameters for the adult). These simplifications lead to some inconsistencies. For example, an effect on assimilation is not consistent with a constant reserve density and will affect the egg costs (Jager et al. 2010).

(c) DEB3

An extended DEB formulation was recently presented by Kooijman et al. (2008), which we will refer to here as DEB3. In contrast to DEBtox, this approach includes explicit calculation of maturation, reserves and egg costs (Kooijman et al. 2008; Sousa et al. 2010). DEB3 ensures complete consistency when parameters change under the influence of toxicants. For example, toxicant stress on some of the primary parameters will affect length at puberty and egg costs. The approach taken to treat toxicant effects in DEB3 was presented by Jager et al. (2010). To allow an estimation of all model parameters from data on growth, survival and reproduction, the DEB3 model is phrased in the dimensions length and time only (see electronic supplementary material).

3. Population models for toxic effects

(a) Simple matrix models

Matrix models (Caswell 2001) are popular to interpret effects of toxicants at the population level. Populations are structured by assigning individuals to classes based on age, size or developmental stage. Individuals within each class are assumed to be exactly the same. Depending on age or size, individuals can progress to the next class or die with a certain probability and produce offspring in some classes. These transfer probabilities and reproduction rates are collected in the matrix form (the ‘projection matrix’) and taken constant in time (in the simplest matrix approaches). Time is taken in discrete steps (the ‘projection interval’), and the projection matrix is applied to current population structure to calculate the number of individuals in each class for the next time point. Using discrete time steps and discrete classes implies that matrix models are particularly suited for species with a discontinuous reproduction pattern, such as birth-pulse populations. Depending on the type of data available, matrix entries may relate in a simple way to observations (e.g. monitoring data), which are usually collected on a discrete time basis. Both the asymptotic population growth rate and the stable stage/age distribution can be calculated analytically (Caswell 2001).

Inevitably, matrix models will introduce errors due to discretization by using a discrete time step and a subdivision in classes, especially for animals with a more-or-less continuous production of offspring in time (birth-flow populations). These continuous time processes need to be discretized, which implies loss of information. The magnitude and relevance of this loss will depend on the particular life cycle of the organism of interest, the number of classes and the choice of projection interval.

(b) Continuous Euler–Lotka

A discretization of population structure into matrix form is not always needed. When continuous age-specific survival and reproduction functions are available (or can be estimated from the observations), the continuous form of the Euler–Lotka equation can be used to calculate the asymptotic population growth rate (e.g. Kooijman & Metz 1984; Jager et al. 2004). An advantage over matrix models is that there is no error from dividing a continuous life cycle into discrete classes and using discrete time. Because survival and reproduction are generally observed at certain points in time only, some sort of interpolation is required. As a crude interpolation, the data in time can be splined (as discussed in §6b), but the use of DEB models provides interpolation (and even extrapolation) possibilities on a more mechanistic basis.

Compared with matrix models, calculating the population growth from the Euler–Lotka equation is less straightforward; the equation is not explicit and requires numerical integration. However, with the current power of desktop computers, this is no longer a strong argument against this method. The Euler–Lotka equation also allows for calculating elasticities and the stable age/stage distribution (De Roos 2008).

4. Description of the dataset

Newly hatched juveniles (2–4 days old) of the earthworm D. octaedra were exposed to copper-spiked natural soil (0–200 mg copper added per kilogram dry soil). Worms were kept for 20 weeks under laboratory conditions with ad libitum food (ground cow dung). Development stage, fresh weight, survival and cocoons were determined regularly (full details in Bindesbøl et al. 2007). The experiments started with, on average, 3 day old juveniles that had been starved after hatching. We accommodated this setup in all re-analysed data by taking t = 0 at the point of hatching, in contrast to the original authors. The DEB3 formulation includes dynamic calculation of reserves, and therefore, we simulated the effects of starvation over the first 3 days. For the KM and DEBtox analyses, we assumed no growth and no exposure for the first 3 days.

Growth data for all treatments revealed an s-shaped pattern when body size is expressed as a length measure (cubic root of body weight). This contradicts DEB theory, which under constant environmental conditions leads to Von Bertalanffy growth. Working in a DEB framework, this discrepancy requires a mechanistic explanation. Similar s-shaped growth curves have been observed in bacterivorous nematodes, which were caused by size-dependent feeding limitation; for juvenile nematodes, the bacterial food supplied in laboratory tests was large relative to the size of their mouthparts, thereby reducing feeding efficiency (Jager et al. 2005). It is not clear whether a similar mechanism is acting for Dendrobaena, but we assumed the same relationship between body size and scaled ingestion rate as applied for nematodes. This assumption required one additional parameter (Lf, see electronic supplementary material) and resulted in a much-improved data fit (figure 1). The goodness of fit across all treatments strongly supports our assumption of an initial feeding limitation in this test situation.

Interestingly, low doses of copper appeared to have a stimulatory effect (hormesis) on growth and reproduction (highest performance at 80 mg kg−1 copper in soil). Based on the same dataset, Klok et al. (2007) found curved relationships between copper and DEB parameters, confirming combined positive and negative effects of copper. The hormetic response does not hamper the data analysis when no dose–response is assumed, as every exposure concentration is treated in isolation. However, in a fully consistent DEB framework, the hormetic response requires a mechanistic explanation; the extra investment in growth and offspring needs to be paid from either an increased feeding rate or a decrease in other expenses (e.g. maintenance costs). Even though copper is an essential micronutrient for earthworms, it appears unlikely that levels were limiting in the test soil (a loamy sand containing 11 mg kg−1 copper). Copper is widely used as a fungicide, so it is more likely that low doses reduce the negative impact of fungi, which may be particularly abundant when manure is used as food source (Klok et al. 2007). A similar explanation was suggested by Svendsen et al. (2005), who observed a stimulatory effect of the antiparasitic drug ivermectin on Lumbricus terrestris. We therefore think it is likely that this low-dose stimulation is an indirect effect through food (quality or quantity) and modelled its effect through the ingestion rate (as indicated by Spurgeon et al. 2003). Quite arbitrarily, we selected the dose with the highest performance as representative control (80 mg kg−1). This means that worms at lower copper concentrations (0 and 40 mg kg−1) were assumed to suffer from food limitation (f < 1, table 2). The observed effects of low copper concentrations on the growth and reproduction curves are completely consistent with this assumption (figure 1), thereby supporting our explanation.

Copper had a large effect on hatching success at the highest test concentration. We described hatching success simply as a linear function of exposure concentration above a threshold (hatching success in control 83.5 per cent, threshold for effect on hatching 151 mg kg−1 and tolerance concentration 77.7 mg kg−1, see electronic supplementary material). There was no dose–response apparent in the survival data (likelihood ratio test using the DEBtox model), and, therefore, mortality in all treatments was considered to be background mortality.

5. Deb analyses of the dataset

All DEB calculations were made with Matlab 7.3, using maximum likelihood optimization (confidence intervals were derived using profile likelihood). Here, we only discuss the analysis with DEB3 in detail. First, we focus on the control response. The three lowest concentrations (0, 40 and 80 mg kg−1 copper) were treated as controls at different food levels. Growth and reproduction at several food levels should be sufficient to identify all DEB3 parameters (Kooijman et al. 2008). The dataset, however, contains limited information on reproduction (only three points in time), and we included an additional parameter for size-dependent food limitation, which requires us to make several simplifications. In this dataset, body length at first reproduction is largely independent of food level, so we fixed k at a value of 1 (Kooijman et al. 2008; Sousa et al. 2010). Additionally, we assumed that the batch of eggs to start the experiment was produced by mothers from ad libitum food conditions (f = 1) and fixed the energy investment ratio (g) to a value of 1 (Kooijman & Bedaux 1996; model behaviour turned out to be rather insensitive to the exact value of this parameter). With these assumptions, the DEB model fitted the control data well. Parameter estimates are given in table 2, and fits are shown in the upper panels of figure 1.

A consequence of these parameter values is that DEB estimates a cocoon development time of less than a week, whereas Holmstrup et al. (1991) report an average hatching time of 92 days for this species at this temperature on filter paper. Cocoon development time is, however, highly variable and depends on season, soil moisture level and temperature (Edwards & Lofty 1977; Holmstrup et al. 1991) and even varies strongly between cocoons produced by the same individual (Svendsen et al. 2005). The development time of 92 days is very long, relative to the life cycle of the organism (longer than the time needed to grow from hatchling to reproductive adult). Therefore, it is likely that the embryo spends substantial time in diapause and that development requires a change in environmental conditions to start, thus ensuring that the juvenile appears under favourable conditions. The work of Holmstrup and co-workers indicates that a sudden increase in temperature triggered hatching. Nevertheless, the origin of this discrepancy requires further study, especially because hatching time (and the influence of stressors on it) can be important for population dynamics.

To estimate direct effects of copper (above 80 mg kg−1), we fixed the physiological parameters to their maximum likelihood estimate (table 2). The effects of copper on growth and reproduction are explained by assuming that copper affects the ingestion rate (figure 1, model fits; table 2, estimates for toxicological parameters). This mode of action was selected because it fitted the data best, although an increase in maintenance costs was only slightly worse. The elimination rate indicates rapid toxicokinetics (not significantly different from instantaneous equilibrium), which is consistent with observations on copper kinetics in another earthworm species (Peijnenburg et al. 1999). It should be noted that taking exposure to 80 mg kg−1 copper as the control could have underestimated toxicity.

The DEB model predicts a slight increase in egg development times with increasing copper concentrations. An increase in hatching time with copper exposure was indeed observed for other earthworm species (Bengtsson et al. 1986; Helling et al. 2000), but the validity of this prediction requires dedicated testing.

For the analyses with DEBtox, we also assumed that the response at 80 mg kg−1 is the representative control and included Lf. In line with the DEB3 analysis, the mode of action is through the ingestion rate (fits and parameter estimates in the electronic supplementary material). The DEBtox fits are very similar to the DEB3 fits in figure 1. Similarity results from the fact that the mode of action does not interfere with maturation and size is a good proxy for maturity in this species (k = 1). Furthermore, this similarity shows that, at least for this exposure scenario, reserves dynamics hardly influence the results.

A KM analysis of this dataset is given in Klok et al. (2007). In comparison to the other two DEB approaches, we additionally made a calculation with the KM formulation including a dose–response relationship between external concentration and the scaled ingestion rate (see electronic supplementary material). Model fits and parameter estimates (with and without feeding limitation for the juveniles) are provided in the electronic supplementary material. With feeding limitation, the growth curves for the KM approach are similar to those of the other DEB models, as the KM growth equation is essentially the same as the DEBtox formula. However, larger deviations are shown for reproduction. KM does not include maturity maintenance, which leads to straighter curves for the cumulative reproductive output in time and a larger size at puberty.

6. Calculating asymptotic population growth rate

In the sections below, we will compare several approaches for population projection based on this dataset. Section 6a will introduce our reference case: an Euler–Lotka calculation based on a DEB3 analysis. In §6b, population projections will be based directly on the experimental data, and in §6c, we discuss the combination of a DEB analysis with matrix calculations. Section 6d will compare the influence of the choice for a DEB approach by comparing our three approaches using the Euler–Lotka equation only.

(a) DEB3 combined with Euler–Lotka

The Euler–Lotka equation is calculated with the reproduction rates and survival probabilities as a function of age, as specified by DEB model parameters (Jager et al. 2004). Because exposure concentration is a parameter in the DEB model, population growth rate (r) can be plotted as a continuous function of concentration (the lines in figure 2). For these population projections, hormesis is treated as an indirect effect of the experimental setup (as explained in §4) and therefore irrelevant for field populations (the response at 80 mg kg−1 copper is taken as a representative control). A hatching time of 92 days is assumed for all treatments and 3 days of initial starvation for the hatchlings. The DEB analysis indicated a shorter hatching time, but a diapause is probably relevant in this species' life cycle. Furthermore, this hatching time allows for a direct comparison with the analyses of Bindesbøl et al. (2007) and Klok et al. (2007), who also assumed 92 days.

Estimated population growth rate using the continuous Euler–Lotka equation with DEB3, integrated over 200 weeks (dashed lines) and 20 weeks (solid line). In plots (a) and (b), symbols represent the results for the matrix calculations: open circles, recalculated from Bindesbøl et al. (2007); filled circles, recalculated from Klok et al. (2007); filled triangles, calculations based on DEB3 analysis. In plot (c), filled squares represent splined data over the 20 weeks of the experiment, used in the Euler–Lotka equation.

For the Euler–Lotka equation, a choice has to be made about the time interval of integration. In principle, one would have to integrate until all animals are dead or until reproduction ceases (although the contribution of reproduction to the population growth rate decreases with age). However, the dataset in this case only provides information about life-history traits for some 20 weeks after hatching. One option is to include observed behaviour only, so integration over the test duration only. This is equivalent to all of the animals dying or ceasing reproduction at the end of the test, which is clearly unrealistic. Another option would be to integrate over a longer time span than the test duration, which would require an extrapolation of survival and reproduction behaviour outside the range of observations. It is also unrealistic to assume that reproduction and survival rates remain unchanged after 20 weeks of the test until the animal dies, but this is what is done in stage-structured matrix models.

For our comparisons, we made two calculations with the Euler–Lotka equation: one integrated over 20 weeks only (solid lines of figure 2) and one over a lifespan that is considerably longer than the test duration (200 weeks), assuming that reproduction and mortality rates remain constant (dashed lines in figure 2). This latter calculation is much more comparable to the calculations made with matrix models. Figure 2a shows that, in this case, the choice of integration range for the Euler–Lotka equation influences r considerably. This implies that the assumptions on the survival and reproduction patterns after 20 weeks do matter (increasing lifespan for more than 200 weeks had very little additional effect).

Without DEB analysis, the experimental data can be used directly in population models, both simple matrix models (e.g. Bindesbøl et al. 2007) and continuous time Euler–Lotka equation. Bindesbøl and co-workers used the life-history traits as observed in each treatment directly in a two-stage model. The population growth factors (λ) they present were recalculated to growth rates and plotted in figure 2a. For the Euler–Lotka projections with DEB3 in the same plot, the hormetic response at low copper concentrations was ignored. This explains the constant r for low concentrations, in contrast to the peaked relationship for the calculations of Bindesbøl et al.

For the Euler–Lotka calculation without DEB analysis, continuous survival and reproduction functions were derived by interpolating the survival and reproduction data in time using splines (Alda Álvarez et al. 2006). Here, we interpolated the data for each treatment using a piecewise cubic Hermite interpolation and for hatching success assumed a linear relationship with threshold as explained earlier. We integrated over a 20 week lifespan only, as it is impossible (or at least unreliable) to extrapolate data beyond the range of the observations without a mechanistic model. These continuous splines are used in the Euler–Lotka equation. Additionally, we included the 3 day initial starvation of the hatchlings (in line with the DEB3 analysis, but in contrast to the results of Bindesbøl et al. and Klok et al.). However, the effect of these additional 3 days on r is negligible. The spline estimates are given as symbols in figure 2c, which shows a fair correspondence to the DEB3 analysis for a 20 week adult lifespan. Note that without a process-based description of the toxic effect, it is not possible to estimate r outside the range of tested concentrations.

(c) KM and DEB3 combined with a matrix model

Klok et al. (2007) analysed this dataset by combining the KM model with a four-stage matrix model formulation (Klok & De Roos 1996). Furthermore, they extended their analysis to estimate variability in r by an approach that takes covariability of the compound DEB parameters into account by using Bayesian statistics. Survival was based on background survival as estimated for field-relevant conditions. Here, we recalculated r with the same approach, but with the survival rate taken from the fit to the data (h0 in table 2). The initial starvation period is ignored, and the hatching success per treatment is used as reported by the original authors. The results are shown in figure 2a. It should be noted that the results are plotted as symbols and not as a continuous function of concentration. This is not a restriction due to the use of matrix models (e.g. Billoir et al. 2007), but results from the DEB approach of Klok and co-workers, which does not assume a dose–response relationship. The differences with the results from Bindesbøl et al. mainly result from different approaches to estimate the start of reproduction (Klok et al. 2007).

Additionally, a matrix calculation was performed using the length at puberty and reproduction rates as resulting from the DEB3 analysis (table 2). A three-stage matrix model was used, ignoring the subadult stage (which is not identified as a stage in DEB) and using a projection interval of 7 days. Results are plotted as symbols in figure 2b. Interestingly, the estimates from the matrix model correspond very closely to the Euler–Lotka calculations using the same DEB3 analysis, integrated over 200 weeks. Only the lowest two concentrations deviate, which is caused by the fact that the hormetic response was ignored for the Euler–Lotka calculations. For the matrix calculations, the life-history characteristics for these two treatments were estimated with the lower food levels specified in table 2, leading to lower values of r.

Both the matrix and Euler–Lotka calculations based on the DEB3 approach produce lower estimates for r than the other two matrix approaches (figure 2a). In DEB3, the age at maturation and reproduction rate are determined by the model parameters (table 2), whereas in the other approaches, it is determined directly from the observations. Clearly, the details of the reproduction curve can have a substantial effect on r.

(d) Comparison of KM, DEBtox and DEB3 with Euler–Lotka

In this section, we compare three different DEB approaches, regarding their prediction for the population response: DEB3, DEBtox and KM (with and without food limitation for the juveniles). For this comparison, we use the Euler–Lotka calculation, integrated over the time period of the observed response (20 weeks). The DEB3 analysis is in each case plotted as a solid line for reference in figure 3. Figure 3a shows that all models yield very similar predictions for the population growth rate under the conditions of the test (f = 1). Apparently, as long as the model provides a reasonable description of the reproduction and survival patterns, the details do not matter.

Estimated population growth rate using the continuous Euler–Lotka equation. (a) Population response based on the fit of several DEB-based approaches (all related to the DEB3 predictions) and an extrapolation to a lower food level. (b) Prediction without presumed artefacts (no food limitation in juveniles and no initial 3 day starvation after hatching) at three food levels. Thick lines represent response at high food (f = 1) and thin lines at limiting food levels (f < 1).

Using a DEB model to analyse the response to toxicants allows for extrapolation, for example, to exposure under lower food availability. This is done by changing the scaled ingestion rate (f) from 1 to an arbitrary lower value. Results are shown by thin lines in figure 3. Figure 3a shows that food limitation leads to a lower population growth rate and also to an increased effect of copper (the same percentage decrease in r is reached at a lower external concentration). If copper acts on assimilation, it acts ‘synergistically’ with food limitation, as predicted earlier (Kooijman & Metz 1984; Jager et al. 2004). Again, all approaches perform remarkably similarly, apart from the KM analysis without food limitation in juveniles, which shows less effect of decreasing food level. The reason for this behaviour is that assuming size-dependent food limitation in juveniles also makes this stage more sensitive to a decrease in the available food, lowering f even more.

As explained in §4, the s-shaped growth curves and the initial starvation after hatching are quite likely specific consequences of the experimental setup and therefore not relevant for field populations. Running the DEB model without these factors leads to higher values of r (figure 3b) and less severe effects of copper. Furthermore, the differences between the predictions of the various models increase. The DEBtox approach is still very similar to the DEB3 analysis, but the KM approach (especially without Lf) deviates more seriously. Thus, correctly interpreting the relevance of the s-shaped growth curve and initial starvation for field situations is very important for predicting effects at the population level. Furthermore, the choice of the model starts to become more important when stronger extrapolation is performed.

Several other factors, which were not included, may need adjustment for a more field-realistic estimation. We did not consider transfer of contaminants from the mother to the egg, or toxic effects on hatching time and hatchling size, which can have an important impact on r. Finally, under field conditions, mortality rates will likely be much higher than those in the laboratory test, e.g. due to disease and predation.

7. Discussion

Simple population models are useful to integrate toxicant effects on all individual endpoints over time into a more meaningful summary statistic: the asymptotic population growth rate (r). When the interest is in r under the specific laboratory test conditions only, it can be calculated without a DEB model or any other explanatory model at the individual level. Data for survival and reproduction in time can simply be splined and r calculated from the Euler–Lotka equation (figure 2), or data can be directly used in a simple two-stage calculation (Bindesbøl et al. 2007). However, when reproduction data are widely spaced in time (such as in our dataset, see figure 1), interpolation is dangerous as small deviations in estimates for the start of reproduction and reproduction rate substantially influence r (figure 2).

A DEB analysis provides a range of advantages over analysis without a mechanistic model, as emphasized in our example.

— DEB models make strong predictions for the life history of individuals, which implies that small experimental deviations are ignored and suggests particular experimental peculiarities that are probably irrelevant for the field situation (in our case, hormesis and the s-shaped growth curves, see §4). Furthermore, reliable estimates for the start of reproduction and reproduction rate can be derived.

— When a dose–response is assumed between the (internal) concentration and a DEB parameter, r can be plotted as a continuous function of time (figures 2 and 3), facilitating interpretation. Furthermore, an NEC can be estimated, which also implies no effects for the population.

— DEB analysis provides insight into underlying processes affected by the toxicant. In this example, copper effects are most likely caused by a reduction in ingestion (equivalent to a reduction in assimilation efficiency).

— Useful extrapolations can be made (figure 3), such as the consequences of experimental artefacts and the interactions between toxicant stress and food limitation.

In this paper, we compared the behaviour of three DEB-based approaches using a single dataset only, which limits the generality of our results. The DEB3 formulation is obviously the most extensive one, allowing for a completely consistent analysis of toxic effects on the energy budget. However, the simpler DEBtox produced almost exactly the same fit to the data and very similar values for the population growth rate. Although the initial starvation and size-dependent feeding limitation would require reserve dynamics for a consistent calculation, the effects are in this case minimal (reserve dynamics is apparently fast). The KM with juvenile food limitation fits the growth curves very well, but provides a poorer fit for the reproduction data (see electronic supplementary material). Nevertheless, the population prediction hardly suffers from this misfit (figure 3). The KM approach without juvenile food limitation provides a relatively poor fit to the body size data, but a rather good fit to the reproduction data. The values for r in figure 3 at f = 1 are therefore hardly distinguishable from the DEB3 analysis. However, the extrapolations to other food levels and to field populations differ strongly because the data are interpreted in a different way. As a consequence, the KM analysis without Lf predicts that body weight will still increase by more than a factor of 3 after 20 weeks of the test, whereas the other DEB analyses show that the final weight is almost achieved by the end of the test duration. In this example, the toxicokinetics for copper appear to be rapid (table 2), which explains why the KM approach (assuming instantaneous steady state) is able to provide a good description of the dataset. Although it is our experience that sublethal effects are often well described by rapid toxicokinetics, this is not always the case.

Ultimately, the choice for a particular DEB model is limited by the available data. A DEB3 analysis requires a strong dataset. Even though the dataset in our example is a fair representation of a typical life cycle toxicity test, it is insufficient to identify all DEB3 parameters accurately. In this case, various DEB approaches lead to similar conclusions about the population effects of copper for this earthworm species. Only for strong extrapolations do differences occur. Clearly, a complex DEB model is not needed for all applications. However, several cases would require a DEB3 approach:

— when the length at puberty varies with stress (e.g. Alda Álvarez et al. 2006), an explicit evaluation of maturation dynamics is needed;

— when food availability varies sharply in time, or in cases of periodic starvation, explicit calculation of reserve dynamics may be needed;

The differences between the population approaches were rather small in our example. When simple matrix modelling was compared with the Euler–Lotka equation on a similar basis, they produced similar results. The discretizations had a small influence on r, which probably relates to the fact that the reproduction rate of the earthworms did not change much with age and the background mortality was assumed to be constant.

In general, the choice for a matrix model or using the Euler–Lotka equation should depend on the life history of the organism, the available data and the problem that needs to be addressed. The standard DEB model with continuous functions for survival and reproduction in time makes the Euler–Lotka equation a more natural choice. However, for birth-pulse populations, matrix models may be the preferred approach. A matrix model may be an acceptable simplification when the reproduction rate changes little with age and size. However, with the current power of personal computers, the complexity of iteratively solving the Euler–Lotka equation is hardly an issue anymore.

It should be stressed that our analysis deals with the asymptotic population growth rate only. The extension to population dynamics under time-varying conditions is not straightforward, as most simple population models deal with a single state variable for the organism. DEB theory can be pragmatically reduced to a single state (Nisbet et al. 2010), but toxic stress will often require an additional one (body residues), as will changing food conditions in the environment (reserves).

8. Recommendations

— The choice for a particular DEB model should be based on available data, observed effect patterns and the extrapolation needs. More complex DEB approaches may not yield a better fit, but could provide more insight into the data and therefore allow for better extrapolation beyond the conditions of the laboratory test. For deriving basic DEB3 parameters, life history at different food levels is very useful.

— Endpoints such as egg size, hatching time and hatchling size are not usually determined in toxicity tests. However, this kind of information may influence population calculations. If these endpoints change under stress, a DEB3 approach is required to interpret them.

— Transfer of chemicals from mother to offspring should be assessed. Due to this transfer, offspring may experience different body residues and therefore different effects over time.

— A field-relevant value for r requires field-relevant background mortality as a function of age, at least up to an age where reproduction still contributes to r.

— When hormesis or s-shaped growth curves (on body length) are observed, the experimental setup should be critically examined (especially the food source, which may be suboptimal for some life stages).

— The choice for a particular matrix model or the Euler–Lotka calculation should depend on the life history of the organism, the available data and the problem that needs to be addressed. For organisms with a more-or-less continuous reproduction pattern, the Euler–Lotka calculation is a natural choice.

Acknowledgements

We wish to thank Anne-Mette Bindesbøl and Martin Holmstrup for providing the original data from their experiments. This study was partly supported by the European Union (EU) Marie Curie Initial Training Network CREAM (http://cream-itn.eu), Contract 238148.

2007Integrating the lethal and sublethal effects of toxic compounds into the population dynamics of Daphnia magna: a combination of the DEBtox and matrix population models. Ecol. Modell.203, 204–214. (doi:10.1016/j.ecolmodel.2006.11.021)

2007Extending a combined dynamic energy budget matrix population model with a Bayesian approach to assess variation in the intrinsic rate of population increase. An example in the earthworm Dendrobaena octaedra. Environ. Toxicol. Chem.26, 2383–2388. (doi:10.1897/07-223R.1)