Qualitative effects of introducing nonlinear birth and death rates for the predator in a predator-prey type model

In this paper, we study how introducing nonlinear birth and death rates for the predator might affect the qualitative behavior of a mathematical model, describing predator-prey systems. We base our investigations on a known model, exhibiting anti-predator behavior. We propose a generalization of the latter by introducing generic birth and death rates for the predator and study the dynamics of the resulting system. We establish existence and uniqueness of positive model solutions, their uniform boundedness, existence, local stability and bifurcations of equilibrium points as well as global stability properties of the solutions. Most of the solution properties are demonstrated numerically and graphically by various numerical examples. Based on the obtained results, we show that the model with nonlinear birth and death rates can describe a much more complex behavior of the predator-prey system than the classical model (i.e., with linear rates) does. Keywords-predator–prey model; Holling type II; generic birth and death rates for the predator; stability analysis; numerical simulation.


I. INTRODUCTION
For many decades, the study of predator-prey type mathematical models has been a main topic in Biomathematics (see, e.g., [2], [14] and the references therein).A classical setting for those kind of models is given by the Gause-type models (see [13] and the references therein): Here, N (t) and P (t) are the densities of prey and predator populations at time t; K is the maximum prey carrying capacity of the environment in the absence of predators; χ and d are positive constants, χ means the ratio between the consumed food and the birth rate of the predator, and d is the death rate of the predator.
For the first equation it is assumed that the prey population grows logistically in the absence of predators and the consumption of preys by the predator is described by a so-called response function F (N, P ).In the literature, this function assumes many different forms (see, e.g., [2], [8]).
For the second equation it holds that the per capita birth rate of the predator is directly proportional to the consumption and the per capita death rate is constant.
In [19], however, it is argued that the assumptions underlying the second equation in (1), are not biologically plausible.The reasoning behind this statement is as follows.
It has been observed for various real world predator-prey systems that for sufficiently small values of the predator's functional response, the predator reproduction will be zero rather than linearly increasing w.r.t. the predator's functional response.Moreover, during times of low prey density, a cessation of predator's breeding may occur.On the other hand, the reproduction rate of predators can achieve a plateau level if their prey-consumption rate becomes sufficiently high.There will always be a limit to the rate at which an individual predator can reproduce.
The predator death rate should also depend on the predator's functional response.The predator needs to consume prey at some minimal rate to avoid death by hunger.If the predator has a sufficiently high prey-consumption rate, then further increases of this rate may have little impact on its short-term chance of death.
Based on these observations, A. J. Terry [19] proposes a generalization of the classical models (1) by introducing generic forms of the birth and death rates for the predator as functions of the consumption and assuming that F = F (N ).Let us denote them by B(F ) and D(F ), respectively.Those functions are defined in [19] with the following assumptions.
The function B(F (N )) = B(F ) is continuously differentiable in N and F , for N ≥ 0 and F ≥ 0. For F ≥ 0, 0 ≤ B(F ) ≤ cF holds for some positive constant c.Also, dB/dF ≥ 0 is fulfilled.Moreover, dB/dF > 0 either for F ∈ [F 1 , ∞) where F 1 is a nonnegative constant, or for F ∈ [F 1 , F 2 ] where F 1 and F 2 are constants with 0 ≤ F 1 < F 2 .
The function D(F (N )) = D(F ) is continuously differentiable in N and F , for N ≥ 0 and F ≥ 0. For F ≥ 0, 0 < d m ≤ D(F ) ≤ d M where d m and d M are constants.Also, for F ≥ 0, dD/dF ≤ 0.
In [19], the following specific model has been studied where F (N ) = aN/(b+N ) is the Holling type-II function.
Allee effect on predator reproduction, based on a generic form for its death rate was considered in [20].
In [12], we have studied a model of the form (2) with a Beddington-DeAngelis functional response . There, we have also slightly changed the original assumptions, made in [19], in order to exclude a few possibilities of degenerate behaviors of the system.The study has shown that introducing non-linear birth and death rates for the predator in the classical Beddington-DeAngelis model does not lead to qualitatively new profiles of the system.The typical behavior is exhibited, i.e. an internal globally stable equilibrium, a limit cycle, or a globally stable boundary equilibrium corresponding to the extinction of the predator population are possible depending on the model parameters [7], [10], [11], [12].
It is, thus, natural to consider whether introducing nonlinear birth and death rates always preserves the dynamics of a model system or, in some cases, it might enrich the dynamics and lead to qualitatively new behavior, depending on the model parameters.In the present work, we show that the latter is indeed the case.
The main goals of the present paper are: 1) to compare the possible qualitative behavior of a particular model with linear birth and death rates versus an analogous model with nonlinear rates; 2) to show, thus, that introducing nonlinear birth and death rates can lead to richer dynamics and, therefore, describe more complex behavior of an ecological system; 3) based on classical methods in the theory of nonlinear dynamical systems, to completely characterize the global dynamics of a model system with generic birth and death rates, when varying the values of the model parameters.The paper is structured as follows.In Section II, we formulate two predator-prey models-with linear and with generic birth and death rates for the predator.In Section III, we study the dynamics of the model with generic birth and death rates, i.e. general properties of the model (positiveness and boundedness, existence and uniqueness of the solutions), conditions for the existence of equilibrium points, local and global asymptotic properties of the solutions.At the end of Section III, we present numerical examples, illustrating the theoretical results in the previous section and showing how they can be easily extended.A comparison between the dynamics of the two models and discussion about the biological relevance of the differences is presented in Section IV.

