Modelling and analysis of a modified May-Holling-Tanner predator-prey model with Allee effect in the prey and an alternative food source for the predator

1.   Introduction

Numerous investigations have extensively studied the interactions between the specialist predator least-weasel (Mustela nivalis) and the field-vole (Microtus agrestis) ^[1,2,3]. These investigations have revealed important findings in the variation of the populations in Fennoscandia from regular cycles to stable populations. This variation is also affected by the fact that the rodent breeding is seasonal, and the duration of the breeding varies with the physical location in Fennoscandia ^[4]. Moreover, cycles may be affected by the existence of an optimum group size and an Allee effect ^[5]. Graham and Lambin ^[6] showed that field-vole survival depends on the variation of weasel predation. They also demonstrated that the proportion of weasels was suppressed in summer and autumn, while the vole population always declined to a low density. Thus, if the predation by weasels is reduced then the field-vole survival increases. In ^[6], the authors assumed that the May–Holling–Tanner model was not appropriate to conduct such studies due to the large number of parameters. Additionally, the authors suggested that there are different variations in these two species that should be considered in order to improve study the results.

The most common mathematical framework for describing this predator-prey interaction is that the weasel and field-vole have a logistic-type growth function and the predator carrying capacity is prey dependent ^[7]. The classic model used to represent the predation between these two species is the May–Holling–Tanner predator-prey model ^[8] which is a modification of original the Leslie–Gower model. Leslie (1948) proposed the following generalisation of the Lotka-Volterra model

where $N(t)$ is used to represent the size of the prey population at time $t$, $P(t)$ is used to represent the size of the predator population at time $t$, $r$ is the intrinsic growth rate for the prey, $s$ is the intrinsic growth rate for the predator, $K$ is the prey environmental carrying capacity, $n$ is a measure of the quality of the prey as food for the predator, $nN$ can be interpreted as a prey dependent carrying capacity for the predator and the term $P/nN$ is known as the Leslie-Gower term. It measures the loss in the predator population due to rarity (per capita $P/N$) of its favourite food.

This model retains the prey equation from the Volterra model, but for the predator it assummes a logistic-like term, in which the predator carrying capacity is directly proportional to prey density. This interesting formulation for the predator dynamics has been discussed in Leslie and Gower in ^[9] and Pielou in ^[10]. Later, May (1974) in ^[8] modified the Leslie-Gower model (1.1) by assuming the hyperbolic functional response (Holling type Ⅱ) in the prey equation, $H(N) = qN/(N+a)$, which occurs in species when the number of prey consumed rises rapidly at the same time that the prey density increases ^[7,8,11]. The parameter $q$ is the maximum predation rate per capita and $a$ is half of the saturated level ^[7]. In particular, it is given by

Model (1.2) is known as Holling–Tanner or as May–Holling–Tanner model which is described by an autonomous two-dimensional system of ordinary differential equations ^[7,12,13].

The original May–Holling–Tanner model (1.2) was studied, among others, by ^[14]. The authors proved that all the solutions which are initiated in $\mathbb{R}^2_+$ are bounded and the equilibrium point $(1, 0)$, which is the scaled equilibrium point $(K, 0)$, is a saddle point. Moreover, by doing a polar blowing-up ^[15] at the origin, Saez and González-Olivares ^[14] showed that there exists a separatrix curve between a hyperbolic and a parabolic sector in the neighbourhood of the origin. The authors also found that there is always one positive equilibrium point in the first quadrant which can be stable; unstable surrounded by a stable limit cycle; or stable surrounded by two limit cycles. Additionally, the global stability of a periodic solution of this predator-prey model was investigated in ^[16] and the analysis of a discrete May–Holling–Tanner model was described in ^[17].

Additional complexity can be incorporated into the predator-prey model (1.2) in order to make it more realistic. We can consider, in the prey population, a mechanism connected with a density-dependent phenomenon (Allee effect ^{[18,19,20,21,22,23]}). The Allee effect may appear due to a wide range of biological phenomena such as reduced anti-predator vigilance; social thermo-regulation; genetic drift; mating difficulty; reduced defence against the predator; and deficient feeding because of low population densities ^[24]. The Allee effect may appear in the field-vole population due to the density-dependent habitat selection showed in ^[25,26]. Moreover, the stabilisation of field-vole populations depends on the variation in reproductive rate and recruitment of the population ^[27]. The influence of the Allee effect upon the logistic-type growth in the prey equation is represented by the inclusion of a multiplier in the form of $N-m$ (1.3). For $0 < m < K$, the per-capita growth rate of the prey population with the Allee effect included is negative, but increasing, for $N\in[0, m)$, and this is referred to as the “strong Allee effect”. When $m\leq0$^*, the per-capita growth rate is positive but increases at low prey population densities and this is referred to as the weak Allee effect ^[20,28]. The per capita growth rate for the logistic growth function and Allee effects (strong and weak) are shown in Figure 1. We observe that the per capita growth rate for a weak Allee effect is positive and growing, when compared to the strong Allee effect. Moreover, the trajectories of a population affected by a weak Allee effect are delayed compared with a strong Allee effect. This effect can be generated by the reduction of the male population ^[29] or the habitat selection ^[25,26].

