Complex dynamics near extinction in a predator-prey model with ratio dependence and Holling type III functional response

In this paper, we analyze a recently proposed predator-prey model with ratio dependence and Holling type III functional response, with particular emphasis on the dynamics close to extinction. By using Briot-Bouquet transformation we transform the model into a system, where the extinction steady state is represented by up to three distinct steady states, whose existence is determined by the values of appropriate Lambert W functions. We investigate how stability of extinction and coexistence steady states is affected by the rate of predation, predator fecundity, and the parameter characterizing the strength of functional response. The results suggest that the extinction steady state can be stable for sufficiently high predation rate and for sufficiently small predator fecundity. Moreover, in certain parameter regimes, a stable extinction steady state can coexist with a stable prey-only equilibrium or with a stable coexistence equilibrium, and it is rather the initial conditions that determine whether prey and predator populations will be maintained at some steady level, or both of them will become extinct. Another possibility is for coexistence steady state to be unstable, in which case sustained periodic oscillations around it are observed. Numerical simulations are performed to illustrate the behavior for all dynamical regimes, and in each case a corresponding phase plane of the transformed system is presented to show a correspondence with stable and unstable extinction steady state.

Introduction

A starting point for many ecological models analyzing predator-prey interactions is provided by models of Gause-Kolmogorov type [1–3]

$\begin{array}{l}\dot{u}=uf(u)-vg\left(u,v\right),\\ \dot{v}=bvg\left(u,v\right)-dv,\end{array}$

where u(t) and v(t) are population densities or abundances of prey and predator, d is the natural death rate of predators, and functions f(u) and g(u, v) describe, respectively, natural per capita growth rate of prey, and the trophic function or functional response of predators [4], which quantifies how efficiently predators are consuming prey, and how prey consumption increases predator reproduction. Function f(·) is often chosen in the form of either a constant (so-called Malthusian growth), or as a monotonically decreasing linear function f(u) = r(1−u/K), which describes intra-specific competition of prey for available resources, and is known as the logistic growth of prey with a linear growth rate r and the carrying capacity K. In the simplest case, the function g(u, v) would only depend on the abundance of prey g(u, v) = g(u), in which case it is known as prey-dependent functional response. Perhaps, the simplest example of prey-dependent functional response is when it is proportional to the number of prey g(u) = au, which yields the classical Lotka-Volterra model [5, 6]. A more realistic representation of interactions between prey and predators is given by a functional response suggested by Holling [7, 8], which accounts for the time predators spend searching for prey, as well as handling (i.e., chasing, killing, and digesting) the prey. Holling proposed three types of functional response g(u), which all satisfy g(0) = 0 and approach some constant value for large u. The difference between different types of Holling functional response is in the behavior for small prey numbers/densities. Type I response is linearly increasing for small prey densities, whereas for large prey numbers it saturates at some constant value; type II and type III are functions that are also saturating at high prey numbers, and are, respectively, concave and sigmoidal. Besides direct impact of predation, prey behavior can also be affected by the fear of predation, which has been explored in a series of recent papers [9–13].

By the early 1990s, a number of ecological studies showed that when studying predator-prey interactions across a variety of scales and biological species, from insects to mammals [14–17], it may be more realistic to consider the functional response that depends on the ratio of prey to predators, which became known as ratio-dependent functional response [18, 19]. Such functional response can effectively model an observation that at higher predator densities, predators would have to compete for and share some of the prey. Formally, this would correspond to writing the functional response g(u, v) in the form g(u/v), which would then yield the following model

$\begin{array}{l}\dot{u}=uf(u)-vg\left(\frac{u}{v}\right),\\ \dot{v}=bvg\left(\frac{u}{v}\right)-dv.\end{array}$

A number of theoretical studies have studied the dynamics of models with ratio dependence [20–26] and various forms of functional response that have usually taken the form of g(z) = a, g(z) = a₁z/(1+a₂z), and $g(z)=a_1{z}^n/(1+a_2{z}^n)$ with n>1, where z = u/v, for Holling type I, II and III responses, respectively [27, 28]. Some additional ecological effects that have been studied in models with ratio dependence include time delays [29–31], as well as spatial dependence [32–34]. One important observation to make here is that due to the very nature of ratio-dependent functional response, the origin, i.e., the point (u, v) = (0, 0) that is characterized by the extinction of both species, may create certain mathematical challenges for analysis from the perspective that either the vector field itself is not defined at that point, or it is defined, but the linearisation is not, which hinders standard stability analysis. To overcome these challenges, one approach that is often used consists in rescaling time with the denominator of functional response, often in combination with Briot-Bouquet transformation to remove the singularity [22, 24, 26].

