Many biological ecosystems exhibit chaotic behavior, demonstrated either analytically using parameter choices in an associated dynamical systems model or empirically through analysis of experimental data. In this paper, we use existing software tools (COPASI, R) to explore dynamical systems and uncover regions with positive Lyapunov exponents where thus chaos exists. We evaluate the ability of the software’s optimization algorithms to find these positive values with several dynamical systems used to model biological populations. The algorithms have been able to identify parameter sets which lead to positive Lyapunov exponents, even when those exponents lie in regions with small support. For one of the examined systems, we observed that positive Lyapunov exponents were not uncovered when executing a search over the parameter space with small spacings between values of the independent variables.
Chaos is known to occur in many natural biological systems involving population dynamics [10, 12, 13, 17]. For example, Becks et al. demonstrated existence of chaotic and non-chaotic states in an experimental study of a microbial predator-prey system with a nutrient source . The authors used daily measurements of the population size of each species to document transitions between chaotic and non-chaotic states after changing the strength of the nutrient source. Another experimental study involving baker’s yeast and a fermenter  also observed chaotic fluctuations in the growth rate of the fungus. The paper demonstrated chaotic behavior through time series analysis.
Even simple dynamical system models are known to exhibit chaotic behavior for certain values of system parameters [19, 28]. However, even such simple systems can possess a wide range of chaotic and non-chaotic behavior, related in a non-trivial way to certain combinations of system parameters and initial conditions. As a result, detecting and analyzing chaotic behavior is far from straightforward and requires a high degree of mathematical sophistication. This can pose a significant challenge to modelers with a mostly biological background and hinder the discussion of chaotic dynamics in biological systems for educational purposes.
In this work, we evaluate the ability of standard computational tools to find subsets of the parameter space which give rise to positive Lyapunov exponents. In particular, we use optimization algorithms to efficiently search the parameter space and return values which define arrangements of dynamical systems giving computationally positive Lyapunov exponents. The framework uses numerical techniques to generate approximate solutions to the dynamical system coupled with optimization algorithms to uncover connections between parameters and highlight regions of chaos.
We pair the results from the optimization with visualization tools available in R to allow investigators to easily establish relationships between parameter choices and identify ranges of parameter values leading to chaotic behavior. The aim is to help biologists develop a better understanding of the causes of chaotic behavior in dynamical systems without having to master the intricate mathematical analysis often required to rigorously detect and analyze chaotic behavior. The main novelty of the approach is the use of optimization algorithms and high-dimensional visualization techniques to assist with identification and understanding of parameters leading to chaotic dynamics.
We validate our computational approach by applying it to population models in the literature. The equations in the systems we considered contain rate-limited functions to model predator/prey dynamics. Systems of this type have been well studied in chemostat environments, both experimentally and analytically [3, 7, 17, 19]. They have been shown to demonstrate chaotic dynamics [7, 6, 14] and can be considered as small-scale studies for larger ecological systems . In particular, we consider the set of coupled differential equations
Given the wide applicability of these models to different physical and ecological systems, we seek to better understand the choices of model parameters and initial conditions leading to chaotic behavior. Let P ⊆ Rm denote all possible valid settings for the m parameters in our system (e.g., one parameter might describe the dilution rate which is the rate at which a nutrient is introduced into the system). Let I ⊆ Rn denote the set of all valid initial conditions for (Y1, Y2, …, Yn). Defining S = P × I, we seek to characterize the regions of S leading to chaotic behavior. A particular region s ∈ S is deemed chaotic if it has negative divergence at its initial conditions (to ensure boundedness) and yields a positive Lyapunov exponent. In our computations we use standard numerical integrators for ordinary differential equations along with popular methods for numerical Lyapunov exponent calculation [8, 31].
In the following sections, we describe the algorithms and visualizations used in our approach, we provide the results from the validation exercises, and we summarize our findings and provide discussion on future research directions.
2Description of the framework
The framework discussed in this paper uses the freely available software tools COPASI (a COmplex PAthway SImulator)  and R  to compute the Lyapunov exponents, evaluate the objective functions, and construct parallel coordinates visualizations of the results. The flowchart in Fig. 1 shows the progression of the algorithm.
Let Y′ = f (t, Y, s) represent a dynamical system with parameters s ∈ S. If the system models population dynamics, for instance, the parameters can contain birth and death rates along with initial conditions. We take S to be the design space in this framework.
The usual test for chaos is calculation of the largest Lyapunov exponent; a positive largest Lyapunov exponent indicates chaos . Lyapunov exponents provide a quantitative measure of the convergence or divergence of nearby trajectories for a dynamical system. The system is said to be chaotic over s ∈ S if the maximum Lyapunov exponent L (s) is positive and the diverence of the system D (s) is negative . The divergence of the system is defined as the trace of the Jacobian associated with the system evaluated at the initial state. That is, the system exhibits chaotic dynamics if L (s) >0 and D (s) <0. Thus, our overall objective is to efficiently search S to find members s in the space for which the maximum Lyapunov exponent L (s) >0 and for which the divergence D (s) <0.
Strategies for computing the Lyapunov spectrum for a dynamical system have been known for some time . At the heart of these algorithms is an efficient computational technique for generating accurate numerical approximations to the system. Our interest is in the inverse problem; that is, given a dynamical system with an associated parameter space, is it possible to use numerical tools to find elements in the space which result in a positive Lyapunov exponent for the system? In addition, we wanted to gain insight into the possible connections between the parameter values leading to these positive exponents. Thus, we used visualization tools to generate graphical representations of these associations.
The algorithm in Wolf, et al.  actually computes the entire Lyapunov spectrum from the discretized solution to the system of ordinary differential equations (ODEs) and the associated linearized equations of motion. The linearized equations of motion are formed from the Jacobian of the system of ordinary differential equations. We use the implementation of this algorithm available in COPASI.
The objective function defined for the optimization routine is to maximize the maximum Lyapunov exponent for a given dynamical system Y′ = f (t, Y, s), over a paramter set s. That is, starting with the definition of the dynamical system and the initial parameter set s0, the optimization algorithm computes the maximum Lyapunov exponent and the divergence of the system associated with that particular parameter set. The optimization problem is constrained by the divergence. That is, the Lyapunov exponent must be positive, and the divergence negative, for this parameter set to be feasible. We note most optimization algorithms minimize objective functions, so our objective function is actually the negative of the maximum Lyapunov exponent.
Based on the information returned from the evaluation of the objective function, the optimization routine then chooses another parameter set s1 and again computes the maximum Lyapunov exponent and divergence for that parameter set. This process continues until the optimization algorithm exhausts its set of function evaluations or determines progress is no longer being made in maximizing (or minimizing) the objective.
We chose to use the genetic algorithm available in COPASI  as the optimization algorithm. Genetic algorithms move through “generations" of design points by evaluating the fitness of members of the generation and selecting members to continue to the next generation (through mutation or cloning), parent offspring for the next generation, or die (i.e., these points are removed from the population) . The algorithm was run multiple times, with different initial choices for parameters, to look for multiple design points.
These results were validated with the simulated annealing optimization routine also available in COPASI. The simulated annealing algorithm randomly selects design points and keeps those which are associated with lower objective function evaluations. In order to make the search global, the algorithm also selects a few random points which increase the value of the objective function. The algorithm moves to a new state or remains in a given state with a specified probability, with an increased probability associated with those states giving lower values of the objective function. The behavior of the algorithm is defined by choosing parameters which define the size of the initial search region and reduce the size of the search space as the algorithm progresses [1, 9].
Let C ⊆ S denote the subset of all parameters and initial conditions leading to chaotic behavior. Now that we have a method for efficiently generating samples from the chaotic regime C, we use an interactive visualization platform to help the user understand the structure of C —a challenging task, due to the potentially complex, high-dimensional nature of this set. Our platform is based on parallel coordinates, a popular means of visualizing high-dimensional data. The plots were generated using the plotly  package in R. The input for the parallel coordinates was obtained from the result files from the optimization task of COPASI. All COPASI and R files can be found in https://github.com/skoshyc/Computational-Framework-Chaotic-Dynamics.
Parallel coordinate plots allow us to understand a number of useful properties by visual inspection and interactive manipulation. For example, in Fig. 2, we see that chaos can be obtained for log (r) ∈ (4.55, 4.6) and values of b ∈ (2.60, 2.63). The purple line depicts one particular combination of the parameters and the Lyapunov exponent. The faint lines are combinations outside the selected range of log (r). By restricting several coordinates at a time, the user can filter an initially large number of sample points down to only a few.
The numerical results presented are intended to demonstrate the ability of the framework to meet the stated goal of this work; i.e., provide researchers with a convenient numerical tool for exploring the parameter space and uncovering regions where parameter choices lead to chaos. We validated the framework on several dynamical systems from the literature. We include results from two of those studies here [14, 19].
Both of the included systems can be used to model predator-prey dynamics in ecosystems. The mathematical analysis of both systems uncovered parameter values which gave chaotic solutions [14, 19]. We sought to replicate these findings by maximizing the Lyapunov exponent over a design space consisting of the initial conditions for the variables and two of the parameters used to define the ODEs. We use this section to describe the equations and the results of our study.
All computations were carried out on a Windows 10 PC with Intel(R) CPU and 4 Cores. The version of COPASI was 4.27 and R was 4.0.0. Both the genetic algorithm and simulated annealing were run multiple times to obtain multiple sample points. Each run had randomized start values which is an option in the optimization task in COPASI. After each run, the maximum Lyapunov exponent along with the optimal parameters are generated. COPASI collects the information from each run and provides the output from all the runs in a text file. The positive Lyapunov exponents and associated parameter set are chosen by R and visualized using the parallel coordinates plot.
We used 200 generations with 20 population members in each genetic algorithm search to generate the results provided below. The algorithm was run 100 times to generate the parallel coordinates plots. For the simulated annealing results, we used the COPASI default parameter values for the starting temperature, the cooling factor, and the tolerance . The start temperature defines the size of the random perturbations on the parameters for the algorithm, and the cooling factor defines the size of the reduction of the random perturbations at each cycle. The algorithm stops when the change in the objective function has been smaller than the tolerance in the last two temperature steps. The simulated annealing algorithm was run 25 times in order to have a defined number of points to compare with the results from the genetic algorithm. The number of repeated runs was limited due to the computational expense of the simulated annealing algorithm.
3.1Double forced Monod system (DFMS)
We first consider one of the dynamical systems presented and analyzed in work by Kot, Sayler, and Schultz . While experimental results are not included in the paper, the authors justify their choice to study this particular system by noting the possibility of obtaining experimental validation of their work. As our primary interest in developing the framework presented here is to aid biologists in their data-based studies of dynamical systems, this particular problem presents an ideal benchmark case.
The dimensionalized equations for the DFMS are given as
Si is the inflowing substrate concentration where i denotes inflow of the substrate into the chemostat. We note the parameter ε changes the dynamics of the inflowing substrate. Smaller values of ε lead to more constant inflows, while larger values of ε introduce more fluctuation into the inflow. Both ε and the frequency of the associated oscillation were important parameters for chaos, indicating nutrition rates play an important role here.
The equations are nondimensionalized by rescaling all variables by the inflow substrate, the prey by its yield constant Y1, and the predator by both yield constants . The dimensionless system is given by
As shown in the referenced paper, the dimensionless system exhibits chaotic behavior for and ε = 0.6, with Si = 115 mg/l, D = 0.1 / h, and initial conditions x (0) =0.42, y (0) =0.4, z (0) =0.42. The remainder of the parameter values are provided in Table 1. Representative manifold plots are given in  and are not duplicated here.
|Yi||μi (h- 1)||Ki (mg/l)|
|Prey (i = 1)||0.4||0.5||8|
|Predator (i = 2)||0.6||0.2||9|
To determine if our framework could recover the values cited in , we used the optimization algorithms to maximize the computed Lyapunov exponent as a function of ε, ω, and the initial states for x (τ), y (τ), and z (τ). A subset of our results is shown using parallel coordinates visualization in Fig. 3. This visualization indicates more positive Lyapunov exponents associated with small values of ω. Note in the rescaled system ω is the angular frequency of the forcing term. Thus, observing positive Lyapunov exponents for small values of ω may indicate to a biologist an interesting connection between periodic forcing and chaotic behavior that they may be able to further study experimentally. Any clustering of the remainder of the parameters chosen for our study is less obvious. Figure 4 contains a manifold plot of the dimensionless system for a randomly chosen set of parameter values from those discovered by our framework; namely, ε = 0.753444, ω = 0.455121, and initial conditions x (0) =4.49844, y (0) =0.428112, and z (0) =0.446482. We note that the strange attractor shown in Fig. 4 is distinct from the one given in the referenced paper.
To further investigate the efficacy of having a targeted search of the parameter space, we also used the simulated annealing approach of COPASI to generate a comparison of Lyapunov exponents with the results of the genetic algorithm. The initial conditions were set to (x, y, z) = (1, 1, 1) to observe the behavior of only the two parameters ε and ω. Each of the algorithms were run 25 times for the comparison. The comparison was repeated twice to account for stochasticity. As is seen in Fig. 5, there are clusters of positive Lyapunov exponents (and hence chaos) in the ranges of ε ∈ (0.18, 0.21) and ω ∈ (0.8, 1).
We sought to better understand these results by computing the maximum Lyapunov exponent over a 100 × 100 discretization of the region ε × ω for ε ∈ (0.18, 0.21), ω ∈ (0.9, 1) using the initial condition (1, 1, 1). The task was completed using the parameter scan option in COPASI. The results are shown in Fig. 6. The graph shows all of the largest Lyapunov exponents at these discrete points were negative, indicating the optimization algorithms were able to locate these small regions where positive Lyapunov exponents were computed. One of the reasons we possibly don’t get the requisite parameters in the 100 × 100 grid search is because the grid is regularly spaced. For example, two of the values associated with the cluster in Fig. 5 are (ε, ω) = (0.180791, 0.955417) , (0.181136, 0.913295). Such numbers do not appear in the regularly spaced grid search.
3.2A three-species food chain
The second example uses a three-species food chain model where chaos was observed for a given set of biologically relevant parameters . The model was a continuous, non-linear system which has biological implications of chaos on a general food web model. Our algorithm was able to find the same parameter values reported in the paper along with other parameter choices which lead to chaos.
This three-species model  is defined using the equations
The authors created a nondimensional system of equations using the following transformations:
For our design space, we considered b1 ∈ (2.5, 4), b2 ∈ (1, 3), and ranges for the initial conditions x (0) ∈ (0.1, 1), y (0) ∈ (0.1, 1), and z (0) ∈ (7.5, 10.5).
The results obtained from this search are provided in Table 3 and Fig. 7. The authors of  note chaos is observed for larger values of b2 (i.e. larger than 2), which we observe in the parallel coordinates plot as well. We also see that most of the positive Lyapunov exponents are found when b1 and b2 are both close to 2.5.
|Lyapunov exponent||b1||b2||x (0)||y (0)||z (0)|
A comparison of simulated annealing and the genetic algorithm was performed to better understand the landscape of the Lyapunov exponents as a function of (b1, b2) with the initial conditions set to (0.5, 0.3, 9). As seen in Fig. 8, the majority of the parameter combinations which show chaos are centered around (b1 = 2.5, b2 = 2.5).
We have provided a computational strategy that allows researchers to efficiently search over a design space including parameter values and initial conditions and discover possible connections between values of these constants and chaotic dynamics. The ability for such a search may reveal chaotic behavior in systems not previously known to have chaotic regime and the existence of parameters and initial conditions not previously known to yield chaotic behavior in studied systems. In one of the models we investigated, our technique found chaotic regions which were not observed using a regularly spaced grid search. This showed that to obtain the chaotic parameters, we would need a finer grid and thus several more iterations. Another brute force technique such as a random search may also miss the chaotic parameters. The use of COPASI and R results in a robust, comprehensive framework for analyzing a given dynamical system without the need for rigorous mathematical analysis.
One of the main strengths of this approach is that it allows the study of system behavior that is virtually impossible to observe in a laboratory environment, thus making it useful to experimentalists. The parallel coordinates plots provide further insight into system dynamics for a wide variety of parameter changes.
Parallel computations will improve the efficiency of this strategy. The genetic algorithm is conceptually easy to parallelize, in that objective function evaluations for a generation can be easily distributed across processors. We have not evaluated the speed-up which would be associated with this modification. We also continue to evaluate the ability of different objective functions to locate parameter values associated with chaos. One of the Lyapunov exponents being positive is the chaos indicator that we employ in this paper. The optimization algorithm which searches for the maximal Lyapunov exponent will find the conditions which give the most positive exponent. As the Lyapunov exponent can be a discontinuous function of the design space, it can be easy for an optimization algorithm based on properties of continuous functions to miss them. Thus, better strategies for finding these points are warranted. If the system has two positive Lyapunov exponents, we are limited in the use of an optimization algorithm as the algorithm will probably give us the parameters which give rise to the more positive exponent thus missing other regions of chaos. However, multiple runs of global optimization algorithms like the genetic algorithm, simulated annealing, evolutionary algorithms in COPASI would give us different parameter combinations of the desired chaotic regions.
The authors express gratitude to Dr. Oleg Yordanov for his guidance and help with appropriate system rescaling; Akshay Galande for discussions about genetic algorithms; and to Dr. Pedro Mendes for his help with troubleshooting COPASI.
http://copasi.org/Support/User_Manual/Methods/Optimization_Methods/Simulated_Annealing/. [Last accessed: Jan 24, 2020].
Kathleen T. , Alligood T.S. and Yorke James A. , Chaos: an introduction to dynamical systems. Springer, 1997.
Banks H.T. , Modeling and Control in the Biomedical Sciences, volume 6 of Lecture Notes in Biomathematics. Springer-Verlag, Heidelberg, 1975.
Banks H.T. and Davidian M. , Introduction to the chemostat. http://www.ncsu.edu/crsc/htbanks/STMA810c/chemostat.pdf, 2009.
Barrio R. and Serrano S. , A three-parametric study of the lorenz model, Physica D: Nonlinear Phenomena 229(1) (2007), 43–51.
Becks L. and Arndt H. , Transitions from stable equilibria to chaos, and back, in an experimental food web, Ecology 89 (2008), 3222–3226.
Becks L. , Hilker F.M. , Malchow H. , Jürgens K. and Arndt H. , Experimental demonstration of chaos in a microbial food web, Nature 435 (2005), 1226–1229.
Benettin G. , Galgani L. , Giorgilli A. and Strelcyn J.-M. , Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems: a method for computing all of them, Meccanica 15(1) (1980), 21–30.
Corana A. , Marchesi M. , Martini C. and Ridella S. , Minimizing multimodal functions of continuous variables with the simulated annealing algorithm, ACM Trans Math Softw 13 (1987), 262–280.
Costantino R. , Desharnais R. , Cushing J. and Dennis B. , Chaotic dynamics in an insect population, Science 275 (1997), 389–391.
Davey Hazel M. , Davey Christopher L. , Woodward Andrew M. , Edmonds Andrew N. , Lee Alvin W. and Kell Douglas B. , Oscillatory, stochastic and chaotic growth rate uctuations in permittistatically controlled yeast cultures, Biosystems 39(1) (1996), 43–61.
Dennis B. , Desharnais R. , Cushing J. and Costantino R. , Transitions in population dynamics: equilibria to periodic cycles to aperiodic cycles, J Anim Ecol 66 (1997), 704–729.
Fussmann G. , Ellner S. , Shertzer K. and Hairston N. , Crossing the Hopf bifurcation in a live predator-prey system, Science 290 (2000), 1358–1360.
Hastings A. and Powell T. , Chaos in a three-species food chain, Ecology 72(3) (1991), 896–903.
Henson M.A. , Dynamic modeling of microbial cell populations, Curr Opin Biotech 14 (2003), 460–467.
Hoops S. , Sahle S. , Gauges R. , Lee C. , Pahle Jürgen , Simus N. , Singhal M. , Xu L. , Mendes P. and Kummer U. , Copasi|a complex pathway simulator, Bioinformatics 22(24) (2006), 3067–3074.
Jost J. , Drake J. , Frederickson A. and Tsuchiya H. , Interactions of Tetrahymena pyrifarmis, Escherichia coli, Azobacter vinelndil, and glucose in a mineral medium, J Bacteriol 113 (1973), 834–840.
Kooi B.W. and Boer M.P. , Chaotic behavior of a predator-prey system in the chemostat, Dynam Cont Dis Ser B 10(2) (2003), 259–272.
Kot M. , Sayler G.S. and Schultz T.W. , Complex dynamics in a model microbial system, B Math Biol 54(4) (1992), 619–648.
Kovárová-Kovar K. and Egli T. , Growth kinetics of suspended microbial cells: From singlesubstrate-controlled growth to mixed-substrate kinetics, Microbiol Mol Biol R 62(3) (1998), 646–666.
Kravchenko L.V. , Strigul N.S. and Shvytov I.A. , Mathematical simulation of the dynamics of interacting populations of rhizosphere microorganisms, Microbiology 73(2) (2004), 189–195.
LorenzEdward N., Deterministic nonperiodic flow, Journal of the atmospheric sciences 20(2) (1963), 130–141.
Michaelis L. and Menten M.L. , Die kinetik der invertinwirkung, Biochemistry Z 49 (1913), 333–369.
Mitchell M. , An Introduction to Genetic Algorithms. MIT Press, Boston, 1995.
Molz F. and Faybishenko B. , Increasing evidence for chaotic dynamics in the soil-plantatmosphere system: A motivation for future research, Procedia Environmental Sciences 19 (2013), 681–690.
R Core Team, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2020.
Sievert C. , Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC, 2020.
Sprott J.C. , Some simple chaotic ows, Phys Rev E 50(2) (1994), 647–650.
Sprott J.C. , Chaos and Time-Series Analysis. Oxford University Press, Oxford, UK, 2003.
Strigul N.S. and Kravchenko L.V. , Mathematical modeling of PGPR inoculation into the rhizosphere, Environ Modell Softw 21 (2006), 1158–1171.
Wolf A. , Swift J.B. , Swinney H.L. and Vastano J.A. , Determining Lyapunov exponents from a time series, Physica D 16 (1985), 285–317.