Introduction

Many environments contain toxic compounds at levels that are too low to induce an acute lethal response, but may impact the dynamics of populations by adversely affecting the physiological performance of individual organisms. We develop a general modeling framework to describe sublethal physiological responses within the context of dynamic energy budget (DEB) theory (Kooijman 2000). Our aim is to develop and test simple formalism that is applicable to a wide range of species and contaminants, the rationale being that there are so many different species, combinations of potentially toxic compounds and habitats, that context-specific models are impractical in all but the best documented studies.

Toxic effect modeling within the framework of DEB theory involves three process-based modules (Billoir et al. 2008; Kooijman et al. 2008a, b; Kooijman and Bedaux 1996; Kooijman 2000). The first is a DEB model that describes the ingestion and assimilation of food by individual organisms and its utilization for maintenance, growth and reproduction (Kooijman 2000; Muller and Nisbet 2000). The second module describes how organisms exchange toxic compounds with the environment. The last module describes how toxic compounds that have been accumulated in the body subsequently affect the performance of organisms through changes in the values of the DEB model parameters.

DEB theory offers substantial advantages over purely descriptive methods, such as those applied to estimate values for the NOEC and LC50, since the interpretation of toxicity assessments from process based models are independent of experimental protocol, i.e. toxicity parameters estimated with the DEB framework are independent of exposure time and choice of toxicant test concentration (Alda Alvarez et al. 2006; Billoir et al. 2008; Jager et al. 2004). Moreover, DEB theory provides a framework within which different endpoints can be compared in a meaningful way, since those endpoints are based on shared processes in DEB theory (Alda Alvarez et al. 2006). Furthermore, with the DEB approach and information about several organisms and toxicant combinations, we can make educated guesses about other combinations via quantitative structure-activity and body-size scaling relationships (Kooijman et al. 2007).

In principle, a toxicant may impact any DEB model parameter. Kooijman and co-workers have argued that many toxicants are likely to have to reach a threshold level, the no-effect body burden, before an organism experiences adverse effects, and that every model parameter has its own threshold level (Alda Alvarez et al. 2006; Kooijman et al. 2008a, b; Kooijman and Bedaux 1996). This may imply that at relatively low contamination levels only one parameter is affected. Kooijman’s approach rests on this assumption, and has been tested with data from standard OECD toxicity tests (Billoir et al. 2008; Kooijman and Bedaux 1996) and a variety of other data (Alda Alvarez et al. 2006; Jager et al. 2004; Lopes et al. 2005). Generally, data are adequately described, but it is often impossible to determine which parameter is affected, as alternative submodels commonly fit equally well.

In this paper we consider an alternative representation of the sublethal effects of toxicants whose effect mechanism may be poorly understood. We recognize that organisms generally spend energy to offset the negative effects of toxicants and that toxicants commonly impair the organism’s ability to acquire resources from the environment. Hence, our toxic effects modules assume that toxicants primarily affect the parameters determining the rates of feeding and maintenance. We compare predictions of our models with a variety of published data. Among the pollutants tested are metals (mercury and copper) and various organic compounds (chlorophenols, toluene, polycyclic aromatic hydrocarbons, tetradifon and pyridine); organisms include mussels, oysters, earthworms, water fleas and zebrafish.

Models

Toxicants affect the vital rates of individual organisms. To model these effects, we need to know how an animal grows and reproduces in absence of toxicants, how it accumulates toxicants, and how toxicant interactions can affect the vital rates. We use three modules to describe these processes. As the first module, we use Kooijman’s DEB model to describe growth and reproduction of an ectothermic heterotroph that does not change its shape during growth (Kooijman 2000, 2001; Muller and Nisbet 2000). Since a detailed discussion of model assumptions and derivation can be found elsewhere (Kooijman et al. 2008a, b; Kooijman 2000), we confine this presentation to the bare essentials. Figure 1 outlines the energy and material flows and state equations, Table 1 lists the equations used in this paper, and Table 2 explains symbols and conventions.

Fig. 1
figure 1

