Study on evolution of a predator–prey model in a polluted environment*

Biological evolution is a common phenomenon in nature. It refers to the process in which an organism interacts with its living environment and its genetic system changes irre- versibly over time, leading to the evolution of corresponding phenotypic characteristics. Some evolution is rapid. Due to the influence of ecological changes and other factors,species selection pressure becomes stronger and rapid evolution occurs. In fishing, larger fish tend to be caught, while slow-growing fish are smaller, which are better able to escape nets and have a higher chance of passing on their genes to the next generation. Some evolution is a long, slow, continuous process in which highly adaptable organisms evolve in their environment through natural selection over generations. In an environment where food was scarce, giraffes with long necks were able to survive by eating leaves from taller trees, while those with short necks were eliminated, so giraffes evolved long necks by natural selection over time.

There are mainly two modelling methods for the evolution of biological phenotypic characteristics: quantitative trait model (QTM) and adaptive dynamic model (ADM). QTM assumes the ecology and evolution occur at the same time scale, and the evolution of quantitative trait is described as a differential dynamical system, which is incorporated into the ecological dynamical system to study the effects of rapid evolution on phenotypic trait and population dynamics. QTM was first proposed by Lande [13], providing a theo- retical basis for the study of the rapid evolution of a biological trait. Next, Saloniemi [25] investigated a coevolutionary quantitative trait model in the predator–prey system, consid- ering the influence of rapid evolution on the dynamics of a predator–prey system. The re- sults showed that rapid evolution promoted the stable coexistence of predator and prey and the existence of periodic solutions. Then the rapid evolution based on the quantitative trait model has been extensively studied [2,5,6,8,11,21–23,26,30]. For example, Fussmann et al. studied the rapid evolutionary response of predator–prey system with sexual predispo- sition trait [5] and genotype trait [6], respectively. Michael et al. in [2] studied a general predator–prey system exhibiting fast evolution in either the predator or the prey. Kaitala et al. in [11] used the quantitative trait evolution model to simulate observational data and concluded that the interaction between predator and prey (bacteria and ciliates) in the microbial system could be best explained as a coevolutionary process, where both the prey and predator evolved. Wang et al. in [30] formulated and explored an eco-evolutionary resource–consumer–predator trophic cascade model incorporating the rapid evolution to study the effects of rapid evolution on the consumer’s body size and to investigate the impact of density-mediate indirect effect on the population dynamics and trait dynamics. ADM assumes that evolutionary population are frequency-dependent and that evo- lutionary and ecological processes occur at different time scales, that is, evolution is a slow process and ecology is a fast process. It provides an important theoretical tool for studying the evolution of species with continuous phenotype trait due to frequency- dependent selection, and its basic framework was proposed by Metz [20], Dieckman [3] and Geritz [7]. ADM is more widely used and is suitable for the study of dynamics such as the evolutionary stability of phenotypic characteristics and evolutionary branching. In recent years, many studies have been conducted on the evolution of biological phenotypic characteristics in predator–prey systems by using adaptive dynamics theory [18,19,24,29, 31,32]. For example, Zu et al. studied the evolutionary response in predator–prey systems with predation rate trait [31] and the prey’s defence ability trait [32], respectively. Wang et al. formulated and investigated a tri-trophic food chain model with foraging effort being the single evolutionary trait and assumed that the predation rate of the top predator to the intermediate predator was the trade-off function of phenotypic characteristics [24]

With the rapid development of industry and agriculture, human production and life lead to the discharge of a large number of toxicants and harmful substances into the surrounding environment, which affect the survival of biological species seriously. With the pollution worsening, there will be obvious changes in morphological characteristics and the density of species. For example, scientists have shown that the body sizes of living polar bears have shrunk compared to their ancestors, the main reason for these changes is environmental pollution, which has affected the evolution of the species. In recent years, much research has been done on the survival and extinction of biological species in a polluted environment by mathematical models [4, 9, 10, 14, 15, 17], especially by impulsive differential equations [9, 10, 14, 15], since in real life, the discharge of some toxins such as waste water, waste gas and waste impurities and the spraying of pesticides are instantaneous and not continuous. Further, some studies have been given to the effects of environmental pollution on the evolutionary changes of species. Liu et al. in [16] considered the effects of continuous discharge pollution on the evolution of phenotypic characteristics of a single species; Veprauskas et al. investigated the dynamics of the daphniid population model with rapid evolution of toxicant resistance trait in a polluted environment [28]. Ackleh et al. developed a discrete-time predator–prey evolutionary model to study how the pest population evolved resistance to the toxicants [1]. In this paper, based on [14], we use QTM and ADM to establish the predator–prey evolution model with the body size as the trait in a pulsed pollution environment to explore the dynamics of fast evolution and slow evolution, respectively, and analyze how the pollution affects the body size of the prey and biodiversity of the species.

