Global Asymptotic Stability and Hopf Bifurcation in a Homogeneous Diffusive Predator-Prey System with Holling Type II Functional Response

Abstract

In this paper, we considered a homogeneous reaction-diffusion predator-prey system with Holling type II functional response subject to Neumann boundary conditions. Some new sufficient conditions were analytically established to ensure that this system has globally asymptotically stable equilibria and Hopf bifurcation surrounding interior equilibrium. In the analysis of Hopf bifurcation, based on the phenomenon of Turing instability and well-done conditions, the system undergoes a Hopf bifurcation and an example incorporating with numerical simulations to support the existence of Hopf bifurcation is presented. We also derived a useful algorithm for determining direction of Hopf bifurcation and stability of bifurcating periodic solutions correspond to j 0 and j = 0, respectively. Finally, all these theoretical results are expected to be useful in the future study of dynamical complexity of ecological environment.

Share and Cite:

Wang, S. , Yu, H. , Dai, C. and Zhao, M. (2020) Global Asymptotic Stability and Hopf Bifurcation in a Homogeneous Diffusive Predator-Prey System with Holling Type II Functional Response. Applied Mathematics, 11, 389-406. doi: 10.4236/am.2020.115028.

1. Introduction

Since C.S. Holling proposed several kinds of functional responses (Holling functional responses) to model the phenomena of predation in 1965; the classical Lotka-Volterra predator-prey system was extended and more realistic [1] [2] [3] [4] in biomathematics. These functional responses describe how predators transform the harvested prey into the growth of itself and were discussed by numbers of researchers [5] - [10]. When the spatial distributions of the two populations are also of interest, the passive dispersal of the populations can be modeled and simulated by diffusive operators [11]. Complicated diffusive predator-prey systems in the form of partial differential equations (PDEs) with Holling type II functional response have been constructed and analyzed in several previous literatures [12] - [18].

For instance, in [12], a continuous diffusive predator-prey model incorporating Holling type II functional response of the predator and a logistic growth of the prey was shown to exhibit temporal chaos at a fixed point in space. Numerical results demonstrated that low diffusion values drive a periodic system into aperiodic behavior with sensitivity to initial conditions. [13] considered the case where densities of predator and prey are both spatially inhomogeneous in a bounded domain subject to homogeneous Neumann boundary condition, and they also studied qualitative properties of solutions to this reaction-diffusion system. They showed that even though positive constant steady state is globally asymptotically stable for the ordinary differential equation (ODE) dynamics, non-constant positive steady states can coexist in a PDE system.

With regard to Hopf bifurcation analysis, [18] carried out Hopf and steady state bifurcation, and the existence of multiple spatially non-homogeneous periodic orbits are showed in particular, while the system parameters are all spatially homogeneous. The global bifurcation theory also suggested the existence of loops of spatially non-homogeneous periodic orbits and steady state solutions. Based on this reference, [19] considered the possibility of the occurrence of Turing patterns and performed detailed Hopf bifurcation analysis in a diffusive predator-prey system with Holling type III functional response. They showed that the system has multiple oscillatory patterns.

Motivated by the reference [18], in this paper, we mainly consider a homogeneous reaction-diffusion predator-prey system with Holling type II functional response with density-dependent predator specific death rate and predator mutual interference:

u t d 1 Δ u = r 1 u ( 1 u K 1 ) α u v a + u m 1 u , t > 0 , x Ω , (1a)

v t d 2 Δ v = α e u v a + u m 2 v d v 2 , t > 0 , x Ω , (1b)

ν u = ν v = 0 , t 0 , x ( Ω ) , (1c)

u ( 0 , x ) = u 0 ( x ) 0 , v ( 0 , x ) = v 0 ( x ) 0 , x Ω . (1d)

Here functions u = u ( t , x ) 0 and v = v ( t , x ) 0 are prey and predator densities, respectively. The one dimensional spatial domain Ω is ( 0 , l π ) ( l > 0 ) . All above positive constants have practically biological considerations. Parameter r 1 is the intrinsic growth rate of prey; K 1 represents the carrying capacity of environment; a is the half-saturation constant; α is the search efficiency of predator for prey; m 1 and m 2 are the mortality rate of prey and predator species, respectively; e is the biomass conversion; d is the intra-specific competition coefficient of predator; d 1 and d 2 are two positive diffusive rates of prey and

predator, respectively. The specific growth term r 1 x ( 1 x K 1 ) governs the increase of prey in the lack of predator. The coupled term x y a + x , named Holling

type II functional response, describes the functional response of predator, which also refers to the change in the density of prey attached per unit time per predator as prey density changes. The square term d y 2 denotes intrinsic decrease and mutual interference of predator. In the absence of diffusion, its corresponding ODEs system is familiar to the Lotka-Volterra system in which populations have the addition of damping terms(or self inhabit) [9]. To describe an environment surrounded by dispersal barriers, we take zero flux at such boundary ( Ω ) [12]. The symbol ν is the outer flux, and no flux boundary condition is imposed, thus the system is closed [18] and we have above Neumann boundary conditions (1c). Bazykin [20] once looked at the ODEs version of above system, and it is also researched by some authors ( [9] and [21], etc.).

For simplicity and convenience in the later, we introduce a new dimensionless change of variables and parameters:

s = r 1 t , u ¯ = u a , v ¯ = v a e , K 1 ¯ = K 1 a , α ¯ = α e r 1 , d ¯ = d a e r 1 , d i ¯ = d i r 1 , m i ¯ = m i r 1 , i = 1 , 2 , (2)

and still denote s , u ¯ , v ¯ , K 1 ¯ , α ¯ , d ¯ , d i ¯ and m i ¯ as t , u , v , K 1 , α , d , d i and m i , i = 1 , 2 . Thus we have following simplified dimensionless system in the form of PDEs:

u t d 1 Δ u = f ( u , v ) : = u ( 1 η u ) α u v 1 + u m 1 u , t > 0 , x Ω , (3a)