^*Note that $m = 0$ is often also called the “weak Allee effect”, however, the behaviour of the model with $m = 0$ is similar to the case of a strong Allee effect ($m > 0$) and thus we will not consider the case of $m = 0$ in this manuscript.

Some species are considered as a generalist and specialist predator and thus it can switch to another available food, although its population growth may still be limited by the fact that its preferred food is not available in abundance. For instance, McDonald et al. ^[30] reveal that weasels' diet involved mainly field-vole, which correspond to $68\%$ of the diet which can be considered as the primary food source. Moreover, the authors also reveal that lagomorphs with $25\%$ of the diet; birds and birds' eggs with $5\%$ of the diet and other types of food with $2\%$ of the diet can be considered as a secondary food source, see Figure 2. This characteristic was not considered in investigations associated with the interaction between the field-vole as a prey and least-weasel as a predator ^{[1,2,7,31,32,33,34,35]}. That is, these studies do not consider that, since the predator is a generalist, it can survive under different environments and utilise a large range of food resources. Instead of adding more species to the model, we assume that these other food sources are abundantly available ^[30] and model this characteristic by adding a positive constant $c$ to the environmental carrying capacity for the predator ^[36]. Therefore, we have made a modification to the prey-dependent logistic growth term in the predator equation, namely $nN$ is replaced by $nN+c$ (1.5). As a result, we observe that the modified May–Holling–Tanner model with an alternative food source (1.6) is no longer singular for $N = 0$.

The aim of this manuscript is to extends the properties of the weasel and field-vole interaction studied in ^[6] by showing the impact of an Allee effect on the prey (1.4) and an alternative food source for the predator (1.6) in the original May–Holling–Tanner model. In addition, the analysis of the modified model with weak Allee effect (i.e., systems (1.4) with $m < 0$) and the modified model considering an alternative food for the predator (i.e., system (1.6)) show a strong connection with the findings presented in ^[3,33]. Study results suggest that the system does not lose the coexistence and oscillation of the population when a constant alternative food source for weasels and a weak Allee effect on the prey are included. Moreover, the survival of the field-vole does not depend on the level of predation when a weak Allee on the prey is included ^[6]. However, we have shown that the May–Holling–Tanner model with strong Allee effect on the prey (i.e., systems (1.4) with $m > 0$) reduce the complexity of the dynamics in the modified model. The modified model with strong Allee effect supports, depending on the initial conditions, extinction and coexistence of both the prey and the predator species.

The mathematical interaction in systems (1.2), considering the modifications separately and together, is study further in Section 2. In Section 3, we conclude the manuscript, presenting the impact of the modification on the dynamic of the modified May–Holling–Tanner models. We also discuss the ecological implications of the inclusion of the strong and weak Allee effects on the prey and of an alternative food for the predator in a real model studied in ^[3,6,33].

2.   Models

The results derived in ^[14] are for the nondimensionalised versions of the modified versions of systems (1.2) which is given by

However, this nondimensionalised version is topologically equivalent to the original version of the model as we used a diffeomorphism preserving the orientation of time. Therefore, the results can be transferred to the original May–Holling–Tanner model. Moreover, we can extend the analysis by using the system parameters presented in ^[3]. In this paper, Turchin and Hanski used the original May–Holling–Tanner model and the authors showed that the oscillation of the field-vole is driven by their interaction with the weasel. This oscillation is obtained by using the system parameters presented in Table 2.

Nowadays, there are several computational methods to find different type of bifurcations. These methods are implemented in software packages such as MATCONT ^[37]. Figure 3 illustrates the Hopf bifurcation which was detected with MATCONT for the parameters $r$, $K$, $a$, $n$ and $c$ fixed, while we will vary $q$ and $s$, see Table 2.

2.1.   Original model

