A phase transition induces chaos in a predator-prey ecosystem with a dynamic fitness landscape

In many ecosystems, natural selection can occur quickly enough to influence the population dynamics and thus future selection. This suggests the importance of extending classical population dynamics models to include such eco-evolutionary processes. Here, we describe a predator-prey model in which the prey population growth depends on a prey density-dependent fitness landscape. We show that this two-species ecosystem is capable of exhibiting chaos even in the absence of external environmental variation or noise, and that the onset of chaotic dynamics is the result of the fitness landscape reversibly alternating between epochs of stabilizing and disruptive selection. We draw an analogy between the fitness function and the free energy in statistical mechanics, allowing us to use the physical theory of first-order phase transitions to understand the onset of rapid cycling in the chaotic predator-prey dynamics. We use quantitative techniques to study the relevance of our model to observational studies of complex ecosystems, finding that the evolution-driven chaotic dynamics confer community stability at the “edge of chaos” while creating a wide distribution of opportunities for speciation during epochs of disruptive selection—a potential observable signature of chaotic eco-evolutionary dynamics in experimental studies.


Supplementary Analysis A Solutions in the absence of evolutionary dynamics
We study the system along the two-dimensional nullclineċ = 0, wherec is a parameter representing some fixed mean trait value in the population. First, we consider the further reduced case of the prey density in the absence of predation,ẏ = y = 0.
When b 1 d 1 < a 1 , the prey density grows exponentially ifc > d 1 /(a 1 − b 1 d 1 ); it decreases exponentially whenc < d 1 /(a 1 − b 1 d 1 ). Because a nonzero predator density only serves to decrease the rate of the prey growth, this constraint also applies to the two-dimensional case (A1, A2). These equations admit two solutions. The first is mutual exclusion, which has associated eigenvalues, As predicted for the one-dimensional problem, the mutual exclusion solution is stable when 15c < d 1 /(a 1 − b 1 d 1 ). However, even ifc > d 1 /(a 1 − b 1 d 1 ), this solution is stable if competition is strong enough that b 1 > a 1 /d 1 .
The second solution to (A1, A2) is an interior point, .
The associated eigenvalues for this solution occur in pairs, ) . (A8) The requirement thatŷ,x ≥ 0 means that this interior solution exists only whenc > d 1 /(a 1 − b 1 d 1 ) (which, as before, entails exponential prey growth in the absence of the predator) and b 1 < a 1 /d 1 , b 2 < a 2 /d 2 . The latter two conditions are equivalent to the non-existence of 20 coexistence solutions under strong competition in the classical Lotka-Volterra predator-prey model. For these conditions, the real parts of (A8) are always positive, and so the interior solution is never stable.
The eigenvalues of (A8) always have nonzero imaginary components if either When either of these conditions is satisfied, cycling is possible in the system.

B Hysteresis and critical points
Under the assumption that evolutionary dynamics are fast enough thatc ≈ c near the maximum of the fitness landscape r(x, y, c,c), then c eq are given by the solutions of the The roots of this equations are intricate expressions; we define first the auxiliary variable β: The roots can now be written in the form, where the latter two equilibria are complex conjugates of one another.

30
The maximum value ofc(t) can be found by determining by the maximum possible value of c eq , which occurs when x = 0 in c For the parameter values used here, max(c eq ) = √ 2/2 ≈ 0.707107 The first turning point, x * , is found by determining the positive value of x at which the two positive equilibria are equal (c , namely ) .

Inserting this value into either c
(3) eq or c (4) eq yields an estimate of c * . For the parameter values used here, x * ≈ 0.4495, c * ≈ 0.3565.

35
The second turning point, x * * , is found by determining the point where the unstable equilibrium c (4) eq first crosses the x axis,

Inserting this equation into c
(3) eq yields the value of c * * , the point to which the equilibrium value ofc jumps when x reaches x * * from above. For the parameter values used here,

