Predator–Prey Models: A Review of Some Recent Advances

: In recent years, predator–prey systems have increased their applications and have given rise to systems which represent more accurately different biological issues that appear in the context of interacting species. Our aim in this paper is to give a state-of-the-art review of recent predator–prey models which include some interesting characteristics such as Allee effect, fear effect, cannibalism, and immigration. We compare the qualitative results obtained for each of them, particularly regarding the equilibria, local and global stability, and the existence of limit cycles.


Introduction
Population dynamics is one of the most widely discussed topics within biomathematics. The study of the evolution of different populations has always been of special interest, starting with populations of a single species, but evolving to more realistic models where different species live and interact in the same habitat. Among them we can find models that study competitive relationships, symbiosis, commensalism or predator-prey relationships.
We will focus on predator-prey systems, and our aim is to give a state-of-the-art review of recent predator-prey models which include Allee effect, fear effect, cannibalism and immigration. We compare the qualitative results obtained for each of them, particularly regarding the equilibria, local and global stability and the existence of limit cycles.
Probably the most famous predator-prey models are Lotka-Volterra systems, which are polynomial differential equations of degree two that were initially proposed by Alfred J. Lotka [1] in 1925 and Vito Volterra [2] in 1926, both in the context of competing species.
At first, Lotka proposed a model for the study of an herbivorous animal species and a plant species, and then he developed the application to the study of the dynamics of a predator-prey system. On his behalf, Volterra considered the same model simultaneously to explain the observations about the percentage of predatory fish caught during the years of World War I, which had increased. This fact seemed confusing, because the fishing effort had been reduced during those years, so Volterra was interested in studying that situation.
To explain this behavior, Volterra stated a simple system of ordinary differential equations. He considered x(t) and y(t) the densities of prey and predators, and assuming linear behaviors and taking positive constants a, b, c and d, he stated the following model x = x(a − by), y = y(−c + dx).
Later, Lotka-Volterra systems were generalized and considered in arbitrary dimension, i.e.,ẋ i = x i a i0 + n ∑ j=1 a ij x j , i = 1, ..., n. (1) Consequently, the applications of these systems started to multiply. New applications on population dynamics had been developed, but also these systems have been used for modeling many natural phenomena, such as chemical reactions [3], plasma physics [4] or hydrodynamics [5] just as other problems from social science and economics [6]. In view of this, it is clear that population models, in particular competition or predator-prey models, are important both in their original field and in their application to many other problems from other areas. Their study is also interesting from a theoretical point of view. Actually, the study of predator-prey models has always been one of the hot spots in biomathematics, so there is a very big number of works on this topic, especially in dimension two, i.e., with one predator species and one prey species. For this reason, we believe that this review may be of particular interest to those researchers involved in modeling this kind of biological and ecological problems, allowing them to have an overview of the work who is currently being carried out.
Therefore, our objective in this paper is to analyze several adaptations of predator-prey models that have been done in recent years to model different real situations that appear in the field of population dynamics.
We have decided to focus on some special issues to which researchers have devoted special attention in recent years, particularly the study of the Allee effect , the influence of fear effect , cannibalism and immigration. We introduce some papers with influences between them, which allow us to understand the importance of these biological aspects, through the study of several models that are completed and improved in successive works.
The Allee effect appears when a population of prey has a really small density, so it makes difficult for them to reproduce or survive. It was first introduced by Allee in [7]. There are several works in the literature that analyze this effect in different population models and conclude that it can have important effects on the system dynamics, including the stabilization or destabilization of a system. It also can cause the solutions of a system to take a much longer time to reach a stable equilibrium point.
Another issue to consider is that in certain ecosystems, prey may feel fear of predators and act accordingly making hunting more difficult for the predators. The theoretical reasonings about the effect of fear behaviors are supported by real experiments. For example, Zanette , White, Allen and Clinchy [8] conducted an experiment on song sparrows during a whole breeding season using electrical fence to eliminate direct predation of both juvenile and adult song sparrows. No direct killing can happen in the experiment; however, broadcast of vocal cues of known predators in the field was employed to mimic predation risk. Two groups of female song sparrows were tested, among which one group was exposed to predator sounds while the other group was not. The authors found that the group of song sparrows exposed to predator sounds produced 40% less offspring than the other group. They believe this is because fewer eggs were laid, fewer eggs were successfully hatched, and fewer nestlings survived eventually, and because there were behavioral changes including less time of adult song sparrows on brood and less feeding to nestlings during breeding period. Similar experiments on other birds and vertebrate species reported the same conclusions [9][10][11], i.e., even though there is no direct killing between predators and prey, the presence of predators cause a reduction in prey population due to anti-predator behaviors.
We analyze different models with the presence of fear effect, combined with different functional responses, with both omnivorous and specialist predators, with hunting cooperation, with Allee effect and with prey refuge.
As we mentioned, we will study ecosystems in which predators feed only on one prey species, which we will call specialists, and others in which predators are omnivorous, i.e., they will be able to feed themselves from other resources. A special case of this is cannibalism, i.e., the act of killing and at least partial consumption of conspecifics. Cannibalism actively occurs in more than 1300 species in the nature [12], and it has been mathematically modeled for some ecosystems, as for example the European regional seas ecosystem model (ERSEM) for the North Sea [13].
In the literature, one of the first contributions to the study of this topic was the work of Kohlmeir [14], who considered cannibalism in the predator in the classic Rosenzweig-McArthur model. In general, cannibalism has been considered primarily in predators, and the results agree that it has stabilizing effects, and can cause the survival of a species that would otherwise be driven to extinction. We have selected two works that illustrate well the different effects of cannibalism when we add it to already well-known models, as the Lotka-Volterra and the Holling-Tanner model. In addition, the inclusion of cannibalism is considered not only in the predator species as usual, but also in the prey species.
Just as cannibalism or omnivorism can be adaptation strategies to the lack or shortage of food, other types of strategies or behaviors can be induced by this lack of food or by the hostility of the habitat, such as migration. In fact, most predator-prey systems in the wild are not isolated, so it is important to consider the effects of the presence of some number of immigrants. We will focus on studying how the inclusion of immigration changes the dynamics of the classic Rosenzweig-MacArthur model and the Lotka-Volterra model.
Although we have chosen the previous topics as the thread of this review, in many of the considered works these ones appear combined with other biological characteristics such as prey refugee, hunting cooperation or with different types of functional responses. Down below we make some considerations and give some references about these topics.
First, we recall that in ecology a functional response is the intake rate of a predator as a function of food density. The most usual classification of functional responses is the one given by C. S. Holling, in which three types of responses are considered. The type I functional response assumes that the consuming and hunting rates are linear up to a maximum where they become constant. This linear increase assumes that the time needed by the predator to process food is negligible and that consuming food does not interfere with searching food. This is the functional response used in the Lotka-Volterra model.
The type II functional response is given by f (x) = ax/(1 + ahx), where x denotes the food or prey density, a is the rate at which the predator encounters food items per unit of food density, and h is the average time spent on processing a food item. This functional response assumes that the consumption rate is decelerated with increased population because of a limited capacity in searching and processing food. This is the functional response considered in the Rosenzweig-MacArthur model, which is a classical model with logistic growth [15].
In type III functional response the consumption rate is more than linear at low levels of resource. This is suitable for instance in a population of predators that must learn how to hunt efficiently. This functional response is given by f (x) = nx p /(a p + x p ), where a and n are positive and p is an integer with p > 1.
These three types of functional responses are widely used in population dynamics models and especially in predator-prey models, not only on ODE models but also on discrete models, diffusion models, stochastic or fractional order models. Some interesting recent works that study this kind of functional responses are [16][17][18][19][20][21].
There are other types of functional responses, for example the one proposed by Taylor, also known sometimes as Holling type IV functional response, which is given by the function where b ≥ 0. Other models consider a functional response of ratio-dependent type, in which the considered function is h(x) = a/x, with a ∈ R.
Regarding the existence of prey refuge, it has also special interest in the study of predator-prey populations, and many scholars have made great achievements in this aspect. There are interesting works on this subject apart from the ones we include in this review, i.e., those related to our main topics. As a suggestion for interested readers we propose the following references [22][23][24][25][26].
Finally, we would like to note that the third characteristic that we have mentioned, hunting cooperation on predators, has been investigated by many authors in mathematical modeling, but most of the time independently of the effect of fear on prey [27][28][29][30][31]. In this review, we focus on the combination of both hunting cooperation and fear effect, which we believe is particularly important due to the ecological evidence. For example, wolves cooperate during hunting and when they are present, elks use anti-predator strategies and avoid areas frequented by the wolves [32]. Elks respond to the presence of wolves by altering foraging, vigilance, habitat selection, patterns of aggregation and sensitivity to environment [33][34][35][36]. In the same way, while lionesses show cooperative hunting behaviors [37], zebras are affected by fear of predation risk, reaching areas where the encounter with lionesses is less frequent [38].
In the following, we will focus on models with Allee effect in Section 2, with fear effect in Section 3, with cannibalism in Section 4 and with immigration in Section 5. We compare the proposed models whenever it is possible, and we include figures that allow us to compare different models and to see how the inclusion of different biological characteristics affects the dynamics. All these figures have been made with the free software R. In Table 1 we summarize the most important characteristics of the models studied.

