Introduction

In their environments, organisms can be exposed to a wide range of naturally occurring or synthetic toxic substances (Moriarty 1999). Physiological or population effects of such chemicals are mostly examined in dose–response studies, where the dosage of a toxicant is varied on a certain biomass of organisms (Moriarty 1999). However, an increasing number of studies demonstrate that toxicity effects may often not only be dependent on the dose, but also on the biomass of exposed organisms. In such cases, toxicity may be alleviated by a high organism biomass because the concentration of a toxicant is reduced to sub-lethal levels due to joint uptake or active detoxification. This effect has been described for a wide range of biological systems such as heavy metal accumulation in organisms (Duxbury and McIntyre 1989; Pickhardt and others 2002), phytotoxins in microbes and plants (for example, Weidenhamer and others 1989; Greig and Travisano 2008; Pollock and others 2008), and for drug treatments of infectious bacteria or cancer (for example, Brook 1989; Kobayashi and others 1992; Brandt and others 2004). In ecotoxicological studies, this mechanism is often referred to as “density-dependent toxicity” (for example, Duxbury and McIntyre 1989; Kobayashi and others 1992; van der Heide and others 2008) or the “dilution effect” (for example, Karimi and others 2007; Pollock and others 2008), whereas it is described as the “inoculum effect” in many pharmaceutical studies (for example, Brook 1989; Kobayashi and others 1992; Brandt and others 2004).

This density-dependent toxicity implies a positive feedback between population density and toxic substance concentration, because an increase in biomass alleviates toxicity, which can in turn further amplify biomass growth. Theory suggests that if such a positive feedback mechanism is strong enough, it could lead to alternative stable states (also called bistability) and hysteresis (Carpenter 2001; Scheffer and others 2001; Scheffer and Carpenter 2003). This implies that environmental changes or disturbances (for example, disease) may push the population beyond a critical threshold, causing a collapse (for example, mass mortality) to an alternative stable state (Scheffer and others 2001). Implications of non-linear response and threshold behavior in biological systems can be profound. Shifts in populations with positive feedbacks are typically hard to predict and recovery of a collapsed system is often difficult.

In this study, we examine whether a positive feedback caused by density-dependent toxicity may cause alternative stable states. First, we show the basic idea in a simple conceptual model, describing a generalized positive feedback system between a population of organisms and a toxic compound. We analyzed two different assumptions of how an established population can alleviate toxicity: the organisms can alleviate toxicity actively (“joint detoxification”) or can take up and store a limited amount of the toxic substance per organism or unit of biomass (“growth dilution”). Next, we created a more realistic model, describing ammonia toxicity in seagrass ecosystems, to further explore our theory for an empirical situation. We analyze whether this model can have alternative stable states using realistic parameter settings based on laboratory experiments and literature.

Methods

Conceptual Models

The basic mechanism of this feedback can be shown in a very simple model. This model describes a system with a population of organisms that cannot exceed a carrying capacity, for instance due to limited nutrients, food or space. Furthermore, we model a toxicant that increases mortality of the population:

$$ {\frac{{{\text{d}}X}}{{{\text{d}}t}}} = a\,\left( {1 - {\frac{X}{{K_{X} }}}} \right)\,X - b\,X\,T $$
(1)

X describes the biomass of the population, a is the maximum net growth rate per unit of time and K X is the carrying capacity. Parameter b is a mortality constant which, multiplied with the concentration of toxicant T represents the mortality rate in population X per unit of time. Although the linear effect of T in the model may often be an oversimplification of reality, simulations comparing linear, Monod and Hill functions, showed that basic model behavior is not sensitive to the choice of the function here in the sense that it consistently produced alternative stable states in a wide parameter range. Moreover, in many cases the response of organisms to toxicants can well be described by this simple relation (Hendriks and others 2005). The second equation describes the change in toxic compounds in the system per unit of time. We assume a constant external input or internal release of the toxicants similar to a chemostat model (for example, Edelstein-Keshet 1988). However, organisms can also reduce the concentration of toxicants in the system:

$$ {\frac{{{\text{d}}T}}{{{\text{d}}t}}} = (T_{\text{in}} - T)\,p - T\,f(X) $$
(2)

where T in is the maximum equilibrium concentration of toxicants that can accumulate in the system and p is the dilution rate (that is, the fraction of the volume replaced per time unit). We assume that the organisms reduce the concentration of toxicants as a function of their biomass. This density-dependent alleviation of toxicity is essential for the feedback.

We analyzed two different assumptions. In the first case, organisms can actively transform the toxic compound into a harmless (or even useful) substance, for instance metabolically or by excretion of compounds that chemically react with the toxicant. This mechanism can result in decreased mortality if the density of organisms if sufficient to reduce ambient toxicant concentrations (“joint detoxification”). Here, f(X) is described by d 1 X, where d 1 is a constant describing the uptake rate of T per unit of X. For the second mechanism, we assume that the toxicant is not converted into a harmless substance, but is stored is the tissues of the organisms (for example, heavy metal storage in fat tissue). In this case, organisms can only take up a limited amount of the toxic substance per unit of biomass and detoxification is thus dependent on growth (“growth dilution”). In this case, f(X) is replaced with the term d 2 a X (1 − X/K x ), where d 2 describes the uptake of T proportional to growth of X. Notably, net uptake of the toxicant in the growth dilution model becomes zero when carrying capacity is reached. In reality, natural mortality and regrowth around carrying capacity in the population would result in some uptake and release dynamics of the toxicant. For simplicity, however, we chose not to include a mortality term in the models, because the general behavior of the model would remain unchanged.

In systems with a relatively high dilution rate p, dynamics of T are much faster than those of the organisms. Therefore, we can assume a quasi steady state (that is, \( {\frac{{{\text{d}}T}}{{{\text{d}}t}}} = 0 \)) without any consequences for the equilibrium density of organisms and the behavior of the model. This assumption simplifies the model to:

$$ {\frac{{{\text{d}}X}}{{{\text{d}}t}}} = a\,\left( {1 - {\frac{X}{{K_{X} }}}} \right)\,X - b\,X\,{\frac{{p T_{\text{in}} }}{p + f(X)}} $$
(3)

The conditions for alternative stable states of this simple model can be determined analytically for the joint detoxification assumption or numerically in case of the growth dilution model (online Appendix 1).

Specific Model of NH x Toxicity in Eelgrass

Recent studies have demonstrated that positive feedbacks are important mechanisms in seagrass ecosystems (van der Heide and others 2007, 2008, 2009). This model, based on empirical data, describes a feedback mechanism between the temperate seagrass Zostera marina (commonly called eelgrass), reduced nitrogen (NH x ) and potentially lethal gaseous ammonia (NH3) in the water layer. We chose this system as a more realistic analysis for our theory because recent research demonstrated that susceptibility of eelgrass to NHx toxicity is highly dependent on vegetation density, indicating that positive feedbacks between eelgrass and reduced nitrogen may lead to alternative stable states in sheltered estuaries with high exposure to NH x . In these systems, high concentrations of NH x may be caused by for instance discharges of waste or river water and degradation of phytoplankton or macro-algal mats (van der Heide and others 2008). Also, because ammonium uptake is well studied in eelgrass, model parameters could be reliably estimated based on these studies and results from our own experimental work.

We based the model on the joint detoxification assumption. In the first place because ammonium is used as a nutrient by the plants and it is therefore metabolized. Secondly, an eelgrass shoot may discard excess nutrients by replacing its leaves without resorbing the nutrients stored in the leaves that are lost (Hemminga and others 1999). Moreover, compared to other vascular plants the lifespan of eelgrass leaves is relatively short, suggesting that excess nitrogen stored in the leaves may be exported from the system through rejection of leaves (Hemminga and others 1999).

In the model, survival of eelgrass is dependent on the concentration of NH3 in the water layer. The equation describing the change in eelgrass shoot density per day (dZ/dt) is similar to Eq. 1:

$$ {\frac{{{\text{d}}Z}}{{{\text{d}}t}}} = r\,\left( {1 - {\frac{Z}{{K_{Z} }}}} \right)\,Z - m\,f\left( {{\text{NH}}_{3} } \right)\,{\text{Z}} $$
(4)