II. MATHEMATICAL MODELS
We consider a model with linear (w.r.t. the functional response) birth and death rates for the predator and study the effects of generalizing them by introducing non-linear functions.It is known that models with non-monotonic functional responses (like Holling type-IV) exhibit richer dynamics (see, e.g., [5], [9], [16], [21] and the references therein).Our investigations are based on a known model exhibiting anti-predator behaviour and proposed in 2015 by Tang and Xiao [18]: where F (N ) = aN/(b + N 2 ), and the positive constant η is the rate of anti-predator behavior of prey to the predator population.
In our work, we compare the behavior of the latter model with a similar model, incorporating generic birth and death rates for the predator.Thus, for convenience, we shall briefly comment on the possible dynamical behavior of the solutions of (3).It was shown in [18] that the following possibilities exist.
• There exist two equilibria-E 0 = (0, 0) (saddle point) and E K = (K, 0) (globally stable).• There exist three equilibria-E 0 = (0, 0) and E K = (K, 0) (saddle points), and one internal equilibrium E 1 = (x 1 , y 1 ), which might be either stable, or unstable (in the latter case a limit cycle exists).• There exist four equilibria-E 0 = (0, 0) (saddle point) and E K = (K, 0) (stable equilibrium), and two internal equilibrium E 1 = (x 1 , y 1 ), which might be either stable, or unstable (in the latter case a limit cycle exists) and E 2 = (x 2 , y 2 ), which is a saddle point.There exist two basins of attraction.Depending on the initial conditions, trajectories might tend either to E K or to E 1 /to a limit cycle if E 1 is stable/unstable.For more details, including conditions for the model parameters leading to each of the cases, see [18].
We shall compare the dynamics of (3) with a model, including generic birth and death rates for the predator.In order to simplify the analysis, we shall consider here a Holling type-II functional response.Thus, the model we study is  (iii) there exist non-negative constants (θ 2 possibly equal to +∞) such that, for Remark 1. Conditions (B)(iv), (D)(iv), and (E) are included in order to make the analysis more tractable, although the results of the present paper could be easily extended if those assumptions were not valid.The main purpose of the present work, however, is to show that nonlinear birth and death rates can enrich the dynamics of a particular predator-prey model system and, thus we decide in favor of a better readability than considering the most general setting.In this way, we want to highlight the important qualitative results.Nevertheless, we believe that the behavior of the functions B and D, as defined by (B)(iv) and (D)(iv), is natural.

III. DYNAMICS OF THE MODEL WITH GENERIC BIRTH AND DEATH RATES
Now, we shall study the dynamics of the model (4) and show that it is able to describe the behavior that (3) does, but for certain values of the model parameters the phase portraits of (4) might become much more complex.

