Weak Allee effect in a predator-prey system involving distributed delays

Weak Allee effect in a predator-prey system involving distributed delays

Paulo C.C. Tabares^I,*; Jocirei D. FerreiraII, ^**; V. Sree Hari Rao^III,**,***

^IUniversidad del Quindío, Facultad de Ciencias Básicas y Tecnologías Armenia, Quindío - Colombia. E-mail: [e-mail]

^IIUniversidade Federal de Mato Grosso, Instituto de Ciências Exatas e da Terra, Pontal do Araguaia - MT - Brasil. E-mail: [e-mail]

^IIIOn leave from Jawaharlal Nehru Technological University, Hyderabad-500 085, India. E-mail: [e-mail]

ABSTRACT

In this paper we study the influence of weak Allee effect in a predator-prey system model. This effect is included in the prey equation and we determine conditions for the occurrence of Hopf bifurcation. The stability properties of the system and the biological issues of the memory and Allee models on the coexistence of the two species are studied. Finally we present some simulations which allow one to verify the analytical conclusions obtained.

Mathematical subject classification: Primary: 34C25; Secondary: 92B05.

Key words: weak Allee effect, population dynamics, Hopf bifurcation, predator-prey model, distributed delay.

1 Introduction

It is well known that in nature some species often co-operate amongst themselves in their search for food or when they try to escape from predators. Allee [15, 16], has studied extensively the aspects of aggregation and associatedco-operative and social characteristics among the members of a species. Inpopulation biology Allee effect refers to a population that has a maximalper capita growth rate at intermediate density. This occurs when the per capita growth rate increases with the increase in the density and decreases when the density passes through a critical value. Clearly this situation is different from the logistic growth in which the per capita growth rate is a decreasing function of the density. Even more, when considering coexistence for a singleprey species with one or more predator species, then the Allee effect zone is modified, because it depends also on the functional response (see in [6]).

The Allee effect is called weak if there exists no critical density population, below which the per capita rate becomes negative. Taking into account thecarrying capacity of the environment with respect to the prey in the per capita growth rate of the population, the weak Allee effect has been modeled in [4] by the following differential equation

in which N(t) denotes the prey population at time t, 0 < ε < 1 is the per capita growth rate of the population and K > 0 is the carrying capacity of the environment.

In this study we propose and analyze the following system of equations

as a model to describe the dynamical behavior of the predator-prey system incorporating a weak Allee effect in the prey populations N(t) and a distributeddelay in the predator populations P(t), at time t. The parameters K, ε, α, γ, β and a are positive. K represents the carrying capacity of the environment with respect to the prey, ε is the intrinsic growth rate of the prey, α is the rate of predation, γ is the predator mortality in the absence of prey, β is the conversion rate and a is the delay parameter. Since ε is small, we have considered a magnification of ε by denoting it as 4 ε in the model equations (2) and in the subsequent simulations. This magnification renders clarity to the simulations carried in this paper.

It is well known in biology that group defence helps decrease (or even prevent) the predation due to the enrichment in the ability of the prey to defend or escape from the predators. In view of the considerations in Freedman and Wolkowicz [11], the model equations (2) may also be viewed as a predator-prey system with group defence exhibited by the prey.

The models considered in earlier studies ([1, 2, 3, 5, 7, 8, 10, 12, 14]) though similar, but our model differs from those in terms of incorporating weak Allee effect as opposed to a logistic functional response in the prey dynamics. Though our study in this paper allows one to replace the density function ae^-at by a more general function G : [0, ∞) → [0, ∞), that satisfies G(t) > 0 and G(t)dt = 1, we prefer to retain the second equation in the system (2) in its present form, as no new ideas would be introduced by this replacement.

The present paper is organized as follows: In Section 2, we determine theequilibria and present the local stability of the equilibrium points that do notdepend on the parameters of the system. The conditions for occurrence ofHopf bifurcation for the parameter dependent equilibrium and the relatedanalysis are discussed in Section 3. The simulation results are carried usingMaple 9 and are described in Section 4. A pseudo Maple code that facilitates the simulations is presented in Section 6. A discussion follows in thefinal section.

2 Equilibria and the linear analysis

In order to determine the equilibria for the system (2), we rewrite (2) in the following form

in which

The change of variables,

transforms the system (3) into