Material and energy flows and primary state variables in a heterotrophic metazoan according to DEB theory (Kooijman 2000). The compositions of reserves and structure are constant, implying that conversion efficiencies for assimilation and growth, and the maintenance rate for a unit of structure are constant. Food is ingested at a rate proportional to the surface area of the organism and the scaled food density (type II functional response). Under non-starvation conditions, a constant fraction of mobilized reserves is used for somatic maintenance and growth, with maintenance having priority over growth; the remaining fraction is used for maturity maintenance and maturation (juveniles) or reproduction (adults), with maturity maintenance having priority over maturation and reproduction

Full size image
Table 1 Equations used in this paper
Full size table
Table 2 Symbols
Full size table

The DEB model distinguishes three types of biomass in a heterotrophic organism: structure, reserves and biomass set aside for reproduction (including sperm and eggs). Structure is defined to be biomass needed to perform vital functions; this type of biomass requires a continuous investment of energy in order to remain functional. Reserves are biomass that is neither structure nor reproductive matter, including reserve compounds stored in tissues, such as lipids and glycogen, and intracellular compounds used for immediate physiological demands. Each type of biomass is assumed to have a constant chemical composition, which implies that it takes a constant amount of energy to produce a new unit of structure (growth), and a constant rate of energy expenditure to maintain the viability of a unit of structure (somatic maintenance). As originally formulated, the model has two state variables: the volume of structure and the amount reserves expressed as a fraction of the structural volume. At constant food densities, this second state variable has a constant value (see Fig. 1). The data we analyze in this paper were measured at food densities that were nearly constant or with food supplied sufficiently frequently that it is reasonable to assume that the reserves adequately buffered the fluctuations in food densities. Then, the dynamics of the model closely follow the dynamics at constant food density (Billoir et al. 2008).

The second module describes the dynamics of the body burden (see Eq. 6 in Table 1). Since DEB theory accounts for a changing lipid content, which depends on feeding and spawning histories, the dynamics of the body burden are potentially complex. However, if ambient conditions are constant, the density of reproductive material is relatively stable, and the partitioning of toxicants is rapid relative to toxicant exchange with the environment, the dynamics reduce to the simple one-compartment model of toxicant exchange (Kooijman and Bedaux 1996). But even this simplified representation is often too complex for analyzing existing data on toxic effects, since toxic effect studies rarely include measurements on body burdens. Therefore, we make use of two properties of this model: the body burden in a previously unexposed organism initially increases approximately linearly and the body burden has an asymptotic maximum (Kooijman and Bedaux 1996). The first property implies that the body burden is approximately proportional to the ambient toxicant concentration and exposure time in relatively short experiments. The second property implies that the ultimate body burden is proportional to the ambient toxicant concentration. Consistent with both properties, we assume that the body burden is proportional to the ambient toxicant concentration when analyzing data lacking measurements of the body burden. Several toxicants, such as cadmium and PCBs, accumulate too slowly to justify this simplifying assumption. We therefore do not use this assumption for data obtained with these types of toxicants.

Organisms may take up toxicants without being appreciably affected in performance. For instance, metals such as copper and cadmium, bind to metallothionins and thereby become toxicologically inactive. As another example, many organic pollutants are altered in the body via Phase 1 and Phase 2 biotransformation pathways (Walker et al. 2006). Biotransformation yields products that are eliminated more rapidly and are often, with notable exceptions, less toxic than the original pollutant. Consistently, we assume that there is a maximum toxicant level in the body and environment up to which an organism will not experience adverse effects and call this level the no-effect body burden and no-effect concentration, respectively. These species-specific threshold values may not differ significantly, in the biological and statistical sense, from zero.

