Analysis of simulated long-term ecosystem dynamics using visual recurrence analysis

Ecological Complexity 1 (2004) 111–125

Analysis of simulated long-term ecosystem dynamics using visual recurrence analysis Lael Parrott∗ Complex Systems Laboratory, Department of Geography, University of Montreal, C.P. 6128 Succursale Centre-ville, Montreal, QC, Canada H3C 3J7 Received 26 September 2003; received in revised form 10 January 2004; accepted 21 January 2004 Available online 2 April 2004

Abstract An ecosystem is a prototypical complex system, exhibiting a self-organised temporal dynamics, the characterisation and description of which is riddled with challenges. In this study, simulated ecosystems in the multi-species model WIST (Weather-driven, Individual-based, Spatially explicit, Terrestrial ecosystem model) are used as surrogates with which to explore and characterise the long-term dynamics of ecosystems subject to disturbances of varying scale and frequency. A visual recurrence analysis is used to study the biomass history in each of the simulated systems. The recurrence analysis illustrates clearly the differences between the dynamic regimes of different systems and the subsequent effects of disturbances on these dynamic regimes. Recurrence analysis may prove to be a useful technique for the characterisation of data from ecosystem models and long-term ecological monitoring programs. © 2004 Elsevier B.V. All rights reserved. Keywords: Ecological modelling; Spatially explicit models; Recurrence analysis; Complex systems; Ecosystem management; Complexity

1. Introduction As the impact of human activities on the Earth’s life support systems and more specifically on the structure and functioning of ecosystems becomes more widespread, natural ecosystems are increasingly subject to persistent and frequent levels of disturbance. Detecting and quantifying the effect of such disturbance poses a considerable challenge: how to characterise the dynamics of an ecosystem in such a way that one can conclude whether an observed change is ∗

Tel.: +1-514-343-8032; fax: +1-514-343-8008. E-mail address: [email protected] (L. Parrott).

of major significance (i.e., indicative of system collapse or complete re-organisation) or simply part of the normal envelope of dynamics? While the occurrence of catastrophic shifts in the states of some ecosystems due to a change in a certain environmental parameter is easy to recognise, with many well documented cases (Scheffer et al., 2001), there may well be a less-dramatic, yet perhaps equally significant, internal re-structuring that occurs in ecosystems subject to persistent and frequent levels of disturbance. Such change may not be immediately discernible until the state of the system has been significantly, and often irreversibly, altered. In such circumstances, it is difficult to judge whether the new

1476-945X/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecocom.2004.01.002

112

L. Parrott / Ecological Complexity 1 (2004) 111–125

system state represents an alternative and equally robust stability domain, in which the disturbed system has found a new self-organised state, or whether the resulting system state is one of lower ecological resilience (sensu Holling, 1996) or of reduced ecological integrity that may ultimately lead to system collapse. From an ecosystem management perspective, understanding and characterising the dynamic regime of an ecosystem is, evidently, of utmost importance in order to be able to track and quantify the effects of disturbance on a given system. This task is complicated by the fact that ecosystems are ever changing systems, whose dynamics exhibit all of the key attributes commonly attributed to a complex system, including: non-linearity, feedback between hierarchical levels of organisation (e.g., between an individual and its community), and self-organisation (Levin, 1998; Parrott, 2002). The result is an ordered-yet-disordered dynamic that is difficult to characterise, and even more difficult to predict: like the weather, the dynamics of an ecosystem may be predictable to some degree over the short term, but is largely unpredictable over longer periods. With the addition of disturbance, the effects of which may be amplified by non-linearities in the system response, the task of prediction (and thus, characterisation) becomes even more difficult. There is, for example, considerable uncertainty regarding the relative importance of different deterministic and stochastic driving forces (e.g., density dependence, environmental and demographic stochasticity and climatic variation) and their interactions on observed dynamics (Higgens et al., 1997; Bjørnstad and Grenfell, 2001). In addition, the extent to which chaos is an appropriate model for ecological dynamics remains unclear. Efforts to detect chaos in ecological data sets by estimating Lyapunov exponents or reconstructing the underlying attractor, for example, has not met with incredible success (Hastings et al., 1993). The challenge of characterisation is compounded not only by the complexity of ecosystem dynamics, but also by a general lack of data that cover the range of temporal and spatial scales over which ecosystem processes operate. To date, long-term ecosystem studies in which key variables have been recorded over a period of decades are rare, and measurements are usually not taken at sufficiently frequent intervals to enable anything but a cursory analysis of year to year trends (an exception may be certain long-term river

flow time series). Thus, while many data sets enable an analysis of short term, rapid ecosystem process (e.g., changes in soil biochemistry over a growing season), the data and methods necessary to characterise the system-level dynamics of ecosystems are essentially non-existent. In this study, simulated ecosystems in the multispecies model WIST (Weather-driven, Individualbased, Spatially explicit, Terrestrial ecosystem model) are used to explore and characterise the long-term dynamics of ecosystems subject to disturbances of varying scale and frequency. The ability of a less conventional method of non-linear time series analysis, known as visual recurrence analysis (Eckmann et al., 1987; Casdagli, 1997), to capture key characteristics of the temporal dynamics is investigated. This method, which is little known to ecologists, has been claimed to detect structure in data that has passed all other tests for randomness and is applicable to shorter time series (Grassberger et al., 1991). The visual recurrence analysis proves to be a useful method of characterising and detecting structure in the dynamics of the different simulated ecosystems.