In this paper we are interested in Holling type III functional response with ratio dependence that has been observed in a number of experimental settings [35–38]. Rather than using a functional form $g(z)=a_1{z}^n/(1+a_2{z}^n)$ given by a ratio of two polynomials, where z = u/v and n>1, we consider an alternative form that can be written as follows,

$\begin{array}{ll}g(z)=a{e}^{-\alpha /z},& (1)\end{array}$

which satisfies the conditions of g(0) = 0, is monotonically increasing, and settles at a constant value as z → ∞. Biological motivation for this form of functional response comes from the studies of plant-parasite interactions, in which plants serve as hosts for insect predators, and hence, current plant population represents a carrying capacity for insect population, in which case the interaction term has the form g(u, v) = a(1−v/u) [39, 40]. Formally, trophic function (Equation 1) is reminiscent of the Ivlev trophic function g(z) = a(1−e^−αz) [41] that represents Holling type II response, as well as of the Ricker model N_t+1 = N_texp[r(1−N_t)] for single-species populations [42]. Making ratio dependence in Equation (1) explicit, we have the functional response g(u, v) being given by

$\begin{array}{ll}g(u,v)=a{e}^{-\alpha v/u}.& (2)\end{array}$

A recent work has used a delayed version of this trophic function to analyze the dynamics of vector-plant interactions in the context of modeling plant mosaic disease [43], while the role of stochastic effects in the time-delayed model has been studied in Blyuss et al. [44]. Of particular interest to us is the analysis of model behavior in the neighborhood of extinction (u, v) = (0, 0), and the effects this has on overall dynamics. While earlier numerical simulations have suggested that extinction is indeed possible in some parameter regimes [44], the question of when exactly this happens, and whether extinction can coexist with other states has not been explored.

The remainder of this paper is organized as follows. In the next section we identify different steady states of the model and show that the origin that corresponds to extinction of both species is a well-defined steady state. Applying Briot-Bouquet transformation, we will then obtain another version of our model, where the original extinction steady state is unfolded as an entire axis in the phase space. Whereas, the origin of the modified model is always unstable, there are two more steady states that also correspond to the extinction steady state of the original model, and whose locations are given by a Lambert W function. We will derive analytical conditions for stability of these steady states and will illustrate regions of their feasibility and stability depending on model parameters. Model dynamics is further explored by numerically computing regions of feasibility and stability for different steady states and periodic orbits of the model depending on predation rate, predator growth rate, and the parameter characterizing functional response. We also demonstrate phase plane of the modified model in each scenario, clearly indicating different steady states that correspond to extinction and coexistence, and complement these by numerical solutions of the original model, which can exhibit such distinct types of behavior, as extinction of both species; survival of prey only, when predator growth rate is not sufficiently high; a regime of bi-stability between these two scenarios, where for the same values of parameters it is the initial values that determine whether both species or only predators will go extinct. Other possibilities include a stable coexistence steady state, where both prey and predators are maintained at some constant level, suggesting that predation is compensated or balanced by the prey growth, a periodic solution around the coexistence steady state, which is reminiscent of oscillations in the standard Lotka-Volterra model, and a regime of bi-stability between a stable coexistence and extinction. The paper concludes with a discussion of results.

Methods

As discussed in the Introduction, we consider a predator-prey model with logistic growth of prey and a Holling type III functional response with ratio dependence, which has the form

$\begin{array}{l}\dot{u}=ru\left(1-\frac{u}{K}\right)-av{e}^{-\alpha v/u},\\ \dot{v}=bv{e}^{-\alpha v/u}-dv.& (3)\end{array}$

In order to simplify the model and to reduce the number of free parameters, we rescale the variables and parameters as follows,

$\begin{array}{c}u=K\hat{u},v=K\hat{v},rt=\hat{t},\\ \\ \frac{d}{r}=\hat{d},\frac{a}{r}=\hat{a},\frac{b}{r}=\hat{b}.\end{array}$