v t d 2 Δ v = g ( u , v ) : = α u v 1 + u m 2 v d v 2 , t > 0 , x Ω , (3b)

ν u = ν v = 0 , t 0 , x ( Ω ) , (3c)

u ( 0 , x ) = u 0 ( x ) 0 , v ( 0 , x ) = v 0 ( x ) 0 , x Ω . (3d)

where η : = 1 K 1 .

Our main contribution in this paper is detailed global asymptotic stability proof and Hopf bifurcation analysis of the system (3). The rest of this paper is organized as follows. In Section 2, we will analyze global stability of trivial equilibria E 0 , E 2 and interior equilibrium E * by using the comparison principle. In Section 3, we firstly give standard stability analysis to show the nonexistence of Turing patterns of this system, then we conduct the Hopf bifurcation analysis to show the existence of oscillatory patterns. The directions of Hopf bifurcation are also performed analytically. Finally, a short summary and some remarks are in Section 4.

2. Global Asymptotic Stability

In this section, we devote to give priori foundations for our system. Firstly, we discuss non-negative equilibria of the system (3) with their sufficient existence conditions. It is obvious to see that this system has following trivial equilibria:

E 0 = ( 0,0 ) and E 2 = ( u 2 ,0 ) , where u 2 is defined as u 2 : = 1 m 1 η . For practical considerations, we omit a singular point E 1 = ( 0, m 2 d ) . The point E 2 is

a desired equilibrium only if 1 > m 1 . Then we make a special effort to derive the existence conditions of an interior equilibrium E * which is denoted by ( u * , v * ) or ( s 1 , s 2 ) . If 0 < u 0 < u 2 , an interior equilibrium E * exists. Meanwhile, we

have u 0 < u * < u 2 and 0 < v * < η 4 α ( 1 + u 2 ) 2 , where u 0 : = m 2 α m 2 [21].

Then we will give analysis of global asymptotic stability at the equilibria E 0 , E 2 and E * . These conclusions can also be extended to a generalized bounded domain Ω 1 with smooth boundary. Here u 0 ( x ) , v 0 ( x ) C 1 ( Ω ¯ ) and positive solutions u , v C 2,1 ( ( T , ) × Ω ) C 1,0 ( [ T , ) × Ω ¯ ) , where T 0 [22].

2.1. Equilibria E0 and E2

Firstly, we consider global asymptotic stability of the trivial equilibrium E 0 by using the comparison principle [22] or [23].

Theorem 1 (Global asymptotic stability at E0) If 1 m 1 , then the equilibrium E 0 is globally asymptotically stable.

Proof. With respect to the Equation (1a), it is obvious to get an inequality

u t d 1 Δ u u η ( u 2 u ) . (4)

By using lemmas in [22] or [23], we have lim sup t max Ω ¯ u = 0 . Thus for any sufficiently small ε > 0 , there must exist T 1 , such that x Ω ¯ , t T 1 , we have

v t d 2 Δ v v [ ( α m 2 ) ε m 2 1 + ε d v ] . (5)

This implies lim sup t max Ω ¯ v = 0 . Thus we complete the proof.

With the same technique at hand, we have following theorem about the axis equilibrium E 2 .

Theorem 2 (Global asymptotic stability at E2) If 1 > m 1 and

α < m 2 ( 1 + 1 u 2 ) , then the equilibrium E 2 is globally asymptotically stable.

Proof. For the Equation (3a), from the inequality (4) we have

lim sup t max Ω ¯ u u 2 . That is to say, for any sufficiently small ε > 0 , we have

u u 2 + ε . Thus from the Equation (3b), a similar inequality

v t d 2 Δ v v [ α ( u 2 + ε ) 1 + u 2 + ε m 2 d v ] . (6)

is derived. By using the lemmas illustrated above, we have lim sup t max Ω ¯ v = 0 , i.e.

an inequality v ε holds. Substitute it into the Equation (3a), we have a new inequality

u t d 1 Δ u u η ( u 2 α ε η u ) . (7)

This implies lim inf t min Ω ¯ u u 2 and positive solutions u converge uniformly

to u 2 in Ω ¯ . Thus we complete the proof.

2.2. Equilibrium E*

Here we consider the global stability of the equilibrium E * . Firstly, we define two functions

φ ( s ) : = ( α m 2 ) s m 2 d ( 1 + s ) , ψ ( s ) : = 1 2 η ( 1 m 1 η + Δ s ) , s 0 (8)

and a discriminant

Δ s = ( 1 m 1 η ) 2 + 4 η ( 1 m 1 α s ) 0.

Notice that φ ( s ) is a monotonic increasing function but ψ ( s ) is a monotonic decreasing function.

From the existence conditions of point E * and inequality (10), for sufficiently small ε > 0 , we have u u 1 ¯ + ε , where u 1 ¯ : = u 2 . Substitute it into the Equation (3b), we obtain an inequality

v t d 2 Δ v v 1 + u 1 ¯ + ε { [ α ( u 1 ¯ + ε ) m 2 ( 1 + u 1 ¯ + ε ) ] d ( 1 + u 1 ¯ + ε ) v } . (9)

This implies

lim sup t max Ω ¯ v v 1 ¯ : = φ ( u 1 ¯ ) (10)

and

u t d 1 Δ u η u 2 + u ( 1 m 1 η ) + ( 1 m 1 ) α ( v 1 ¯ + ε ) 1 + u u . (11)

For any sufficiently small ε > 0 , suppose now that the numerator in the right hand side of the inequality (11) has two different real roots u v 1 ¯ , + ε > 0 and u v 1 ¯ , ε < 0 , i.e.

1 m 1 α v 1 ¯ > 0 , (12)

where u v 1 ¯ , + ε = ψ ( v 1 ¯ + ε ) . Hence we have a positive lower bound

lim inf t min Ω ¯ u u 1 _ : = ψ ( v 1 ¯ ) . (13)

Similarly, we have lim inf t min Ω ¯ v v 1 _ : = φ ( u 1 _ ) . This positive lower bound enquires that