C Calculation of Global Lyapunov Exponents
Lyapunov exponents were calculated numerically using the "renormalization" algorithm originally described by Bennettin et al. 1,2 First, a long trajectory (T = 40, 000) was generated from an arbitrary initial condition in order to sample a large range of points on the attractor.
Then, N locations on this attractor were randomly chosen, and the Lyapunov exponents were calculated for trajectories originating at each point using the "renormalization" algorithm.
The algorithm depends on several parameters: the renormalization time K, the integration time per renormalization T , and the integration timestep dt. The total time sampled to generate a single estimate of the global Lyapunov exponent is given by K T dt.
In order to test for ergodicity, the Lyapunov exponent was calculated for many short runs,

55
The resulting Lyapunov spectra are given in Table A. For each set of Lyapunov exponents, the Kaplan-Yorke fractal dimension (D KY ) may be directly calculated.
In general, the shorter simulation runs yielded a wide distribution of estimates for the exponents, primarily due to some initial conditions producing trajectories that remain stuck within the "metastable" slow dynamics on the rim of the teacup attractor. While the dis-  Short Runs (N=500) Long runs (N=100)

D Appropriateness of mean trait gradient dynamics
The gradient dynamics model used herein assumes that 1) the predator-prey dynamics depend solely on the mean of the trait distribution 2) the mean trait dynamics depend only on the current mean trait value and no high-order moments of the trait distribution 3) the additive genetic variance (V ) remains constant over long timescales. This results in the form of the dynamical equations used throughout the paper, In this section we note some of the underlying assumptions of these equations, and comment on their applicability.

D.1 Underlying assumptions of the mean trait evolution equation 75
First we consider the accuracy of assuming that the right-hand side of (A19) depends only on the current values of x(t), y(t),c(t). Following the derivation originally given by Lande,4 for an infinitely large population the mean fitness is given bȳ the dependence of each term on time has been suppressed from the notation; p(c) here represents a snapshot distribution of trait values in the prey population at given time, where the denominator represents the total prey population size. Taking the gradient of (A20) and inserting (A21), We now assume that the prey trait distribution has the form of a perturbed normal distribution, where the substitution V = h 2 V c has been performed; V thus represents the only the additive genetic variation (and not total variation) in the population. The first term in this 90 series is derived in Lande's original paper (which assumes κ 3 and all higher cumulants are zero); it is equivalent to the standard gradient dynamics model. 5,6 In general, the second "correction" term in this equation (which arises from a departure from Gaussianity in the trait distribution) cannot be solved analytically.
Comment on the assumption of constant additive genetic variance 95 As noted by previous investigators, even when an observed phenotypic distribution has strongly non-Gaussian form it is often possible to transform it into a distribution with Gaussian form due to the underlying additivity of the random processes that create the trait distribution. 4, 7 By the same token, a Gaussian with a first-order correction term can be used to describe nontrivial distributions that exhibit skew, given an appropriate coordinate 100 transform. 8,9 Additionally, it has observed that the moments of genetic distributions tend to remain fixed over time, justifying the assumption of holding V c and κ 3 constant in some cases.
However, even if the trait distribution has non-stationary V c , the additive genetic variance V may nonetheless remain stationary. Constant additive genetic variance is a common as-105 sumption in models in which gene selection is weak compared to selection on phenotypes. 10,11 The absence of net directional selection in the model varies. 14, 15 Examples include cases in which there is a constant degree of environmental heterogeneity, as well as fixed mating preferences and mutation rates among the population that serve to enforced a fixed degree of overall variation even in the absence of strong selection 115 forces. 10,11,16,17 For this reason, even as the prey population evolves, the additive genetic variation (as determined by the realized heritability observed in response to selection) may stay fixed. Because the dynamical equations used here do not specify the dynamics of reproduction or mutation, but rather just the fitness landscape for traits and the mean trait, the mean trait dynamics model may be most appropriate for populations that have been observed experimentally to maintain nearly-constant heritability values [18][19][20][21] However, for some cases-such as a prey population that fully speciates, directional selection that occurs for extended epochs, or mating and selection that deplete additive genetic variation--more advanced models of phenotypic evolution have been developed, which relax the assumption of constant genetic variation. 10, 16, 22 125 Numerical estimation of the potential contribution of non-Gaussianity in the trait distribution In order to determine the potential error in the dynamics introduced by neglecting the second term in (A23), the relative size of this term is computed ex post facto for simulations of the system generated for the case when κ 3 = 0 (the case used in the main text). The 130 relative contribution of the correction terms is determined using the ratio of the first and second terms in (A23) y, c,c) ) dc This quantity should be as small as possible for the gradient dynamics approximation to remain accurate. It has units of 1/κ 3 , and thus sets the maximum acceptable κ 3 in the trait distribution. In order to be conservative, V c may be set to the same value as V in Depending on the initial conditions, a distribution resulting from an exhaustive numerical 160 simulation of many individuals may never reach a κ 3 as high as the upper bound presented above. However, if circumstances arise in which the the additional term in (A23) has an overall effect on the long-term dynamics, this effect would likely be to stabilize the dynamics the the system and lead to transient chaos (in which the dynamics eventually exit the chaotic attractor and seek a stable equilibrium or limit cycle). This is because the negative sign of 165 additional term in (A23) dampens the dynamics when the fitness landscape has sharp peaks and valley. Interestingly, a constant fitness landscape (r(x, y, c,c) = r 0 (x, y)) causes the second term in (A23) to equal zero.