A. General properties of the model
First, we shall prove that the model (4) possesses some standard properties that we would expect from a predator-prey model, namely uniqueness and positiveness of the solutions for positive initial conditions.Proposition 1.The positive cone R 2 + := {(N, P ) : N ≥ 0, P ≥ 0} is a positively invariant set for the model (4).
Proof: From the first equation in (4) we obtain and, therefore, Analogously, from the second equation in (4), we have If N (0) = 0 it follows that N (t) = 0 for every t > 0 and if P (0) = 0, then P (t) = 0 for every t > 0. Using the uniqueness of the solutions, it follows that the coordinate semi-axes are positively invariant.Hence, it is obvious that if the initial conditions are non-negative, the solutions are also non-negative for all positive times t.
Proposition 2. The non-negative solutions of the model (4) are uniformly bounded above.
Proof: The statement follows from the dissipativeness of the system and a standard result (cf.[6, pp. 17-18]).
In order to study the dynamics of the system (4), it is sufficient to consider only initial conditions in the set {(N, P ) : 0 < N < K, 0 < P }.

B. Existence of equilibrium points
The equilibrium points of the model (4) are the solutions of the algebraic system Let us note that we are interested only in the equilibria having non-negative components.
Obviously, the boundary equilibria E 0 = (0, 0) and E K = (K, 0) always exist.If internal (i.e. with positive components) equilibria exist, they should satisfy where 0 < N < K is a solution of the equation and E (N ), as defined in (E), is a monotonically increasing function that has either one, or three inflection points (see Fig. 1).We shall now study the conditions for the existence of internal equilibrium points of the model (4).First, the following proposition is obvious.
, then there is one internal equilibrium E 1 ; (iii) if N 0 < K and ηK > E (K), then there are two internal equilibria E 1 and E 2 .
Proof: The internal equilibria of (4) correspond to the positive solutions of equation ( 6).We shall examine the graphs of the left-hand side and the right-hand side of the latter (see Fig. 2).Let l = l(N ) be a tangent to the graph of E (N ), passing through the origin.Let (N 0 , E (N 0 )) be the point of tangency.Then, N 0 is a zero of the function ϕ(N First, we will show that there exists exactly one positive zero of ϕ(N ).Let C 1 and C 2 be such points that Taking into account the above results, it follows that the tangent l(N ) = E (N 0 )N is uniquely defined.Moreover, N 0 > N inf l holds.
Then, if η is greater than the slope E (N 0 ), equation ( 6) has no positive solutions.If η < E (N 0 ) holds true, there exist two crossing points of the two graphs.Taking into account that we are only interested in those with abscissa less than K, the statement of the proposition is easy to be verified.We illustrate the cases in Fig. 3.
The case when the function E (N ) has three inflection points, can be treated in a similar way.Using the same geometrical approach it can be shown that the model ( 4) can have up to four equilibrium points.We shall discuss this case in more details on the basis of the numerical experiments in Section III-E.

C. Local stability and bifurcations of the equilibrium points
In this section, we shall study the conditions for stability of the equilibria.As we have shown in the previous section, the model (4) has two boundary equilibria-E 0 = (0, 0) and E K = (K, 0) and it may have up to four internal equilibria.Let us denote them by E i = (N i , P i ), i = 1, . . ., n, where n is the number of internal equilibria, and The variational matrix of the system (4) is Proposition 5.The equilibrium point E 0 = (0, 0) is a saddle point for all positive values of the parameters in the model (4).
Proof: Taking into account conditions (B) and (D) as well as the variational matrix evaluated at E 0 , where the positive constant D 2 is defined in (D), the proposition follows directly.
Proposition 6.The boundary equilibrium Proof: The proposition follows from the variational matrix of (4) evaluated at E K , which has the form , then E i is an unstable equilibrium (a node or a focus).
Proof: The variational matrix, evaluated at E i , has the form Thus, if η > E (N i ), the real parts of the eigenvalues have opposite signs and the equilibrium is a saddle point.If η < E (N i ), then the real parts of the eigenvalues have the same signs.If they are positive, i.e. if the trace τ Otherwise, it is locally asymptotically stable.Taking into account (5), for the trace τ (E i ) we obtain The condition τ (E i ) < 0 is equivalent to and, thus, the proposition is proved.
For further use, denote Remark 2. Proposition 7 implies that the internal equilibria with even indices are saddle points.The ones with odd indices might be locally asymptotically stable or unstable nodes or foci, depending on the sign of τ (E i ), i.e. whether N i > N c or N i < N c (in the latter case it is assumed that K > b).
Remark 3. If E (N ) has exactly one inflection point, then the following possibilities for the model (4) exist.
(a) There are no internal equilibria; E 0 is a saddle point, E K is locally asymptotically stable.(b) There is one internal equilibrium that can be either asymptotically stable, or unstable; E 0 and E K are saddle points.(c) There are two internal equilibria -E 1 is either asymptotically stable, or unstable, and E 2 is a saddle point; E 0 is a saddle point, E K is asymptotically stable.
It follows from (6) that the components N i of the equilibria E i = (N i , P i ), i = 1, 2, (when they exist) do not depend on the parameter K. Thus, we can consider K as a bifurcation parameter.
The above propositions imply that when K = N 1 or K = N 2 , transcritical bifurcations of the equilibria occur at E K leading to the appearance of E 1 or E 2 respectively and exchange of stability.
Consider the case when the model (4) possesses two internal equilibria, E i = (N i , P i ), i = 1, 2. Proposition 8. Let E 1 = (N 1 , P 1 ) be the internal equilibrium point of the model ( 4) and K = 2N 1 + b.Then a Hopf bifurcation occurs at E 1 .
Proof: For K = 2N 1 + b (or equivalently when N 1 = N c ) the trace τ (E 1 ) of the Jacobian matrix J(E 1 ) becomes equal to zero, i.e.
thus, J(E 1 ) possesses a pair of pure imaginary complex conjugate eigenvalues.Further, The latter means that there exists a constant > 0 equivalently, when N c − < N 1 < N c ) then a limit cycle appears.We shall see below (cf.Theorem 2) that the limit cycle is stable.