The original May–Holling–Tanner predator-prey model with functional response Holling type Ⅱ (1.2) was studied mathematically, for instance, in ^[14,38] and used to model real predator-prey interactions ^[1,2,7,39]. The bifurcation diagram of the original May–Holling–Tanner model (1.2) created with MATCONT ^[37] and two typical phase portraits are shown in the first row of Figure 3. The Hopf curve^† for the May–Holling–Tanner model shows that if the predation rate $q$ and the intrinsic growth of the predator $s$ are located in the hatched green area, then system (1.2) supports the stabilisation of the predator and the prey populations. If the predation rate is located in the blue area, then there are conditions for which the system supports the oscillation of both populations, i.e., there is a stable limit cycle surrounding a positive equilibrium point in the first quadrant. Moreover, there are conditions in the $(q, s)$ parameter space (near the Hopf curve) for which populations can coexist and also oscillate simultaneously, i.e., there are two limit cycles surrounding a positive equilibrium point in the first quadrant. However, for the sake of the presentation, this region is not included in the bifurcation diagram. Instead, we refer to ^[14]. Note that the extinction of (one of) the species is not possible in this model.

^†The Hopf curve is where the stability of an equilibrium point switches and, in some cases, a periodic solution arises.

2.2.   Allee effect (strong and weak)

The May–Holling–Tanner model with a strong Allee effect ($m > 0$) was studied in ^[40]. The authors proved that system (1.4) with $m > 0$ has three equilibrium points in the $x$-axis and at most two positive equilibrium points, see the left panel of Figure 6. The authors also proved that one of the positive equilibrium points is always a saddle. While, the other can be stable, unstable or stable surrounded by an unstable limit cycle. For $m < 0$ there are two equilibrium points on the $x$-axis and at most three positive equilibrium points, see the right panel of Figure 6. The stability of the equilibrium points of system (1.4) considering $m < 0$ is studied partially in Appendix A.

First, we recall that as we used a diffeomorphism preserving the orientation of time the nondimensionalised version (A.2) is topologically equivalent to the original version of the model (1.4). Therefore, the results can be transferred to the modified May–Holling–Tanner model. Moreover, we can extend the analysis by using the system parameters $(r, K, a, n, c, m) = (4,150, 6, 0.025, 0.01, \pm15)$ (see Table 2). The bifurcation diagram of the May–Holling–Tanner model with weak Allee effect (1.4) ($m < 0$) and two typical phase portraits are shown in the third row of Figure 3. System (1.4) with $m < 0$ is not discussed in detail in this manuscript, but the bifurcation diagram can be obtained by using the techniques presented in ^[40]. If the parameters $(q, s)$ are located in the hatched green region, then system (1.4) has one positive equilibrium point which is stable, and thus both the predator and the prey population can coexist. If the parameters $(q, s)$ are located in the blue region, then the equilibrium point is unstable surrounded by a stable limit cycle and thus both the predator and the prey populations oscillate.

However, if the parameters $(q, s)$ are located in the grey region, then system (1.4) has three equilibrium points. These equilibrium points can be a saddle, a stable and/or an unstable node. Therefore, system (1.4) supports the coexistence and/or oscillation of the predator and prey populations. Note that the dynamics of system (1.4) with $m < 0$, outside the grey region, are qualitatively similar to the dynamics of system (1.2) as shown in row one in Figure 3 since there is a large region in the $(q, s)$ parameter space in which there is always one positive equilibrium point, i.e., the hatched green region.

Additionally, the Allee effect, both strong and weak, decreases the range of the value of the predation rate for which both populations can coexist, i.e., $\tilde{q}_2 < q_1$ (see Figure 3). Moreover, as in the original May–Holling–Tanner model, the predator and prey populations never become extinct in system (1.4) and thus there are always conditions in the system parameters where both populations coexist (see hatched green area in Figure 3) and/or oscillate (see blue and grey areas in Figure 3).

The bifurcation diagram of the modified May–Holling–Tanner model with strong Allee effect (1.4) ($m > 0$) studied in ^[40] and two typical phase portraits are shown in the fourth row of Figure 3. The Hopf curve for the model shows that if $(q, s)$ are located in the green area, system (1.4) supports the stabilisation of both populations and also the extinction of both populations. However, if the predation rate is located in the red hatched area in Figure 3 ($q > \tilde{q}_1$), then system (1.4) does not have positive equilibrium points in the first quadrant and thus the system only supports the extinction of both the predator and the prey populations. In contrast, if the predation rate and the intrinsic growth of the predator $(q, s)$ are located in the solid red area, then system (1.4) has two positive equilibrium points which are a saddle and an unstable node, and thus again we observe the extinction of the predator and prey populations.