Study of Allee Effect
The Allee effect is the biological phenomenon of correlation between the population size or density and the growth rate. Roughly speaking, it occurs when the population of a species has a really small density, so it makes difficult for them to reproduce or survive.
There are two types of Allee effect depending on the nature of density dependence at low densities, the strong Allee effect and the weak Allee effect. The distinction between these two types is based on whether or not the population exhibits a critical population size or density. A population with a weak Allee effect has a reduced per capita growth rate at lower population density or size, but even at this low population size or density, the population will always exhibit a positive per capita growth rate. On the other hand, a population with a strong Allee effect has a critical population size or density under which the population growth rate becomes negative. Then when the population density or size is below this threshold, the population will be drive to extinction.
None of these types of Allee effect were considered in the first classical models as the Lotka-Volterra model. This phenomenon was first introduced by Allee in [7]. There are several works in the literature that analyze this effect in different population models [39][40][41] and conclude that it can have important effects on the system dynamics. These effects include stabilizing or destabilizing a system by changing the local stability of some singularities, as can be seen in [42][43][44][45][46] or cause the solutions of the system to take a much longer time to reach a stable equilibrium point [44,47].
Mathematically, the Allee effect is usually represented modifying the growth function. The most common modification is to use a multiplicative factor, in which case the equation for a single species is given byẋ where r is the intrinsic growth rate and K the environmental carrying capacity. It is said that the Allee effect is strong if m > 0 and weak if m ≤ 0. For the strong Allee effect, the per capita growth rate is negative for x ∈ (0, m), which implies the extinction of the population.
In this section, we will focus first on some Lotka-Volterra systems subject to an Allee effect and then we will study the inclusion of this effect in more general population models.
First, in [48] the author begins by studying the following Lotka-Volterra systeṁ where the variables x and y represent the number of individuals of the prey population, and the number of individuals of the predator population, respectively. We will maintain this notation throughout all the review. The equilibrium points for this system are (0, 0), (1, 0) and (x * , y * ) = (r/(a + r), r/(a + r)). The most interesting singularity from a biological point of view is the one in the positive quadrant, which represents the coexistence of both species. Based on the eigenvalues of the Jacobian matrix, the author concludes that the equilibrium on the axes is unstable and the positive equilibrium is locally asymptotically stable. The interest of this work lies in the comparison with the case in which Allee effect is taken into account. For that, the following system is considered:ẋ where α(x) = x/(β + x) represents the Allee effect and β > 0 is called the "Allee effect constant". We understand that as β increases, the existing Allee effect on the population is stronger, and the population growth of that species decreases. Based on biological facts the function α(x) must verify some conditions: α (x) must be positive for any positive x because the Allee effect decreases as density increases and lim x→∞ α(x) = 1 because the Allee effect vanishes at high densities. The authors also state that α(0) must be zero because reproduction is not possible without any individuals, but as the net growth rate is xα(x), the condition α(0) = 0 is not necessary as xα(x) is always zero if x = 0, regardless of the value of α(0).
The proposed function α(x) = x/(β + x) verifies this set of conditions, but maybe it would be interesting to analyze different functions representing the Allee effect.
For system (5) there exist again three equilibria. Two of them are the same as in the system with no Allee effect, (0, 0) and (1,0), and the third one acquires a new expression (x * , y * ) = ((r − aβ)/(a + r), (r − aβ)/(a + r)). It is important to recall that now, the equilibrium (x * , y * ) has only biological meaning if r − aβ is positive.
The stability of the equilibria is studied by means of the eigenvalues. The author says that the equilibria (0, 0) and (1, 0) are unstable, but the proof for the origin is not given. Since it is a linearly zero equilibrium, i.e., the linear part of the system (5) in the origin is the identically zero matrix, it would be convenient to give a more detailed analysis. The equilibrium (x * , y * ) is locally asymptotically stable.
Numerical simulations are also carried out in the article. The outstanding conclusions of [48], obtained from both analytical and numerical considerations, are the following: The Allee effect may cause there to be no equilibrium in the positive quadrant, while in the classical model this singularity always exists, and furthermore it is locally asymptotically stable. In the case that the Allee effect does not affect the existence of equilibria, it also has consequences on the dynamics of the model; on the one hand, the system takes much longer time to reach the state of equilibrium, and on the other hand, the values of population density of both prey and predators in that point are reduced with respect to the classical model.
In [49] the study of the case in which the Allee effect affects predator population rather than prey is carried out, based on the reason that higher levels in the food chain are more likely to become extinct. The proposed system is analogous to the one in [48], considering the same function α, but in this case as a function of the density of predators: In this case, for the study of the stability, different results are combined: the study of the trace and determinant of the Jacobian matrix of the system, the Bendixon-Dulac criterion and the property of persistence. This notion of persistence or permanence appears frequently in ecology as it has special importance on biological systems because it is a condition ensuring the long-term survival of all species. Considering a general systeṁ on the positive region R n + , where x = (x 1 , ..., x n ) and conditions ensuring the global existence and uniqueness of solutions in forward time hold, it is called permanent or persistent if there exist m, M ∈ (0, ∞) such that given any x ∈ intR n + , there is a t x such that From a biological point of view, it is reasonable to expect that a permanent system will have a coexistence equilibrium in intR n + , as can be seen in [50]. Because of the importance of this property, the first objective addressed by the authors in [49] is the proof of the persistent property of system (6), because they rely on that to prove the local stability. They conclude that the system is persistent under condition r > a.
The persistence property makes it possible to ensure that the equilibria on the axes, (0, 0) and (1, 0) are unstable, as no positive solution can approach to these equilibria. For the equilibrium point in the positive quadrant (x * , y * ) = (r(1 − x)/a, r(1 − x)/a), they conclude that it is locally asymptotically stable, since the real part of the eigenvalues is negative. Locally, the stability of the equilibria is the same as the obtained in [48] for the system with Allee effect on the prey species.
Here, beyond the local study of equilibria, a simple argument is used to prove that the asymptotic stability of the point (x * , y * ) is global. They prove that there are no limit cycles, which is easily deduced applying the Bendixon-Dulac criterion, and that is enough because we know (x * , y * ) is the only stable equilibrium. This global result is not proved for the model in [48].
After the stability analysis, numerical simulations are carried out in [49]. They show that if r > a, the Allee effect on the predator species has no influence on the existence and stability of the positive equilibrium and on the final density of the predator and prey species. This result differs from the obtained in [48]. This is showed in Figure 1. A similar behavior is also obtained in some cases with r < a, as for example in the case given in Figure 2.  Numerical simulations also show that the system with Allee effect on the predators takes much longer time to reach its positive equilibrium when β increases its value, which coincides with the results in the case of Allee effect on the prey species.
In [51] the previous model is combined with the choice of a more realistic model for the growth of prey, instead of the logistic one used previously. The proposed model iṡ In particular, they take the Beverton-Holt function [52] as the birth rate. In the same way that in [49], the authors prove that the system is permanent. In this case, the condition under which the system (7) is permanent is: As the authors say, this condition can be changed into other which is simpler but more restrictive Analogous to the previous system, there exist three equilibria: (0, 0), (u * , 0) and (x * , y * ), and the conclusions about their stability are similar: if the condition that makes the system permanent holds, then the equilibrium points on the axes are unstable, and the positive equilibrium, which appears if a 1 > b 1 d 1 , is always locally asymptotically stable, and globally asymptotically stable if the condition for the permanence holds.
Here an extinction property is proved: if This kind of property was not proved in the systems studied in [48,49]. At last, the numerical simulations in this case are only oriented to give some examples for the validation of the conclusions about the stability, but the influence of the Allee effect is not analyzed.
We will now consider the inclusion of Allee effect in more general predator-prey models. In [40] the authors study the global dynamics of a model with increasing and bounded functional response which includes a strong Allee effect. They consider a systeṁ where the functions f and g verify the following hypothesis: The function g is the functional response of predator and g(x) f (x) is the net growth rate of the prey. The conditions about g include Holling type I, II and III functional responses, but no type IV.
The hypothesis on f reflect that in the absence of the predator, the prey has a strong Allee effect growth. The Equation (10) is dimensionless and the conversion efficiency of prey biomass into predator biomass is 1. The parameter b is the survival threshold or sparsity constant of the prey. The parameter d is the mortality rate of predator but also λ, which is the stationary prey population density coexisting with predator, can be thought as a measure of the predator mortality as it increases with d. As pointed out by Bazykin [53], it is natural to regard λ as a measure of how well the predator is adapted to the prey. In [40] the parameter λ plays an important role as a bifurcation parameter.
A complete theoretical analysis of the model is carried out, including an analysis of the phase portrait and Hopf bifurcation. The system is transformed into other of Liénard type. The existence and uniqueness of limit cycles and heteroclinical orbits are studied, which is important for some ecological phenomena. The first quadrant is shown to be invariant.
There exist four possible equilibria: (0, 0), (1, 0), (b, 0) and (λ, f (λ)). This latest equilibrium only exists if b < λ < 1, and in any other case the system has three boundary equilibria, there is no coexistence and no limit cycles. In that case the origin is globally asymptotically stable if 0 < λ ≤ b, and if λ ≥ 1 the regions of attraction of the origin and (1, 0) divide the first quadrant.
If b < λ < 1 then the origin is a stable node, (b, 0) is a saddle with stable manifold on the x-axis and unstable manifold on the prey nullcline y = f (x). The equilibrium (1, 0) is a saddle with stable manifold on the x-axis and unstable manifold on the region y > f (x).
The authors also analyze the existence and multiplicity of limit cycles, and for some conditions the uniqueness of limit cycle is proved.
The system has a Hopf bifurcation when λ = λ in (λ, f (λ)). The Hopf bifurcation can be supercritical and backwards (respectively, subcritical and forward) if the first Lyapunov coefficient is negative (respectively positive). Hopf bifurcation is backwards if there exists a periodic orbit of small amplitude for each λ ∈ (λ − ε, λ) (similarly with forward), it is supercritical if the periodic solutions which bifurcate are orbitally asymptotically stable (respectively, subcritical if periodic solutions which bifurcate are unstable).
The key bifurcation parameter λ is the location of the predator nullcline x = λ. When λ is below the prey extinction threshold b, the prey will become extinct and consequently the predator will also starve to death. The authors show that there exists a threshold λ * < λ such that if λ < λ * , then the predator invasion leads to the extinction of both species. This already studied phenomenon is known as overexploitation [54], and mathematically it means the global stability of the equilibrium (0, 0). The threshold λ * is at the unique point where a point-to-point heteroclinic orbit loop exists.
When λ moves across λ = λ * , a limit cycle appears. It can be established if initially the prey is over the threshold b, and the predator is under the overexploit threshold to make a successful invasion. There is another change when it moves past λ, as it is the value at which the Hopf bifurcation takes place. Then, the alternative stable state of limit cycle switches to a coexistence equilibrium point, and the oscillation dies down to the stable equilibrium. Finally, when λ > 1, the alternative stable state becomes the equilibrium in which there are only prey, and the predator species is not able to establish itself and will become extinct. For these cases, higher initial predator density will always lead to the overexploitation.
The authors apply these results to special cases determined by different functional responses g(u), and they study concrete models, setting the expressions of both functions f and g. In particular, they apply their results to a cubic model with Holling type II functional response, a cubic model with linear functional response and to the Boukal-Sabelis-Berec model [55].
After analyzing this more general model, we study another one that cannot be written in the form (10). In [39] the studied model is based on the Leslie-Gower model, i.e., the predator growth is given by a logistic-type equation and the environmental carrying capacity of predators is a function on the prey population, f (x) = nx. The Leslie-Gower model has only one positive equilibrium which is asymptotically stable. Adding the prey population affected by Allee effect and the lineal functional response, the authors propose the following system:ẋ The parameter m is the Allee threshold or minimum of viable population, they consider m ∈ [0, K), resulting in a weak Allee effect when m = 0 and in a strong Allee effect in any other case.
The system (11) has always the boundary equilibria (K, 0) and (m, 0). The positive By changing the variables by x = Ku, y = nKv and rescaling the time τ = trK/(u + c/(nK)), the authors obtain the following Kolmogorov system, which is topologically equivalent and a continuous extension of system (11) Here M = m/K, Q = qn/r, C = c/nK and S = s/rK. The equilibrium points of system (12) are (0, 0), (M, 0), (1, 0) and (u * , u * ). For this system the authors obtain that for certain values of the parameters there are two positive equilibria (u * 1 , u * 1 ) and (u * 2 , u * 2 ), only one positive equilibrium or none.
In the case with strong Allee effect, i.e., with M > 0, the equilibrium (M, 0) is a hyperbolic repellor point and (1, 0) is a hyperbolic saddle point. The point (0, 0) has a hyperbolic and a parabolic sector, so there exists a separatrix curve in the phase plane that divides the behavior of the trajectories. This implies that the model is highly sensitive to the initial conditions. Then, this point is an attractor point for certain trajectories and a saddle point for others. In the proof of this result the blowup technique is applied, which can be found in [56,57].
When the system has two positive equilibria, (u * 1 , u * 1 ) is a saddle point and (u * 2 , u * 2 ) can be a stable point or an unstable equilibrium surrounded by a stable limit cycle. In the transition from one case to another, the equilibrium (u * 2 , u * 2 ) is at least a two-order weak focus and a stable limit cycle is generated by a Hopf bifurcation.
There exist values of the parameters for which there is a homoclinic curve determined by the stable and unstable manifold of the point (u * 1 , u * 1 ) and there is an unstable limit cycle that bifurcates from the homoclinic orbit surrounding point (u * 2 , u * 2 ). For some values of the parameters the equilibria (u * 1 , u * 1 ) and (u * 2 , u * 2 ) collapse and give rise to a unique equilibrium point in the interior of the positive quadrant. This point can be a non-hyperbolic attractor node, a non-hyperbolic repellor node or a cusp, in which chase there exists a unique trajectory which attains the equilibrium. When this point is a repellor, the origin is a global attractor.
In the case with weak Allee effect, i.e., with M = 0, the equilibrium points (0, 0) and (1, 0) are a non-hyperbolic saddle and a saddle, respectively. Moreover, there exists a unique positive equilibrium if and only if Q < 1. This equilibrium can be a local attractor or a repellor, and in the transition from one case to another, the equilibrium is a weak focus of order one. The equilibrium can be also a repellor focus surrounded by a limit cycle or a repellor node (in which case the origin is globally asymptotically stable). This work shows that the Allee effect may destabilize a predator-prey model, since the equilibrium points of the system could be changed from stable to unstable.
In the previous works a simple Allee effect was considered, but mathematically it is possible to represent the Allee effect with more than one modification on the system simultaneously. In [58] two models are compared, one with a unique Allee effect and other with a double Allee effect on prey. The model with double Allee effect is given bẏ where n < K for biological reasons and m ∈ (−K, K). The Allee effect is expressed once in the factor α(x) = x − m, and a second time in the term r(x) = rx/(x + n) which can be interpreted as an approximation of a population dynamics where the differences between fertile and non-fertile are not explicitly modeled. The one with a unique Allee effect is given bẏ We recall that this last system is a particular case of (10).
For system (13) there exist always three equilibria, two of them in the boundary. In the case with strong Allee effect there exists another boundary equilibrium.
The authors prove the existence of a subset of parameter values for which there are two limit cycles when the Allee effect is either strong or weak in system (13). This is a different result by comparing with system (14), generated by the double Allee effect, because system (14) has a unique limit cycle generated by Hopf bifurcation, surrounding the unique positive equilibrium.
In system (13) there exists a separatrix curve determined by the unstable manifold of equilibrium point (m, 0) or the origin. Then, there are trajectories close this separatrix with different ω-limit sets for the same values of the parameters, so the system is highly sensitive to initial conditions. For a same set of parameters different behaviors are possible: the extinction of two populations, the coexistence, and the oscillation of both populations. This is in accordance with the results in [39].
The separatrix curve disappears when m < 0 and in this case the origin is a saddle point in both systems (13) and (14). We give an example in Figure 3 of how the systems (13) and (14) show an analogous behavior but in the system with double Allee effect the predator density in the coexistence equilibrium is smaller than in the case with simple Allee effect.  In Figure 4 we compare the phase portraits of systems (13) and (14) with the same values of the common parameters. We note that for the system (13) there exist two limit cycles surrounding the positive equilibrium, while it does not occur for the system with simple Allee effect.  With similar techniques, the influence of the Allee effect is studied in many other works, combined with different particular biological conditions. For example, the influence of the inclusion of different characteristics in Gauss-type predation models has been studied. In particular, in [59] the competition among predators is included, and in [60,61] the hyperbolic functional response, which is a particular case of Holling type II functional response, is considered. The authors study both the strong and weak Allee effect. A generalization of this last work is given in [62], with different functional responses. In this last one, the results show a marked difference between cases of strong and weak Allee effect, in contrast with [59] where there exist minor differences among both cases.
Instead of assuming a functional response of type I, as in the previous models, now we are presenting a model with type II functional response. This is the case of the work of D. Sen, S. Petrovskii, S. Ghorai and M. Banerjee, [63]. They also consider generalist predators, i.e., the ones that feed not only on the species of prey. The predator species is subject to an Allee effect, and this model is motivated by the real-world example of the phytoplankton-zooplankton system. Although the zooplankton growth is known to be subject to the Allee effect, there is considerable evidence that the phytoplankton growth is logistic [64].
They show that increasing the carrying capacity of preys helps the population to survive in a coexistence equilibrium. Their model can produce bi-stable dynamics for a reasonable range of parameter values. The obtained dynamics is quite rich and there exist different types of bifurcation. The system they consider is: which has a Holling Type II functional response. Taking x = N/k 1 and y = P (r 2 /(k 2 r 1 )) they provide a dimensionless system: where α = (1/(hk 1 )) k 2 /(r 1 r 2 ), δ = 1/(ahk 1 ), η = √ r 2 k 2 /r 1 , l = m r 2 /(k 2 r 1 ) and α 1 = e/(r 1 h).
This system has the trivial equilibrium E 0 = (0, 0) and three boundary equilibria E 1 = (1, 0), E 2 = (0, l) and E 3 = (0, η). The system can also present up to four positive equilibria E * j obtained as the intersection points of non-trivial prey and predator nullclines in the first quadrant: Obviously, the number of interior equilibria depends on the value of the parameters, shape and position of the nullclines and position of the points of intersection of the nullclines on the positive axes with respect to the boundary equilibria.
Stability conditions for all boundary equilibria are obtained by applying standard linear stability analysis, but for the positive equilibria the authors make a graphical approach, since the expressions of E * j are not available explicitly. For the boundary equilibria they prove that E 0 is always a saddle; E 1 is a saddle if α 1 > (1 + δ)ηl and locally asymptotically stable if α 1 < (1 + δ)ηl; E 2 is a saddle if δ/α < l and an unstable node if δ/α > l; E 3 is locally asymptotically stable if δ/α < η and it is a saddle if δ/α > η. They also give a description of the basin of attractions of the equilibria.
Regarding bifurcations appearing on the system, the authors proof the existence of saddle node and transcritical bifurcations through which the positive equilibria appear or disappear and the stability of the boundary equilibria changes. Their results are the following: The point E 1 undergoes a transcritical bifurcation when δ passes the threshold δ TC1 = α 1 /(ηl) − 1 for α 1 > ηl. E 1 undergoes another transcritical bifurcation at the thresholds δ TC2 = αl. At last, E 3 undergoes a transcritical bifurcation at δ TC3 = αη.
They also prove analytically the existence of saddle-node bifurcations, more precisely: the system (16) undergoes a saddle-node bifurcation when the curve y = p(x) touches y = q 1 (x) at E(x, y) such that x is a double root of h(x) = 0 with x ∈ (0, 1) and y ∈ (0, l). System (16) undergoes another saddle-node bifurcation at α = α SN2 when y = p(x) touches y = q 2 (x) at E(x, y) such that x is a double root of h(x) = 0 and y = p(x) = q 2 (x) with 0 < x < 1. The expressions for p(x), q 1 (x), q 2 (x) and h(x) are: As the expressions of the positive equilibria are difficult to find, the authors cannot give the analytical conditions for Hopf bifurcation, but they find the numerical threshold value of the parameter α at which the system (16) undergoes this kind of bifurcation. For this they use the software Maple. They prove numerically that the first Lyapunov coefficient is negative, so the bifurcating limit cycle is stable and the bifurcation is supercritical.
They also study numerically the so-called Bogdanov-Takens bifurcation. It is a local co-dimension 2 bifurcation which occurs at the point where the bifurcation curves corresponding to the saddle-node bifurcation and the Hopf bifurcation intersect when we draw them together in a two parametric bifurcation plane. The authors give the threshold value of this bifurcation.
The study of the complete dynamics of the system (16) carried out in [63] shows with a well-founded analysis that the dynamics of the predator-prey system can be very complex also in systems with only two species. The authors show that the dynamics of the generalist predator-prey system assuming that the predator is subjected to Allee effect is much more complicated than in the case of the specialist predator. The prey-predator system with a generalist predator has more sensible biological properties, because the extinction of prey does not lead to the extinction of predator. Although the specialist predator system has some unrealistic properties, the generalist predator system apparently does not. One remarkable property of the specialist case is its counter-intuitive response to the system enrichment described by an increase in either its linear per capita growth rate or its carrying capacity. The complexity is also increased because the inclusion of the Allee effect in the predator causes different types of bifurcation.
The work done in [63] leaves some open questions as if the complexity of the system may increase considering Holling type III or ratio-dependent functional response instead of the Holling type II functional response.

