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

In this paper, we investigate the effects of pollution on the body size of prey about a predator–prey evolutionary model with a continuous phenotypic trait in a pulsed pollution discharge environment. Firstly, an eco-evolutionary predator–prey model incorporating the rapid evolution is formulated to investigate the effects of rapid evolution on the population density and the body size of prey by applying the quantitative trait evolutionary theory. The results show that rapid evolution can increase the density of prey and avoid population extinction, and with the worsening of pollution, the evolutionary traits becomes smaller gradually. Next, by employing the adaptive dynamic theory, a long-term evolutionary model is formulated to evaluate the effects of long-term evolution on the population dynamics and the effects of pollution on the body size of prey. The invasion fitness function is given, which reflects whether the mutant can invade successfully or not. Considering the trade-off between the intrinsic growth rate and the evolutionary trait, the critical function analysis method is used to investigate the dynamics of such slow evolutionary system. The results of theoretical analysis and numerical simulations conclude that pollution affects the evolutionary traits and evolutionary dynamics. The worsening of the pollution leads to a smaller body size of prey due to natural selection, while the opposite is more likely to generate evolutionary branching.


Introduction
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 irreversibly 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 theoretical 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, considering the influence of rapid evolution on the dynamics of a predator-prey system. The results 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][22][23]26,30]. For example, Fussmann et al. studied the rapid evolutionary response of predator-prey system with sexual predisposition 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 evolutionary 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 frequencydependent 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.

Ecological model in a polluted environment and its dynamic analysis
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 0 (t) is the toxicant concentration in the organism at time t; r 1 (x) is a trait-dependent intrinsic growth 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 1 (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 1 (x) > 0. For simplicity, we take l 1 (x) = l 1 x and l 1 > 0; r 2 is the intrinsic growth rate of the predator; l 2 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]: where c, v and β are positive constants. c is the maximum competition coefficient between prey species with traits x i and x j ; v reflects the intensity of competition. The larger value of v, the stronger the competition; f (x i , x j ) reflects the effect of the species with trait x j on the species with trait x i . In model (1), since the phenotypic trait is the same, then is the predation rate using the following function: and it reaches its maximum when x = x 0 , β 0 is the maximum predation rate, and σ 2 β is the phenotypic variance of the trait x to x 0 ; α > 0 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 periodically at fixed times. Therefore, c 0 (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; kc e (t) represents the species's net uptake of toxicants from the environment; −gc 0 (t) and −mc 0 (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 discharge. Since c 0 (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 ( c 0 (t), c e (t)) T , and for any solutions 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. (7), if β(t) 0 for any t ∈ R and τ 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 From Lemma 2 we know the ultimate boundedness of solutions of system (7), that is, there is a constant M > 0 such that for any positive solution u(t) of system (7), Considering the one-dimensional subsystem of system (6) holds, then system (8) has a unique τ -periodic solution P 20 (t) > 0. For any positive solution P (t) of system (8), we have P (t) → P 20 (t) as t → ∞. Call (0, P 20 (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).
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 20 (t)) is linearly unstable, that is, In the following study, we assume system (6) is uniform persistent, that is, conditions (9) and (10) hold.
Further, we can deduce the following conclusions.
is the solution of system (6). Denote holds, (N * , P * ) satisfies where Since system (6) is uniform persistent, we have Let t → ∞, then equations (13) becomes which is equivalent to Solve the above equations, we have (9) and (10) it is obvious that r 1 (x) − l 1 (x)c 0 > 0 and r 2 > 0 hold, then we have P * > 0.