Proof:
The proof follows ideas from Theorem 4 in [18].
It is straightforward to see that under the above conditions the Jacobian matrix J(E b ) possesses a double zero eigenvalue: The coordinate change X 1 = N − N 0 , Y 1 = P − P 0 translates the point E b into the origin (0, 0), Taylor expansions of the right-hand side functions about (0, 0) yield In the above system, O(|(X 1 , Y 1 )| 3 ) indicates smooth functions containing terms of order at least 3 in X 1 and Y 1 .The next coordinate change translates the latter system into the following one: Obviously, the coefficient 2 is strongly positive.The sign of the coefficient in front of X 2 Y 2 depends on the sign of the term Since both coefficients A 1 and A 2 are nonzero, this finishes the proof.
The existence of the cusp bifurcation at E b = E c means that two new equilibrium points, E 1 and E 2 , are 'born' when the model parameters are varied.

D. Global behavior of the solutions 1) No internal equilibria or one internal equilibrium:
Theorem 1.If the system (4) has no internal equilibria, then the boundary equilibrium E K is a globally asymptotically stable equilibrium point of the system.Proof: Since (4) has no internal equilibria, it follows that ηN > E (N ) for every Then the following holds true: Therefore, the asymptotic behavior of the trajectories of ( 4) is determined by their behavior on the positively invariant set {(N, P ) : N > 0, P = 0}.
Thus, E K is a globally asymptotically stable equilibrium point.
Proof: In this case, all the equilibrium points are unstable, see Remark 3(b).
We shall use the Poincaré-Bendixson Theorem [6] and the Butler-McGehee Lemma [4] to obtain the desired result.First, let us note that the ω-limit set of no trajectory can consist of one point only because the system has no locally stable equilibria.
Let us suppose that E 0 = (0, 0) is in the ωlimit set of an orbit γ.Then, from the Butler-McGehee Lemma it follows that the ω-limit set has a non-empty and non-trivial intersection with the stable manifold of E 0 , i.e. with S = {(N, P ) : N = 0, P ≥ 0}.Since the ω-limit set is invariant with respect to the trajectories of the system, it follows that it contains the whole set S, but this is a contradiction with the fact that the ω-limit set is bounded [6, p. 47].Therefore, E 0 cannot be in the ω-limit set of any trajectory.Now, we shall prove the same for the equilibrium point E K = (K, 0).Similarly, if we assume that E K is in the ω-limit set of a trajectory, using the fact that the ω-limit set is closed, it should also contain one of the following two sets: The first one, however, is unbounded, and the second one contains the point E 0 , which cannot be in the ω-limit set of any trajectory.As we have already noted, the ω-limit set cannot contain only the point E 1 and, therefore, it contains no equilibria.Thus, the ω-limit set of every trajectory is a periodic orbit.Theorem 3. If the system (4) has exactly one internal equilibrium, E 1 = (N 1 , P 1 ), which is locally asymptotically stable, i.e.N 1 > N c = (K − b)/2, then E 1 is globally asymptotically stable.
Proof: We shall construct a Lyapunov function for the system (4), using ideas, proposed in [1] and [3].Let us define where θ is a parameter that is to be defined.
The derivative of V (N, P ) along the trajectories of ( 4) is We want to choose θ in such a way that V is negative for every N = N 1 .First, let N > N 1 hold true.Then, we have E (N ) − ηN > 0 and θ should satisfy The numerator of the latter is positive for all for all N ∈ (0, N 1 ).
It is easy to see that the right-hand side of the latter is non-positive for N ∈ [K −b−N 1 , N 1 ), see Biomath 6 (2017), 1703167, http://dx.doi.org/10.11145/j.biomath.2017.03.167 Fig. 4. Therefore, it is sufficient to find a positive θ, such that the following condition is satisfied: Let us introduce the notation .
It remains to be proven that max Using conditions (B) and (D), for N > N 1 we obtain Then, in order to prove (9), it is sufficient to show that max Let us assume the contrary.Then, taking into account that W (K − b − N 1 ) = 0, there exists β > 0, such that W (N ) = β has at least three distinct roots in the interval [0, K] and, thus, the function has at least four distinct roots in [0, K].From the Rolle's Theorem it follows that W (N ) has at least one root in [0, K].After some calculations, we obtain which is in a clear contradiction with the latter statement.This concludes the proof of the theorem.
2) Two internal equilibria: be the two internal equilibrium points of (4).In order to study the dynamic behavior of the system in this case, first, we shall study the structure of the nullclines of ( 4).
The non-trivial nullcline, defined by the first equation is the parabola Obviously, N c is the point, where the parabola P = P (N ) takes its maximum.The non-trivial nullclines, defined by the second equation, are Then, the nullclines divide the positive cone into six regions, denoted by I, II, . .., VI, as shown in Fig. 5.
First, we shall study the case, when the internal equilibrium E 1 = (N 1 , P 1 ) is locally asymptotically stable, i.e. when N 1 > N c , see Remark 3(c).
Then, the stable manifold of the saddle point E 2 divides the positive cone into two disjoint invariant sets Ω in and Ω out .We define Ω in to be the set, containing E 1 , see Fig. 5. Let us denote the stable manifold of E 2 with W S (E 2 ).Theorem 4. Let the equilibrium point E 1 be asymptotically stable, i.e.N 1 > N c .Then, for every trajectory γ(t), originating in Ω in it holds true that γ(t) → E 1 as t → ∞.Proof: Since no equilibria exist in region I (as denoted in Fig. 5), then all trajectories originating in Ω in enter the set Ω in := {(N, P ) ∈ Ω in : N ≤ N 2 }.Moreover, Ω in is a positively invariant set for (4) because the vector field on the boundary with region I points towards Ω in .Now, analogously to the proof of Theorem 3, it can be shown that V (N, P ) (as defined in the same proof) is a Lyapunov function for (4) in Ω in .Thus, we obtain the statement of the theorem.
Theorem 5. Let the equilibrium point E 1 be asymptotically stable, i.e.N 1 > N c .Then, for every trajectory γ(t), originating in Ω out it holds true that γ → E K as t → ∞.
Proof: Since the solutions of (4) are uniformly bounded and there are no internal equilibria in Ω out , taking into account the structure of the vector field (see Fig. 5), all trajectories, originating in Ω out , eventually enter region VI.
The vector field on the boundary of region VI points towards it (see Fig. 5), thus, region VI is a positively invariant set for the system (4).There are no internal equilibria in it and, hence, no periodic orbits.Therefore, from the Poincaré-Bendixson Theorem it follows that every trajectory originating in this set tends to the only equilibrium point E K .This concludes the proof.Now, we shall study the asymptotic behavior of the model solutions when the internal equilibrium E 1 = (N 1 , P 1 ) is unstable, i.e. when K > b and Let us denote the nullcline N = N 2 with n.Theorem 6.Let the equilibrium E 1 = (N 1 , P 1 ) be asymptotically unstable, i.e.K > b and N 1 < N c hold true.If W S (E 2 ) ∩ n = {E 2 } (see Fig. 6(b)), then the ω-limit set of every trajectory starting in Ω in is a periodic orbit.
Proof: In this case the equilibrium E 2 = (N 2 , P 2 ) is a saddle point.Denote by λ the positive eigenvalue of J(E 2 ), the eigenvector corresponding to λ is Then the slope k 1 of the tangent line to the unstable manifold of E 2 is equal to The slope m of the tangent line to the nullcline P = P (N ) from ( 5) at the point < 0, and obviously k 1 < k 2 < 0 holds true.The latter inequality ensures that the branch of the unstable manifold of E 2 with P > P 2 lies above the nullcline P = P (N ).Following some ideas from the proof of Theorem 5.1 in [17] it is straightforward to see that when the orbit of the unstable manifold of E 2 leaves E 2 , it emerges consecutively in regions II, III by crossing the nullcline N = N 1 above N 1 , then in regions IV and V (see the direction field in Figure 6b), staying inside Ω in and finally crosses the nullcline Denote by G the curve, consisting of the orbit of the unstable manifold of E 2 and the part of the nullcline P = P (N ) between N Q and N 2 , and let G in be the "inside of G".Further, G in is compact and positively invariant; moreover it contains only one equilibrium point E 1 .From the Poincaré-Bendixson theorem, the ω-limit set of every trajectory, starting at a point in G in which is different from E 1 , is a periodic orbit in G in .Finally, the dissipativeness of (4) implies that every trajectory of the system starting in Ω in enters G in .This proves the theorem.Theorem 7. If the equilibrium E 1 is asymptotically unstable (i.e.K > b and N 1 < N c ) and if W S (E 2 ) ∩ n = {E 2 } (see Fig. 6(b)), then every trajectory of (4), originating in Ω out , converges to E K .
Proof: The proof follows the same ideas as the proof of Theorem 5.
The case, when W S (E 2 ) ∩ n = {E 2 }, is graphically presented in Fig. 6(a).The stable manifold of E 2 does not enclose a positively invariant set and all the trajectories converge to E K .
3) More than two internal equilibria: It can be shown that when more than two internal equilibria exist, there are many possibilities for the dynamics of the system (4).We shall comment on those possibilities on the basis of a couple of numerical examples in the next section, since the reasoning behind the analytic study mainly copies what has been said so far.