Δ v 1 ¯ > max { 0 , 2 η u 0 1 + m 1 + η } . (14)

With the same procedure at hand, we have new bounds u 2 ¯ : = ψ ( v 1 _ ) , .

Without loss of generality, for these positive lower and upper bound sequences, we have following iteration relations

v i _ = φ ( u i _ ) , v i ¯ = φ ( u i ¯ ) , u i _ = ψ ( v i ¯ ) , u i + 1 ¯ = ψ ( v i _ ) (15)

and comparison relations

0 < u i _ u i + 1 _ u i + 1 ¯ u i ¯ , 0 < v i _ v i + 1 _ v i + 1 ¯ v i ¯ .

It is obvious to see that four limitations lim i v i _ , lim i v i ¯ , lim i u i _ and lim i u i ¯

are all exist with the help of mathematical analysis. But we denote them as v _ , v ¯ , u _ and u ¯ , respectively. Notice that they satisfy following equations:

v _ = φ ( u _ ) , v ¯ = φ ( u ¯ ) , u _ = ψ ( v ¯ ) , u ¯ = ψ ( v _ ) , (16)

or

f ( u _ , v ¯ ) = 0 , f ( u ¯ , v _ ) = 0 , g ( u _ , v _ ) = 0 , g ( u ¯ , v ¯ ) = 0 , (17)

thus u _ = u ¯ is equivalent to v _ = v ¯ from above last two equations. If u _ = u ¯ holds or v _ = v ¯ holds, we know that u _ = u ¯ = u * and v _ = v ¯ = v * due to the existence condition of point E * .

Case 1. If condition d ( 1 m 1 + η ) ( α m 2 ) 2 holds, substitute the equation v _ = φ ( u _ ) into the equation u _ = ψ ( v ¯ ) , then we have

( 1 m 1 ) ( α m 2 d v _ ) η ( m 2 + d v _ ) v ¯ ( α m 2 d v _ ) 2 = 0.

Similarly, we have

( 1 m 1 ) ( α m 2 d v ¯ ) η ( m 2 + d v ¯ ) v _ ( α m 2 d v ¯ ) 2 = 0.

Let the two equations to subtract each other, we derive a contradiction! Thus we have a theorem.

Theorem 3 (Global asymptotic stability at E*) Suppose 0 < u 0 < u 2 , if conditions (12), (14) and d ( 1 m 1 + η ) ( α m 2 ) 2 hold, then E * is globally asymptotically stable.

Case 2. Substitute the equation v _ = φ ( u _ ) into the equation u i _ = ψ ( v i ¯ ) , then we have

d ( 1 η u _ m 1 ) ( 1 + u _ ) = α ( α u ¯ 1 + u ¯ m 2 ) .

Similarly, we have

d ( 1 η u ¯ m 1 ) ( 1 + u ¯ ) = α ( α u _ 1 + u _ m 2 ) .

Let the above two equations to divide each other, we derive p ( u _ ) = p ( u ¯ ) , where p ( s ) = η ( α m 2 ) ( u 2 s ) ( s u 0 ) is a quadratic function. If there exist an index i 0 , such that

[ u i 0 _ , u i 0 ¯ ] [ u 0 , u 0 + u 2 2 ] or [ u 0 + u 2 2 , u 2 ] ,

for instance, i 0 = 1 , i.e. u 1 _ u 0 + u 2 2 or

Δ v 1 ¯ α η α m 2 , (18)

then the equilibrium E * is also globally asymptotically stable.

Theorem 4 (Global asymptotic stability at E*) Suppose 0 < u 0 < u 2 , if conditions (12) and (18) hold, then E * is globally asymptotically stable.

3. Hopf Bifurcation

In this section, we concentrate on the Hopf bifurcation analysis. Firstly, we define a real-valued Sobolev space

X : = { ( u , v ) | u , v H 2 ( Ω ) ; ν u = ν v = 0 , x ( Ω ) }

and the complexification of X [18] as

X C : = X i X = { x 1 + i x 2 | x 1 , x 2 X } .

Denote the complex-valued L 2 ( Ω ) inner product , on space X C as U , V : = Ω U V d x , where column vectors U , V X C , and U : = U T ¯ is the Hermitian conjugate (or adjoint) vector of U, thus we notice that the space X C equipped with inner product , is a Hilbert space. It is easy to verify the “linear” relation λ U , V = λ ¯ U , V .

3.1. Nonexistence of Turing Instability

In this section, we consider the nonexistence of Turing instability of above positive constant steady state E * . Firstly, we recall the corresponding ODEs system of (3) again:

u t = f ( u , v ) , v t = g ( u , v ) , (19)

and the Jacobian matrix of the system (19) at E * reads

J ( E * ) = [ α u * v * ( 1 + u * ) 2 u * K 1 α u * 1 + u * α v * ( 1 + u * ) 2 d v * ] : = [ A B C D ] . (20)

Then we denote some notations δ : = u * , A 1 : = t r [ J ( E * ) ] = A + D and A 2 : = det [ J ( E * ) ] = A D B C , where A 1 and A 2 are the trace and determinant of the matrix J ( E * ) , respectively. If the parameter δ satisfies

( δ + 1 ) 2 α v * η , (21)

conditions A 1 < 0 and A 2 > 0 hold, we know that E * is asymptotically stable in the ODEs system (19).

The linearized operator of the system (3) at steady state E * is

L ^ ( E * ) = [ A + d 1 Δ B C D + d 2 Δ ] . (22)

Suppose now that Φ = ( ϕ ψ ) is an eigenvector of operator L ^ ( E * ) corresponding to an eigenvalue μ , i.e. L ^ ( E * ) Φ = μ Φ or equations

( A μ ) ϕ + d 1 Δ ϕ + B ψ = 0, C ϕ + ( D μ ) ψ + d 2 Δ ψ = 0, (23)

{ ε i j } 0 i < , 1 j m i is a basis and ϕ , ψ have following expressions [19]:

ϕ = i = 0 j = 1 m i a i j ε i j , ψ = i = 0 j = 1 m i b i j ε i j . (24)

