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

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

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)=a1zn/(1+a2zn) 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)=a1zn/(1+a2zn) 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,

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

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.

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

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

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

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

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

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

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

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

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

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₂.

Linearisation of the system (Equation 5) near any steady state E^=(u^,z^) gives the Jacobian

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

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.

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.

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.

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.

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.

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

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.