2. Model description WIST stands for weather-driven, individual-based, spatially explicit, terrestrial ecosystem model. The model was originally created to explore the engineering (specifically design and control) of materially closed ecosystems having a structure and composition patterned after natural systems. An earlier description of the model’s dynamics, particularly the emergence of complex spatio-temporal patterns, can be found in Parrott and Kok (2000). Detailed descriptions of the animal and plant sub-models are given in Parrott and Kok (2002) and Parrott and Kok (2001). The basic functioning of the model will therefore be only briefly described here. WIST is “object-based”, in the sense that each component of the modelled ecosystem is represented as a discrete entity whose state is updated at each simulation time step according to a number of deterministic rules. Thus, in accordance with the bottom-up modelling approach (Kawata and Toquenaga, 1994), the global level behaviour of the ecosystem arises as the collective result of many individually simulated en-

L. Parrott / Ecological Complexity 1 (2004) 111–125

tities. Ecosystem components represented by entities include: plant and animal objects, terrain cells, and an atmosphere. The model is largely mechanistic, and is based on a representation of mass flows within and between different ecosystem objects. A complete mass balance is performed during a simulation, such that all mass exchanges, however miniscule, are accounted. The main categories of objects and sources of mass flows in the modelled ecosystem are shown in Fig. 1. Plant objects may represent individual plants or “lumps” of plants in the case of grasses and herbs (covering a maximum area equal to one terrain cell, see below). The growth of each plant object over its life cycle is modelled in a mechanistic fashion in which each plant photosynthesises according to available light and temperature conditions and then allocates photosynthates to its various parts (leaf, stem, root, etc.) according to seasonal requirements. While plant form is not explicitly represented, light is limited according to a simple three-layer canopy model that takes into consideration the leaf area index of plants in the canopy in the calculation of available radiant energy for any given plant. Growth may be limited if there is not enough available water or nutrients to synthesise photosynthates. The plant sub-model can be parameterised to represent plants of different species, including trees, bushes, grasses and herbs. Animal objects represent individual animals (currently only mammalian species have been modelled) that interact with one another and with their environment. Animals must eat in order to meet their metabolic requirements (calculated as a function of body mass, current activity, gestation, etc.). Ingested food is digested into lean or fatty mass and indigestible portions are excreted. A starving animal may metabolise its fatty mass reserves as required to stay alive. Animals may prey upon plants or other animals according to feeding preferences. The community food web is a directed graph of trophic relationships in which each species in the community, as well as detritus, carrion and seeds are represented by nodes. Edges connect two nodes for which there is a predator prey relationship. The strength of each relationship is given a relative value between 0 and 1, which is intended to represent the delectability of a prey item, where a low value means that the predator will eat the prey species very rarely (generally only if nothing else exists to eat) and a value close to 1 represents nearly exclusive pre-

113