2.3.   Alternative food

The model with alternative food for predators corresponds to system (1.6). This model was studied in ^[13,41,42] where the authors presented the stability of the equilibrium points and conditions for periodic solutions. The bifurcation diagram of the modified May–Holling–Tanner model (1.6) and two typical phase portraits are shown in the second row of Figure 3, see also ^[42,44]. The Hopf curve for the modified May–Holling–Tanner model shows that if the predation rate $q$ and the intrinsic growth of the predator $s$ are located in the green area, then system (1.6) supports the stabilisation of both populations and also the extinction of only the prey population. If the predation rate is located in the hatched brown area in Figure 3 ($q > q_2$), then system (1.6) does not have a positive equilibrium point in the first quadrant and thus the extinction of the prey population is observed, while the predator population stabilises at $(0, c)$ (the equilibrium point $(0, c)$ is related to the alternative food source for the predator).

Moreover, when the parameters $(q, s)$ are located in the solid brown area, then system (1.6) has two positive equilibrium points in the first quadrant, namely a saddle and an unstable node, and the extinction of the prey population is again observed for all initial conditions. When the parameters $(q, s)$ are located in the hatched brown area, then system (1.6) does not have equilibrium points in the first quadrant and the extinction of the prey population is observed for all initial conditions. Additionally, in ^[42,44] the authors have shown that there is a region, near to the Hopf bifurcation, in which both populations oscillate and coexist at the same time, i.e., there are two limit cycles surrounding a positive equilibrium point in the first quadrant. For the sake of the presentation, we did not include this region in the bifurcation diagram.

2.4.   Allee effect (strong and weak) and alternative food

The May–Holling–Tanner model with Allee effect on the prey and alternative food for predators is given by system (1.7). This model was studied partially in ^[43] and in Appendix B. The bifurcation diagram of the modified May–Holling–Tanner model with Allee effect on the prey and alternative food for the predator (1.7) and two typical phase portraits are shown in the last row of Figure 3. System (1.2) with Allee effect is not discussed in detail in this manuscript, but the bifurcation diagram can also be obtained by using the techniques presented in ^[40]. The Hopf curve for the system (1.7) with Allee effect and alternative food shows that system (1.7) supports the coexistence of both populations, the extinction of the prey population only (green area) and the extinction of the prey population (solid and hatched brown regions). That is, the behaviour is similar to the modified May–Holling–Tanner model (1.6), see the second row of Figure 3.

In Figure 3, we show the bifurcation diagrams (first columns) in which different regions represent different behaviours. In the solid green region, sets of initial conditions exist that lead to the coexistence of both species, while there are also open sets of initial conditions that lead to the extinction of both species. In the hatched green regions, all initial conditions lead to the coexistence of both populations. In the solid blue regions, all initial conditions lead to the oscillation of both populations. In the solid and hatched red regions, all initial conditions lead to the extinction of both species, while in the solid and hatched brown regions, all initial conditions lead to the extinction of the prey only. In the right column, the orange region (dark and light) represents the basin of attraction of the stabilisation of both populations, the light blue (light grey) region represents the basin of attraction of the extinction of both (only the prey) populations and the yellow region represent the oscillation of both populations.

Figure 4 shows the time series of the predator and prey populations in the original May–Holling–Tanner model (1.7), the May–Holling–Tanner model with alternative food for the predator (1.6) and the May–Holling–Tanner model with Allee effect on the prey (1.4). The system parameters are as in Table 2 and $(q, s) = (700, 1.25)$. Note that these values lie inside the solid or hatched green regions in Figure 3. Moreover, different predator and prey initial densities are taken into account to highlight the different types of behaviours. If the parameters $(q, s)$ are located in the blue region of Figure 3, then the original model supports oscillation of the solutions and it does not depend on the initial conditions, see the first column of Figure 4. That is, every initial density evolves to the coexistence of the species. If the parameters $(q, s)$ are located in the green region of Figure 3, then the May–Holling–Tanner model with strong Allee effect ($m > 0$) supports both coexistence and extinction of both species. That is, there exists an initial density which evolves to the coexistence of the species, while there are also initial conditions which evolve to the extinction of both species (see second column of Figure 4).

