Schooling behavior driven complexities in a fear-induced prey–predator system with harvesting under deterministic and stochastic environments

In a region under consideration, let at any instant (t>0), x and y represent the prey and predator population densities, respectively. The rate of change of each model species density at time t is made on the following assumptions:

1. 1.

   Prey population grow logistically in the absence of predator with birth rate r, which is affected by the fear ((f_1)) of predator (when predators are around).

2. 2.

   There is a reduction in the rate of prey density change due to three types of death, namely, natural death with the rate (d_1), fear related death⁵ with the level of fear (f_2) and over crowding death with the rate (d_2).

3. 3.

   Also, the rate of change of prey density decreases due to predation of predator population following a predator-dependent functional response describing both predatory and prey schooling behaviors¹⁰. Response function is expressed in functional form describing as (zeta (x, y)=frac{cxy}{1+chxy}), where c denotes the rate of consumption and h represents handling time of predator for one prey.

4. 4.

   Predator population survive in the system by consuming prey population only. They grow with conversion efficiency (c_1) of prey biomass into predator biomass.

5. 5.

   Predator population harvested from the system which reduces its rate of density. We consider a nonlinear harvesting term (Michaelis-Menten type) given by, (H(y)=dfrac{qEy}{p_1E+p_2y}). Here, parameters q and E, respectively, represent the catchability rate and harvesting effort. It is easy to observe that (Hrightarrow frac{q}{p_1}y) as (Erightarrow infty) for a fixed value of y. Also, (Hrightarrow frac{q}{p_2}E) as (yrightarrow infty) for a fixed value of E. Therefore, at higher effort levels, (p_1) is proportional to the stock level-catch rate ratio and at higher levels of stock, (p_2) is proportional to the effort level-catch rate ratio.

6. 6.

   Lastly, we assume that the predator population experience natural as well as over crowding related death with the rates (d_3) and (d_4), respectively.

Keeping all these above assumptions in mind, we formulate the following prey–predator model:

System (1) is to be analyzed with the initial conditions (x(0),y(0)>0). All the model parameters are assumed to be positive constants and their hypothetical values that we used for numerical calculations are as follows:

In Table 1, we have provided system’s equilibria, sufficient conditions of their existence and stability. Mathematically, it is difficult to determine the existence of coexistence (interior) equilibrium point(s) by the given nullclines. So, we visualize it numerically (see Fig. 1). It is apparent from the figure that on increasing the value of E, number of coexistence equilibrium points reduces and after a certain range there is no coexistence equilibrium point.

Transcritical bifurcation

From Table 1, it is clear that the equilibrium (E_0) is stable if (r<d_1), which is opposite to the existence condition of (E_1). That is, equilibrium (E_0) is stable whenever the equilibrium (E_1) does not exist, and hence these two equilibria are related via transcritical bifurcation. Using Sotomayor theorem³⁸, we can easily prove that model system (1) experiences transcritical bifurcation at the trivial equilibrium point (E_0) as the growth rate of prey crosses a critical value (r^{[TB]}=d_1).

We visualize the transcritical bifurcation graphically in Fig. 2. It is clear from the figure that when the value of r is less than (d_1=0.1), the equilibrium point (E_0) only exists and it is stable. If we increase the value of r ((r>r^{TB}=d_1=0.1)) then (E_0) becomes unstable and the equilibrium point (E_1) exists and becomes stable.

Hopf bifurcation

One of the most common dynamics in interacting population dynamics is oscillating behavior, which implies that there is a Hopf bifurcation. By local changes in equilibrium properties, Hopf bifurcation describes when a periodic solution appears or disappears. In this section, we study the Hopf bifurcation through the coexistence equilibrium (E^*) with respect to the model parameter E. Discussion for the existence of Hopf bifurcation is as follows:

As it is easy to follow, we verify Hopf bifurcation numerically. We have considered the parameters value same as in (2) except (c=0.1) and E. At (E=E^{[HB]}=0.1196559641), the trace of the Jacobian matrix at (E^*(2.618402886, 2.352228027)) is zero and determinant, (Det(J_{E^*})=0.4474794791>0). The value of (dfrac{d(Tr(J_{E^*}))}{dE}Big |_{E=E^{[HB]}}=-0.02965188514ne 0). Therefore, the transversality condition for Hopf bifurcation is also satisfied at (E=E^{[HB]}). Thus, these results confirm that the system (1) experiences a Hopf bifurcation² around (E^*(2.618402886, 2.352228027)).

Moreover, we obtain Lyapunov number (L_1=-0.04728284756pi <0) at (E^{[HB]}=0.1196559641). This implies that system (1) goes through a supercritical Hopf bifurcation at (E^*(2.618402886, 2.352228027)).