The rest of this paper is organized as follows. In the next section, the ecological model of a predator–prey system in a pulsed polluted environment is formulated, and its dynamics are given. In Section 3, we establish a quantitative trait model with the body size of the prey as a continuous phenotypic trait to discuss the effects of rapid evolution and pollution on the population density and the evolutionary trait. In Section 4, an adaptive dynamics model is established to investigate the evolutionary dynamics of slow evolution, and the effects of pollution on evolutionary trait and stability are analyzed theoretically and numerically. We give our conclusions in the last section.

In this section, based on [14], an ecological predator–prey model with the phenotypic trait (the body size of prey) in a polluted environment with pulse toxicants input at fixed times is established as follows:

where N (t) and P (t) are the densities of the prey and predator at time t; c₀(t) is the toxicant concentration in the organism at time t; r₁(x) is a trait-dependent intrinsic grow function of the prey, which exhibits a trade-off between the body size and the intrinsic growth of prey, and is assumed to be a continuous and decreasing function because big individuals generally absorb more toxin and imply a low survival and fecundity; l₁(x) denotes the decrease rate function of prey due to the absorption of the toxin. In general, the amount of toxin absorbed by the organism is related to its body size, and the larger it is, the more toxins it absorbs, thus we assume l₁′ (x) > 0. For simplicity, we take l₁(x) = l₁x and l₁ > 0; r₂ is the intrinsic growth rate of the predator; l₂ denotes the decrease rate of predator due to the absorption of toxin; f (x, x) is the intraspecific competition coefficient, where we adopt the following phenotypic dependent asymmetric competition functions proposed by Kisdi [1]:

and it reaches its maximum when is the maximum predation rate, and is the phenotypic variance of the trait x to denotes the conversion rate; a > 0 denotes the interspecific competition coefficient of predators

In this paper, we assume that the toxicants in the environment are discharged period- ically at fixed times. Therefore, c₀(t) satisfies the following dynamical system in Liu et al. [14]:

where c_e(t) is the concentration of toxicants in the environment at time t; k_ce(t) repre- sents the species’s net uptake of toxicants from the environment; gc₀(t) and mc₀(t) represent the egestion and depuration rates of toxicants in the species, respectively; hc_e(t) is the toxicants loss rate from the environment itself by volatilization and so on; b is the discharge amount of toxicants at time t = n τ ; τ is the impulse period of toxicants dischage

Since c₀(t) and c_e(t) are the concentrations of toxicants, they should satisfy the following inequalities:

From Liu et al. [14] we know it is necessary for (4) to be true that

hold. In the following, we assume that conditions (5) always hold true.

Lemma 1.(See [14].) System (3) has a unique positive τ-periodic solution and for any solutions (c₀(t), c_e(t))^T of system (3), we have

for

Therefore, the dynamics of model (1) is equivalent to the following limit system

Now we give the dynamics of the nonautonomous system

where α(t) and β(t) are τ -periodic continuous functions defined on R, τ > 0

Lemma 2. (See [27].) In system (7), if β(t) > 0 for any t ∈ R and R τ 0 β(t)dt > 0, then system (7) has a unique globally asymptotically stable nonnegative τ -periodic solution u ∗ (t), that is, for any positive solution u(t) of system (7), we have u(t) → u ∗ (t) as and if R τ 0 α(t) dt 6 0, then

From Lemma 2 we know the ultimate boundedness of solutions of system (7), that is, there is a cnstant M > 0 such that for any positive solution u(t) of system (7), limt→∞ sup u(t) ≤ M

Considering the one-dimensional subsystem of system (6)

holds, then system (8) has a unique τ -periodic solution P₂₀(t) > 0. For any positive solution P (t) of system (8), we have P (t)P₂₀(t) as t. Call (0, P₂₀(t)) the semitrivial periodic solution of system (6).