where the prime represents the derivative respect to s, x = (n, p, q) ∈ ³ and K ∈ (0, ∞). It is easy to see that system (4) has the following equilibria

(n₁, p₁, q₁) = (0, 0, 0) , (n₂, p₂, q₂) = (1, 0, 1) and

The following result provides sufficient conditions for the local stability and equilibria (n₁, p₁, q₁) and (n₂, p₂, q₂).

Proposition 2.1. For the system (4) the equilibrium (n₁, p₁, q₁) is unstablefor all K > 0. On the other hand, the equilibrium (n₂, p₂, q₂) is locallyasymptotically stable if K < and unstable if K > .

Proof. We establish this result by considering the Jacobian matrix of system (4), which may be written as

The Jacobian matrix J(n₁, p₁, q₁) has eigenvalues 0, - and - and byreduction procedure to center manifold ([13]), we obtain = 4n² and from this it follows that the equilibrium (n₁, p₁, q₁) is unstable for all K > 0.

Similarly the eigenvalues associated with J(n₂, p₂, q₂) are given by -1, - and . Thus (n₂, p₂, q₂) is locally asymptotically stable provided K < and unstable when K > .

Remark 2.2. The instability of the origin (n₁, p₁, q₁) may be interpreted asthe non vanishing of the species simultaneously in increasing time. On theother hand the instability of the equilibrium (n₂, p₂, q₂) implies that the carrying capacity of the environment K and the conversion rate β of the predators do not support the predator population which vanishes in time. This conclusion is the same as observed in ([9]). The equilibrium is biologically meaningful only if

The condition (6) means that the dependence on time of the predator growth rate is positive, at least when assumes the value K (see [9]). Moreover, 0 < < 1 implies 0 < < β, provided that the mortality rate γ of the predator species with respect to the maximum capacity of prey K, must be less than the conversion rate β. However, for large enough K, it is always trueand necessarily (for small populations of prey), this does not ensure the survival of the predator species.

We shall discuss the stability of this equilibrium in the next section.

3 Stability analysis of the equilibrium (n₃, p₃, q₃)

In this section we study the stability properties of the equilibrium (n₃, p₃, q₃) = . We study the local codimension one Hopf bifurcation which occurs in the system (4) and also determine the direction of the bifurcation.

The Jacobian matrix in (n₃, p₃, q₃) is given by

and the characteristic equation associated with J(n₃, p₃, q₃) is

The following proposition is a direct consequence of the Routh-Hurwitz stability criterion (see [8]).

Proposition 3.1. The critical point (n₃, p₃, q₃) is locally asymptotically stable for the system (4) if the following inequalities are satisfied

Remark 3.2. The condition (8) coincide with (6). The same is observed in [8]. From (8) and (9), it follows that

Clearly this inequality highlights the contribution of the weak Allee effect, i.e. the refinement in the relationship between β, γ and K.

A rearrangement of terms in (10) yields

and this inequality is trivially satisfied if the right hand side of (11) is either zero or negative, in which case the equilibrium (n₃, p₃, q₃) is locally asymptotically stable for all a > 0. The case of interest for us would be when

In this case, (n₃, p₃, q₃) loses its stability at a = a₀ (see Fig. 1). We study the local codimension one Hopf bifurcation which occurs in the system (4) for the state variables (n, p, q) ∈ ³ and the parameter a ∈ ₊.

To simplify the calculations we introduce the notation b = and rewrite conditions (8)-(10) respectively as

At a = a₀ the eigenvalues of system (4) are given by

in which

We have the following result

Theorem 3.3. With the notation and conditions (13)-(17) the real part of the eigenvalue λ₁(a) that is Re( λ(a)) satisfies the conditions that

Proof. Note that λ₁(a₀) = i ω and define the function

obtained from the right hand side of (7) in which b = .

Since F( λ₁(a₀), a₀) = F(i ω, a₀) = 0 and i ω is a simple root of F( λ, a₀); by the implicit function theorem exist a unique one λ = λ(a) in a neighborhood of a₀ such that

where

and

Thus

and

Since

substituting

in Re we have

where

G( ε,·) = 256 ε²γ⁴b⁶ - 512 ε²γ³b⁵ + 16 ε γ³b⁴ + 384 ε²γ²b⁴ - 32 ε γ²b³