dation. In addition to trophic relationships, a number of simple social rules are included. For example, reproduction is currently modelled after mammals, thus young animals stay with and are fed by their mother. Also, prey will flee from predators, and burrowing animals must select a place to build a burrow, and the choice is affected by the presence of other members of their species (allowing for density-dependent rules). The terrain or landscape is modelled as a matrix of cells. Each terrain cell is a rectangular column composed of soil (various layers of decomposing material mixed with inorganic substrate) and water (saturated and unsaturated). The area of the terrain cells is a model parameter and determines the spatial resolution of the model (a resolution of 10 m × 10 m is currently used). While all cells have the same surface area, the height may differ according to landscape topography. During each time step, organic material in each terrain cell is decomposed, rainwater may infiltrate, and subsurface water may flow in or out to neighbouring cells. The atmosphere of the ecosystem is also modelled as a discrete entity. The atmosphere is represented as mixture of water, oxygen, carbon dioxide and nitrogen gases in the same proportions as they are found in the Earth’s atmosphere. Animals and plants exchange gas with the atmosphere due to respiration and photosynthesis. The model ecosystem also contains a number of hypothetical mass storage chambers (“material storage”, Fig. 1) in which about 60% of the total system mass resides. These chambers are intended to mimic the buffering and regulating effects of the enormous material reserves available to the Earth’s natural ecosystems (e.g., in the ocean and the Earth’s crust). Material is moved in and out of these chambers as required to maintain the partial pressure of the atmospheric gases at their set values. These chambers also serve as the source of water for rainfall and as the basin into which surplus surface runoff is collected (to avoid flooding events). The model does not include a climate sub-model. Rather, an artificial weather generator is used to generate a realistic series of temperature, rainfall and radiation values, at a temporal resolution equal to the simulation time step. Alternatively, historic weather data from a real site may be used in place of the simulated values. The various biological activities in the ecosystem (photosynthesis, decomposition, etc.)

114

L. Parrott / Ecological Complexity 1 (2004) 111–125

O CO HO N

ATMOSPHERE

O CO HO

MATERIAL STORAGE

O CO HO

HO BIOLOGICAL COMPONENTS

N

(Living Organisms)

HO (rain)

HO

seeds

whole and part dead organisms

POND

N

TERRAIN

INORGANIC NITROGEN

SEEDS

HO SATURATED AND UNSATURATED WATER

seeds

DECOMPOSING PLANT & ANIMAL MATTER

partially decomposed material

ORGANIC MATTER

N

O CO HO

N

Fig. 1. Main avenues of mass flow in the modelled ecosystem. The biological components category includes all individual plants and animals, whose regular activities of respiration, photosynthesis, litter fall, water and nutrient uptake, food consumption and defecation, and death result in regular mass exchanges with the atmosphere and the terrain (as well as between organisms in the case of predation). The terrain is modelled as a regular grid of 3D cells filled with soil and water. Decomposition of material is a first-order kinetic process that results in the release of gas to the atmosphere as well as the release of nutrients (nitrogen) into the soil. Soil nitrogen is also augmented via the activities of nitrogen fixing microorganisms that capture nitrogen from the atmosphere. The atmosphere, which contains nitrogen, oxygen, carbon dioxide and water, is modelled as a perfectly mixed gas. The partial pressures of each of these components is maintained constant via the controlled exchange of mass with the material storage chambers. Aside from serving to regulate the atmosphere, the material storage chambers serve as a source of rainwater, as well as a reservoir into which water may be pumped to avoid the occurrence of floods. The material storage chambers serve to imitate the role played by the material “buffers” available to natural ecosystems (oceans, atmosphere, Earth’s crust), thereby avoiding catastrophic material shortages in the much smaller modelled system.

L. Parrott / Ecological Complexity 1 (2004) 111–125

are thus functions of the current weather conditions provided by the weather generator or data source. Weather is, therefore, a forcing function of the ecosystem, and there is no feedback between climate and ecosystem-level processes. WIST allows for the parameterisation of any number of different plant and animal species, the initial population sizes of which may be defined for a given simulation. To date, the model has been parameterised to represent a hypothetical temperate climate ecosystem (annual temperature range between 2 and 38 ◦ C) composed of up to about 25 plant and animal species interacting on a gently sloping 500 m × 500 m terrain. The time step used is 10 min, or 144 iterations per simulated day. Simulations based on these model configurations have demonstrated the capability of the model to generate emergent spatial heterogeneity on the terrain (starting from random species distributions), emergent and persistent quasi-cyclic temporal population dynamics, as well as many other non-trivial structures and dynamics that are characteristic of complex systems (Parrott and Kok, 2000).

3. Description of simulations performed A series of simulations were run using two sets of initial conditions, referred to respectively as “prairie” and “herbivore” scenarios. The two scenarios differed only in their initial species compositions. The prairie scenario was based on a set of initial conditions that consisted of only herbaceous plant species and decomposers. The herbivore scenario had the same initial species composition as the prairie scenario, plus three additional herbivore species intended to represent small rodents (“rabbits”, “squirrels” and “voles”). All other initial conditions (initial configuration of the terrain and atmosphere as well as weather parameters) were identical. Simulations were run for a period of 100 years with a time step of 10 min. 3.1. Simulations without disturbance The population dynamics and biomass histories of these simulated systems are shown in Figs. 2(a, b) and 3(a–c). From these figures, it can be observed that both systems undergo a transient period of about 3 years, during which the system re-organises from the pre-