Next, we give the definition and the sufficient and necessary conditions of the uniform persistence for system (6).

Definition 1. System (6) is said to be uniform persistent if there exist constants M > m > 0 such that

And

hold for any positive solution (N (t), P (t)) of system (6)

From Theorem 2 of [27] we obtain

Lemma 3. If condition (9) holds, then system (6) is uniform persistent if and only if the semitrivial periodic solution (0, P₂₀(t)) is linearly unstable, that is

In the following study, we assume system (6) is uniform persistent, that is, condi- tions (9) and (10) hold.

Further, we can deduce the following conclusions.

Theorem 1. Suppose (N (t), P (t)) is the solution of system (6). Denote

If

holds, (N^∗ , P^∗ ) satisfies

where

Proof. From system (6) we have

Since system (6) is uniform persistent, we have

Let . → ∞, then equations (13) becomes

which is equivalent to

Solve the above equations, we have

thus (N^∗, P^∗) satisfies (12). From (11) we know holds true, so N^∗ > 0. From (9) and (10) it is obvious that hold, then we have P^∗ > 0.

Remark 1. In the case of no pollution, the persistent condition of system (6) is a / β(x) > r₂ / r₁(x). Condition (11) in Theorem 1 is actually equivalent to the persistent condition of system (6) under which the toxicant concentration in the organism is averaged, and also can be written as In the following discussion, we assumehold true.

In this section, we study the evolutionary dynamics of the body size of prey by using the method of quantitative trait model. Assume that the evolution process is density dependent, and the evolution occurs rapidly enough such that the ecological processes and evolutionary processes take place on the same time scales. The eco-evolutionary model with the rapid evolution of trait (the body size of prey) is as follows:

where (1/N )(dN/dt) is the fitness of prey, G_x is the genetic variance of the trait x, and G_x > 0.

In the following, by numerical simulations we study the evolutionary dynamics of population density and the body size of prey in the quantitative trait model (14) and further discuss the effects of an ecosystem on evolutionary traits and the feedback mechanism of the traits on the ecosystem with the rapid evolution.

The plots in Figs.1,2(a) and 3(a) show that the rapid evolution of traits promotes the survival of prey, increases the prey density and avoids prey population from extinction, but has little influence on the density of the predator population. It is also found that in the polluted environment, the rapid evolution drives the body sizes of prey to be smaller (see Fig.3).

Next, we study the effects of pollution on the evolutionary trait. The plots in Fig.4 show that pollution affects the phenotypic traits of the prey. With the increase of the

discharge amount . of toxicants, the body size of prey decreases gradually (see Fig.4(a)), while with the prolongation of the impulsive discharge period, the body size of prey increases (see Fig.4(b)). So it can be seen that the worsening pollution will result in a smaller body size of prey.

In this section, suppose that the evolutionary singularity is continuously stable or an evolutionary branching point. We discuss the effects of pollution on the phenotype trait of the prey species theoretically and numerically.

Firstly, the effects of the discharge amount b of toxicants emissions on the evolution- ary singularity x^∗ are discussed. Since there is no explicit analytical solution for x^∗, it is impossible to solve the derivative of x∗ concerning b directly. By the definition of evolutionary singular strategy we know D(x^∗) = 0, which is given by (19). Considering x^∗ as a function of b, by the derivative of an implicit function we have

Since we only consider the case that the evolutionary singular strategy x∗ is con- tinuously stable or an evolutionary branching point, both of which means x∗ must be convergently stable and So the sign of dx^∗/db is determined by the sign of numeratorBy equation (19) we obtain

Since

By (2) we get

If x₀ < x^∗, then β′(x^∗) < 0, and we have ∂D(x^∗)/∂b < 0, thus dx^∗/db < 0. So, with the increase of the discharge amount b of toxicants, the continuously stable trait x^∗ will always decrease. This indicates that the worsening of pollution causes a continuously stable singular strategy to change, and the evolutionary trait tends to be smaller to adapt to the environment. If the evolutionary singular strategy experiences evolutionary branching, the trait value of the branching point will also decrease with an increase in b, suggesting that evolution first makes the prey species smaller and then splits into two different subspecies with different traits.

Similarly, the effect of the pulse discharge period τ on the evolutionary singular strategy x^∗ is discussed. Considering x^∗ as a function of τ , calculate the derivative of x concerning τ , we have