Substitute them into above Equation (23), we have following algebraic linear equations

[ A μ d 1 λ i B C D μ d 2 λ i ] ( a i j b i j ) = 0. (25)

Set ( a i 0 , j 0 b i 0 , j 0 ) 0 or above linear equations have nonzero solutions for index

i 0 , j 0 , then the determinant of the equation must be zero, i.e. we have a characteristic equation

μ 2 P i 0 μ + Q i 0 = 0 , (26)

where

P i 0 = A 1 ( d 1 + d 2 ) λ i 0 < 0 , Q i 0 = ( A d 1 λ i 0 ) ( D d 2 λ i 0 ) B C .

With the condition (21) at hand, we know Q i 0 > 0 , i.e. R e ( μ ) < 0 , then the steady state E * is asymptotically stable in PDEs system (3). Turing instability will not occur. See following theorem to summarize above analysis.

Theorem 5 (Turing instability) Suppose that conditions 0 < u 0 < u 2 and (21) hold, then E * is asymptotically stable in the PDEs system (3) and Turing instability phenomenon will not occur.

3.2. Existence of Hopf Bifurcation and Spatial Periodic Patterns

In this section, with the help of standard Hopf bifurcation theory, we will prove the existence of spatially homogeneous and non-homogeneous periodic patterns of the system (3). Here we choose δ = u * as a bifurcation parameter (or equivalently d as a bifurcation parameter). Firstly, take the linear transformation for later use: u ^ = u δ , v ^ = v v * , and still denote u ^ , v ^ as u, v, we have a new system

u t d 1 Δ u = F ( u , v ) : = f ( u + δ , v + v * ) , v t d 2 Δ v = G ( u , v ) : = g ( u + δ , v + v * ) , (27)

where smooth functions are

F ( u , v ) = ( u + δ ) [ 1 η ( u + δ ) ] α ( u + δ ) ( v + v * ) 1 + u + δ m 1 ( u + δ ) , G ( u , v ) = α ( u + δ ) ( v + v * ) 1 + u + δ m 2 ( v + v * ) d ( v + v * ) 2 .

Take the operator of (27) near ( 0,0 ) which determines the eigenvalues of linearized operator L ^ ( δ ) : = L ^ ( E * ) :

L ^ n ( δ ) = [ A ( δ ) d 1 n 2 l 2 B ( δ ) C ( δ ) D ( δ ) d 2 n 2 l 2 ] , (28)

where

A ( δ ) = δ ( 1 2 δ η m 1 η ) 1 + δ , B ( δ ) = α δ 1 + δ , C ( δ ) = 1 δ η m 1 1 + δ , D ( δ ) = m 2 α δ 1 + δ ,

then we derive its characteristic equation det [ L ^ n ( δ ) β I ] = 0 or an algebraic equation

β 2 β T n ( δ ) + D n ( δ ) = 0, (29)

where coefficients are:

T n ( δ ) = A 1 ( δ ) ( d 1 + d 2 ) n 2 l 2 , D n ( δ ) = A 2 ( δ ) + d 1 d 2 n 4 l 4 A ( δ ) d 2 n 2 l 2 D ( δ ) d 1 n 2 l 2 .

From the theorem 5, we know that any potential Hopf bifurcation occurs

when ( δ + 1 ) 2 < α v * η or

δ < 1 η m 1 2 η . (30)

We shall identify possible Hopf bifurcation values δ H > 0 under this condition accompanied by 1 η m 1 > 0 .

Suppose now that the bifurcation parameter δ satisfies the condition (30), and A ( δ ) ± i w ( δ ) are a pair of conjugate eigenvalues of L ^ n ( δ ) , i.e.

A ( δ ) = T n ( δ ) 2 or

A ( δ ) = δ η δ 2 m 1 δ α δ 2 ( 1 + δ ) δ η 2 + m 2 2 ( d 1 + d 2 ) n 2 2 l 2 (31)

with

A ( δ ) = 1 4 η δ 2 η δ 2 m 1 α η 2 ( 1 + δ ) 2 . (32)

Hence the transversality condition is satisfied as long as δ δ * (if δ * > 0 exist) or δ * > 0 can not holds, where

δ * = η + 1 m 1 α 2 η 1.

Let A ( δ ) = 0 , we obtain potential Hopf bifurcation points

δ j , ± H = 1 m 1 η α + m 2 ( d 1 + d 2 ) j 2 l 2 ± Δ j 4 η , j = 0,1,2, (33)

with a discriminant

Δ j = [ 1 m 1 η α + m 2 ( d 1 + d 2 ) j 2 l 2 ] 2 8 η [ ( d 1 + d 2 ) j 2 l 2 m 2 ] .

Note that the Hopf bifurcation at δ 0, ± H occurs without any restriction on l and δ 0, H is always non-positive. Now we only need to verify that D n ( δ j , ± H ) 0 , for instance, D n ( δ j , ± H ) > 0 holds forever.

Recall the condition (30) again, we have A ( δ ) > 0 , A ( δ ) ( A _ , A ¯ ) ; B ( δ ) < 0 , B ( δ ) ( B _ , B ¯ ) ; C ( δ ) > 0 , C ( δ ) ( C _ , C ¯ ) and D ( δ ) < 0 , D ( δ ) ( D _ , D ¯ ) . If there is a positive lower bound (or a local positive minimum) L > 0 such that A 2 ( δ ) L > 0 , for instance, A 2 ( δ ) > A ¯ D _ B ¯ C _ > 0 , then

D n ( δ j , ± H ) L + d 1 d 2 t 1 2 A ¯ d 2 t 1 , t 1 : = n 2 l 2 . (34)

Suppose further that the discriminant of above quadratic function in righthand side satisfies following condition

A ¯ 2 4 L < d 1 d 2 , (35)

then D n ( δ j , ± H ) > 0 since the quadratic function is always positive. Summarizing discussions above, we obtain following theorem.