- γ²b² - 128 ε²γb³ + 20 ε γb² + 2 γb + 16 ε²b² - 4 εb - 1.

From (14) the numerator in the right hand side of (19) is positive. To complete the proof we need to show that G( ε,·) ≠ 0. To this end, consider the equation G() = 0, and observe that

are its roots. Since G^"( ε,·) = 32b²(2 γb - 1)⁴ > 0, G has a minimum in the interval ( ε₁, ε₂). Rewriting (14) as

we have

So, G( ε,·) < 0 which implies Re < 0, and this proves the Theorem.

Now we introduce the new bifurcation parameter µ by

(see Fig. 1). In this case, the equilibrium (n₃, p₃, q₃) is locally asymptotically stable for µ < 0 and looses its stability at µ = 0, that is, µ = 0 is the new bifurcation point. Further,

3.1 Projection to the center manifold

At beginning of this section we show that the eigenvalues curve of system (4) satisfies the transversality conditions of the Hopf's Theorem. In the following, we will restrict system (4) to the center manifold, that is, we will project the vector field associated to system (4) to the center manifold. So, the ODE system given by (4) will behave as a system on plane.

To reach this purpose, we consider the change of variables b = toconvert the equilibrium (n₃, p₃, q₃) to the form . Furthermore, we move this equilibrium to origin introducing the variables

Thus, with the new variables we carried out system (4) to the equivalent system

At µ = 0 (i.e. at a = a₀) we have

The eigenvectors of the above system associated with the eigenvalues λ₀(a₀) = and λ_1,2(a₀) = ±i ω are

respectively.

We take s₀ and the real and imaginary parts of s_1,2 as a new basis for ³. So, with the change of coordinate

where

we transform system (22) to an equivalent system

for x = (x₁, x₂, x₃) and

with

Now we project the system (23) to the center manifold M which is characterized by the fact that it is tangent to the x₁x₂ plane at the origin (the eigenspace corresponding to λ_1,2 = ±i ω) and is locally invariant with respect to the flow of system (23). Note that this center manifold can be parameterized by the variables x₁ and x₂ in the form x₃ = h(x₁, x₂) where h : ²→ admits a Taylor expansion of the form

and, the Center Manifold Theorem implies that h must satisfy

where h ∈ C^k, k ∈ with k > 3.

If (x₁(s), x₂(s), x₃(s)) is a solution of system (23) near the origin withinitial conditions starting on M then the flow of system (23) remains locally on M for all time, that is x₃(s) ≡ h(x₁(s), x₂(s)). Consequently

Using system (23) the above identity turns out to be

since

Omitting the terms of third and higher order we obtain

A rearrangement of terms yields

Equating the coefficient to zero and solving the resulting system we have

where

F = 3 γ⁴b³ - 6 γ³b² + 3 γ²b + εD.

Thus,

where

Finally, to get to the restriction of system (23) to the µ = 0 section of the center manifold M, we consider the new change of coordinates

y₁ = x₁, y₂ = x₂, y₃ = x₃ - h(x₁, x₂) .

In this new coordinate system, the equation of the µ = 0 section is y₃ = 0 , that is the µ = 0 section is the y₁, y₂ plane. By this coordinate transformation, system (23) becomes

Restricted to the center manifold M we can write y₃ = 0 everywhere and taking into that ₃ = 0. Hence, the system (23) restricted to the center manifold is given by

The calculations may be summarized in the next proposition, which provides the direction of the Hopf bifurcation.

Proposition 3.4. Consider the one-parameter family of differential equations (4). The first Lyapunov coefficient associated with the equilibrium (n₃, p₃, q₃) is given by

Proof. Consider the two dimensional ODE system

where

f₂(x, y) = a₁x² + a₂xy + a₃y²,

f₃(x, y) = b₁x³ + b₂x²y + b₃xy² + b₄y⁴,

g₂(x, y) = c₁x² + c₂xy + c₃y²

g₃(x, y) = d₁x³ + d₂x²y + d₃xy² + d₄y⁴.

In accordance with Bautin's Lemma (see [8]) the first Poincaré-Lyapunov coefficient for system (26) is given by

Comparing (26) with (24) and using (27), we obtain G₄ as given by (25).

□

With the results obtained in this section we can enunciate the following Theorem.