E. Numerical examples
In this section, we shall show some numerical examples that illustrate the behavior of the solutions.For the particular expressions of B(F ) and D(F ), in the numerical examples we shall follow [19].Define where all the parameters are positive and, furthermore, A 1 < A 2 , θ 1 < θ 2 , and As discussed in [19], these functions are constructed such that they capture several important E 2 The stable manifold of E2 crosses N = N2 twice.Fig. 6.Nullclines for the model (4) in the case of two unstable internal equilibria; the thick dashed line denotes the stable manifold of the saddle point E2.ideas (Fig. 7).B(F ) is zero for all sufficiently small values of F , which represents the idea that it needs a certain level of energy intake before a predator can reproduce.Also, if the functional response is sufficiently large, the reproduction rate Similarly, the death function D(F ) represents the idea that a predator will suffer its greatest risk of mortality if its prey consumption rate is little or nothing.Also, there is some minimal death rate that is reached once the consumption rate reaches a certain threshold.
1) Zero or one internal equilibrium: In the numerical experiments here we shall use the following values of the parameters in the birth and death rate functions: In this case, the function E (N ) has one inflection point.
Example 1. First, we illustrate the case, when the model ( 4) has no internal equilibria.We choose the model parameters to be r = 3, a = 5, b = 1, K = 6.For those values of the model parameters, the point N 0 , as defined in Proposition 4 is approximately equal to 3.7662.Then E (N 0 ) ≈ 0.4516.We choose η = 1 > E (N 0 ).All the trajectories tend to the boundary equilibrium E K , see Fig. 8. Example 2. Our next example is for the case when the model (4) has one unstable internal equilibrium.Let all the model parameters, except η, are the same as in Example 1.We choose η = 0.3 < E (N 0 ).For those values we have for the trace τ (E 1 ) ≈ 0.431 > 0 and, thus, the internal equilibrium E 1 is unstable and the trajectories are periodic, see Fig. 9(a).2) Two internal equilibria: Here, we shall use again the values (10) for the model parameters; this means that the function E (N ) possesses one inflection point.
Example 4. Increasing K (with respect to Example 3) to K = 22 and leaving the other parameters unchanged, another internal equilibrium appears.In this case the internal equilibria are E 1 = (12, 4.0909) and E 2 = (16.5587,2.778).
Since N c = 9 < 12, E 1 is stable and according to Theorem 4 and Theorem 5, E 1 and E K are stable and the separatrix is the stable manifold of the saddle point E 2 , see Fig. 10(a).
Example 5. Increasing K further (with respect to Example 4) to K = 31 and K = 32, E 1 becomes unstable.In the first case K = 31, the trajectories originating in Ω in tend to a limit cycle.The ones that originate in Ω out tend to E K as t → ∞, see Fig. 10(b).
In the second case K = 32, E K becomes globally stable, see Fig. 10(c).
3) More than two internal equilibria: In the next examples, we consider the following values for the parameters in B(F ) and D(F ): and we shall use η = 1.375.In this case, as can be seen from Fig. 11, the function E (F ) has three inflection points and the model (4) can have up to four internal equilibria, depending on the value of K.
First, when all the internal equilibria with odd indices are locally asymptotically stable, they are globally stable within their basins of attraction, determined by the stable manifolds of the saddle points.In this case, all trajectories converge to a stable equilibrium point, as we illustrate in the following two examples.From Proposition 7, it can be easily seen that both E 1 and E 3 are locally asymptotically stable.Now, analogous arguments to the ones in Section III-D for the vector field, shown in Fig. 12(a), can show that the system (4) has two stable equilibria-E 1 and E 3 -and the separatrix is the stable manifold of the saddle point E 2 , see Fig.When there exists an unstable internal equilibrium, there are many possible dynamics of the system.We shall illustrate the principal possibilities in the case of four internal equilibria.If E 1 is unstable and E 3 is locally stable, there exist two possibilities.
Example 8. Increasing K with respect to the previous example and choosing K = 8.3, the equilibrium point E 1 becomes unstable and a limit cycle appears.Depending on the initial condition, the trajectories either tend to the limit cycle, or to the equilibrium point E 3 , or to the boundary Example 10.In the case of K = 9.4, the equilibrium point E 3 becomes also unstable and a limit cycle appears, enclosing it.All trajectories tend either to this limit cycle, or to the boundary equilibrium E K , depending on their initial condition, see Fig. 16(a) and Fig. 16(b).
The case of three internal equilibria can be studied similarly.When the carrying capacity K is increased, the following possibilities appear: • E 1 and E 3 are both stable; there exist two basins of attraction, corresponding to both the stable equilibria; • E 1 is unstable and E 3 is stable; a limit cycle appears, enclosing E 1 ; trajectories tend to either the limit cycle or E 3 ; • E 1 is unstable and E 3 is stable; all trajectories  tend to E 3 as t → ∞; As can be seen from the discussion so far, model (3) can have up to two internal equilibria.Comparing to the results, derived in Section III-D2  and illustrated in Sections III-E1 and III-E2, we can conclude that when (4) does not have more than two internal equilibria, it can exhibit the same qualitative behavior as the model with linear birth and death rates for the predator.Model (4), however, can also have more than two internal equilibria, which leads to qualitatively new dynamics that are neglected in the classical models, but can have serious impact on the possibilities for persistence and control of a real system.The nonlinearity in the birth and death rates leads to new "regimes" of long-term persistence of a predator-prey system.
The main difference from biological point of view, however, is in the way the predator-prey system reacts to large perturbations.Let us consider the case, presented in Example 6.Let us assume that for some period of time the population fluctuates around the equilibrium E 1 .A large perturbation at time t P (i.e.leaving the basin of attraction of E 1 ) can lead to a future stabilization of the predator-prey system around another equilibrium, E 3 , see Figure 17.Such a behavior cannot be described by the classical models, i.e. with linear birth and death rates, but has nevertheless been observed in real systems [15].Also, the rich dynamics of the system give more possibilities from a control point of view.For example, if the system has two or more stable internal equlibria, one can artificially design the necessary perturbations in the system so that it is stabilized around the best (in some sense) possible equilibrium.
Another possible scenario is the following.Let the system be characterized with periodic oscillations of the two populations, as depicted in Fig. 18.Large perturbations at time t P can stabilize the system around a stable equilibrium. .Stabilization of a predator-prey system after a large perturbation in the case of a limit cycle and a stable internal equilibrium.
Taking into account the above examples, it is evident that a principle biological implication of the Biomath 6 (2017), 1703167, http://dx.doi.org/10.11145/j.biomath.2017.03.167 enriched dynamical behavior of the system lies in the way the system reacts to large perturbations.A perturbation in a predator-prey system, described by a model with nonlinear birth and death rates for the predators, may lead to a change in the long term behavior of the system, while still preserving the persistence of both the organisms.
On the other hand, as our analytical studies so far show, generalizing a model by including generic birth and death rates is not necessarily connected with too big technical difficulties for its analysis and, thus, most of the classical models could be generalized in a similar manner.It is up to be studied what implications this will have on the dynamics the resulting models can exhibit.