The model (3) can then be rewritten in the form

$\begin{array}{l}\dot{u}=u(1-u)-av{e}^{-\alpha v/u},\\ \dot{v}=bv{e}^{-\alpha v/u}-dv,& (4)\end{array}$

where we have dropped hats for notational convenience. It is straightforward to show that this system is well-posed in that for arbitrary non-negative initial conditions, its solutions will remain non-negative and bounded for all t≥0. For any values of parameters, the system (Equation 4) has a prey-only steady state E = (1, 0), which is stable for b<d and unstable for b>d [44]. If b>d and αb−adln(b/d)>0, this system also has a coexistence steady state E^* = (u^*, v^*), with

$\begin{array}{l}{u}^{*}=\frac{1}{\alpha b}[\alpha b-ad\ln (\frac{b}{d})],\hfill \\ v^{*}=\frac{1}{\alpha b}\ln (\frac{b}{d})[\alpha b-ad\ln (\frac{b}{d})].\hfill \end{array}$

This steady state is stable, provided the following condition is satisfied [44]:

$-1+dln(\frac{b}{d})\left[\frac{a}{\alpha b}(2-ln(\frac{b}{d}))-1\right]<0.$

In an earlier paper [44], we focused on the analysis of coexistence steady state E^* = (u^*, v^*) under effects of maturation time delay in predators, and also investigated the role of stochasticity.

In the neighborhood of the point (u, v) = (0, 0), in the first quadrant the term ve^−αv/u is well-defined and positive, and in the limit (u, v) → (0, 0), we have

$0\le \left|v{e}^{-\alpha v/u}\right|\le \left|v\right|{\displaystyle \underset{(u,v)\to (0,0)}{\overset{}{\to }}}0,$

which shows that (0, 0) is indeed a steady state of the system (Equation 4), which biologically describes extinction of both species. Since the system (Equation 4) is not differentiable at the point (0, 0), to explore its dynamics close to this steady state, we employ the Briot-Bouquet transformation [45] that in our case consists in introducing a new variable z, such that v = zu, and the system (Equation 4) then transforms into

$\begin{array}{l}\dot{u}=u(1-u)-azu{e}^{-\alpha z},\\ \dot{z}=z(b+az)e^{-\alpha z}-z(d+1-u).& (5)\end{array}$

The Briot-Bouquet transformation blows up the origin into an entire z-axis [46]. Hence, we should investigate equilibria of the system (Equation 5) that lie on the z-axis.

For any parameter values, the transformed system (Equation 5) has a steady state E₀ = (0, 0). Besides this steady state, as long as the following condition is satisfied

$\begin{array}{ll}\frac{\alpha (d+1)e^{-b\alpha /a}}{a}\le \frac{1}{e},& (6)\end{array}$

the system (Equation 5) can have up to two more steady states E₁ = (0, z₁) and E₂ = (0, z₂), where

$\begin{array}{l}{z}_1=-\frac{b}{a}-\frac{1}{\alpha }{W}_0[-\frac{\alpha (d+1)e^{-b\alpha /a}}{a}],\hfill \\ z_2=-\frac{b}{a}-\frac{1}{\alpha }{W}_{-1}[-\frac{\alpha (d+1)e^{-b\alpha /a}}{a}],\hfill \end{array}$

and W₀(·) and W₋₁(·) are two branches of Lambert W function. Figure 1 illustrates the dependence of these two functions on their argument, and we note that the smallest possible value of W₀(x) is attained at x = −1/e and is equal to W₀(−1/e) = −1, which suggests that the steady state E₁ is only biologically feasible when 0 < αb/a < 1. Furthermore, since W₋₁(x) ≤ W₀(x) < 0 for x < 0, this implies that whenever both E₁ and E₂ exist, we have z₁ ≤ z₂.

FIGURE 1

Figure 1. Lambert functions W₀(x) and W₋₁(x).

Linearisation of the system (Equation 5) near any steady state $\hat{E}=(\hat{u},\hat{z})$ gives the Jacobian