With r as the maximum net growth rate (day−1), K Z as the carrying capacity (shoots m−2), and m as the maximum mortality rate (day−1). The toxic effect of NH3 is described by the function f(NH3). To estimate the toxicity effect of NH3 in eelgrass, we recalculated and analyzed experimental data of Van der Heide and others (2008) (online Appendix 2). Our analyses revealed that toxicity by NH3 in eelgrasses is best described by a Hill-curve (Figure 1). This is an equation that is typically used to describe toxicity in organisms. The function expresses a sigmoid toxicity effect in the organism in response to increasing exposure to a toxicant (Hill 1910):

$$ M = i + M_{\max } \,{\frac{{{\text{NH}}_{3}^{n} }}{{{\text{NH}}_{3}^{n} + {\text{H}}_{{{\text{NH}}_{3} }}^{n} }}} $$
(5)

Here M describes the fraction of leaf tissue mortality in eelgrass after 5 days of exposure to NH3, i is the background leaf mortality at zero exposure. \( {\text{H}}_{{{\text{NH}}_{3} }} \) is the half-saturation constant (mmol m−3) and n is a dimensionless exponent determining the slope of the curve. To describe the effect of NH3 in our model, we adopted the part of Eq. 5 describing the relative effect of NH3 on eelgrass mortality:

$$ f\left( {{\text{NH}}_{3} } \right) = {\frac{{{\text{NH}}_{3}^{n} }}{{{\text{NH}}_{3}^{n} + {\text{H}}_{{{\text{NH}}_{3} }}^{n} }}} $$
(6)

In water, the total NH x concentration is made of the sum of NH3 and ammonium (NH4 +). NH3 and NH4 + are in equilibrium and the balance between these compounds is determined by the pH of the water. The concentration of NH3 in the water can be calculated from the pH and the total concentration of reduced nitrogen in the water layer:

$$ {\text{NH}}_{3} = {\frac{{k_{a} \,{\text{NH}}_{x} }}{{k_{a} + 10^{{ - {\text{pH}}}} }}} $$
(7)

where k a is the dimensionless dissociation constant of NH x in water with a salinity of 16 PSU at 20°C. The change of NH x in the water layer is described by the second differential equation:

$$ {\frac{{{\text{dNH}}_{x} }}{{{\text{d}}t}}} = \left( {{\text{NH}}_{x} {\text{in}} - {\text{NH}}_{x} } \right)\,R - U_{\max } \,{\frac{{{\text{NH}}_{x} }}{{{\text{NH}}_{x} + {\text{H}}_{{{\text{NH}}_{x} }} }}}\,f\left( Z \right) $$
(8)

With NH x in as the NH x concentration of the water flowing into the meadow and R as the dilution rate of the water inside the meadow. U max is the maximum uptake rate of NH x by eelgrass (mmol g dry weight−1 day−1) and \( {\text{H}}_{{{\text{NH}}_{x} }} \) is the half-saturation constant for NH x uptake (mmol m−3). Finally, f(Z) is a function describing the conversion from eelgrass shoot density (shoot m−2) to the amount of dry weight biomass per unit of volume:

$$ f\left( Z \right) = {\frac{Z}{C}}\,{\text{Dw}}_{Z} $$
(9)

Here C is the height of the canopy (m) and Dw Z is the dry weight of one eelgrass shoot (g).

Figure 1
figure 1

Response of eelgrass shoots to ammonia at varying concentrations after 5 days of exposure.

Bifurcation Analysis

We analyzed the stability of the equilibria of the model at varying settings of key parameters. Critical thresholds were determined by a numerical procedure. The key parameter was increased in small steps, after which the model was run to stabilize to its equilibrium. Next, this analysis was also performed backwards, by decreasing the key parameter in small steps. These analyses were combined to construct bifurcation plots of various parameters. We determined unstable equilibria by making the quasi steady state assumption (\( {\frac{{{\text{d}}T}}{{{\text{d}}t}}} = 0 \)) and plotting equilibria for different values of the control parameters in GRIND for MATLAB.