Theorem 3.5. If (12) holds and G₄ < 0 (respectively G₄ > 0), then there exists a δ > 0 such that for each a ∈ (a₀ - δ, a₀) (respectively a ∈ (a₀, a₀ + δ)), the system (4) (or equivalently system (3)) has a unique orbitally asymptotically stable (respectively unstable) periodic orbit around the equilibrium (n₃, p₃, q₃).

4 Simulation results

In this section we present some simulations carried in Maple 9 to verify the veracity of the analytical results obtained for system (4).

The equilibrium point (0,0,0) is always unstable; to illustrate this case we made the phase portrait in the case by choosing ε = 0.5, K = 1.4, α = 0.5, γ = 0.7, β = 0.8 y a = 0.25 (see Fig. 2).

The equilibrium point (1,0,1) is stable if βK - γ < 0 and unstable if βK - γ > 0. So, to see the stability of this equilibrium we take ε = 0.5, K = 1.4, α = 0.5, γ = 0.7, β = 0.3 and a = 0.25 (see Fig. 3-(a)); on the other hand, the instability is obtained when ε = 0.5, K = 1.4 , α = 0.5, γ = 0.7, β = 0.8 and a = 0.25 (see Fig. 3-(b)).

Finally, we illustrate the stability properties of the equilibrium (n₃, p₃, q₃).If we consider the following values of the parameters ε = 0.2, K = 1.4, α = 0.5, γ = 0.7, β = 0.8 and a = 0.985, the equilibrium (n₃, p₃, q₃) is asymptotically stable (see Fig. 4-(a)). Note that these values yield (n₃, p₃, q₃) = (0.625, 0.2678571428, 0.625) and the conditions (8), (9) and (10) are equivalent to 0.42 > 0, -0.28 < 0 and a = 0.985 > a₀ = 0.925 respectively. Furthermore, decreasing the value of a up to a₀ we see in Figure 4-(b) that the equilibrium (n₃, p₃, q₃) is a weak stable focus.

The parameters given in the above paragraph we obtain Figures 5 (a) and (b) illustrating the occurrence of a periodic orbit around the equilibrium (n₃, p₃, q₃) at a = 0.825 and a = 0.414628 respectively. Note that in this case we have G₄≈ -2.595557406 which shows numerically that the periodic orbit generated by the Hopf bifurcation is supercritical and a₀ = 0.925 is the bifurcation point. Finally, in Figure 6 we give a phase portrait for system (24),the restriction of system (4) to the center manifold.

5 Maple Code

In this section we include the main Maple code used to simulate Figure 2. The same code is appropriately modified for simulating other figures

1 > restart:

2 > with(plots):with(DEtools):with(plottools):with(RandomTools):

3 > unprotect(gamma):epsilon:=0.5:K:=1.4:alpha:=0.5:gamma:=0.7:

4 beta:=0.8:a:=0.25:

5 > sist:={diff(x(t),t)=4*x(t)^2*(1-x(t))-alpha*K/epsilon

6 *x(t)*y(t),diff(y(t),t)=-gamma/epsilon*y(t)+1/(epsilon*b)

7 *y(t)*z(t),diff(z(t),t)=a/epsilon*(x(t)-z(t))}:

8 > amin:=0.05:amax:=0.4:

9 > cond:=[[x(0)=Generate(float(range=amin..amax,digits=3)),

10 y(0)=Generate(float(range=amin..amax,digits=3)),

11 z(0)=Generate(float(range=amin..amax,digits=3))]]:

12 > for i from 1 by 1 to 12 do

13 cond:=[op(cond),[x(0)=Generate(float(range=amin..amax,

14 digits=3)),y(0)=Generate(float(range=amin..amax,

15 digits=3)),z(0)=Generate(float(range=amin..amax,

16 digits=3))]]:end do:

17 > xmax:=0.6:ymax:=0.6:zmax:=0.6:

18 > opc:=x=-0.1..xmax+0.07,y=-0.1..ymax+0.07,z=0..zmax+0.07,

19 stepsize=.05,axes=box,linecolor=black,color=grey:

20 > tiempo:=t=-15..15:x1:=0:y1:=0:z1:=0:

21 > espfase:=DEplot3d(sist,{x(t),y(t),z(t)},tiempo,cond,opc,

22 thickness=2,labels=[n,p,q],orientation=[145,55]):