$J=\left|\begin{array}{l}1-2\widehat{u}-a\widehat{z}{e}^{-\alpha z}-a\widehat{u}{e}^{-\alpha \widehat{z}}(1-\alpha \widehat{z})\hfill \\ \widehat{z}(b+a\widehat{z})(1-\alpha \widehat{z})e^{-\alpha \widehat{z}}+az{e}^{-\alpha \widehat{z}}-(d+1-\widehat{u})\hfill \end{array}\right|.$

At E₀, the two eigenvalues are 1 and b−d+1, suggesting that this steady state is always unstable and is either a saddle, or an unstable node. At the steady states E₁ and E₂, we have the Routh-Hurwitz conditions for stability in the form

$1<a\hat{z}{e}^{-\alpha \hat{z}}<\alpha (d+1)\hat{z}.$

Results

In order to illustrate the effects of different parameters on feasibility and stability of the two steady states E₁ and E₂ that correspond to extinction in the original model, we plot in Figure 2 maximum real part of the characteristic eigenvalues as determined by the Jacobian (Equation 7). We observe that the steady state E₁ only exists for sufficiently small levels of predator growth rate b, and either low enough predation rate a, or sufficiently small rate α For any parameter combination, where this steady state is biologically feasible, it is unstable. In contrast, starting with rather small values of b, the steady state E₂ is biologically feasible for an entire range of values of a or α and is stable for smaller values of a, or for larger values of α, and the region of stability shrinks with increasing predator growth rate b.

FIGURE 2

Figure 2. Regions of feasibility and stability of steady states E₁ (A,C) and E₂ (B,D). Steady states are feasible in a colored region, and biologically infeasible in the white region. Color code denotes the largest real part of characteristic eigenvalues. Parameter values are a = d = 1 (A,B), α = d = 1 (C,D).

Under the inverse Brio-Bouquet transformation, we can interpret regions of stability of the steady state E₂ of the transformed system (Equation 5) as parameter regions, where the extinction steady state E₀ = (0, 0) of the original system (Equation 4) is stable. Combining these results with those on stability of the prey-only equilibrium E and the coexistence equilibrium E^*, we now illustrate in Figure 3 an entire range of dynamical scenarios that can be exhibited in the model.

FIGURE 3

Figure 3. Regions of feasibility and stability of different steady states of the model (Equation 4). The prey-only steady state E = (1, 0) is stable to the left of the dashed black line at b = 1, and unstable to the right of that line. The extinction steady state E₀ is stable above red line in plot (A) and below red line in plot (B). In region 1, there is bi-stability between a stable steady state E and the extinction steady state E₀; in region 2 only the steady state E is stable; in region 3, only the coexistence steady state E^* is stable; in region 4, the coexistence steady state E^* is unstable, and there is periodic orbit around this state; in region 5, there is bi-stability between the extinction steady state E₀ and the coexistence steady state E^*; in region 6, only extinction steady state E₀ is stable. Parameter values are α = d = 1 (A), a = d = 1 (B).

For values of b smaller than 1, which biologically describes a scenario, where the rate of growth in the population of predators is too small to be able to maintain their population as based on the available prey, the coexistence steady state E^* is not biologically feasible, while the prey-only steady state E is stable. In this case, the system either goes to this prey-only steady state as the only stable equilibrium for parameter values in region 1, or it exhibits a bi-stability between this steady state, and a stable extinction steady state E₀ in region 2. This is shown in detail in Figure 4, where, for convenience, for each parameter combination, in the right column we plot numerical solution of the model (4), and in the left column we plot the corresponding phase plane of the transformed system (5) to demonstrate how changes in stability of the extinction steady state E₀ affect global dynamics. Since in region 1 there is bi-stability between two steady states, both of which are characterized by the absence of predators, this suggests that it is rather the initial conditions of the system that determine whether prey will survive (though, initially, prey population will decline due to predation), or will also be driven to extinction. Interestingly, even though in this parameter regime, predator population goes to zero, this does not automatically guarantee the survival of the prey, which is affected not only by the values of parameters, but also by the initial presence of sufficiently many prey to maintain its stable population.

FIGURE 4