115

scribed starting state (including rather arbitrary initial population sizes and a random spatial distribution of individuals) into a more persistent state, and that this organisation entrains a considerable loss of species. After this period, in the prairie scenario, the two remaining grass species become established over the territory and ultimately reach the system’s carrying capacity, as indicated by the biomass which increases steadily until a relatively constant annual average is obtained. The regular fluctuation of biomass is due to the annual cycle of warm season growth and winter dormancy. In the herbivore scenario, the same two grass species ultimately persist, but their population dynamics is changed markedly by the introduction of predators. The predators have the effect of moderating the amplitude of the annual cycles, as well as introducing a clearer disparity between the two species, due most likely to their differing ability to recover from grazing (given equal availability, the herbivores do not prefer one grass species over another). In addition, a longer term predator–prey cycle with an irregular period of about 20 years emerges spontaneously from the dynamics of this system, and can been seen in the population dynamics of all species. The extinction of the rabbit species at year 24 is most likely due to competitive exclusion by the smaller and more fecund vole species with which it competes directly for food. 3.2. Simulations with disturbance After establishing the baseline behaviour of these two scenarios, a series of additional simulations were run in order to ascertain the response of these two systems to disturbance. For the purposes of this study, a “disturbance” consists of a “mowing” event in which all of the individuals (both plants and animals) in a certain area of the ecosystem are killed by a “mower”, and their biomass is left in place to decompose. A moving event is a singular occurrence. After mowing, the system proceeds as usual and the mowed area is re-colonised as a natural consequence of seed dispersal and migration from other unaffected regions. While mowing does not specifically represent any particular disturbance-type known to natural ecosystems, it is reasonably similar to the effect that a fire, flood, or extremely ravaging parasite might cause to a system. For each scenario, a series of six additional simulations were performed, each with a different

116

L. Parrott / Ecological Complexity 1 (2004) 111–125

Fig. 2. Plant population and total living biomass dynamics for selected simulations based on the “prairie” scenario with disturbance frequencies of 0, 10, 40, and 70 events per 100 years. This is a single trophic-level system in which the resulting system dynamics is dominated by two persistent herbaceous species. Both species (dashed line and solid line in graphs a, c, e and g) are perennial grasses. The number 1 in d corresponds to a key mowing event that is also evident in the recurrence plot of Fig. 4(b).

L. Parrott / Ecological Complexity 1 (2004) 111–125 Fig. 3. Plant and animal population and total living biomass dynamics for selected simulations based on the “herbivore” scenario with disturbance frequencies of 0, 10, 40, and 70 events per 100 years. This is a two trophic level system in which the final dynamics is dominated by a predator–prey relationship between one herbivorous species (“voles”, represented by the dash-dot line in graphs b, e, h, and k) and two herbaceous species. Both herbaceous species (dashed line and solid line in graphs a, c, e and g) are perennial grasses. The two less successful animal species are a “squirrel” and a “rabbit” species (represented by the solid line and the long dash line respectively in graphs b, e, h, and k). Numbers on the plots correspond to key events in the dynamics that are also evident in the recurrence plots (Fig. 5; see text).

117

L. Parrott / Ecological Complexity 1 (2004) 111–125

Fig. 3. (Continued ).

118

L. Parrott / Ecological Complexity 1 (2004) 111–125

probability of occurrence of mowing. The probabilities used were: 10, 25, 40, 55, 70, or 95% per year, where a 70% probability means that, on average, there should be 70 mowing events occurring in a 100-year period. Thus, for each simulated year, it would first be determined whether or not to mow according to the given probability, and then, if a moving were to occur, the location and surface area to be mowed would be selected randomly. According to this scheme, the “size” of the mowing events (in this case area mowed) is distributed in time according to an exponential distribution, whereby small events are very common and large events happen less frequently. The population and biomass dynamics for the two scenarios subjected to increasing disturbance frequencies are shown in Figs. 2(c–h) and 3(d–l). Based on a visual inspection of these figures it can be observed that, while the overall system composition is not altered by the mowing activity (i.e., the species composition remains as before), the dynamics does appear to undergo certain significant changes. First, as might be expected, total living biomass in the system tends to decrease with increased mowing frequency, since the vegetation never manages to re-establish itself between events. For lower disturbance frequencies, however, the systems do seem to recover. For example, in the case of the herbivore scenario with a 10% probability of disturbance, the system demonstrates a typical overshoot before re-establishing previous biomass levels after the large mowing event at year 47 (Fig. 3(f)). Higher frequency disturbance in the herbivore system has the effect of breaking up the long-term predator prey cycles. Lastly, the grass species which is normally less abundant in the undisturbed herbivore system seems to be favoured by the opportunities presented by the disturbance, since it becomes increasingly dominant as the disturbance probability exceeds 70%.

