A phase transition induces chaos in a predator-prey ecosystem with a dynamic fitness landscape
In many natural ecosystems , at least one constituent species evolves quickly enough relative to its population growth that the two effects become interdependent ., This phenomenon can occur when selection forces are tied to such sudden environmental effects as algal blooms or flooding 1 , or it can arise from more subtle , population-level effects such as overcrowding or resource depletion 2 ., Analysis of such interactions within a unified theory of “eco-evolutionary dynamics” has been applied to a wide range of systems—from bacteria-phage interactions to bighorn sheep 3—by describing population fluctuations in terms of the feedback between demographic change and natural selection 4 ., The resulting theoretical models relate the fitness landscape ( or fitness function ) to population-level observables such as the population growth rate and the mean value of an adapting phenotypic trait ( such as horn length , cell wall thickness , etc ) ., The fitness landscape may have an arbitrarily complex topology , as it can depend on myriad factors ranging from environmental variability 5 , 6 , to inter- and intraspecific competition 7 , 8 , to resource depletion 9 ., However , these complex landscapes can be broadly classified according to whether they result in stabilizing or disruptive selection ., In the former , the landscape may possess a single , global maximum that causes the population of individuals to evolve towards a state in which most individuals have trait values at or near this maximum 10 ., Conversely , in disruptive selection , the fitness landscape may contain multiple local maxima , in which case the population could have a wide distribution of trait values and occupy multiple distinct niches 11 ., In eco-evolutionary models , the shape of the fitness landscape may itself depend on the population densities of the interacting species it describes ., Specifically , the concept that the presence of competition can lead a single-peaked fitness landscape to spontaneously develop additional peaks originates in the context of “competitive speciation” first proposed by Rosenzweig 12 ., This is formalized in genetic models in which sympatric speciation is driven by competitive pressures rather than geographic isolation 13 ., Competition-induced disruptive selection has been observed in natural populations of stickleback fish 14 , microbial communities 15 , and fruit flies 16 , 17 ., Here , we model eco-evolutionary dynamics of a predator-prey system based on first-order “gradient dynamics” 10 , 18 , a class of models that explicitly define the fitness in terms of the population growth rate r , which is taken to depend only on the mean value of the trait across the entire population , c ¯ 19 ., Despite this simplification , gradient dynamics models display rich behavior that can account for a wide range of effects observed in experimental systems—in particular , recent work by Cortez and colleagues has shown that these models can result in irregular cycles and dynamical bifurcations that depend on the standing genetic variation present in a population 20 , 21 ., In our model , gradient dynamics cause the prey fitness landscape to change as a result of predation , and we find that the resulting dynamical system exhibits chaotic dynamics ., Chaos is only possible in systems in which three or more dependent dynamical variables vary in time 22 , and previously it has been observed in predator-prey systems comprising three or more mutually interdependent species , or in which an external environmental variable ( such as seasonal variation or generic noise ) is included in the dynamics 23 , 24 ., Here we show that evolution of just one species in a two-species ecosystem is sufficient to drive the ecosystem into chaos ., Moreover , we find that chaos is driven by a density-dependent change of the fitness landscape from a stabilizing to disruptive state , and that this transition has hysteretic behavior with mathematical properties that are strongly reminiscent of a first-order phase transition in a thermodynamical system ., The resulting dynamics display intermittent properties typically associated with ecosystems poised at the “edge of chaos , ” which we suggest has implications for the study of ecological stability and speciation ., Adapting the notation and formulation used by Cortez ( 2016 ) 21 , we use a two-species competition model with an additional dynamical variable introduced to account for a prey trait on which natural selection may act ., The most general fitness function for the prey , r , accounts for density-dependent selection on a prey trait c ,, r ( x , y , c ¯ , c ) ≡ G ( x , c , c ¯ ) - D ( c , c ¯ ) - f ( x , y ) , ( 1 ), where x = x ( t ) is the time-dependent prey density , y = y ( t ) is the time-dependent predator density , c is a trait value for an individual in the prey population , and c ¯ = c ¯ ( t ) is the mean value of the trait across the entire prey population at time t ., r comprises a density-dependent birth rate G , a density-independent death rate D , and a predator-prey interaction term f , which for simplicity is assumed to depend on neither c nor c ¯ ., Thus the trait under selection in our model is not an explicit predator avoidance trait such as camouflage , but rather an endogenous advancement ( i . e . , improved fecundity , faster development , or reduced mortality ) that affects the prey’s ability to exploit resources in its environment , even in the absence of predation ., The continuous-time “gradient dynamics” model that we study interprets the fitness r as the growth rate of the prey: 19 , 25, x ˙ = x r ( x , y , c ¯ , c ) | c → c ¯ ( 2 ), y ˙ = y ( f ( x , y ) - D ˜ ( y ) ) ( 3 ), c ¯ ˙ = V ∂ r ( x , y , c ¯ , c ) ∂ c | c → c ¯ ., ( 4 ) Eq ( 2 ) is evaluated with all individual trait values c set to the mean value c ¯ because the total prey population density is assumed to change based on the fitness function , which in turn depends on the population-averaged value of the prey trait c ¯ 21 ., The timescale of the dynamics in c ¯ are set by V , which is interpreted as the additive genetic variance of the trait 10 ., While Eq ( 2 ) depends only on the mean trait value c ¯ , the full distribution of individual trait values c present in a real-world population may change over time as the relative frequencies of various phenotypes change ., In principle , additional differential equations of the form of Eq ( 4 ) could be added to account for higher moments of the distribution of c across an ensemble of individuals , allowing the gradient dynamics model to be straightforwardly extended to model a trait’s full distribution rather than just the population mean ., However , here we focus on the case where the prey density dynamics x ˙ depend only on the mean trait value to first order , and we do not include differential equations for higher-order moments of the prey trait value distribution ., The use of a single Eq ( 4 ) to describe the full dynamics of the trait distribution represents an approximation that is exact only when the phenotypic trait distribution stays nearly symmetric and the prey population maintains a constant standing genetic variation V 10 ., However , V may remain fixed even if the phenotypic variance changes , a property that is observed phenomenologically in experimental systems , and which may be explained by time-dependent heritability , breeding effects , mutation , or other transmission effects not explicitly modeled here 26–29 ., More broadly , this assumption may imply that gene selection is weak compared to phenotype selection 30 , 31 ., S1D Appendix further describes the circumstances under which V remains fixed , and also provides a first-order estimate of the magnitude of error introduced by ignoring higher-order effects ( such as skewness ) in the trait distribution ., The results suggest that these effects are small for the parameter values ( and resulting range of x and y values ) used here , due in part to limitations on the maximum skewness that a hypothetical trait distribution can achieve on the fitness landscapes studied here ., In S1D Appendix , we also compare the results presented below to an equivalent model in which a full trait distribution is present , in which case Eq ( 2 ) becomes a full integro-differential equation involving averages of the trait value over the entire prey population ., Detailed numerical study of this integro-differential equation is computationally prohibitive for the long timescales studied here , but direct comparison of the contributions of various terms in the velocity field suggests general accuracy of the gradient dynamics model for the fitness landscapes and conditions we study here ., However , in general the appropriateness of the gradient dynamics model should be checked whenever using Eq ( 4 ) with an arbitrary fitness function ., Fig 1A shows a schematic summarizing the gradient dynamics model , and noting the primary assumptions underlying this formulation ., Next , we choose functional forms for f , G , D , and D ˜ in Eqs ( 2 ) and ( 3 ) ., We start with the assumption that , for fixed values of the trait c an d its mean c ¯ , the population dynamics should have the form of a typical predator-prey system in the absence of evolutionary effects ., Because the predator dynamics are not directly affected by evolutionary dynamics , we choose a simple form for predator growth consisting of a fixed death rate and a standard Holling Type II birth rate , 32, f ( x , y ) = a 2 x y 1 + b 2 x ( 5 ), D ˜ ( y ) = d 2 ( 6 ), The predator birth rate f saturates at large values of the prey density , which is more realistic than the standard Lotka-Volterra competition term xy in cases where the prey density is large or fluctuating 22 ., A saturating interaction term ensures that solutions of the system remain bounded for a wider range of parameter values , a necessity for realistic models of long-term interactions 33 ., For the prey net growth rate ( Eq ( 1 ) , the fitness ) in the absence of the predator , we use the following functional forms ,, G ( x , c ¯ , c ) = a 1 c ¯ 1 + b 1 c ¯ ( 1 - k 1 x ( c - c ¯ ) ) ( 7 ), D ( c , c ¯ ) = d 1 ( 1 - k 2 ( c 2 - c ¯ 2 ) + k 4 ( c 4 - c ¯ 4 ) ) ., ( 8 ), The first term in Eq ( 7 ) specifies that the prey population density growth rate r | c → c ¯ depends only on a primary saturating contribution of the mean trait to the birth rate G . In other models a similar effect is achieved by modifying the mean trait evolution Eq ( 4 ) , such that extremal values of the trait are disadvantaged 21; alternative coupling methods based on exponential saturation would be expected to yield similar qualitative results 19 ., However , the additional series terms in Eqs ( 7 ) and ( 8 ) ensure that the any individual’s fitness r may differ from the rest of the population depending on the difference between its trait value c and the population mean c ¯ ., Because the functional form of this difference is unknown , its contribution expressed as second-order truncation of the series difference of the form r ( c , c ¯ ) = r ˜ | c → 0 + ( r ˜ ( c ) - r ˜ ( c ) | c → c ¯ ) ( where r ˜ represents an unscaled fitness function ) ., This ensures that when c ˙ = 0 or c = c ¯ , the system reduces to a standard prey model with a Holling Type II increase in birth rate in response to increasing mean trait value 25 ., In the results reported below , we observe that all dynamical variables remain bounded as long as parameter values are chosen such that the predator density does not equilibrate to zero ., This is a direct consequence of our use of saturating Holling Type II functional forms in Eqs ( 7 ) and ( 8 ) , which prevent the fitness landscape from increasing without bound at large c , c ¯ and also ensure that the predator and prey densities do not jointly diverge ., That the dynamics should stay bounded due to saturating terms is justified by empirical studies of predator-prey systems 34 , 35; moreover , other saturating functional forms are expected to yield similar results if equivalent parameter values are chosen 33 , 36 ., The nonlinear dependence of the mortality rate Eq ( 8 ) on the trait is based on mechanistic models of mortality with individual variation 19 , 37 , 38 ., The specific choice of a quartic in Eq ( 8 ) allows the fitness function r to have a varying number of real roots and local maxima in the domain c , c ¯ > 0 , affording the system dynamical freedom not typically possible in predator prey models with constant or linear prey fitness—in particular , for different values of k2 , k4 the fitness landscape can have a single optimal phenotype , multiple optimal phenotypes , or no optimal intermediate values ., Because any even , continuous form for the fitness landscape can be approximated using a finite number of terms in its Taylor series around c = 0 , our choice of a quartic form simply constitutes truncation of this expansion at the quartic order in order to include the simplest case in which the fitness function admits multiple local maxima—for this reason , a quartic will always represent the leading-order series expansion of a fitness landscape with multiple local maxima ., Below , we observe numerically that ∣ c - c ¯ ∣ < 1 , ex post facto justifying truncation of the higher order terms in this series expansion ., However , if the trait value c was strictly bounded to only take non-zero values on a finite interval ( as opposed to the entire real line ) , then a second-order , quadratic fitness landscape would be sufficient to admit multiple local maxima ( at the edges of the interval ) 14 ., However , the choice here of an unbounded trait value c avoids creating boundary effects , and it has little consequence due to the steep decay of the quartic function at large values of |c| , which effectively confines the possible values of c ¯ accessible by the system ., In physics , similar reasons—unbounded domains , multiple local optima , and continuity—typically justify the use of quartic free energy functions in minimal models of systems exhibiting multiple energetic optima , such as the Ginzberg-Landau free energy used in models of superconducting phase transitions 39 ., We note that the birth rate Eq ( 7 ) contributes a density-dependent term to the fitness function even in the absence of predation ( y = 0 ) 21 ., Unlike the death rate function , the effect of the individual trait value on this term is directional: the sign of c - c ¯ determines whether birth rates increase or decrease ., As the population density x increases , the effect of these directional effects is amplified , consistent with the observed effect of intraspecific competition and crowding in experimental studies of evolution 40 , 41 ., The chaotic dynamics reported below arise from this density-dependent term because the term prevents the Jacobian of the system ( 2 ) , ( 3 ) and ( 4 ) from having a row and column with all zeros away from the diagonal; in this case , the prey trait ( and thus evolutionary dynamics ) would be uncoupled from the rest of the system , and would thus relax to a stable equilibrium ( as is necessary for a first-order single-variable equation ) ., In that case , c ¯ would essentially remain fixed and the predator-prey dynamics would become two-dimensional in x and y , precluding chaos ., For similar reasons , density-dependent selection has been found to be necessary for chaos in some discrete-time evolutionary models , for which chaotic dynamics require a certain minimum degree of association between the fitness and the trait frequencies 42 ., Inserting Eqs ( 5 ) , ( 7 ) and ( 8 ) , into Eq ( 1 ) results in a final fitness function of the form, r ( x , y , c ¯ , c ) = a 1 c ¯ 1 + b 1 c ¯ ( 1 - k 1 x ( c - c ¯ ) ) - d 1 ( 1 - k 2 ( c 2 - c ¯ 2 ) + k 4 ( c 4 - c ¯ 4 ) ) - a 2 x y 1 + b 2 x ., ( 9 ), This fitness landscape is shown in Fig 1B , for typical parameter values and predator and prey densities used in the numerical results below ., Depending on the current predator and prey densities , the local maximum of the system can appear in two different locations , which directly affects the dynamics described in the next section ., Inserting Eq ( 9 ) into Eqs ( 2 ) , ( 3 ) and ( 4 ) results in a final form for the dynamical equations ,, x ˙ = x ( a 1 c ¯ 1 + b 1 c ¯ - a 2 y 1 + b 2 x - d 1 ) ( 10 ), y ˙ = y ( y a a 2 x 1 + b 2 x - d 2 ) ( 11 ), c ¯ ˙ = c ¯ V ( ( 2 k 2 d 1 ) - ( 4 k 4 d 1 ) c ¯ 2 - ( a 1 k 1 ) x 1 + b 1 c ¯ ) ., ( 12 ), Due to the Holling coupling terms , the form of these equations qualitatively resembles models of vertical , tritrophic food webs—the mean trait value c ¯ affects the growth rate of the prey , which in turn affects the growth rate of the predator 24 , 32 , 43 ., The coupling parameter ya introduces asymmetry into the competition when ya ≠ 1; however , it essentially acts as a scale factor that only affects the amplitude of the y cycles and equilibria rather than the dynamics ., Additionally , because the predator-prey interaction term Eq ( 5 ) is unaffected by the trait , our model contains no triple-product x y c ¯ interaction terms , which typically stabilize the dynamics ., For our analysis of the system ( 10 ) , ( 11 ) and ( 12 ) , we first consider the case where evolution proceeds very slowly relative to population dynamics ., In the case of both no evolution ( V = 0 ) and no predation ( y = 0 ) , the prey growth Eq ( 10 ) advances along the one-dimensional nullcline y ˙ , c ¯ ˙ = 0 , y = 0 ., Depending on whether the fixed mean trait value c ¯ exceeds a critical value ( c ¯ † ≡ d 1 / ( a 1 - b 1 d 1 ) ) , the prey density will either grow exponentially ( c ¯ > c ¯ † ) or collapse exponentially ( c ¯ < c ¯ † ) because the constant c ¯ remains too low to sustain the prey population in the absence of evolutionary adaptation ., The requirement that c ¯ > c ¯ † carries over to the case where a predator is added to the system but evolutionary dynamics remain fixed , corresponding to a two dimensional system advancing along the two-dimensional nullcline c ¯ ˙ = 0 ., In this case , as long as c ¯ > c ¯ † , the prey density can exhibit continuous growth or cycling depending in the relative magnitudes of the various parameters in Eqs ( 10 ) and ( 11 ) ., The appearance and disappearance of these cycles is determined by a series of bifurcations that depends on the values of c ¯ and b1 , b2 relative to the remaining parameters a1 , a2 , d1 , d2 ( S1A Appendix ) ., In the full three-variable system ( 10 ) , ( 11 ) and ( 12 ) , c ¯ passes through a range of values as time progresses , resulting in more complex dynamics than those observed in the two-dimensional case ., For very small values of V , the evolutionary dynamics c ¯ ˙ are slow enough that the system approaches the equilibrium predicted by the two-variable model with c ¯ constant ., The predator and prey densities initially grow , but the prey trait value does not change fast enough for the prey population growth to sustain—eventually resulting in extinction of both the predator and prey ., However , if V takes a slightly larger value , so that the mean trait value can gradually change with a growing prey population density ( due to the density-dependent term in Eq ( 10 ) ) , then the population dynamics begin to display regular cycling with fixed frequencies and amplitudes ( Fig 2A , top ) ., This corresponds to a case where the evolutionary dynamics are slow compared to the ecological dynamics , but not so slow as to be completely negligible ., Finally , when V is the same order of magnitude as the parameters governing the ecological dynamics , the irregular cycles become fully chaotic , with both amplitudes and frequencies that vary widely over even relatively short time intervals ( Fig 2A , bottom ) ., Typically , the large V case would correspond to circumstances in which the prey population develops a large standing genetic variation 10 , 44 ., That the dynamics are chaotic , rather than quasi-periodic , is suggested by the presence of multiple broad , unevenly-spaced peaks in the power spectrum 45 ( Figure A in S1E Appendix ) , as well as by numerical tabulation of the Lyapunov spectrum ( described further below ) ., Due to the hierarchical coupling of Eqs ( 10 ) , ( 11 ) and ( 12 ) , when plotted in three-dimensions the chaotic dynamics settle onto a strange attractor that resembles the “teacup” attractor found in models of tritrophic food webs 24 , 46 ( Fig 2B ) ., Poincare sections though various planes of the teacup all appear linear , suggesting that the strange attractor is effectively two-dimensional—consistent with pairings of timescales associated with different dynamical variables at different points in the process ( Figure B in S1E Appendix ) ., In the “rim” of the teacup , the predator density changes slowly relative to the prey density and mean trait value ., This is visible in a projection of the attractor into the x - c ¯ plane ( Fig 2B , bottom inset ) ., However , in the “handle” of the teacup , the mean trait value varies slowly relative to the ecological dynamics ( c ¯ ˙ ≈ 0 ) , resulting in dynamics that qualitatively resemble the two-dimensional “reduced” system described above for various fixed values of c ¯ ( Fig 2B , top inset ) ., The structure of the attractor suggests that the prey alternately enters periods of evolutionary change and periods of competition with the predator ., A closer inspection of a typical transition reveals that this “two timescale” dynamical separation is responsible for the appearance of chaos in the system ( Fig 3A ) ., As the system explores configuration space , it reaches a metastable configuration corresponding to a high mean trait value c ¯ , which causes the prey density to nearly equilibrate to a low density due to the negative density-dependent term in Eq ( 10 ) ., During this period ( the “rim” of the teacup ) , the predator density gradually declines due to the lack of prey ., However , once the predator density becomes sufficiently small , the prey population undergoes a sudden population increase , which triggers a period of rapid cycling in the system ( the “handle” of the teacup attractor ) ., During this time , the predator density continuously increases , causing an equivalent decrease in the prey density that resets the cycle to the metastable state ., The sudden increase in the prey population at low predator densities can be understood from how the fitness function r ( from Eq ( 9 ) ) changes over time ., Fig 3B shows a kymograph of the log-scaled fitness Eq ( 9 ) as a function of individual trait values c , across each timepoint and corresponding set of ( x , y , c ¯ ) values given in panel A . Overlaid on this time-dependent fitness landscape are curves indicating the instantaneous location of the local maximum ( black ) and minimum ( white ) ., By comparing panels A and B , it is apparent that the mean trait value during the “metastable” period of the dynamics stays near the local maximum of the fitness function , which barely varies as the predator density y changes ., However , when y ( t ) ≈ 0 . 25 , the fitness function changes so that the local minimum and local maximum merge and disappear from the system , leading to a new maximum spontaneously appearing at c = 0 ., Because V is large enough ( for these parameters ) that the gradient dynamics occur over timescales comparable to the competition dynamics , the system tends to move rapidly towards this new maximum in the fitness landscape , resulting in rapidly-changing dynamics in x and c ¯ ., Importantly , because of the symmetric coupling of the prey fitness landscape r to the prey density x , this rapid motion resets the fitness landscape so that the maximum once again occurs at the original value , resulting in a period of rapid cycling ., The fitness landscape at two representative timepoints in the dynamics is shown in Fig 3C ., That the maxima in the fitness Function ( 9 ) suddenly change locations with continuous variation in x , y is a direct consequence of the use of a high-order ( here , quartic ) polynomial in c to describe the fitness landscape ., The quartic represents the simplest analytic function that admits more than one local maxima in its domain , and the number of local maxima is governed by the relative signs of the coefficients of the ( c 2 - c ¯ 2 ) and ( c 4 - c ¯ 4 ) terms in Eq ( 9 ) , which change when the system enters the rapid cycling portion of the chaotic dynamics at t = 500 in Fig 3A ., This transition marks the mean prey trait switching from being drawn ( via the gradient dynamics ) to a single fitness peak at an intermediate value of the trait ceq ≈ 0 . 707 to being drawn instead to one of two peaks: the existing peak , or a new peak at the origin ., Thus the metastable period of the dynamics corresponds to a period of stabilizing selection: if the fitness landscape were frozen in time during this period , then an ensemble of prey would all evolve to a single intermediate trait value corresponding to the location of the global maximum ., Conversely , if the fitness landscape were held fixed in the multipeaked form it develops during a period of rapid cycling , given sufficient time an ensemble of prey would evolve towards subpopulations with trait values at the location of each local fitness maximum—representing disruptive selection ., That the fitness landscape does not remain fixed for extended durations in either a stabilizing or disruptive state—but rather switches between the two states due to the prey density-dependent term in Eq ( 9 ) — underlies the onset of chaotic cycling in the model ., Density-dependent feedback similarly served to induce chaos in many early discrete-time ecosystem models 23 ., However , the “two timescale” form of the chaotic dynamics and strange attractor here is a direct result of reversible transitions between stabilizing and disruptive selection ., If the assumptions underlying the gradient dynamics model do not strictly hold—if the additive genetic variance V slowly varies via an additional dynamical equation , or if the initial conditions are such that significant skewness would be expected to persist in the phenotypic distribution , then the chaotic dynamics studied here would be transient rather than indefinite ., While the general stability analysis shown above ( and in the S1 Appendix ) would still hold , additional dynamical equations for V or for high-order moments of the trait distribution would introduce additional constraints on the values of the parameters , which would ( in general ) increase the opportunities for the dynamics to become unstable and lead to diverging predator or prey densities ., However , in some cases these additional effects may actually serve to stabilize the system against both chaos and divergence ., For example , if additional series terms were included in Eq ( 8 ) such that the dependence of mortality rate on c ¯ and c had an upper asymptote 25 , then c ¯ ˙ = 0 would be true for a larger range of parameter values—resulting in the dynamical system remaining planar for a larger range of initial conditions and parameter values , precluding chaos ., The transition between stabilizing and disruptive selection that occurs when the system enters a period of chaotic cycling is strongly reminiscent of a first-order phase transition ., Many physical systems can be described in terms of a free energy landscape , the negative gradient of which determines the forces acting on the system ., Minima of the free energy landscape correspond to equilibrium points of the system , which the dynamical variables will approach with first-order dynamics in an overdamped limit ., When a physical system undergoes a phase transition—a qualitative change in its properties as a single “control” parameter , an externally-manipulable variable such as temperature , is smoothly varied—the transition can be understood in terms of how the control parameter changes the shape of the free energy landscape ., The Landau free energy model represents the simplest mathematical framework for studying such phase transitions: a one-dimensional free energy landscape is defined as a function of the control parameter and an additional independent variable , the “order parameter , ” a derived quantity ( such as particle density or net magnetization ) with respect to which the free energy can have local minima or maxima ., In a first-order phase transition in the Landau model , as the control parameter monotonically changes the relative depth of a local minimum at the origin decreases , until a new local minimum spontaneously appears at a fixed nonzero value of the order parameter—resulting in dynamics that suddenly move towards the new minimum , creating discontinuities in thermodynamic properties of the system such as the entropy 47 ., First-order phase transitions are universal physical models , which have been used to describe a broad range of processes spanning from superconductor breakdown 48 to primordial black hole formation in the early universe 49 ., In the predator-prey model with prey evolution , the fitness function is analogous to the free energy , with the individual trait value c serving as the “order parameter” for the system ., The control parameter for the transition is the prey density , x , which directly couples into the dynamics via the density-dependent term in Eq ( 7 ) ., Because the fitness consists of a linear combination of this term in Eq ( 7 ) and a quartic landscape Eq ( 8 ) , the changing prey density “tilts” the landscape and provokes the appearance of the additional , disruptive peak visible in Fig 3C ., The appearance and disappearance of local maxima as the system switches between stabilizing and disruptive selection is thus analogous to a first-order phase transition , with chaotic dynamics being a consequence of repeated increases and decreases of the control parameter x above and below the critical prey densities x* , x** at which the phase transition occurs ., Similar chaotic dynamics emerge from repeated first-order phase transitions in networks of coupled oscillators , which may alternate between synchronized and incoherent states that resemble the “metastable” and “rapid cycling” portions of the predator-prey dynamics 50 ., The analogy between a first-order phase transition and the onset of disruptive selection can be used to study the chaotic dynamics in terms of dynamical hysteresis , a defining feature of such phase transitions 47 ., For different values of x , the three equilibria corresponding to the locations of the local minima and maxima of the fitness landscape , ceq , can be calculated from the roots of the cubic in Eq ( 12 ) ., The resulting plots of ceq vs x in Fig 4 are generated by solving for the roots in the limit of fast prey equilibration , c ¯ → c e q , which holds in the vicinity of the equilibria ( S1B Appendix ) ., The entry into the transient chaotic cycling occurs when x increases gradually and shifts ceq with it; x eventually attains a critical value x* ( x* ≈ 0 . 45 for the parameters used in the figures ) , causing ceq to jump from its first critical value c* to the origin ( the red “forward” branch in Fig 4 ) ., This jump causes rapid re-equilibration of c ¯ ( t ) , resulting in the rapid entry into cycling observable in Fig 3A ., However , x cannot increase indefinitely due to predation; rather , it decreases until it reaches a second critical value x** , at which point ceq jumps back from the origin to a positive value ( the blue “return” branch in Fig 4; x** = 0 . 192 for these parameter values ) ., This second critical point marks the return to the metastable dynamics in Fig 3A ., This asymmetry in the forward and backwards dynamics of x lead to dynamical time-irreversibility ( hysteresis ) and the jagged , sawtooth-like cycles visible in the dynamics of the full system ., Because the second jump in ceq is steeper , the parts of the trajectories associated with the “return” transition in Fig 3A appear steeper ., Additionally , the maximum value obtained by c ¯ ( t ) anywhere on the attractor , c e q m a x , is determined by the limiting value o
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 .
Evolution is usually thought to occur very gradually , taking millennia or longer in order to appreciably affect a species survival mechanisms ., Conversely , demographic shifts due to predator invasion or environmental change can occur relatively quickly , creating abrupt and lasting effects on a species survival ., However , recent studies of ecosystems ranging from the microbiome to oceanic predators have suggested that evolutionary and ecological processes can often occur over comparable timescales—necessitating that the two be addressed within a single , unified theoretical framework ., Here , we show that when evolutionary effects are added to a minimal model of two competing species , the resulting ecosystem displays erratic and chaotic dynamics not typically observed in such systems ., We then show that these chaotic dynamics arise from a subtle analogy between the evolutionary concept of fitness , and the concept of the free energy in thermodynamical systems ., This analogy proves useful for understanding quantitatively how the concept of a changing fitness landscape can confer robustness to an ecosystem , as well as how unusual effects such as history-dependence can be important in complex real-world ecosystems ., Our results predict a potential signature of a chaotic past in the distribution of timescales over which new species can emerge during the competitive dynamics , a potential waypoint for future experimental work in closed ecosystems with controlled fitness landscapes .
ecology and environmental sciences, predator-prey dynamics, population dynamics, systems science, mathematics, population biology, thermodynamics, computer and information sciences, ecosystems, dynamical systems, free energy, community ecology, physics, population metrics, ecology, predation, natural selection, trophic interactions, biology and life sciences, physical sciences, population density, evolutionary biology, evolutionary processes