Saddle-node bifurcation

Sotomayor’s theorem says that saddle-node bifurcation may occur at coincident equilibrium points depending on threshold value of the bifurcation parameters. Since the analytical result is difficult to follow, we verify the existence of saddle-node bifurcation of system (1) numerically for given set of parameters value (2) except (c=0.15) and E. At (E=E^{[SN]}=2.6670276565), system (1) has a coincident equilibrium (E^*(12.07536421, 1.327137944)). The corresponding Jacobian matrix (J_{E^*}(E=E^{[SN]})= left( begin{matrix} -0.9247512109 &{} -10.89563235 0.08585790896 &{} 1.011635174 end{matrix} right)) has one zero eigenvalue. For the matrices (J_{E^*}) and (J_{E^*}^T) , the eigenvectors V and W associated with the zero eigenvalue can be found as follows: (V=[-7.625518027 quad 0.6469271584]^T) and (W=[-0.06805575397 quad -0.7326946399]^T). This yields,

Based on Sotomayor’s theorem³⁸, we can conclude that system (1) experiences a saddle-node bifurcation at the equilibrium point (E^*) when the parameter E crosses the critical value (E=E^{[SN]}=2.6670276565).

Bogdanov–Takens bifurcation

A number of one-dimensional bifurcations such as transcritical, Hopf and saddle-node have been studied in earlier segments. Each of these bifurcations belongs to one parametric bifurcation. In the present segment, we will discuss two parametric bifurcation, i.e., of codimension two bifurcation. Bogdanov–Takens bifurcation (BT) is a such type of bifurcation. When Hopf and saddle-node bifurcation curves collide, a bifurcation of this type occurs at the vicinity of colliding point. In this case, the critical values of two bifurcation parameters give zero eigenvalues of multiplicity two in the Jacobian matrix of the system.

One can used Kuznetsov’s³⁹ technique to achieve the standard form of BT-bifurcation. Due to model complexity, we did not find any explicit expression of BT-bifurcation. So, we verify it numerically, which is discussed below.

System (1) undergoes BT-bifurcation for the set of parameters given in (2) except c and E. We consider c and E as bifurcation parameters and observe that saddle-node and Hopf bifurcation curves collide with each other at ((E, c)=(1.42829255, 0.105462715)=left( E^{[BT_1]}, c^{[BT_1]}right),) and (,(E, c)=(7.345497850, 0.2137076372)=left( E^{[BT_2]}, c^{[BT_2]}right)), and the instantaneous equilibrium points are (11.12253608, 1.485272279) and (13.778434, 1.13304881), respectively. Therefore, system (1) undergoes through two BT-bifurcation points: one for (E^*=(11.122536, 1.4852723)) at (left( E^{[BT_1]}, c^{[BT_1]}right) =(1.42829255, 0.105462715)) and another for (E^*=(13.778434, 1.13304881)) at (left( E^{[BT_2]}, c^{[BT_2]}right) =(7.345, 0.213)).

Numerical results of the deterministic system (1)

To explore the rich dynamics of the system (1), we perform extensive numerical simulations in this section. For numerical simulations, we choose a set of biologically feasible hypothetical parameter values as given in (2).

A system may exhibit rich dynamical behavior if its dynamics are investigated in different bi-parameter spaces. So, here we plot two bi-parametric bifurcations (Figs. 3a, 5a). In Fig. 3a, we plot saddle-node, Hopf and Bogdanov–Takens bifurcations in (E) (-) (c) plane, which divide the entire region into six different regions ((R_1) (-) (R_6)) of different dynamical behaviors. Complete phase portraits are drawn for every region in order to understand their dynamics transparently, Fig. 3b–i. It is apparent from Fig. 3a that the saddle-node curve divides the whole region into two parts. Among these two parts one part is region (R_1), which has no interior equilibrium point (see Fig. 3b) and other part is sum of rest regions ((R_2) (-) (R_6)), that contain two interior equilibrium points (see Fig. 3c–i). Next, for lower and higher values of E, there exist Hopf bifurcation curves, which divide two portions of the region into two subregions that are ((R_2)/(R_3)) and ((R_5)/(R_6)). In one side of the Hopf curve, one of the two interiors is stable spiral and other one is saddle, Fig. 3c,i (correspond to the regions (R_2) & (R_6), respectively). On the other hand, as the value of E and c crosses Hopf curve, stable spiral point becomes unstable spiral (stable limit cycle) and saddle one remain same, Fig. 3d,h (correspond to the regions (R_3) & (R_5), respectively). Similar as Hopf curve, there exist two Homoclinic curves for lower and higher values of E, which divide a subregion into two another subregions ((R_3)/(R_4)) and ((R_4)/(R_5)). As the values of (E, c) closest to the Homoclinic curve system gives larger amplitude-limit cycle and on the curve it gives maximum amplitude-limit cycle, Fig. 3e,g. Further, if the values of (E, c) crosses the Homoclinic curve, stable limit cycle destroys and the equilibrium point becomes unstable (see Fig. 3f corresponds to the region (R_4)). Also, it is clear from Fig. 3a that there exist two BT-points ((BT_1=(1.428292550, 0.1054627152)) and (BT_2=(7.345497850, 0.2137076372))), where aforementioned three bifurcation curves meet with each other. Therefore, the neighbourhood of BT-points exhibits complex dynamics of the system.

