Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 1, Number 2, June 2020, pp.181-206 https://doi.org/10.5206/mase/10739 ON MECHANISMS OF TROPHIC CASCADE CAUSED BY ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS YANG WANG AND XINGFU ZOU Abstract. Motivated by a recent field study [Nat. Commun. 7(2016), 10698] on the impact of fear of large carnivores on the populations in a cascading ecosystem of food chain type with the large carnivores as the top predator, in this paper we propose two model systems in the form of ordinary differential equations to mechanistically explore the cascade of such a fear effect. The models are of the Lotka-Volterra type, one is three dimensional and the other four dimensional. The 3-D model only considers the cost of the anti-predation response reflected in the decrease of the production, while the 4-D model considers also the benefit of the response in reducing the predation rate, in addition to the cost by reducing the production. We perform a thorough analysis on the dynamics of the two models. The results reveal that the 3-D model and 4-model demonstrate opposite patterns for trophic cascade in terms of the dependence of population sizes for each species at the co-existence equilibrium on the anti-predation response level parameter, and such a difference is attributed to whether or not there is a benefit for the anti-predation response by the meso-carnivore species. 1. Introduction Predator-prey interactions have attracted the great attention of both ecologists and mathematical biologists, not only because of their vast existence in nature but also because of their diversified forms and rich consequences in the real world. Mathematically, if only considering direct interaction through predation, a classic predator-prey model can be generally described by a system of ordinary differential equations of the form:   du dt = f1(u(t)) −p(u(t),v(t))v(t), dv dt = f2(v(t)) + cp(u(t),v(t))v(t). (1.1) See, e.g., [16, 17, 28, 29]. Here u(t) and v(t) are the populations of the prey and predator respec- tively, f1(u) and f2(v) denote growth functions of the prey and the predator respectively, p(u,v) is the functional response which accounts for the predation rate and biomass transfer from the prey to the predator after predation, and the constant c explains the efficiency in biomass transfer. Predator-prey ODE models of the above form only consider interactions of two species with direct effect reflected by the predation term. However, since 1990s, more and more ecologists have realized the existence of indirect effects (e.g., fear effect), and observed impact of such effects; see, e.g. [1, 14, 21, 26]. Recent field experiments have found the presence of predator itself, can have significant influence on prey’s population through changes in reproduction [15, 35], habitat selection [4, 27] and physiology [3, 5, 34]. In contrast, as far as mathematical modeling is concerned, indirect effects have been largely (if Received by the editors 25 May 2020; accepted 16 June 2020; online first 18 June 2020. 2000 Mathematics Subject Classification. 34K20, 92B05. Key words and phrases. Predator-prey, food chain, anti-predation strategy, stability. Supported by NSERC of Canada (RGPIN-2016-04665). 181 182 Y. WANG AND X. ZOU not all) ignored in those existing models describing predator-prey interactions and those on conservation and management of the ecosystem. Motivated by the field study in [35] which observed an as high as 40% decrease in prey’s reproduction rate when the prey perceived a risk of predator coming from the playback of the predator’s voice, [30] formulated a mathematical model in the form of the following ordinary differential equations   du dt = f(k,v(t)) r0 u(t) −du(t) −au2(t) −p(u(t)) v(t), dv dt = cp(u(t)) v(t) −mv(t). (1.2) Here u is the population of the prey species and v is the population of predator species, and the prey’s growth follows a logistic growth with the intrinsic growth rate being split into the net growth rate r0 and natural death rate d, given by r0 −d. To mimic the scenario of the field experimental study in [35] in which predation actually did not occur due to the use of electronic fence, in (1.2) the fear effect is only incorporated into the production term by the function f(k,v(t)) accounting for a cost. The term au2 = (au)u reflects the self-limiting mechanism of u (due to intra-species competition) and c is the biomass transform efficiency constant. The function p(u) is the functional response which is assumed to depend on the prey population only. Analysis of (1.2), both analytical and numerical, have revealed some interesting dynamics that would have not occurred without considering the fear effect. Since [30], there have been some follow-up modelling works that extend the model (1.2) to accom- modate various aspects of fear effect. For example, [31] considered age structure and discussed different impacts of fear effect on different age stages; [32] explored the fear effect reflected thought the dispersals and its impact on the pattern formation. [8] incorporated an extra food source for the predator in (1.2) and also added a white noise to the death rates of the prey and predator, and analyzed the resulting stochastic model. [18] further considered a digestion delay in addition to the cost of fear, white noise in the death rates, and extra food for the predator but ignored the benefit of the anti-predation response. More recently, in addition to a digestion delay and a cost in prey’s production due to the anti-predation response, [33] further incorporated a benefit term from the anti-predation response, and explored the joint impact of the fear effect and the digestion delay on the population dynamics for both predator and prey. In the real world, it is quite common that a species can predate on one species and in the mean time, it can be a prey of other species. This leads to occurrence of food chains between multiple species, consisting of cascading predator-prey interactions. One may naturally ask how a fear effect arising from one or more species in the chain will affect the dynamics of the whole cascaded populations? In a more recent field experimental study [25], the authors tested a meso-predator cascade by manipulating the large carnivores playback, which resulted in a decrease in the population of meso-carnivore and increase in the population of its prey. That is, the fear effect on the top species in that chain of three species actually affects every species in the ecosystem explicitly or implicitly. To better understand the above mentioned propagation of the fear effect from top layer to the bottom layer in that food chain ecosystem reported in [25], we incorporate a fear effect on the top species into a Lotka-Volterra type food chain model, formulated by the following system of ordinary differential equations ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 183   dN1 dτ = N1 (R1 −a11N1 −a12N2) , dN2 dτ = N2 (R2 −a22N2 −a23N3 + a21N1) , dN3 dτ = N3 (B(α) −D −a33N3 + a32N2) , N1(0) ≥ 0 , N2(0) ≥ 0 , N3(0) ≥ 0. (1.3) Here, N3 is the population of the meso-carnivore (e.g., racoon) which is affected by the large carnivore’s (e.g., wolf, bear) playback; N2 is the prey (e.g., crab) of the meso-carnivore N3, and N1 is the prey of N2. Each species is assumed to follow a logistic growth with growth rates R1, R2 and B(α)−D respectively. The population of the large carnivores does not appear in the system because in the field study [25], only their voices are played, and hence they only have fear (indirect) effect on the meso-carnivores represented by B(α), where the net birth (production) rate B(α) depends on a parameter α standing for the meso-carnivore’s anti-predation response level. By its biological meaning, B(α) is assumed to satisfy B′(α) < 0, B(0) = B3 > 0 and lim α→∞ B(α) = 0. (1.4) The constant D is the natural death rate of N3, aii (i = 1, 2, 3) are the intra-species competition coefficients. a12 and a23 are the predation rates, while a21 and a32 are the conversion rate of the biomass from N1 to N2 and from N2 to N3 respectively; thus, a12/a21 and a23/a32 actually account for the efficiency of biomass transfer from the predations. We point out that this type of three dimensional food chain models have been intensive and extensively studied by some researchers, see, e.g., [9, 11, 12, 13] and the references therein. As far as fear effect in food chains is concerned, two recent papers [19, 20] have also followed line of [30] to consider fear effect in food chain of three species; their scenario is different from ours: they considered other types of functional responses, they incorporated fear effects in the bottom and middle species, and they assumed that the top and middle species are specialist predators. The first goal of this paper is to explore the dynamics of (1.3), particularly, the impact of the meso-carnivore’s anti-predation response level α on the dynamics. In (1.3), only the cost for the meso-carnivore’s anti-predation response is considered. For such a response, in addition to cost, there should also be a benefit (see, e.g., [6, 26]), typically reflected by the decrease in the chance of being predated. A response strategy is expected to seek balance between cost and benefit to achieve certain optimality. In order to also add the benefit into the interplay in the above model system (1.3), we need the term of predation on the meso-carnivore by a large carnivore, which inevitably requires us to add the the population of that large carnivore into the system. This leads to the following four dimensional food chain model  dN1 dτ = N1 (R1 −a11N1 −a12N2) , dN2 dτ = N2 (R2 −a22N2 −a23N3 + a21N1) , dN3 dτ = N3 (B(α,N4) −D3 −a33N3 + a32N2 −a34(α)N4) , dN4 dτ = N4 (−D4 + c̄a34(α)N3) , N1(0) ≥ 0 , N2(0) ≥ 0 , N3(0) ≥ 0 , N4(0) ≥ 0, (1.5) where N4 is the population of the restored top predator (large carnivore) which is assumed to be a specialist predator with the mortality rate D4. Now the net growth function of N3 depends not only 184 Y. WANG AND X. ZOU on the anti-predation response level α but also on the population of its predator N4. The function a34(α) denotes the encounter rate between N3 and N4 which is affected by the protective behaviours of N3 species characterized by its dependence on the anti-predation response level α. By their biological meanings of B(α,N4) and a34(α), they are assumed to satisfy the following conditions: { B(α,N4) is decreasing in α and N4, B(0,N4) = B(α, 0) = B3 > 0, lim α→∞ B(α,N4) = lim N4→∞ B(α,N4) = 0, (1.6) and a34(α) is decreasing, a34(0) = a0 > 0, lim α→∞ a34(α) = 0. (1.7) Finally, c̄ is the efficiency of biomass transform. So in this model, we consider both the complex multi-trophic predator-prey structure and the trade-off from anti-predation response. The remainder of this paper is organized as follows. In Section 2, we analyze the model system (1.3). We establish the well-posedness of system (1.3) and find the condition for existence and stability for all its equilibrium solutions. We also discuss the relationship between the anti-predation level and the final population size. In the end, some numerical examples, together with some discussions, are given to demonstrate our results. In Section 3, we investigate the dynamics of the four dimension model (1.5), including the existence and stability of equilibria as well as the continuously dependence between the final population size with respect to the anti-predation level. We also discuss the difference of the results from those for (1.3) in Section 2. In addition, we also present some numerical examples to illustrate that different functional response functions may lead to slightly different dynamical behaviour of the solution. In Section 4, we summarize our main results and discuss their biological implications. We also discuss some possible future projects along this direction of anti-predation response in predator-prey interactions. 2. Analysis of the model without large carnivores In this section, we analyze the three-species model (1.3). We first show the well-posedness of the system (1.3). Then we find all equilibrium solutions and discuss their stability in terms of the parameter values and anti-predation strategy level. As we mentioned in last section, this kind of food chain models have been studied in literatures, and thus, some technical results can be found in existing researches, e.g., [9, 11, 12, 13] and their references. But we need to associate the results to the new parameter α, the anti-predation response level of the species N3 to shy light on influence of the fear effect for this model. 2.1. Preliminaries. For mathematical simplification, we first non-dimensionalize the model (1.3). Let t = R1τ, x = a11N1 R1 , y = a12N2 R1 , z = a23N3 R1 , then model (1.3) becomes   dx dt = x(1 −x−y), dy dt = y(k −d1y −z + β1x), dz dt = z(f(α) −d2z + β2y), (2.1) ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 185 where k = R2 R1 , d1 = a22 a12 , d2 = a33 a23 , β1 = a21 a11 , β2 = a32 a12 , f(α) = B(α) −D R1 , f(0) = B3 −D R1 , lim α→∞ f(α) = − D R1 . By the basic theory of ODE systems, we can easily show that the initial value problem for (2.1) has a unique solution; moreover, the solution is nonnegative (positive) with nonnegative (positive) initial conditions because each equation in (2.1) is of Gaussian type. Now we show that the solution to system (2.1) is bounded. From the first equation in system (2.1) and by the non-negativity of y(t), we have dx dt = x(1 −x−y) ≤ x(1 −x). By the comparison theorem [24], we can obtain limt→∞ sup x(t) ≤ 1. Therefore, for any �1 > 0, there holds x(t) ≤ 1 + �1 for large t. Incorporating this estimate for large t into the second equation in (2.1) results in dy dt = y(k −d1y −z + β1x) ≤ y(k + β1(1 + �1) −d1y), for large t. Applying the comparison theorem again, we then obtain lim t→∞ sup y(t) ≤ k + β1(1 + �1) d1 . Since �1 > 0 is arbitrary small, the above inequality actually implies lim t→∞ sup y(t) ≤ k + β1 d1 . Incorporating the above inequality into the third equation in (2.1) and by the same argument, we can obtain lim t→∞ sup z(t) ≤ max ( f(α)d1 + β2(k + β1) d1d2 , 0 ) . Combining the above, we have proved that the solution (x(t),y(t),z(t)) to (2.1) is bounded. 2.2. Existence and Stability of the boundary equailibria. In this section, we find all boundary equilibrium solutions and give the condition for their existence and stability. For the result to be biologically meaningful, we are only interested in equilibria with nonnegative components. An equilibrium of (2.1) solves the following system  x(1 −x−y) = 0, y(k −d1y −z + β1x) = 0, z(f(α) −d2z + β2y) = 0. There are seven boundary equilibrium solutions. E0 = (0, 0, 0) is the trivial equilibrium solution which always exists. E1 = (1, 0, 0), E2 = (0,k/d1, 0) and E3 = (0, 0,f(α)/d2) are the equilibria representing the scenario that only one species survives; E1 and E2 always exist, while E3 exists only when f(α) > 0. 186 Y. WANG AND X. ZOU There are also other three possible equilibria corresponding to the scenario of two species coexisting, and they are given by E12 = ( 1 − k + β1 d1 + β1 , k + β1 d1 + β1 , 0 ) , E13 = ( 1, 0, f(α) d2 ) , E23 = ( 0, d2k −f(α) d1d2 + β2 ,k − d1d2k −d1f(α) d1d2 + β2 ) . By the nonnegative requirement, E12 exists when k < d1, E13 exists when f(α) > 0 and E23 exists when d2k > f(α) > −β2k/d1. The local stability of an equilibrium is obtained by linearization at the equilibrium. The Jacobian matrix at equilibrium (x∗,y∗,z∗) is given by J (E = (x∗,y∗,z∗)) =   1 − 2x∗ −y∗ −x∗ 0 β1y ∗ k − 2d1y∗ + β1x∗ −z∗ −y∗ 0 β2z ∗ f(α) − 2d2z∗ + β2y∗   . (2.2) At E0 = (0, 0, 0), the Jacobian is given by J(E0) =   1 0 0 0 k 0 0 0 f(α)   , therefore E0 is unstable, as there are positive eigenvalues λ1 = 1 and λ2 = k. Similarly, at E1 = (1, 0, 0), the Jacobian is given by J(E1) =   −1 −1 0 0 k + β1 0 0 0 f(α)   , therefore E1 is unstable, as there is a positive eigenvalue λ = k + β1. At E2 = (0,k/d1, 0), the Jacobian is given by J(E2) =   1 − k d1 0 0 β1k d1 −k − k d1 0 0 f(α) + β2k d1   , therefore E2 is asymptotically stable if and only if 1 −k/d1 < 0 and f(α) < −β2k/d1. ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 187 At E3 = (0, 0,f(α)/d2), the Jacobian is given by J(E3) =   1 0 0 0 k − f(α) d2 0 0 β2f(α) α2 −f(α)   , hence, E3 is unstable since there is a positive eigenvalue λ = 1. At E12, the Jacobian becomes J(E12) =   k + β1 d1 + β1 − 1 k + β1 d1 + β1 − 1 0 β1k + β 2 1 d1 + β1 − d1k + d1β1 d1 + β1 − k + β1 d1 + β1 0 0 f(α) + β2k + β2β1 d1 + β1   , thus, E12 is asymptotically stable if and only if f(α) < − β2k + β2β1 d1 + β1 . At E13 = (1, 0,f(α)/d2), the Jacobian reduces to J(E13) =   −1 −1 0 0 k + β1 − f(α) d2 0 0 β2f(α) α2 −f(α)   , therefore E13 is asymptotically stable if and only if f(α) > d2(k + β1). At E23 = ( 0, d2k −f(α) d1d2 + β2 ,k − d1d2k −d1f(α) d1d2 + β2 ) , the Jacobian is given by J(E23) =   1 − d2k −f(α) d1d2 + β2 0 0 β1d2k −β1f(α) d1d2 + β2 − d1d2k −d1f(α) d1d2 + β2 − d2k −f(α) d1d2 + β2 0 β2k − β2d1d2k −β2d1f(α) d1d2 + β2 −d2k + d1d 2 2k −d1d2f(α) d1d2 + β2   , therefore E23 is asymptotically stable if and only if f(α) > d2k −d1d2 −β2. We summarize the above analysis in the following proposition. Proposition 2.1. For system (2.1), the following statement hold. (i) Equilibrium E0 and E1 always exist and are unstable. (ii) Equilibrium E2 always exists; it is asymptotically stable if and only if k > d1 and f(α) < −β2k/d1. 188 Y. WANG AND X. ZOU (iii) When f(α) > 0, E3 exists; it is unstable. E13 exists and is asymptotically stable if and only if f(α) > d2(k + β1). (iv) When k < d1, E12 exists; it is asymptotically stable if and only if f(α) < −β2(k +β1)/(d1 +β1). (v) When d2k > f(α) > −β2k/d1, E23 exists; it is asymptotically stable if and only if f(α) < d2k −d1d2 −β2. 2.3. Existence and stability of a positive equilibrium solution. There is a unique positive equi- librium solution E∗ = (x∗,y∗,z∗) if d2(k + β1) > f(α) > max ( d2k −d1d2 −β2, −β2(k + β1) d1 + β1 ) . (2.3) Indeed, E∗ = (x∗,y∗,z∗) is given by   x∗ = 1 −y∗, y∗ = d2(k + β1) −f(α) d2(d1 + β1) + β2 , z∗ = k + β1 − (d1 + β1)y∗. (2.4) The conditions in (2.3) directly comes from the formulas in (2.4). The Jacobian matrix at the positive equilibrium can be simplified as J(E∗) =   −x∗ −x∗ 0 β1y ∗ −d1y∗ −y∗ 0 β2z ∗ −d2z∗   , Thus, the corresponding characteristic equation is λ3 + (d1y ∗ + d2z ∗ + x∗)λ2 + (d1d2y ∗z∗ + β1x ∗y∗ + β2y ∗z∗ + d1x ∗y∗ + d2x ∗z∗)λ + (β1d2 + d1d2 + β2)x ∗y∗z∗ = 0. where   a1 = d1y ∗ + d2z ∗ + x∗ > 0, a2 = d1d2y ∗z∗ + β1x ∗y∗ + β2y ∗z∗ + d1x ∗y∗ + d2x ∗z∗ > 0, a3 = (β1d2 + d1d2 + β2)x ∗y∗z∗ > 0, and a1a2 −a3 = d21d2(y ∗)2z∗ + d1d 2 2y ∗(z∗)2 + β1d1x ∗(y∗)2 + β2d1(y ∗)2z∗ + β2d2y ∗(z∗)2 + d21x ∗(y∗)2 + 2d1d2x ∗y∗z∗ + d22x ∗(z∗)2 + β1(x ∗)2y∗ + d1(x ∗)2y∗ + d2(x ∗)2z∗ > 0. By Routh-Hurwitz criterion, the positive equilibrium is locally asymptotically stable. Moreover, we can prove that it is actually globally asymptotically stable as long as it exists (i.e., (2.3) holds). Theorem 2.1. If (2.3) holds, then the positive equilibrium E∗ is globally asymptotically stable. Proof. Consider the Lyapunov function V (x,y,z) = β1β2(x−x∗ −x∗ ln x x∗ ) + β2(y −y∗ −y∗ ln y y∗ ) + (z −z∗ −z∗ ln z z∗ ), then dV dt = −β1β2(x−x∗)2 −β2(y −y∗)2 − (z −z∗)2 ≤ 0 ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 189 and dV dt = 0 if and only if (x,y,z) = (x∗,y∗,z∗). By LaSalle’s Invariant Principle, we conclude E∗ is globally asymptotically stable. � For readers’ convenience, we summarize the analytical results on the dynamics of the model (2.1) obtained above in the following Table 1. Table 1. Condition of existence and stability of the equilibria in model (2.1) Equilibrium solution Existence Stability E0 Always Unstable E1 Always Unstable E2 Always d1 < k and f(α) < −β2k/d1 E3 f(α) > 0 Unstable E12 d1 > k f(α) < −β2(k + β1)/(d1 + β1) E13 f(α) > 0 f(α) > d2(k + β1) E23 d2k > f(α) > −β2k/d1 f(α) < d2k −d1d2 −β2 E∗ d2(k + β1) > f(α) > max ( d2k −d1d2 −β2,− β2(k + β1) d1 + β1 ) GAS From Theorem 2.1 , we know that the positive equilibrium is always globally asymptotically stable as long as it exists (i.e., (2.3) holds), implying that the populations of all three species will converge to co-existence state at the respect levels x∗, y∗ and z∗. Thus, it is worthwhile to investigate how the response strength α will affect these levels. Indeed, direct calculations give  dy∗ dα = −f′(α) d2(d1 + β1) + β2 > 0, dx∗ dα = − dy∗ dα < 0, dz∗ dα = −(d1 + β1) dy∗ dα < 0. That is, within the range of α that guarantees (2.3), with the increase of the anti-predation response strength α, the final population sizes of the meso-carnivore, its prey and the prey’s prey will decrease, increase and decrease respectively, demonstrating an alternative pattern for the fear effect in the cascade, which was observed in the field study [25]. Therefore, the model (1.3) does provide a mechanism that can explain the phenomenon of trophic cascade caused by a fear of large carnivores reported in [25]. 2.4. Numerical simulations. In this subsection, we present some numerical simulations to illustrate the analytical results obtained above. For this purpose, we choose a particular form for the function B(α) given by B(α) = R3 1 + cα , 190 Y. WANG AND X. ZOU and this sends (1.3) to the following system:  dN1 dτ = N1 (R1 −a11N1 −a12N2) , dN2 dτ = N2 (R2 −a22N2 −a23N3 + a21N1) , dN3 dτ = N3 ( R3 1 + cα −D −a33N3 + a32N2 ) , N1(0) ≥ 0 , N2(0) ≥ 0 , N3(0) ≥ 0, (2.5) We fix the parameters R1 = 1, a11 = 1, a12 = 0.4, R2 = 1, a22 = 0.2, a23 = 0.5, a21 = 0.5, R3 = 3, D = 1, a33 = 0.5, a32 = 0.05, c = 0.4, (2.6) and demonstrate how α impacts the population dynamics. In this case, k = R2/R1 = 1 and d1 = a22/a12 = 1/2 so k > d1. According to Proposition 2.1, there are three threshold values for α, denoted by α∗1, α ∗ 2 and α ∗ 3, which are given by  f (α∗1) = d2(k + β1), f (α∗2) = d2k −d1d2 −β2, f (α∗3) = − β2k d1 . Using the parameter values in (2.6), we obtain α∗1 = 0.5, α ∗ 2 ≈ 2.954545455 and α∗3 = 7.5. By Proposition 2.1, when 0 < α < 0.5, the equilibrium E13 is stable (as demonstrated in Figure 1-(a) for α = 0.4); when 0.5 < α < 2.954545455, the coexistence equilibrium E∗ is stable (as demonstrated in Figure 1-(b) for α = 2); when 2.954545455 < α < 7.5 (destroying (2.3), hence E∗ no longer exists), the equilibrium E23 is stable (as demonstrated in Figure 1-(c) for α = 5); when α > 7.5, the equilibrium E2 is stable (as demonstrated in Figure 1-(d) for α = 10). The bifurcation diagram with respect to α is given in Figure 2. Now, we change R1 to R1 = 3 and a33 to a33 = 0.1 and keep other parameters the same as in (2.6). Then k = R2/R1 = 1/3 and d1 = a22/a12 = 1/2 leading to the scenario of k < d1. According to Proposition 2.1, there are two threshold values for α, denoted by α1 and α2, which are given by  f (α1) = d2(k + β1), f (α2) = − β2(k + β1) d1 + β1 . Using the parameter values in (2.6), we obtain α1 = 2.5 and α2 ≈ 8.409090911. By Proposition 2.1, when 0 < α < 2.5, the equilibrium E13 is stable (as demonstrated in Figure 3-(a) for α = 2); when 2.5 < α < 8.409090911, the coexistence equilibrium E∗ is stable (as demonstrated in Figure 3-(b) for α = 6); when α > 8.409090911, the equilibrium E12 is stable (as demonstrated in Figure 3-(c) for α = 10). The bifurcation diagram with respect to α is given in Figure 4. We find that depending on the difference between k and d1, we have two kinds of bifurcation. In both cases, when the anti-predation response α passes a threshold, it leads to a transcritical bifurcation and we can always observe a meso-predator cascade inside the coexistence region when increasing α: the population of meso-predator is decreasing, the population of its prey is increasing and the population of the prey’s prey (bottom prey) is decreasing. However, this pattern will be dramatically changed when we restore large carnivores instead of only manipulating their playback to induce fear. In the next section, we will model the case when we also introduce the large carnivores back into the food chain, leading to a 4-D model. ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 191 0 20 40 60 80 100 time 0 1 2 3 4 p o p lu a ti o n Solution of the three dimension system = 0.4 level-1 level-2 level-3 (a) 0 20 40 60 80 100 time 0 0.5 1 1.5 2 2.5 p o p lu a ti o n Solution of the three dimension system = 2 level-1 level-2 level-3 (b) 0 20 40 60 80 100 time -1 0 1 2 3 4 5 p o p lu a ti o n Solution of the three dimension system = 5 level-1 level-2 level-3 (c) 0 20 40 60 80 100 time 0 1 2 3 4 5 6 p o p lu a ti o n Solution of the three dimension system = 10 level-1 level-2 level-3 (d) Figure 1. Population dynamics of (2.5) when d1 < k. (a) When 0 < α = 0.4 < 0.5, the equilibrium E13 is stable, (b) 0.5 < α = 2 < 2.954545455, the coexistence equilibrium E∗ is stable, (c) when 2.954545455 < α = 5 < 7.5, the equilibrium E23 is stable, (d) when α = 10 > 7.5, the equilibrium E2 is stable. 3. Model with restoring large carnivores In this section, we analyze (1.5) which has the population of large carnivores incorporated together with a benefit in preventing predation of the meso-carnivore by the large carnivores, in addition to the cost in the meso-carnivore’s production. Parallel to Section 2, we first establish the well-posedness of the 4-D model (1.5), discuss all possible equilibrium solutions and find the condition for their existence and stability in terms of the parameter values and anti-predation strategy level. 3.1. Well-posedness. For mathematical simplification, we still apply non-dimensionalization for our model (1.5) which is a natural expansion of the non-dimensionalization for (1.3) in Section 2, given by t = R1τ, x = a11N1 R1 , y = a12N2 R1 , z = a23N3 R1 , w = a0N4 R1 , 192 Y. WANG AND X. ZOU 0 2 4 6 8 10 0 1 2 3 4 5 p o p lu a ti o n Bifurcation diagram when k>d 1 level-1 level-2 level-3 (a) 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 p o p lu a ti o n Meso-predator cascade in the coexistence region level-1 level-2 level-3 (b) Figure 2. (a) Bifurcation diagram of (2.5) when k > d1, (b) In the region when the coexistence equilibrium E∗ is stable, we can observe a meso-predator cascade. with a0 = a34(0). Then model (1.5) becomes  dx dt = x(1 −x−y), dy dt = y(k −d1y −z + β1x), dz dt = z (f(α,w) −d2z + β2y −p(α)w) , dw dt = w(−m + cp(α)z) (3.1) where k = R2 R1 , d1 = a22 a12 , d2 = a33 a23 , β1 = a21 a11 , β2 = a32 a12 , m = D4 R1 , c = c̄a0 a23 . For the transformed and rescaled functions f(α,w) and p(α), the conditions (1.6) and (1.7) are trans- formed to   f(α,w) = B(α,N4) −D3 R1 , f(α, 0) = f(0,w) = f0 = B3 −D3 R1 , ∂f ∂α < 0, ∂f ∂w < 0, lim α→∞ f(α,w) = lim w→∞ f(α,w) = −D3 R1 , p(α) = a34(α) a0 , dp dα < 0, p(0) = 1 and lim α→∞ p(α) = 0. (3.2) ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 193 0 20 40 60 80 100 time 0 2 4 6 8 p o p lu a ti o n Solution of the three dimension system = 2 level-1 level-2 level-3 (a) 0 20 40 60 80 100 time 1 2 3 4 5 6 p o p lu a ti o n Solution of the three dimension system = 6 level-1 level-2 level-3 (b) 0 20 40 60 80 100 time 0 1 2 3 4 5 6 p o p lu a ti o n Solution of the three dimension system = 10 level-1 level-2 level-3 (c) Figure 3. Population dynamics of (2.5) when d1 > k. (a) When 0 < α = 2 < 2.5, the equilibrium E13 is stable, (b) 2.5 < α = 6 < 8.409090911, the coexistence equilibrium E∗ is stable, (c) when 8.409090911 < α = 10, the equilibrium E12 is stable. By the fundamental theory of ODEs, we easily see that the initial value problem associated with (3.1) has a unique solution; moreover, the solution is nonnegative (positive) with nonnegative (positive) initial conditions. Next we show that the solution to system (3.1) is bounded. Firstly, by the same argument as in subsection 2.2, we can obtain the same estimates for x(t) and y(t): lim t→∞ sup x(t) ≤ 1 and lim t→∞ sup y(t) ≤ k + β1 d1 . For z and w, we consider P = cz + w. Then we have dP dt = cz (f(α,w) −d2z + β2y) −mw. 194 Y. WANG AND X. ZOU 0 2 4 6 8 10 0 5 10 15 20 p o p lu a ti o n Bifurcation diagram when k k, (b) In the region when the coexistence equilibirum E∗ is stable, we can observe a meso-predator cascade. Now for any given � > 0, there holds y(t) ≤ (k + β1)(1 + �)/d1 for large t. Then for some µ ∈ (0,m) and large t, we have dP dt + µP = cz (µ + f(α,w) −d2z + β2y) − (m−µ)w, ≤ cz ( µ + f0 + β2(k + β1)(1 + �) d1 −d2z ) =: cz(K −d2z) (3.3) where K = µ + f0 + β2(k + β1)(1 + �) d1 . If K ≤ 0, since z(t) ≥ 0, we can conclude that dP dt + µP ≤ 0, leading to P ≤ P0e−µt, where P0 = cz0 + w0 comes from the initial condition. Thus, lim t→∞ sup P(t) ≤ 0. If K > 0, then dP dt + µP ≤ cK 2 4d2 . By the comparison theorem [24], we can obtain P ≤ P0e−µt + ( 1 −e−µt ) cK2 4µd2 . which implies lim t→∞ sup P(t) ≤ cK 2 4µd2 . Thus, in both cases, P is bounded, implying that z and w are both bounded. Therefore, all four components of the solution are bounded. ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 195 3.2. Existence and Stability of the boundary equailibrium solutions. As in Section 2, we first analyze the boundary equilibria of the model (3.1). Again, we are only interested in the equilibria with nonnegative components. By solving the system   x(1 −x−y) = 0, y(k −d1y −z + β1x) = 0, z (f(α,w) −d2z + β2y −p(α)w) = 0, w(−m + cp(α)z) = 0, we find that there are eleven possible boundary equilibria and they are described below. E0 = (0, 0, 0, 0) is the trivial equilibrium solution which always exists, E1 = (1, 0, 0, 0), E2 = (0,k/d1, 0, 0) and E3 = (0, 0,f0/d2, 0) are the equilibria that correspond to the case of only one species surviving: E1 and E2 are always exist while E3 exists only when f0 > 0. There are also four possible equilibria accounting for the scenario of two species coexisting and they are E12 = ( d1 −k d1 + β1 , k + β1 d1 + β1 , 0, 0 ) , E13 = ( 1, 0, f0 d2 , 0 ) , E23 = ( 0, d2k −f0 d1d2 + β2 , kβ2 + d1f0 d1d2 + β2 , 0 ) , E34 = ( 0, 0, m cp(α) ,w ) where w satisfies the equation f(α,w) − d2m cp(α) −p(α)w = 0. (3.4) Due to the requirement of non-negativity, E12 exists when k < d1, E13 exists when f0 > 0, and E23 exists when d2k > f0 > −β2k/d1. For the existence and uniqueness of E34, we denote F1(w) = f(α,w) − d2m cp(α) −p(α)w. Then it is obvious that F1(w) is a decreasing function with respect to w and limw→∞F1(w) = −∞. Thus, the sufficient and necessary condition for E34 to exist is F1(0) > 0, that is f0 > d2m/cp(α). There are three possible equilibria for the case of three species coexisting, given by E123 = ( d1d2 −d2k + β2 + f0 β1d2 + d1d2 + β2 , β1d2 + d2k −f0 β1d2 + d1d2 + β2 , β1β2 + β1f0 + β2k + d1f0 β1d2 + d1d2 + β2 , 0 ) , E134 = ( 1, 0, m cp(α) ,w ) , E234 = ( 0, cp(α)k −m d1cp(α) , m cp(α) , ŵ ) , where w̄ is as in (3.4) and ŵ satisfies the equation f(α,ŵ) − d2m cp(α) + β2cp(α)k −β2m d1cp(α) −p(α)ŵ = 0. (3.5) Thus, E123 exists when d2(k + β1) > f0 > max ( d2k −d1d2 −β2,− β2(k + β1) d1 + β1 ) , 196 Y. WANG AND X. ZOU The condition for E134 to exist is the same as for the existence of E34, that is, when f0 > d2m/cp(α). Condition for E234 to exist is cp(α)k −m > 0 and f0 > d2m cp(α) − β2cp(α)k −β2m d1cp(α) . In order to discuss the local stability, we calculate the Jacobian matrix at equilibrium E = (x∗,y∗,z∗,w∗) as J(E) =   1 − 2x∗ −y∗ −x∗ 0 0 β1y ∗ k − 2d1y∗ + β1x∗ −z∗ −y∗ 0 0 β2z ∗ J33 J34 0 0 cp(α)w∗ −m + cp(α)z∗   , (3.6) where J33 = f (α,w ∗) − 2d2z∗ + β2y∗ −p(α)w∗, and J34 = z ∗ ∂f(α,w) ∂w ∣∣∣∣ w=w∗ −p(α)z∗ < 0 for z∗ > 0, w∗ > 0. At E0 = (0, 0, 0, 0), the Jacobian becomes J(E0) =   1 0 0 0 0 k 0 0 0 0 f0 0 0 0 0 −m   , thus, E0 is unstable. At E1 = (1, 0, 0, 0), the Jacobian reduces to J(E1) =   −1 −1 0 0 0 k + β1 0 0 0 0 f0 0 0 0 0 −m   , hence E1 is unstable. ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 197 At E2 = (0,k/d1, 0, 0), the Jacobian is given by J(E2) =   1 − k d1 0 0 0 β1k d1 −k − k d1 0 0 0 f0 + β2k d1 0 0 0 0 −m   , so E2 is asymptotically stable if and only if 1 −k/d1 < 0 and f0 < −β2k/d1. At E3 = (0, 0,f0/d2, 0), the Jacobian is now J(E3) =   1 0 0 0 0 k − f0 d2 0 0 0 β2f0 α2 −f0 J34 0 0 0 −m + cp(α)f0 d2   , therefore E3 is unstable. At E12, the Jacobian is reduced to J(E12) =   k −d1 d1 + β1 k −d1 d1 + β1 0 0 β1k + β 2 1 d1 + β1 − d1k + d1β1 d1 + β1 − k + β1 d1 + β1 0 0 0 f0 + β2k + β2β1 d1 + β1 0 0 0 0 −m   , therefore E12 is asymptotically stable if and only if f0 < −(β2k + β2β1) d1 + β1 . 198 Y. WANG AND X. ZOU At E13, the Jacobian is J(E13) =   −1 −1 0 0 0 k + β1 − f0 d2 0 0 0 β2f0 α2 −f0 J34 0 0 0 −m + cp(α)f0 d2   , consequently, E13 is asymptotically stable if and only if md2 cp(α) > f0 > d2(k + β1). At E23, the Jacobian is given by J(E23) =   1 − d2k −f0 d1d2 + β2 0 0 0 β1d2k −β1f0 d1d2 + β2 − d1d2k −d1f0 d1d2 + β2 − d2k −f0 d1d2 + β2 0 0 kβ22 + d1f0β2 d1d2 + β2 − kβ2d2 + d1f0d2 d1d2 + β2 J34 0 0 0 −m + cp(α) (kβ2 + d1f0) d1d2 + β2   , therefore E23 is asymptotically stable if and only if f0 > d2k −d1d2 −β2 and p(α) < m (d1d2 + β2) c (kβ2 + d1f0) . At E34 = ( 0, 0, m cp(α) ,w ) , the Jacobian becomes J(E34) =   1 0 0 0 0 k − m cp(α) 0 0 0 β2m cp(α) −d2m cp(α) J34 0 0 cp(α)w 0   , therefore E34 is unstable. ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 199 At E123 = (X1,Y1,Z1, 0), the Jacobian is given by J(E123) =   −X1 −X1 0 0 β1Y1 −d1Y1 −Y1 0 0 β2Z1 −d2Z1 0 0 0 0 −m + cp(α)Z1   , where X1, Y1 and Z1 denote the three positive components in E123. In Section 2, for the three dimensional model (2.1), we have proved the principle 3 × 3 sub-matrix of J(E123) only has negative eigenvalues. Therefore E123 is asymptotically stable if and only if −m + cp(α)Z1 < 0, that is p(α) < m (β1d2 + d1d2 + β2) c (β1β2 + β1f0 + β2k + d1f0) . At E134 = ( 1, 0, m cp(α) ,w ) , the Jacobian is given by J(E134) =   −1 −1 0 0 0 k + β1 − m cp(α) 0 0 0 β2m cp(α) − d2m cp(α) J34 0 0 cp(α)w 0   , Since J34 < 0, E123 is asymptotically stable if and only if k + β1 − m/cp(α) < 0, that is p(α) < m/c (k + β1). At E234, the Jacobian becomes J(E234) =   1 − cp(α)k −m d1cp(α) 0 0 0 β1 (cp(α)k −m) d1cp(α) − d1 (cp(α)k −m) d1cp(α) − cp(α)k −m d1cp(α) 0 0 β2m cp(α) − d2m cp(α) J34 0 0 cp(α)ŵ 0   , the lower 3 × 3 principle sub-matrix can be written as A =   −d1Y2 Y2 0 β2Z2 −d2Z2 J34 0 cp(α)ŵ 0   , 200 Y. WANG AND X. ZOU where Y2 = cp(α)k −m d1cp(α) and Z2 = m cp(α) . Then the characteristic polynomial of matrix A is λ3 + (d1Y2 + d2Z2)λ 2 + (−cp(α)ŵJ34 + d2Z2d1Y2 + β2Z2Y2) λ− cp(α)J34d1ŵY2 = 0, where   a1 = d1Y2 + d2Z2 > 0, a2 = −cp(α)ŵJ34 + d2Z2d1Y2 + β2Z2Y2 > 0, a3 = cp(α) −J34d1ŵY2 > 0, and a1a2 −a3 = −cp(α)J34d2ŵZ2 + d21d2Y 2 2 Z2 + d1d 2 2Y2Z 2 2 + β2d1Y 2 2 Z2 + β2d2Y2Z 2 2 > 0. Therefore, by the Routh-Hurwitz criterion, the sub-matrix A only has negative eigenvalues. Thus, E234 is asymptotically stable if and only if 1 − cp(α)k −m d1cp(α) < 0. 3.3. Existence and stability of the positive equilibrium. By solving the system  1 −x−y = 0, k −d1y −z + β1x = 0, f(α,w) −d2z + β2y −p(α)w = 0, −m + cp(α)z = 0, (3.7) we can find the expression of the possible positive equilibrium solution E∗ = (x∗,y∗,z∗,w∗) where  x∗ = cp(α) (d1 −k) + m cp(α) (β1 + d1) , y∗ = cp(α) (β1 + k) −m cp(α) (β1 + d1) , z∗ = m cp(α) , and w∗ is given by the equation f(α,w∗) −d2z∗ + β2y∗ −p(α)w∗ = 0. Therefore, there exists a unique positive equilibrium if and only if the following inequalities hold  cp(α) (d1 −k) + m > 0 cp(α) (β1 + k) −m > 0, f0 − d2m cp(α) + cp(α)β2 (β1 + k) −β2m cp(α) (β1 + d1) > 0, (3.8) ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 201 The Jacobian of the positive equilibrium is given by J(E∗) =   −x∗ −x∗ 0 0 β1y ∗ −d1y∗ −y∗ 0 0 β2z ∗ −d2z∗ J34 0 0 cp(α)w∗ 0   , The corresponding characteristic equation is λ4 + a1λ 3 + a2λ 2 + a3λ + a4 = 0, where   a1 = d1y ∗ + d2z ∗ + x∗ > 0, a2 = −J34cp(α)w∗ + d1d2y∗z∗ + β1x∗y∗ + β2y∗z∗ + d1x∗y∗ + d2x∗z∗ > 0, a3 = −J34cp(α)d1w∗y∗ −J34cp(α)w∗x∗ + β1d2x∗y∗z∗ + d1d2x∗y∗z∗ + β2x∗y∗z∗ > 0, a4 = −J34β1cp(α)w∗x∗y∗ −J34cp(α)d1w∗x∗y∗ > 0. Calculation gives a1a2 −a3 = −J34cp(α)d2w∗z∗ + d21d2(y ∗)2z∗ + d1d 2 2y ∗(z∗)2 + β1d1x ∗(y∗)2 + β2d1(y ∗)2z∗ + β2d2y ∗(z∗)2 + d21x ∗(y∗)2 + 2d1d2x ∗y∗z∗ + d22x ∗(z∗)2 + β1(x ∗)2y∗ + d1(x ∗)2y∗ + d2(x ∗)2z∗ > 0 and a1a2a3 −a21a4 −a 2 3 = d2x ∗z∗ (cp(α)w∗J34 + β1x ∗y∗) 2 + d1d2y ∗z∗ (β1x ∗y∗ + cp(α)w∗J34) 2 + positive terms > 0. Therefore, by the Routh-Hurwitz criterion, the positive equilibrium is locally asymptotically stable. For readers’ convenience, we summarize the results obtained above about the existence and stability of all the equilibria in Table 2. As was done to rescaled model (1.5) in Section 2, we can also examine the relationship how the population size for each species at the stable positive equilibrium depends on the anti-predation level α. Indeed, we can calculate to obtain dz∗ dα = −m cp2(α) dp dα > 0. By using implicit differentiation on the system (3.7), we can also determine dx∗ dα > 0 and dy∗ dα < 0. However, we are not able to determine the sign of dw∗ dα . From the above discussion, we see that after incorporating the benefit obtained by the meso-predator’s anti-predation response in reducing the predation by the large carnivores, the final population sizes of the meso-predator, its prey and its prey’s prey are increasing, decreasing, and increasing respectively with respect to the response strength α. This alternating pattern is totally opposite to the one obtained in Section 2 on this context. 202 Y. WANG AND X. ZOU Table 2. Condition of existence and stability of the equilibria in model (3.1) Equilibrium solution Existence Stability E0 Always Unstable E1 Always Unstable E2 Always 1 −k/d1 < 0 and f0 < −β2k/d1 E3 f0 > 0 Unstable E12 k < d1 f0 < − β2k + β2β1 d1 + β1 E13 f0 > 0 md2 cp(α) > f0 > d2(k + β1) E23 d2k > f0 > −β2k/d1 f0 > d2k −d1d2 −β2 and p(α) < m (d1d2 + β2) c (kβ2 + d1f0) E34 f0 > d2m/cp(α) Unstable E123 d2(k + β1) > f0 > max ( d2k −d1d2 −β2,− β2(k + β1) d1 + β1 ) p(α) < m (β1d2 + d1d2 + β2) c (β1β2 + β1f0 + β2k + d1f0) E134 f0 > d2m/cp(α) p(α) < m c (k + β1) E234 cp(α)k −m > 0 and f0 > d2m cp(α) − β2cp(α)k −β2m d1cp(α) . 1 − cp(α)k −m d1cp(α) < 0 E∗ Condition (3.8) Stable 3.4. Numerical simulations. In this part, we present some numerical simulations to illustrate the analytical results obtained above. To this end, we choose B(α,N4) = R3 1 + c1αN4 and a34(α) = 1 1 + c2α in (1.5), leading to the following system   dN1 dτ = N1 (R1 −a11N1 −a12N2) , dN2 dτ = N2 (R2 −a22N2 −a23N3 + a21N1) , dN3 dτ = N3 ( R3 1 + c1αN4 −D3 −a33N3 + a32N2 − N4 1 + c2α ) , dN4 dτ = N4 ( −D4 + cN3 1 + c2α ) , N1(0) ≥ 0 , N2(0) ≥ 0 , N3(0) ≥ 0 , N4(0) ≥ 0, (3.9) We fix the parameters R1 = 3, a11 = 1, a12 = 0.4, R2 = 1, a22 = 0.2, a23 = 0.5, a21 = 0.5, R3 = 3, D3 = 1, a33 = 0.5, a32 = 0.05, c1 = 0.4,D4 = 0.1,c = 0.5,c2 = 0.2, (3.10) and illustrate how α impacts the population dynamics. For the above set of parameter values, k = 1/3 < d1 = 1/2, the bifurcation diagram with respect to α is given in Figure 5. There is a transcritical bifurcation between E∗ and E123 where the critical value α∗ is given by f0 − d2m cp(α∗) + cp(α∗)β2 (β1 + k) −β2m cp(α∗) (β1 + d1) = 0. Using the parameter values in (3.10), we can solve this equation to obtain α∗ ≈ 97.77777780. We can observe a trophic cascade in Figure 5 inside the coexistence region with respect to increment of ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 203 0 50 100 150 0 1 2 3 4 5 6 p o p lu a ti o n Bifurcation diagram p( )=1/(1+0.2 ) level-1 level-2 level-3 level-4 Figure 5. Bifurcation diagram of (3.9) when k < d1. 0 50 100 150 0 1 2 3 4 5 p o p lu a ti o n Bifurcation diagram p( )=1/(1+0.2 ) level-1 level-2 level-3 level-4 Figure 6. Bifurcation diagram of (3.9) when k > d1. α. Contrary to the previous section, this cascade shows an increasing population in odd level and decreasing population in even level. Next, we change R1 from R1 = 3 to R1 = 1 in (3.10) and keep the same values for other parameters. We then have k = R2/R1 = 1 and d1 = a22/a12 = 1/2 so that k > d1. For this case, we observe more complicated dynamical behaviours: there are three critical values for α, denoted by α1, α2 and α3, which are given by   cp(α1) (d1 −k) + m = 0 cp(α2) (β1 + k) −m = 0, md2 cp(α3) = f0. Using the parameter values in (3.10), we obtain α1 = 20, α2 = 70 and α3 = 95. In figure 6, when 0 < α < α1, E234 is stable; when α1 < α < α2, E ∗ is stable; when α2 < α < α3, E134 is stable; when α3 < α, E13 is stable. We can also observe trophic cascade in this case as is shown in Figure 6. 204 Y. WANG AND X. ZOU 0 5 10 15 0 1 2 3 4 5 6 p o p lu a ti o n Bifurcation diagram p( )=1/(1+2 ) level-1 level-2 level-3 level-4 (a) 0 1 2 3 4 5 0 1 2 3 4 5 6 p o p lu a ti o n Meso-predator cascade in the coexistence region level-1 level-2 level-3 level-4 (b) Figure 7. (a) Bifurcation diagram of (3.9) when c2 = 2, (b) In the region when the coexistence equilibrium E∗ is stable, we can observe a meso-predator cascade with non-monotonically change on top predator. In the last two cases, we can see that population size of large carnivores at the positive equilibrium is monotonically decreasing. We point out that this is dependent on the choice of the benefit reflecting term p(α) = a34(α)/a0. To see this, we change c2 from 0.2 to 2 (corresponding to a more significant benefit to N3 and disadvantage to N4), then, as is shown in Figure 7-(a), we can see a transcritical bifurcation from E∗ to E123 at threshold value α̃ = 10.43750000, and before this value, w ∗(α) is not monotone: it increases first and then decreases. Figure 7-(b) is an enlargement of (a) in which, one can more clearly see that while the lower trophic still follow the cascade, the top predator (w∗) increases first and then decreases with respect to the increment of α. Thus, the response function p(α) has an impact of the trophic cascade. 4. Conclusion and discussions A recent experimental study [25] in fields observed a phenomenon of trophic cascade in a food chain population system consisting of three species, i.e., meso-carnivore on top, its prey in the middle and the prey’s prey in the bottom, caused by the fear of virtual large carnivores which is implemented by playback of the large carnivores. This phenomenon, together with some recent works on fear effect in two species predator-prey models, has motivated us to theoretically explore the mechanisms for such trophic cascade in this paper. To this end, we have proposed two models, (1.3) and (1.5), with (1.3) directly corresponding to the scenario of field study in [25], and (1.5) being an extension to include a benefit in the meso-carnivore from the anti-predation response, in addition to a cost, as in (1.3). In order to incorporate the benefit term into the model, we have to add the population of the large carnivores into the interplay, making (1.5) a 4-D system. We have thoroughly analyzed the two models, using the approach of dynamical systems. For each of the two models, we have obtained complete structure of the equilibria, and established their stabil- ity/instability in terms of the model parameters, in the form of thresholds for certain parameters. For model (1.3) our results show that an anti-predation response at lower level is beneficial to the top and ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 205 bottom species (N3 and N1); while at higher level, it is beneficial to the middle species (N2) but disad- vantageous to N3 and N1, confirming the phenomenon of trophic cascade reported in that experimental study. For model (1.5), our results show that there are now three threshold values for the response level α, distinguishing ranges for α accounting for various combinations of co-existence among the four species. Particularly, within certain range of parameters, the model also demonstrates the phenomenon of trophic cascade but with an opposite alternating pattern for the three species (the meso-carnivore, its prey and its prey’s prey): increasing the response level α is beneficial to N3 and N1 but disadvantageous to N2. This change is attributed to the effect of the benefit of the anti-predation response in reducing the predation and its balance with cost of such a response. Our results can have ecological implications as they may suggest practical strategies of manage- ment/control for maintaining biodiversity. For example, in some ecosystems, populations of some meso-predators have been observed to increase significantly due to the loss/extinction of larger carni- vores, and this has in turn put some pressure on the meso-predators’ preys for their survivals. Our results on model (1.3) suggest that by creating certain virtual situations (e.g, vocal) mimicking the large carnivore predators, one may expect to reduce the populations of the meso-predators, and conse- quently relax the pressure on the meso-predators’ preys. On the other hand, if the large (top) predators of the meso-predators are present, their predation on the meso-predators poses a major threat to the meso-predators. In such a case, by out results on model (1.5), creating the aforementioned virtual situ- ations may stimulate the meso-predators to increase their anti-predation response level, which will then reduce the predation risk by the top predators. This way, the benefit of anti-predation response of the meso-predators in reducing the predation risk may outplay the cost in production, and thus, enhance the survival probability of the meso-predators. Such a net benefit in the meso-predators can then be passed on to the lower level species in an alternating fashion. Therefore, the risk events such as fear effect in some species in an ecosystem may actually offer a management tool in shifting the structure of ecosystem and help conserve the biodiversity. Note that in our model, we have used the mass action or Holling Type I functional responses as the predation mechanism. For some species between which the predation involves foraging, this mechanism is not suitable and other types of functional responses should be adopted. It would be interesting and worthwhile to investigate the population dynamics in models like (1.3) and (2.1) with such replaced functional responses. We also point out that in our models, we have only considered fear effect of meso-carnivore species N3 against the large carnivores N4. Such fear effect may also exist in N2 against N3 and in N1 against N2. Modeling fear effect in those or in all levels would also be interesting but would be very challenging mathematically. We remark that for predator-prey interactions between two species only, the recent works mentioned in the introduction may also suggest some possible extensions and expansions of the two models in this paper. For example, one may also incorporate age structure, spatial structure, digestion delay, extra food, stochastic noise, as was done in [8, 18, 31, 32, 33]. Efforts on all these lines will greatly enhance our understanding of predator-prey interactions, and enrich the theory in this area. References [1] R. Boonstra, D. Hik, G.R. Singleton and A. Tinnikov, The impact of predator-induced stress on the snowshoe hare cycle, Ecol. Monogr. 68(1998) , 371-394. [2] C. Castillo-Chavez and H. R. Thieme, Asymptotically autonomous epidemic models, In: O. Arino et al: Mathematical Population Dynamics: Analysis of Heterogeneity, I. Theory of Epidemics (pp. 33-50). Wuerz, Canada, 1995 [3] M. Clinchy, L. Y. Zanette, R. Boonstra, J. C. Wingfield and J. N. M. Smith, Balancing food and predator pressure induces chronic stress in songbirds, Proc Biol Sci. 271(2004), 2473-2479. 206 Y. WANG AND X. ZOU [4] S. Creel, J. Winnie Jr., B. Maxwell, K. Hamlin and M. Creel.,Elk alter habitat selection as an antipredator response to wolves, Ecology 86(2005), 3387-3397. [5] S. Creel, D. Christianson, S. Liley and J. A. Winnie, Predation risk affects reproductive physiology and demography of elk, Science 315(2007), 960. [6] S. Creel and D. Christianson, Relationships between direct predation and risk effects, Trends Ecol Evolut 23(2008), 194-201. [7] W. Cresswell, Predation in bird populations, J Ornithol 152(2011), 251-263. [8] D. Das and G. P. Samanta, Modeling the fear effect on a stochastic prey-predator system with additional food for the predator, J. Phys. A: Math. Theor. 51 (2018): 465601 (26pp). [9] H. Freedman and P. Waltman, Persistence in models of three interacting predator-prey populations, Math. Biosci. 68 (1984), 213-231. [10] G. F. Gauze, The Struggle for Existence, The Williams and Wilkins company, 1934. [11] A. Hastings and T. Powell, Chaos in a three-species food chain, Ecology, 72(3) (1991), 896-903. [12] S. Hsu, S. Ruan and T. Yang, Analysis of three species Lotka-Volterra food web models with omnivory, J. Math. Anal. Appl. 426 (2015), 659-687. [13] N. Krikorian, The Volterra model for three species predatorprey systems: boundedness and stability, J. Math. Biol. 7 (2) (1979), 117-132. [14] S. L. Lima, Nonlethal effects in the ecology of predator-prey interactions, Bioscience 48(1998), 25-34. [15] S. L. Lima, Predators and the breeding bird: behavioural and reproductive flexibility under the risk of predation, Biol Rev 84(2009), 85-513. [16] A. J. Lotka, Analytical Note on Certain Rhythmic Relations in Organic Systems, Proc. Natl. Acad. Sci. U.S.A. 6(1920), 410-415 . [17] A. J. Lotka, Elements of Physical Biology, Williams and Wilkins (1925). [18] S. Mondal, A. Maiti and G. P. Samanta, Effects of Fear and Additional Food in a Delayed Predator-Prey Model, Biophysical Reviews and Letters, 13 (2018), 157-177. [19] P. Panday, N. Pal, S. Samanta and J. Chattopadhyay, Stability and bifurcation analysis of a three-species food chain model with fear, Int. J. Bifurc. Chaos 28 (2018): 1850009. [20] P. Panday, N. Pal, S. Samanta and J. Chattopadhyay, A Three Species Food Chain Model with Fear Induced Trophic Cascade, Int. J. Appl. Comput. Math 5 (2019): 100. [21] B. K. Peckarsky, C.A. Cowan, M.A. Penton and C. Anderson,Sublethal consequences of stream-dwelling predatory stoneflies on mayfly growth and fecundity, Ecology 74(1993), 1836-1846. [22] O. J. Schmitz, A. P. Beckerman and K. M. OBrien. Behaviorally mediated trophic cascades: effects of predation risk on food web interactions, Ecology 78(1997), 1388-1399. [23] O. J. Schmitz, V. Krivan and O. Ovadia. Trophic cascades: the primacy of trait-mediated indirect interactions, Ecology Letters 7(2004), 153-163. [24] H. L. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, Amer. Math. Soc. Math. Surveys and Monographs Vol. 41, 1995. [25] J. P. Suraci, M. Clinchy, L. M. Dill, D. Roberts and L. Zanette, Fear of large carnivores causes a trophic cascade, Nat Commun. 7(2016): 10698. [26] R. Tollrian and C. D. Harvell, The evolution of inducible defenses: current ideas, in The Ecology and Evolution of Inducible Defenses (Tollrian, R. and Harvell, C.D., eds), Princeton University Press, (1999), 306-321. [27] M. Valeix, A. J. Loveridge, S. Chamaille-Jammes, Z. Davidson, F. Murindagomo, H. Fritz, and D. W. Macdonald, Behavioral adjustments of African herbivores to predation risk by lions: spatiotemporal variations influence habitat use, Ecology 90(2009), 23-30. [28] V. Volterra, Variazioni e fluttuazioni del numero d’individui in specie animali conviventi, Mem. Acad. Lincei Roma. 2(1926), 31-113. [29] V. Volterra, Variations and fluctuations of the number of individuals in animal species living together, in Chapman, R. N. Animal Ecology. McGraw-Hill, 1931. [30] X. Wang, L. Y. Zanette and X. Zou, Modelling the fear effect in predator-prey interactions, J Math. Biol. 73(2016 ), 1179-1204. [31] X. Wang and X. Zou, Modeling the Fear Effect in Predator-Prey Interactions with Adaptive Avoidance of Predators, Bull Math Biol 79(2017), 1325-1359. [32] X. Wang and X. Zou, Pattern formation of a predator-prey model with the cost of anti-predator behaviours, Math. Biosci. Eng. 15(2018), 775-805. ANTI-PREDATION RESPONSE IN FOOD CHAIN SYSTEMS 207 [33] Y. Wang and X. Zou, On a predator prey system with digestion delay and anti-predation strategy, J Nonlinear Sciences, 30(2020), 1579-1605. [34] L. Y. Zanette, Synergistic effects of food and predators on annual reproductive success in song sparrows, Proc. Biol. Sci. 270(2003), 799-803. [35] L. Y. Zanette, A. F. White, M. C. Allen and M. Clinchy, Perceived predation risk reduces the number of offspring songbirds produce per year, Science 334 (6061)(2011), 1398-1401. Corresponding author. Department of Applied Mathematical Sciences, University of Western Ontario London, ON, Canada N6A 5B7 E-mail address: ywan342@uwo.ca Department of Applied Mathematical Sciences, University of Western Ontario London, ON, Canada N6A 5B7 E-mail address: xzou@uwo.ca