The third module specifies sublethal toxic effects. Although the biochemical modes of action of toxicants are manifold, toxicants can affect only a small number of DEB parameters, of which eight are relevant to this paper (see Table 2). Because we strive for maximum simplicity, we will consider only parameters that can be clearly identified as targets of toxicants: the maximum specific searching rate, {F m}, the maximum specific feeding rate, {J XAm }, and the maintenance rate, [p M ]. As a result, compound parameters that contain those parameters, i.e. the von Bertalanffy growth rate and the asymptotic length, are also affected by toxicants. The abilities to search for and handle food items decline as motor skills are impaired by neurotoxins, for instance. Since effects on the ability to search for food coincides with those on food handling, {F m } and {J XAm } are likely to be affected in a comparable way by toxicants. Consistently, we assume that the food searching and handling parameters decline hyperbolically with the effective body burden, i.e. the body burden in excess of the no-effect burden. Furthermore, we assume that the toxicity parameters for food searching and handling are identical. Because the ratio of these two parameters is defined as the half-saturation constant in the functional response (see Eq. 1), this implies that the latter is not affected by toxicants. The corresponding toxic effect functions are listed in Table 1 (see Eqs. 5, 6). Maintenance expenses are likely to increase in the presence of toxicants threatening the physiological integrity of an organism. Organisms regulate their physiological activities in order to create a stable internal environment. When this stability is at stake, organisms try to compensate for the disruptive effects of toxicants. We assume that the maintenance requirements of an animal increase linearly with the effective body burden (see Eq. 7 in Table 1).

We use maximum likelihood methods to fit our model data (Bolker 2008; Hilborn and Mangel 1997; McCullagh and Nelder 1989). To simplify calculations, we assume that a measured value is the sum of a deterministic process, as defined by our model, plus a normally distributed error term with a variance that is potentially a function of the independent variable. With these assumptions, the maximum likelihood estimates for the parameters are simply the values for which the weighted residuals are minimized. To determine confidence intervals, we use the (biased) maximum likelihood estimator for the variance. We usually need to make additional assumptions in order to use maximum likelihood methods. For data sets with unknown standard deviations we assume that the coefficient of variation is constant; we weigh those data with the reciprocal of the squared mean of the dependent variable and the number of observations (Kooijman et al. 2008a, b; Kooijman and Bedaux 1996). Furthermore, some data sets, notably growth trajectories and cumulative reproductive effort, do not meet the requirement of independence. We assume that this violation has only minor impact on the estimation of parameter values and confidence intervals.

We consider four submodels in our tests. Submodel 1 assumes that only feeding is affected; submodel 2 assumes that only maintenance is affected; submodel 3 assumes that the toxicant has the same effect on both maintenance and feeding; and submodel 4 assumes that both maintenance and feeding are impacted but by different amounts. Because these four submodels form a web of nested models, we can use the likelihood ratio test to evaluate which submodel best fits the data (Bolker 2008). Among submodels 1, 2 and 3, the submodel with the lowest negative log likelihood value fits the data best; as the deviance is approximately chi-square distributed with one degree of freedom, submodel 4 fits the data best if its associated negative log likelihood value is at least 1.92 (95% confidence level) or 3.32 (99% confidence level) less than the negative log likelihood values of other submodels. The likelihood ratio test cannot be used to analyze data representing mean values without known standard deviations, since standard deviations are required to compute likelihood values. For those data, we assumed before that the standard deviation is proportional to the calculated mean of the dependent variable. Therefore, we can still use the methodology to calculate pseudo-likelihood values in order to rank submodels in terms of fitting success, but the significance of this ranking remains unknown.

Test against experimental data

We test the toxicity submodels with published data about effects of metallic and lipophilic compounds on mussels, oysters, earthworms, zebra fish and water fleas. The results are presented in an increasing order of biological complexity. First, we examine the toxicant effects on a primary energy flows, i.e. the feeding and respiration rate. Then, we study the effects on growth, which is assumed to be a consequence of changes in the primary flows. Finally, we examine the effects on both growth and reproduction together.

Toxic effects on feeding and respiration are examined with data on similarly sized blue mussels, Mytilus edulis exposed to pentachlorophenol (PCP) (Widdows and Donkin 1991). Because both rates are clearly affected by PCP (see Fig. 2), we fit the data to submodel 3 and 4 only. Given that PCP is an uncoupler of the respiratory chain, it is noteworthy that the respiration rate is an increasing rather than decreasing function of the body burden of PCP; mitochondrial activity depends on the availability of NADH, which, in contrast to standard in vitro conditions, may be in limited supply in the cell. We neglect the no-effect body burden, since the data indicate that this value is relatively small, and exclude data obtained with mussels with extremely low filtering rates, i.e. those obtained at the highest PCP level, as those mussels switch to anaerobic metabolism (Wang and Widdows 1993), disqualifying oxygen consumption as a proper measure for total metabolism. Submodel 3 appropriately describes the effects of PCP on feeding and respiration (see Fig. 2). The difference in log likelihood values between submodel 3 and 4 is 2.79 (see Table 3), implying that the effects on maintenance and feeding differ significantly at the 0.05 but not at the 0.01 level.