The Influence of Fear Effect
In a biological system where predators and prey coexist, the behavior of each of these species can affect both populations. The way of life of the prey affects the evolution of the ecosystem, for example when prey, aware of the presence of their predators, feel fear and act accordingly making hunting more difficult for the predators. This fact drives prey to avoid direct predation, which can increase short-term survival of prey, but may cause a long-term decrease in prey population as a consequence [65].
There exist experiments that support the theoretical reasonings about the effect that anti-predator behaviors may have, as the ones mentioned in the introduction [8][9][10][11]. All these experiments on birds and vertebrate species reported the same kind of conclusions: even though there is no direct killing between predators and prey, the presence of predators cause a reduction in prey population due to anti-predator behaviors.
We analyze different models with the presence of fear effect. First, in [66] the influence of fear effect in two different systems is studied. These two systems consider respectively a Holling type I functional response and a Holling type II functional response. In this model predators feed only on prey, so in [67] this model is generalized into one in which predators are omnivorous. Then we present a study of the impact of fear on a Leslie-Gower model with hunting cooperation [68]. As we said in the introduction, the combination of the fear effect and hunting cooperation is supported by much ecological evidence. Next, we present two examples where the fear effect is combined with Allee effect, the first one in [69] which has a multiplicative Allee effect and the second in [70] with an additive Allee effect. Finally, in [71] the influence of prey refuge is studied, including it on a model with Holling type II functional response. The prey refuge is also considered in [72], where it is included in the model proposed in [69] and combined with an Allee effect.
We start with the work of X. Wang, L. Zannete and X. Zou in [66] where they propose the following model,ẋ which modifies a classic predator-prey model multiplying the production term by a factor f (s, y) related to the cost of anti-predator defense due to fear. The parameter s represents the level of fear and by the biological meaning the following conditions are assumed: In [66] two types of functional responses are considered: a Holling type I functional response or lineal response, g(x) = xp, and a Holling type II functional response, g(x) = xp/(1 + qx).
It is assumed that predators feed exclusively on this species of prey. Consequently, only two boundary equilibria appear, since there could never be a population of predators in the absence of prey, and therefore the equilibrium on the vertical axis which appears on the classical model does not exist in this case.
The authors prove that if r < d then both populations become extinct, independently of the fear effect and the particular predation mechanism, because x (t) ≤ (r − d)x, and by a comparison argument and the theory of asymptotically autonomous systems [73], they obtain that x(t) and y(t) tend to zero when t tends to infinity. Therefore, it is enough to consider the case when r > d.
Taking the linear functional response, if r ∈ (d, d + an/cp) then there is not positive equilibrium, i.e., there is not coexistence of both species, and the equilibrium on the boundary, E 1 = ((r − d)/a, 0) is globally asymptotically stable, so predators become extinct. If r > d + an/cp there exists also a positive equilibrium, E * , which is globally asymptotically stable. We note that the threshold value d + an/cp does not depend on the parameter s, so it seems that fear effect has no influence on the dynamics in this case.
With the Holling type II functional response, they consider that the term related to the fear effect is f (s, y) = 1 1 + sy .
If cp ≤ nq, then there is not positive equilibrium and the boundary equilibrium E 1 is globally asymptotically stable. If cp > nq, applying a transformation on the variables t, x and y they provide the Kolmogorov systeṁ where a i > 0 for i = 1, . . . , 6, whose analysis allows them to determine the stability of the equilibrium points of the original system, obtaining the following results. If (r − d)(cp − nq) > an there exist only two equilibrium points, the origin which is unstable and E 1 which is stable. If (r − d)(cp − nq) > an the origin and E 1 are unstable and there exist a positive equilibrium E * which stability depends on the parameters. E * is locally asymptotically stable if or and E * is unstable if When the birth rate of prey is not large enough to support oscillations, there exists a globally asymptotically stable positive equilibrium, and when the birth rate is large enough to support oscillations, the positive equilibrium can be also locally asymptotically stable if prey is sensitive enough to perceive potential dangers and show an anti-predator behavior motivated by fear. So, the conclusion is that in some cases, the fear effect can stabilize the system. It is proved that the asymptotical stability of E * is global if (21) holds and In addition, in the case with Holling type II functional response, it is possible the existence of limit cycles. In this case, the positive equilibrium becomes unstable, and a limit cycle appears as a result of a Hopf bifurcation. The impact of the fear effect on this Hopf bifurcation is also studied in [66].
Although the authors choose f (s, y) = 1/(1 + sy) for the theoretical analysis of the model, they also conduct some numerical simulations with another decreasing functions, in particular: f (s, y) = e −sy and f (s, y) = 1 1 + s 1 y + s 2 y 2 .
The behavior observed in the simulations regarding the existence of Hopf bifurcation with these other functions is similar to that obtained for the first chosen function. In general, the same results are obtained for a general monotone decreasing function of y.
The model analyzed in this paper assumes that perceived predation risks only reduce the birth rate and survival of offspring, and does not include the possible impact on the mortality rate of adult prey. They also do not include the case where fear affects intraspecific competition. The reason for not incorporating them in the models is that there is not yet experimental evidence for both effects. However, in [8,74] they argue that fear can increase the mortality rate of adults due to long-term physiological impacts. There are also theoretical arguments in [65] that the effect of fear may change the strength of intraspecific competition due to the complexity of the food web.
In the previous model, predators feed only on the prey species, but recently, Z. Zhu, R. Wu, L. Lai and X. Yu in [67] have proposed a model where predators are omnivorous. As we saw in Section 2 for the case of populations with Allee effect, the consideration of omnivorous or generalist predators affects the dynamics of the system. Here the authors assume predators can feed on other species, particularly when the main prey species becomes extinct. This characteristic is incorporated on the system with lineal functional response considered in [66], by imposing a logistic growth of predators in absence of prey obtaining:ẋ Depending on the values of the parameters, the system can exhibit four equilibria. Two of them are the boundary equilibriums that also appeared on the system without omnivorous predators, E 0 and E 1 , and now it appears also a new boundary equilibrium E 2 = (0, n/d 1 ), because predators can survive in absence of prey. Here, the equilibrium points E 0 and E 1 are unstable, so the dynamic behavior of both systems are quite different. Furthermore, the equilibrium E 2 is globally asymptotically stable under certain conditions which shows that the predator species can be permanent despite the extinction of the prey species. When the positive equilibrium exists, it is globally asymptotically stable, as in the case with no omnivorous predators.
In addition, the authors conclude that the density of both prey and predator species decrease when the fear parameter s increases. In Figure 5 we show an example of how the variation of the fear parameter affects the solutions and the final population density of each species in the equilibrium point. It is concluded that fear effect can cause the extinction of the species, in opposition to the conclusions of [66], where the fear did not affect the dynamics. In Figures 6 and 7 we give examples of the different behaviors obtained for systems (25) and (19) (with linear response) with the same values of the common parameters, so it is clear that the inclusion of the omnivorism affects the dynamics.  Considering the importance of the predator predation strategy and the prey avoidance strategy it is convenient to analyze how both prey and predators adjust their behaviors to obtain maximum benefits and increase their biomass for each of them. In [68], S. Pal, N. Pal, S. Samantha and J. Chattopadhyay consider a modified Leslie-Gower predator-prey model where predators cooperate during hunting and due to fear of predation risk prey populations show anti-predator behavior. The impact of hunting cooperation and fear effect on the dynamics of the system are analyzed.
The authors consider the following system with Holling functional response of type II, fear effect and Berec's encounter-driven functional response [27]: Please note that in absence of hunting cooperation (p = 0), the previous model coincides with the modified Leslie-Gower model.
It is proved that solutions of system (26) starting at positive conditions are bounded if b > d and q > (b − d + ε 1 ) with ε 1 > 0. Moreover, system (26) The model has the following boundary equilibria: E 0 = (0, 0), E 1 = (b − d, 0) and E 2 = (0, r/q). It can also present a positive equilibrium E * = (x * , y * ). The equilibrium points E 0 and E 1 are both unstable. The equilibrium E 2 is locally asymptotically stable if the cooperation strength α is greater than the threshold value or if the fear strength β is greater than the threshold value bpq 3 r 2 (αr + q) + dpq 2 r − q r .
Using the Routh-Hurwitz criterion, the positive equilibrium E * is locally asymptotically stable if Regarding the existence of bifurcations, this system presents a Hopf bifurcation at α = α h from the equilibrium point E * if and only if D(x * (α), y * (α)) > 0 and d dα T(x * (α), y * (α)) = 0 in α = α h .
To obtain the stability and direction of the Hopf equilibrium, the first Lyapunov coefficient l is calculated, applying the results on [75]. They obtain that this bifurcation can be supercritical or subcritical depending on the parameters.
We note that when l = 0 the system exhibits a generalized Hopf bifurcation, which is also called Bautin bifurcation, at which the interior equilibrium has a pair of purely imaginary eigenvalues. The Bautin bifurcation point separates branches of subcritical and supercritical Hopf bifurcation in the parameter plane.
The Bogdanov-Takens bifurcation from the equilibrium point E * is also analyzed. The employed techniques are given in [76].
Some of the important conclusions of this work are the following. In the absence of fear, if we increase the strength of hunting cooperation, then the system becomes unstable around the positive equilibrium and undergoes multiple Hopf bifurcations, the first of them is supercritical and the second one is subcritical. The coexisting equilibrium can lose the stability via supercritical Hopf bifurcation with an intermediate value of hunting cooperation and regain stability via subcritical Hopf bifurcation if the strength of hunting cooperation is large.
The existence of multiple limit cycles and bistability via subcritical Hopf bifurcation has been observed by the authors. The solutions tend to a coexisting equilibrium or oscillate periodically depending on the initial population size.
When the system shows limit cycle oscillations around the interior equilibrium, increasing the strength of hunting cooperation or cost of fear, then the population oscillations will be replaced by a stable focus.
It is concluded that fear factor has a stabilizing effect, and by including it, the stable system remains stable, and the oscillating system becomes stable by excluding the existence of periodic orbits. It is observed that fear factor has more stabilizing effect than hunting cooperation and makes the system more robust.
For low level of fear, if the cooperation level is increased then the system switches its dynamics from stable to limit cycle oscillation and from limit cycle oscillation to stable focus. When the strengths of hunting cooperation and fear factor are very high, prey population goes to extinction due to overexploitation by the predators. In general, increasing both hunting cooperation and fear factor, the density of both species populations eventually will decrease.
In conclusion, the system with hunting cooperation and fear phenomena exhibits rich dynamical behaviors so these two characteristics are interesting factors with an important role in determining the long-term population dynamics.
Inspired by [66], S. K. Sasmal considered in [69] a model with a reduction of prey growth rate due to the fear effect and with the prey species affected by Allee effect. This model is given byẋ where s represents the effect of fear and m the Allee effect. The authors work on a dimensionless system obtained through variables changes, and obtain the following results for the original system. The system presents three boundary equilibria which always exist, E 0 = (0, 0), E 1 = (m, 0) and E 2 = (K, 0), and a positive equilibrium E * which appears under condition m < n aα < s.
The origin is always locally asymptotically stable. The boundary equilibrium E 1 is an unstable node if m > n/aα and it is saddle if m < n/aα. The boundary equilibrium E 2 is a saddle if n/aαs < 1 and it is locally asymptotically stable if n/aαs > 1. The equilibrium E * is locally asymptotically stable if (s + m)/2 < n/aα < s and it is unstable if m < n/aα < (s + m)/2.
We note that the origin is asymptotically stable, so there exists always some initial conditions that will lead to extinction of both species. Due to strong Allee effects in prey, when the predator invasion is large enough, then predator grows very fast to drive prey below its survival threshold density, and extinction of both the populations is certain.
The analysis carried out concluded that the stability of the positive equilibrium does not depend on the fear effect, the same showed in [66], but here the cost of fear changes the density of the predator species on the coexistence equilibrium: as the value of s increases, the density of predators on the positive equilibrium decreases.
The system presents a subcritical Hopf bifurcation at n/aα = (s + m)/2, and changing the values of m and n can produce bistability or stable oscillatory coexistence of both species. The author observed that modifications on the value of s can only change the density of predator species at the positive equilibrium, but it does not affect the stability.
Motivated by this work, L. Lai, Z. Zhu and F. Chen, in [70] replace the usual Allee effect with the additive Allee effect, which they take from [77]. They study the systeṁ where the term 1/(1 + sy) represents the fear effect and the term m/(x + b) represents the additive Allee effect. In particular, the values of s and m represent, respectively, the intensity of fear effect and Allee effect. Depending on the parameters the system can present two or three boundary equilibrium points. The origin E 0 is always an equilibrium for system (29). It is a stable node if b < m or b = m = 1, a saddle node if b = m = 1 and a saddle is b > m. When b ∈ (0, 1) and b < (b + 1) 2 /4 < m there are no more boundary equilibria. Now we consider the different cases depending on the existence of other boundary equilibria besides E 0 and their local stability.
Besides the boundary equilibria, the system (29) has a positive equilibrium , , and if also E 3 exists and is unstable, the asymptotical stability is global. This allows us to know under what conditions the system always evolves towards the coexistence of both species.
The authors give conditions for the existence of saddle-node bifurcation and transcritical bifurcation of boundary equilibriums. They apply the Sotomayor's theorem, which can be found in [78,79]. The saddle-node bifurcation occurs from E 1 at m = (b + 1) 2 /4 if b < 1 and αa(1 − b) − 2n = 0. The transcritical bifurcation occurs from E 0 at b = m when m < 1. They also study the supercritical Hopf bifurcation of the positive equilibrium that occurs when b = √ m − x * . In general, the additive Allee effect gives rise to a more complex dynamics than the multiplicative Allee effect. As an example, for a fixed set of parameters we show in Figure 8a comparation between the behaviors of system (29) and system (27). In the system with multiplicative Allee effect, (27), both species become extinct while in system with additive Allee effect, (29), for the same common parameters and varying the additive Allee effect coefficient b we obtain different dynamics. Regarding the influence of fear effect, this work conclude that it only affects the final density of predator species. Now we will focus on some models that include a prey refuge. First, in [71], H. Zhang, Y. Cai, S. Fu and W. Wang, investigated a Holling-II predator-prey model incorporating the fear effect and a prey refuge. Their results show that the impact of prey refuge and fear effect is significant. They based their study on the model proposed by T. K Kar in [22], i.e.,: on which they add the fear effect by multiplying by the factor f (s, y) = 1/1 + sy, obtaining the modelẋ where s > 0 refers to the level of fear, and b = α/K. This system represents what they call the anti-predator defense due to fear, and for proposing it, they are also inspired by the results on [66]. The origin E 0 = (0, 0) is always an equilibrium for system (31) and it is always a saddle point. There is also a boundary equilibrium E 1 = (α/b, 0) that always exists. This equilibrium is locally asymptotically stable if cβ − aγ ≤ 0, or cβ − aγ > 0 and 1 − bγ/α(cβ − aγ) < m < 1, so under these conditions there exists a region from which the system always goes to the extinction of the predator species. E 1 is unstable if cβ − aγ > 0 and 0 ≤ m < 1 − bγ/α(cβ − aγ). We note that these conditions do not depend on the parameter s, i.e., they do not depend on the fear effect. When E 1 is unstable, it appears a coexistence equilibrium . This positive equilibrium is locally asymptotically stable if and The asymptotical stability is global if 0 < cβ − aα ≤ aα and the first condition in (32) holds. At last, E * is unstable if the second condition in (32) holds and 0 < s < S * .
The existence of Hopf bifurcation and limit cycles is also studied in this paper. If the second condition in (32) holds, then the system has a Hopf bifurcation from E * at s = S * , and if 0 < s < S * then it is proved that there exists a limit cycle. The Poincaré-Bendixon theorem is used in the proof of this result.
There are some important conclusions, some that agree with other works and others that yield different results. When the coexistence equilibrium exists, the increase of the fear parameter s causes the decrease of the predator density, but does not induce the extinction. These facts agree with the results in [66,69].
Moreover, the fear effect can increase the stability of the system, because it can prevent the existence of limit cycles. When the parameter m related to the prey refuge has a value on the interval [0, 1 − b(aγ + cβ)/aα(cβ − aγ)), the increase of the parameter s changes the local stability of E * from unstable to asymptotically stable, and the limit cycle vanishes, so fear makes the system evolve to the coexistence of both species. This result does not agree with the ones in [69].
Regarding the prey refuge, when m ∈ [0, 1 − 2bγ/α(cβ − aγ)), with the increase of s the predator density changes from increase to decrease. When m reaches a high risk threshold of the prey refuge, the predator goes to extinct.
We have seen how incorporating different biological issues, such as the Allee effect and the prey refuge, can change the qualitative behaviors. Now we continue combining some of these biological characteristics, in particular, the next model incorporates at the same time the Allee effect and the prey refuge. Motivated by the work of Sasmal, in [72] the authors analyze the influence of prey refuge on the model raised in [69]. They study the systemẋ where η ∈ (0, 1) is the prey refuge constant and so ηx(t) is the capacity of the refuge at time t. If we compare with the previous model included in this review, we are adding the Allee effect in addition to fear effect and prey refuge.
For system (33) the origin is always an equilibrium, particularly, a stable node. There are also two boundary equilibria, E 1 = (K, 0) and E 2 = (m, 0). E 1 can be a saddle, an attracting saddle node or locally asymptotically stable, depending on the relation between the parameter η and n/aαb. In a similar way, E 2 can be a repelling saddle-node or a saddle depending on the relation between the parameter η and n/aαm.
Under certain conditions there exist also a positive equilibrium E * . As they consider the prey refuge as a parameter, they arrive to a threshold condition for the stability of the system. It is also proved that the system has a supercritical Hopf bifurcation.
It is concluded that the prey refuge has an important role on the dynamics of the system. Comparing the results with the ones obtained in [69], the dynamic is more complex in this case. The influence of the prey refugee may not be positive for the population, as shown in Figure 9. For a fixed set of parameters, in system (27) without prey refugee the coexistence of both species is possible, but in system (33) it is not possible. The fear effect and Allee effect do not affect the final density of prey on the coexistence equilibrium, but they can cause a decrease in the predator population.