Figure 4. Phase planes (left column) of the transformed system (Equation 5) and the corresponding numerical solution (right column) of the rescaled original system (Equation 4), illustrating different dynamical regimes listed in Figure 3A. Blue/red dots in phase plane plots indicate stable/unstable steady states. Parameter values are α = d = 1. (A,B) Case 1, a = 6, b = 0.5, bi-stability between prey-only steady state E and the extinction steady state E₀. (C,D) Case 2, a = 3, b = 0.5, prey-only steady state E is the only stable state. (E,F) Case 3, a = 1.5, b = 1.4, only the coexistence steady state E^* is stable.

As the value of b exceeds 1, the prey-only equilibrium E becomes unstable, while the coexistence steady state E^* becomes feasible and stable, as shown in Figures 4E,F (region 3). As the value of predation rate a increases, as long it stays below the boundary of stability of the extinction steady state E₀, the coexistence steady state E^* loses its stability through a supercritical Hopf bifurcation, giving rise to stable periodic solutions around this steady state, as illustrated in Figures 5A,B (region 4). Once the stability boundary of the extinction steady state is crossed, for smaller values of b and sufficiently high values of predation rate a, stable coexistence steady state E^* co-exists with a stable extinction steady state E₀ shown in Figures 5C,D (region 5). Biologically, this means that for the same parameter values, whether or not the system evolves toward coexistence of prey and predators, or they drive each other to extinction, is again determined by the choice of initial conditions. In region 6, which corresponds to high values of both predation rate and the growth rate of predators, the only stable state of the system is that of extinction, as illustrated in Figures 5E,F.

FIGURE 5

Figure 5. Phase planes (left column) of the transformed system (Equation 5) and the corresponding numerical solution (right column) of the rescaled original system (Equation 4), illustrating different dynamical regimes listed in Figure 3A. Blue/red dots in phase plane plots indicate stable/unstable steady states. Parameter values are α = d = 1. (A,B) Case 4, a = 3.3, b = 1.5, periodic solution around the unstable coexistence steady state E^*. (C,D) Case 5, a = 6, b = 1.05, bi-stability between the extinction steady state E₀ and the coexistence steady state E^*. (E,F) Case 6, a = 6, b = 1.4, only the extinction steady state E₀ is stable.

Discussion

In this article, we have studied the dynamics of a predator-prey model with Holling type III functional response and ratio dependence. Earlier numerical results for the same model with maturation delay suggested the possibility of extinction of both prey and predator populations that was, however, not investigated at that time [44]. Whereas, extinction steady state is a feasible steady state of the model for any values of parameters, due to ratio dependence in the functional response, the model could not be linearised in the neighborhood of extinction steady state in the standard manner. To overcome this difficulty, we have transformed the system using a Briot-Bouquet transformation, which resulted in the unfolding of the extinction steady state into three distinct steady states, whose stability could now be studied. Two of those steady states are unstable for any values of parameters, while the third is stable for sufficiently large predation rate to have an effect on prey who are reproducing logistically, and for sufficiently small values of the parameter α characterizing functional response. In both cases, we note that the region of stability of extinction steady state is shrinking for higher rates of predator fecundity, which can be explained by the fact that if for the same predation rate predators are better able to reproduce, this increases the possibility of a stable coexistence, where both populations are maintained at some steady level. Regions of feasibility of two steady states of the transformed model that correspond to extinction in the original model are determined explicitly in terms of system parameters, since they are given in the form of principal and n = −1 branches of the, in general, complex-valued Lambert W function.

By combining the results on stability of extinction steady state with conditions for stability of prey-only and coexistence equilibria, we are able to obtain a full picture of system dynamics in different parameter regions, identifying situations where only one of those steady states is stable, as well as cases of bi-stability, where extinction coexists with either a prey-only, or coexistence equilibria. In the parameter region, where the system exhibits sustained periodic oscillations around an unstable coexistence equilibrium, resulting from a supercritical Hopf bifurcation of this steady state, the extinction steady state is unstable, hence, no bi-stability is observed. While coexistence equilibrium is only biologically feasible, when the prey-only steady state is unstable, in contrast, the extinction steady state exists for all values of parameters and can coexist with either of those steady states. In the case of coexistence of steady states, for the same parameter values, depending on the initial conditions, the system will approach one of the two coexisting steady states. This observation is important from the perspective of analyzing real ecological data, as it provides clues about the possibility of different types of dynamics for the same or very similar values of parameters, some of which may be difficult to accurately measure. Earlier papers that focused largely on ratio dependence in models with Holling type II functional response [22, 24, 26] primarily showed either extinction of both species, or their stable coexistence. In contrast, our model with Holling type III response exhibits significant parameter regions, where prey-only equilibrium is stable either by itself, or in combination with extinction steady state. On the other hand, in the case when coexistence steady state is unstable, and there is a periodic solution around it, it is the only attractor in the model.

