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

In the present study, we have modified the traditional May-Holling-Tanner predator-prey model used to represent the interaction between least weasel and field-vole population by adding an Allee effect (strong and weak) on the field-vole population and alternative food source for the weasel population. It is shown that the dynamic is different from the original May-Holling-Tanner predator-prey interaction since new equilibrium points have appeared in the first quadrant. Moreover, the modified model allows the extinction of both species when the Allee effect (strong and weak) on the prey is included, while the inclusion of the alternative food source for the predator shows that the system can support the coexistence of the populations, extinction of the prey and coexistence and oscillation of the populations at the same time. Furthermore, we use numerical simulations to illustrate the impact that changing the predation rate and the predator intrinsic growth rate have on the basin of attraction of the stable equilibrium point or stable limit cycle in the first quadrant. These simulations show the stabilisation of predator and prey populations and/or the oscillation of these two species over time.


Introduction
Numerous investigations have extensively studied the interactions between the specialist predator least-weasel (Mustela nivalis) and the field-vole (Microtus agrestis) [22,23,43]. 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 [14]. Moreover, cycles may be affected by the existence of an optimum group size and an Allee effect [31]. Graham and Lambin [21] 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 [21], 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 [42]. The classic model used to represent the predation between these two species is the May-Holling-Tanner predator-prey model [32] 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 [28] and Pielou in [35]. Later, May (1974) in [32] modified the Leslie-Gower model (1) by assuming the hyperbolic functional response (Holling type II) 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 [32,39,42]. The parameter q is the maximum predation rate per capita and a is half of the saturated level [42]. In particiular, it is given by Model (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 [1,5,42].
The original May-Holling-Tanner model (2) was studied, among others, by [38]. The authors proved that all the solutions which are initiated in 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 [16] at the origin, Saez and González-Olivares [38] 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 by [30] and the analysis of a discrete May-Holling-Tanner model was described in [12].
Additional complexity can be incorporated into the predator-prey model (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 [2,3,10]). 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 [40]. The Allee effect may appear in the field-vole population due to the density-dependent habitat selection showed in [27,34]. Moreover, the stabilisation of field-vole populations depends on the variation in reproductive rate and recruitment of the population [36]. 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 (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 ∈ [0, m), and this is referred to as the "strong Allee effect". When m ≤ 0 1 , the per-capita growth rate is positive but increases at low prey population densities and this is referred to as the weak Allee effect [10,13]. 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 [27,34].
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 Alternative food: Allee effect and alternative food: Figure 2: Main and secondary diet of a specialist predator least weasel (mustela nivalis) studied in [33].
et al. [33] 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 [4, 17, 22-25, 42, 46]. 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 [33] and model this characteristic by adding a positive constant c to the environmental carrying capacity for the predator [8]. Therefore, we have made a modification to the prey-dependent logistic growth term in the predator equation, namely nN is replaced by nN + c (5). As a result, we observe that the modified May-Holling-Tanner model with an alternative food source (6) is no longer singular for N = 0.
The aim of this manuscript is to examine and review the mathematical results of the original May-Holling-Tanner model with an Allee effect on the prey (4), an alternative food source for the predator (6) and both modifications together (7). In particular, by using real system parameter values presented in [24,43], we aim to understand the consequent change in dynamics produced by the Allee effects (strong or weak) and/or the presence of an alternative food source. This manuscript also extends the properties of the weasel and field-vole interaction studied in [21] by showing the impact of these modifications in the original Holling-Tanner model. In addition, the analysis of the modified models (systems (4), (6) and (7)) shows a strong connection with the findings presented in [24,43]. Study results suggest that the system does not lose the coexistence and oscillation of the population when a constant alternative food source for weasels is included. Moreover, the survival of the field-vole still depends on the level of predation when strong and weak Allee on the prey is included [21].
The mathematical interaction in systems (2), considering the modifications separately and together, is reviewed further in Section 2. In Section 3, we conclude the manuscript, summarising the results and discussing 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 [21,24,43].

Models
The results derived in [38] are for the nondimensionalised versions of the modified versions of systems (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 [43]. 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 [15]. Figure 6 illustrates the Hopf bifurcation which was detected with MATCONT in the (q, s)-plane with parameter values Π = (r, K, a, n, c) = (4, 150, 6, 0.025, 0.01) where the parameters r, K, a, n and c are fixed, while we will vary q and s, see Table 2.

Original model
The original May-Holling-Tanner predator-prey model with functional response Holling type II (2) was studied mathematically, for instance, in [26,38] and used to model real predator-prey interactions [22,23,37,42]. The bifurcation diagram of the original May-Holling-Tanner model (2) and two typical phase portraits are shown in the first row of Figure 6. The Hopf curve 2 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 (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 [38]. Note that the extinction of (one of) the species is not possible in this model.

Allee effect (strong and weak)
The May-Holling-Tanner model with a strong Allee effect (m > 0) was studied in [6] in which the authors simplified the system by introducing a change of variable and time rescaling. This is given by the function φ :Ω × R → Ω × R, where φ(u, v, τ ) = (N, P, t) = (Ku, nKv, τ u(u + a/K)/rK). Moreover, by setting Note that system (10) with a strong Allee effect (M > 0) has three equilibrium points in the u-axis ((0, 0), (M, 0) and (1, 0)) and at most two positive equilibrium points [6], see the left panel in Figure 3. The stability of the equilibrium points of system (10) considering M > 0 was studied in [6,11,44,45]. The authors proved that the equilibrium point (0, 0) is always stable, (M, 0) is always unstable and (1, 0) is always a saddle point. Moreover, there are at most two positive equilibrium points: one is always a saddle (P 1 ) and the other (P 2 ) can be stable, unstable or stable surrounded by an unstable limit cycle.
For M < 0 there are two equilibrium points on the u-axis ((0, 0) and (1, 0)) and at most three positive equilibrium points, see the right panel in Figure 3. The stability of the equilibrium points of system (10) considering M < 0 is studied partially in Appendix A. The equilibrium points on the axis (1, 0) and (0, 0) are again a saddle point and unstable respectively. The stability of the positive equilibrium point(s) is described in detail in Appendix A and summarised in Table 3. The dark red circle represents a saddle equilibrium point, the grey circles represent an equilibrium point which can be stable or unstable, the black circle represents an equilibrium point which is always stable, the orange circle represents an equilibrium point which is always unstable and the green circle represents an equilibrium point which is the collision of two equilibrium points.  (10) First, we recall that as we used a diffeomorphism preserving the orientation of time the nondimensionalised version (10) is topologically equivalent to the original version of the model (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 (8). The bifurcation diagram of the May-Holling-Tanner model with weak Allee effect (4) (m < 0) and two typical phase portraits are shown in the third row of Figure 6. System (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 [6]. If the parameters (q, s) are located in the hatched green region, then system (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 (4) has three equilibrium points. These equilibrium points can be a saddle, a stable and/or an unstable node. Therefore, system (4) supports the coexistence and/or oscillation of the predator and prey populations. Note that the dynamics of system (4) with m < 0, outside the grey region, are qualitatively similar to the dynamics of system (2) as shown in row one in Figure 6 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.q 2 < q 1 (see Figure 6). Moreover, as in the original May-Holling-Tanner model, the predator and prey populations never become extinct in system (4) and thus there are always conditions in the system parameters where both populations coexist (see hatched green area in Figure 6) and/or oscillate (see blue and grey areas in Figure 6).
The bifurcation diagram of the modified May-Holling-Tanner model with strong Allee effect (4) (m > 0) studied in [6] and two typical phase portraits are shown in the fourth row of Figure 6. The Hopf curve for the model shows that if (q, s) are located in Figure 4: The intersection of the predator nullcline (blue curve) and the prey nullcline (red curve) in the May-Holling-Tanner model with alternative food for predator by changing the parameter Q. The purple circle represents a saddle equilibrium point, the grey circle represents an equilibrium point which can be stable or unstable, the black circle represents an equilibrium point which is always stable, the orange circle represents an equilibrium point which is always unstable and the green circle represents an equilibrium point which is the collision of two equilibrium points. the green area, system (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 6 (q >q 1 ), then system (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 (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.

Alternative food
The model with alternative food for predators corresponds to system (6). This model was studied in [5,9,18] where the authors presented the stability of the equilibrium points and conditions for periodic solutions. We can also simplify the system by introducing a change of variable and time rescaling presented in [5,9,18]. It is given by the function χ : with Q and S defined in (9), into system (6) we obtain System (12) has three equilibrium points on the axis namely (0, 0) which is always stable, (1, 0) which is always a saddle point and (0, C) which can be a saddle point if C < A/Q (see right panel of Figure 4), saddle node if C = A/Q (see middle panel of Figure 4) and attractor if C > A/Q (see left panel of Figure 4). While, there are at most two positive equilibrium points (see also left panel of Figure 4) one is always a saddle and the other can be stable, unstable or stable surrounded by an unstable limit cycle.
The bifurcation diagram of the modified May-Holling-Tanner model (6) and two typical phase portraits are shown in the second row of Figure 6, see also [18,20]. 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 (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 6 (q > q 2 ), then system (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 (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 (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 [18,20]  The purple circle represents a saddle equilibrium point, the grey circle represents an equilibrium point which can be stable or unstable, the black circle represents an equilibrium point which is always stable and the orange circle represents an equilibrium point which is always unstable.
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.

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 (7). This model was studied partially in [7] and in Appendix B . In [7] the authors only proved the stability of the equilibrium points of the system considering strong Allee effect on the prey (m > 0). The system can be simplified by introducing a change of variable and time rescaling presented in [7]. It is given by the function ζ :Ψ × R → Ψ × R, where ζ(u, v, τ ) = (N, P, t) = (Ku, nKv, τ (u + a/K)(u + c/(nK))/rK) [7]. By setting the parameters A, C, Q, S and M defined in (9) and (11), into system (7) we obtain The equivalent system has four equilibrium points on the axis namely (0, 0) which is a saddle, (1, 0) which is a saddle, (M, 0) which is unstable and (0, C) which is stable, see Appendix B. While, there are at most two positive equilibrium points P 1 which is always a saddle and P 2 which can be stable, unstable or stable surrounded by an unstable limit cycle (see Figure 5).
The bifurcation diagram of the modified May-Holling-Tanner model with Allee effect on the prey and alternative food for the predator (7) and two typical phase portraits are shown in the last row of Figure 6. System (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 [6]. The Hopf curve for the system (7) with Allee effect and alternative food shows that system (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 (6), see the second row of Figure 6.
In Figure 6, 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 7 shows the time series of the predator and prey populations in the original May-Holling-Tanner model (7), the May-Holling-Tanner model with alternative food for the predator (6) and the May-Holling-Tanner model with Allee effect on the prey (4). The system parameters Π are as in (8) and (q, s) = (700, 1.25). Note that these values lie inside the solid or hatched green regions in Figure 6. 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 6, then the original model supports Figure 6: The phase plane and bifurcation diagrams of modified May-Holling-Tanner models ( (4), (6), (7) and (2)) for the system parameters Π (8) fixed and created with the numerical bifurcation package MATCONT [15]. oscillation of the solutions and it does not depend on the initial conditions, see the first column of Figure 7. 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 6, 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 7).
If the parameters (q, s) are located in the grey region of Figure 6, 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).

Conclusion
In this manuscript, the May-Holling-Tanner predator-prey model with functional response Holling type II (2) was reviewed. We also considered that the prey population is affected by strong and weak Allee effect (3) and the predator population has an alternative food source (5). Using a diffeomorphism we extend the analysis of a topologically equivalent system (4), (6), (7) and (2).
We have also extended the analysis by using the system parameters presented in [23,42], 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 6. 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 8 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 (4) is reduced, see the top right panel of Figure 8. In particular, the Hopf curve of the model with strong Allee effect shows that system (4) can have up to two positive equilibrium points which can collapse at q =q 1 (respectively for m < 0 at q =q 2 ). Therefore, for q <q 1 system (4) presents two positive equilibrium points and thus the coexistence of both populations is observed [6]. In all panels we used the system parameters Π (8) fixed and created with the numerical bifurcation package MATCONT [15].
In contrast, the Hopf curve of the model with week Allee effect (m < 0) shows that if q >q 2 then system (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 (4) with m < 0 is similar with the dynamics of system (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 (6) and the system parameters Π (8) and c fixed, see the bottom left panel of Figure 8. 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 (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 [18,20]. Note that if q > q 2 the equilibrium point related to the alternative food is globally stable. Therefore, system (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 (7), i.e. a strong Allee effect on the prey and alternative food for the predator, and the system parameters Π (8), m > 0 and c fixed, see the bottom right panel of Figure 8. In particular, the modified system (7) can have up to two positive equilibrium points which can collapse at q =q 1 . Therefore, for q <q 1 system (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 [6,7,38] 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.

A.1 Stability of positive equilibrium points:
To determinate the stability of the equilibrium points of system (10) with M < 0 we consider the trace and the determinant of the Jacobian matrix of system (10) [6]. Then, we consider the stability of these positive equilibrium points P 1 , P 2 and P 3 of system (10) with M < 0 in the interior of Φ. Proof. Evaluating u 2 (2u − (1 − A + M )) − AM at u = G gives: Then, Therefore, the sign of the trace (17), and thus the behaviour of (G, G) depends on the parity of f (G) − S. Proof. Evaluating u 2 (2u − (1 − A + M )) − AM at u = G gives: Then, the sign of the determinant and the trace (17), and thus the behaviour of P 1 depends on the parity of h(G) and f (G)−S.
Then, the sign of the determinant depends on the parity of G − u 2 . Then, the trace (17), and thus the behaviour of P 2 depends on the parity of f (u 2 ) − S. Proof. Evaluating Then, the sign of the determinant depends on the parity of u 3 − G. Then, the trace (17), 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.
To analyse the positive solution(s) of system (13), which is given by we consider three cases of the intersection of the nullcline of system (13). The u-nullclines of system (13) 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 (19) are: such that M < u 1 ≤ E ≤ u 2 < 1, where E = (1 − A + M + G)/2. In particular,

B.1 Stability of positive equilibrium points
To determine the nature of the equilibrium points we compute the Jacobian matrix J(u, v) of (13) J(u, v) = u(u + C)g (u) + (2u + C)(g(u) − Qv) −Qu(u + C) Sv(A + C + 2u − v) S(C + u − 2v)(u + A) , Next, we consider the stability of the positive equilibrium points P 1 and P 2 of system (13). The Jacobian matrix (21) of system (13) at these equilibrium points becomes J(u, u + C) = u(u + C)g (u) −Qu(u + C) S(u + C)(u + A) −S(u + C)(u + A) , Here, the determinant and the trace of the Jacobian matrix (22) are given by det(J(u, u + C)) = Su(u + A)u + C) 2 (Q − g (u)) tr(J(u, u + C)) = (u + C)(u + A) ug (u) u + A − S Therefore, the sign of the determinant the trace (23), 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 (13) be such that ∆ > 0 (23). Then, the equilibrium point P 1 is a saddle point.
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 2 g (u 2 )/(u 2 + A) − S.
Next we analyse the case when ∆ = 0 (23) then the equilibrium points P 1 and P 2 collapse such that u 1 = u 2 = E = (1 − A + M + G)/2. Therefore, system (13) Therefore, the behaviour of the equilibrium point (E, E + C) depends on the parity of Eg (E)/(E + A) − S.