D.2 Analysis of the direct dynamics of the full trait distribution
We can further assess the accuracy of the gradient dynamics model by eliminating (A19) 170 and instead writing the predator-prey system in terms of a full integro-differential equation that depends on the full trait distribution. If x(t, c) denotes the density of prey with trait value c, then the mean trait value becomes where ∫ x(c)dc represents the total prey density across all trait values. In this case, the system becomes a system of two coupled ordinary differential equations, the first of which depends on the integral (A24) Full numerical solution of this system of equations is difficult for long periods due to the requirement that the integral term (A24) be evaluated at every timepoint that the numerical 175 integrator computes the numerical derivative constituting the right hand side of (A25). For this reason, small errors in the computation of the mean trait value accumulate quickly, leading the numerical integrator to converge prohibitively slowly to allow direct comparison of numerical solutions to (A17),(A18),(A19) to those of (A25),(A26) over long timescales.
A simple method of determining the accuracy of the gradient dynamics approximation 180 instead relies on the observing that both the integro-differential equation, and the gradient dynamics approximation, cast the predator prey coevolution problem in terms of first-order dynamics. For this reason, comparing the time evolution of the velocity field (as a function of x, y, andc) under the two formulations can be used to compare how close the dynamics of the two models would be expected to be, even in the absence of full numerical integration 185 (in which small errors in the velocity fields can accumulate over time). We thus use the following algorithm to determine the effect having a distribution of trait values, x(t, c), has on the dynamics relative to the gradient dynamics approximation: 1. Using the gradient dynamics model (the three-dimensional system (A17),(A18),(A19)), we simulate a long trajectory in the system that resides on the chaotic attractor, Because the predator dynamics depend only on the total prey density ( ∫ x(t, c)dc) then the trait-averaged predator density has the simple formȳ id (t) =ẏ id (t). This produces 210 a time series of trait-averaged velocity values under the integro-differential equation, v id (t) ≡ (x id (t),ȳ id (t)).
6. The two time series of velocity values under the gradient dynamics (v gd (t), from step 2) and integro-differential (v id (t), from step 5) models are compared using the average squared error, which penalizes cases in which the instantaneous velocity greatly differs 215 between the two models. A large squared error suggests that the full integro-differential equation would produce different dynamics at many points on the attractor, resulting in qualitatively different dynamics between the two models that would grow larger as  curve). Figure D shows the similarity score between the velocity time series of the gradient dynamics model and the velocity time series of the integro-differential model, computed as a function of the variance parameter used in the ID model. Similarity is calculated as one minus the mean squared difference between the normalized velocity time series for the two models, with each calculated for a very long trajectory residing on the chaotic attractor.

235
The peak in similarity appears near the point at which the variances of the two models are roughly the same, suggesting general agreement between the two models. More precise peak determination could be achieved with greater computational resources.