Fig. 2
figure 2

Toxic effects of PCP on the rates of feeding (circles, solid line) and respiration (squares, broken line) in the mussel Mytilus edulis (data from Widdows and Donkin 1991). The model fit assumes similar toxic effects on feeding and maintenance (C K,M  = C K,F  = 73.6(4.8) nmol g−1)

Full size image
Table 3 Negative likelihood values
Full size table

Toxicants affect growth via reduced assimilation and increased maintenance rates. We test these indirect effects with data on two bivalves (Mytilus galloprovencialis and Crassostrea gigas), an earthworm (Lumbricus rubellus) and a fish (Brachydanio rerio). The growth of larvae of M. galloprovencialis and C. gigas exposed to several levels of mercury starting was followed for 10 days (Beiras and His 1995). However, larvae of bivalves typically do not exhibit a von Bertalanffy growth pattern (Bayne 1976), possibly because their capacity to catch food increases with size or time. For simplicity, we assume that larvae need a 2 day period to gain optimal food-catching capabilities and discard earlier data. Despite the fact that mussels and oysters are closely related, growth of mussel larvae is described best by submodel 1 (mercury affects feeding only) and worst by submodel 2 (mercury affects maintenance only), whereas the growth of oyster larvae is best described by submodel 2 and worst by submodel 1 (see Fig. 3a, b; Table 3). This paradox indicates that the quality of the data is insufficient to discriminate between submodels or that the submodels differ little in describing growth patterns. The latter indication is supported by the estimates of DEB parameter, which differ little among submodels.

Fig. 3
figure 3

Growth of Mytilus galloprovencialis (a) and Crassostrea gigas (b) larvae with mercury (data from Beiras and His 1995) and growth of Lumbricus rubellus with copper (data from Ma, cited in Klok and deRoos AM 1996). Ambient mercury concentrations in panel A and B are 0 (circle), 5 (plus), 10 (square), 20 (asterisk) and 40 nM (diamond); cupper densities in soil (dryweight) in panel C are 0.21 (circle), 0.95 (plus), 2.28 (square) and 5.70 (asterisk) μmol/g. Error bars represent standard deviations. The curves represent, from top to bottom, model fits at increasing mercury or copper levels. The fit in panel A assumes toxic effects on feeding (L ∞ = 350 μm; L 0 = 81.2(0.5) μm; r B  = 21.5(1.2) a−1; C K,F  = 63.1(11.4) nM); the fit in panel B assumes toxic effects on maintenance (L ∞ = 350 μm; L 0 = 66.6(0.7) μm; r B  = 19.7(1.1) a−1; nC K,M  = 27.7(4.5) nM); the fit in panel C assumes toxic effects on feeding (L ∞ = 1.4(0.1) g−1/3; L 0 = 0.2(0.0) g−1/3; r B  = 3.4(0.4) a−1; C K,F  = 10.6(0.9) μmol g−1 C NEC,F  = 0.5(0.3) μmol g−1)

Full size image

Consistent with the fact that copper is an essential nutrient, the growth of the earthworm L. rubellus (data from Ma, cited in Klok and deRoos 1996) is not affected by relatively low levels of dietary copper (see Fig. 3c; note that the shape of the curves differ from the previous example because a weight rather than a length measure is used). However, at higher levels, copper strongly reduces growth. The submodel that assumes an effect on feeding only fits the data best (see Table 3), but DEB parameter estimates and curve fits obtained with the other two submodels (results not shown) are quite similar.