There are several directions, in which the work presented in this paper could be extended. In an earlier paper [44], we studied the effects of stochasticity and maturation delays for the same model, and some of the simulations there also indicated the possibility of observing extinction in the time-delayed model. Unfortunately, the standard Briot-Bouquet transformation, as used in this paper, would not work for the analysis of stability of extinction steady state in the time-delayed model, because even after the transformation, the singularity at the point (0, 0) would remain, and some alternative approach for studying stability of extinction equilibrium in a model with time delay needs to be developed. Another feature of many ecological models is the so-called stochastic extinction, where for sufficiently small population densities, species can go extinct, even though deterministically they could recover to some sustained levels at a later stage. It would be interesting and important to investigate how (in)stability of extinction steady state is affected by stochasticity, and how it impacts the regimes of bi-stability between extinction and other equilibria. This could be achieved through numerical exploration of basins of attraction of different steady states in the stochastic model in a manner similar to how it was done by Fatehi et al. [47, 48] in the context of modeling autoimmunity arising from immune response to a viral infection.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

KB and YK: conceptualization and analysis. KB, YK, and OB: numerical simulations and write-up. All authors contributed to the article and approved the submitted version.

Funding

OB has been supported by the Cancer Research UK (EDDCPJT\100022).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Freedman HI. Deterministic Mathematical Models in Population Ecology. New York, NY: Marcel Dekker, Inc. (1980).

Google Scholar

2. Gause GF. The Struggle for Existence. Baltimore, MD: Williams and Wilkins (1934).

Google Scholar

3. Kolmogorov AN. A qualitative study of mathematical models of population dynamics. (In Russian) Problemy Kibernetiki. (1972) 25:100.

Google Scholar

4. Solomon ME. The natural control of animal populations. J Anim Ecol. (1949) 18:1–35. doi: 10.2307/1578

CrossRef Full Text | Google Scholar

5. Lotka AJ. Elements of Physical Biology. Baltimore, MD: Williams and Wilkins (1925).

Google Scholar

6. Volterra V. Fluctuations in the abundance of a species considered mathematically. Nature. (1926) 188:558–60. doi: 10.1038/118558a0

CrossRef Full Text | Google Scholar

7. Holling CS. The components of predation as revealed by a study of small-mammal predation of the European Pine Sawfly. Can Entomol. (1959) 91:293–320. doi: 10.4039/Ent91293-5

CrossRef Full Text | Google Scholar

8. Holling CS. Some characteristics of simple types of predation and parasitism. Can Entomol. (1959) 91:385–98. doi: 10.4039/Ent91385-7

CrossRef Full Text | Google Scholar

9. Mondal B, Roy S, Ghosh U, Tiwari PK. A systematic study of autonomous and nonautonomous predator-prey models for the combined effects of fear, refuge, cooperation and harvesting. Eur Phys J Plus. (2022) 137:724. doi: 10.1140/epjp/s13360-022-02915-0

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Roy S, Tiwari PK, Nayak H, Martcheva M. Effects of fear, refuge and hunting cooperation in a seasonally forced eco-epidemic model with selective predation. Eur Phys J Plus. (2022) 137:528. doi: 10.1140/epjp/s13360-022-02751-2

CrossRef Full Text | Google Scholar

11. Roy S, Tiwari PK, Nayak H, Martcheva M. An ecoepidemic seasonally forced model for the combined effects of fear, additional foods and selective predation. J Biol Syst. (2022) 30:285–321. doi: 10.1142/S0218339022500103

CrossRef Full Text | Google Scholar

12. Hossain M, Pal S, Tiwari PK, Pal N. Bifurcations, chaos, and multistability in a nonautonomous predator-prey model with fear. Chaos. (2021) 31:123134. doi: 10.1063/5.0067046

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Tiwari PK, Verma M, Pal S, Kang Y, Misra AK. A delay nonautonomous predator-prey model for the effects of fear, refuge and hunting cooperation. J Biol Syst. (2021) 29:927–69. doi: 10.1142/S0218339021500236

