Complex dynamics in an eco-epidemiological model with the cost of anti-predator behaviors

In this paper, we investigate the complex dynamics of a predator–prey model with disease in the prey, which is characterized by the reduction of prey growth rate due to the anti-predator behavior. The value of this study lies in two aspects: Mathematically, it provides the existence and the stability of the equilibria and gives the existence of Hopf bifurcation. And epidemiologically, we find that the influence of the fear factor is complex: (i) The level of the population density decreases with the increasing of the fear factor; (ii) the cost of fear can destabilize the stability and benefit the emergency of the periodic behavior; and (iii) the high level of fear can induce the extinction of the predator. These results may enrich the dynamics of the eco-epidemiological systems.

ogy due to its universal existence and importance [17]. It is a central topic to study the mechanisms driving predator-prey systems in ecology and evolutionary biology [35]. The research results of prey-predator models are helpful to understand the interactions of different species within a fluctuating natural environment [4,5,10].
Predator-prey interactions are modified by considering various ecological factors [15,26]. In the previous models, predators only affect the density of prey population through direct killing [18]. A new view is that the presence of predators can also change the behavior and physiological characteristics of the prey to a certain extent, so as to affect the amount of the prey population, and its influence degree is even greater than that of direct killing [7,20,28,37,40]. And Zanette et al. [40] experimentally showed that there is about 40% reduction in offspring production of the song sparrows due to the fear from predators. Wang et al. [35] first introduced the fear factor F(k, w), which accounts for the cost of anti-predator defense due to the fear on a predator-prey system. Their insightful work has drawn many scholars' attention and developed some predator-prey models incorporating the fear effect [24,25,30,33,34,36,41].
On the other hand, in the natural world, species do not exist alone, and more and more studies show that the competition of resources and space among populations, as well as infectious diseases, affects population growth. Therefore, the combination of the predator-prey model and the epidemic model, i.e., eco-epidemiological model, has a very high research ability [9,11,12]. Simply speaking, the so-called ecoepidemiological models are demographic models accounting for interactions among different populations in which a disease also spreads [21]. Recently, ecoepidemiological models have drawn much attention of scientists [2,3,6,13,16,29,32,38,39]. At the same time, complex dynamics, such as bistability, quasiperiodicity and chaos, has been found in many ecoepidemiological models [2]. But little attention has been paid so far to the study of the dynamics of ecoepidemiological model incorporating fear factor.
And there comes a question: How does the fear factor affect the persistence and extinction of the predatorprey system?
The main goal of this article is to investigate how the fear factor affects the predation dynamics through studying the stability of the equilibria and Hopf bifurcation of the model.
The rest of this paper is organized as follows: In Sect. 2, we give the details of the model derivation. In Sect. 3, we provide the existence of the equilibria of the model. In Sect. 4, we give the local stability of the equilibria and the existence of Hopf bifurcation. In Sect. 5, we numerically discuss the influence of the cost of fear on the dynamics of the eco-epidemic model. Finally, in Sect. 6, we conclude the ecological significance of the analytical results.

Model derivation
Following [27], we construct an eco-epidemic model of one predator and one prey with a disease in prey population incorporating the fear factor from the predator. First of all, we give some assumptions as follows: (i) Only the prey population is infected by an infectious disease. (ii) In the presence of disease, total prey population is divided into two subclasses, namely susceptible prey population (u) and infected prey population (v). (iii) In the absence of predation, the infected prey population dies before having capability of reproducing. Consider the simple birth-death process of the prey u(t) modeled by where a is the birth rate and d the natural death rate. For the positive growth of the prey, we assume r := a − d > 0. And by introducing the cost of fear according to [35], we obtain: In addition, the susceptible prey population feels the population pressure of both susceptible and infected prey population, i.e., intraspecific and interspecific competition. (iv) The disease in prey population is transmitted vertically from susceptible to infected prey population according to law of mass action at constant rate of infection β. Infected prey population do not recover from the disease. (v) The predator consumes only the infected prey and does not consume the susceptible prey. This is in accordance with the fact that the infected prey are less active and can be caught easily [38], or the experimental observations that some fish (trematode Euhaplorchis californiensis) infect the brain of killifish (Fundulus parvipinnis), and since upper surface contains more oxygen, infected fish come closer to water surface and are eaten by piscivorous birds strikingly about 31 times on average that of non-parasitized fish [1,14,19]. In addition, Peterson and Page [23] observed that wolf often attacks successfully on the moose infected by Echinococcus granulosus.
Based on the assumptions above, the model structure is shown in the transfer diagram of Fig. 1.
And the model can be presented by the following set of ordinary differential equations: with initial conditions where u(t), v(t) and w(t) are the density of the uninfected (susceptible) prey, the infected prey and the b: the density dependent on death rate due to intraspecies competition; β: the contact rate between the susceptible and the infected prey; p: the attack rate on the infected prey; δ: the death rate of the infected prey including natural death rate and disease-induced mortality; c: the conversion coefficient; μ: the natural death rate of the predator; F(k, w): the fear factor. Here, the parameter k represents the level of fear which drives anti-predator behavior of the prey.
Particularly, for convenience of the analysis, similar to that in [35], we adopt the following form for the fear effect term F(k, w): Then, model (1) can be rewritten as follows:

Existence of the equilibria
Following [8,31], we obtain the basic reproduction number R 0 for model (3) as follows: Easy to verify that model (3) has the following trivial and semi-trivial equilibria: (1) Trivial equilibrium: E 0 = (0, 0, 0); (2) Axial equilibrium: bR 0 (δR 0 +r ) and R 0 > 1. Next, we mainly study the existence of the positive equilibrium E * = (u * , v * , w * ) of model (3). Any positive equilibrium (u * , v * , w * ) of model (3) satisfies: and w * = w * (k) is the positive root of the following equation: where has a unique positive root: The discussion of the sign of ξ then we obtain the existence of the endemic equilibrium of model (3).

Theorem 1 Assume r cp > bμ. If
and which indicates that w * is a strictly decreasing function with respect to k. That is, the value of the predator equilibrium w * decreases as k increases. And from (10), we know that there exists a k * such that i.e., if the level of the fear k is larger than k * , the value of the predator w * will tend to 0, and w(t) will almost go to extinction. And we can call k * the high risk value of the fear factor.

Remark 2
In the case without the fear factor, i.e., k = 0, model (3) can be rewritten as It is clear that model (12) has trivial equilibriumÊ 0 = (0, 0, 0), axial equilibriumÊ 1 = r b , 0, 0 and planar And if the conditions in Theorem 1 hold, there is a unique interior equilibrium Comparing the results of the values of equilibria of model (3) with those of model (12), we can see that the fear factor has no influence on the values of the trivial, axial or planar equilibria of model (3). What's more, the cost of fear k does not affect the existence of the interior equilibrium E * of model (3), and it can decrease the values of u * and w * as k increases, but it has no influence on the value of v * directly.

Local stability
Firstly, we directly give the results of stability of trivial and axial equilibria E 0 and E 1 as follows. The proof is standard, and hence, we omit it here. Theorem 2 For model (3), Secondly, we present the stability of the planar equilibrium E 2 of model (3).
The proof of Theorem 3 is complex and given in Appendix B.
Thirdly, we show the stability of the interior equilibrium E * of model (3) by setting (13).
Proof The Jacobian matrix of model (3) around E * = u * , μ cp , w * is given as Hence, the characteristic equation of J is given as where Set Thus, by the Routh-Hurwitz criterion, we can arrive at the conclusion.

Hopf bifurcation
In this subsection, we take k as the bifurcation parameter. The characteristic equation of model (3) at E * = (u * , v * , w * ) is (15), and A i (k)(i = 1, 2, 3) are defined as in (16). There exists a positive constant K H such that (13).

Theorem 5 Model (3) possesses a Hopf bifurcation around E
Proof (II) the other eigenvalue is negative at k = k H .
Next, we verify the transversality condition Putting λ(k) = p 1 (k) + i p 2 (k) into (15), we have Taking derivative of both sides of (20) with respect to k, we obtain Separating the real and imaginary parts, we have where It follows from (21) that We know that when the Hopf bifurcation occurs, the characteristic Eq. (15) has a pair of purely imaginary roots of the form ±i √ A 2 , that is, In order to show that the transversality condition (19) holds, we consider whether D 2 D 4 + D 1 D 3 is zero at k = k H in the following two cases. Therefore, Therefore, From (9), one can obtain that Ψ = − pb (1 + kw * ) w * bδ 2 kμ R 0 2 + w * cdδ kprR 0 + w * bδ kμ r R 0 + r 2 cp 2 w * + cδ pr 2 bδ 2 kμ R 0 2 + cdδ kprR 0 + bδ kμ r R 0 + 2 r 2 ckp 2 w * + cδ kpr 2 + r 2 cp 2 k < 0, which implies that Based on the results above, we obtain: Therefore, conditions (I) and (II) are satisfied and this completes the proof.
For the sake of understanding the complex theoretical results above further, in Fig. 2, we give a bifurcation diagram of model (3) in R 0 − k parameter space. One can see that the fear factor can destabilize the coexistence equilibrium E * of model (3) and will benefit the occurrence of limit cycle oscillation. More precisely, if k < k H , below the Hopf bifurcation, E * is always stable (see Theorem 4); if k = k H , model (3) undergoes a Hopf bifurcation (see Theorem 5) and there is a limit cycle around E * , while if k > k H , E * is unstable (see Theorem 4).

Numerical simulations
In this section, we verify the correctness of the theoretical results by numerical simulation. What's more, we mainly focus on the impact of the level of fear factor k in model (3). As an example, we take the parameters of model (3) as follows: Then, we get: Thus, the conditions of Theorem 1 are satisfied. In addition, model (3) and model (12) have the same boundary equilibria: trivial equilibriumÊ 0 = (0, 0, 0), axial equilibriumÊ 1 = (4, 0, 0) and planar equilibrium E 2 = (0.4, 0.6, 0). Next, we take some different values of k to understand the impact of the fear factor on the resulting dynamics of the system (3).
Example 2 We adopt k = 0.4 < k H ; with the cost of the fear factor, model (3) has a unique positive equilibrium E * = (1.120, 0.375, 0.360). In this case, we obtain −acδ pR 0 kw * = 0.06072 > 0, which means that E * is locally asymptotically stable (see Fig. 4). In addition, the dynamics is similar to that in Fig. 3, and the difference between these two cases is that the value of E * of model (3) is smaller thanÊ * of model (12), which is induced by the cost of fear factor.