The growth of the zebrafish B. rerio is decreased by benzo(k)fluoranthene and by a mixture of six polycyclic hydrocarbons (PAHs) (see Fig. 4; data from Hooftman et al. 1993). Due to a substantially increased mortality at higher benzo(k)fluoranthene concentrations, the first three data points in Fig. 4a have the most influence on the fit. The ratio of PAH concentrations in the mixture is constant, so the concentration of the mixture can be directly used in the toxicity functions, assuming that the PAHs work additively. Submodel 4 (different effects on feeding and maintenance) fits both data sets significantly better than other submodels (see Table 3). The second best fit is obtained assuming effects on maintenance only with the single PAH and effects on feeding only with the mixture of PAHs, which again indicates that the data cannot discriminate among submodels. This is also suggested by the DEB the estimates of DEB parameters and curve fits (results not shown), which are comparable among the submodels.

Fig. 4
figure 4

The length of Brachidanio rerio grown in an intermittent flow-through system for 37 days after fertilization with benzo(k)fluoranthene (a) or a mixture of 6 polycyclic aromatic hydrocarbons at a highest concentration of 1 nM each (b). Error bars represent standard deviations. Curves represent model fits assuming toxic effects on maintenance (parameter values in panel A are L ∞ = 4.89(0.04) cm; L 0 = 0.40 cm; r B  = 0.15 a−1; C K,M  = 0.92(0.05) nM; C NEC,M  = 0.92(0.05) nM; and in panel B L ∞ = 3.62(0.01) cm; L 0 = 0.40 cm; r B  = 0.15 a−1; C K,M  = 0.44(0.00) nM); C NEC,M  = 0.10 nM) (data from Hooftman et al. 1993; Hooftman and Ruiter 1992)

Full size image

The water flea Daphnia magna produces less and needs a longer time to produce the first clutch of eggs when exposed to tetradifon (see Fig. 5a, b, c; data from Villarroel et al. 2000) In order to fit these data, we need to estimate the length increment from the time the water fleas reach reproductive maturity until the time they release the first clutch. This length increment varies a bit among treatments, but in the fits we assume it has a constant value in order to avoid undesirable mathematical complexity as well as biological complexity due to unknown degrees of synchrony of molts among individuals. The submodels assuming effects on feeding only and similar effects on feeding and maintenance are nearly indistinguishable (see Table 3). The submodel assuming different effects on assimilation and maintenance fits the data marginally better, but this improvement is not significant.

Fig. 5
figure 5

Production in Daphnia magna grown with tetradifon (a, b and c, data from Villarroel et al. 2000) and pyridine with ambient concentrations 0 (circle), 25 (plus), 50 (square) and 100 mg l−1 (asterisk) (d and, e data from Santojanni et al. 1998). The curves of length, reproduction and time to first brood with tetradifon represent model fits assuming toxic effects on feeding; error bars represent standard deviations (L ’ = 0.87(0.06) mm; L 0 = 0.90 mm; L p  = 2.10 mm; L ∞ = 5.30 mm; r B  = 30.8(1.2) a−1; g = 1; f = 1; p r  = 0.78(0.07) eggs mm−3; C K,F  = 3.40(0.28) mg l−1). The curves of length and cumulative reproduction with pyridine represent model fits assuming toxic effects on maintenance (L 0 = 0.90 mm; L p  = 2.10 mm; L ∞ = 5.30 mm; r B  = 25.8(0.8) a−1; g = 1; f = 0.7; p r  = 0.42(0.04) eggs mm−3; C K,M  = 141.5(03.9) mg l−1)

Full size image

A similar trend is observed when D. magna is exposed to pyridine (see Fig. 5d, e; data from Santojanni et al. 1998). The submodel assuming effects on maintenance only fits the data best (see Table 3), but the submodel assuming similar effects on maintenance and assimilation also describes the data well (fits not shown).

Population consequences