CrossRef Full Text | Google Scholar

14. Akcakaya HR. Population cycles of mammals: evidence for a ratio-dependent predation hypothesis. Ecology. (1992) 62:119–42. doi: 10.2307/2937172

CrossRef Full Text | Google Scholar

15. Arditi R, Ginzburg LR, Akcakaya HR. Variation in plankton densities among lakes: a case for ratio-dependent predation models. Am Nat. (1991) 138:1287–96. doi: 10.1086/285286

CrossRef Full Text | Google Scholar

16. Gutierrez AP. Physiological basis of ratio-dependent predator-prey theory: the metabolic pool model as a paradigm. Ecology. (1992) 73:1552–63. doi: 10.2307/1940008

CrossRef Full Text | Google Scholar

17. Arditi R, Perrin N, Saïah H. Functional responses and heterogeneities: an experimental test with cladocerans. Oikos. (1991) 60:69–75. doi: 10.2307/3544994

CrossRef Full Text | Google Scholar

18. Arditi R, Ginzburg LR. Coupling in predator-prey dynamics: ratio-Dependence. J Theor Biol. (1989) 139:311–26. doi: 10.1016/S0022-5193(89)80211-5

CrossRef Full Text | Google Scholar

19. Arditi R, Ginzburg LR. How Species Interact: Altering the Standard View on Trophic Ecology. Oxford: Oxford University Press (2012).

Google Scholar

20. Beretta E, Kuang Y. Global analyses in some delayed ratio-dependent predator-prey systems. Nonlinear Anal. (1998) 32:381–408. doi: 10.1016/S0362-546X(97)00491-4

CrossRef Full Text | Google Scholar

21. Jost C, Arino O, Arditi R. About deterministic extinction in ratio-dependent predator-prey models. Bull Math Biol. (1999) 61:19–32. doi: 10.1006/bulm.1998.0072

CrossRef Full Text | Google Scholar

22. Xiao D, Ruan S. Global dynamics of a ratio-dependent predator-prey system. J Math Biol. (2001) 43:268–90. doi: 10.1007/s002850100097

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Haque M. Ratio-dependent predator-prey models of interacting populations. Bull Math Biol. (2009) 71:430–52. doi: 10.1007/s11538-008-9368-4

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Arancibia-Ibarra C, Aguirre P, Flores J, Heijster P. Bifurcation analysis of a predator-prey model with predator intraspecific interactions and ratio-dependent functional response. Appl Math Comput. (2021) 402:126152. doi: 10.1016/j.amc.2021.126152

CrossRef Full Text | Google Scholar

25. Aguirre P, Flores JD, González-Olivares E. Bifurcations and global dynamics in a predator-prey model with a strong Allee effect on the prey, and a ratio-dependent functional response. Nonlinear Anal. (2014) 16:235–249. doi: 10.1016/j.nonrwa.2013.10.002

CrossRef Full Text | Google Scholar

26. Pal PJ, Saha T. Qualitative analysis of a predator-prey system with double Allee effect in prey. Chaos Solitons Fractals. (2015) 73:36–63. doi: 10.1016/j.chaos.2014.12.007

CrossRef Full Text | Google Scholar

27. Tyutyunov YV, Titova LI. From Lotka-Volterra to Arditi-Ginzburg: 90 years of evolving trophic functions. Biol Bull Rev. (2020) 10:167–85. doi: 10.1134/S207908642003007X

CrossRef Full Text | Google Scholar

28. Dawes JHP, Souza MO. A derivation of Holling's type I, II and III functional responses in predator-prey systems. J Theor Biol. (2013) 327:11–22. doi: 10.1016/j.jtbi.2013.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Fan YH, Li WT. Permanence for a delayed discrete ratio-dependent predator-prey system with Holling type functional response. J Math Anal Appl. (2004) 299:357–74. doi: 10.1016/j.jmaa.2004.02.061

CrossRef Full Text | Google Scholar

30. Pal PJ, Mandal PK, Lahiri KK. A delayed ratio-dependent predator-prey model of interacting populations with Holling type III functional response. Nonlinear Dyn. (2014) 76:201–20. doi: 10.1007/s11071-013-1121-3

