Non-linear analysis indicates chaotic dynamics and reduced resilience in model-based Daphnia populations exposed to environmental stress.

In this study we present evidence that anthropogenic stressors can reduce the resilience of age-structured populations. Enhancement of disturbance in a model-based Daphnia population lead to a repression of chaotic population dynamics at the same time increasing the degree of synchrony between the population's age classes. Based on the theory of chaos-mediated survival an increased risk of extinction was revealed for this population exposed to high concentrations of a chemical stressor. The Lyapunov coefficient was supposed to be a useful indicator to detect disturbance thresholds leading to alterations in population dynamics. One possible explanation could be a discrete change in attractor orientation due to external disturbance. The statistical analysis of Lyapunov coefficient distribution is proposed as a methodology to test for significant non-linear effects of general disturbance on populations. Although many new questions arose, this study forms a theoretical basis for a dynamical definition of population recovery.

Introduction

Recovery, resilience and disturbance

Ecosystems are increasingly affected by human activity. Structural as well as chemical impacts threaten ecological populations, communities, the connected ecosystem functions, processes and services [1]. The questions how the constituents of ecosystems react to those impairments and whether they can recover are matter of many ecological research studies (e.g., [2], [3], [4], [5], [6], [7], [8]).

The term recovery is often stressed in ecology and monitoring. In its broadest sense recovery is defined as the recurrence of a system to a preceding status after disturbance, often quantified as the time to achieve this pre-disturbance status. Thereby the concept of recovery is used to quantify how severely a system is affected by a disturbance but also to assess whether a system is enduringly negatively affected.

Recovery constitutes a part of the much debated concept of ecological stability which is connected to various other ecological aspects as resilience [9], [10]. Like for recovery different definitions of the term resilience exist (e.g., [11], [12]). Engineering resilience is defined as the time that a system takes to return to equilibrium following perturbation [13], whereas ecological resilience is defined as the amount of perturbation a system can withstand before it moves into a different basin of attraction or stability domain [14]. Often resilience is simply defined as the capacity of the system to return to a given state after a disturbance [15]. In this last case recovery and resilience are used synonymously (resilience  =  rate of recovery).

As resilience reflects the ability of a community or a population to withstand and recover from disturbance it is assumed a fundamental property of ecological systems and thus has been proposed as an indicator of ecological health and integrity [16]. Within this context disturbance is broadly defined as any relatively discrete event that disrupts ecosystem, community, or population structure and changes resources, substrate availability or the physical environment [17].

It is generally accepted that most species communities are constantly exposed to natural disturbance resulting in changes of community structure and life history characteristics of the constituting species. From the idea that many of the environmental parameters that determine resilience of populations and communities to natural disturbance may also influence the response to chemical stressors [16] it can be concluded that exposure to contaminants or other anthropogenic impacts alters the ability to withstand and recover from natural disturbance via reduced resilience.

Two different types of disturbance can be distinguished. On the one hand pulse disturbance [18] results in prompt change in species abundance within a community (e.g. forest fire or chemical spill that reduces densities of certain species). Effects of pulse disturbances are generally assumed to recede after some time so that the system can return to its initial equilibrium state. On the other hand press disturbance causes permanent change in species abundance and often leads to the loss of some taxa and establishment of an alternative community (e.g. increased temperature associated with global climate change or continuous input of toxic material into a system). Until the elimination of the stressor a community under continuous press disturbance is generally not accepted to return to its original condition [16].

Dynamics of age-structured populations

The fact that age structure is an important feature of natural populations has been deduced in many studies based on population models [19], [20], [21]. Many animal species, like crustaceans, have life-histories composed of a sequence of stages within which their individuals' characteristics are broadly similar to those of other individuals in the same stage and markedly different from those of individuals in other stages [22]. Within these complex life-cycles the vital rates (rates of survival, development, and reproduction) almost always depend on individuals' age, size, or development stage. It has become widely recognized that this diversity can have a significant influence upon dynamics at the population level and that in order to understand the dynamics of most populations it is necessary to take this diversity adequately into account [23], [24].

Such variation can be described with the help of time discrete nonlinear age-structured population models (for examples see [25], [26], [27], [28], [29], [30]) which have been widely studied in the biomathematical literature [31] (e.g., [32], [33], [34], [35], [36], [37], [38]). A large number of theoretical analyses and numerous applications to ecological systems have proven the relevance of age-structured population models (e.g. [22], [28], [31], [39], [40], [41], [42], [43], [44], [45], [46]).

