We used a standard three-species model with a linear trophic link (N1: resource, N2: prey, and N3: predator). Nonlinear functional responses, namely Holling’s type II, were assumed for the prey grazing on the resource and the predator preying on the prey. Additional nonlinearity was implemented for intraspecific competition at the basal level of resources, with an autotroph population. The dynamics of the abundance of the three species on the time scale T are
$${\frac{d{N}_{1}}{dT}=rN}_{1}\left(1-\frac{{N}_{1}}{K}\right)-\frac{{A}_{2}{N}_{1}{N}_{2}}{1+{A}_{2}{h}_{2}{N}_{1}}$$
1
$$\frac{d{N}_{2}}{dT}={N}_{2}\left(\frac{{{c}_{2}A}_{2}{N}_{1}}{1+{A}_{2}{h}_{2}{N}_{1}}-{D}_{2}\right)-\frac{{A}_{3}{N}_{2}{N}_{3}}{1+{A}_{3}{h}_{3}{N}_{2}}$$
2
$$\frac{d{N}_{3}}{dT}={N}_{3}\left(\frac{{{c}_{3}A}_{3}{N}_{2}}{1+{A}_{3}{h}_{3}{N}_{2}}-{D}_{3}\right)$$
3
,
where r and K are the intrinsic rate of natural increase and the carrying capacity of the resource, A is the consumption efficiency; h is the inverse maximum consumption rate; c is the conversion coefficient; and D is the mortality. All subscripts denote the different trophic levels.
We simplified the system using Eq. (1)–(3), according to Hastings and Powell (1991), by standardizing the biomass for all species by the resource carrying capacity K (\({n}_{i}={N}_{i}/K\), i = 1,2 and 3) and by rescaling time with the intrinsic rate of natural increase r of the resource (\(t=rT\)),
$${\frac{d{n}_{1}}{dt}=n}_{1}\left(1-{n}_{1}\right)-\frac{{a}_{2}{n}_{1}{n}_{2}}{1+{a}_{2}{h}_{2}{n}_{1}}$$
4
$$\frac{d{n}_{2}}{dt}={n}_{2}\left(\frac{{{c}_{2}a}_{2}{n}_{1}}{1+{a}_{2}{h}_{2}{n}_{1}}-{d}_{2}\right)-\frac{{a}_{3}{n}_{2}{n}_{3}}{1+{a}_{3}{h}_{3}{n}_{2}}$$
5
$$\frac{d{n}_{3}}{dt}={n}_{3}\left(\frac{{{c}_{3}a}_{3}{n}_{2}}{1+{a}_{3}{h}_{3}{n}_{2}}-{d}_{3}\right)$$
6
,
where \({a}_{j}={A}_{j}/r\) and \({d}_{j}={D}_{j}/r\) (j = 2 and 3). The original timescale was multiplied by r so that the intrinsic growth rate on the standardized time became 1, and the biomass of the three species was scaled according to the carrying capacity K of the producer (species 1), ensuring the minimum number of model parameters.
In the above tri-trophic system, Eq. (1)–(3) or Eqs. (4)–(6), is known to have heteroclinic cycles, exhibiting chaotic behavior within a specific range of model parameters (Fig. 1; Hastings and Powell 1991; McCann and Yodzis 1994).
We introduced an anti-predator defence trait for the prey and focused on the effect of its evolutionary change on the ecological dynamics of the community. We excluded any predator traits that would facilitate predation efficiency because the evolutionary rate of the predator was likely much slower than that of the prey, allowing the coevolutionary process to be negligible. This was followed by the prey having a much higher metabolic rate and a much shorter generation time than the predator. We observed no altered results in terms of eco-evolutionary dynamics with an alternative model that included predator evolution (simulation data not shown). Therefore, the results of the present study may not hold for cases in which coevolutionary instability significantly affects eco-evolutionary dynamics.
The quantitative trait z of the prey is assumed to be an anti-predator trait that determines vulnerability to predation by the predator. The following predation efficiency function was used in place of a3 in Equations (5) and (6).
$${a}_{3}\left(z\right)={a}_{max}\text{e}\text{x}\text{p}\left\{-{\left(\frac{z}{s}\right)}^{2\gamma }\right\}$$
7
,
where amax is the maximum predation efficiency when the trait is at an optimal state 0, s measures the scale of the trait and, and γ is the slope parameter of the association between the scaled trait value and predation efficiency (γ is set to unity for all the simulations presented in this article).
The anti-predator trait entailed a cost on fitness as long as it increased the tolerance against predation. The fitness component of the prey that indicated a negative selection of the prey trait was proportional to the following Gaussian cost function: \(\text{e}\text{x}\text{p}\left(-\omega {z}^{2}/2\right)\), in which ω is the fitness cost of the anti-predator trait and was measured in the context of predation not being directly related. We excluded the demographic effect of the cost of the trait on the prey, because preliminary simulations have demonstrated that such a contribution played a lesser role in comparison to the effect of prey trait evolution and is likely to have almost no effect on the behavior of the whole system with the ecological parameter values set in this study.
The rate of evolution of a single trait per unit time is generally the selection gradient, which is the slope of the log mean fitness \(\stackrel{-}{W}\) against the mean trait value \(\stackrel{-}{z}\), multiplied by the additive genetic variance of the trait G, that is, \(\text{G}dln\stackrel{-}{W}/d\stackrel{-}{z}\) (Lande 1982). Applying the selection gradient approximated as \({\left(dlnW/dz\right)}_{z=\stackrel{-}{z}}\) (Abrams et al 1993), the fitness of the prey is \({W}_{2}=\text{e}\text{x}\text{p}\left(\frac{d{n}_{2}}{{n}_{2}dt}-\frac{\omega }{2}{z}^{2}\right)\), which, along with Eq. (5), specifies the evolutionary rate as:
$$\frac{dz}{dt}=-G\left\{\omega z+\frac{{a}_{3}^{\text{'}}{n}_{3}}{{\left(1+{a}_{3}{h}_{3}{n}_{2}\right)}^{2}}\right\}$$
8
,
where \({a}_{3}^{\text{'}}\) is the differential coefficient of \({a}_{3}\) by z: \({a}_{3}^{\text{'}}=-{a}_{3}\frac{2\gamma }{s}{\left(\frac{z}{s}\right)}^{2\gamma -1}\).
Numerical procedures
Numerical evaluations of the community model and trait evolution based on Eqs. (4)–(6) and (8), was undertaken to examine the effects of trait evolution on the dynamical properties of the trajectories and periodicity of demographic fluctuations. We focused on how the evolvability or the evolution limit imposed by the fitness cost and evolutionary rate could affect the stability and long-term periodicity of community dynamics.
Parameterization
The set of baseline ecological parameter values used in the simulations, i.e., a2 = 2, amax=0.1, c2 = c3 = 1, d2 = 0.2 (varying from 0.1 to 0.3) and d3 = 0.01, is known to generate chaotic dynamics in the linear tri-trophic food chain model (Hastings and Powell 1991).
We examined the interaction effect between trait evolution and community dynamics to affect the stability and dynamical properties of the system and elucidated the driving factors by manipulating two parameters, namely the cost coefficient ω and the genetic variance G. The cost coefficient limits the extent of trait evolution by selection pressures arising from predation, and genetic variance determines the process rate of evolution relative to the population dynamics of the prey. The parameters related to trait evolution were set to γ = 1, s = 2 and G = 0.01 (a reference value).
The effects of trait evolution on community ecological dynamics may depend on the demographic state of the community without evolution. To characterize the dynamic properties of the community, we used the process rate of the prey population relative to the predator population. This is because previous studies have indicated that the bifurcation in ecological dynamics of the tri-trophic food chain is largely determined by the relative process rate (Hastings and Powell 1991; McCann and Yodzis 1994). McCann and Yodzis (1994) further simplified this system (Eqs. 4–6) in terms of the metabolic rate per unit biomass, which is likely directly associated with or indicated by the mortality rate of the two species, \({d}_{2}\) and \({d}_{3}\), and the assimilation rate relative to the metabolic rate, \({y}_{2}=\frac{{c}_{2}}{{d}_{2}{h}_{2}}\) and \({y}_{3}=\frac{{c}_{3}}{{d}_{3}{h}_{3}}\). Discrepancies in the metabolic rate between the prey and the predator may lead to heteroclinic cycles and chaotic dynamics. Parameters \({y}_{2}\) and \({y}_{3}\) indicate the potential reproductive capacity of the species scaled by the death rate, or the metabolic rate with unlimited consumption and are referred to as the ecological scope. In our simulations, we set the ecological scope to biologically realistic values for the two trophic levels, with two for the prey and five for the predator, followed by other related parameters (Appendix 1). The death rate (metabolic rate) d2 and the maximum consumption rate h2−1 of the prey were manipulated at the same rate so that the ecological scope of the prey was kept constant. This treatment changed the relative process rate between the prey and predator within the constraint of the ecological scope.
Simulation methods
Numerical integrations were performed according to the fourth-order Runge–Kutta method from t = 0 to 105 with \(3\bullet {10}^{5}\) steps using Mathcad 15 (Mathsoft) for each parameter setting. To find the equilibrium points of the dynamical systems, we used the built-in function of Mathematica (Wolfram). “NSolve” was used for communities without evolution and “FindRoot” was used for communities with evolution.
Wavelet transform
The nonlinear system comprises four variables with community dynamics at three trophic levels and the evolutionary dynamics of the prey trait, which exhibited complex and chaotic behaviors depending on the parameter values. To objectively and quantitatively present the fluctuation patterns, the simulated sample paths were wavelet-transformed to extract periodicities at multiple time scales. The wavelet transform to time-series data Xt is \(\text{W}\left(\sigma ,\tau \right)={\left\{\sum _{t}{\left|\psi \left(t-\tau ,\sigma \right)\right|}^{2}\right\}}^{-\frac{1}{2}}\sum _{t}\psi \left(t-\tau ,\sigma \right){X}_{t}\), where \(\psi \left(t,\sigma \right)\) is the mother function, t is the shift parameter (the point of transform in the time series), and \(\sigma\) is the scale parameter, which was set to powers of 2 (2h; h = 1, 2, 3,…) in this study. For \(\psi \left(t,\sigma \right)\), we used the following Mexican hat function: \(\psi \left(t,\sigma \right)\propto \left(1-\frac{{\left(2\pi /\sigma \right)}^{2}{t}^{2}}{2}\right)exp\left[-\frac{{\left(2\pi /\sigma \right)}^{2}{t}^{2}}{4}\right]\).