and arising from the majority of biological functions. For each of the simulations, the biomass was recorded once per simulated week, giving a total of 5214 values for each 100-year series. 4.1. Recurrence analysis A recurrence plot is a relatively simple graphical way of discerning recurring structure in a time series that cannot be readily detected using other methods (Casdagli, 1997; Zbilut et al., 1998; Gao and Huaqing, 2000). The method was originally proposed by Eckmann et al. (1987) and has since been applied to the analysis of various theoretical systems in physics as well as to physiological time series (see Weber and Zbilut, 1994 for a good example). The method is of potential interest for the analysis of ecological data since it is robust to the presence of transients and non-stationarity, both of which are typically present in ecological data sets and which limit the applicability of most statistical methods (Legendre and Legendre, 1998). The procedure is briefly summarised here; for more precise details the reader is referred to the original references. A recurrence plot is essentially a graphical representation of the correlation integral of a time series. To construct the plot, the original time series (x1 · · · xt , where t is the total number of observations made at a fixed frequency) must first be embedded in an M-dimensional phase space, according to Takens method of embedding (Hilborn, 2000), whereby each xi is replaced by a vector (vi ) as follows: vi = (xi , xi+d , xi+2d , . . . , xi+(M−1)d )

The overall objective of the analysis was to characterise the temporal dynamics of the two simulated scenarios with and without the introduced disturbances. The total living biomass (kg dry biomass/m2 ) was selected as the variable of study, since any structure detected in the biomass time series should be the result of biological activity occurring at all temporal scales

(1)

where i is the time index, M the embedding dimension, and d the time delay for embedding (d ≥ 1). The result is a series of vectors V = {v1 , v2 , v3 , . . . , vt−(M−1)d }

4. Analysis of temporal system dynamics

119

(2)

that specify the coordinates of consecutive points in an M-dimensional space. Thus, in the simplest case, for M = 2 and d = 1, the t–1 vectors would be: V = {(x1 , x2 ), (x2 , x3 ), (x3 , x4 ), . . . , (xt−1 , xt )). (3) Next, the Euclidean distance between each vector (vi ) and every other vector (vj ) is computed, resulting in a triangular matrix of all distances. A vi × vj recurrence plot is then constructed by placing a dot at point

120

L. Parrott / Ecological Complexity 1 (2004) 111–125

(i, j) if vi is sufficiently close to vj (i.e., vi and vj are within an arbitrary distance ε). Recurrence plots (M = 2, d = 1 week, ε = 0.025 kg) for the prairie and herbivore ecosystem simulations with increasing disturbance frequencies, are shown in Figs. 4 and 5. Interpretation of the patterns that appear in a recurrence plot has been the subject of several studies, notably Weber and Zbilut (1994) and Casdagli (1997). One can imagine the temporal evolution of the time

series as following the diagonal from the lower left corner (where distances between the first few vectors in the series are plotted) to the upper right corner (where distances between the final vectors in the series are plotted). Changes along this diagonal represent changing dynamics in time. Patterns in the recurrence plot are, in general, indicative of structure in the time series, and thus a completely random time series will appear as an evenly scattered mix of dots with no

Fig. 4. Recurrence plots of the prairie scenario for selected disturbance frequencies: (a) 0, (b) 10, (c) 40, and (d) 70 showing the progressive break-up of structure in the time series as the system is increasingly disturbed. The embedding dimension used is 2 and the delay is 1. Values on the x- and y-axes refer respectively to the i- and j-indices of the biomass history. Numbered locations on the plot refer to numbered events in the temporal dynamics (see Fig. 2). Recurrences in the time series appear as dark points. Note that the plots are symmetrical due to the fact that the matrix is triangular.

L. Parrott / Ecological Complexity 1 (2004) 111–125

121

Fig. 5. Recurrence plots of the herbivore scenario for selected disturbance frequencies: (a) 0, (b) 10, (c) 40, and (d) 70 in which cycles (rectangular structures) corresponding to the emergent predator–prey cycle are apparent. The embedding dimension used is 2 and the delay is 1. Values on the x- and y-axes refer, respectively, to the i- and j-indices of the biomass history. Recurrences in the time series appear as dark points. Numbered locations on the plot refer to numbered events in the temporal dynamics (see Fig. 3). Note that the plots are symmetrical due to the fact that the matrix is triangular.