and the sign of dx^∗/dτ is determined by the sign of numerator ∂D(x^∗)/∂τ . By equa- tion (19) we have

Similar to the above analysis, we can get if x₀ < x^∗, then dx^∗/dτ > 0. So, with the longer pulse discharge period, the continuously stable trait x∗ will get larger. This indi- cates that the amount of toxicants in the environment will decrease, which is conducive to the survival of the population, thus making the body size of the prey population larger. If the evolutionary singularity experiences an evolutionary branching point, its trait value will also increase with the prolongation of the impulse discharge period τ , indicating that evolution will first make the prey species larger and then split into two subspecies with different traits.

The influence of pollution on the phenotype trait is further analyzed through numerical simulations. Select b as a bifurcation parameter to simulate the adaptive dynamics of model (18). Figure 7(a) shows that the evolutionary singular strategy x^∗ becomes smaller with the increase of the discharge amount b of toxicants in the environment. When the trait x of the mutant is smaller than x₀, the evolutionary singular strategy remains continuously stable; when the trait x of the mutant is larger than x₀, the evolutionary singular strategy changes from evolutionary branching to continuously stable with the increase of b.

We can also select τ as a bifurcation parameter for analysis. Figure 7(b) shows that evolutionary singular strategy x^∗ strictly increase with the prolongation of the period τ , which means that the increase of τ leads to larger body size of prey and that when the trait x of the mutant are smaller than x₀, the evolutionary singularity strategy remains

continuously stable; when the trait x of the mutant is larger than x₀, the evolutionary sin- gular strategy changes from continuous stable to only convergent stable with the increase of τ .

Therefore, it can be seen from the above analysis that pollution affects the body size of prey and the stability of the evolutionary singularity strategy no matter the exact relationship between the trade-off function r₁(x) and x. The worsening of the pollution (the increase of the discharge amount b of toxicants or shortening the pulse period τ ) promotes evolutionary stability but discourages the presence of evolutionary branching, and drives the body size of prey species to be smaller, while the opposite is more likely to generate evolutionary branching and promote species diversity.

In this paper, the evolutionary response of the phenotypic characteristics of the prey population in the predator–prey model under a pulsed pollution environment was ex- plored. We assumed that the prey species contained only a single phenotypic trait (body size) and that the absorption function of the prey to the toxin, the growth rate of the prey species, the competitive intensity between the prey population, and the predation rate of predators were all related to this trait. Moreover, it was assumed that there was a trade-off relationship between the growth rate of prey species and its body size, the absorption function of the prey species to toxin was proportional to the body size of prey, the interspecies competition function between the prey population adopted the phenotypic dependent asymmetric competition function proposed by Kisdi [12], and the predation rate of the predator to the prey was taken as a power exponential function of the body size of the prey. We discussed the evolutionary dynamics from the rapid evolution and the long-term evolution process by applying the quantitative trait evolutionary theory and adaptive dynamic theory, respectively. Our results showed that pollution affected the evolutionary dynamics and the trait value. The worsening of the pollution led to a smaller body size of prey due to natural selection, while the opposite was more likely to generate evolutionary branching and promote biological diversity.

Recently, there are many researchers on the evolutionary dynamics of single-population model, predator–prey model, competition model, food chain model, infectious disease model, but the influence of pollution on evolution and the trait value is less considered in the polluted environment. Liu et al. [16] investigated the effects of toxicants on the evolving trait of prey species in a continuous polluted environment. They concluded that an increase in the strength of toxicants uptake resulted in a decrease in the trait value of singular evolutionary strategy and promoted the evolutionary stability. In contrast, low levels of toxicants emissions could lead to biological diversity. Our conclusions are consistent with those obtained in [16], which is a good explanation for why the body size of prey species tends to get smaller in a polluted environment.

Finally, we would like to point out that our studies have some limitations. In this paper, we only considered a deterministic evolutionary predator–prey model to study the effects of pollution on the single evolutionary trait and evolutionary dynamics. The growth of species in the natural world is inevitably affected by the environmental fluctuations, so more interesting work is to analyze the influence of the stochastic disturbance on evolutionary changes in adaptive trait dynamics. At the same time, the evolutionary model incorporating more well selected traits, deserves further investigations. We leave these in future work.

This research was supported by the National Natural Science Foundation of China (grant No. 11371030)and the Natural Science Foundation of Liaoning Province (grant No. 20170540001)