Remark 1.
In the case of no pollution, the persistent condition of system (6) is a/β(x) > r 2 /r 1 (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 assume ar 1 (x) > r 2 β(x) and al 1 (x) > l 2 β(x) hold true.

Quantitative trait model and its dynamics analysis
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   https://www.journals.vu.lt/nonlinear-analysis discharge amount b 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. 4 Adaptive dynamic model and its dynamics analysis 4.1 Adaptive dynamic model In this subsection, the adaptive dynamics theory [20] is applied to study the long-term evolutionary response of the body size of prey. We assume that evolution is a slow process and ecology is a fast process, that is, evolution and ecology occur at different time scales. Suppose that evolution is a long-term evolutionary behavior, and rare mutant prey population N mut with a slightly different trait y appears in the community of resident prey population N with trait x, and they compete with each other. Then the residentmutant predator-prey model is described by The instantaneous growth rate of rare mutant with trait y is The long-term invasion fitness function is the average function of the instantaneous growth rate of mutant species, then the invasion fitness function [20] is as follows: where N * and P * are given by equation (12). h(y, x) is the long-term average exponential growth rate of the mutant with respect to trait x and h(x, x) = 0. The sign of invasion fitness function determines whether the mutant can successfully invade or not. If the invasion fitness is positive, that is, h(y, x) > 0, then the mutant can invade the resident and will coexist with the resident (see Fig. 5(a)) or replace the resident and become the new resident species. Otherwise, if h(y, x) < 0, the mutant cannot invade the resident (see Fig. 5(b)). It can be seen from equation (15) that the invasion success of a mutant species is related to the intrinsic growth rate of mutant species and resident species, the decrease rate of prey due to the absorption of toxin, the intraspecific and interspecific competition intensity between mutant species and resident species, the discharge amount and discharge period of toxicants and the predation rate of the predator. If the mutations are rare and the mutant trait y is very close to the resident trait x, then the invasion fitness function can be approximately calculated as and r 1 (x) = dr 1 (y)/dy| y=x , l 1 (x) = dl 1 (y)/dy| y=x , f (x, x) = ∂f (y, x)/∂y| y=x . Substitute N * , P * of (12) into the above equation (16), and let then equation (16) can be reduced to the following equation: D(x) is the value of invasion fitness function's partial derivative to mutant trait y at the resident trait x, which is called the invasion fitness gradient and determines the direction of evolutionary change. If D(x) > 0, then only the mutants with trait y > x can invade and replace the resident, whereas if D(x) < 0, then only the mutants with trait y < x can invade residents. If the mutations are rare and random, then the evolution process of residents species with trait x can be well approximated by the following canonical equation of adaptive dynamics proposed by Dieckmann and Law [14]: where , µ is the probability that a newly born prey species is a mutant, σ 2 is the variance of the mutation distribution, µ, σ 2 are constants, 1/2 means that only half of the mutant prey species are at a selective disadvantage and are doomed to be eliminated.

Dynamics analysis
In this subsection, we first analyze the stability of evolutionary traits of adaptive dynamic model (18) to explore the evolutionary dynamics of a biological trait in the presence of mutant prey species. For evolutionary dynamics, it is very important to identify the point x * , where the invasion fitness gradient function D(x * ) = 0. Such trait values are called evolutionary singular strategies or evolutionary singular points. By equation (17), an evolutionary singular strategy x * satisfies the following equation: It is difficult to solve x * analytically from (19) since the trade-off function r 1 (x) does not have the explicit expression. Now we give the stability analysis of evolutionary singular strategy x * by using critical function analysis [20].
Critical function r 1crit (x) is a continuously differentiable function and has the same slope value as the trade-off function r 1 (x) at x * . From equation (19) the slope of the critical function at evolutionary singular strategy x * is so the critical function r 1crit (x) satisfies the following differential equation: In general, the evolutionary singular strategy is the point x * at which the critical function r 1crit (x) is tangent to the trade-off function r 1 (x) (see Fig. 6(a)).
An evolutionary singular strategy x * is called to be locally convergently stable or an attractor (CS) if the resident prey species can be invaded by the mutants whose phenotype trait are closer to x * , which indicates that in the neighborhood of evolutionary singularity x * , the invasion fitness function h(y, x) > 0, that is, if the residents with trait x can be invaded by the mutants with y and x < x * , then x < y < x * and D(x) > 0; if the residents with trait x can be invaded by the mutants with y and x > x * , then x * < y < x and D(x) < 0. So the sufficient condition of the local convergence stability of the evolutionary singularity point x * is as follows: An evolutionary singular strategy x * is called a locally evolutionarily stable strategy (ESS) if any mutant species with trait y = x * fails to evolve and the evolution stops. The fitness function h(y, x * ) as a function of trait y attains a maximum at y = x * . So from equation (15) the sufficient condition of the local evolutionary stability of the evolutionary singularity point x * satisfies where f (x * , x * ) = ∂ 2 f (y, x)/∂y 2 | y=x=x * . An evolutionary singular strategy x * is called a continuously stable strategy (CSS) if it is both convergently stable and evolutionarily stable. The continuously stable strategy indicates that selection pressure drives the evolution of trait towards evolutionary singularity x * and then comes to a halt, and cannot be invaded by the prey species with other traits (A in Fig. 6(b), x * = 0.235).
An evolutionary singular strategy is called an evolutionary branching point (BP) if it is convergently stable but not evolutionarily stable. The presence of evolutionary branching indicates that selection pressure drives the evolution of traits towards the evolutionary singularity, and any mutant near the evolutionary singularity can invade. Over a long period of evolution, the monomorphic resident prey species will split into two prey subspecies with different traits, which promotes diversity of species. (C in Fig. 6(b)), x * = 0.79).
An evolutionary singular strategy x * is called an evolutionary repellor if it is neither convergently stable nor evolutionarily stable ((B in Fig. 6(b), x * = 0.401)).
For the evolutionary singular strategy x * of adaptive dynamic model (18), we have the following theorem.
Theorem 2. If inequalities (20) and (21) are both true, then the evolutionary singular strategy x * is continuously stable; if inequality (20) is true, while inequality (21) is not true, then the evolutionary singular strategy x * is an evolutionary branching point; if neither of inequalities (20) and (21) is true, then the evolutionary singular strategy x * is an evolutionary repellor.

Effects of pollution on the evolutionary trait
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 evolutionary 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 continuously stable or an evolutionary branching point, both of which means x * must be convergently stable and ∂D(x * )/∂x < 0. So the sign of dx * /db is determined by the sign of numerator ∂D(x * )/∂b. By equation (19) we obtain Since al 1 (x) > l 2 β(x) and f (x * , x * ) < 0, we know By (2) we get If x 0 < 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 equation (19) we have Similar to the above analysis, we can get if x 0 < x * , then dx * /dτ > 0. So, with the longer pulse discharge period, the continuously stable trait x * will get larger. This indicates 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 0 , the evolutionary singular strategy remains continuously stable; when the trait x of the mutant is larger than x 0 , 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 0 , the evolutionary singularity strategy remains continuously stable; when the trait x of the mutant is larger than x 0 , the evolutionary singular 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 1 (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.

Conclusions
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 explored. 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.