V. CONCLUSION
Our investigations are based on a known predator-prey model exhibiting anti-predator behavior, proposed and studied in [18].We propose a generalization of the model by introducing generic functions for the birth and death rates of the predator as suggested in [19] and already exploited by the authors in [12].The generic birth and death rates are defined under general assumptions.Although the functional response in our model is taken to be of Holling type II (a monotone function) the new model exhibits very rich and complex dynamics.We establish existence and uniqueness of positive model solutions, their uniform boundedness, existence, local stability and bifurcations of equilibrium points as well as global stability properties of the latter.A variety of numerical examples is considered to confirm the theoretical studies and to demonstrate the many forms of coexistence between predator and prey as well as possible extinction of the predator population.
In particular, allowing for the predator's birth and death rates to be nonlinear w.r.t. the functional response, the system might have more than two internal equilibrium points, which is not valid in the original model, studied in [18].This results in the possibility of having two or more stable internal equilibria or having a stable internal equilibrium and a stable limit cycle.The basins of attraction of the respective ω-limit sets are defined by the stable manifolds of the internal saddle points.
From biological point of view, the findings of our study show that models with nonlinear birth and death rates may describe behavior of the biological system that is more complex than what the classical models could describe, but is, nevertheless, relevant to the real-world food webs.One possible example concerns the case when the model (4) has two stable internal equilibrium points.It corresponds to a behavior that has been observed in some biological systems-after a large perturbation in a stable system, the population begins to fluctuate around an equilibrium, which is different than the original one.
A natural next step in the study of those models is their validation with respect to experimental data, which should give further insight on the ways the theoretical results are applied in real systems.