CrossRef Full Text | Google Scholar

31. Wang X, Peng M, Lu X. Stability and Hopf bifurcation analysis of a ratio-dependent predator-prey model with two time delays and Holling type III functional response. Appl Math Comput. (2015) 268:496–508. doi: 10.1016/j.amc.2015.06.108

CrossRef Full Text | Google Scholar

32. Li WT, Wu SL. Traveling waves in a diffusive predator-prey model with Holling type-III functional response. Chaos Solitons Fractals. (2008) 37:476–86. doi: 10.1016/j.chaos.2006.09.039

CrossRef Full Text | Google Scholar

33. Apreutesei N, Dimitriu G. On a prey-predator reaction-diffusion system with Holling type III functional response. J Comput Appl Math. (2010) 235:366–79. doi: 10.1016/j.cam.2010.05.040

CrossRef Full Text | Google Scholar

34. Guin LN, Mandal PK. Spatial pattern in a diffusive predator-prey model with sigmoid ratio-dependent functional response. Int J Biomath. (2014) 7:1450047. doi: 10.1142/S1793524514500478

CrossRef Full Text | Google Scholar

35. Morozov A. Emergence of Holling type III zooplankton functional response: bringing together field evidence and mathematical modelling. J Theor Biol. (2010) 265:45–54. doi: 10.1016/j.jtbi.2010.04.016

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Schenk D, Bersies LF, Bacher S. An experimental test of the nature of predation: neither prey- nor ratio-dependent. J. Anim. Ecol. (2005) 74:86–91. doi: 10.1111/j.1365-2656.2004.00900.x

CrossRef Full Text | Google Scholar

37. Sarnelle O, Wilson AE. Type III functional response in D aphnia. Ecology. (2008) 89:1723–32. doi: 10.1890/07-0935.1

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Kratina P, Vos M, Bateman A, Anholt BR. Functional responses modified by predator density. Oecologia. (2009) 159:425–33. doi: 10.1007/s00442-008-1225-5

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Holt J, Jeger MJ, Thresh JM, Otim-Nape GW. An epidemiological model incorporating vector population dynamics applied to African cassava mosaic virus disease. J Appl Ecol. (1997) 34:793–806. doi: 10.2307/2404924

CrossRef Full Text | Google Scholar

40. Blyuss KB, Basir FA, Tsygankova VA, Biliavska LO, Iutynska GO, et al. Control of mosaic disease using microbial biostimulants: insights from mathematical modelling. Ricerche di Matematica. (2020) 69:437–55. doi: 10.1007/s11587-020-00508-6

CrossRef Full Text | Google Scholar

41. Ivlev VS. Experimental Ecology of the Feeding of Fishes. New Haven, CT: Yale University Press (1947).

Google Scholar

42. Ricker WE. Stock and recruitment. J Fisher Res Board Can. (1954) 11:559–623. doi: 10.1139/f54-039

CrossRef Full Text | Google Scholar

43. Basir FA, Kyrychko YN, Blyuss KB, Ray S. Effects of vector maturation time on the dynamics of Cassava Mosaic Disease. Bull Math Biol. (2021) 83:87. doi: 10.1007/s11538-021-00921-4

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Blyuss KB, Kyrychko SN, Kyrychko YN. Time-delayed and stochastic effects in a predator-prey model with ratio dependence and Holling type III functional response. Chaos. (2021) 31:073141. doi: 10.1063/5.0055623

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Zhang T, Ding T, Huang W, Dong Z. Qualitative Theory of Differential Equations. Providence: American Mathematical Society (1992).

Google Scholar

46. Dumortier F, Llibre J, Artés J. Qualitative Theory of Planar Differential Systems. Berlin: Springer-Verlag (2006).

Google Scholar

47. Fatehi F, Kyrychko SN, Ross A, Kyrychko YN, Blyuss KB. Stochastic effects in autoimmune dynamics. Front Physiol. (2018) 9:45. doi: 10.3389/fphys.2018.00045

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Fatehi F, Kyrychko YN, Blyuss KB. Stochastic dynamics in a time-delayed model for autoimmunity. Math Biosci. (2020) 322:108323. doi: 10.1016/j.mbs.2020.108323

PubMed Abstract | CrossRef Full Text | Google Scholar