Theorem 6 (Hopf bifurcation) Suppose that 0 < u 0 < u 2 and points δ j , ± H > 0 exist with 1 > η + m 1 , the parameters satisfy L > 0 and (35), then the system (3) undergoes a Hopf bifurcation at δ = δ j , ± H ( δ * ) , and the bifurcating periodic solutions can be parameterized (see the Formula (2.32) in [18]). Furthermore, we have:

(1) The bifurcation periodic solutions from δ = δ 0, + H are spatially homogeneous, which coincides with the periodic solutions of the corresponding ODEs system;

(2) The bifurcation periodic solutions from δ = δ j , ± H ( j 0 ) are spatially non-homogeneous.

Example From the existence condition of point E * and the condition (30), we have

A ( δ ) < A ¯ : = m 2 α ( 1 m 1 η ) , B ( δ ) < m 2 , C ( δ ) > η , D ( δ ) > m 2 α , (36)

then the positive lower bound is

L = m 2 ( 1 m 1 η ) ( m 2 α ) α + m 2 η .

Here we let l = 50 , d 1 = 5 , d 2 = 1 , K 1 = 17 , m 1 = 15 / 17 , m 2 = 1 and α = 15 , from some complicated calculations, the Hopf bifurcation values are δ 0, + H 0.07168660 and δ 1, + H 0.07150236 .

3.3. Direction of Hopf Bifurcation

Under the given conditions in above subsection, by the center manifold theorem and the normal form theory [24], the system (3) has a series of periodic solutions. In this subsection, we consider the direction and stability of spatially non-homogeneous periodic solutions at δ = δ j , ± H corresponding to j = 0 and j 0 , respectively. Here we obey the framework in references [24] (Chapter 5), [18] and [19], and only need to calculate R e [ c 1 ( δ j , ± H ) ] . For convenience, we denote δ j = δ j , ± H , w j = w ( δ j ) , A = A ( δ j ) , B = B ( δ j ) , C = C ( δ j ) and D = D ( δ j ) , where

w j 2 = α δ j ( 1 δ j η m 1 ) ( 1 + δ j ) 2 ( m 2 α δ j 1 + δ j d 2 j 2 l 2 ) 2 . (37)

To summarize above discussions, we firstly give following Hopf bifurcation theorem for our diffusive predator-prey system.

Theorem 7 (Direction of Hopf bifurcation) For the diffusive system (3), suppose that the theorem 6 holds, then Hopf bifurcation at point δ = δ j is supercritical (subcritical) if following number

σ ( δ j ) : = R e [ c 1 ( δ j ) ] A ( δ j ) < 0 ( > 0 ) . (38)

Moreover, we have:

(1) The bifurcating (spatially homogeneous) periodic solutions are stable (unstable) at δ 0 if R e [ c 1 ( δ 0 ) ] < 0 ( > 0 ) ;

(2) The bifurcating periodic solutions are all unstable at δ j ( j 0 ) .

3.3.1. The General Case: j 0

For the operators L ^ j ( δ j ) , we take an eigenvector q = ( a j b j ) φ j ( x ) and a “conjugate” vector q * = ( a j * b j * ) φ j ( x ) , such that q * , q = 1 and q * , q ¯ = 0 , where φ j ( x ) : = cos j x l and

a j = 1 , b j = A d 1 j 2 l 2 i w j B , a j * = i B b j w j l π , b j * = i B w j l π .

For later use, from functions F ( u , v ) and G ( u , v ) , we obtain partial derivatives evaluated at ( 0,0 ; δ j ) as follows:

F u u = 2 η ( 1 + δ j ) 3 + 2 α s 2 ( 1 + δ j ) 3 , F v v = 0 , F u v = α ( 1 + δ j ) 2 ; F u u u = 6 α s 2 ( 1 + δ j ) 4 , F v v v = 0 , F u u v = 2 α ( 1 + δ j ) 3 , F u v v = 0 ; G u u = 2 η F u u , G v v = 2 d , G u v = F u v ; G u u u = F u u u , G v v v = 0 , G u u v = F u u v , G u v v = 0 ,

and vectors in the form of symmetric B, C functions (see [25]):

B ( q , q ) = ( c j d j ) φ j 2 , B ( q , q ¯ ) = ( e j f j ) φ j 2 and C ( q , q , q ¯ ) = ( g j h j ) φ j 3 , where coefficients are

c j = F u u + 2 F u v b j , d j = G u u + 2 G u v b j + G v v b j 2 , e j = F u u + F u v ( b j ¯ + b j ) , f j = G u u + G u v ( b j ¯ + b j ) + G v v | b j | 2 , g j = F u u u + F u u v ( 2 b j + b j ¯ ) , h j = g j .

Notice that the integrals Ω φ j 3 d x = 0 , it is straightforward to drive relations

q * , B ( q , q ) = q * , B ( q , q ¯ ) = q * ¯ , B ( q , q ) = q * ¯ , B ( q , q ¯ ) = 0. (39)

So far, we have H 20 = B ( q , q ) and H 11 = B ( q , q ¯ ) . From following inverse operators

[ 2 i w j I L ^ 2 j ( δ j ) ] 1 = 1 α 1 + i α 2 [ 2 i w j D + 4 d 2 j 2 l 2 B C 2 i w j A + 4 d 1 j 2 l 2 ] , (40a)

[ 2 i w j I L ^ 0 ( δ j ) ] 1 = 1 α 3 + i α 4 [ 2 i w j D B C 2 i w j A ] , (40b)

L ^ 2 j ( δ j ) 1 = 1 α 5 [ D 4 d 2 j 2 l 2 B C A 4 d 1 j 2 l 2 ] , (40c)

L ^ 0 ( δ j ) 1 = 1 α 6 [ D B C A ] , (40d)

where

α 1 = 3 ( A D B C ) + 12 d 1 d 2 j 4 l 4 , α 2 = 2 w j [ A D + 4 ( d 1 + d 2 ) j 2 l 2 ] , α 3 = 4 w j 2 + A D B C , α 4 = 2 w j ( A + D ) , α 5 = ( A 4 d 1 j 2 l 2 ) ( D 4 d 2 j 2 l 2 ) B C , α 6 = A D B C ,