Cannibalism
Cannibalism is the act of killing and at least partial consumption of conspecifics. This actively occurs in more than 1300 species in the nature [12]. In the literature, cannibalism has been principally considered in the predator species, and in general the results agree that this cannibalism stabilizes the system and causes persistence in a population doomed to go extinct. The first contribution in this sense was the work of Kohlmeir [14], who considered cannibalism in the predator in the classic Rosenzweig-McArthur model. As an example, we present the work of H. Deng, F. Chen, Z. Zhu and Z. Li [80] and after we present the first work that considers the effect of prey cannibalism in a two species unstructured ODE model [81].
In [80], the authors consider a modification of the original Lotka-Volterra model which includes the cannibalism in the predator species. The system they propose is the followinġ where c 1 < c and cy 2 /(y + d) denotes the cannibalism of the predator.
The authors obtain the equilibria of the system and study their local and global stability. The results are compared with the ones obtained for the original system without cannibalism. There are two equilibria that always exist, the origin E 0 and the boundary equilibrium E 1 = (b/α, 0), and another boundary equilibrium This positive equilibrium E * is globally asymptotically stable whenever it exists.
The conclusion is that the presence of cannibalism modifies the stability of the system, but not always in the same direction. When in the original system the boundary equilibrium in which there are only prey is globally asymptotically stable, the introduction of cannibalism (in a suitable quantity) benefits the system, i.e., it allows the two species to coexist without becoming extinct, since a globally asymptotically stable positive equilibrium appears. If the parameter related to cannibalism increases too much, this positive equilibrium disappears and another boundary equilibrium appears, with the result that prey becomes extinct but even so, the predators manage to survive on cannibalism.
On the other hand, when the original Lotka-Volterra system is in conditions in which the coexistence equilibrium exists and is globally asymptotically stable, as the parameter associated with cannibalism increases the population of prey increases and the population of predators decreases. In this way, when a certain value of this parameter is reached, cannibalism causes the extinction of the predators.
In [81], A. Basheer, E. Quansah, S. Bhowmick and R. D. Parshad consider the the first ODE model with two species and prey cannibalism. They compare the simple Holling-Tanner model with ratio-dependent functional response in prey [82], i.e., the systeṁ with the systems obtained by adding cannibalism in the predator or in the prey species. For the original system, if β < 1/(1 − α) then there exists a positive equilibrium (x * , y * ) = (1 − β/(1 + αβ), β(1 − β/(1 + αβ))), and it is stable whenever it exists. By considering the inclusion of predator cannibalism they obtain the systeṁ Here the food source of predators is represented by γx + cy, where cy is the cannibalism term. There is a gain of energy to the cannibalistic predators from the act of cannibalism. This results in an increase in reproduction which is represented by adding a β 1 y term in the predation equation, assuming β 1 > β. For system (36) there is a positive According to the results obtained in [80], for this equilibrium the following result is proved: For a parameter set such that the positive equilibrium of system (35) is unstable, if γ < β/β 1 then there exists a value of the cannibalism rate c for which the positive equilibrium of system (36) is stable. This implies that in that case the equilibrium can be stabilized via predator cannibalism. This does not occur in the case γ > β/β 1 . Now, with the inclusion of prey cannibalism instead of predator cannibalism, the authors study the systemẋ The cannibalism term C(x) = cx(x/(x + d))) is added in the prey equation. The gain of energy to the cannibalist prey and the consequent increase in reproduction is modeled by adding a c 1 x term to the prey equation. It is assumed that c 1 < c as it takes depredation of several prey by the cannibal to produce one new offspring. For system (37) the positive equilibrium is with m = β/(αβ + 1) < 1 and under some of the following feasibility conditions: The authors answer two questions about the effect of prey cannibalism. First, can prey cannibalism destabilize the positive equilibrium? They prove that there always exists a cannibalism rate c such that the positive equilibrium of system (35) is stable while the positive equilibrium of the system (37) is unstable, so it has a destabilizing effect. The following question is, can prey cannibalism stabilize the positive equilibrium? In this case, it is proved that if m + d + x < 1 + c 1 , then no amount of prey cannibalism can stabilize the unstable positive equilibrium of the system (35). This corresponds with the first feasibility condition, under which the result is rigorously proven. The authors conjecture the same result for the second feasibility conditions, based on their numerical study. Additionally, by numerical analysis, the authors note that even with stochasticity in the cannibalism parameter, prey cannibalism is unable to stabilize the dynamics.
Moreover, the authors study the existence of limit cycles. In particular, they show that for system (37), if the conditions hold, then there exists at least one limit cycle. Please note that the uniqueness was not proved.
From an ecological point of view, these results show that in a predator-prey community with cyclical oscillations, introducing cannibalism in the prey would just maintain the behavior, and the population cannot be driven to a stable steady state.
Cannibalism on prey had previously been considered in multicomponent structured model via integer differential equations [83], structured discrete model [84] or three species structured ODE models [85].
Additionally, on mathematical literature about cannibalism, other types of predator-prey models with inclusion of cannibalism have been investigated, as three species ODE models incorporating stage structure, two species PDE models, discrete models or integro-differential equations models. [14,[86][87][88][89].

Models with Immigration
Most predator-prey systems in the wild are not isolated, so it is important to consider the effects of the presence of some number of immigrants. We will focus on studying the inclusion of immigration in the Rosenzweig-MacArthur model and in the Lotka-Volterra model.
In the literature one can find works of a very different nature that consider immigration. On many occasions, species of prey and predators migrate instantaneously after facing a hostile situation. However, after going through one of these situations, the individuals may encounter some barrier that causes delayed migration. In this sense, some models based on delay equations can be found, [90][91][92][93]. Some other studies consider the species distributed between different patches, so systems in R n with n > 2 are studied [94][95][96][97].
Taking into account the hypothesis that spatial interaction contributes to the regulation of predator-prey populations, in [98] the authors constructed a model for spatially interacting populations, adding constant immigration of prey to the Rosenzweig-MacArthur predator-prey model with a Holling type II functional response. They consider the systeṁ where the prey immigration is given by c ≥ 0. In the case c = 0, i.e., for the Rosenzweig-MacArthur model, the dynamics are well known. The origin is always an equilibrium and it is the point (K, 0). There exists a positive stable equilibrium when µ > d and a(µ + d)/(µ − d) > K > ad/(µ − d), which bifurcates at K = a(µ + d)/(µ − d). When K > a(µ + d)/(µ − d) a limit cycle emerges. When the carrying capacity K increases, the cycle brings one of both populations closer to zero.
Although E * exists, the system does not always have limit cycles, in fact there are no limit cycles if µ > d and rλ 2 /(rλ + c) < K ≤ λ. There exists a unique limit cycle if µ > d, When E * exists, it is globally asymptotically stable if and only if condition (39) does not hold. This allows the authors to state some conclusions from the point of view of ecology. They assume the hypothesis µ > d, which they consider reasonable for the maintenance of the predator population. If K ≤ λ the model without immigration has no positive equilibrium, and the boundary equilibrium is asymptotically stable, so the predators become extinct. By incorporating immigration into the model and increasing the value of c, a positive equilibrium appears when K = rα 2 /(rλ + c). This equilibrium is asymptotically stable, so authors conclude that immigration allows the coexistence.
If K > λ there exist a positive equilibrium (for any value of c). Depending on the values of the parameters, the system can present a globally asymptotically stable equilibrium, in which case the two population coexist, or it can have a limit cycle, in which case the two populations have an oscillatory behavior.
For system (38) the authors give necessary and sufficient conditions for the asymptotic stability of the positive equilibrium and the uniqueness of limit cycles. They use the Liénard's theory, after transforming the system into one of Liénard type by means of a standard way, which is detailed in the proof of the result. Their results clarify that the prey constant immigration dampens the large amplitude of oscillations emerging around the positive equilibrium. This agrees with the general idea that immigration stabilizes predator-prey populations.
In [99] the authors consider the modified system with little predator and prey immigration. They analyze the asymptotic stability of the predator-prey system but adding an immigration factor c 1 (x) into the prey population and c 2 (x) into the predator population, so they study the systemẋ where the immigration functions c i can be modeled into two ways, as c i (z) = c i or as c i (z) = c i /z for z = 0. The first expression represents a constant number of immigrants that arrive into a population, which can occur if the habitat is attractive for its quality. If we choose the second expression the number of individuals arriving depends on the population, for high populations densities the number of immigrants is lower. We recall that the case with c i < 0 can be also considered. Three different models are considered, which can be written in general aṡ where h and α are the functional response coefficients. Please note that for h = 0 and α = 0 we have a type I functional response, for h = 0 and α = 0 we have type II functional response and for h = 0 and α > 0 a type III functional response.
For each one of these three models, four different types on small immigration are considered, determined by the definition of c 1 and c 2 . Respectively they consider: c 1 (x) = c 1 and c 2 (y) = 0; c 1 (x) = c 1 /x and c 2 (y) = 0; c 1 (x) = 0 and c 2 (y) = c 2 ; and at last, c 1 (x) = 0 and c 2 (y) = c 2 /y.
For the model with type I functional response, i.e., the linear model, the authors find that without immigrants the positive equilibrium is unstable and there exist limit cycles. In this case, the coexistence equilibrium becomes stable with any of the four immigration types.
For the model with type II functional response, the positive equilibrium is unstable without immigrants, but there are no limit cycles. In this case, the coexistence equilibrium becomes locally asymptotically stable with any of the four immigration types.
Finally, for the model with type III functional response, the positive equilibrium is locally asymptotically stable if rb − hnr > 0. This equilibrium remains locally asymptotically stable regardless of the type of immigration considered. This work concludes that very small immigration into either prey or predator population acts as a stabilizing factor to the Lotka-Volterra system. Adding positive immigration will average out all fluctuations in both the populations of prey and predator. Note also that a positive immigration factor is enough to change the qualitative dynamics of the population. This paper may imply that cyclic populations can be stabilized by adding a few immigrations into them. As an example, for a fixed set of parameters, in Figure 10 we compare the behavior of system (40) with and without constant immigration. In most natural populations there are at least a few immigrants over time, and these small number of immigrating individuals are sufficient for obtaining asymptotic stability in the predator-prey system. Please note that the case with immigration in both species is not studied in [99], but we have carried out some simulations and the results obtained are similar to the obtained for the cases with immigration in only one species, as can be seen in Figure 10d.

Conclusions
Population models, and in particular predator-prey models, are a widely studied topic in biomathematics. We have reviewed some recent works focusing on the modelization of Allee effect, fear effect, cannibalism and immigration. Some of the stated models show a very rich dynamics which is much more realistic that the one of the first historical models.
The inclusion of Allee effect in a predator-prey system can have important effects on the dynamics; it can stabilize or destabilize the system or cause that the solutions of the system take a much longer time to reach a stable equilibrium point as seen in [39,48]. We note that the Allee effect can produce different consequences depending on whether it is considered in the prey species or in the predator species, see the works of Merdan [48] and Guan, Liu and Xie [49] and the comparative Figures 1 and 2. Additionally, limit cycles can appear in systems with Allee effect, as analyzed in the work of Wang, Xi and Wei [40], and even two limit cycles can appear as proved in [58].
The conclusions about the influence of fear effect are not the same in all the studied models. Considering a system with linear functional response and specialist predators, in [66] the authors show that fear does not affect the dynamics, but considering omnivorous predators, in [67] it is shown that the increase of the fear constant decreases the population density of the species in the singular point and can even cause extinction. Different behaviors for these two systems are compared in Figures 6 and 7. Additionally, for the system considered in [66], it is proved that if we change the functional response from lineal to Holling type II, then the fear affects the dynamics: a value large enough for the fear parameter can help to maintain the asymptotic stability of a singular point.
The stabilizing effect of fear is also showed when combined with other biological characteristics. In [68] the fear effect is combined with hunting cooperation, and it is shown that when the system has oscillating behaviors, the increase of hunting cooperation or the increase of fear effect produces the stabilization. Additionally, a stabilizing effect appears when combining the fear effect with prey refugee, as shown in [71].
In contrast, in [69] we see that the stabilizing effect does not appear when fear is combined with Allee effect. In this case, the only consequence is a reduction of the population density in the positive singular point. Similar conclusions are obtained in [70] with a modified Allee effect.
Cannibalism had been principally considered in the predator species, and in that case the results agreed that cannibalism stabilizes the system and causes persistence in a population doomed to go extinct, but we have seen that considering cannibalism in prey species can produce the opposite result, as it can destabilize the system [81].
Regarding the presence of immigration, we can observe that it has a stabilizing effect, reducing the amplitude of oscillations when they exist [98]. Even very low immigration can cause cyclic populations to be stabilized, as proved in [99]. Furthermore, this stabilization occurs whether immigration occurs on prey or predators.
We notice that the same methods and results are used in most papers for the study of the stability of the equilibria and the qualitative behavior of the system. Principally, the authors use the standard linear stability analysis, they analyze the sign of the real part of the eigenvalues of the Jacobian matrix of the system, sometimes using the Routh-Hurwitz criterion.
It is important to remember that for non-hyperbolic equilibria, the computation of the eigenvalues does not allow us to conclude on the stability. We have observed that in such cases the equilibria were not usually studied in the papers reviewed. Some other frequently used tools are the well-known Poincaré-Bendixon theorem, the Bendixon-Dulac criterion [57] and the computation of Lyapunov coefficients [100].
In most works, the authors make variable changes to obtain systems that are easier to study, for example by writing them down in a dimensionless form. In many cases, the systems are transformed into systems of Kolmogorov type, which are a generalization of Lotka-Volterra systems (1) to an arbitrary degree. These systems are often used in biomathematics and their theoretical study is currently progressing [101][102][103][104][105][106][107][108]. The transformation into Liénard type systems is also employed in some cases.
There are some open problems and some characteristics that would be interesting to investigate more in detail. In the following we list some problems which researchers may find interesting to address in the future, which are proposed or motivated by the works included in this review.

•
Study new functions that can model the Allee effect, i.e., functions representing the lowest population growth when population density is very low. This suggestion of considering different expression to model the same effect is also raised for some of the authors working on this topic, as for example in [48]. We recall that according to biological facts, a function α(x) representing the Allee effect must verify that α (x) is positive for any positive x because the Allee effect decreases as density increases and lim x→∞ α(x) = 1 because the Allee effect vanishes at high densities. • Among the systems that consider the Allee effect, the one proposed in [63] is the one that shows a richer dynamics. As the authors say, it is possible that by changing the functional response and taking one of Holling type III, the behavior is even more realistic, so it would be interesting to study it. • It would be interesting to consider models in which fear affects intraspecific competition, because as stated in [66], there are arguments in support of this. It would be also important to acquire experimental evidence of this phenomenon to model it as realistic as possible. • In [81] it is proved that cannibalism can lead to limit cycle dynamics, but it would be interesting to determine the uniqueness or non-uniqueness of limit cycles for cannibalistic populations, and answer the question of how many limit cycles can appear. • The effect of cannibalism on predator and prey simultaneously has not been studied as far as we know. However, what happens when both are considered, and predator cannibalism has a stabilizing effect while prey cannibalism has a destabilizing effect? A problem of this type is proposed in [81] for Holling-Tanner models, but we think this question can be studied in any general model including cannibalism. • In [99] system (41) was studied when there exists small immigration in prey or in predator, but the case with immigration in both species was not studied. It would be interesting to study the effect of immigration in both species simultaneously. • From a theoretical point of view, it would be especially interesting to carry out a study of the singular points of the systems for the values of the parameters for which they are non-hyperbolic, as in most articles these cases are omitted. • In [63] the existence of Hopf bifurcations and Bogdanov-Takens bifurcations is proved numerically. It would be interesting to find an analytical proof of these results.
• There are lots of effects and characteristics of the species and the environments that would be interesting to study, as for example, the consequences of anomalous diffusion and heterogeneous environments, see [109].
Author Contributions: All the authors have contributed in a similar way to this paper. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors are partially supported by the Ministerio de Ciencia e Innovación, Agencia Estatal de Investigación (Spain), grant PID2020-115155GB-I00 and the Consellería de Educación, Universidade e Formación Profesional (Xunta de Galicia), grant ED431C 2019/10 with FEDER funds. The first author is also supported by the Ministerio de Educacion, Cultura y Deporte de España, contract FPU17/02125.