If the parameters $(q, s)$ are located in the grey region of Figure 3, then the May–Holling–Tanner model with weak Allee effect ($m < 0$) supports the stabilisation of both populations to two different equilibrium points. That is, there exists an initial density which evolves to a coexistence point, while there are also initial conditions which evolve to a lower coexistence point (see third column of Figure 1). Finally, the modified May–Holling–Tanner model with an alternative food source for the predator supports, depending on the initial conditions, both oscillation and the extinction of the prey species only. That is, there exists as initial density which evolves to a limit cycle which represents the oscillation of the species, while there are also initial conditions which evolve to the extinction of the prey and the stabilisation of the predator population (see fourth column of Figure 1).

3.   Conclusion

In this manuscript, the May–Holling–Tanner predator-prey model with functional response Holling type Ⅱ (1.2) was reviewed. We also considered that the prey population is affected by strong and weak Allee effect (1.3) and the predator population has an alternative food source (1.5). Using a diffeomorphism we extend the analysis of a topologically equivalent system (1.4), (1.6), (1.7) and (1.2).

We have also extended the analysis by using the system parameters presented in ^[2,7], see Table 2. In addition, we have shown that the predation rate per capita (parameter $q$) impacts the number of equilibrium points in the first quadrant while the intrinsic growth rate of the predator (parameter $s$) changes the behaviour of the positive(s) equilibrium point(s), see Figure 3. Therefore, it is clear that the self-regulation depends on the values of these two parameters. Consequently, we have shown that the original May–Holling–Tanner model with a strong Allee effect ($m > 0$), a weak Allee effect ($m < 0$) and with both modifications, i.e., Allee effects (strong and weak) for the prey and alternative food for the predator, can mean extinction and coexistence of prey and predator populations, depending on particular conditions. This is because the stable manifold of a saddle equilibrium point or an unstable limit cycle act as a separatrix between the basin of attraction of the origin and a stable equilibrium point in the first quadrant. However, the original May–Holling–Tanner model with a weak Allee effect ($m < 0$) and with alternative food for the predator are present with the coexistence or oscillation of the population. This is due to the existence of at least one positive equilibrium point which can be either stable, unstable surrounded by a stable limit cycle or stable surrounded by two limit cycles.

In Figure 5 we observe that the introduction of a density dependent phenomenon in the prey population and/or an alternative food source for the predator population have an impact on the predation rate (parameter $q$). The predation rate at which the Hopf curve of the Holling–Tanner model with Allee effect (1.4) is reduced, see the top right panel of Figure 5. In particular, the Hopf curve of the model with strong Allee effect shows that system (1.4) can have up to two positive equilibrium points which can collapse at $q = \tilde{q}_1$ (respectively for $m < 0$ at $q = \tilde{q}_2$). Therefore, for $q < \tilde{q}_1$ system (1.4) presents two positive equilibrium points and thus the coexistence of both populations is observed ^[40].

In contrast, the Hopf curve of the model with week Allee effect ($m < 0$) shows that if $q > \tilde{q}_2$ then system (1.4) with $m < 0$ has always one positive equilibrium point which is stable or unstable surrounded by a stable limit cycle and thus the coexistence or oscillation of both population is observed. Note that the dynamic of system (1.4) with $m < 0$ is similar with the dynamics of system (1.2) since in both systems there is always one positive equilibrium point. Consequently, there is a coexistence and/or oscillation of predator and prey population.

We also consider the Hopf curve of the Holling–Tanner model with alternative food for predators (1.6) and the system parameters $(r, K, a, n, c, m) = (4,150, 6, 0.025, 0.01)$ (see Table 2) fixed, see the bottom left panel of Figure 5. In particular, the Hopf curve for the original May–Holling–Tanner model, considering an alternative food for predators, shows that the positive equilibrium points collapse at $q = q_2$. Therefore, for $q < q_2$ system (1.6) present also two positive equilibrium points. One of them is always a saddle point and the other can be stable or unstable surrounded by stable limit cycle and thus the coexistence or oscillation of both population is also observed ^[42,44]. Note that if $q > q_2$ the equilibrium point related to the alternative food is globally stable. Therefore, system (1.6) supports the stabilisation of only the predator population and the extinction of the prey population.

Next, we consider the Hopf curve of the May–Holling–Tanner model with both modification (1.7), i.e., a strong Allee effect on the prey and alternative food for the predator, and the system parameters $(r, K, a, n, c, m) = (4,150, 6, 0.025, 0.01, \pm15)$ (see Table 2) fixed, see the bottom right panel of Figure 5. In particular, the modified system (1.7) can have up to two positive equilibrium points which can collapse at $q = \hat{q}_1$. Therefore, for $q < \hat{q}_1$ system (1.7) support the coexistence of both population or the extinction of only the prey population.