Results

Conceptual Models

Figure 2A and B shows phase planes of the joint detoxification and growth dilution model, respectively, based on the default settings presented in Table 1. Both the graphs show two stable equilibria and one unstable equilibrium (saddle point). Whereas toxicant concentrations show a straightforward decrease with increasing biomass in the joint detoxification model, toxicant levels in the growth dilution model increase again when X nears its carrying capacity. This is because the population growth and therefore also the detoxification rate is highest halfway to the population’s carrying capacity. Next, we analyzed the sensitivity of both models to varying values of the maximum toxicant concentration (T in) in a one-dimensional bifurcation plot (Figure 2C, D). The results demonstrate that both models can have alternative stable states, one without organisms and one with a population that can alleviate toxicity. The systems collapse to a bare state when organism density is pushed below the critical threshold (Figure 2C, D, dashed lines).

Figure 2
figure 2

Analyses of the conceptual models. A and B Nullclines at default settings of the “joint detoxification” model and “growth dilution” model, respectively. The closed dots represent stable equilibria in the models, the open dots are unstable saddle points. The separatrix indicates the critical boundary above which the population can survive and develop to the equilibrium. C indicates a colonized state, B is a bare state. C and D Bifurcation analyses of the “joint detoxification” model and “growth dilution” model, respectively, with varying values of T in. Solid lines represent stable equilibria, whereas the dashed line indicates unstable equilibria. Dots indicate bifurcation points (F fold bifurcation, T transcritical bifurcation), arrows show the direction of change. Note that all parameter settings for the “joint detoxification” model and the “growth dilution” model were identical, except for parameter b, which were set at 0.2 and 1.5, respectively.

Table 1 Variables and Default Parameter Settings of the Conceptual Model

A more thorough bifurcation analysis of the joint detoxification model shows that the conditions for alternative stable states in this model are relatively simple. After reducing the number of model parameters to 2 by non-dimensionalization (online Appendix 1), conditions for alternative stable state can be summarized in a simple 2D plot. This plot shows all parameter combinations at which alternative stable states occur (Figure 3A). It appears that there are two prerequisites that determine whether the feedback is strong enough to cause alternative stable states. First, the equilibrium concentration of toxicant without organisms (T in) should be able to prevent colonization of the organism [that is, its effect on the organisms (b T in) should be higher than the maximum growth rate of the population (a)]. The second prerequisite is that the effect of a full-grown population (K X ) on the toxic substances should be strong enough to let the concentration of the toxicant decrease, that is, the refresh rate of the toxic substance (p) should be less than the maximum effect of the organisms (K X d). Note that this means that the chances for alternative equilibria increase if the turnover rate of the toxic substance (p) is low. If these two prerequisites are met, there is a range of T in with alternative stable states. With increasing carrying capacity (or decreasing turnover rate), this range increases.

Figure 3
figure 3

Two-dimensional plots of the scaled conceptual models (see online Appendix 1). On the axes are the two combined parameters: the scaled toxic load \( \alpha = {\frac{b}{a}}T_{\text{in}} \) and the scaled carrying capacity of X, which is \( \beta =\frac{d}{p}K_{X} \) for the “joint detoxification” model and \( \gamma = \frac{d}{p}a\,K \) for the “growth dilution” model. The figures give all parameter combinations where we get alternative stable states. (C colonized, C/B alternative states, B bare only, F (solid line) fold bifurcation, T (dashed line) transcritical bifurcation).

The growth dilution model is too complex for a similar analytical bifurcation analysis. However, we did this analysis numerically, showing very similar results (Figure 3B). Although the range for bistability in this model is narrower when compared to the joint detoxification model, there is still a rather large parameter space with alternative stable states. Moreover, the qualitative effect of scaled toxicity and carrying capacity is remarkably similar.

Specific Model of NH x Toxicity in Eelgrass

