MODEL ANALYSIS OF A SIMPLE AQUATIC ECOSYSTEMS WITH SUBLETHAL TOXIC EFFECTS

. The dynamic behaviour of simple aquatic ecosystems with nutrient recycling in a chemostat, stressed by limited food availability and a toxicant, is analysed. The aim is to ﬁnd eﬀects of toxicants on the structure and functioning of the ecosystem. The starting point is an unstressed ecosystem model for nutrients, populations, detritus and their intra-and interspeciﬁc interactions, as well as the interaction with the physical environment. The fate of the toxicant includes transport and exchange between the water and the populations via two routes, directly from water via diﬀusion over the outer membrane of the organism and via consumption of contaminated food. These processes are modelled using mass-balance formulations and diﬀusion equations. At the population level the toxicant aﬀects diﬀerent biotic processes such as assimilation, growth, maintenance, reproduction, and survival, thereby changing their biological functioning. This is modelled by taking the parameters that described these processes to be dependent on the internal toxicant concentration. As a consequence, the structure of the ecosystem, that is its species composition, persistence, extinction or invasion of species and dynamics behaviour, steady state oscillatory and chaotic, can change. To analyse the long-term dynamics we use the bifurcation analysis approach. In ecotoxicological studies the concentration of the toxicant in the environment can be taken as the bifurcation parameter. The value of the concentration at a bifurcation point marks a structural change of the ecosystem. This indicates that chemical stressors are analysed mathematically in the same way as environmental (e

recycling in a chemostat, stressed by limited food availability and a toxicant, is analysed.The aim is to find effects of toxicants on the structure and functioning of the ecosystem.The starting point is an unstressed ecosystem model for nutrients, populations, detritus and their intra-and interspecific interactions, as well as the interaction with the physical environment.The fate of the toxicant includes transport and exchange between the water and the populations via two routes, directly from water via diffusion over the outer membrane of the organism and via consumption of contaminated food.These processes are modelled using mass-balance formulations and diffusion equations.At the population level the toxicant affects different biotic processes such as assimilation, growth, maintenance, reproduction, and survival, thereby changing their biological functioning.This is modelled by taking the parameters that described these processes to be dependent on the internal toxicant concentration.As a consequence, the structure of the ecosystem, that is its species composition, persistence, extinction or invasion of species and dynamics behaviour, steady state oscillatory and chaotic, can change.To analyse the long-term dynamics we use the bifurcation analysis approach.In ecotoxicological studies the concentration of the toxicant in the environment can be taken as the bifurcation parameter.The value of the concentration at a bifurcation point marks a structural change of the ecosystem.This indicates that chemical stressors are analysed mathematically in the same way as environmental (e.g.temperature) and ecological (e.g.predation) stressors.Hence, this allows an integrated approach where different type of stressors are analysed simultaneously.Environmental regimes and toxic stress levels at which no toxic effects occur and where the ecosystem is resistant will be derived.A numerical continuation technique to calculate the boundaries of these regions will be given.Introduction.The individual level is the natural starting point for studying toxicological effects of biological systems, for the individuals suffer from toxic stress.However research is not primarily concerned with the fate of a specific individual, but with that of populations or higher levels of organization, such as communities or ecosystems, in which the physical environment is also important.In [13,15,14,12,17] the step from the individual to the population level is made for species with multiple life-stages, such as Daphnia and fish, which suffer sublethal toxic stress in a constant environment.An individual-based, physiologically structured population model was used, mathematically leading to a partial differential equation.
In this paper we study the sublethal effects of toxicants on the functioning and structure of aquatic ecosystems.The functioning of ecosystems involve the transformation by the biota of material between organics and inorganic pools through decomposition, nutrient mineralization, assimilation, and production.Structural properties of an ecosystem are the species composition and its dynamical state, such as equilibrium or oscillatory behaviour.Besides ecological and environmental stresses there is stress due to the toxicant.Since the life history of the species involved is assumed simple (binary fission), the populations are mathematically treated as a sort of super-individual and mathematically described by ordinary differential equations.The toxicological exposure is formulated at the population level.
We consider populations that compose the ecosystem modelled by a simplified version of the Dynamic Energy Budget (deb) model [23].For an overview of different mathematical model formulations and their analysis of the unstressed ecosystems the reader may consult [18] and references therein.In our model, food is ingested at a rate given by the Holling type II functional response.A fixed portion of the ingested food is assimilated.The assimilates are used for maintenance, proportional to population size, and the remaining part for growth.
There is a long history of interest in how the properties of food webs can influence aspects of ecosystem functioning [26].Biodiversity is the variation of taxonomic life forms within a given ecosystem.Vertical diversity is the functional diversity across trophic levels along the food chain.A phenomenon related to vertical diversity is the paradox of enrichment.In mathematical models of food chains oscillatory dynamics is generally predicted when nutrients at the bottom of the food chain are abundant.In ecotoxicological studies [6] the predator-prey interaction is important, where the effects of bioaccumulation of toxicants is via trophic interactions.Horizontal diversity is the taxonomic and functional diversity within trophic levels.Much redundancy among species at one tropic level implies that species can be lost without a change of the ecosystem functioning (the redundancy hypothesis).On the other hand, if all species differ to some extent in their resources (niche complementarity), then extinction of one species can have effects.Loss of species may then lead to fewer utilised niches, stronger competition, and lower process rates, thus affecting ecosystem functioning.In mathematical models horizontal diversity is related to competitive exclusion: the number of species is not more than the number of nutrients.
Toxicological effects on individuals and populations are described by the processbased debtox [24] approach.Potential effects on physiological processes such as assimilation, maintenance, growth, and mortality are the possible mode of actions of the toxicant.The rates of these physiological processes are not constant but depend on the internal toxicant concentration.Here, this approach is extended to derive the effects on the ecosystem behaviour.Multiple nutrients and a toxicant are supplied into and removed from a spatially well-mixed chemostat at a constant rate.The kinetics of the toxicant in the organism is modelled with a first order onecompartment model where two processes are involved describing the uptake rate (from water and food) and elimination rate.When toxic food is consumed there is an additional intake rate assumed proportional to the food consumption rate.
In [29,3,7,27] bioaccumulation in food webs is studied.In these papers the toxicological effects are dealt with separately from the ecological effects.An underlying assumption is that the ecosystem is in equilibrium.We follow the mass-balance and process based deb approach [24,23,19] which allows the analysis of non-equilibrium dynamical states.In the effect module physiological parameters depend on the internal toxicant concentration.
We analyse a simple ecosystem model consisting of nutrients and populations (food chain and food web arrangements) in a chemostat.Control parameters are the dilution rate, the nutrient supply concentration and the toxicant concentration in the inflow.We call these parameters control parameters, since in our model formulation these three entities determine the environment of the ecosystem.To analyse the long-term dynamics of the resulting model we use the bifurcation analysis approach.In ecotoxicology studies the concentration of the toxicant in the environment can be taken as the bifurcation parameter.Then the value of the concentration at a bifurcation point marks a structural change of the ecosystem, e.g.extinction of a species.This indicates that chemical stressors are analysed mathematically in the same way as environmental and ecological stressors.Hence, this approach allows an integrated approach where different types of stressors are analysed simultaneously.Environmental regimes will be derived where, for a given toxic stress, no toxic effects occur and where the structural properties remain unchanged after exposure.
In a previous paper [19] we studied similar issues for a nutrient-prey-predator model.Here we extend this system by nutrient cycling and mortality of the organisms that compose the populations.The nutrients are recycled by mineralization of the dead organisms and of the products formed in physiological processes such as assimilation and maintenance.Similar ecosystem models with nutrient recycling are discussed in [4,10].To study vertical effects of toxic stress on ecosystem functioning, we introduce here a top-predator consuming the predator in the nutrientprey-predator system studied in the previous paper [19].We showed therein that the feed-back mechanism, where the toxicant is removed from the reactor as internal concentration with the removed biota while the same internal concentration affects the dynamics of the biota, is very important.Horizontal effects are studied by analysing simple food webs of which the ecological functioning has been discussed in [20,21].
2. Formulation of the model.We formulate models for a two-nutrient-twoprey-predator-top-predator system with nutrient recycling where toxicant uptake is from the water.In the inflow of the chemostat besides two nutrients, with concentrations N ri , i = 1, 2, a toxicant with concentration c r enters the reactor.Via the outflow the nutrient, toxicant, and populations are removed from the reactor, all with the same rate D, called the dilution rate.However, there is also an extra sink of toxicant, namely the toxicant stored in the organisms that leave the reactor.
The three parameters N r1 + N r2 , D, and c r are control parameters that determine the external forcing of the system.
Both prey populations are specialists and feed on their nutrient.The predator is a generalist and consumes both prey populations.Finally the top-predator consumes the predator population.For each population, products of assimilation, maintenance, and mortality are added to detritus pools in the water.The labile (easily degraded) and refractory (not easily degraded) detritus are mineralised and regenerated as nutrients.
The state variables are the nutrient density, N i (t), the biomass of the prey population, R i (t), the biomass of the predator, P (t), the biomass of the top-predator, F (t), the two variants of the detritus D L (t) (labile) and D R (t) (refractory) the toxicant concentration in the water c W (t), and the internal toxicant concentrations c R i (t), c P (t), c F (t), i = 1, 2. The state variables and parameters of the system are listed in Table 1.
2.1.Ecosystem model.The mass-balance equations for the ecosystem read The parameter m p is the maintenance rate and h p is the per capita mortality rate both for the tree trophic levels, where p ∈ {R i , P, F }, i = 1, 2. In the unstressed system these parameters are constant (m p = m p0 and h p = h p0 ) and with toxic stress they depend on the internal toxicant concentration in the population defined below.
Recycling of nutrients is via two forms of detritus.The labile compounds are products of the maintenance (excretion) process D L , and the refractory compounds are products from the assimilation (defecation) and mortality (death) processes D R (see also [10]).Both types of detritus are remineralised with different rates α L and α R where α L > α R .The partition factor of the recycled products between the two target nutrients is denoted by β i .
Table 1.State variables and parameters.The control parameters that can be experimentally manipulated are D ∈ (0, 0.5) h −1 , N ri ∈ (0, 150) mg dm −3 and c r ∈ (0, 1) mg dm −3 .m mass of toxicant [mg], l length of environment [dm] (l 3 is dimension of the volume of the reactor), L length of organism [dm] and t time [h].Here we use only the dimensionless versions of the state-variables, parameters and equations.p, q ∈ {Ri, P, F }, i = 1, 2.

Variable Description Dimension N i
Nutrient mass density Emission of the toxicant is via the in-flowing water with concentration c r and dilution rate D. We define the total toxicant concentration in the reactor as the sum of the dissolved toxicant in the reactor, both in the water and the biota where c W is the concentration in the water and c p , p ∈ {R i , P, F }, i = 1, 2 the internal concentration in the populations.
The mass-balance equation for the toxicant reads where we use the property that the internalised toxicant and the toxicant dissolved in the water leave the reactor with the same dilution rate.

2.3.
Exposer model.The one-compartment model for the internal toxicant concentration for the prey, predator and top-predator and the transport equation for the toxicant through the reactor read where k pu are the uptake rates and k pa the elimination rates, p ∈ {R i , P, F }, i = 1, 2 of the toxicant by the populations.We assume that the efficiencies for the intake of the toxicant equals the assimilation efficiency for both predator-prey interactions.
As in [12,7,19], we assume that the uptake and elimination rates are much faster than other physiological population rates.Then only the classical bioconcentration factor bcf, being the ratio of the uptake rate and elimination rate (bcf W p = k pu /k pa ), p ∈ {R i , P, F }, i = 1, 2, is involved and c p = c W bcf W p .The total toxicant concentration c T , given by (2), reads now of which the dynamics is described by (3).
2.4.Concentration-effect relationship.Direct toxic effects generally reduce population abundances by increasing mortality, increasing costs for growth or maintenance and decreasing assimilation efficiency.In the process-based debtox approach the populations are affected via a parameter alteration depending on the internal toxicant concentration: the concentration-effect relationship.In principle all population processes can be affected: assimilation, maintenance, growth, reproduction, and survival.Here we assume that the maintenance rate and the hazard rate depend on the toxicant concentrations for the prey c Ri and predator c P and top-predator are c F given by The subscripted plus operator is defined as x + = max(0, x).This illustrates that the nec parameters c pM 0 and c pH0 are threshold concentration values.These expressions are substituted in the expressions for m p and h p , p ∈ {R i , P, F }, i = 1, 2 in (1).
(1), ( 2), ( 4), and ( 6) form together the integrated model where ecological and toxicological quantities are modelled simultaneously.If not stated otherwise we consider the food chain configuration whereby 3. Bifurcation analysis.The long-term dynamics behaviour of ecosystems is described by the steady solutions of the systems of odes that described their dynamics.See [9,30,25] for an introduction to bifurcation analysis, [1,18] for applications to ecosystem models and [28] for chemostat systems.The structural stability is studied with respect to free or bifurcation parameters.Steady solutions can be equilibria where the right-hand sides of the systems of odes are zero, periodic solutions and chaotic solutions.For the composition of the ecosystem at equilibrium, the equilibrium values of the state variables are decisive.Biodiversity is measured as the number of species with positive abundance.Here we focus on equilibria, the dependency of their abundances on parameters and the associated transcritical and Hopf bifurcations.
With a bifurcation analysis, parameters are varied and critical points where the stability of a system changes and mark regions in the parameter space where the long-term dynamics change qualitatively.In principle this can be done for all parameters.In practice, depending on the purpose of the analysis, a limited number of parameters are chosen as so-called free or bifurcation parameters.For our purpose these are three parameters, namely the control parameters N r , D, and c r that determine the loading of the system (N r for ecology, c r for toxicity and D both).Computer packages such as auto [5] and locbif [16], can be used to calculate the bifurcation curves.
The results of the bifurcation analysis of the ecosystem are given in bifurcation diagrams.These diagrams directly give information about the qualitative long-term dynamic behaviour (equilibrium, oscillatory or chaotic) at specific environmental conditions (nutrient input, dilution rate and toxicant input).
3.1.Discontinuous Jacobian matrix.The expressions in (6) show that the affected physiological parameter is continuous with respect to the internal concentration but that its derivative with respect to this concentration is discontinuous at the nec concentration.The dynamic behaviour in time is well defined also when the internal concentration crosses the nec concentration.The right-hand sides of the odes satisfy the Lipschitz condition with respect to the state variables.Hence, we have uniqueness of the initial value problem.
When c p = c pM 0 or c p = c pH0 , p ∈ {R, P, F } (see Eq. ( 6)) the partial derivatives of right-hand side functions of the ode system with respect to the state variables are discontinuous with respect to bifurcation parameters.At this point the Jacobian matrix is discontinuous.As a consequence, when a parameter is varied the eigenvalues can jump from one half of the complex plane into the other half without a smooth transition as with the Hopf bifurcation.As a result the bifurcation curves at which the long-term dynamics changes qualitatively do not have to be smooth curves at points where the interior concentration crosses the nec.
We show first that the transcritical bifurcations T C j , j = 1, 2, 3, are not affected when they describe invasion of a species, for instance the predator into the nutrientprey system.The partial derivatives of (1c) with respect to the state variable (here the predator population P ) and substitution of the equilibrium value (here P * = 0, F * = 0) reveals that only one term is non-zero, namely the diagonal term . This expression equals the per capita growth rate of the predator population when it is rare.The parameter values where this rate equals zero marks a transcritical bifurcation (see also [28]).Since the rate is continuously dependent on the parameters, this curve is also continuous with respect to the internal concentration at levels across the nec level.This holds in a similar way for invasion of the prey population R and the top-predator population F .

3.2.
Calculation of the nec-isocline.The environmental regime where the applied toxic stress has no effects is determined for equilibrium states.The points where c p = c pM 0 or c p = c pH0 , p ∈ {R i , P, F } form a codim-1 curve in the N r , D bifurcation diagram.This curve is denoted by I j , where j = 1, 2, 3 is the food chain length after invasion and is referred to as nec-isocline.
In equilibrium we have c T = c r where c r is assumed to be given.Then the nec-isocline I 3 is fixed by the relationship where c W = c F M 0 /bcf W F when for instance c F = c F M 0 .To continue this curve, the equilibria of the ecosystem (1) are continuated with N r as a bifurcation parameter where D is implicitly fixed by (7).
A standard bifurcation package can be used when we treat D as a state variable instead of a parameter of which the dynamics is described by When N r is varied as a bifurcation parameter the equilibrium values for D describe the wanted curve.When instead of the toxicant concentration c r , the dilution rate D is fixed, a similar procedure can be used by continuing system (1) together with and N r as bifurcation parameter.
4. Analysis of food chain systems with recycling.In this section we analyse the long-term dynamics of the ecosystems with increasing food chain length.In Table 1 and Table 2 we give the parameter values for the food chain.The physiological parameter values are those for a bacterium-ciliate-carnivorous ciliate system.

4.1.
Effects with constant dilution rate.In this subsection we keep the dilution rate D = 0.02 constant and vary both the nutrient and toxicant concentrations in the influent: N r and c r .The results are shown in Figure 1.Only equilibrium long-term dynamics is discussed.For the nutrient-prey system the transcritical bifurcation T C ± 1 and the tangent bifurcation T 1 (Fig. 1A) occur.Bi-stability occurs between the tangent bifurcation T 1 and the transcritical bifurcation T C + 1 .
Table 2. Parameter set for the two-nutrient-two-prey-predator system.These values were also used in [21].Observe that the maximum ingestion rates I R 1 P and I R 2 P are taken equal.For all calculations we assumed that both nutrients are supplied with the same density:  1B.For the nutrient-prey-predator-top-predator system the same bifurcations denoted by T C − 3 and H − 3 , Figure 1C, occur.The curves I j with j = 1, 2, 3 the food chain length, are the nec-isoclines calculated by continuing system (1) together with (9).These curves describe properties at the population and the ecosystem level and they are related to the No Effect Concentration (nec) at the individual level [23,24].Due to deteriorating effects of the toxicant the bifurcation curves for the longer food chains are at higher nutrient levels in the diagram than for the unstressed systems with c r = 0 on the horizontal axis.
Figure 1D is the superposition of the three diagrams for the different food chain lengths.From this figure we conclude that the single population system is less resistant to toxic stress than a food chain.For instance, for inflow concentrations higher than, say, four there are two end-points.Suppose at the initial state with N r = 50 there is a prey population.Now this system is exposed to the toxicant.Then depending on the initial conditions the prey population can go extinct.In a second situation before exposure to the toxicant a predator is introduced into the reactor.In this case the nutrient-prey-predator system establishes itself.When this system is exposed now to the toxicant with c r = 4 the system is resistant.

4.2.
No effect regimes.We are interested in a region in the two-parameter diagram with N r and D as bifurcation parameters, where the toxic stress has no effects: the dynamic system is equivalent to the unstressed system and consequently posseses the same dynamics.The boundaries of the regions in the two-parameter diagram where no effects occur are the nec-isoclines calculated by continuing system (1) together with (8).For the nutrient-prey system, in Figure 2A the curve is denoted by I 1 .On the right-hand side of this curve I 1 the internal concentration is lower than the nec and therefore the biomass abundances equal those in the unstressed system, while on the left-hand side the abundances are affected.We call the region at the right-hand side the No Effect Region (ner).
For the prey-predator system, in Figure 3A the nec-isocline curve I 2 connects the two end-points of the two branches of the Hopf bifurcation curve H − 2 .Observe that this curve for the nutrient-prey-predator system terminates at the point where that curve for the nutrient-prey system intersects the transcritical bifurcation curve  N r and dilution rate D as bifurcation parameters for the ecosystem with toxicity stress in the chemostat system c r = 1. A. The solid line is the nec-isocline.The bifurcations belonging to the nutrient-prey system are also shown in Figure 2. B. A two-parameter bifurcation diagram for the nutrient-prey-predator system.Between the two thick lines there is a structural change of the ecosystem due to toxic stress.

T C −
2 .This is in agreement with our finding in Section 3.1, namely that the transcritical bifurcations T C j , j ∈ {1, 2, 3} are not affected when they describe invasion of a species, for instance the predator into the nutrient-prey system j = 2.For the nutrient-prey-predator-top-predator system Figure 4A is similar to that of the nutrient-prey-predator system.The nec-isocline is labeled I 3 .
Notice that because the nec values are the same for the three populations there is a single nec-isocline curve.In general, when the nec values differ for different trophic levels these curves have to be calculated for each value of the nec values.These curves are plotted jointly in the diagram, the cross-section of all regions equals the ner's.

Resistance conditions.
Here regions in the bifurcation diagrams are considered where an increase of the toxic stress does not lead to a change of the structure of the ecosystem.These regions are bounded by bifurcation curves for c r = 0, the control case, and for the stressed case, here c r = 1.Mathematically structural changes are related to bifurcations and biologically population extinction or periodic instead of equilibrium dynamics.
In [8] the term resistance, is defined as: "staying essentially unchanged despite the presence of disturbance."When "essentially unchanged" means "the same longterm dynamic behaviour as the reference state (stable equilibrium, limit cycle of even chaotic attractor)" and "presence of disturbance" as "the toxic stress," then this definition applies to the regions in the bifurcation diagrams where increase of the toxic stress does not lead to a change in the structure of the ecosystem.
Figure 2B shows that the resistance region is bounded by the transcritical bifurcation curve T C ± 1 where c r = 1.The effects of toxic stress for the region in the parameter space between the two transcritical bifurcation curves for c r = 0 and c r = 1 (light grey) is dramatic.There is a large region with bi-stability.Figure 3B shows that the toxic effects on the the nutrient-prey-predator system is rather small.Only at low nutrient enrichment there is effect.From Figure 4B we conclude that this is also the case for the nutrient-prey-predator-top-predators system.In  N r and dilution rate D as bifurcation parameters for the ecosystem with toxicity stress in the chemostat system c r = 1. A. The solid line is the nec-isocline.The bifurcations belonging to the nutrient-prey and nutrient-prey-predator system are also shown in Figure 3. B. A two-parameter bifurcation diagram for the nutrient-prey-predator-toppredator system with toxicity stress, c r = 0 and c r = 1, in the chemostat.Between the two thick lines there is a structural change of the ecosystem due to toxic stress.
this case the resistance region is bounded by the transcritical bifurcation curve T C 2 and T C 3 for c r = 0 and c r = 1 (light grey).There is a region where the predator and subsequently the top-predator go extinct when the toxic stress is increased from c r = 0 to c r = 1.In Figure 4B this region is the small (dark grey) area.
Notice that in order to determine the region of resistance, the diagram for the stressed system is compared with that for the unstressed system.The property of resistance depends on both situations, unstressed and stressed, while with the calculation of the no effect region only the stressed situation is considered.The end-points of the two branches of the supercritical Hopf bifurcations are connected by the nec-isocline I 2 .
5. Pseudo-Hopf bifurcations.In this section we investigate the consequences of a discontinuous Jacobian matrix with respect to a parameter.Figure 5 is a detail of the two-parameter diagram for a stressed prey-predator food chain.
In Figure 6 bifurcation diagrams are shown with D as free parameter for different N r values: N r = 72 (Fig. 6A), N r = 72.2 (Fig. 6B), N r = 72.5 (Fig. 6C) and N r = 73.5 (Fig. 6D).Below the point I 2 at which the Jacobian is discontinuous with respect to D, the equilibrium is unstable.Above that point there is bi-stability of the stable equilibrium and a stable limit cycle.With N r = 72 (Fig. 6A) this limit cycle originates from the lowest Hopf bifurcation point and becomes unstable at a tangent bifurcation point for this limit cycle denoted by T c .For N r = 72.2 (Fig. 6B) we passed the cusp bifurcation C in Figure 5 where two tangent bifurcations T c collide.Now between two tangent bifurcations three limit cycles coexist whereby the intermediate limit cycle is unstable and the other two are stable.
With N r = 72.5 two tangent bifurcations disappear again but in another manner than at the cusp bifurcation.With respect to the situation at N r = 72.2, the stable limit cycles that originate at the Hopf bifurcation with the lowest nutrient inflow are now connected with the Hopf bifurcation with the highest nutrient inflow.Between the intermediate Hopf bifurcation and the point I 2 there is a stable equilibrium and an unstable limit cycle which forms the separatrix.Just above the point I 2 a stable limit cycle originates that merges with the unstable limit cycle at the tangent bifurcation T c .For N r = 73.5only two Hopf bifurcations are present with a stable limit cycle in the enclosed interval.
We conclude that point I 2 acts as a subcritical pseudo-Hopf bifurcation.Figure 5 shows that the tangent bifurcation of the limit cycle undergoes a cusp bifurcation point C and this explains why the bifurcations in the diagrams presented in Figure 6 possess different patterns.
6. Analysis of food web systems with recycling.
6.1.Unstressed food web system.For the unstressed food web model described by (1) where m p = m p0 and h p = h p0 , p ∈ {R i , P, F }, i = 1, 2 given in Table 2, the calculated bifurcation diagram is shown in Figure 7A.Only curves for the nutrients-preys-predator system are shown, not of the sub-systems.The transcritical bifurcation curve T C 2 and the Hopf bifurcation curve H − 2 are the only important bifurcation curves.On the right-hand side of the T C 2 curve the predator can invade the nutrients-preys system at stable equilibrium.In the region on the right-hand side of the Hopf bifurcation curve, the positive equilibrium is unstable and the system shows oscillatory behaviour.
Qualitatively the diagram resembles that of the unstressed food chain nutrientprey-predator system.This is explained as follows.It is easy to show that if both prey populations are identical there exists an equivalent food chain model described by where the saturation constant is twice that for the nutrient-prey interactions.6.2.Stressed food web system.Figure 7D gives the diagram for the stressed food web system with c r = 1.Only for low dilution rates and nutrient levels the internal toxicant concentration becomes larger than the nec value.In Figure 7B and Figure 7C the diagrams are shown for food chains (no top-predator) with the prey R 2 and prey R 1 population.respectively.These diagrams show the extreme cases when the applied toxic load is lethal for one of the two prey populations.Comparing Figure 7D with Figure 7B and Figure 7C shows that the effects of toxic stress on the food web model resembles that for nutrient-prey-predator food-chain models.These results suggest redundancy when ecosystem functioning is considered while structure and biodiversity are affected.
7. Discussion and conclusions.The ecosystem model ( 1) is a straightforward extension of the model for the nutrient-prey-predator system in [19], but now with nutrient recycling via two forms of detritus, another nutrient and prey population and a top-predator.Comparison with the results reported there and obtained here show an essential effect of recycling for low dilution rates.
In [6] many direct and indirect effects on the ecosystem level reported in the literature are classified.Direct effects lead to reduced population abundances due to negative effects on the population growth.They are measured in single species experiments.In multi-species ecosystems a toxicant can exert indirect effects on resistant populations by a number of mechanisms, such as cascading effects in food chains by altering trophic or competitive interactions.Indirect effects enhance, mask, or spuriously indicate direct toxic effects.In our approach ecological, physical/environmental, and toxicological factors affect ecosystem process rates.As a result direct and indirect effects on the ecosystem are taken into account as well as the feed-back mechanisms between the concentration of dissolved toxicant and the ecosystem state.For instance, despite high toxicant concentrations in the influent (ten times the nec-values), the water toxicant concentrations in the reactor need not be high.This is the result of the fact that the supplied toxicant-free (unstressed) nutrient is converted in the reactor into biomass which takes up the toxicant and removes it partly from the reactor together with the removed individuals.In large regions of the bifurcation diagram the concentration of dissolved toxicant, c W , in the effluent is low even when c r is high in the influent.In the literature generally the concentration toxicant in the water c W is assumed constant and known in advance, see [7].
In [24] direct and indirect effects on the individual or population level are discussed for the reproduction.Since processes, assimilation (i.e. the combination of feeding and digestion), growth, maintenance, mortality and reproduction are interlinked, changes of any of these processes can result in changes in the reproduction.Here the organisms propagate by binary fission and maintenance and mortality as modes of action for all populations, are directly affected.
In [11] the bifurcation technique is used to calculate the effects of toxicity on ecosystem functioning.In this model the death rate depends linearly on the toxicant concentration.The resulting expression equals the concentration-effect relationship (6b) where the nec concentration equals zero.The toxicant concentration in the system is taken constant.Hence no emission, transport, and exposer models are formulated.Therefore the feed-back, that is the dependency of the toxicant concentration on the ecosystem state of the system, cannot be considered.
Bifurcation analysis techniques for the assessment of toxic effects on ecosystem functioning fits well in the modelling philosophy that simplified models will not provide precise forecasts for real ecosystems.It enables us to understand a range of potential outcomes by dealing with extreme situations.By reformulation of the set of equations where a state variable is treated as a parameter and a parameter as a variable, continuation techniques can be used to determine under which toxic stress levels or environmental conditions the ecosystem is resistant.Open systems (chemostat reator), as well as closed systems (batch reactors) without mass exchange with the environment, can be analysed with this approach.This facilitates the use of the same modelling formalisms for indoor laboratory and outdoor natural field systems.

AFigure 1 .Figure 2 .
Figure1.Two-parameter bifurcation diagrams with total nutrient inflow Nr and toxicant inflow concentration cr as bifurcation parameters for stressed food chains where the dilution rate is fixed D = 0.02.The curves I j , with j = 1, 2, 3 the chain length, are nec-isocline curves .A. Nutrient-prey system, B. nutrient-prey-predator system, C. nutrientprey-predator-top-predator system and D. superposition of the diagrams.

Figure 5 .
Figure 5. Bifurcation diagram for region where the two branches of the Hopf bifurcations H −2 terminate due to a discontinuous Jacobian matrix and which are connected by the nec-isocline.Point C is a cusp bifurcation for the tangent bifurcation of limit cycle Tc.

Figure 6 .
Figure 6.One-parameter bifurcation diagrams for the prey-predator system with dilution rate D as free parameter, where A. N r = 72, B. N r = 72.2,C. N r = 72.5 and D. N r = 73.5.The equilibrium values or the extreme values of the nutrient density in the chemostat N are plotted.Points H −2 are Hopf bifurcation points where the equilibrium becomes unstable.Point I2 is the point where the internal concentration equals the nec-value and where the Jacobian matrix is discontinuous.Point T c is a tangent bifurcation point for this limit cycle.

AFigure 7 .
Figure 7. Two-parameter bifurcation diagram with total nutrient inflow and dilution rate D as bifurcation parameters for different compositions of the food web.For the stressed systems the toxicant concentration in the inflow is cr = 1.The drawn line is the nec-isocline curve.A. Unstressed food web system c r = 0, with N r1 = N r2 and dilution rate D as bifurcation parameters.B. Stressed system, c r = 1, with only the competitor R 2 , and N r2 and dilution rate D as bifurcation parameters.C. Stressed system, cr = 1, with only the prey R1, and Nr1 and dilution rate D as bifurcation parameters.D. Stressed system, cr = 1, with prey R 1 and competitor R 2 , nutrient inflow N r1 = N r2 and dilution rate D as bifurcation parameters.