In summary, the bifurcation diagrams of the models studied in ^[14,40,43] are often qualitatively similar but their solutions behave quantitatively differently. In other words, it is observed that there are two groups of modes which support equivalent ecological behaviour due to the addition of the modifications in the May–Holling–Tanner model. That is, an alternative food source for the predator, a strong Allee effect or the combination of both modifications support the coexistence and the extinction of at least one of the species. In contrast, the original model and the model with weak Allee effect ($m < 0$) do not support the extinction of one of the species.

Acknowledgements

This project was fully supported by Universidad de las Américas (UDLA) Chile. Especially, authors would like to give special thanks to Jaime Garcia for his support in the edition of this manuscript.

Conflict of interest

The authors declare that no conflict of interest.

Appendices

A.   The modified May–Holling–Tanner model with weak Allee effect (A.2) ($M < 0$)

The May–Holling–Tanner model with a strong Allee effect ($m > 0$) was studied in ^[40] in which the authors simplified the system by introducing a change of variable and time rescaling. This is given by the function $\phi :\breve{\Omega}\times\mathbb{R}\rightarrow \Omega\times\mathbb{R}$, where $(N, P, t) = \phi(u, v, \tau) = (Ku, nKv, \tau u(u+a/K)/rK)$. Moreover, by setting

into system (1.4) we obtain

To analyse the positive solution(s) of system (A.2) we consider three cases of the intersection of the nullcline of system (A.2) and $M < 0$. The $u$-nullclines of system (A.2) are $u = 0$ and $v = (u+A)(1-u)(u-M)/Q$, while the $v$-nullclines are $v = 0$ and $v = u$. Hence, the equilibrium points for system (A.2) are $(0, 0)$, $(1, 0)$ and the point(s) $(u^*, v^*)$ with $v^* = u^*$ and where $u^*$ is determined by the solution(s) of

where $H(A, M) = M+1-A$ and $L(A, M, Q) = A(M+1)-Q-M$. System (A.2) always has one positive equilibrium point, which we denote by $P_1 = (u_1, u_1)$ where $0 < u_1 < 1$. Factoring out $(u-u_1)$ from the equation resultant (A.3) of the intersection of the nullcline of system (A.2) ^[40], leaves as with second order polynomial equation

Therefore, if $1-A+M-u_1 > 0$ then the roots of (A.4) are given by

such that $0 < u_2 < u_3 < 1$.

A.1. Stability of positive equilibrium points:

To determinate the stability of the positive equilibrium points $P_1$, $P_2$, and $P_3$ of system (A.2) with $M < 0$ we compute the Jacobian matrix (A.6) at these equilibrium points

with $i = 1, 2, 3$. Thus, the trace and the determinant of the Jacobian matrix (A.6) are

with

Then, we consider the stability of these positive equilibrium points $P_1$, $P_2$ and $P_3$ of system (A.2) with $M < 0$ in the interior of $\Phi = \{\left(u, v\right), \ 0 < u\leq1, \ v\geq0\}$.

Lemma A.1. Let the system parameters of (A.2) be such that $1-A+M\leq0$ or $1-A+M > 0$ and $A(M+1)-Q-M\geq0$. Then system (A.2) has only one positive equilibrium point $(u_1, u_1)$ which can be

(i) a repeller if $S < f(u_1)$,

(ii) an attractor if $S > f(u_1)$,

with $f(u_*)$ defined in (A.8).

Proof. Evaluating $u^2(2u-(1-A+M))-AM$ at $u = u_1$ gives:

Then,

(ⅰ) if $1-A+M < 0$ then $u_1^2(2u_1-(1-A+M))-AM > 0$. Hence $\det(J(P_1)) > 0$ since $M < 0$,

(ⅱ) if $1-A+M = 0$ then $2u_1^3-AM > 0$. Hence $\det(J(P_1)) > 0$,

(ⅲ) if $1-A+M > 0$ then $u_1^2(2u_1-(1-A+M))-AM = u_1^3+(A(1+M)-Q-M)u_1-2AM > 0$. Hence $\det(J(P_1)) > 0$ since $A(M+1)-Q-M\geq0$.

Therefore, the sign of the trace (A.7), and thus the behaviour of $(u_1, u_1)$ depends on the parity of $f(u_1)-S$.

Lemma A.2. Let the system parameters of (A.2) be such that $1-A+M > 0$, $A(M+1)-Q-M < 0$ and $\Delta > 0$ (A.5). Then the equilibrium point $P_1$ is