Example 3
We adopt k = 0.7615187994 = k H and get: −acδ pR 0 kw * = 1 × 10 −10 ≈ 0, From Theorem 5, we learn that the cost of the fear factor k can destabilize E * from a stable spiral sink In other words, the system exhibits periodic phenomenon. The numerical results are shown in Fig. 5. Comparing Fig. 4 with Fig. 5, one can see that there are two different implications induced by the fear factor k: The first is the decrease of values of u * and w * of E * , and the second is that the stability of E * converts from stable into unstable.

Example 4
We adopt k = 1 > k H and obtain −acδ pR 0 kw * = −0.02628 < 0, which means that E * is unstable and model (3) exhibits periodic phenomenon. The numerical results can be found in Fig. 6, similar to that in Fig. 5. The difference between them is the decrease of value of E * from (0.922,0.375,0.261) to (0.843,0.375,0.222), which is induced by the influence of the fear factor.

Example 5
We analyze the effect of high level of fear factor on the population dynamics of model (3).
As examples, we adopt four different values of k = 10, 100, 1000 and 50000 to show the influence of the fear factor. The simulation results are shown in Fig. 7; where as k increases, the amplitude of w(t) is becom-ing smaller and smaller (see Fig. 7(a), 7(b) and 7(c)). And when k = 50000, 0 < w(t) < 1 × 10 −4 ≈ 0, according to (11), we can claim that the value of the predator w(t) will tend to 0, and w(t) will almost go to extinction (see Fig. 7(d)). In addition, at the same time, u(t) and v(t) oscillate periodically.

Conclusions and discussion
In this paper, we divide the prey population into two subclasses: the susceptible and the infected. To study the influence of the fear factor on the population dynamics, we establish an eco-epidemiological model incorporating the cost of fear in the susceptible prey induced by the predator. Mathematically, we give the conditions for existence (see Theorem 1) and local asymptotic stability (see Theorems 2,3 and 4) of the equilibria of model (3). In addition, we show that the fear factor can destabilize the stability of E * and the model exhibits the Hopf bifurcation (see Theorem 5). Biologically, we summarize our main findings as follows.
(a) If the level of fear factor k < k H , the coexistence equilibrium E * of model (3) is stable (see Fig. 4), which is the same asÊ * of model (12), the case without the fear factor (see Fig. 3). These results agree with those in Refs. [24,25,[33][34][35]41]. (b) If k = k H , the fear effect can destabilize E * and will benefit the occurrence of periodic oscillation (see Fig. 5). In other words, the system has periodic behavior and the model undergoes a limit cycle. These results agree with those in Refs. [24,[33][34][35]41]. (3) is unstable and the system exhibits periodic behavior. These results are similar to those in [24]. (d) If k > k * , the solutions of model (3) will tend to E 2 = (u * , v * , 0), i.e., the predator w will go to extinction, and the fear factor can prevent the occurrence of the limit cycle oscillation and increase the extinction of the system. This can be seen as a population example of "kill 1000 enemies with 800 soldiers killed" extended from "Sun Tzu's Art of War," one of the earliest military works in the world. As we know, k reflects the level of fear factor which drives anti-predator behavior of the prey, and the fear factor is induced by the predator w. Biologically speaking, if the fear factor is very large, the prey will respond to perceived predation risk and show a variety of anti-predator responses (such as changing habitat, forging behavior, vigilance and physiological changes). When the anti-predator response is the key factor in the predator-prey system, there is not enough prey for the predator, which leads to the extinction of the predator. Hence, for the coexistence of system (3), k needs to be less than the high risk value of the fear factor k * ; otherwise, the predator w(t) will almost go extinct, and the prey population will be permanent (see Remark 1). In the meantime, model (3) degenerates to the following limit epidemic model: which has two boundary equilibria E 0 = (0, 0), E 1 = ( r b , 0) and a unique coexistence equilibrium E * = ( δ β , rβ−δ b β (b+β) ). Consider Theorem 1 and Remark 1 again, and we know that the value of u * (k) in E * of model (3) is related to the value of δ β in E * as follows: Furthermore, if rβ > bδ, E * is stable, while if rβ < bδ, E * is unstable. With the parameters set (22), the phase portraits of system (23) are shown in Fig. 8, in which the boundary points E 0 = (0, 0) and E 1 = (4, 0) are saddle points, and the coexistence equilibrium E * = (0.4, 0.6) is globally stable.
To the best of our knowledge, the extinction of the predator w in eco-epidemiological model (3) induced by the fear factor seems to be the first reported case.