then we have

W 20 = 1 2 { [ 2 i w j I L ^ 2 j ( δ j ) ] 1 φ 2 j + [ 2 i w j I L ^ 0 ( δ j ) ] 1 } ( c j d j ) : = ( λ μ ) φ 2 j + ( τ χ ) , W 11 = 1 2 [ L ^ 2 j ( δ j ) 1 φ 2 j + L ^ 0 ( δ j ) 1 ] ( e j f j ) : = ( λ ˜ μ ˜ ) φ 2 j + ( τ ˜ χ ˜ ) . (41)

These calculations of inverse operators in W 20 and W 11 are restricted to the subspaces spanned by the eigen-modes 1 and φ 2 j . Precisely, we have

λ = 1 2 ( α 1 + i α 2 ) [ ( 2 i w j D + 4 d 2 j 2 l 2 ) c j + B d j ] , μ = 1 2 ( α 1 + i α 2 ) [ C c j + ( 2 i w j A + 4 d 1 j 2 l 2 ) d j ] , τ = 1 2 ( α 3 + i α 4 ) [ ( 2 i w j D ) c j + B d j ] , χ = 1 2 ( α 3 + i α 4 ) [ C c j + ( 2 i w j A ) d j ] ,

λ ˜ = 1 2 α 5 [ B f j ( D 4 d 2 j 2 l 2 ) e j ] , μ ˜ = 1 2 α 5 [ C e j ( A 4 d 1 j 2 l 2 ) f j ] , τ ˜ = 1 2 α 6 ( B f j D e j ) , χ ˜ = 1 2 α 6 ( C e j A f j ) ,

and (41) yield

B ( W 20 , q ¯ ) = ( λ F u u + ( λ b j ¯ + μ ) F u v λ G u u + ( λ b j ¯ + μ ) G u v + μ b j ¯ G v v ) φ 2 j φ j + ( λ F u u + ( λ b j ¯ + μ ) F u v τ G u u + ( τ b j ¯ + χ ) G u v + χ b j ¯ G v v ) φ j , B ( W 11 , q ) = B ( W 20 , q ¯ ) { λ λ ˜ , μ μ ˜ , τ τ ˜ , χ χ ˜ , b j ¯ b j } . (42)

Since the integrals of φ j are

Ω φ j 2 d x = 1 2 l π , Ω φ 2 j φ j 2 d x = 1 4 l π , Ω φ j 4 d x = 3 8 l π , (43)

by some calculations, we obtain

q * , B ( W 20 , q ¯ ) = l π 4 { a j * ¯ [ λ F u u + ( λ b j ¯ + μ ) F u v ] + b j * ¯ [ λ G u u + ( λ b j ¯ + μ ) G u v + μ b j ¯ G v v ] } + l π 2 { a j * ¯ [ τ F u u + ( τ b j ¯ + χ ) F u v ] + b j * ¯ [ τ G u u + ( τ b j ¯ + χ ) G u v + χ b j ¯ G v v ] } (44)

and

q * , B ( W 11 , q ) = q * , B ( W 20 , q ¯ ) { λ λ ˜ , μ μ ˜ , τ τ ˜ , χ χ ˜ , b j ¯ b j } , q * , C ( q , q , q ¯ ) = 3 8 l π ( a j * ¯ b j * ¯ ) g j . (45)

Finally, from the Formula (2.31) in [18] or page 47 in [24], we have

R e [ c 1 ( δ j ) ] = R e [ q * , B ( W 11 , q ) ] + 1 2 R e [ q * , B ( W 20 , q ¯ ) ] + 1 2 R e [ q * , C ( q , q , q ¯ ) ] = p , q = u , v ( F p q f p q + G p q g p q ) + F u u u f u u u + F u u v f u u v , (46)

where coefficients f u u , , g u v , are

f u u = 1 8 l 2 w j { [ ( 2 λ ˜ + 2 τ R + 4 τ ˜ + λ R ) w j + 2 A ( τ I + ( 1 / 2 ) λ I ) ] l 2 2 j 2 ( τ I + ( 1 / 2 ) λ I ) d 1 } , f v v = 0 , f u v = 1 8 l 4 w j B { ( ( 2 τ I + λ I ) w j 2 + ( ( 2 μ ˜ + 4 χ ˜ + μ R + 2 χ R ) B 4 A ( τ R + ( 1 / 2 ) λ R ) ) w j 2 A ( ( ( 1 / 2 ) μ I χ I ) B + A ( τ I + ( 1 / 2 ) λ I ) ) ) l 4 + 4 j 2 d 1 ( ( τ R + ( 1 / 2 ) λ R ) w j + ( ( 1 / 4 ) μ I ( 1 / 2 ) χ I ) B + A ( τ I + ( 1 / 2 ) λ I ) ) l 2 2 j 4 ( τ I + ( 1 / 2 ) λ I ) d 1 2 } ,

g u u = B ( 2 τ I + λ I ) 8 w j , g v v = 1 8 w j l 2 { ( ( 4 χ ˜ μ R 2 χ R + 2 μ ˜ ) w j A ( μ I + 2 χ I ) ) l 2 + j 2 d 1 ( μ I + 2 χ I ) } , g u v = 1 8 w j l 2 { ( ( 4 τ ˜ λ R + 2 λ ˜ 2 τ R ) w j A λ I 2 A τ I + B ( μ I + 2 χ I ) ) l 2 + j 2 d 1 ( 2 τ I + λ I ) } , f u u u = 3 16 , f u u v = 1 16 l 2 B [ ( 6 A 3 B ) l 2 + 6 d 1 j 2 ] ;

and decompositions of real part and imaginary part are: Γ R = R e ( Γ ) , Γ I = I m ( Γ ) , Γ = λ , μ , τ and χ , where