It is observed that the system preserves bistability for the regions (R_2), (R_3), (R_5) and (R_6). For a fixed set of parameter values, depending on the initial population size the trajectories initiating inside the green regions converge to the stable interior (stable limit cycle) whereas the trajectories initiating in the red regions converge to the predator-free equilibrium point, Fig. 4. In Fig. 5a, we plot another bi-parametric bifurcation in (f_1) (-) (q) plane, which is divided into four subregions ((R_1) (-) (R_4)) of different dynamics by different bifurcation curves. Phase portraits of every region are plotted (see Fig. 5b–f). We observe that there is no interior for diagonally higher values of both of the parameters (f_1) and q whereas two interiors exist when the values of (f_1) and/or q are lower, Fig. 5a. Detailed discussions of the figure are same as previous one (Fig. 3).

Next, to observe the effect of important model parameters explicitly, we plot one parameter bifurcation diagrams while other parameters are fixed, Fig. 6. It is clear form Fig. 6a,d,e that for lower values of the parameters r, c and (c_1), there is no interior. As the value of these parameters surpasses saddle-node bifurcation point there exist two interior one of which always remain saddle and another one switches its stability (stable spiral (Rightarrow) unstable spiral) through Hopf bifurcation point. Therefore, up to a certain range of these parameters have destabilizing effect. On the other hand, for lower values of the parameters (f_1), (f_2) and E there exist two interior one of them is always saddle and another one switches its stability (unstable spiral (Rightarrow) stable spiral) through Hopf bifurcation point (see Fig. 6b,c,f). But, when the value of these parameters crosses saddle-node bifurcation point, no interior equilibrium point exists. Note that up to a certain range of these parameters, (f_1) and (f_2) have stabilizing effect whereas E has both destabilizing as well as stabilizing effects.

Interpretation of the deterministic results in the context of biology

For the deterministic system, we mainly investigate different types of bifurcation results. They are Transcritical, Hopf, Saddle-node and Bogdanov–Takens. In the context of a biological system, a bifurcation occurs when a small smooth change made to the parameter values (the bifurcation parameters) of a system causes a sudden ’qualitative’ or topological change in the behavior of the system. Bifurcations mainly describe changes in the stability and/or existence of fixed points (equilibrium points) and the bifurcation parameters act as a control parameter. Changes in the control parameter eventually changed the qualitative behavior of the system. Above mentioned bifurcations of the considered system (1) were found for the model parameters r, (f_1), (f_2), c, (c_1), q and E. Therefore, all these ecological parameters (i.e., growth rate of prey, fear, capture rate, harvesting effort) act as control parameters of the proposed system. In Transcritical bifurcation, there is a critical value of the bifurcation parameter from which two equilibrium points exchange their stability. Moreover, on the one side of the critical value one of the two equilibrium points exist while on the other side both equilibrium points exist. Therefore, we have a critical value of an ecological parameter from which species persistency and stability can be described. In case of Hopf bifurcation, there is a critical value of a parameter from which densities of model species either fluctuate from fixed stable densities or vice versa. In the case of Saddle-node bifurcation, coexistence equilibrium points exist or diminish in pair by varying bifurcation parameter i.e., from the critical value of the bifurcation parameter coexistence equilibrium points exists or diminish in pair. Lastly, BT-bifurcation, the point where saddle-node curve and Hopf curve meet is known as BT-point. Thus, all these phenomena occur around the BT-point, i.e., near the BT-point complex dynamics (like species persistency, stability, extinction) occurs. Thus, studying the bifurcation in prey–predator systems has great importance from environmental perspective. In order to maintain an ecological balance between prey and predator populations, identification of bifurcation parameters plays a crucial role in determining an effective control strategy.

Source: Ecology - nature.com