evident pattern. Isolated recurrences (shown as single blackened dots) are due to chance recurrences in the dynamics, whereas groups of recurrences represent existing structure in the time series and, for deterministic systems, represent points in time where the system is revisiting the same region of its attractor. Squares formed by lines parallel and perpendicular to the main diagonal represent cycles. Readers are

referred to the cited literature for examples of recurrence plots for many well known dynamical systems. The recurrence plots shown in Figs. 4 and 5 illustrate a number of immediately striking features: (I) The dynamics of the two simulated ecosystems (prairie and herbivore) can be easily differentiated by comparing Fig. 4(a) with Fig. 5(a), in

122

L. Parrott / Ecological Complexity 1 (2004) 111–125

which quite different recurrence patterns are evident. In Fig. 4(a), the paling of the plot away from the main diagonal during the first 2500 i- and j-indices is indicative of a non-stationary 50-year trend in the data, which becomes less pronounced during the subsequent 50 years. The large number of recurrences (dark points) in the latter part of Fig. 4(a) (i- and j-indices 3000–5000) is due to the relatively stable cycling of the biomass during the final 50 years of the non-disturbed prairie scenario. The white horizontal and vertical bands that criss-cross the dark areas correspond to single years or seasons that are somewhat different from the others, enabling rapid identification of non-average periods in the system’s evolution. In Fig. 5(a), the various oval patterns correspond to different cycles, as described below in (IV). (II) A sequential break-up of observable structure with increasing disturbance frequency can be noted in both Figs. 4 and 5. The recurrence plots for higher disturbance frequencies appear much more random than for the non-disturbed cases, suggesting that there is a fundamental un-structuring of the system’s normal dynamics with increased disturbance frequencies. The dark points in Fig. 4(c) and (d) are probably not indicative of any particular structure in the time series. In these two cases, recurrences occur each time the system repeats a crash-and-recover cycle due to a mowing event. The events are so frequent, however, that the distribution of recurrences in the plots is almost random. The same is true for Fig. 5(c) and (d), although a little more structure seems to be apparent, due perhaps to the added effects of herbivores on the biomass cycles.