The impact of toxicants on the long-term population growth rate can be calculated from data about reproduction and survival (see Eq. 9 in Table 1). This was one of the original motivations for the development of DEB theory (Kooijman and Metz 1984). Thus, DEB models can be applied to analyze data from life cycle test in order to determine effects of pollutants on the long term population growth rate (Billoir et al. 2007; Jager et al. 2004; Klok et al. 2007; Lopes et al. 2005). We have analyzed reproduction in Daphnia in the presence of two toxicants, tetradifon and pyridine, but the survival characteristics are unknown. In order to visualize the impact of the different toxicity submodels we have analyzed in this paper, we make a simplistic assumption about mortality in those experiments: the background mortality is independent of age and toxicant action. With this assumption, the Lotka–Euler equation (see Eq. 9), and the estimated parameter values listed (see Fig. 6), we calculate the long term population growth rate as a function of the tetradifon or pyridine concentration and as a function of the food density and scale the results to the long term population growth rate calculated without pollutants: r(C,f)/r(0,f).

Fig. 6
figure 6

The long term population growth rate of Daphnia magna declines as a function of the food density and the tetradifon (a) or pyridine (c) concentration; the growth rate is scaled to the growth rate in absence of pollutants: r(C,f)/r(0,f). Population growth rates are estimated with the parameters of the submodel assuming similar effects on feeding and maintenance (see Table 3; negligible background mortality) and scaled to the estimated population growth rates in absence of toxicants. The tetradifon (b) and pyridine (d) concentration at which the population growth rate is reduced by 50% (EC50) increases with food density in the submodel assuming similar effects on feeding and maintenance (solid line), effects on feeding only (broken line) and assuming effects on maintenance only (dotted line)

Full size image

The scaled population growth rate, assuming similar effects on assimilation and maintenance, declines as a function of the toxicant concentration and the food density (see Fig. 6a, c). The scaled population growth rate declines fairly gradually until close to a critical combination of toxicant concentration and food density at which point it drops rapidly to zero, a property anticipated in the seminal paper by Kooijman and Metz (1984). Below this critical combination, the animals do not reach the threshold size needed for reproduction. The results obtained with the other toxicity submodels yield numerically comparable results (results not shown); the main difference is that these submodels predict different critical combinations of toxicant concentrations and food densities. As a result, the difference among submodels is most pronounced at low population growth rates. Figure 6b and d illustrate how the EC50 for the population growth rate increases with increasing food densities. The predictions differ somewhat among submodels, but given the uncertainty in parameter values, those differences may not be biologically significant.

Discussion

The central idea in modeling toxicant effects in a DEB framework is that contaminants impact the fluxes of energy between the organism and its environment and/or within the organism by changing parameter values. In this paper, we focused on toxic impact on the parameters quantifying the rates of feeding and maintenance. We have considered four toxic effect mechanisms: feeding only is affected; maintenance only is affected; assimilation and maintenance are both affected in a similar way (toxicity parameters have the same values in both effect modules); and assimilation and maintenance are both affected but with different values for the toxicity parameters. We specifically excluded toxic effects related to survival, offspring viability and resource partitioning from this paper; these effects are important but are outside the scope of this work. Background mortality and diminished off-spring viability due to toxicant action have been modeled by Kooijman and Bedaux (Kooijman and Bedaux 1996), but other, potentially more relevant mortality mechanisms, such as an increased susceptibility to infectious diseases, parasites or predators, have yet to be investigated in a DEB modeling framework, as does the change in partitioning of resources due to endocrine disruptors.

Our most striking conclusion is that the data analyzed in this study are generally well described by any of the toxic effect submodels. Our analysis of data on growth and reproduction in the presence of pollutants indicate that these type of data generally do not have the discriminatory power needed to distinguish between submodels, since the parameters of each submodel can be tweaked to yield satisfactory descriptions of production data. Only direct measurement, e.g. the data on respiration and feeding by mussels in the presence of pentachlorophenol (see Fig. 2), shows convincingly that any compound affects multiple processes. Because in that case the feeding rate declines concomitantly with an increase in the respiration rate as the toxicant concentration increases, it is obvious that both feeding and maintenance are affected. A second striking finding was that similar toxicant-organism combinations were not necessarily best described by the same toxicity submodel. The growth of mussel and oyster larvae with mercury are best described by different submodels: effects on assimilation only and effects on maintenance only, respectively (see Table 3). Also the fish growth data with either a single PAH or a mixture of PAHs yields different best fitting submodels. Thus, in order to distinguish between submodels, we need to include other types of measurements, such as feeding and respiration rates. In the absence of such information we propose as a default the submodel in which effects on maintenance and assimilation are both included, but are described by a single parameter. This submodel never gave a poor fit to data, and introduces the same number of parameters as the single-effect submodels. Furthermore, the parameters in this submodel can be tweaked such that projections of the population growth rate are comparable to those of other submodels (see Fig. 6), which suggests that identification of the ‘correct’ submodel is relatively unimportant for describing population level effects under constant food conditions.

