CHAOS 20, 043120 共2010兲 Regularity analysis of an individual-based ecosystem simulation Abbas Golestani1 and Robin Gras2 1 School of Computer Science, University of Windsor, Windsor, Ontario N9B 3P4, Canada Department of Biological Sciences, School of Computer Science, University of Windsor, Windsor, Ontario N9B 3P4, Canada 2 共Received 18 June 2010; accepted 19 October 2010; published online 18 November 2010兲 We analyze the results of a large simulation of an evolving ecosystem to evaluate its complexity. In particular, we are interested to know how close to a stochastic or a deterministic behavior our simulation is. Four methods have been used for this analysis: Higuchi fractal dimension, correlation dimension, largest Lyapunov exponent, and P&H method. Besides, we use a surrogate data test to reach a final decision about analysis. As we expect, our results show that there is a deterministic and chaotic behavior in ecosystem simulation. © 2010 American Institute of Physics. 关doi:10.1063/1.3514011兴 We analyze the results of a large simulation of an evolving ecosystem to evaluate its complexity. This simulation is a logical description of how a simple ecosystem performs. Because in reality, biologists do not have much data regarding variation of ecosystem during the time, analyzing logical simulation can help biologists in better understanding of long-term behavior of ecosystem. Therefore, biologists can have access to the valuable data in order to use them as input to their analysis about speciation mechanisms, ecosystem stability, the emergence of a trait, etc. As many scientists believe that natural phenomena have to be considered as deterministic and chaotic systems, it is very important that a simulation used to model such a phenomenon generate nonrandom and complex chaotic patterns. For this purpose, we have used four different methods to evaluate some patterns generated by our simulation. With each of them, we have been able to show clearly that the behavior of our simulation is deterministic and chaotic. I. INTRODUCTION Nonlinear signal processing is an important research area for many applications. Specifications and identifications of nonlinear signals can help us to detect nonlinear behavior of dynamical systems. The discrimination of stochastic and deterministic behaviors of nonlinear time series is a basic topic in nonlinear dynamic fields.1 This specification has attracted researchers for a long time.2 Some studies have shown that when a fractal dimension is a function of the state-space dimension, the deterministic properties of a system can be determined by low values of fractal dimension and the stochastic properties of a system can be determined by high values of fractal dimension.3–5 Among different versions of fractal dimension that are used for complexity estimation, the Higuchi fractal dimension is a precise and applicable mechanism to estimate self-similarity that also gives a stable value for the fractal dimension.6–9 Besides, this method is considered to be the standard for representing the irregularity of time series.10–12 The Higuchi fractal dimension can be calculated directly in the time do1054-1500/2010/20共4兲/043120/13/$30.00 main and is therefore simple and fast. It has also been proved to be able to estimate samples as short as 150–500 data points reliably.13 In the past few years, fractal analysis techniques have gained increasing attention in medical signal and image processing. For example, analyses of electroencephalographic data and other biosignals are among its applications.14 The fractal complexity of the signal in the time domain 共calculated using Higuchi’s algorithm兲 seems to be the simplest method.15–25 The same method may also be used in other biomedical applications.15,26,27 Besides, we also used the correlation dimension using Gaussian kernel algorithm 共GKA兲 method28,29 for the analysis, and the results were similar to Higuchi fractal dimension. We used these criteria to judge whether an ecosystem simulation30 behavior is stochastic or deterministic. The simulation has already presented promising results showing coherent behaviors of the whole simulation with the emergence of strong correlation patterns also observed in existing ecosystems. Therefore, biologists can have access to the valuable data in order to use them as input to their analysis. Obviously, this amount of data is not available for biologists in reality, and it could be a very useful tool to test some hypothesis on several open questions about speciation mechanisms, ecosystem stability, the emergence of a trait, etc. To understand if our simulation is able to generate complex emerging behaviors, we would like to evaluate its level of complexity and predictability. Most of the scientists believe that the chaotic behavior can be observed in many natural systems, such as the weather,31 so natural phenomena have to be considered as a chaotic system.32 Because of the similarity between chaotic and stochastic signals, distinguishing these two types is difficult. In this paper, we try to evaluate the level of complexity of an ecosystem simulation. So, to overcome this problem, we used the largest Lyapunov exponent and an efficient method, which is based on the Poincaré section and the Higuchi fractal dimension, the P&H method.33 So, using four different methods in analyzing could be more reliable to judge about the results. 20, 043120-1 © 2010 American Institute of Physics Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-2 A. Golestani and R. Gras An important point is that a simple chaotic system can produce a time series that passes most tests for randomness. Conversely, a purely random system with a nonuniform power spectrum can masquerade for chaos. Thus, we validate our conclusion about whether an ecosystem simulation behavior is deterministic and chaotic or not by applying the tests to surrogate data designed to mimic the statistical properties of original data using Higuchi fractal dimension, correlation dimension using GKA method, Lyapunov exponent, and P&H method as statistics. A brief overview of our method is as follows: • The Higuchi fractal dimension, correlation dimension, Lyapunov exponent, and P&H method value are extracted from the original time series and its surrogates. • For Higuchi fractal dimension and correlation dimension, random and deterministic distributions are constituted, respectively, as H0 and H1 hypotheses by means of the histograms of surrogate and original data. For Lyapunov exponent and P&H method, nonchaotic and chaotic distributions are constituted, respectively, as H0 and H1 hypotheses by means of the histograms of surrogate and original data. • A desired confidence level is selected regarding the value of standard deviation of the distribution. Typically, the interval out of ␮ ⫾ 4␴ is considered as proper to be the confidence interval. • Finally, the criterion value of the original data is compared with the confidence interval of the distribution. If it lies on the confidence interval, the corresponding time series is considered to have deterministic behavior when the Higuchi’s fractal dimension and correlation dimension are used and have chaotic behavior when Lyapunov exponent and P&H method are used. This paper is organized as follows. In Sec. II, we present our ecosystem simulation. In Sec. III, we explain the different analyzing methods such as Higuchi’s fractal dimension, correlation dimension, Lyapunov exponent, P&H method, and the surrogate data test method. Finally, in Sec. IV, we present the results of applying these methods to the ecosystem simulation data. II. INDIVIDUAL-BASED EVOLVING PREDATOR-PREY ECOSYSTEM SIMULATION USING FUZZY COGNITIVE MAPS AS A BEHAVIOR MODEL In this section, the main parts of the already existing evolving agent-based predator/prey ecosystem are briefly introduced. The comprehensive description of this simulation has been proposed in Ref. 30. This simulation is a logical description of how a simple ecosystem performs. In this simulation, complex adaptive agents 共or simply individuals兲 are either a prey or a predator and a world is implemented as a 1000⫻ 1000 matrix of cells. Our global aim was to develop a generic platform that was able to simulate complex ecosystems with “intelligent” agents interacting and evolving in a large and dynamic environment. For that, we have chosen to use an individual-based model and to implement the predator-prey paradigm. In a run Chaos 20, 043120 共2010兲 of the simulation, time is divided into discrete steps 共as in a digital clock兲 called time steps, composed of a set of common phases 共see Sec. II B 3兲. Each agent’s behavior is modeled by a fuzzy cognitive map.34 This enables an agent to evaluate 共with memory兲 its environment 关for example, its distance to food or to a possible mate, or its energy level兲 and its internal state 共for example, fear, hunger, or curiosity兲, and to choose among several possible actions 共for example, escaping, eating, or breeding兴. Furthermore, agents’ behaviors are allowed to evolve 共at the phylogenetic scale兲 over the time steps of the simulation.35 A. General concepts founding our simulation 1. Individual-based modeling In the area of ecosystem simulation, individual-based modeling is a bottom-up approach allowing for the consideration of the traits and behavior of individual organisms. Instead of modeling an entire ecosystem as a whole, individual-based models aim to “treat individuals as unique and discrete entities.”36 By modeling organisms with varying characteristics 共such as age, mating preferences, and role in the ecosystem兲, the properties of the system that the individuals represent can begin to emerge from their interactions. While the use of individual behavior has been included in many models during recent decades, the individual-based modeling approach is increasing as the cost to purchase and operate a machine capable of running time-consuming simulations reduces. This has brought interesting benefits, for example, in forest ecology, fish recruitment modeling, and spatial heterogeneity depiction.37 Yet, a few attempts have been made to simulate a complete and complex ecosystem. An example of such a system is the platform Echo, which also includes an evolutionary mechanism. However, the organisms in Echo are very simple and are not provided with any behavior model. The system closed to what we are interested to simulate is the Avida simulation.38 It has nevertheless some limitations: the individuals do not move and are quite limited in number, there is no explicit representation of species, and more importantly, there is a fix fitness function which means that the system is in fact an optimization process and therefore cannot be considered as a complex system. 2. Predator-prey models An important property we wanted to integrate into our simulation was the fact that agents would have to develop efficient behaviors to be able to survive in their environment. We have therefore chosen a predator-prey model in which behaviors of prey and predators have to evolve simultaneously to give them abilities to survive. A direct consequence of this choice is that our ecosystem encompasses individuals belonging to two trophic levels. Some other predator-prey models have already been proposed, such as the one in Ref. 39, for example. Yet, this particular agent model is dedicated to represent schooling behaviors, and the evolution is an offline mechanism using a genetic algorithm.35 Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-3 Regularity analysis of a simulation Chaos 20, 043120 共2010兲 3. Fuzzy cognitive maps In general, fuzzy cognitive maps 共FCMs兲 are a weight graph aiming to represent the causal relationship between concepts and to analyze inference patterns 共the final states of the system after convergence兲. In our simulation, the FCM is not only the base for describing and computing the agent behaviors, but also the platform for modeling the evolutionary mechanism and the speciation events. The individuals perform one action each during a time step based on their perception of the environment. The FCM, called a map in our system, is used to model the agent behaviors 共structure of the graph兲 and to compute the next action of the agent 共dynamics of the map兲. Formally, a FCM is a graph which contains a set of nodes C, each node Ci being a concept, and a set of edges I, each edge Iij representing the influence of the concept Ci on the concept C j. A positive weight associated with the edge Iij corresponds to an excitation of the concept C j from the concept Ci, whereas a negative weight is related to an inhibition 共a null value which means there is no influence of Ci on C j兲. An activation level ai is associated to each concept. A FCM allows to compute the new activation levels of the concepts of an agent based on its perception and on the current activation levels of its concepts. This computation is called the dynamics of the map and is a normalized matrix product. A map contains three kinds of concepts: sensitive, internal, and motor. The activation level of a sensitive concept is computed by a fuzzification of the information coming from the environment 共see Fig. 1兲. The activation level of the motor concept is used to determine what the next action of the agent will be, and a defuzzification of its value can be used to determine the amplitude of the action. Finally, the internal concept’s activation levels correspond to the levels of intensity of the internal states of the agent and affect the computation of the dynamics of the map. B. Description of our evolving ecosystem simulation In our simulation, each agent possesses its own genome. Agents can mate with genetically similar individuals and produce offspring that inherit a modified combination of the genomes of their parents. We have also implemented a speciation mechanism based on the idea of genotypic pool. Species can emerge from the evolution of individuals or get extinct if all their members die. 1. Agents Each agent possesses its own FCM that models its behavior. This FCM contains sensitive concepts such as predator_close 共prey only兲, prey_close 共predator only兲, food_close, mate_close, and energy_low; internal concepts such as hunting 共predator only兲, fear, hunger, sexual_need, curiosity, and satisfaction; and motor concepts such as escape 共prey only兲, search_for_prey 共predator only兲, search_for_food, socialize, eat, and breed. It includes also links and weights representing the mutual influences of these concepts. It is important to notice that the activation level of each motor concept depends on a complex and nonlinear combination of excitatory and inhibitory influences from both sensitive and internal concepts. Therefore, this behavior model enables the FIG. 1. 共Color online兲 A simple fuzzy cognitive map for detection of foe and decision to evade with its corresponding matrix L and 0 for “foe close,” 1 for “foe far,” 2 for “fear,” and 3 for “evasion,” and the fuzzification and defuzzification functions. representation of very complex phenomena. It is also worth noting that the activation levels of an agent’s concepts are never reset during its life. As the computation of the activation level of a concept involves its previous value, the evaluation of the current state of an agent incorporates all its past states. It means that an agent has a memory of its own past that will influence its future states. Eventually, it appears that an agent’s behavior dynamically depends on a complex combination of the information it currently receives from its environment, its current internal state, and all the past states it went through during its life. The FCM of an agent 共or more precisely the set of edges of its FCM as well as their weights兲 represents also its genome. It is transmitted to the agent’s offspring after combination with one of the other parent and after possible mutations. Links between concepts can appear or disappear during this process, so that the structure and complexity of the maps can change during the evolutionary process. Each agent also possesses several physical characteristics such as its maximum ages, its minimum age for mating, its maximum speeds, its vision distance, and its levels of energy. Energy is provided to individuals by the resources 共grass or meat兲 they find in their environment. An agent consumes some energy each time it performs an action, proportionally to the complexity 共number of edges兲 of its FCM. If it uses all its energy, an individual dies. In our simulation, each individual possesses its proper map which contains around 30 concepts and hundreds of edges. 2. Species Our simulation implements a speciation mechanism directly related to the genotypic cluster definition proposed in Ref. 40. Our mechanism accounts for the gradualism and the fuzziness of the speciation process. In this simulation, a species is a set of individuals associated with the average of the genetic characteristics of its members. The average map is computed based on the FCM matrices of all individuals that Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-4 A. Golestani and R. Gras FIG. 2. 共Color online兲 Population of prey and predator agents. Chaos 20, 043120 共2010兲 FIG. 3. 共Color online兲 Evolution of the prey and predator species. 4. Evolution are members of a species. It is considered that an individual belongs to a species if the difference between its matrix and the average matrix of the species is below a speciation threshold; the threshold is the same for all species.41 Our speciation method begins by finding the individual in a species S with the greatest distance from the species. If this distance is greater than a predefined threshold for speciation, a two-means clustering is performed to split the initial species into two species. Otherwise, species S remains unchanged. 3. World Our simulation takes place in a virtual toric world composed of 1000 cells in both dimensions. Each cell can contain an unlimited number of agents. The initial numbers of prey and predators in the world are parameters of the simulation. Each cell can also contain resources, providing agents with energy: some grass, eatable only by prey, and some meat 共dead prey actually兲 edible only by predators. In order to create a competition for resources, the amount of resources available in a cell has been limited 共by a parameter of the simulation兲. We have modeled some growth and diffusion mechanisms in order to create a realistic evolution of the dispersion and of the amount of grass. Predators have two modes of nutrition: hunting and scavenging. When a predator kills a prey, new meat units are added in the corresponding cell, which can be eaten by the same or other predators. At each time step, the values of the states of all the parameters in the model are updated. The successive phases of the update process for each agent are as follows: perception of the environment, computation of all concepts, application of their action, and update of the energy level. Then there is an update of the lists of agents, species, and cells of the world. An agent has a quite short lifespan 共in terms of the number of time steps兲 and performs only a few dozens of actions during its life. This enables us to obtain a high level of population renewal, which is an important criterion for studying an evolutionary process. In our simulation, evolution stems from several mechanisms: mating, mutation, and speciation. Since species membership of the agents is evaluated at each time step, births and deaths of individuals influence the general species composition.35 Thus, a species can emerge or disappear at any time step. This enables us to model the evolution of populations of individuals sharing important genetic properties. Due to our species model, species evolution is derived directly from individual evolution. This one can occur when a breeding event happens. If the mating is successful, the two parents give birth to a unique offspring. This one inherits a combination of the genomic information of its parents with possible mutations. Moreover, some new edges can be created 共based on a probability of apparition which is a parameter of the simulation兲 and some old edges can be removed. The apparition of new edges is a very important mechanism in the sense that new influences between concepts can emerge during the evolutionary process. This allows the apparition of more complex and potentially more adaptive behaviors. If they show a selective advantage, such behaviors will be preserved 共and thus transmitted through generations兲 by the process of natural selection, inherent to the interaction of the individuals with their environments. As a counterpart, the possibility for edges to disappear is also fundamental. When the complexity 共i.e., the number of edges兲 of the FCM grows, this increases the energy needs of the agent, which then needs a more efficient behavioral model for being able to obtain this energy. Thus, the influence links between concepts are somehow “tested” by the evolutionary process and are removed if they appear to be not beneficial enough. This allows agents 共at the phylogenetic scale兲 to react to changes in the environment and to balance the interest of a complex behavioral model with its energy cost. A unique breeding event generating a mutated offspring cannot produce a very different behavior model. It is the accumulation of neutral mutations during several generations that allows the apparition of new individual behaviors and then of new species. Figure 2 shows the population of prey and predator agents after each time step 共generation兲. Similarly, Fig. 3 depicts the evolution of the number of prey and predator species in successive time steps. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-5 Regularity analysis of a simulation Chaos 20, 043120 共2010兲 A recent execution of the simulation produced approximately 40 000 generations 共time steps兲 in 36 days. In the final generation, there are 190 000 preys, 22 500 predators, 49 prey species, and 19 predator species, which are values quite stable during the whole execution of the simulation. We want to understand how complex and predictable our simulation is. Dynamic studies of nonlinear systems allow us to describe the specification of biological processes.42 Calculation tools of chaos theory such as fractal dimension or signal presentation in state space are suitable for these processes.13 III. ANALYSIS METHODS A. Higuchi fractal dimension method In Higuchi’s method,6 which is used for fractal dimenk sion calculation, we must construct new time series xm from the input time series x共1兲 , x共2兲 , . . . , x共N兲 as follows: 再 冉 ⌊ ⌋ 冊冎 k xm = x共m兲,x共m + k兲,x共m + 2k兲, . . . ,x m + N−m ·k k , 共1兲 m = 1,2, . . . ,k, where both m and k are integers and ⌊ ⌋ is Gauss’s notation. m is the initial time and k is the interval time. For a time interval k, we get k sets of new time series 共for example, if k = 3 there are three time series: x31, x32, and x33兲. The length of the curve Lm共k兲 is defined as follows: Lm共k兲 = 兺i=1兩x共m + ik兲 − x共m + 共i − 1兲k兲兩共N − 1兲 . N−m k k ⌊ ⌋ 共2兲 N is the number of samples and N − 1 / ⌊N − m / k⌋k is the normalization factor. The length of the curve for the time interval k, L共k兲, is defined as the average value over k sets of Lm共k兲. If L共k兲 ⬀ k−D, then the curve is a fractal with dimension D. Time series can be identified by low values of Higuchi fractal dimension and the stochastic properties of a time series can be determined by high values of Higuchi fractal dimension. In current implementation, the range is between 0 and 1 in order that values close to 0 specify deterministic time series and values close to 1 specify random time series. B. Correlation dimension Correlation dimension is an extension of the usual notion of dimension to objects with a fractional dimension.43 In dimensions one, two, three, or more, it is easily established, and intuitively obvious, that a measure of volume V共␧兲 共e.g., length, area, volume, and hypervolume兲 varies as V共␧兲 ⬀ ␧d, where ␧ is a length scale 共e.g., the length of a cube’s side or the radius of a sphere兲 and d is the dimension of the object. For a general fractal, it is natural to assume a relation like V共␧兲 ⬀ ␧d holds true, in which case its dimension is given by d⬀ log V共␧兲 , log ␧ and therefore 共3兲 log V共␧兲 . ␧→0 log ␧ d ⬇ lim 共4兲 Let 兵zt其 be an embedding of a time series in Rd. We therefore define the correlation function CN共␧兲 by CN共␧兲 = 冉冊 N 2 −1 兺 I共储zi − z j储 ⬍ ␧兲. 共5兲 0ⱕi⬍jⱕN Here, I共X兲 is the indicator function, which has a value of 1 if condition X is satisfied and 0 otherwise, and 储 · 储 is the usual distance function in Rd. N is the total number of points. The correlation dimension d is then defined as the slope of the line log共CN共␧兲兲 versus log共␧兲 at small scales where ␧ → 0, CN共␧兲 ⬀ ␧d , 共6兲 log CN共␧兲 . d = lim lim ␧→0 N→⬁ log ␧ Although correlation dimension has been introduced as above, for reliable estimation of correlation dimension, we compute correlation dimension using GKA method.44 Gaussian kernel algorithm takes an entirely different approach to modeling the noise in a signal. We return to the original definition of correlation dimension, i.e., Eqs. 共5兲 and 共6兲. Since the observations xn are contaminated by noise, one cannot know zn precisely.43 Therefore, the computation of I共兩zi − z j兩 ⬍ e兲 in Eq. 共5兲 is actually somewhat fuzzy. Rather than adopting contemporary density estimation to improve the estimate of the distribution, one can model this uncertainty by replacing the hard indicator function I共 · 兲 with a continuous one. The choice 共implied by its title兲 of the Gaussian kernel algorithm is the Gaussian basis function exp− 储zi − z j储2 / 4␧. Details of this algorithm are described in Ref. 45 and an efficient implementation of this technique is presented by Yu.44 Correlation dimension estimates from the GKA suggest the complexity of the system attractor 共the number of active degrees of freedom兲.43 The correlation dimension gives an estimate of system complexity. The GKA separates the data into purely deterministic and stochastic components.28 C. P&H method This method is a new and efficient method, which is based on the Poincaré section and the Higuchi fractal dimension, and it could detect the random signals from deterministic signals.33 Besides, this method can be used to detect chaotic behavior in a signal. Recently, this method has been applied to some biomedical data.42 P&H method is composed of several steps that can be summarized as follows. The first step is to intersect the time series trajectory with the Poincaré section. This intersection induces a set of points that indicate dynamic flow. This intersection leads to a series 共P兲 specified by these points. Applying the Higuchi fractal dimension to the P series yields to a vector L共k兲. This vector is a basis for decision making about series specification. Figure 4 gives an example of a Poincaré section for a 3-dimensional flow 共⌫兲. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-6 A. Golestani and R. Gras Chaos 20, 043120 共2010兲 FIG. 4. Intersection between the flow 共⌫兲 and the Poincaré section 共S兲 generating the set of points P = 共P0 , P1 , P2兲. P0, P1, P2, and so on come from intersections of the flow with the Poincaré section 共see Fig. 5兲. This intersection leads to a set of points P = 共P0 , P1 , P2 , . . .兲 that indicate the dynamic flow behavior.46–48 Many of the time series features, such as periodicity and quasiperiodicity of the original curve 共⌫兲, could still be extracted from the set of points. One-dimensional signals must be embedded according to the embedding theorem.49 Then the calculating procedure of Higuchi fractal dimension is applied to S series except the last level, 再 冉 ⌊ ⌋ 冊冎 k Pm = P共m兲, P共m + k兲, P共m + 2k兲, . . . , P m + N−m ·k k 共7兲 m = 1,2, . . . ,k. For each k , Sm Lm共k兲 = , Lm共k兲 is equal to 兺i=1兩P共m + ik兲 − P共m + 共i − 1兲k兲兩共N − 1兲 . N−m k k ⌊ ⌋ 共8兲 N is the number of samples and N − 1 / ⌊N − m / k⌋k is the normalization factor. Thus, Lm共k兲 is obtained that indicates series of S specification. Now we consider the figure with log共Lm共k兲兲 as the vertical axis and k as the horizontal axis. This has been used to propose the typical criteria for detecting stochastic signals from deterministic signals. For random time series, this value decreases according to Fig. 6共b兲. As shown in Sec. III A, Lm共k兲 is obtained by summing approximately 共N / k兲 terms. If we ignore the normalization factor, then as k increases, the number of summed k time series were random series, terms decreases. If the xm then the positive subtraction of consequent terms yields a value that is also random. So in the value of Lm共k兲, 共N / k兲 random terms are summed; in the value of Lm共k + 1兲, N / k + 1 random terms are summed. Therefore, 共N/k兲X → Lm共k兲, 共9兲 共N/k + 1兲X → Lm共k + 1兲. As k increases, the number of summands decreases. Therefore, the value of Lm共k兲 decreases. Based on this specification, this irregularity criterion has been proposed to distinguish between random time series and deterministic time FIG. 5. 共Color online兲 Intersection of time series trajectory and Poincaré section. series.33 For random time series, a decreasing pattern 关Fig. 6共b兲兴 is observed; for chaotic time series, a nondecreasing pattern 关such as in Fig. 6共a兲兴 is observed. For chaotic time series, it has been shown that because of the stretching and folding property,50 there is at least one k value by which the value of Lm共k + 1兲 is greater than Lm共k兲 and the Lm共k兲 vector indicates a zigzag pattern.33 According to this property, L共k兲 is greater for any odd values of k than for the previous even values of k. As it is apparent in Fig. 6, the P&H method generates a decreasing monotonous pattern for all stochastic signals. For chaotic time series, because of the stretching and folding property, the L共k兲 vector indicates a zigzag pattern. There are different chaotic signals based on route to chaos 共perioddoubling or quasiperiodic or intermittent兲.50 This method has been tested to examples from all types of route to chaotic signals.33 The P&H method output can be easily mapped to negative and positive numbers for stochastic and chaotic time series, respectively. While other criteria do not have the ability of identifying discrete maps such as Henon map,51 P&H method shows satisfying results.3,33 D. Lyapunov exponent Most experts would agree that chaos is the aperiodic, long-term behavior of a bounded, deterministic system that demonstrates sensitive dependence on initial conditions. For that purpose, we must quantify the sensitivity.50 Lyapunov exponents quantify the exponential divergence of initially close state-space trajectories and estimate the amount of chaos in a system.52 A bounded dynamical system with a positive largest Lyapunov exponent is chaotic.50 Suppose a one-dimensional function Xn+1 = f共Xn兲. Imagine two nearby initial points at X0 and X0 + ⌬X0, respectively. After one application of the function 共f共Xn兲兲, the points are separated by Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-7 Regularity analysis of a simulation Chaos 20, 043120 共2010兲 FIG. 6. 共Color online兲 Applying of P&H method over Lorenz time series and random time series. ⌬X1 = f共X0 + ⌬X0兲 − f共X0兲 ⬵ ⌬X0 f ⬘共X0兲, where f ⬘ = df / dX. Now define the local Lyapunov exponent ␭ at X0 such that e␭ = 兩⌬X1 / ⌬X0兩 or ␭ = ln兩⌬X1/⌬X0兩 = ln兩f ⬘共X0兲兩. 共11兲 To obtain the largest Lyapunov exponent, an average of the above equation is taken over many iterations, N−1 1 ␭ = lim 兺 ln兩f ⬘共Xn兲兩. N N→⬁ n=0 N 共10兲 共12兲 The largest Lyapunov exponent determines the average exponential rate of separation of two nearby initial conditions or the average stretching of the space. A positive value shows chaos.50 The different methods that have been proposed for computing Lyapunov exponents from time series can be divided in two classes: Jacobian-based methods and direct methods. Direct methods directly estimate the divergent motion of the reconstructed states without fitting a model to the data.53 The method, which has been used in this paper, was proposed by Sato et al.54 and Kurths and Herzel.55 The average exponential growth of the distance of neighboring orbits is studied on a logarithmic scale via the prediction error, p共k兲 = 冉 冊 1 储X共n + k兲 − X共nn + k兲储 兺 log2 储X共n兲 − X共nn兲储 , Nts n=1 共13兲 where X共nn兲 is the nearest neighbor of X共n兲 共nn is a neighbor index of n which could be n − 1 , n + 1 , . . .兲. The dependence of the prediction error p共k兲 on the number of time steps k may be divided into three phases.53 Phase I is the transient where the neighboring orbit converges to the direction corresponding to the largest Lyapunov exponent. During phase II, the distance grows exponentially until it exceeds the range of validity of the linear approximation of the flow. Then phase III begins where the distance increases slower than exponentially until it decreases again because of folding in the state space. In phase II, a linear segment with slope ␭1 appears in the p共k兲 versus k diagram. This allows an estimation of the largest Lyapunov exponent ␭1.56 Figure 7 gives an example to determine the largest Lyapunov exponent ␭1 of data by this method.53 With some modifications, there are several implementations of Lyapunov exponent in case that time series is nonstationary.57 E. Surrogate data test method A simple chaotic system may produce an output time series, such that it passes randomness testing. In contrast, a Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-8 A. Golestani and R. Gras Chaos 20, 043120 共2010兲 N/2 Xn ⬵ 2␲mn 2␲mn a0 + bm sin , + 兺 am cos 2 m=1 N N am = 2␲mn 2 , Xn cos 兺 N n=1 N N 共14兲 N bm = 2 2␲mn 兺 Xn sin N , N n=1 N/2 S m = a m2 + b m2 ⇒ Y n = FIG. 7. Prediction error p for experimental data vs the number of time steps k. The slope of the solid line in the intermediate range of k gives the largest Lyapunov exponent ␭1 = 0.16. completely random process with nonuniform power spectrum 共correlated noise兲 may reshape to a chaotic system. To solve the problem of separating such systems from each other, time series are tested by means of surrogate data. It is rather simple to generate a time series with similar probability distribution as the original data. It is enough to make random sampling in time series. Randomizing samples preserve the probability distribution, while it does not preserve the power spectrum and autocorrelation function of the time series. Even if the original data are not uncorrelated and white, surrogate data will have these characteristics. To constitute surrogate data which have the same power spectrum of Xn, a surrogate data time series Y n with the same Fourier coefficients but with random phase is generated, 冉 冊 mn a0 + rm , + 兺 冑Sm sin 2␲ 2 m=1 N where rm is a random number between 0 and 1 共0 ⱕ rm ⱕ 1兲. In this method, the power spectrum does not change, but probability distribution becomes a Gaussian distribution.58 A simple solution for getting the same probability distribution is to sort original and random phase data samples separately from smallest to largest. Then the smallest value of surrogate data is assigned to the smallest value of original data, and so on, until the end of the samples. Finally, we will have a random time series with the identical distribution function of original time series and with only slightly altered power spectrum.3,53 IV. RESULT This part is an examination of simulation’s population data and it will be shown whether these time series have a stochastic behavior or followed a specific order. For a better understanding, the methods have been applied to the simulation’s population data, random time series, and Lorenz time series, which is one of the well-known chaotic time series. The Higuchi fractal dimension and correlation dimension are used to show the deterministic behavior of the simulation data, and the largest Lyapunov exponent and P&H method are used to show the chaotic behavior of simulation data. FIG. 8. 共Color online兲 The results of hypothesis testing using 24 surrogate data sets for random time series 共left兲 and Lorenz time series 共right兲 using Higuchi fractal dimension. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-9 Regularity analysis of a simulation FIG. 9. 共Color online兲 The results of hypothesis testing using 24 surrogate data sets over simulation’s population time series; 共a兲 prey and 共b兲 predator using Higuchi fractal dimension. As we stated in Sec. III E, in surrogate data test method, several time series are generated randomly from original data, but with the same characteristics of original data, such as power spectrum and probability distribution. Then an appropriate statistic 共Higuchi fractal dimension, correlation dimension, largest Lyapunov exponent, and P&H method兲 is extracted from both the original and surrogate data. Then we should decide if the difference between the values of original data and surrogate data is statistically meaningful or not. For this purpose, one can generate several surrogate time series with different random phases. Then the mean and the standard deviation for surrogate data sets are computed. Finally, the value of original data is compared to the mean to see if it is significantly far from the mean considering the value of standard deviation. Figures 8–13 show the results of hypothesis testing on several time series by using Theiler algorithm Chaos 20, 043120 共2010兲 FIG. 10. 共Color online兲 The results of hypothesis testing using 24 surrogate data sets over simulation’s population time series; 共a兲 prey and 共b兲 predator series using correlation dimension II for generating surrogate data sets. The results of hypothesis testing have been obtained by using 24 surrogate data sets. A. Higuchi fractal dimension To determine if the Higuchi fractal dimension of a typical time series presents a stochastic or deterministic behavior, a hypothesis testing procedure is adopted. Several surrogate data are generated from the original time series, but with different random phases. These surrogate signals have a completely different structure in comparison to the original data, but similar histogram and power spectrum. The larger the number of generated surrogates is, the more accurate randomness distributions for hypothesis testing are. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-10 A. Golestani and R. Gras Chaos 20, 043120 共2010兲 FIG. 11. 共Color online兲 The results of hypothesis testing using 24 surrogate data sets for random time series 共left兲 and Lorenz time series 共right兲 using largest Lyapunov exponent. The hypothesis under H0 models randomness and the hypothesis H1 implies determinism. These two hypotheses are specified in more details. H0: time series has a stochastic behavior. In other words, the time series is completely irregular without any special structure. H1: time series corresponds to a deterministic system. In other words, the extracted statistic represents a more regular structure in comparison to the pure random surrogates. Surrogate data methods are used to construct a probability distribution function under H0. As we stated before, in surrogate data methods, several time series are generated randomly from original data, but with the same characteristics of original data, such as power spectrum and probability distribution. Then a Higuchi fractal dimension statistic is extracted from both the original and surrogate data. Then, it should be decided if the difference between the values of original data and surrogate data is statistically meaningful or not. For this purpose, one can generate several surrogate time series with different random phases. Then the mean and the standard deviation for surrogate data sets are computed. Finally, a confidence level for acceptance or rejection of the hypothesis H0 can be determined. Typically, the intervals out of ␮ ⫾ 4␴ are proper to be defined as a confidence interval. If the value of the statistic lies on the confidence level interval, the H0 hypothesis is rejected; otherwise, H0 is accepted. To show the clear distinction between the results obtained with stochastic and deterministic time series, we apply this method to random time series, Lorenz time series, and simulation’s population data. Figure 8 shows the results of hypothesis testing on several well-known systems. The red sign indicates the value of the statistic extracted from the original data and the blue one indicates the value of the statistic extracted from the surrogate data. As it can be seen in Fig. 8 共left兲, the hypothesis testing result is in the central range of the distribution for random time series. Therefore, hypothesis H0 is accepted which means that the time series is random. On the other hand, for Lorenz time series, Fig. 8 共right兲, the red signal indicates the value of the statistic ex- tracted from the original data, which is completely separable from its imposter random distribution. It should be noticed that due to the change of abscise scale compared to Fig. 8 共left兲, the surrogate distribution looks like a line, but it still corresponds to the distribution Higuchi fractal dimension value for 24 surrogate time series. Therefore, hypothesis H0 is rejected which means that the time series is deterministic. Having observed the result of this test for well-known time series, the comparison with the results obtained with the simulation time series leads to a clear interpretation. According to Fig. 9, the red signals indicate the value of the statistic extracted from the original data both for prey and predator time series. These values are completely separable from the surrogate distribution. It means that the red signals are in the confidence interval, so the null hypothesis is rejected 共random behavior rejected兲 and the deterministic behavior of simulation’s population data can be concluded. This criterion, therefore, shows that the simulation’s population data have a deterministic behavior. B. Correlation dimension using GKA method To determine if the correlation dimension 共using GKA兲 of a typical time series presents a deterministic or stochastic behavior, the same hypothesis testing procedure is adopted. As in the previous section 共see Sec. IV A兲, the hypothesis under H0 models a stochastic behavior and the hypothesis H1 implies a deterministic behavior. The correlation dimension statistic is also extracted from both the original and surrogate data. It should be decided if the difference between the values of original data and surrogate data is statistically meaningful or not. For this purpose, the same process as in the previous section 共see Sec. IV A兲 is applied to the random time series, Lorenz time series, and simulation’s population data. The same results as in the previous section are observed for the random and Lorenz time series 共data not shown兲. In Fig. 10, the red signals indicate the value of the statistic extracted from the original data for both prey and predator time series for GKA embedding dimension m = 5. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-11 Regularity analysis of a simulation FIG. 12. 共Color online兲 The results of hypothesis testing using 24 surrogate data sets over simulation’s population time series; 共a兲 prey and 共b兲 predator series using largest Lyapunov exponent. These values are completely separable from the surrogate distribution. It means that the red signals are in the confidence interval, so the null hypothesis is rejected 共stochastic behavior rejected兲. This criterion, therefore, also confirms that the simulation’s population data have a deterministic behavior. C. Lyapunov exponent To determine if the largest Lyapunov exponent of a typical time series presents a nonchaotic or chaotic behavior, the same hypothesis testing procedure is adopted. The hypothesis under H0 models a nonchaotic behavior and the hypothesis H1 implies a chaotic behavior. The Lyapunov exponent statistic is also extracted from both the original and surrogate data. It should be decided if the difference between the values of original data and surro- Chaos 20, 043120 共2010兲 FIG. 13. 共Color online兲 The results of hypothesis testing using 24 surrogate data sets over simulation’s population time series; 共a兲 prey and 共b兲 predator using P&H method. gate data is statistically meaningful or not. For this purpose, the same process as in the previous section 共see Sec. IV A and IV B兲 is applied to the random time series, Lorenz time series, and simulation’s population data. Figure 11 shows the results obtained as reference for the random and Lorenz time series. As it can be seen in Fig. 11 共left兲, the hypothesis testing result is falling into imposter distribution for random time series. Therefore, hypothesis H0 is accepted and it means that the time series is nonchaotic. On the other hand, for Lorenz time series, Fig. 11 共right兲, the red sign indicates the value of the statistic extracted from the original data, which is completely separable from its imposter distribution. Therefore, hypothesis H0 is rejected and it means that the time series is chaotic. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-12 A. Golestani and R. Gras In Fig. 12, the red signals indicate the value of the statistic extracted from the original data for both prey and predator time series. These values are completely separable from the surrogate distribution. It means that the red signals are in the confidence interval, so the null hypothesis is rejected 共nonchaotic behavior rejected兲. This criterion, therefore, shows clearly that the simulation’s population data have a chaotic behavior. D. P&H method To determine if the P&H method value of a typical time series presents a nonchaotic or chaotic behavior, the same hypothesis testing procedure as in Sec. IV C is adopted. The hypothesis under H0 models a nonchaotic behavior and the hypothesis H1 implies a chaotic behavior. The same results as in the previous section are observed for the random and Lorenz time series 共data not shown兲. In Fig. 13, the red signals indicate the value of the statistic extracted from the original data both for prey and predator time series. This value is completely separable from the surrogate distribution. It means that the red signals are in the confidence interval, so the null hypothesis is rejected 共nonchaotic behavior rejected兲. This last criterion confirms that the simulation’s population data have a chaotic behavior.59 V. CONCLUSION The purpose of this paper is the examination of the stochastic and deterministic behavior of signals that are produced by an ecosystem simulation. First, we gave a brief overview of the existing simulation. It is an individual-based evolving predator-prey ecosystem simulation, which uses fuzzy cognitive maps as a behavior model. We are interested to understand how complex and predictable our simulation is. The variations of the time series associated to ecosystem simulation are the result of complex simulation mechanisms. In this study, we focus on the time series corresponding to the variation of the number of prey and predator individuals. To understand how close our system is to the random or chaotic processes, we examine whether a chaotic behavior exists in these signals. To enforce our result, we use four different methods: Higuchi fractal dimension, correlation dimension, largest Lyapunov exponent, and P&H method. For each of them, in order to obtain a statistically significant evaluation, we apply the surrogate test method on 24 samplings of the considered data. According to the results obtained after applying these different methods, all of them providing clear predictions, we can conclude that the behavior of simulation is deterministic. This has been shown by Higuchi method and correlation dimension. Also among various cases of deterministic behavior, we show that the behavior of simulation is chaotic. This has been shown by Lyapunov exponent and P&H methods. It is important to show the existence of a complex behavior in our simulation because any attempt to model a realistic system needs to have the capacity to generate patterns as complex as the ones that are observed in real systems. In this regard, we will be interested to analyze other Chaos 20, 043120 共2010兲 versions of the current simulation or even other different simulations. This will give us the possibility to better understand what important factors lead to complex systems with behavior showing an important level of regularity but still very hard to predict. Using the same process, it is also possible to examine other time series generated by our simulation. It can be useful to determine among the huge amount of data generated by the simulation the ones that are more interesting to analyze more deeply. ACKNOWLEDGMENTS This work was supported by the NSERC 共Grant No. ORGPIN 341854兲, the CRC 共Grant No. 950-2-3617兲, and the CFI 共Grant No. 203617兲, and was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network 共SHARCNET: www.sharcnet.ca兲. 1 J. Hubbard and B. West, Differential Equations: A Dynamical Systems Approach: Ordinary Differential Equations 共Springer, Berlin, 1991兲, Vol. 5. 2 P. Grassberger, T. Schreiber, and C. Schaffrath, Int. J. Bifurcation Chaos Appl. Sci. Eng. 1, 521 共1991兲. 3 A. H. Omidvarnia and A. Nasrabadi, Fractals 16, 129 共2008兲. 4 A. Brandstäter, J. Swift, and H. L. Swinney, Phys. Rev. Lett. 51, 1442 共1983兲. 5 P. Berge, Y. Pomeau, and C. Vidal, Order within Chaos: Towards a Deterministic Approach to Turbulence 共Wiley, New York, 1984兲. 6 T. Higuchi, Physica D 31, 277 共1988兲. 7 T. Higuchi, Physica D 46, 254 共1990兲. 8 L. Telesca, G. Colangelo, V. Lapenna, and M. Macchiato, Chaos, Solitons Fractals 18, 385 共2003兲. 9 W. Klonowski, E. Olejarczyk, and R. Stepien, Attractors, Signals, and Synergetics, Frontiers of Non-linear Dynamics 共Pabst Science, Lengerich, 2002兲, Vol. 1. 10 W. Klonowski, “Everything you wanted to ask about EEG but were afraid to get the right answer,” Nonlinear Biomed. Phys. 3, 2 共2009兲. 11 A. Anier, T. Lipping, S. Melto, and S. Hovilehto, “Higuchi fractal dimension and spectral entropy as measures of depth of sedation in intensive care unit,” Proceedings of the 26th IEEE EMBS Conference, 2004, pp. 526–529. 12 F. Cervantes-De la Torre, C. G. Pavía-Miller, A. Ramirez-Rojas, and F. Angulo-Brown, “Time evolution of the fractal dimension of electric selfpotential time series,” Nonlinear Dynamics in Geosciences 共SpringerVerlag, New York, 2007兲, pp. 407–418. 13 A. Anier, T. Lipping, S. Melto, and S. Hovilehto, Proceedings of the 26th Annual International Conference of the IEEE IEMBS APOS ‘04, 2004, Vol. 1, pp. 526–529. 14 A. Accardo, M. Affinito, and M. Carrozzi, Biol. Cybern. 77, 339 共1997兲. 15 W. Klonowski, “Chaotic dynamics applied to signal complexity in phase space and in time domain,” Chaos, Solitons Fractals 14, 1379 共2002兲. 16 S. Spasić, “Surrogate data test for nonlinearity of the rat cerebellar electrocorticogram in the model of brain injury,” Signal Process. 90, 3015 共2010兲. 17 X. Li, Z. Deng, and J. Zhang, “Function of EEG temporal complexity analysis in neural activities measurement,” Advances in Neural Networks 共Springer-Verlag, Berlin, 2009兲, Vol. 5551, pp. 209–218. 18 J. Virkkala, S. L. Himanen, A. Värri, and J. Hasan, “Fractal dimension of EEG in sleep onset,” Proceedings of the Third European Interdisciplinary School on Nonlinear Dynamics for System and Signal Analysis, Warsaw, Poland, 2002. 19 M. T. R. Peiris, R. D. Jones, P. R. Davidson, P. J. Bones, and D. J. Myall, “Fractal dimension of the EEG for detection of behavioural microsleeps,” Proceedings of the 27th International Conference on IEEE Engineering in Medicine and Biology Society, Shanghai, China, 2005. 20 J. Z. Liu, Q. Yang, B. Yao, R. W. Brown, and G. H. Yue, “Linear correlation between fractal dimension of EEG signal and handgrip force,” Biol. Cybern. 93, 131 共2005兲. 21 B. S. Raghavendra, D. N. Dutt, H. N. Halahalli, and J. P. John, “Complexity analysis of EEG in patients with schizophrenia using fractal dimension,” Physiol. Meas. 30, 795 共2009兲. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions 043120-13 22 Regularity analysis of a simulation S. Georgiev, Z. Minchev, C. Christova, and D. Philipova, “EEG fractal dimension measurement before and after human auditory stimulation,” BIOAUTOMATION 12, 70 共2009兲. 23 R. Sullivan, T. Holden, G. Tremberger, Jr., E. Cheung, C. Branch, J. Burrero, and G. Surpris, “Fractal dimension of breast cancer cell migration in a wound healing assay,” Int. J. Biol. Life Sci. Technol. 共IJLST兲 6, 3 共2010兲. 24 C. Gómez, A. Madiavilla, R. Hornero, D. Abásalo, and A. Fernández, “Use of Higuchi’s fractal dimension for the analysis of MEG recordings from Alzheimer’s disease patients,” Med. Eng. Phys. 31, 306 共2009兲. 25 T. L. A. Doyle, E. L. Dugan, B. Humphries, and R. U. Newton, “Discriminating between elderly and young using a fractal dimension analysis of centre of pressure,” Int. J. Med. Sci. 1, 11 共2004兲. 26 W. Klonowski, Nonlinear Biomed. Phys. 1, 5 共2007兲. 27 M. Moradi, P. Abolmaesumi, and A. Phillip, Proceedings of the 28th Annual International Conference of the IEEE EMBS APOS ‘06, 2006, pp. 2400–2403. 28 B. S. Raghavendra and D. Narayana Dutt, “Nonlinear dynamical characterization of heart rate variability time series of meditation,” Int. J. Life Sci. Technol. 共IJLST兲 1, 167 共2010兲. 29 Y. Zhao and M. Small, “Deterministic propagation of blood pressure waveform from human wrists to fingertips,” Lect. Notes Comput. Sci. 3177, 142 共2004兲. 30 R. Gras, D. Devaurs, A. Wozniak, and A. Aspinall, “An individual-based evolving predator-prey ecosystem simulation using fuzzy cognitive map as behavior model,” Artif. Life 15, 423 共2009兲. 31 R. Sneyers, “Climate chaotic instability: Statistical determination and theoretical background,” Environmetrics 8, 517 共1997兲. 32 L. Romanelli, M. A. Figliola, and F. A. Hirsch, “Deterministic chaos and natural phenomena,” J. Stat. Phys. 53, 991 共1988兲. 33 A. Golestani, M. R. Jahed Motlagh, K. Ahmadian, A. H. Omidvarnia, and N. Mozayani, “A new criterion for distinguish stochastic and deterministic time series with the Poincaré section and fractal dimension,” Chaos 19, 013137 共2009兲. 34 B. Kosko, “Fuzzy cognitive maps,” Int. J. Man-Mach. Stud. 24, 65 共1986兲. 35 D. Devaurs and R. Gras, “Species abundance patterns in an ecosystem simulation studied through Fisher’s logseries,” Simulation Modelling Practice and Theory 共Elsevier, Amsterdam, 2010兲, pp. 100–123. 36 V. Grimm, “Ten years of individual-based modelling in ecology: What have we learned and what could we learn in the future?,” Ecol. Modell. 115, 129 共1999兲. 37 D. L. DeAngelis and W. M. Mooij, “Individual-based modeling of ecological and evolutionary processes,” Annu. Rev. Ecol. Evol. Syst. 36, 147 共2005兲. 38 R. E. Lenski, C. Ofria, T. C. Collier, and C. Adami, “Genome complexity, robustness, and genetic interactions in digital organisms,” Nature 共London兲 400, 661 共1999兲. 39 C. R. Ward, F. Gobet, and G. Kendall, “Evolving collective behavior in an artificial ecology,” Artif. Life 7, 191 共2001兲. Chaos 20, 043120 共2010兲 40 J. Mallet, “A species definition for the modern synthesis,” Trends Ecol. Evol. 10, 294 共1995兲. 41 A. Aspinall and R. Gras, “K-means clustering as a speciation method within an individual-based evolving predator-prey ecosystem simulation,” The International Computer Science Conference—Active Media Technology 共AMT ‘10兲, 2010. 42 A. Golestani, A. Ashouri, K. Ahmadian, M. Jahed, and M. A. Doostari, “Irregularity analysis of iris patterns,” WORLDCOMP ‘08 Congress 共IPCV兲, 14–17 July 2008, pp. 691–695. 43 M. Small, “Applied nonlinear time series analysis: Applications in physics, physiology and finance,” Nonlinear Science Series A 共Singapore, World Scientific, 2005兲, Vol. 52. 44 D. Yu, M. Small, R. G. Harrison, and C. Diks, “Efficient implementation of the Gaussian kernel algorithm in estimating invariants and noise level from noisy time series data,” Phys. Rev. E 61, 3750 共2000兲. 45 C. Diks, “Estimating invariants of noisy attractors,” Phys. Rev. E 53, R4263 共1996兲. 46 A. Tsonis, Chaos—From Theory to Applications 共Plenum, New York, 1992兲. 47 N. B. Tufillaro, Poincare Map, 1997. 48 S. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering 共Addison-Wesley, Reading, MA, 1994兲. 49 F. Takens, Dynamical Systems and Turbulence, Lecture Notes in Mathematics Vol. 898 共Springer, Berlin, 1981兲, p. 366. 50 J. Clinton Sprott, Chaos and Time Series Analysis 共Oxford University Press, New York, 2003兲. 51 P. Grassberger and I. Procaccia, Physica 共Amsterdam兲 9D, 189 共1983兲. 52 M. T. Rosenstein, J. J. Collins, and C. J. De Luca, “A practical method for calculating largest Lyapunov exponents from small data sets,” Physica D 65, 117 共1993兲. 53 U. Parlitz, “Nonlinear time-series analysis,” Nonlinear Modeling— Advanced Black-Box Techniques, edited by J. A. K. Suykens and J. Vandewalle 共Kluwer, Dordrecht, 1998兲, pp. 209–239. 54 S. Sato, M. Sano, and Y. Sawada, “Practical methods of measuring the generalized dimension and largest Lyapunov exponent in high dimensional chaotic systems,” Prog. Theor. Phys. 77, 1 共1987兲. 55 J. Kurths and H. Herzel, “An attractor in solar time series,” Physica D 25, 165 共1987兲. 56 M. Dämmig and F. Mitschke, “Estimation of Lyapunov exponent from time series: The stochastic case,” Phys. Lett. A 178, 385 共1993兲. 57 A. Ray and R. Chowdhury, “On the characterization of non-stationary chaotic systems: Autonomous and non-autonomous cases,” Physica A 389, 5077 共2010兲. 58 J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, and J. D. Farmer, “Testing for nonlinearity in time series: The method of surrogate data,” Physica D 58, 77 共1992兲. 59 C. Merkwirth, U. Parlitz, I. Wedekind, and W. Lauterborn, TSTOOL, Version 1.11, 2002, downloadable from http:// www.physik3.gwdg.de/tstool/index.html. Downloaded 11 Apr 2011 to 137.207.140.84. Redistribution subject to AIP license or copyright; see http://chaos.aip.org/about/rights_and_permissions