Dynamical behaviour of an ecological system with Beddington-DeAngelis functional response

0.5 MB·English

by  Sahabuddin Sarwardi

Dynamical behaviour of an ecological system with
Beddington-DeAngelis functional response
5 Sahabuddin Sarwardi( ) , Md. Reduanur Mandal( ) and Nurul Huda Gazi( )
1 † ∗ † †
0 (†) Department of Mathematics, Aliah University, IIA/27, New Town
2 Kolkata - 700 156, West Bengal, India.

n email: [email protected]
a
J

0
2
Abstract

]
S The objective of this paper is to study the dynamical behaviour systematically of an
D ecologicalsystemwithBeddington-DeAngelisfunctionalresponsewhichavoidsthecriticism
occurred in the case of ratio-dependent functional response at the low population density
.
h
of both the species. The essential mathematical features of the present model have been
t
a analyzedthoroughlyintermsofthelocalandtheglobalstabilityandthebifurcationsarising
m
insomeselectedsituationsaswell. Thethresholdvaluesforsomeparametersindicatingthe
[ feasibility and the stability conditions of some equilibria are also determined. We show

that the dynamics outcome of the interaction among the species are much sensitive to the
1
system parameters and initial population volume. The ranges of the significant parameters
v
7 under which the system admits a Hopf bifurcation are investigated. The explicit formulae
4 fordeterminingthe stability,directionandotherpropertiesofbifurcatingperiodicsolutions
8 are also derived with the use of both the normal form and the central manifold theory (cf.
4
Carr[1]). Numericalillustrations areperformedfinally inorderto validate the applicability
0
of the model under consideration.
.
1
0
5 Mathematics Subject Classification: 92D25, 92D30, 92D40.
1
: Keywords: Ecological model; Stability; Hopf bifurcation; Limit Cycle; Center manifold; Nu-
v
merical Simulation.
i
X
r
a
1 Introduction
Mathematical model is an important tool in analyzing the ecological models. Ecological prob-
lems are challenging and important issues from both the ecological and the mathematical point
ofview (cf. AndersonandMay [2], Beretta andKuang[3], Freedman[4], Hadeler andFreedman
[5], Hethcote et al. [6], Ma and Takeuchi [7], Venturino [8], Xiao and Chen [9]). The dynamic
relationship between predator and its prey has long been and will continue to be one of the
dominant themes in both ecology and mathematical ecology due to its universal existence and
importance. The most common method of modelling that ecological interactions consists of
two differential equations with simple correspondence between the consumption of prey by the
admissible predator and their population growth. The traditional predator-prey models have
∗Authorto whom all correspondence should beaddressed
1
been studied extensively (cf. Cantrell and Cosner [10], Cosner et al. [11], Cui and Takeuchi
[12], Huo et al. [13] and Hwang [14]), but those are questioned by several biologists. The
most crucial element in these models is the “functional response”–the expression that describes
the rate at which the number of prey consumed by a predator. Modifications were limited to
replacing the Malthusian growth function, the predator per capita consumption of prey func-
tions such as Holling type I, II, III functional responses or density dependent mortality rates.
These functional responses depend only on the prey volume x, but soon it became clear that
the predator volume y can influence this function by direct interference while searching or by
pseudo interference (cf. Curds and Cockburn [15], Hassell and Varley [16] and Salt [17]). A
simple way of incorporating predator dependence in the functional response was proposed by
Arditi and Ginzburg [18], who considered this response function as a function of the ratio x/y.
The ratio-dependent response function produces richer dynamics than that of all the Holling
types responses, but it is often criticized that the paradox occurred at the low densities of both
populations size. Normally one would expect that the population growth rate decrease when
both the populations fall bellow some critical volume, because food-searching effort becomes
very high. For some ecological interaction ratio-dependent model give the negative feed back.
Thus, the Lotka-Volterra type predator-prey model with the Beddington-DeAngelis functional
response has been proposed and well studied. Keeping these in mind, the proposed model can
be expressed as follows:
dx′1 =rx 1 x′1 c1x′1x′2
 dddxtt′2 =−δ′11(cid:0)x′2−+ ka1c(cid:1)+1e−x1′1x+a′1bx11+′2xx′2′1+b1x′2 (1.1)

0 0
with the initial conditions x (0) = x > 0 and x (0) = x >0. The functions x (t), x (t) are
′1 ′1 ′2 ′2 ′1 ′2
the volumes of prey and predator at any time t. All the system parameters are assumed to be
positive and have their usualbiological meanings. Thefunctionalresponse c1x′1x′2 insystem
a1+x′1+b1x′2
(1.1) was introduced by Beddington [19] and DeAngelis et al. [20] as a solution of the observed
problem in the classic predator-prey theory. It is similar to the well-known Holling type-II
functional response but has an extra term b x in the denominator which models mutual inter-
1 2
ference between predators. It represents the most qualitative features of the ratio-dependent
models, but avoids the “low-densities problem”, which usually the source of controversy. It
can be derived mechanistically from considerations of time utilization (cf. Beddington [19]) or
spatial limits on predation.
The present study under consideration has been carried out sequentially in the latter sections
as follows: The basic assumptions and the model formation are proposed in Section 2. Section
3 deals with some preliminary results. The equilibria and their feasibility are rightly given in
Section 4. The local analyses of the system around the boundary as well as interior equilibria
are discussed in Section 5. The global analysis of the system around the interior equilibrium
is studied at length in Section 6. Simulation results are reported in Section 7 while a final
discussion and interpretation of the results of the present study in ecological terms are rightly
included in the concluding Section 8.
2
2 Model formulation
Firstly we replaced the logistics growth function rx (1 x1) of the prey species by the modified
1 − k
quasi-linear growth function rx (1 x1 ) = r( k )x = r x (r < r) in order to make the
1 − x1+k k+x1 1 ′ 1 ′
model free from any axial equilibrium. Which fits better for some special type of ecosystem,
whereof environmental carrying capacity varies w.r.t. its prey volume, i.e., carrying capacity
is always greater than its present prey volume. In the present model we introduce one more
predator species in the model (1.1) to make it one step closure to reality. Thus, our final model
is extended to the following form:
dx1 = rx (1 x1 ) c1x1x2 c2x1x3
dt 1 − x1+k − a1+x1+b1x2 − a2+x1+b2x3
 dx2 = δ x + c1e1x1x2 , (2.1)
 ddxt3 = −δ1x2+ a1c+2ex21x+1bx13x2
dt − 2 3 a2+x1+b2x3


where x is the population volume of the two prey species and x , x are the population
1 2 3
volumes of the predator species at any time t. It is assumed that all the system parameters
are positive constants. Here r and k are the growth rate and the half-saturation constant for
the prey species, δ , δ are the first and second predators death rate respectively. c , c are
1 2 1 2
the respective search rates of the first and second predator on the prey species, c1, c2 are the
a1 a2
maximum number of prey that can be eaten by the first and second predator per unit time
respectively; 1 , 1 being their respective half saturation rates while e , e are the conversion
a1 a2 1 2
factors, denoting the number of newly born first and second predator for each captured prey
species respectively (0 < e , e < 1). The parameters b and b measure the coefficients
1 2 1 2
of mutual interference among the first and second predator species respectively. The terms
c1x1x2 and c2x1x3 denote the respective predator responseson the firstand second prey
a1+x1+b1x2 a2+x1+b2x3
species. This type of predator response function is known as Beddington-DeAngelis response
function (cf. Beddington [19] and DeAngelis et al. [20]).
3 Some preliminary results
3.1 Existence and positive invariance
Letting, x = (x ,x ,x )t, f : R3 R3, F = (f ,f ,f )t, the system (2.1) can be rewritten as
1 2 3 1 2 3
→
x˙ = f(x). Here f C (R) for i= 1,2,3, where f = rx (1 x1 ) c1x1x2 c2x1x3 ,
i ∈ ∞ 1 1 −x1+k −a1+x1+b1x2 −a2+x1+b2x3
f = δ x + c1e1x1x2 and f = δ x + c2e2x1x3 . Since the vector function f is a smooth
2 − 1 2 a1+x1+b1x2 3 − 2 3 a2+x1+b2x3
function of the variables (x ,x ,x ) in the positive octant Ω0 = (x ,x ,x ) : x > 0,x >
1 2 3 1 2 3 1 2
{
0,x >0 , the local existence and uniqueness of the solution of the system (2.1) hold.
3
}
3.2 Persistence
If a compact set D Ω0 = (x ,x ,x ) : x > 0, i = 1,2,3 exists such that all solutions of
1 2 3 i
⊂ { }
(2.1) eventually enter and remain in D, the system is called persistent.
Proposition 3.1. The system (2.1) is persistent if the conditions: (i)r > δ +δ , (ii) x >
1 2 11
a2δ2 , (iii)x > a1δ1 are satisfied.
c2e2 δ2 12 c1e1 δ1
− −
3
Proof. Weusethemethodofaverage Lyapunovfunction(cf. GardandHalam [21]), considering
a function of the form
V(x ,x ,x )= xγ1xγ2xγ3,
1 2 3 1 2 3
where γ ,γ and γ are positive constants to be determined. We define
1 2 3
V˙
Π(x ,x ,x ) =
1 2 3
V
rx c x c x c e x
1 1 2 2 3 1 1 1
= γ r +γ δ +
1 2 1
− x +k − a +x +b x − a +x +b x − a +x +b x
(cid:16) 1 1 1 1 2 2 1 2 3(cid:17) (cid:16) 1 1 1 2(cid:17)
c e x
2 2 1
+γ δ + .
3 2
− a +x +b x
(cid:16) 2 1 2 3(cid:17)
We now prove that this function is positive at each boundary equilibrium. Let γ = γ, for
i
i = 1, 2, 3. In fact at E , we have Π(0,0,0) = r δ δ > 0 from the condition (i). Moreover,
0 1 2
− −
from condition (ii) and (iii), we find the values of Π at E and E respectively,
1 2
e c x
Π(x ,x ,0) =γ δ + 2 2 11 > 0,
11 21 − 2 a +x
(cid:16) 2 11(cid:17)
e c x
Π(x ,0,x ) = γ δ + 1 1 12 > 0.
12 32 − 1 a +x
(cid:16) 1 12(cid:17)
Hence, there always exists a positive number γ such that Π > 0 at the boundary equilibria.
Hence V is an average Lyapunov function and thus, the system (2.1) is persistent.
Since the system is uniformly persistent, there exists σ > 0 and τ > 0 such that x (t) > σ, for
i
all t > τ, i = 1, 2, 3.
3.3 Boundedness
Boundedness implies that the system is consistent with biological significance. The following
propositions ensure the boundedness of the system (2.1).
Proposition 3.2. The prey population is always bounded from above.
Proof. Before proving that the prey population is bounded above, we need to prove that the
predator populations x and x areboundedabove. To prove this result, consideringthesecond
2 3
sub equation of the system (2.1) and one can obtain the following differential inequality:
dx
2
(δ c e )x .
1 1 1 2
dt ≤ − −
Integrating the above differential inequality between the limits 0 and t, we have x (t)
2
≤
x (0)e (δ1 c1e1)t. Thus, if (δ c e ) > 0, then it is obviously found a positive number τ
2 − − 1 1 1 1
−
there exists a positive constant m such that x (t) m , for all t τ . By using the similar
1 2 1 1
≤ ≥
argument, one can obtain that, if (δ c e ) > 0, then corresponding to a positive number τ
2 2 2 2
−
there exists a positive constant m such that x (t) m , for all t τ . Both the results can
2 3 2 2
≤ ≥
be written unitedly as x > m = min(m ,m ) for all t > τ = max(τ ,τ ), i = 2,3, with the
i 1 2 3 1 2
additional condition min(δ c e , δ c e ) > 0.
1 1 1 2 2 2
− −
4
Now from the first sub-equation of (2.1), the following inequality is found
dx (e +e )σ rk x k rm (e +e )σ
1 1 2 1 1 2
− − x .
dt ≤ (cid:0) km (cid:1) (cid:0)(e +e )σ rk (cid:1) − 1
(cid:16) 1 2 (cid:17)
−
Hence, by using standard but simple argument, we have
krm (e +e )kσ rm rk
1 2
limsupx (t) − = w, where < σ < .
1
≤ (e +e )σ rk e +e e +e
t + 1 2 1 2 1 2
→ ∞ −
Proposition 3.3. The solutions of (2.1) starting in Ω0 are uniformly bounded with an ultimate
bound.
Proof. Considering the total environment population χ = x + x2 + x3. Using the theorem
1 e1 e2
on differential inequality (cf. Birkhoff and Rota [22]) and following the steps of Haque and
Venturino [23], Sarwardi et al [24], boundedness of the solution trajectories of this model is
established. In particular,
x x (r+1)k+w
2 3
limsup x + + = M, where ρ = min(1,δ ,δ ), (3.1)
1 1 2
e e ≤ ρ
t + (cid:16) 1 2(cid:17)
→ ∞
with the last bound is independent of the initial condition.
Hence, all the solutions of (2.1) starting in R3 for any θ > 0 evolve with respect to time in the
+
compact region
x x
Ω¯ = (x ,x ,x ) R3 : x + 2 + 3 M +θ . (3.2)
(cid:26) 1 2 3 ∈ + 1 e e ≤ (cid:27)
1 2
4 Equilibria and their feasibility
The equilibria of the dynamical system (2.1) are given as follows:
1. The trivial equilibrium point E (0,0,0) is always feasible.
0
2. (a) The first boundary equilibrium point is E (x ,x ,0). The component x is a root
1 11 21 11
of the quadratic equation l x2 + (l + l k + rkb e )x + l k = 0, where l = (δ c e ),
1 11 2 1 1 1 11 2 1 1 − 1 1
l = a δ . If l < 0, then the quadratic equation in x possesses a unique positive root and
2 1 1 1 11
consequently x21 = (c1e1−δb11)δx111−a1δ1. The feasibility of the equilibrium E1 is maintained if the
condition x > b1δ1 is satisfied.
11 c1e1 δ1
−
2. (b) The second boundary equilibrium point is E (x ,0,x ). The component x is the root
2 12 32 12
of the quadratic equation m x2 +(m +m k+rkb e )x +m k = 0, where m = (δ c e ),
1 12 2 1 2 2 12 2 1 2− 2 2
l = a δ . If m < 0, then the quadratic equation in x possesses a unique positive root and
2 2 2 1 12
consequently x32 = (c2e2−δb22)δx212−a2δ2. The feasibility of the equilibrium E2 is maintained if the
condition x > b2δ2 holds.
12 c2e2 δ2
−
3. The interior equilibrium point is E (x ,x ,x ), where the first component x is the root
1 2 3 1
∗ ∗ ∗ ∗ ∗
5
of the following quadratic equation:
n x2 +(n +n k+rkb b e e )x +n k = 0, (4.1)
1 1 2 1 1 2 1 2 1 2
∗ ∗
where n = b e (δ c e )+b e (δ c e ) and n = b e δ a +b e δ a .
1 1 1 2 2 2 2 2 1 1 1 2 2 2 1 1 1 1 2 2
− −
Case I: Let n < 0. In this case there exists exactly one positive root of the quadratic equation
1
(4.1) irrespective of the sign of (n +n k+rkb b e e ).
2 1 1 2 1 2
Case II: Let n > 0. In this case there are two possibilities: (i) if n +n k+rkb b e e > 0,
1 2 1 1 2 1 2
then there is no positive solution and (ii) if n +n k +rkb b e e < 0, then there exists two
2 1 1 2 1 2
positive roots or no positive root.
In this presentanalysis we consider theCase I.Under this assumption thenext two components
of the interior equilibrium can be obtained as x2∗ = (c1e1−δb11)δx11∗−a1δ1, x3∗ = (c2e2−δb22)δx21∗−a2δ2.
The feasibility of this important equilibrium point E is confirmed under the condition x >
1
∗ ∗
max a1δ1 , a2δ2 . Moreover, the positivity condition of second and third components of
c1e1 δ1 c2e2 δ2
n − − o
the interior equilibrium ensures the impossibility of the Case II.
Remark: The feasibility and existences conditions of both the planer equilibria E and E
1 2
immediately implies the existence of the unique feasible interior equilibrium point E . But the
∗
existence of the unique feasible interior equilibrium point E implies three possibilities: (i) E
1
∗
exists and E does not exist, (ii) E exists and E does not exist, (iii) existence of both.
2 2 1
5 Local stability and bifurcation
The Jacobian matrix J(x) of the system (2.1) at any point x = (x ,x ,x ) is given by
1 2 3
rk2 c1x2(a1+b1x2) c2x3(a2+b2x3) c1x1(a1+x1) c2x1(a2+x1)
(x1+k)2 − (a1+x1+b1x2)2 − (a2+x1+b2x3)2 −(a1+x1+b1x2)2 −(a2+x1+b2x3)2
J(x) =  c1e1x2(a1+b1x2) δ + c1e1x1(a1+x1) 0 .
3×3 (a1+x1+b1x2)2 − 1 (a1+x1+b1x2)2
 c(2ae22+x3x(1a+2b+2bx23x)32) 0 −δ2+ (ca22e+2xx11(+a2b2+xx31))2
(5.1)
Its characteristic equation is ∆(λ)= λ3+k λ2+k λ+k = 0, where k = tr(J), k = M and
1 2 3 1 2
−
k = det(J); M being the sum of the principal minors of order two of J.
3
−
Note that the conditions for occurrence of Hopf bifurcation are that there exists a certain
bifurcation parameter r = r such that C (r ) = k (r )k (r ) k (r ) = 0 with k (r ) > 0 and
c 2 c 1 c 2 c 3 c 2 c
−
d (Re(λ(r))) = 0, where λ is root of the characteristic equation ∆(λ) = 0.
dr |r=rc 6
5.1 Local analysis of the system around E , E , E
0 1 2
Stability: The eigenvalues of the Jacobian matrix J(E ) are r, δ and δ . Hence E is
0 1 2 0
− −
unstable in nature (saddle point). Let J(E ) = (ξ ) and J(E ) = (η ) . Using the
1 ij 3 3 2 ij 3 3
× ×
RouthHurwitz criterion, it can be easily shown that the eigenvalues of the matrices J(E ) and
1
J(E2)willhavenegative realpartsifftheconditionse1x11+x21 > k(1−b1b1e1)−a1 ande2x12+x32 >
k(1−b2b2e2)−a2 respectively. HencetheequilibriaE1 andE2 arelocally asymptotically stableunder
6
theconditions e1x11+x21 > k(1−b1b1e1)−a1 ande2x12+x32 > k(1−b2b2e2)−a2 respectively (cf. Section
4 of of Sarwardi et al. [25]).
Bifurcation: Since the equilibrium point E is a saddle in nature, hence, there is no question
0
of Hopf bifurcation around this equilibrium. In order to have Hopf bifurcation around the
equilibria E , E , it is sufficient to show that the coefficient of λ in the quadratic factor of the
1 2
characteristic polynomial of J(E ) (k = 1,2) is zero and the constant term is positive. The
k
conditionsforwhichannihilationofthelineartermsinthequadraticfactorsofthecharacteristic
polynomials of J(E ) and J(E ) can be made possible are ξ +ξ = 0 and η +η = 0.
1 2 11 22 11 33
For detailed analysis, interested readers are referred to Appendix A of Haque and Venturion
[26]. The parametric regions where Hopf bifurcations occur around E and E are respectively
1 2
established bytheequality constraints e1x11+x21 = k(1−b1b1e1)−a1 ande2x12+x32 = k(1−b2b2e2)−a2.
5.2 Local analysis of the system around the interior equilibrium
Proposition 5.1. The system (2.1) around E is locally asymptotically stable if the condition
∗
(i) k < min a +b x , a +b x is satisfied.
1 1 2 2 2 3
{ ∗ ∗}
Proof. LetJ(x ) =(J ) is theJacobian matrix atthe interior equilibriumpointE = x of
ij 3 3
∗ × ∗ ∗
c1x1∗x2∗ k (a1+b1x2) c2x1∗x3∗ k (a2+b2x3)
thesystem(2.1). ThecomponentsofJ(x )areJ = − + − ,
∗ 11 (x1∗+k)((cid:0)a1+x1∗+b1x2∗(cid:1)) (x1∗+k)((cid:0)a2+x1∗+b2x3∗(cid:1))
J = c1x1∗(a1+x1∗) < 0, J = c2x1∗(a2+x1∗) < 0, J = c1e1x2∗(a1+b1x2∗) > 0, J =
12 −(a1+x1∗+b1x2∗)2 13 −(a2+x1∗+b2x3∗)2 21 (a1+x1∗+b1x2∗)2 22
b1c1e1x1∗x2∗ < 0, J = 0, J = c2e2x3∗(a2+b2x3∗) >0, J = 0, J = b2c2e2x1∗x3∗ < 0.
−(a1+x1∗+b1x2∗)2 23 31 (a2+x1∗+b2x3∗)2 32 33 −(a2+x1∗+b2x3∗)2
Then the characteristic equation of the Jacobian matrix J(x ) can be written as
∗
λ3+k λ2+k λ+k = 0, (5.1)
1 2 3
where k = tr(J) = (J +J +J ), k = M +M +M = (J J J J )+J J +
1 11 22 33 2 11 22 33 11 22 21 12 22 33
− − −
(J J J J ), k = det(J) = J J J J J J J J J , and C = k k k =
11 33 31 13 3 11 22 33 12 21 33 31 13 22 2 1 2 3
− − − − − −
(J +J ) J (J +J +J )+((cid:0)J J J J ) +J J (J +J(cid:1) ).
11 22 33 11 22 33 11 22 21 12 13 31 11 33
− −
(cid:0) (cid:1)
Itisclearthatk > 0ifJ < 0,i.e., k < min a +b x , a +b x andconsequentlyC > 0.
1 11 1 1 2 2 2 3 2
{ ∗ ∗}
HencetheRouth-Hurwitzconditionissatisfied forthematrix J ,i.e., allthecharacteristic roots
∗
of J are with negative real parts. So the system is locally asymptotically stable around E .
∗ ∗
Theorem 5.2. The dynamical system (2.1) undergoes Hopf bifurcation around the interior
equilibrium point E whenever the critical parameter value r = r contained in the domain
c
∗
dC
D = r R+ :C (r ) = k (r )k (r ) k (r )= 0 with k (r )> 0 and 2 = 0 .
HB c ∈ 2 c 1 c 2 c − 3 c 2 c dr |r=rc 6
n o
Proof. The equation (5.1) will have a pair of purely imaginary roots if k k k = 0 for
1 2 3
−
some set of values of the system parameters. Let us now suppose that r = r be the value of
c
r satisfying the condition k k k = 0. Here only J contains r explicitly. So, we write the
1 2 3 11
−
equation k k k = 0 as an equation in J to find r as follows:
1 2 3 11 c
−
h J2 +h J +h = 0, (5.2)
1 11 2 11 3
where h = J +J , h = J2 +J2 J J J J , h = (J +J )J J J J J
1 22 33 2 − 22 33− 13 31− 12 21 3 22 33 22 33− 13 31 33−
7
J J J .
12 21 22
Thus, J = 1 ( h h2 4h h ) = J .
11 2h1 − 2± 2− 1 3 1∗1
p
Or,
(x +k)2 c x (a +b x ) c x (a +b x )
1 1 2 1 1 2 2 3 2 2 3
r = ∗k2 (cid:20)J1∗1+ (a +∗x +b x ∗)2 − (a +∗x +b x ∗)2(cid:21) = rc. (5.3)
1 1 1 2 2 1 2 3
∗ ∗ ∗ ∗
Using the condition k k k = 0, from equation (5.1) one can obtain
1 2 3
−
(λ+k )(λ2 +k ) = 0, (5.4)
1 2
which has three roots λ = +i√k , λ = i√k , λ = k , so there is a pair of purely
1 2 2 2 3 1
− −
imaginary eigenvalues i√k . For all values of λ, the roots are, in general, of the form λ (r) =
2 1
±
ξ (r)+iξ (r), λ (r)= ξ (r) iξ (r), λ (r) = k (r).
1 2 2 1 2 3 1
− −
Differentiating the characteristic equation (5.1) w.r.t. r, we have
dλ λ2k˙ +λk˙ +k˙
1 2 3
=
dr −3λ2+2k λ+k |λ=i√k2
1 2
k˙ k k˙ +ik˙ √k
3 2 1 2 2
= −
2(k ik √k )
2 1 2
−
k˙ (k˙ k +k k˙ ) √k (k k˙ +k k˙ k k˙ k )
3 1 2 1 2 2 1 3 2 2 1 1 2
= − +i −
2(k2 +k ) 2k (k2+k )
1 2 2 1 2
dC2 √k k˙ k √k dC2
= dr +i 2 2 1 2 dr . (5.5)
−2(k2+k ) 2k − 2k (k2+k )
1 2 h 2 2 1 2 i
Hence,
d dC2
Re(λ(r)) = dr = 0. (5.6)
dr |r=rc −2(k2 +k )|r=rc 6
(cid:0) (cid:1) 1 2
d(Re(λ(r)))
Using the monotonicity condition of the real part of the complex root = 0 (cf.
dr |r=rc 6
Wiggins [27], pp. 380), one can easily establish the transversality condition dC2 = 0, to
dr |r=rc 6
ensure the existence of Hopf bifurcation around E .
∗
6 Global analysis of the system around the interior equilibrium
6.1 Direction of Hopf bifucation of the system (2.1) around E
∗
In this Section we study on the direction of Hopf bifucation around the interior equilibrium.
From the model equations (2.1), we have
x˙ = f(x), (6.1)
8
rx (1 x1 ) c1x1x2 c2x1x3
1 − x1+k − a1+x1+b1x2 − a2+x1+b2x3
where x = (x1,x2,x3)t, f = (f1,f2,f3)t =  −δ1x2+ a1c+1ex11x+1bx12x2 .
δ x + c2e2x1x3
 − 2 3 a2+x1+b2x3 
Here, at x = x , f = 0. Let y = (y ,y ,y ) = (x x , x x , x x ). Putting in
1 2 3 1 1 2 2 3 3
∗ − ∗ − ∗ − ∗
equation (6.1), we have
y˙ = J(x ,x ,x )y+φ, (6.2)
1 2 3
∗ ∗ ∗
where the components of nonlinear vector function φ= (φ ,φ ,φ )t are given by
1 2 3
φ = fi y2+fi y2+fi y2+2fi y y +2fi y y +2fi y y +h.o.t., i = 1,2,3. (6.3)
i x1x1 1 x2x2 2 x3x3 3 x2x3 2 3 x3x1 3 1 x1x2 1 2
The coefficients of nonlinear terms in y , i = 1,2,3 are given by
i
f1 = 2rk2 + 2c1x2(a1+b1x2) + 2c2x3(a2+b2x3), f1 = 2b1c1x1(a1+x1), f1 = 2b2c2x1(a2+x1),
x1x1 −(x1+k)3 (a1+x1+b1x2)3 (a2+x1+b2x3)3 x2x2 (a1+x1+b1x2)3 x3x3 (a2+x1+b2x3)3
f1 = 0, f1 = c1a1(a1+x1+b1x2)+2b1c1x1x2, f1 = c2a2(a2+x1+b2x3)+2b2c2x3x1; f2 =
x2x3 x1x2 − (a1+x1+b1x2)3 x3x1 − (a2+x1+b2x3)3 x1x1
2c1e1x2(a1+b1x2), f2 = 2b1c1e1x1(a1+x1), f2 = 0, f2 = a1c1e1(a1+x1+b1x2)+2b1c1e1x1x2,
− (a1+x1+b1x2)3 x2x2 − (a1+x1+b1x2)3 x3x3 x1x2 (a1+x1+b1x2)3
f2 = 0, f2 = 0; f3 = 2c2e2x3(a2+b2x3), f3 = 0, f3 = 2b2c2e2x1(a2+x1),
x2x3 x3x1 x1x1 − (a2+x1+b2x3)3 x2x2 x3x3 − (a2+x1+b2x3)3
f3 = 0, f3 = a2c2e2(a2+x1+b2x3)+2b2c2e2x1x3, f3 = 0.
x2x3 x3x1 (a2+x1+b2x3)3 x1x2
Let P be the matrix formed by the column vectors (u2,u1,u3), which are the eigenvec-
tors corresponding to the eigenvalues λ = i√k and λ = k of J(x ,x ,x ), then
1,2 2 3 1 1 2 3
± − ∗ ∗ ∗
J(x1 ,x2 ,x3 )u2 = i√k2u2, J(x1 ,x2 ,x3 )u1 = i√k2u1, and J(x1 ,x2 ,x3 )u3 = k1u3.
∗ ∗ ∗ ∗ ∗ ∗ − ∗ ∗ ∗ −
Thus,
0 1 1
 J21√k2 J21J22 J21 
P = −J222+k2 −J222+k2 J−22+k1 = (pij)3 3.
 J31√k2 J31J33 J31  ×
−J323+k2 −J323+k2 J−33+k1
Letus makeuseof thetransformation y = Pz, soas thesystem (6.2)is reducedto thefollowing
one
0 i√k 0
2
−
z˙ = P−1J(x1 ,x2 ,x3 )Pz+P−1φ = i√k2 0 0 z+P−1φ. (6.4)
∗ ∗ ∗
0 0 k
1
 − 
Here P 1 = AdjP = (q ) , where
− detP ij 3 3
×
q = 1 ( J21J22J31 J21J33J31 ), q = 1 ( J31 J31J33 ),
11 detP (J222+k2)(J33+k1) − (J222+k2)(J22+k1) 12 detP J323+k1 − J323+k2
q = 1 ( J21 + J21J22 ), q = 1 ( J21√k2J31 J21√k2J31 ),
13 detP J−21+k1 J222+k2 21 detP (J222+k2)(J22+k1) − (J323+k2)(J22+k1)
q = 1 √k2J31 , q = 1 √k2J31, q = 1 ( J21√k2J31J33 J21√k2J31J22 ),
22 detP (J323+k2) 23 detP −(J323+k2) 31 detP (J222+k2)(J323+k2) − (J323+k2)(J222+k2)
q = 1 √k2J21, q = 1 √k2J21 .
32 detP −(J222+k2) 33 detP (J222+k2)
The system (6.4) can be written as
d z 0 √k z
dt (cid:18) z1 (cid:19) = (cid:18)√k −0 2(cid:19)(cid:18) z1 (cid:19)+F(z1,z2,z3), (6.5)
2 2 2
dz
3
= k z +G(z ,z ,z ). (6.6)
1 3 1 2 3
dt −
9
On the center-manifold (cf. Carr [1], Kar [28])
1
z = b z2+2b z z +b z2 (6.7)
3 2 11 1 12 1 2 22 2
(cid:0) (cid:1)
Therefore,
0 √k z
z˙3 = (cid:0) b11z1+b12z2 b12z1+b22z2 (cid:1)(cid:18)√k2 −0 2(cid:19)(cid:18)z12(cid:19) = pk2b12z12+pk2(b22−b11)z1z2−pk2b12z22
(6.8)
Using (6.4) and (6.6), we have
z˙ = k z +q φ +q φ +q φ (6.9)
3 1 3 31 1 32 2 33 3
−
From the equations (6.8) and (6.9), we have
k b z2+ k (b b )z z k b z2
2 12 1 2 22 − 11 1 2− 2 12 2
p 1 p p 1
= k (b z2+2b z z +b z2)+q f1 p z +p z +p (b z2+2b z z +b z2) 2
−2 1 11 1 12 1 2 22 2 31 x1x1{ 11 1 12 2 132 11 1 12 1 2 22 2 }
(cid:2)
1 1
+f1 p z +p z +p (b z2+2b z z +b z2) 2 +f1 p z +p z +p (b z2
x2x2{ 21 1 22 2 232 11 1 12 1 2 22 2 } x3x3{ 31 1 32 2 332 11 1
1
+2b z z +b z2) 2 +2f1 p z +p z +p (b z2+2b z z +b z2) p z +p z
12 1 2 22 2 } x3x1{ 31 1 32 2 332 11 1 12 1 2 22 2 }{ 11 1 12 2
1 1
+p (b z2+2b z z +b z2) +2f1 p z +p z +p (b z2+2b z z +b z2)
132 11 1 12 1 2 22 2 } x1x2{ 11 1 12 2 132 11 1 12 1 2 22 2 }
1
p z +p z +p (b z2+2b z z +b z2)
×{ 21 1 22 2 232 11 1 12 1 2 22 2 }
(cid:3)
1
+q [f2 p z +p z +p (b z2+2b z z +b z2) 2+f2 p z +p z
32 x1x1{ 11 1 12 2 132 11 1 12 1 2 22 2 } x2x2{ 21 1 22 2
1 1
+p (b z2+2b z z +b z2) 2+f2 p z +p z +p (b z2+2b z z +b z2) 2
232 11 1 12 1 2 22 2 } x3x3{ 31 1 32 2 332 11 1 12 1 2 22 2 }
1 1
+2f2 p z +p z +p (b z2+2b z z +b z2) p z +p z +p (b z2
x2x3{ 21 1 22 2 232 11 1 12 1 2 22 2 }{ 31 1 32 2 332 11 1
1
+2b z z +b z2) +2f2 p z +p z +p (b z2+2b z z +b z2) p z +p z
12 1 2 22 2 } x3x1{ 31 1 32 2 332 11 1 12 1 2 22 2 }{ 11 1 12 2
1 1
+p (b z2+2b z z +b z2) +2f2 p z +p z +p (b z2+2b z z +b z2)
132 11 1 12 1 2 22 2 } x1x2{ 11 1 12 2 132 11 1 12 1 2 22 2 }
1
p z +p z +p (b z2+2b z z +b z2)
×{ 21 1 22 2 232 11 1 12 1 2 22 2 }
(cid:3)
1
+q f3 p z +p z +p (b z2+2b z z +b z2) 2 +f3 p z +p z +p
33 x1x1{ 11 1 12 2 132 11 1 12 1 2 22 2 } x2x2{ 21 1 22 2 23
(cid:2)
1 1
(b z2+2b z z +b z2) 2+f3 p z +p z +p (b z2+2b z z +b z2) 2
×2 11 1 12 1 2 22 2 } x3x3{ 31 1 32 2 332 11 1 12 1 2 22 2 }
1 1
+2f3 p z +p z +p (b z2+2b z z +b z2) p z +p z +p (b z2
x2x3{ 21 1 22 2 232 11 1 12 1 2 22 2 }{ 31 1 32 2 332 11 1
1
+2b z z +b z2) +2f3 p z +p z +p (b z2+2b z z +b z2) p z
12 1 2 22 2 } x3x1{ 31 1 32 2 332 11 1 12 1 2 22 2 }{ 11 1
1 1
+p z +p (b z2+2b z z +b z2) +2f3 p z +p z +p (b z2
12 2 132 11 1 12 1 2 22 2 } x1x2{ 11 1 12 2 132 11 1
1
+2b z z +b z2) p z +p z +p (b z2+2b z z +b z2) .
12 1 2 22 2 }{ 21 1 22 2 232 11 1 12 1 2 22 2 }
(cid:3)
10