The eelgrass model was parameterized to describe a sheltered estuary, where water mixing between the eelgrass meadow and its surroundings is limited (Table 2). Water flowing into the seagrass bed has a NH x concentration of 100 μmol l−1 NH x , a value comparable to various measurements in the field (for example, Hauxwell and others 2001; Brun and others 2002). In these systems pH can vary strongly. At night, pH is generally around 8 whereas pH can rise up to 9 or even 10 during the day, due to photosynthesis of algae and seagrass itself (Choo and others 2002; Feike and others 2007; van der Heide and others 2008), hence leading to higher NH3 concentrations. For simplicity, we assumed an average pH of 8.5 for our model system.

Table 2 Variables and Default Parameter Settings of the Eelgrass Model

The nullclines of this model at default parameter settings are presented in Figure 4A. Similar to the conceptual model, the graph shows one unstable equilibrium and two stable points. Depending on the initial conditions, the meadow will either develop towards carrying capacity or collapse to a bare state. A bifurcation analysis on the NH x concentration of the water flowing into the eelgrass meadow (NH x in) reveals that alternative stable states are present over a wide range of realistic concentrations, from 75 to over 158 μmol l−1 (Figure 4B). We analyzed the interactive effects of NH x in, pH, and dilution rate R, because these parameters are often variable in the field. Results demonstrate that the effect of the NH x concentration of the incoming water is highly dependent on both pH and the dilution rate of the water inside the meadow (Figure 5A, B). The analysis shows that alternative stable states are present at pH values higher than 7.9 (Figure 5A). Below pH 7.9, the toxicity of NH x is too low as only little NH x is present as toxic NH3. Therefore, the meadow tolerates extremely high concentrations of NH x in the incoming water. Sensitivity to NH x exposure increases strongly with rising pH levels, as the NH4 +/NH3 equilibrium shifts towards NH3. At pH 10, alternative stable states exist between NH x concentrations of 10 and 55 μmol l−1 in the water flowing into the meadow. Figure 5B demonstrates the interactive effects of NH x in and water dilution rate R. No alternative stable states are present when the concentration of NH x is below 75 μmol l−1, because these concentrations are not lethal for the eelgrass plants at pH 8.5 (compare the first prerequisite of the conceptual model). The effect of NH x becomes dependent on both NH x input concentrations and the turnover rate R, when NH x concentrations of the incoming water rise above the 75 μmol l−1 threshold. Alternative stable states exist far beyond NH x concentrations of 500 μmol l−1 for NH x in when R drops below 1 day−1.

Figure 4
figure 4

Analyses of the eelgrass model. A Nullclines of the model at default settings. B Bifurcation analysis of the model with varying NH x concentrations in the incoming water (NH x in). See Figure 2 for the meaning of symbols used.

Figure 5
figure 5

Two-dimensional bifurcation analyses of the eelgrass model. A Bifurcation analysis with varying pH and NH x concentrations in the incoming water (NH x in). B Bifurcation plot with varying dilution rates of the water in the meadow (R) and NH x concentrations in the incoming water (NH x in). Solid lines represent fold bifurcations, whereas the dashed lines indicate transcritical bifurcations. B indicates a bare state; C/B indicates the area where alternative stable states occur. Left of the dashed lines (indicated with C), eelgrass presence is the only stable state.

Discussion

We show in both a conceptual and a more realistic model that “density-dependent toxicity,” a positive feedback mechanism between a population of organisms and a toxic compound may lead to bistability in biological systems. Organisms may alleviate adverse effects of the toxicant by actively lowering ambient concentrations through either “joint detoxification” or “growth dilution.” Joint detoxification is a mechanism where the toxicant is actively broken down by the exposed organisms. The population can maintain itself, provided that its biomass is sufficient to reduce toxicant concentrations to a level where organism growth may equalize or exceed mortality. Growth dilution is a mechanism where the toxicant is not broken down, but is stored in the organism’s tissues. Because these tissues are only able to store a limited amount of toxicants, they will become saturated. In this case, reduction of the toxicant is dependent on population growth rather than the biomass present in the system.