Even though these models are of simple structure it was shown that they display all kinds of dynamical behaviour, like equilibria, periodic cycles, and unpredictable chaotic dynamics [40], [47], [48], [49], [50], [51], [52], [53], [54], [55]). The occurrence of this complex dynamic behaviour is related to the fact that minor changes in parameter or initial values can strikingly change the dynamic behaviour of systems [31]. Model dynamics were found to be very sensitive to changes in life-history parameters as well as environmental perturbations [28]. In more general terms of system dynamics it was shown that trends and fluctuations in populations are determined by complex interactions between extrinsic forcing and intrinsic dynamics [46].

Daphnia magna as model organism

Daphnia magna is a Crustacean in the order of Cladocera. In terms of assessing the effects of human contaminants D. magna exhibits an ideal test organism for this purpose as it constitutes a central component in aquatic food web structure [56], [57], [58] and is known to be sensitive to a multitude of xenobiotics (organophosphates, heavy metals, organochlorines, pyrethroids etc.) [59], [60], [61]. Due to its ease of culture and its clonal method of reproduction [56], [62] it is commonly used as a model organism in ecological risk assessment of pesticides, biocides, and industrial chemicals [29].

Within a Daphnia population three age classes (roughly equal to size classes) representing the different life stages can be distinguished, defined by means of age and fertility: neonates (less than 24 h old), juveniles (at least 24 hours old but non-mature), and adults (mature). These age classes differ in terms of their toxicological sensitivity in mortality, reproduction and developmental time. Effects within the population structure due to xenobiotics emerge from these age class specific sensitivities, e.g. reproduction rate.

Such a daphnid population exhibits a simple example of an ecological system with multi-stage interactions of age-classes an, n = {1, 2, 3}. Due to its reduced dimensionality and the differential interactions within it, the population comes close to the classical example of a 3-dimensional dynamical system like the Lorenz attractor. This suggests to analyse the effects of disturbance (effect  =  altered intrinsic rate of natural increase) of this system in terms of non-linear dynamics.

DCA as a model substance

3,4-Dichloroaniline (DCA, CAS-Nr.: 95-76-1, EEC-Nr.: 202-448-4) is an intermediate used for the synthesis of plant protection products (e.g., linuron, diuron, propanil) and as an azo dye for polyester fabrics. Release into the environment mainly occurs during use of plant protection agents after biotransformation or simply as an impurity of these substances [63]. DCA is known to affect survival and reproduction in D. magna already in low concentrations, resulting in reduced offspring number and increased number of aborted eggs in chronic toxicity tests [64], [65], [66]. In higher concentrations DCA has a lethal effect, causing increased mortality rates in acute toxicity tests [67]. Thus a strong effect of DCA on Daphnia population dynamics is obvious [68] and it exhibits a useful model substance to investigate effects on population level.

Non-linear time series analysis

The autocorrelation function allows for determining the non-whiteness of data and detecting periodic components in data. Because the autocorrelation function and the power spectrum are related by the Fourier transform, they are mathematically equivalent, and contain the same information [69].