(III) The response of the system to mowing events is observable, with the most obvious case being the large mowing event at year 47 for the herbivore disturbance frequency 10 scenario (Figs. 3(f) and 5(b), point #5). Here, a large gap (horizontal and vertical white band with no recurrences) is evident in the recurrence plot after i- and j-indices 2440 (∼47 years×52 weeks), during which time the system is restructuring. This period of restructuring takes about 2 years for the herbivore system, after which the system seems to return to a dynamical regime that is (at least qualitatively) similar to that before the incident. A similar response is observable for the herbivore scenario with 40% disturbance at year 52 (Figs. 3(i) and 5(c), point #6) and for the prairie scenario with 10% disturbance at year 45 (Figs. 2(d) and 4(b), point #1). In this case, the dark rectangle formed by recurrences between i-indices 3000 and 4000, and between j-indices 2000 and 2300 in Fig. 4(b) corresponds to the system’s recovery at approximately year 57. Here, the years 57–76 are shown to be very similar to the period just before the clearcut (years 38–46). The white band that follows (between i-indices 3000 and 4000, and between j-indices 2300 and 3000) corresponds to the recovery period between years 46 and 56. Further along, the dark region of recurrences (between 3000 and 4000 on both axes) is indicative of the regular dynamics that is established between years 57 and 76. This dynamic is subsequently interrupted at year 76 due to a small mowing event (causing a shift in the total biomass values and therefore a break in the repetitive recurrence pattern).

Table 1 Quantitative recurrence analysis for the scenarios shown in Figs. 4 and 5 Simulation Herbivore: disturbance Herbivore: disturbance Herbivore: disturbance Herbivore: disturbance Prairie: disturbance 0 Prairie: disturbance 10 Prairie: disturbance 40 Prairie: disturbance 70

0 10 40 70

See text for explanations.

RR

DET

L

ENTR

0.0179 0.0128 0.0121 0.0169 0.0168 0.0154 0.0079 0.0124

0.9489 0.938 0.9273 0.9413 0.9321 0.9363 0.9324 0.9428

8.7386 8.026 7.1796 7.7955 10.0918 10.023 8.8977 9.246

15.1649 14.7533 14.8006 15.1925 14.9045 14.7918 13.9393 14.5392

L. Parrott / Ecological Complexity 1 (2004) 111–125

123

(IV) A number of other events in the biomass histories have more or less detectable traces in the recurrence plot patterns. The extinction of the rabbit species in the non-disturbed herbivore scenario (Fig. 3(b), point #1), for example, corresponds to a node point on the recurrence plot diagonal (Fig. 5(a)). Also, cycles and other repeating structures in the temporal dynamics appear as rectangles in the recurrence plots. Of particular note are the rectangles marked 2–4 in Fig. 5(a), with corresponding cycles marked in Fig. 3(c). While analysis of the structures in recurrence plots is largely a visual exercise, some authors have proposed ways that these structures may be quantified. Four common measures include: (1) the percentage of recurrence (dark) points (RR); (2) the percentage of determinism (DET), which is the percentage of all recurrence points that form part of line segments parallel to the diagonal (when the length of the line segment exceeds a certain threshold value, usually 2); (3) the average length of these line segments (L); and (4) the Shannon entropy (ENT) of the distribution of the lengths of line segments parallel to the main diagonal. These statistics have been calculated for the plots shown in Figs. 4 and 5 and are given in Table 1. While some authors have reported marked changes in determinism for different dynamical states (e.g., Weber and Zbilut, 1994), in the examples shown here, none of the calculated statistics is particularly revealing. Lastly, the effect of certain recurrence analysis parameters on the plot structure was investigated. In the plots shown in Figs. 4 and 5, an embedding dimension of 2 was chosen arbitrarily for the sake of simplicity, and to emphasise observable structures. In Fig. 6, the effect of changing the embedding dimension (M) for the non-disturbed herbivore system is shown. It is clear that, for this system, the structure deteriorates rapidly as the embedding dimension increases. Choice of a suitable embedding dimension for non-linear analyses has been the subject of many studies and there is no general consensus on how best to choose an appropriate value for empirical data (Hilborn, 2000). Based on Fig. 6, it is clear that for the system studied here, recurrence analysis can be very sensitive to the chosen embedding dimension. In Fig. 7, the effect of chang-

Fig. 6. Effect of changing the embedding dimension on the recurrence analysis of the herbivore scenario without disturbances (Fig. 5(a)). (a) M = 4; (b) M = 8; and (c) M = 16.

ing the time delay (d), equivalent to changing the data sampling frequency, is shown. Here, the structure in the plot remains quite clear, even for d = 8. This is an important outcome, since it suggests that real-world sampling frequencies of once per month (d = 4) or per every 2 months (d = 8) may be sufficient to capture the key elements of the temporal dynamics in ecosystems.

124

L. Parrott / Ecological Complexity 1 (2004) 111–125

Fig. 7. Effect of changing the delay, or sampling frequency, on the recurrence analysis of the herbivore scenario without disturbances (Fig. 5(a)). (a) d = 2, bi-monthly sampling; (b) d = 4, monthly sampling; and (c) d = 8, sampling once every 2 months.

5. Discussion and conclusion For the simulated ecosystems studied here, the recurrence analysis provided an effective way of differentiating between two systems having different dy-

namical regimes; time series from similar ecosystems will show similar patterns in a recurrence plot. In addition, the effect of disturbances on the ecosystem dynamics was immediately discernable via a breakdown of the structure and pattern normally visible in the recurrence plots. Lastly, many of the pattern changes evident in the recurrence plots could be re-traced to known events that occurred in time (example, species losses or mowing events). Thus, even without an a priori knowledge of the history of a system, it may be possible to use a recurrence analysis to detect significant points in its dynamic evolution that correspond to changes in its stability regime, or movements away from a previously established dynamic state. While essentially a visual analysis, the use of recurrence plots to analyse ecological time series may prove to be a useful qualitative method of revealing underlying structure in an ecosystem’s dynamics. Shifting patterns or other changes in a recurrence plot may permit early detection of more subtle change or restructuring that is occurring in an ecosystem. While not a replacement for other methods of analysis, recurrence analysis may lend itself well to hypothesis testing regarding changes in an ecosystem’s dynamics, where a recurrence plot of an historical time series serves as the control or null case to which recurrence plots of current time series may be compared. Such an approach may be particularly useful when other methods of non-linear analysis cannot be applied due to technical limitations such as non-stationarity, cyclicity or short time series lengths. The use of the model WIST made it possible to have sufficiently long time series that were relatively complex and that could be used for comparative analysis purposes. Models play an integral role in the study of complex systems and, while a model cannot replace observations of a real system, it may often serve as surrogate system, particularly when experimental work is impracticable (Schank, 2001). Models have the advantage that they provide complete experimental control over the system of study, and the number of replicates and analyses performed is limited only by available computing power. Models are of particular use in ecology where extensive data collection is an expensive and time-consuming enterprise. Preliminary analysis work with models such as WIST may, therefore, help to identify appropriate methods of characterisation that can be adapted to real situa-

L. Parrott / Ecological Complexity 1 (2004) 111–125

tions. Obviously, the time series studied here, system biomass sampled weekly for a 100-year period, does not exist for any real ecosystem. The method of recurrence analysis is, however, equally applicable to other time series that might realistically be measured during a long-term monitoring program. The method of visual recurrence analysis presented here, while not new, has rarely, if ever, been applied to ecological data. Methods such as these, developed specifically with the objective of describing and characterising the dynamics of complex systems, in conjunction with modelling and simulation studies, may prove to be key tools in elucidating the nature of the inherently complex dynamics of ecological systems. Acknowledgements Recurrence analysis was done in the MATLAB environment (Mathworks Inc.) using version 4.5 of the Cross Recurrence Toolbox created by the Nonlinear Dynamics Group at the University of Potsdam. The author gratefully acknowledges N. Marwan for providing access to this toolbox. Acknowledgements are also due to R. Kok who was involved in the development of earlier versions of the WIST model. A working version of WIST can be obtained by contacting L. Parrott ([email protected]). Financial support for this project was provided by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Canadian Foundation for Innovation.

References Bjørnstad, O., Grenfell, B., 2001. Noisy clockwork: time series analysis of population fluctuations in animals. Science 293, 638–643. Casdagli, M.C., 1997. Recurrence plots revisited. Physica D 108, 12–44.

125

Eckmann, J.P., Oliffson Kamphorst, S., Ruelle, D., 1987. Recurrence plots of dynamical systems. Europhys. Lett. 4 (9), 973–977. Gao, J., Huaqing, C., 2000. On the structures and quantification of recurrence plots. Phys. Lett. A 270, 75–87. Grassberger, P., Schreiber, T., Schaffrath, C., 1991. Nonlinear time sequence analysis. Int. J. Bifurcation Chaos 1 (3), 521–547. Hastings, A., Hom, C., Ellner, S., Turchin, P., Godfray, H.C., 1993. Chaos in ecology: is mother nature a strange attractor? Annu. Rev. Ecol. Syst. 24, 1–33. Higgens, K., Hastings, A., Sarvela, J., Botsford, L., 1997. Stochastic dynamics and deterministic skeletons: population behavior of Dungeness crab. Science 276, 1431–1435. Hilborn, R., 2000. Chaos and Non-Linear Dynamics: An Introduction for Scientists and Engineers. Oxford University Press, London, 650 pp. Holling, C.S., 1996. Engineering resilience versus ecological resilience. In: P. Schulze (Ed.), Engineering Within Ecological Constraints. National Academy, Washington, DC, pp. 31–44. Kawata, M., Toquenaga, Y., 1994. From artificial individuals to global patterns. Trends Ecol. Evol. 9 (11), 417–421. Legendre, P., Legendre, L., 1998. Numerical Ecology, second. Elsevier, Amsterdam, 853 pp. Levin, S., 1998. Ecosystems and the biosphere as complex adaptive systems. Ecosystems 1, 431–436. Parrott, L., 2002. Complexity and the limits of ecological engineering. Trans. ASAE 45 (5), 1697–1702. Parrott, L., Kok, R., 2002. A generic, individual-based approach to modelling higher trophic levels in simulation of terrestrial ecosystems. Ecol. Modell. 154, 151–178. Parrott, L., Kok, R., 2001. A generic primary producer model for use in ecosystem simulation. Ecol. Modell. 139 (1), 75–99. Parrott, L., Kok, R., 2000. Use of an object-based model to represent complex features of ecosystems. In: Proceedings of the Third International Conference on Complex Systems, New England Complex Systems Institute, Nashua, NH, 21–26 May 2000 (InterJournal manuscript #371). Schank, J., 2001. Beyond reductionism: refocusing on the individual with individual-based modelling. Complexity 6 (3), 33–40. Scheffer, M., Carpenter, S., Foley, J., Folkes, C., Walker, B., 2001. Catastrophic shifts in ecosystems. Nature 413, 591–596. Weber, C.L., Zbilut, J.P., 1994. Dynamical assessment of physiological systems and states using recurrence plot strategies. J. Appl. Physiol. 76 (2), 965–973. Zbilut, J., Guiliani, A., Webber, C., 1998. Recurrence quantification analysis and principal components in the detection of short complex signals. Phys. Lett. A 237, 131–135.