λ R = ( ( B d j D c j ) l 2 + 4 j 2 c j d 2 ) α 1 + 2 l 2 c j w j α 2 2 l 2 ( α 1 2 + α 2 2 ) , λ I = ( ( B d j + D c j ) l 2 4 j 2 c j d 2 ) α 2 + 2 l 2 c j w j α 1 2 l 2 ( α 1 2 + α 2 2 ) , μ R = ( ( A d j + C c j ) l 2 + 4 j 2 d 1 d j ) α 1 + 2 l 2 d j w j α 2 2 l 2 ( α 1 2 + α 2 2 ) , μ I = ( ( A d j C c j ) l 2 4 j 2 d 1 d j ) α 2 + 2 l 2 d j w j α 1 2 l 2 ( α 1 2 + α 2 2 ) ,

τ R = B α 3 d j D α 3 c j + 2 α 4 c j w j 2 ( α 3 2 + α 4 2 ) , τ I = B α 4 d j + D α 4 c j + 2 α 3 c j w j 2 ( α 3 2 + α 4 2 ) , χ R = A α 3 d j + C α 3 c j + 2 α 4 d j w j 2 ( α 3 2 + α 4 2 ) , χ I = A α 4 d j C α 4 c j + 2 α 3 d j w j 2 ( α 3 2 + α 4 2 ) .

3.3.2. The Special Case: j = 0

In this subsection, we consider the special case: j = 0 . Similarly, we take two

vectors q = ( a 0 b 0 ) and q * = ( a 0 * b 0 * ) , where

a 0 = 1 , b 0 = i w 0 A B , a 0 * = w 0 + i A 2 w 0 l π , b 0 * = i B 2 w 0 l π .

Suppose that B ( q , q ) = ( c 0 d 0 ) , B ( q , q ¯ ) = ( e 0 f 0 ) and C ( q , q , q ¯ ) = ( g 0 h 0 ) , where

c 0 = F u u + 2 F u v b 0 , d 0 = G u u + 2 G u v b 0 + G v v b 0 2 , e 0 = F u u 2 A B F u v , f 0 = G u u + G u v ( b 0 ¯ + b 0 ) + G v v | b 0 | 2 , g 0 = F u u u + F u u v ( 2 b 0 + b 0 ¯ ) , h 0 = g 0 ,

it is straightforward to see H 20 = H 11 = 0 and W 20 = W 11 = 0 , thus we have

q * , B ( W 11 , q ) = q * , B ( W 20 , q ¯ ) = 0 (47)

and (see the Formula (2.31) in [18])

c 1 ( δ 0 ) = i 2 w 0 q * , B ( q , q ) q * , B ( q , q ¯ ) + 1 2 q * , C ( q , q , q ¯ ) . (48)

From following calculations of inner product:

q * , B ( q , q ) = l π ( a 0 * ¯ c 0 + b 0 * ¯ d 0 ) , q * , B ( q , q ¯ ) = l π ( a 0 * ¯ e 0 + b 0 * ¯ f 0 ) , q * , C ( q , q , q ¯ ) = l π ( a 0 * ¯ b 0 * ¯ ) g 0 , (49)

we have the real part

R e [ c 1 ( δ 0 ) ] = 1 4 w 0 2 B 2 1 i + j 3 f i j A i B j , (50)

where some coefficients f i j are

f 30 = 2 F u v 2 + F u v G v v G v v 2 , f 21 = 3 F u u F u v + 3 G u v G v v , f 12 = F u u 2 F u u G u v F u v G u u G u u G v v 2 G u v 2 , f 11 = 2 F u u v w 0 2 , f 10 = 2 F u v 2 w 0 2 + F u v G v v w 0 2 G v v 2 w 0 2 , f 03 = F u u G u u + G u u G u v , f 02 = F u u u w 0 2 F u u v w 0 2 , f 01 = F u u F u v w 0 2 + G u v G v v w 0 2 ,

and coefficients f i j unlisted here are zero.

4. Summary and Remarks

In summary, with the framework of homogeneous reaction-diffusion systems, we have considered global asymptotic stability and Hopf bifurcation in a homogeneous diffusive predator-prey system with Holling type II functional response subject to Neumann boundary conditions, which is also an extended version of the predator-prey system in [18]. Some sufficient results were obtained to ensure that the equilibria of this system were globally asymptotically stable and Hopf bifurcation could occur. In the Example R e [ c 1 ( δ j , + H ) ] is negative while σ ( δ j , + H ) is positive due to the non-existence of δ * , j = 0,1 . That is to say, the bifurcation directions are subcritical at δ j , + H ( j = 0,1 ) ; the bifurcating periodic solutions are stable at δ 0, + H .