(i) a saddle point if $h(u_1) < 0$

(ii) a repeller if $h(u_1) > 0$ and $S < f(u_1)$,

(iii) an attractor if $h(u_1) > 0$ and $S > f(u_1)$.

with $f(u_*)$ and $h(u_*)$ defined in (A.8).

Proof. Evaluating $u^2(2u-(1-A+M))-AM$ at $u = u_1$ gives:

Then, the sign of the determinant and the trace (A.7), and thus the behaviour of $P_1$ depends on the parity of $h(u_1)$ and $f(u_1)-S$.

Lemma A.3. Let the system parameters of (A.2) be such that $1-A+M > 0$, $A(M+1)-Q-M < 0$ and $\Delta > 0$ (A.5). Then the equilibrium point $P_2$ is

(i) a saddle point if $u_2 > u_1$

(ii) a repeller if $u_2 < u_1$ and $S < f(u_2)$,

(iii) an attractor if $u_2 < u_1$ and $S > f(u_2)$,

with $f(u_*)$ defined in (A.8).

Proof. Evaluating $u^2(2u-(1-A+M))-AM$ at $u = u_2 = (1-A+M-u_1-\sqrt{\Delta})/2$ gives:

Then, the sign of the determinant depends on the parity of $u_1-u_2$. Then, the trace (A.7), and thus the behaviour of $P_2$ depends on the parity of $f(u_2)-S$.

Lemma A.4. Let the system parameters of (A.2) be such that $1-A+M > 0$, $A(M+1)-Q-M < 0$ and $\Delta > 0$ (A.5). Then the equilibrium point $P_3$ is

(i) a saddle point if $u_3 < u_1$

(ii) a repeller if $u_3 > u_1$ and $S < f(u_3)$,

(iii) an attractor if $u_3 > u_1$ and $S > f(u_3)$,

with $f(u_*)$ defined in (A.8).

Proof. Evaluating $u^2(2u-(1-A+M))-AM$ at $u = u_3 = (1-A+M-u_1+\sqrt{\Delta})/2$ gives:

Then, the sign of the determinant depends on the parity of $u_3-u_1$. Then, the trace (A.7), and thus the behaviour of $P_3$ depends on the parity of $f(u_3)-S$.

B.   The modified May–Holling–Tanner model with alternative food for the predator and Allee effect

The May–Holling–Tanner model with Allee effect on the prey and alternative food for predators is given by system (1.7). The system can be simplified by introducing a change of variable and time rescaling presented in ^[43]. It is given by the function $\zeta :\breve{\Psi}\times\mathbb{R}\rightarrow \Psi\times\mathbb{R}$, where $\zeta(u, v, \tau) = (N, P, t) = (Ku, nKv, \tau(u+a/K)(u+c/(nK))/rK)$ ^[43]. By setting the parameters $A$, $C$, $Q$, $S$ and $M$ defined in (A.1), into system (1.7) we obtain

To analyse the positive solution(s) of system (B.1) we consider three cases of the intersection of the nullcline of system (B.1). The $u$-nullclines of system (B.1) are $u = 0$ and $v = (u+A)(1-u)(u-M)/Q$, while the $v$-nullclines are $v = 0$ and $v = u+C$. Hence, the equilibrium points for system (B.1) are $(0, 0)$, $(1, 0)$, $(M, 0)$ and the point(s) $(u^*, v^*)$ with $v^* = u^*+C$ and where $u^*$ is determined by the solution(s) of

Define the functions $g(u) = (u+A)(1-u)(u-M)$ and $h(u) = Q(u+C)$ and observe that $\lim\limits_{u \rightarrow \pm \infty} g(u) = \mp \infty$ and $g(0) = -AM < 0$. So, (B.2) will always have a single negative real root, which we denote by $u = -G$ where $G > 0$. Factoring out $(u+G)$ from (B.2) leaves us with second order polynomial

and thus we get $G(M+Q-A(M+1)+G(G-A+M+1))-AM-CQ = 0$. Note that since $Q > 0$, then we obtain that $A < G < C$ or else, $A > G > C$. In addition, we get that $Q = (G+1)(G+M)(G-A)/(C-G)$. The solution of (B.3) are:

such that $M < u_1\leq E \leq u_2 < 1$, where $E = (1-A+M+G)/2$. In particular,

1. if $\Delta < 0$ (B.4), then (B.1) has no equilibrium points in the first quadrant;