Biomath 6 (Fig. 5 .
Fig. 5. Nullclines for the model (4) in the case of one stable internal equilibrium and an internal saddle point; the thick dashed line separates the positive cone into two disjoint sets-Ωin and Ωout.
The stable manifold of E2 does not cross N = N2 twice.

Fig. 7 .
Fig. 7.The generic birth and death rate functions

Fig. 13 .
Fig. 13.Dynamics of system (4) in the case of four internal equilibrium points.

8 .
E1 is unstable and E3 is locally asymptotically stable; the thick dashed line separates the basin of attraction of EK .

Fig. 14 .
Fig. 14.Dynamics of system (4) in the case of one stable and three ustable internal equilibrium points.Example 8.

Fig. 15 .
Fig. 15.Dynamics of system (4) in the case of one stable and three ustable internal equilibrium points.Example 9.

10 .
E1 and E3 are both unstable; the thick dashed line separates the basin of attraction of EK .
Example 10.Dynamics of the model in the vicinity of E1 and E3.A limit cycle exists around E3. Two trajectories depictedwith a dashed and a solid line.The black squares denote the corresponding initial conditions.

Fig. 16 .
Fig.16.Dynamics of system (4) in the case of four unstable internal equilibrium points.

Fig. 17 .
Fig.17.Stabilization of a predator-prey system after a large perturbation in the case of two stable internal equilibria.

Fig. 18
Fig.18.Stabilization of a predator-prey system after a large perturbation in the case of a limit cycle and a stable internal equilibrium.
• E 1 and E 3 are both unstable; a stable limit cycle appears around E 3 .