In Section 2, we induced global asymptotic stability theorems but neglected critical cases due to the used lemmas, which need to be considered further. In Subsections 3.1 and 3.2, combing the phenomenon that Turing instability will not occur, more sufficient conditions could be used to ensure asymptotic stability and existence of Hopf bifurcation, such as the conditions of Theorem 7 in reference [21], but the condition (36) is well-done for the Hopf bifurcation analysis. In Subsection 3.3, similar to the references listed above, we derived a useful algorithm for determining direction of Hopf bifurcation and stability of bifurcating periodic solutions. Furthermore, in [18] and [19], interior equilibria E * are all analytically and easily solvable, but the interior equilibrium E * in our system can not be solved easily. The methods in this paper are forward guidance for other complicated reaction-diffusion consumer-resource (predator-prey) systems, even some general reaction-diffusion systems in other fields. Finally, in some extent, it is our expectancy that these conclusions can provide theoretical support for more complex problems in biomathematics.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No.31570364) and the National Key Research and Development Program of China (Grant No.2018YFE0103700). We thank the Editor and the referee for their works.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1] Holling, C.S. (1965) The Functional Response of Predator Density and Its Role in Mimicry and Population Regulations. Memoirs of the Entomological Society of Canada, 45, 3-60.
https://doi.org/10.4039/entm9745fv
[2] Pei, Y.Z., Chen, L.S., Zhang, Q.R. and Li, C.G. (2005) Extinction and Performance of One-Prey Multi-Predators of Holling Type II Function Response System with Impulsive Biological Control. Journal of Theoretical Biology, 235, 495-503.
https://doi.org/10.1016/j.jtbi.2005.02.003
[3] Misha, P. and Raw, S.N. (2019) Dynamical Complexities in a Predator-Prey System Involving Teams of Two Prey and One Predator. Journal of Applied Mathematics and Computing, 61, 1-24.
https://doi.org/10.1007/s12190-018-01236-9
[4] Huang, J.C., Ruan, S.G. and Song, J. (2014) Bifurcations in a Predator-Prey System of Leslie Type with Generalized Holling Type III Functional Response. Journal of Differential Equations, 257, 1721-1752.
https://doi.org/10.1016/j.jde.2014.04.024
[5] Zhang, C.M., Chen, W.C. and Yang, Y. (2006) Periodic Solutions and Global Asymptotic Stability of a Delayed Discrete Predator-Prey System with Holling II Type Functional Response. Journal of Systems Science and Complexity, 19, 449-460.
https://doi.org/10.1007/s11424-006-0449-x
[6] Rosenzweig, M.L. and MacArthur, R. (1963) Graphical Representation and Stability Conditions of Predator-Prey Interactions. The American Naturalist, 97, 209-223.
https://doi.org/10.1086/282272
[7] Hsu, S.B. and Huang, T.W. (1995) Global Stability for a Class of Predator-Prey Systems. SIAM Journal on Applied Mathematics, 55, 763-783.
https://doi.org/10.1137/S0036139993253201
[8] Freedman, H.I. (1979) Stability Analysis of a Predator Prey System with Mutual Interference and Density Dependent Death Rate. Bulletin of Mathematical Biology, 41, 67-78.
https://doi.org/10.1016/S0092-8240(79)80054-3
[9] Wang Y.Q. and Jing, Z.J. (1999) Multiple Limit Cycles and Global Stability in Predator-Prey Model. Acta Mathematicae Applicatae Sinica, 15, 206-219.
https://doi.org/10.1007/BF02720497
[10] May, R.M. (1973) Stability and Complexity in Model Ecosystems. Monographs in Population Biology, Princeton University Press, Princeton.
[11] Du, Y.H. and Shi, J.P. (2007) Allee Effect and Bistability in a Spatially Heterogeneous Predator-Prey Model. Transactions of the American Mathematical Society, 359, 4557-4593.
https://doi.org/10.1090/S0002-9947-07-04262-6
[12] Pascual, M. (1993) Diffusion-Induced Chaos in a Spatial Predator-Prey System. Proceedings of the Royal Society of London B, 251, 1-7.
https://doi.org/10.1098/rspb.1993.0001
[13] Pang, P.Y.H. and Wang, M.X. (2003) Qualitative Analysis of a Ratio-Dependent Predator-Prey System with Diffusion. Proceedings of the Royal Society of Edinburgh, 133, 919-942.
https://doi.org/10.1017/S0308210500002742
[14] Yang, W.S. (2014) Dynamics of a Diffusive Predator-Prey Model with General Nonlinear Functional Response. The Scientific World Journal, 2014, Article ID: 721403.
https://doi.org/10.1155/2014/721403
[15] Li, S.B., Wu, J.H. and Dong, Y.Y. (2018) Effects of a Degeneracy in a Diffusive Predator-Prey Model with Holling Functional Response. Nonlinear Analysis: Real World Applications, 43, 78-95.
https://doi.org/10.1016/j.nonrwa.2018.02.003
[16] Peng, R. and Shi, J. (2009) Non-Existence of Non-Constant Positive Steady States of Two Holling Type-II Predator-Prey System: Strong Interaction Case. Journal of Differential Equations, 247, 866-886.
https://doi.org/10.1016/j.jde.2009.03.008
[17] Ko, W. and Ryu, K. (2006) Qualitative Analysis of a Predator-Prey Model with Holling Type-II Functional Response Incorporating a Prey Refuge. Journal of Differential Equations, 231, 534-550.
https://doi.org/10.1016/j.jde.2006.08.001
[18] Yi, F.Q., Wei, J.J. and Shi, J.P. (2009) Bifurcation and Spatiotemporal Patterns in a Homogeneous Diffusive Predator-Prey System. Journal of Differential Equations, 246, 1944-1977.
https://doi.org/10.1016/j.jde.2008.10.024
[19] Wan, A.Y., Song, Z.Q. and Zheng, L.F. (2016) Patterned Solutions of a Homogenous Diffusive Predator-Prey System of Holling Type-III. Acta Mathematicae Applicatae Sinica (English Scries), 32, 1073-1086.
https://doi.org/10.1007/s10255-016-0628-z
[20] Bazykin, A.D. (1998) Nonlinear Dynamics of Interacting Populations. World Scientific, Singapore.
https://doi.org/10.1142/2284
[21] Wang, S.T., Yu, H.G., Dai, C.J. and Zhao, M. (2020) The Dynamical Behavior of a Certain Predator-Prey System with Holling Type II Functional Response. Journal of Applied Mathematics and Physics, 8, 527-547.
https://doi.org/10.4236/jamp.2020.83042
[22] Ye, Q.X., Li, Z.Y., Wang, M.X. and Wu, Y.P. (2011) Introduction to Reaction-Diffusion Equations. 2nd Edition, Science Press, Beijing.
[23] Wang, M.X. and Pang, P.Y.H. (2009) Qualitative Analysis of a Diffusive Variable-Territory Prey-Predator Model. Discrete Continuous Dynamical Systems, 23, 1061-1072.
https://doi.org/10.3934/dcds.2009.23.1061
[24] Hassard, B.D., Kazarinoff, N.D. and Wan, Y.H. (1981) Theory and Application of Hopf Bifurcation. Cambridge University Press, Cambridge.
[25] Kuznetsov, Y.A. (1998) Elements of Applied Bifurcation Theory. 2nd Edition, Springer-Verlag, New York.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.