23 > display(espfase,tickmarks=[2,2,2],axes=boxed);

Discussion

In this paper we have studied a Lotka-Volterra predator-prey system introducing the weak Allee effect in the prey dynamics. A distributed delay is considered for the predator population. Aggregation and cooperation are commonamong various species in the nature. These characteristics are observedespecially when the species search for food or when they try to escape from predators. The model studied is formulated along these lines. It is observed that this model exhibits interesting dynamics of the interacting populations due to the presence of the weak Allee effect different from the logistic growth models. Stability analysis of the equilibrium states is carried out. Conditions are derived for the occurrence of Hopf bifurcation of the equilibria of this model.

Farkas [8] has considered a predator-prey system similar to that of the present model with logistic type functional response and with weak Allee effect. The functional response considered in our model is different from that in [8]. Theorem 3.5 is an analogue of Theorem 7.3.1 (see [8]). It is observed thatthough the weak Allee effect drives the system to instability, the parameters associated with the memory will play a neutralizing role which is evident from the presence of periodic orbits.

Acknowledgments. The authors express their gratitude to late professor Miklós Farkas for the stimulating discussions. The authors wish to thank the anonymous reviewers for their constructive suggestions, which have led to the improvement of our earlier presentation.

Received: 01/XII/10.

Accepted: 19/IV/11.

#CAM-302/10.

• [1] A. Farkas, M. Farkas and L. Kaitar, On Hopf bifurcation in a predator-preymodel Diferential Equations: qualitative theory, pp. 283-290, Amsterdam, North-Holland (1987).
• [2] A. Farkas, M. Farkas and G. Szabo, Multiparameter bifurcation diagrams in predator-prey models with time lag Journal of Mathematical Biology, 26 (1988), 93-106, Springer-Verlag.
• [3] A. Bocsó and M. Farkas, Political and economic rationality leads to velcro bifurcation Applied Mathematics and Computation, 140 (2003), 381-389.
• [4] D. Boukal and L. Berec, Single-species models of the Allee effect: Extinction boundaries, sex ratios and mate enconunters Journal of Theoretical Biology, 218 (2002), 375-394.
• [5] J.M. Cushing, Two species competition in a periodic environment Journal Mathematical Biology, 10 (1980), 385-400.
• [6] K. Kiss and J. Tóth, n-dimensional ratio-dependent predator-prey systems with memory Differential Equations and Dynamical Systems, 17(1-2) (2009), 1111-1141.
• [7] L. Dai, Nonconstant periodic solutions in predator-prey systems with continuous time delays. Mathematical Bioscience, 53 (1981), 149-157.
• [8] M. Farkas, Periodic Motions Springer-Verlag, New York (1994).
• [9] M. Farkas, Dynamical Models in Biology Academic Press, New York (2001).
• [10] N. MacDonald, Time delay in a prey-predatormodels. II. Bifurcation theoryMathematical Bioscience, 33 (1977), 227-234.
• [11] H.I. Freedman and G.S.K. Wolkowicz, Predator-Prey Systems with Group Defence: The Paradox of Enrichment Revisited Bulletin of Mathematical Biology, 48(516) (1986), 493-508. Printed in Great Britain.
• [12] N. MacDonald, Time lags in biological models Lecture Notes in Biomathematics 27, Springer-Verlag, Berlin (1978).
• [13] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos Springer-Verlag, Heidelberg (1990).
• [14] V. Volterra, Lecons sur la théorie mathématique de la lutte pour la vie Gauthier-Villars, Paris (1931).
• [15] W.C. Allee, Animal aggregations Quart. Rev. Biol., 2 (1931), 367-398.
• [16] W.C. Allee, Animal aggregations: A study in General Sociology Chicago Univ. Press, Chicago (1993).

*

Research supported by Master program of Biomathematics, UNIQUINDIO, Colombia

**

Research supported by FAPEMAT - Fundação de Amparo à Pesquisa do Estado de Mato Grosso, Grant Number 462997/2009.

***

Research supported by a grant from the Foundation for Scientific Research and Technological Innovation - A division of Sri Vadrevu Seshagiri Rao Memorial Charitable Trust, Hyderabad India.

Publication Dates

• Publication in this collection
  06 Jan 2012
• Date of issue
  2011

History

• Received
  01 Dec 2010
• Accepted
  19 Apr 2011