Our approach of modeling toxic effects differs from a previous one based on DEB theory (Kooijman et al. 2008a, b; Kooijman and Bedaux 1996; Kooijman 2000). The latter assumes that every parameter is subject to toxicant action beyond a threshold level in structural mass and that the threshold of the most sensitive parameter is substantially lower than other thresholds. This implies that in practice only a single parameter would be affected by toxicants. As an approximation of a potentially complex relationship, the toxic effect on this parameter is assumed to change linearly with the toxicant density in structural mass. Because this linear relationship implies that the rates of feeding and assimilation can become negative, we propose to model toxicant effects on parameter values that decrease with increasing toxicant levels with a hyperbolic function, such as in Eq. 6 (see Table 1). A hyperbolic relationship is consistent with the synthesizing unit concept, a more recent development in DEB theory, and can be derived from first principles (Zonneveld 1998). Furthermore, at least in the original publications (Kooijman and Bedaux 1996), sublethal toxic effects are assumed to affect compound rather than primary DEB parameters (as in this paper), which obscures the relationship between effect and process. For instance, in the growth cost submodel of Kooijman and Bedaux (Kooijman and Bedaux 1996), either toxicants increase the costs for both structural biomass formation and maintenance with identical linear toxic effect functions (otherwise their maintenance cost submodel would be in effect as well), or toxicants decrease the maximum reserve density and the maximum rates of assimilation and feeding with identical hyperbolic (cf. Eq. 6 in Table 1) toxic effect functions (otherwise their assimilation submodel would be in effect as well; note that the latter implies a nonlinear toxic effect function). In addition to an unclear relationship between effect and process, using compound parameters as the target of toxicants is inconsistent with the assumption that no-effect threshold levels differ substantially among DEB parameters, as multiple primary parameters are involved.

Toxicants often elicit a wide variety of sub-organismal responses. For instance, in water fleas cadmium induces changes in the expression of genes involved in immune and stress responses, digestion, metabolism and signal transduction (Soetaert et al. 2007). It would be unlikely that such a wide range of responses could be captured in an effect on a single DEB parameter, with the possible exception of the primary maintenance rate parameter, since it can be argued that toxicants generally threaten the physiological integrity of organisms. Moreover, even if a toxicant targets a single process, it is not necessarily true that only one DEB parameter will be affected, since a given process may be involved in more than one energy flux, with the added complication that the relative level of involvement is a function of resource availability (Hanegraaf and Muller 2001; Jager and Kooijman 2005).

In conclusion, it is not feasible to aspire to a single mechanistic model that describes the sublethal effects of every toxicant on every species. We think that mechanistic models can only be formulated when the molecular mode of action of a toxicant is highly specific. With notable exceptions (Jager and Kooijman 2005), this seems impossible given the virtually endless list of toxic pollutants and the complex physiological interactions that ensue after exposure. However, it is feasible to incorporate toxic effects in a DEB modeling framework, because molecular toxic effects can translate in changes in only a limited number of parameter values. The resulting toxicity parameters are less ambiguous, more accurate and ecologically more relevant than classic measures, such as the NOEC and EC x , which cannot be applied to estimate population level consequences of toxicants. Furthermore, the impact of multiple stressors can be delineated in a DEB modeling framework. For example, factors proposed as possible causes of population decline in the Northern Right whale include contaminants (PCBs) and food shortage due to the North Atlantic Oscillation; both factors can be incorporated in a DEB model (Klanjscek et al. 2007). We therefore recommend including the estimation of DEB toxicity parameters in protocols for the testing of toxic substances.