Based on Takens' theorem [70] time-delay embedding allows for the reconstruction of the system's phase space and its attractors [71], [72], [73]. Scalar observations x(n) and their time delays x(n + kT); k = 1,2,…,m−1 are used to make m-dimensional vectors y(n) = [x(n),x(n+T),…,x(n+(m−1)T) whose components provide a coordinate system in which one can identify the attractor structure associated with the observations (T = time delay)[74].

Recurrence plots are a graphical tool for measuring the time constancy of dynamical systems [75]. They represent the recurrence of the phase space trajectory to a certain state, which is a fundamental property of deterministic dynamical systems [76], [77]. The analysis of recurrence plots allows for the detection of periodicity and regression [75] as well as typical transitions, e.g., bifurcation points occurring in complex systems [78].

A Poincaré section (also called a surface of section, Poincaré map or first recurrence map) [79] creates a projection of intersection points of an attractor and a hyper-plane in phase space [77]. It differs from a recurrence plot in that space, not time, determines the plot points. The hyper-plane is adjusted to cut the attractor transversal, thereby depicting the orientation of the n-dimensional attractor as a cloud of points in (n-1)-dimensional space.

Lyapunov exponents (λ) are a measure for exponential divergence of initially nearby trajectories in embedded phase space and are thus a hint for chaotic dynamics within a dynamical system [80]. If a dynamical system exhibits at least one Lyapunov exponent larger than zero it is assumed to show sensitivity to initial conditions, a prerequisite for chaotic dynamics [81], [82]. Calculating maximal global Lyapunov coefficients instead of local coefficients has the advantage of averaging out local divergence rates [83].

Aims

In this paper the authors present an analysis of the dynamics of a model-based age-class structured population of Daphnia magna exposed to different concentrations of a chemical stressor (3,4-Dichloroaniline, DCA). To support the understanding of the effects of chemical disturbance the study focused on the non-linear characteristics of these dynamics and the question whether there is evidence that disturbance may lead to non-reversible alterations of population dynamics. This was achieved by quantification of periodic and chaotic parts of population dynamics, using methods from non-linear analysis and concepts from dynamical systems theory.

The hypothesis to be tested in this study was that anthropogenic stressors can alter the resilience of an age-structured population and thus have a negative influence on the ability to respond to further natural disturbance, thereby enhancing the risk of population extinction.

The Lyapunov was tested as an indicator to detect disturbance thresholds leading to a discrete change in the investigated population dynamics. Interpreting the results in terms of resilience, the study was intended to contribute to a dynamical definition of population recovery.

Materials and Methods

Model-based study

In principle there are three ways to reconstruct dynamics from non-linear systems: (1) analysis of system equations, (2) analysis of empirical data and (3) analysis of model-based data. In stress ecology simulation models are often used to predict recovery time [5], [29], [30], [84], [85], [86], [87], [88].

Analysis of system equations is more a theoretical task and not applicable in this context. Analysis of empirical data is often not optimal in the case at hand as the time series available is too short to analyze non-linear dynamics sufficiently and detect chaotic dynamics. Additionally empirical data are noisy and to find deterministic chaos the deterministic skeleton must be extracted from the data (which requires long data series as well). Thus, in this study we used a validated Daphnia population model (details see below) to simulate the empirical data. By this means we were able to (1) extend the observational time frame, as the model is validated to predict Daphnia population dynamics, (2) extract the deterministic skeleton, eliminating the confounding noise (we do this by setting chance constant) and (3) take into account the diversity within the complex life-cycle to understand the dynamics of Daphnia magna [23], [24].

Simulated data sets

The time discrete data used in this study were generated by use of an individual based population model for Daphnia magna (IDamP, [29]). The model predicts the population dynamics based on individual life cycles which include feeding rate, growth, development, reproduction and survival processes for which food conditions and population density (crowding effects) are the main drivers. The life cycle process is represented as descriptive regression models, based on a large dataset from D. magna life cycle tests.

IDamP reliably predicts population growth and reproduction for D. magna under laboratory conditions using the algae Desmodesmus subspicatus as a food source. It was recently shown to be able to extrapolate the individual effects of DCA to the population level as well as for other toxicants [30], [89], [90]. Besides mortality the only sub-lethal effects for DCA considered in the model is the inhibition of reproduction. The interactions between individuals are due to crowding and competition for food. To reduce the model stochasticity all individual properties, except a maximum lifetime assigned randomly to the initial individuals and the daphnids born during the simulation, were switched off. This maximum lifetime value was not changed anymore during the simulation. Raw data from the simulations can be obtained from the authors on individual request.

Five concentration levels for DCA were used in the simulations (0 µg/l, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l), held constant over the simulation time frame of 365 days thus representing a press disturbance. Based on experimental laboratory data concentration-response relationships for DCA on the individual level had been calibrated in an integrated toxicity module in the IDamP model and validated on population level before [30].

As the DCA concentrations are far below the acute mortality concentrations for daphnids (48 h LC50: 0.23 mg/l [63], [67]), it can be assumed that the effect will mainly occur on the adult daphnid reproduction rate.

A Monte-Carlo set of 100 simulation runs was conducted using the same parameter set as published in [30]. Differences in the 100 simulation runs resulted from differences in initial conditions. Before starting the simulations 5 neonate and 3 adult daphnids were introduced into the system. These 8 individuals were assigned random but fixed values for maximum lifetime.

As the chaotic parts of the population dynamics would be lost in non-linear analysis when based on mean values from these simulation runs, non-linear dynamical characteristics were calculated for each simulation run separately. Statistics were subsequently calculated from the resulting distributions of non-linear characteristics.

The overall population predicted by the model was additionally divided into three age classes (size classes): neonates (<1.4 mm), juveniles (≥1.4 mm and <2.6 mm) and adults (≥2.6 mm). In this way the data structure allowed for statistical analysis of population age class dynamics.

Modelling non-linear dynamics

Data analysis was based on the time series data described above. To exclude the effects from transient dynamics, the first 50 data points were omitted from all analyses. First, 3-dimensional phase space was plotted from the raw data of neonate, juvenile and adult daphnids. Secondly, exploratory analysis was performed on the time series data by use of linear autocorrelation functions. Finally, after embedding the time series, non-linear analysis was applied by means of recurrence plots, Poincaré sections and global maximum Lyapunov exponents. Data visualisation and analysis was performed by use of the software package Tisean [83], [91] in MATLAB (Ver. 8.2.0.701, release R2013b) on a Linux-based system (Debian testing Jessie). The MATLAB source code S1 is accessible from the supporting information.

In a 3-dimensional plot neonate, juvenile, and adult Daphnia abundances were plotted for each DCA-concentration. Trajectories were assessed in terms of compactness and orientation in phase space. Autocorrelation analysis was performed with function corr from the Tisean package to assess the amount of linear structures within the data. To find the best unfolding of the time series in the reconstructed phase space a series of time-delay embeddings was performed by means of the function delay from the Tisean package and the optimal parameters were chosen. Based on this assessment, finally all subsequent analyses were conducted with embedding dimension m = 2 and delay d = 1 (except the Poincaré section, which was performed with m = 3). Recurrence analysis [78] was performed to quantify the number and duration of recurrences of the system in its state space. We used the function recurr from the Tisean package, using m = 2 as embedding dimension, d = 1 as delay and  = 1 as neighbourhood. Poincaré sections were created with the function poincare from the Tisean package applying an embedding dimension m = 3 and delay d = 1 to find a geometric depiction of the trajectories in a lower-dimensional space. Quantification of attractor characteristics was done by calculation of global Lyapunov exponents. Maximum global Lyapunov coefficients for all 100 simulation runs per age class and concentration level were calculated with the function lyap k from the Tisean package [80], [81] using an embedding dimension m = 2 and delay d = 1, revealing a distribution of Lyapunov coefficients. Values larger than 1 were assumed to be a strong evidence for chaotic dynamics within the time series.

Statistical analysis

From the distributions of 100 Lyapunov coefficients mean values and standard deviations per age class and treatment level were calculated and tested for significant differences by means of a Wilcoxon signed-rank test (function wilcoxon.test in R, two.sided, conf.level = 0.95, [92]).

Results

All results are described for the overall population as well as the age classes (neonates, juveniles and adults). Figures are only shown for the overall population. The according figures for the age classes can be found in the supporting information.

The mean Daphnia abundance from all 100 simulation runs (including the 95% confidence interval) for the overall population in the control treatment is exemplarily shown as time series data in Fig. 1 (abundances for neonates, juveniles and adults can be found in the supporting information, Fig. S1). It was assumed that the population and its age classes existed under equilibrium conditions under control conditions from t = 50 on.

Figure 1

Daphnia abundance from all 100 simulation runs for the overall population in control.

Solid line: mean abundance, dashed line: 95% confidence interval.

Population dynamics were quite comparable for control and DCA concentrations lower than 40 µg/l (Fig. 2). For the highest DCA concentration level the dynamics were obviously different.

Figure 2

3-dimensional phase spaces from raw data of simulation run 1 for neonates, juveniles and adults.

From top left to bottom right, row-wise: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l.

In general the autocorrelation functions (Fig. 3) showed that with increasing DCA concentration no significant change in structure occurred. Comparable periodic dynamics of fluctuating correlation coefficients were observed up to 20 µg/l. This effect was especially strong for juveniles and adults, while in the neonate age class the effect was clearly less pronounced (see Fig. S2 in supporting information). In contrast the autocorrelation function for 40 µg/l showed a very strong correlation structure with clearly increased amplitude and decreased frequency.

Figure 3

Autocorrelation functions of simulation run 1 for the overall population.

From top to bottom: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp.

From time-delay embedding (Fig. 4) it was concluded, that the trajectories for the overall population settled in the middle of the phase space. Especially for neonates and juveniles (see Fig. S3 in supporting information) the trajectories settled within a small part of the phase space nearer to the origin for the highest DCA concentration (40 µg/l).

Figure 4

Embedded time series of simulation run 1 for the overall population.

From top left to bottom right, row-wise: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. (d = 1, m = 2).

Increased pattern formation was observed with increasing DCA concentrations in the recurrence plots (Fig. 5), especially for 40 µg/l. This was also observed for the age classes (see Fig. S4 in supporting information). For neonates a very strong lineal pattern occurred in the highest concentration level, while for juveniles and adults this pattern formation was less pronounced. Especially for juveniles additional circular structures were observed in the highest concentration level.

Figure 5

Recurrence plot of simulation run 1 for the overall population.

From top left to bottom right, row-wise: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. (d = 1, m = 1,2, ε = 1).

Poincaré sections did not show large differences in dynamics for treatments from control up to 20 µg/l (Fig. 6). For 40 µg/l the attractor orientation in phase space abruptly changed in location (shift towards origin) and cross sectional area decreased. This effect was especially pronounced for juveniles and adults (see Fig. S5 in supporting information). For neonates the attractor was always located near the phase space origin.

Figure 6

Poincaré sections of simulation run 1 for the overall population.

From top left to bottom right, row-wise: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. (d = 1, m = 3, C = 1).

Table 1 shows mean values and standard deviations for maximum Lyapunov coefficients calculated from the Lyapunov distribution for the overall population in the treatment levels. All maximum Lyapunov coefficients for all treatment levels were larger than zero (Fig. 7). Lyapunov coefficients of DCA treatments up to 20 µg/l were only slightly smaller than for control, while for 40 µg/l a large decrease was observed.

Table 1

Mean values and standard deviations for maximum Lyapunov coefficients per treatment level for the overall population with n = 100.

	Control	2.5 µg/l	5 µg/l	10 µg/l	20 µg/l	40 µg/l
mean	0.503	0.477	0.470	0.467	0.459	0.261
stdev	0.0471	0.0492	0.0485	0.0528	0.0507	0.0336

Figure 7

Boxplots for maximum Lyapunov coefficients for DCA treatments for the overall population.

Boxes show mean value (horizontal line), 25 and 75% quantile (box edges), 5 and 95% quantile (whiskers) and outliers>2 sd from mean (circles).

Table 2 shows the results from the Wilcoxon signed-rank tests on statistical differences between mean maximum Lyapunov coefficients for all treatment levels. The mean maximum Lyapunov coefficients for 40 µg/l were significantly different from all other treatment levels. Also, all mean maximum coefficients for the treatment levels were significantly different from the control.

Table 2

p-values from Wilcoxon signed-rank tests on statistical differences between mean maximum Lyapunov coefficients per treatment level for the overall population.

	2.5 µg/l	5 µg/l	10 µg/l	20 µg/l	40 µg/l
Control	4.853E-4 *	3.489E-6 *	1.137E-6 *	2.267E-9 *	† *
2.5 µg/l		0.2401	0.1183	6.919E-3 *	† *
5 µg/l			0.5842	0.09201	† *
10 µg/l				0.316	† *
20 µg/l					† *

: p<2.2e − 16, significant differences with p<0.05 are marked by an asterisk.

Discussion

From the time series data in Fig. 1 it was concluded that under control conditions (no disturbance) the population and its age classes existed under equilibrium conditions. The existence of equilibrium is prerequisite in the attempt to define recovery because only if population structure is relatively stable, the population can return to defined pre-disturbance conditions if given sufficient time following the disturbance [16].

Generally the system seemed to converge to a quasi-stable state in the long run (see Fig. 2). For concentrations lower than 40 µg/l the system seemed to settle in a relatively small phase space region. This region occurred well bounded (a prerequisite for chaos), but the trajectory seemed to follow a complex attractor, exhibiting considerable population fluctuations. For the highest DCA concentration level fluctuations were much smaller compared to lower DCA concentrations. The trajectory settled within a much smaller phase space region, nearly approaching the zero abundance level, especially for the adult age class. These dynamics clearly exhibited an increased extinction risk for the overall population if this effect would have lasted too long (although this does not pose a problem to the age-class itself in the first place, as the classes are refilled by aging individuals).

For the overall population (as well as for all age classes) in the control and all treatment levels a sinusoid like oscillation was observed (see Fig. 3). Whether unpredictability results from a random process or is based on non-linear effects (deterministic chaos) can not be assessed by linear autocorrelation analysis as it is unable to distinguish between these two sources of fluctuation. However, since individual properties, except random individual maximum lifetime, were switched of to reduce the stochasticity of the model it is fare assumption that the unpredictability resulted exclusively from the non-linear effects.

For the highest concentration level the amplitude of the periodic oscillations clearly increased while the frequency decreased. This argues for an increased formation of long-range correlations within the data [93] for high concentration levels, enhancing the long-term predictability. It was thus concluded that the linear predictability of system dynamics increased with increasing disturbance.

As already observed for the raw data phase space, the system was forced towards an asymptotically stable equilibrium point for all DCA concentrations in embedded phase space (see Fig. 4). The fluctuations within the trajectory's settlement area decreased for the overall population and especially for the juvenile daphnids for the highest concentration level, reflecting the population's decreasing ability to respond to changes in environmental conditions (decreasing resilience).

Additionally, as the equilibrium point moved towards the phase space origin with increasing DCA concentration, the risk of extinction also increased.

The generally increasing pattern formation in the recurrence plots revealed, that the system dynamics increasingly tended to return to formerly engaged states when the disturbance increased (see Fig. 5).

While the more uniform spreading of points in the lower concentration levels suggested a high portion of chaotic dynamics (especially for the neonate daphnids), the lineal pattern found for in the highest concentration level suggested a high amount of periodic dynamics within the time series [94]. In contrast, the circular structures observed for juveniles in the highest concentration level suggested the occurrence of quasi-periodic dynamics, like superimposed harmonic oscillations. These structures can be suspected to originate from the realistic and model-built in time delay between juvenile response, their maturation and the according effects on reproduction.

For the juvenile daphnids a mixture-like pattern was observed in the highest concentration level. It seems like the system settled at a transition between quasi-periodic and periodic dynamics. Such a transition has already been observed in a population model for Tribolium [95], in which the varied parameter of disturbance was the rate of larvae developing to adults, resembling the effect on age class dependent mortality of juveniles within the Daphnia model. As already deduced from the autocorrelation function, from our findings it was concluded, that the predictability of the system dynamics increased with increasing disturbance.

The Poincaré sections illustrated clearly, that the attractor determining the system dynamics was more or less stable from the control up to 20 µg/l DCA in the overall population as well as all age classes (see Figs. 6 and S5). The abruptly changing location of the attractor (esp. for juveniles) suggested that the attractor became unstable between 20 µg/l and 40 µg/l.

As all Lyapunov coefficients (5 treatment levels ×3 age classes ×100 simulation runs = 1500 values) were larger than zero, and together with the fact that the dynamics were clearly bounded, there was strong evidence for chaotic dynamics within this model-based population of Daphnia magna. Consulting the Lyapunov coefficient as a measure of rate of chaotic dynamics within the population development, this rate decreased for DCA treatments compared to control (see Fig. 7). From the biological point of view this seems reasonable as positive feedback (influence of exponential growth) is likely to be smaller in the presence of a disturbance affecting reproduction.

The decrease of the coefficients indicated that higher amounts of disturbance decreased the potential for exponential increase of trajectory distances. This was already obvious from the former analyses, but the statistics on the coefficients quantified and validated this. It suggested some sort of concentration-response relationship, thereby approving the Lyapunov coefficient as an indicator to detect non-linear effects of population dynamic disturbance due to environmental stress. Additionally, the significant decrease of the Lyapunov coefficient revealed a true and important emergent effect which could not have been expected from the raw population abundance pattern.

Conclusions

Implications

The decreasing Lyapunov coefficients with increasing DCA concentration levels resembled a concentration-response relationship. This is an important finding as it is directly connected to the important species trait of reproductive capacity that strongly triggers the daphnids' population resilience (see discussion below). Unfortunately it is not straightforward to quantify the effect as the Lyapunov coefficient and the treatment level exhibit a reciprocal relationship. Furthermore, the upper effect bound is difficult to define, because Lyapunov coefficients approaching zero or values below zero mean a discrete change in system dynamics (phase transition towards a limit cycle or a stable equilibrium point respectively), not a steady change as required for a continuous concentration-response relationship. Thus the concentration-response relationship would exhibit discontinuities, difficult to interpret. In addition, Poincaré sections suggested that between 20 µg/l and 40 µg/l the attractor became unstable and a distinct transition to another attractor (basin of attraction) occurred. This would mean an extra discontinuity within the system's reaction to the increased treatment, which is a further clear contradiction to the concept of concentration-response relationships.

Nevertheless, a no-effect threshold value can easily be inferred from the distribution of Lyapunov coefficients. A simple test on significant differences between the mean values of controls and treatment levels (see Tab. 1) can be conducted. In the present study this threshold turned out to be smaller than 2.5 µg/l for all three age classes as well as the overall population (see Tab. 2). Such a threshold could serve as a relevant indicator to assess non-linear effects in population time-series data.

The population dynamics can generally be decomposed into a deterministic (periodic plus chaotic) and a stochastic part, the latter not being considered in this study. The authors are well aware about the diverse discussion about the problem how to distinguish stochastic fluctuations from chaotic dynamics (e.g. [96], [97], [98], [99], [100], [101]). To reduce sensitivity on stochastic fluctuations we used the function lyap k from the Tisean package which is said to suppress statistical fluctuations by excluding reference points which have no more than nfmin (default = 10) neighbours closer than ε [83].

The control level was found to exhibit a larger amount of chaotic dynamics than the treatments. Nevertheless, this does not mean that the controls were not predictable and thus not suitable for assessment. The dynamics produced by the IDAMP-model appeared clearly bounded and represented the naturally occurring fluctuations within an ecologically stable population correctly, esp. under undisturbed conditions [30].

Natural populations as well as theoretical populations from ecosystem models are known to exist in a mixture of deterministic (periodical, ordered) and chaotic (irregular) dynamics [102], [103]. In our study it was shown that the ratio between the two deterministic parts of population dynamics (periodic vs. chaotic) increased with increasing disturbance resulting in the hypothesis that the capability for resilience decreased with increasing disturbance. As can be seen from Figs. 4 and 6 the dynamics are clearly bounded and the population got more and more trapped within a very small phase space region. From these findings it can be supposed that the population's flexibility to react to further external disturbance was reduced.

As can be concluded from the Poincaré maps used in this study, the DCA treatment seemed to destabilize the attractor above 20 µg DCA/l. Consulting the definition of ecological resilience, it can be assumed that this effect was caused by altering the phase space stability landscape, in the way that the system settled within a new basin of attraction. The fact that the system attractor became unstable could mean more than the gradually changing ratio of chaotic and periodic dynamics with increasing disturbance (see autocorrelation function and recurrence plot). It means a discrete, qualitative change in system dynamics (phase transition), which can be attributed to a shift in demographic parameters (like survival or fecundity) and has already been reported before for insect populations [104].

But, generally it is assumed that a press disturbance causes sustained alterations within a system's dynamics. A system under press disturbance is generally not expected to return to an initial configuration [16]. At least for the highest treatment tested it was supposed that there was an increased chance, that the system would not be able to return to the dynamics typical for less disturbed conditions simply by virtue of its own strength when the disturbance was removed. Instead it would need a further ‘disturbance’ to force the system back to the attractor which was observed for up to 20 µg/l due to the altered stability landscape. Recovery in this context would mean a reconfiguration of the stability landscape to a state typical for less disturbed systems. The problem about this definition within the experimental context is that for treatment levels up to 20 µg/l no change of the stability landscape was visible at all, whereas for 40 µg/l a change was found indeed. As a reconfiguration step (switching off the press disturbance) was not included within the experimental design, it was only possible to observe the potential for recovery, not an actual recovery event within the populations. But, from the findings in this study it was concluded that the chance for a change in attractor orientation increased with increasing treatment level.

Outlook

The study revealed that recovery and resilience in age-class structured populations exposed to environmental stress are complex matters and have to be assessed rather deeply to come to relevant conclusions.

Chaos was observed besides periodicity in the model-based age-structured Daphnia population. This was reported before for different models, populations under laboratory conditions and even for multi-species communities under laboratory and semi field conditions (e.g., [38], [95], [104], [105], [106], [107]). Our study connects this fact to the idea of a dynamical definition of recovery. In this study the amount of chaotic dynamics turned out to be the largest in the undisturbed control treatments, pointing out the importance of the natural mechanism of population stabilization via fluctuation in terms of resilience, defined as the flexibility to changes in environmental conditions (dynamical plasticity). The unpredictable part of population dynamics in the control treatments showed that the specification of recovery in terms of the return to the undisturbed control state maybe a questionable definition, as the theory behind the recovery concept is usually a static, not a dynamical one. Using Lyapunov coefficients as indicators it was shown that a disturbance reduces the amount of chaoticity or non-linear dynamics within the population. In our opinion this effect gives strong evidence that anthropogenic stressors can reduce the resilience of age-structured population and foster linear dynamics and resulting periodicity, reflecting equilibrium conditions.

Reduced resilience has a negative influence on populations' ability to respond to further disturbances and makes them more susceptible to extinction due to climatic changes, for example. This is an increasingly important aspect in ecological research which has been reported for many systems and their population, some of them highly endangered (e.g., [108], [109]).

These results stand in contrast to some other studies which instead suggest the increase of chaos and according decrease of stability in populations exposed to environmental stress [110]. But, in this study the constant press disturbance in combination with the discrete reproduction cycle of the daphnids results in a situation comparable to a periodic forcing which has already been shown to result in phase transitions in the dynamics of time discrete non-linear ecosystem models, either from periodic behaviour to chaos [111] or from chaotic dynamics to periodicity, the latter resulting in a chaos control mechanism [112].

The study at hand supports the hypothesis that chaotic oscillations reduce the degree of synchrony among populations or subdivisions of populations (decorrelation effect), thereby stabilizing the system and reducing the probability of extinction (chaos-mediated survival), an interpretation already been published for meta-populations before [113]. Thus, we think that our findings add an important aspect to the open discussion about resilience and chaotic population dynamics in the presence of environmental disturbance, which is still an exciting area of research [114].

The quantitative change in system dynamics and the two constituting components can effectively be visualised by linear methods (autocorrelation analysis) and non-linear techniques (recurrence plot, Poincaré section, Lyapunov exponents). Lyapunov coefficient proved to be a useful indicator to detect disturbance thresholds quantifying non-linear effects of contaminants within the Daphnia population. Even a concentration-response like relationship was observed for Lyapunov coefficients in this study. The modelling of non-linear effects is thus regarded as important by the authors, as non-linearities in species populations cannot be predicted by linear techniques in the long run.

The authors propose their methodology as a way to test for statistically significant non-linear effects of general disturbance on dynamics in long-term population data sets, especially when assessing the risk that a disturbance (like introducing chemicals in ecosystems) persistently alters the dynamics. It allows assessing an important emergent ecological effect on the population level which cannot be retrieved by standard assessment methods, neither based on individual endpoints like filtration rates nor based on population endpoints like growth rates.

The aim of this study was not to propose alternative definitions for the terms stability or recovery, nor to discuss the appearance of non-linear dynamics in population models. To come to a time explicit statement about recovery, like in the approach of engineering resilience, a larger number of simulation studies have to be conducted, so that a mathematical probability can be computed for the fact whether a population dynamic time series returns to the control dynamics after a given time.

Further studies will be done in the future to find a suitable statistical representation for a recovery statement, elucidate the questions about a lasting altered fitness landscape induced by disturbance, calculate mathematical probabilities for this and use finer increments in the disturbance concentration levels to detect a possible point of attractor instability more precisely.

Daphnia abundance from all 100 simulation runs. Solid line: mean abundance, dashed line: 95% confidence interval. Top panel: neonates, mid panel: juveniles, bottom panel: adults.

(TIFF)

Click here for additional data file.

Autocorrelation function of simulation run 1. From top to bottom: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. Top panel: neonates, mid panel: juveniles, bottom panel: adults.

(TIFF)

Click here for additional data file.

Embedded time series of simulation run 1. From top left to bottom right: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. (d = 1, m = 2). Top panel: neonates, mid panel: juveniles, bottom panel: adults.

(TIFF)

Click here for additional data file.

Recurrence plot of simulation run 1. From top left to bottom right, row-wise: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. (d = 1, m = 1,2, ε = 1). Top panel: neonates, mid panel: juveniles, bottom panel: adults).

(TIFF)

Click here for additional data file.

Poincaré sections of simulation run 1. From top left to bottom right, row-wise: control, 2 µg/l, 5 µg/l, 10 µg/l, 20 µg/l, 40 µg/l, resp. (d = 1, m = 3, C = 1). Top panel: neonates, mid panel: juveniles, bottom panel: adults.

(TIFF)

Click here for additional data file.

Non-linear time series analysis using the Tisean package.

(DOC)

Click here for additional data file.