Our eelgrass model suggests that density-dependent toxicity may indeed be important in real ecosystems. Although the model is somewhat more complicated, its essence is identical to our conceptual joint detoxification model. Sudden die-off events caused by high reduced nitrogen (NH x ) loads, combined with a high pH may be prevented by joint uptake if shoot density of the meadow is high enough. This mechanism fails if shoot densities are pushed below a certain threshold, resulting in a shift to a bare state. This illustrates that density-dependent toxicity can have important implications for toxicity research and management in ecosystems. For toxicity research in laboratory and field studies, our results indicate that it may be very important to choose realistic population densities instead of working with standardized biomass or densities. Moreover, our simulations also show that it is not sufficient to measure ambient toxicant levels to assess ecosystem health. At high population densities, toxicant levels may be low due to detoxification mechanism, whereas the toxicant load may actually be very high. Therefore, ecosystem monitoring should focus on determining the toxicant load in such cases. Finally, it should also be noted that when such an ecosystem collapses, it may not only affect the community structure directly. After the collapse many associated species may now also experience toxicity effects because toxicant levels will increase dramatically.

Although we studied only one example, density-dependent toxicity is most likely an important mechanism in a wide range of biological systems. Joint detoxification has also been reported in for example isoetid macrophytes. In these vegetations, ammonium toxicity can be prevented because ammonium concentrations in the pore water are actively lowered, not only by uptake, but also by density-dependent oxidation of ammonium to nitrate due to high radial oxygen loss of the roots (Smolders and others 2002). Toxic effects of sulfide in salt-marshes (Webb and others 1995; Webb and Mendelssohn 1996), seagrasses (for example, Goodman and others 1995; Pedersen and others 2004) or sulfate-rich freshwater wetlands (Lamers and others 1998; Armstrong and Armstrong 2001; van der Welle and others 2006) may be prevented in a similar way. In these systems, sulphide can be oxidized to harmless sulfate if oxygen loss by the root system is sufficiently high.

A second possible mechanism for density-dependent toxicity, growth dilution, may for instance reduce toxic effects of heavy metals; toxicants that cannot be broken down. The dilution effect increases tolerance of microbes to heavy metal exposure (Duxbury and McIntyre 1989). In aquatic ecosystems, accumulation of toxic metals in the trophic chain of food webs has been shown to be reduced with increasing concentrations of phytoplankton (Pickhardt and others 2002) or even with increasing nutritional quality of the algae (stoichiometric dilution) (Karimi and others 2007).

Although our analyses suggest that the mechanism presented in this study may lead to alternative stable states in many biological systems, it should be noted that dynamics in our models are described in a simplified manner. This implies that the models may disregard or oversimplify processes that might in reality be important. These can include factors that weaken the positive feedback as well as processes that enhance it. In general, processes strengthening the feedback may include symbiosis or natural selection leading to more resistant individuals (Brook 1989), whereas factors weakening it can include limitation of resources (for example, nutrients, water) (Weidenhamer 1996), competition with other species (Weidenhamer 1996), or disease (van der Heide and others 2007). More specifically, a factor that could weaken the positive feedback in our eelgrass model is that high shoot densities may imply a higher photosynthetic activity, leading to a higher pH in conditions with low flow rates and poor mixing. High pH in turn may lead to increased ammonia toxicity (van der Heide and others 2008). Additionally, low flow rates may result in temporary spikes in ambient NH x concentrations at the end of the growing season due to natural die-off of eelgrass itself. In contrast to systems with higher flow rates, export of seagrass litter in sheltered embayments may be slow. This may cause decomposing litter to temporarily accumulate in the system, resulting in increased release of NH x . Still, despite the fact that the described mechanism might not lead to hysteresis in all biological systems, either due to interfering factors or simply because the feedback mechanism is too weak, density-dependent toxicity may still be important. In such systems, feedbacks may still cause a strong nonlinear response of organism density to changing toxicant loads (for example, a sharp sigmoid response), leading to unexpectedly sudden community shifts and ecosystem management problems (Scheffer and others 2001).

In summary, we present a feedback mechanism between organisms and toxic compounds that may potentially lead to bistability in biological systems. Adverse effects of toxicants that are being produced in or come into the system may be prevented by actively lowering ambient concentrations through either rapid joint detoxification or growth dilution. The presented general mechanism may be important in a wide range of biological systems.