2. if $\Delta > 0$ (B.4), then (B.1) has two equilibrium points $P_{1, 2} = (u_{1, 2}, u_{1, 2}+C)$ in the first quadrant; and

3. if $\Delta = 0$ (B.4), then (B.1) has one equilibrium point $(E, E+C)$ of order two in the first quadrant,

Additionally, in Figure 7 we observe that by changing the parameter $Q$ and if $M < 0$ (weak Allee effect), then system (B.1) also has no equilibrium points, two equilibrium points or one equilibrium point of order two in the first quadrant.

B.1. Stability of positive equilibrium points

To determine the nature of the equilibrium points we compute the Jacobian matrix $J(u, v)$ of (B.1)

with $g(u) = (u+A)(1-u)(u-M)$ and $g'(u) = (1-u)(u-M)+(u+A)(1-u)-(u+A)(u-M)$.

Lemma B.1. The equilibrium point $(0, 0)$ is a saddle point, $(1, 0)$ is a saddle point, $(M, 0)$ is an unstable point and $(0, C)$ is a stable point.

Proof. The Jacobian matrix (B.5) evaluated at $(0, 0)$ and $(1, 0)$ gives

Similarly, the Jacobian matrix (B.5) evaluated at $(M, 0)$ gives

Finally, the Jacobian matrix (B.5) evaluated at $(0, C)$ gives

Since $M < 1$, it follows that $(0, 0)$ and $(1, 0)$ are a saddle point, $(M, 0)$ is an unstable point and $(0, C)$ is a stable point in system (B.1).

Next, we consider the stability of the positive equilibrium points $P_{1}$ and $P_{2}$ of system (B.1). The Jacobian matrix (B.5) of system (B.1) at these equilibrium points becomes

Here, the determinant and the trace of the Jacobian matrix (B.6) are given by

Therefore, the sign of the determinant the trace (B.7), and thus the stability of the positive equilibrium points, depends on the sign of $Q-g'(u_{1, 2})$ and $ug'(u)/(u+A)-S$ respectively. This gives the following results.

Lemma B.2. Let the system parameters of (B.1) be such that $\Delta > 0$ (B.7). Then, the equilibrium point $P_1$ is a saddle point.

Proof. Evaluating $Q-g'(u)$ at $u_1$ gives:

Hence, $\det(J(P_1)) < 0$ and $P_1$ is thus a saddle point.

Lemma B.3. Let the system parameters be such that $\Delta > 0$ (B.7). Then, the equilibrium point $P_2$ is:

1. unstable if $S < \dfrac{D}{4\left(1+A+M+G+\sqrt{\Delta}\right)}$; and

2. stable if $S > \dfrac{D}{4\left(1+A+M+G+\sqrt{\Delta}\right)}$,

with $D = \left(1-A+M+G+\sqrt{\Delta}\right)\left(1+A-M-G-\sqrt{\Delta}\right)\left(1-A-M+G+\sqrt{\Delta}\right)$+$2\left(1+A+M+G+\sqrt{\Delta}\right)\left(A-G-\sqrt{\Delta}\right)$.

Proof. Evaluating $Q-g'(u)$ at $u_2$ gives:

Hence, $\det(J(P_2)) > 0$. Evaluating $ug'(u)/(u+A)-S$ at $u_2$ gives

Therefore, the sign of the trace, and thus the behaviour of $P_2$ depends on the sign of $u_2g'(u_2)/(u_2+A)-S$.

Next we analyse the case when $\Delta = 0$ (B.7) then the equilibrium points $P_1$ and $P_2$ collapse such that $u_1 = u_2 = E = (1-A+M+G)/2$. Therefore, system (B.1) has one equilibrium point of order two in the first quadrant given by $(E, E+C)$.

Lemma B.4. Let the system parameters be such that $\Delta = 0$ (B.4). Then, the equilibrium point $(E, E+C)$ is:

1. a stable saddle-node if $S < \dfrac{H}{4(1+A+M+G)}$; and

2. an unstable saddle-node if $S > \dfrac{H}{4(1+A+M+G)}$,

with $H = (1-A+M+G)(1+A-M-G)(1-A+G-M)+2(1+M+G+A)(A-G)$.

Proof. Evaluating $Q-g'(u)$ at $u = E$ gives

Therefore, $\det(J(E, E)) = 0$. Next, evaluating $ug'(u)/(u+A)-S$ at $E$ gives

Therefore, the behaviour of the equilibrium point $(E, E+C)$ depends on the parity of $Eg'(E)/(E+A)-S$.