Nonlinear Analysis: Real World Applications 12 (2011) 377–387
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications
journal homepage: www.elsevier.com/locate/nonrwa
Stability and bifurcation in two species predator–prey models
I. Kusbeyzia,c,∗
, O.O. Aybara,c
, A. Hacinliyana,b,c,d
a
Department of Information Systems and Technologies, Yeditepe University, 26 August Campus, Kayisdagi Street, Istanbul, 34755, Turkey
b
Department of Physics, Yeditepe University, 26 August Campus, Kayisdagi Street, Istanbul, 34755, Turkey
c
Department of Mathematics, Gebze Institute of Technology, Kocaeli, Turkey
d
Department of Physics, Bogazici University, Istanbul, Turkey
a r t i c l e i n f o
Article history:
Received 19 April 2010
Accepted 10 June 2010
Keywords:
Hopf bifurcation
Lotka–Volterra model
Normal form
a b s t r a c t
Changes in the number and stability of equilibrium points in the Lotka–Volterra model
as well as some of its generalizations that lead to qualitative changes in the behavior of
the system as a function of some of its parameters are studied by bifurcation analysis.
A generalization involving a cubic interaction as proposed by Nutku is shown to change
the stability properties in a simple way and in particular cases introduce additional
stability. Numerical methods and the approach provided by approximate techniques near
equilibrium points are used in the analysis.
© 2010 Elsevier Ltd. All rights reserved.
1. Introduction
The predator–prey problem attempts to model the relationship in the populations of different species that share the
same environment where some of the species (predators) prey on the others. The prey is assumed to exhibit linear
growth given by a positive parameter. Predator species consume preys with a nonlinear interaction with another set of
parameters that determine the rate of competition between predators. The natural death rate of the predator is assumed
to be linear and given by a negative parameter. One of the earliest implementations, the Lotka–Volterra model serves
as a starting point of more advanced models in the analysis of population dynamics. Because of its unrealistic stability
characteristics [1], stability analysis of the model and its generalizations have recently gained much attention. To understand
the behavior of a nonlinear system one can analyze the existence and stability of equilibrium points. As parameters are
varied changes in the number and stability of equilibrium points lead to bifurcation. Numerical methods are usually
employed to perform this analysis [2]. Approximate techniques near equilibrium points, such as the normal form method [3]
provides a complementary approach for our study. The well-known generalizations of the Lotka–Volterra model include
the addition of polynomial interactions [4], non-monotonic response functions [5], time delayed [6] and diffusion effected,
time delayed [7] non-monotonic interactions. Nutku has proposed a generalization where an additional cubic rather than a
quadratic interaction is involved. In this work, bifurcation properties of quadratic and cubic local and monotonic interactions
are studied. The Nutku generalization is shown to introduce additional stability in a simple way.
2. The classic Lotka–Volterra model and normal forms
The classic mathematical model of predator–prey systems was first developed [8] by Alfred Lotka and Vito Volterra to
model the cyclic changes in the populations of predatory and prey fish in the Adriatic Sea during World War I and also to
∗ Corresponding author at: Department of Information Systems and Technologies, Yeditepe University, 26 August Campus, Kayisdagi Street, Istanbul,
34755, Turkey. Tel.: +90 2165780000; fax: +90 2165780938.
E-mail addresses: ikusbeyzi@yeditepe.edu.tr (I. Kusbeyzi), oaybar@yeditepe.edu.tr (O.O. Aybar), ahacinliyan@yeditepe.edu.tr (A. Hacinliyan).
1468-1218/$ – see front matter © 2010 Elsevier Ltd. All rights reserved.
doi:10.1016/j.nonrwa.2010.06.023
378 I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387
describe chemical reactions in which chemical concentrations oscillate. This model also serves as a first approximation to
more complex models of physical systems such as resonantly coupled laser models [8,9]. The classic Lotka–Volterra model
is given by the following system of differential equations, where dots denote time derivatives.
˙x = ax − bxy,
˙y = −cy + dxy. (1)
Here a, b, c, d > 0 are parameters characterizing the predator–prey environment, x(t) and y(t) are the prey and predator
populations respectively. The equilibrium points are (0, 0) and (c
d
, a
b
). Linearized eigenvalues about the first equilibrium
point (henceforth to be referred to as ‘‘eigenvalues for’’ for brevity) (0, 0) are {a, −c}; they do not ordinarily satisfy resonance
conditions and the origin is a saddle point. On the other hand the eigenvalues for the equilibrium point (c
d
, a
b
) are {±i
√
ac};
these purely imaginary values satisfy the resonance conditions and the system can be expanded into a resonant normal
form. Moving this equilibrium point to the origin gives the following system [5–7]:
˙x = −
bc
d
y − bxy,
˙y =
ad
b
x + dxy. (2)
We can diagonalize the linear part using the transformation x = −ib
d

c
a
(x1 − x2) and y = x1 + x2. This gives the following
system, where the replacement x1 → x, x2 → y has been done following the substitution;
˙x = −i
√
acx −
b
2

1 + i

c
a

(x2
− y2
),
˙y = i
√
acy +
b
2

1 − i

c
a

(x2
− y2
). (3)
These last two steps will not be henceforth shown for brevity.
If we consider x and y as the complex conjugates of one another, this can be rewritten as:
˙z = −i
√
acz −
b
2

1 + i

c
a

(z2
− z
2
) (4)
and its complex conjugate.
Definition 1. The definition and the notation used for normal forms of 2x2 autonomous systems are summarized below.
Let the system of differential equation ˙x = F(x) be given where xϵR2
and FϵR2
→ R2
. To consider the system near its
equilibrium point x = x0 such that F(x0) = 0, this point is moved to the origin by x = x − x0. Taylor expansion near x = 0
gives
˙x = Ax + F2(x, y) + · · · . (5)
Diagonalizing A or bringing it to the Jordan canonical form by a linear transformation produces
˙x = Jx +F2(x, y) + · · · . (6)
Then a sequence of near identity transformation of the form
x = u + F2(u) + · · · (7)
are applied to simplify the system.
In the rest of this work, u and v will refer to the variables in the near identity transformation.
Eigenvalues of the linearized matrix of coefficients around this equilibrium point make the system resonant in every odd
order of the normal form expansion starting from the third order and also yield a Hopf bifurcation. For example,
˙u = i
√
acu

−1 +
b2
6a

1
c
+
1
a

uv +
17b4
432a2

1
c
+
1
a
2
(uv)2

+ · · ·
˙v = i
√
acu

1 −
b2
6a

1
c
+
1
a

uv −
17b4
432a2

1
c
+
1
a
2
(uv)2

+ · · · (8)
is the fifth order normal form.
These equations have the first integral uv = constant and both parameters have periodic solutions that can be
generalized to all orders. Furthermore all the solutions are periodic, causing the predator and prey populations to cycle
I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387 379
and oscillate around this equilibrium. This model is popular as a simple predator–prey system since it shows cycling of
predator–prey populations [10]. However, the model inappropriately describes realistic predator and prey populations
which typically reach characteristic steady state values. Since it does not exhibit structural stability, the simulated values
for the population cycle endlessly without settling down [5–7,11]. Thus, it fails to describe realistic predator and prey
populations which typically reach characteristic steady state values. Realistic models should predict a single closed orbit, or
perhaps finitely many, but not a continuous family of neutrally stable cycles. For this reason, several generalizations of the
model has been proposed by adding additional information that offer more realistic descriptions [12–14].
The integrability of the simple Lotka–Volterra model has been shown [11,15] where elimination of the time gives
dx
dy
=
(a − yb)x
(−c + xd)y
. (9)
The integral is x d + y b − c log x − a log y = constant.
3. Generalized quadratic Lotka–Volterra model with polynomial interaction
The simplest polynomial generalization involves the addition of a quadratic self-coupling suggested by the Verhulst
model. The relative coefficient between the linear and quadratic term can be set to one without any loss of generality by
rescaling the variables, and if necessary, the parameters. In all cases, the quadratic self-coupling does not possess structural
stability, as in the ordinary Lotka–Volterra model. This generalization can introduce one or two stable equilibrium points in
the linearized approximation for appropriate values of the parameters.
There are four different versions of the quadratic self-interaction model, that depends on the relative sign between the
linear and quadratic terms in each equation. By changing the signs of variables and the coefficient of the coupling term,
one can transform the unphysical region of one version to the physical region of another one. Since each form has different
bifurcation properties, we will examine each form and clearly indicate the transformation leading to it from the basic form
given below.
3.1. The case of a supercritical and a subcritical transcritical bifurcation
The basic form is
˙x = ax(1 − x) − bxy,
˙y = −cy(1 − y) + dxy. (10)
The equilibrium points are (0, 0), (1, 0), ((a−b)c
ac−bd
, a(c−d)
ac−bd
) and (0, 1) which have the eigenvalues {a, −c}, {−a, −c + d},
{ac(a−b−c+d)±
√
∆
2(−ac+bd)
} and {a − b, c} in that order where
∆ = ac(ac(a − b − c + d)2
+ 4(a − b)(c − d)(ac − bd)).
The change to a quadratic self-interaction leaves the equilibrium point at the origin unchanged. We no longer have an
unconditional Hopf bifurcation point. This is to be expected since the self-interaction has a form characteristic of transcritical
rather than Hopf bifurcation. The point ((a−b)c
ac−bd
, a(c−d)
ac−bd
) gives Hopf bifurcation under the condition, a < b and c < d. The
system has a stable equilibrium point within physical values if d < c [12].
The system possesses an algebraic integral of the form
xc1
yc2
(1 + c3x + c4y) (11)
where c1, c2, c3, c4 are given expressions involving the equation parameters and provided that ac ̸= bd and
ac(a − b − c + d)
bd − ac
= 0.
3.2. Special case
Using a = c = 1 and b = d = 2 gives:
˙x = x(1 − x) − 2xy,
˙y = −y(1 − y) + 2xy. (12)
This system has four equilibrium points (0, 0), (1, 0), (1
3
, 1
3
) and (0, 1) with the corresponding eigenvalues {±1}, {±1}, {± i√
3
}
and {±1}. The bifurcation diagram is shown in Fig. 1. The first two and the last equilibrium points are saddle points and they
are unstable. Nevertheless the third equilibrium point has pure imaginary eigenvalues satisfying the resonance condition,
380 I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387
Fig. 1. Family of limit cycles for system (12).
the system can be expanded into a resonant normal form. Moving the third equilibrium point to the origin and transforming
the new system into complex form, one has:
˙z = i
z
√
3
− i
√
3z
2
(13)
and its complex conjugate.
Normal form expansion to the fifth order gives:
˙u =
iu
√
3
(1 − 6uv − 24(uv)2
),
˙v =
iv
√
3
(−1 + 6uv + 24(uv)2
). (14)
Except for the difference in coefficients this normal form is exactly the same as that for the ordinary Lotka–Volterra case
and possesses the identical integral uv = constant [12–14].
When each of the parameters is allowed to become variable one by one, and a bifurcation analysis is carried out,
the system exhibits a subcritical, a supercritical Hopf bifurcation and a transcritical bifurcation for each parameter.
The supercritical Hopf bifurcation point (1
3
, 1
3
) has the eigenvalues {± i√
3
} and leads to a limit point cycle with double
characteristic exponents µ = 1 and nonzero normal form coefficients, hence intermittency is observed. A family of
limit cycles surround the supercritical Hopf equilibrium point (1
3
, 1
3
); their period and extent continues to grow until the
transcritical bifurcation point at (0, 1) is reached. All of the limit cycles lie within a right triangle whose vertices are the
equilibrium points (0, 0), (1, 0) and the transcritical bifurcation point. The median between the equilibrium points that do
not show transcritical bifurcation leads to the subcritical Hopf point at (0.5, 0).
This system admits as a special case of (11), the first integral xy(−1 + x + y) = constant and the triangle that forms the
boundary consists of the lines x = 0, y = 0 and x + y = −1 that forms the degenerate cases of the first integral with the
constant set to zero [13,14,11].
3.3. The case of two supercritical transcritical bifurcations
Unphysical regions of the model given by [16] can be mapped to physical regions of two models where either the predator
or prey self-coupling acquires the positive sign. This will impose the same asymptotic behavior for both self-couplings.
However the modified systems lose Hopf bifurcation property in the physical region.
Changing the signs of y and b yields the system
˙x = ax(1 − x) − bxy,
˙y = −cy(1 + y) + dxy. (15)
The equilibrium points are (0, 0), (1, 0), ((a+b)c
ac+bd
, a(−c+d)
ac+bd
) and (0, −1). Eigenvalues associated with the equilibrium point
(0, 0) are {a, −c} so that this is a saddle node. Eigenvalues associated with the equilibrium point (1, 0) are {−a, d − c}.
If d > c this equilibrium is a saddle node, otherwise it is stable. The eigenvalues associated with the equilibrium point
(c(a+b)
ac+bd
, a(d−c)
ac+bd
) form a complex conjugate pair and the corresponding eigenvectors are complex conjugates of one another.
If
ac(a+b−c+d)
−2(ac+bd)
< 0 then this equilibrium is stable. Eigenvalues associated with the equilibrium point (0, −1) are {a + b, c}.
I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387 381
This equilibrium is unstable. Since the stability of the equilibrium points can change with the parameter values, bifurcation
may, in principle, occur [8,15–17].
Lemma 1. For the system (15), the condition for the existence of an algebraic integral in the form xy(−1 + x + y) = constant
coincides with the Hopf bifurcation condition.
Proof. The required condition for Hopf bifurcation to occur is ac(a + b − c + d) = 0. The system possesses an algebraic
integral of the form xc1 yc2 (1 + c3x + c4y) provided that ac ̸= −bd and
ac(−a − b + c − d)
ac + bd
= 0
which coincides with the required condition for the third equilibrium point to give Hopf bifurcation. Otherwise only
transcritical bifurcation is observed in the physical region.
3.4. The case of two subcritical transcritical bifurcations
On the other hand changing the signs of x and the parameter d yields the system
˙x = ax(1 + x) − bxy,
˙y = −cy(1 − y) + dxy (16)
where this time there is a stable equilibrium point at (−1, 0) with the eigenvalues {−a, −c − d}. The origin is again a saddle
point with the persistent eigenvalues {a, −c}. The remaining two equilibrium points are (c(b−a)
ac+bd
, a(c+d)
ac+bd
) and (0, 1) which
have the eigenvalues {ac(a−b−c−d)±
√
∆
2(−ac−bd)
} and {a − b, c} respectively where
∆ = ac(ac(−a + b + c + d)2
+ 4(a − b)(c + d)(ac + bd)).
Although the third equilibrium point seems to promise a Hopf bifurcation, the required condition
c(c+d)(b+c+d)
b+c
< 0 indicates
that this is impossible for physical parameter values [8,15–17].
3.5. The case of a subcritical and a supercritical transcritical bifurcation
Combining both transformations (changing the signs of x, y, b and d) yields the system
˙x = ax(1 + x) − bxy,
˙y = −cy(1 + y) + dxy (17)
which has an identical stability analysis with a slight modification in the denominator of the eigenvalue of the third
equilibrium point and in the eigenvalues of the last equilibrium point as compared to the previous system. The condition
for Hopf bifurcation is also different, namely,
c(c+d)(c−b+d)
c−b
< 0. It can occur for physically acceptable values of the
parameters [15–17].
4. Cubic model with xk
y interaction
Nutku has suggested a model that includes a cubic rather than quadratic self-interaction. Unlike the classical
Lotka–Volterra system or the system with quadratic coupling [8,15–17], the cubic self-coupling can lead to a structurally
stable system if the cubic terms have the negative sign [18–20]. If one insists on structural semi-stability, the interaction
coupling should be of the form xy in both equations or xk
y in the prey and xyk
in the predator equation where k is a positive
integer and k ≤ 2. A cubic interaction of this form has not been studied before. The possibility of xy2
in the prey and x2
y can
be reduced to the quadratic coupling case with the transformation u = x2
v = y2
.
It is always possible to set the relative coefficient in the self-coupling term to 1 by rescaling the variables and to change
the relative sign between the linear and cubic term in either equation. Details will be given in Lemma 4 for the xy coupling,
other couplings can also be treated in a similar way.
The suggested system is:
˙x = ax(1 − x2
) − bxk
y,
˙y = −cy(1 + y2
) + dxyk
. (18)
Its generalization would involve self-couplings of the form ax(1 ± x2
) and −cy(1 ∓ y2
). To find the equilibrium points, we
need to solve the equations x(a(1 ± x2
) − bxk−1
y) = 0 and y(c(1 ∓ y2
) − xdyk−1
) = 0.
Lemma 2. From the factorization of the equations above, it is clear that for all appropriate values of k, five of the roots are of the
form (0, 0), (0 ± 1) or (0, ±i), (±1, 0) or (±i, 0).
382 I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387
Proof. Setting the factor x = 0 into each of the two factors in the y equation and the factor y = 0 into each of the two factors
in the x equation. The remaining roots are obtained by solving the two nonzero factors simultaneously: a(1 ± x2
) = bxk−1
y
and c(1 ∓ y2
) = bxyk−1
. In spite of the fact that these equations depend on k, it will be seen that for k ≤ 2, the solution will
involve a quartic equation in a single variable, so that there will be nine equilibrium points.
Returning to system (18), the Jacobian matrix is
[
a − 3ax2
− bkyxk−1
−bxk
dyk
−c − 3cy2
+ dkxyk−1
]
.
The simple equilibrium points are (0, 0), (0, ±i), (±1, 0). If k = 1, the eigenvalues are {−a, −c}, {a ± ib, 2c} and
{−2a, −c ± d} and if k = 2, eigenvalues are {−a, −c}, {a, 2c} and {−2a, −c} in the given order.
The remaining roots for k = 1 and k = 2 will now be analyzed. The quartic equation for the remaining four equilibrium
points for k = 1 is
(1 − x2
)2
=
b2
a2

xd
c
− 1

(19)
which leads to the following cases.
Lemma 3. The quartic equation above can have two or zero positive real roots.
Proof. Since the left-hand side is (1−x2
)2
, its intersection with the right-hand side b2
a2 (xd
c
−1) cannot lead to negative roots
for physical values of the parameters. Hence, two of the roots are complex. As far as the remaining roots are concerned, the
equation can be rewritten as:
x4
− 2x2
−
b2
a2
xd
c
+
b2
a2
+ 1 = 0 (20)
and we have the following cases:
i. two real and positive roots if c < d,
ii. two complex roots if d < c,
iii. If c = d, the equation is
(1 − x2
)2
=
b2
a2
(x − 1)
and this gives a degenerate case and a root of 1.
We turn to the system for k = 2. The equilibrium points are (0, 0), (0, ±i), (±1, 0) and four nontrivial points. The
remaining four equilibrium points satisfy the following quadratic equation for t = x2
.
bc
ad
t +
ca
bd
(1 − t2
) = t(1 − t). (21)
The roots of this equation are
x1,2 = ±

−c + 2a2c + abd −

c2 − 4(ac)2 − 2abcd + (abd)2
2a(ac + bd)
x3,4 = ±

−c + 2a2c + abd +

c2 − 4(ac)2 − 2abcd + (abd)2
2a(ac + bd)
.
The cubic self-coupling is of interest since the system undergoes pitchfork bifurcation rather than transcritical
bifurcation [17–19] irrespective of the interaction coupling. Although further self-couplings retain similar results in the
linearized eigenvalue expansions they are not of interest for structural stability reasons.
We now give examples of the cases mentioned in Lemma 3.
4.1. Case I
According to the choice of the parameters in the predator equation the generalized cubic system has two real roots one
between zero and one, the other greater than one.
In this case the generalized Lotka–Volterra system has additionally two complex equilibrium points.
I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387 383
Fig. 2. Bifurcation diagram for system (22).
Example. As an example, we use a = c = 0.5 and b = d = 1 so that we have the system
˙x = 0.5x(1 − x2
) − xy,
˙y = −0.5y(1 + y2
) + xy (22)
which has nine equilibrium points (0, 0), (0, ±i), (±1, 0), (−1.35 ± 1.52i, 0.73 ± 2.06i) and the two real positive roots
(0.55, 0.34), (2.15, −1.81), the quartic has four complex roots. Varying the predator’s growth rate gives the bifurcation
diagram in Fig. 2. There are two pitchfork bifurcation points (BP) each of which bifurcate one new branch. In addition the
main bifurcation curve gives rise to two saddle node bifurcation points (LP). Each branch of the bifurcation point with x ̸= y
illustrates the phenomenon of symmetry breaking. Where each of the two branches meet, four cusp bifurcation points (CP)
are observed signaled by a zero eigenvalue and vanishing quadratic coefficient for the saddle node bifurcation [15,17–19].
4.2. Case II
Using, a = c = 1.5 and b = d = 1 we have the system
˙x = 1.5x(1 − x2
) − xy,
˙y = −1.5y(1 + y2
) + xy (23)
which has nine equilibrium points (±1, 0), (0, 0), (0, ±i), (−1.049 ± 0.415i, 0.106 ± 1.308i), (1.049 ± 0.177i, −0.105 ∓
0.558i). A local bifurcation analysis around these equilibrium points gives no bifurcation points but only a closed orbit.
4.3. Case III (degenerate case)
This is the degenerate case where x = 1 [21].
˙x = x(1 − x2
) − xy,
˙y = −y(1 + y2
) + xy. (24)
Equilibrium points are (0, 0), (0, ±i), (±1, 0), (1, 0), (1.205, −0.453), (−1.102 + 0.665i, 0.226 + 1.467i), (−1.102 −
0.665i, 0.226 − 1.467i). A local bifurcation analysis similar to Case I is obtained for the equilibrium point (1.205, −0.453).
5. The generalized cubic Lotka–Volterra model with simple polynomial interaction
Lemma 4. The following systems of generalized cubic Lotka–Volterra models with simple polynomial interaction can be
transformed into each other.
Proof. Let us consider the system
˙x = ax(1 − x2
) − bxy,
˙y = −cy(1 + y2
) + dxy. (25)
384 I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387
Using the transformation x = ix and d = −id and dividing the first equation by i gives the first varied version of the system
˙x = ax(1 + x2
) − bxy,
˙y = −cy(1 + y2
) + dxy. (26)
Using the transformation of y = iy and b = −ib and dividing the second equation by i we have the second varied version
which is
˙x = ax(1 − x2
) − bxy,
˙y = −cy(1 − y2
) + dxy. (27)
The third varied version is obtained by carrying out both transformations and divisions mentioned above.
˙x = ax(1 − x2
) − bxy,
˙y = −cy(1 − y2
) + dxy. (28)
The original system and its three varied forms mentioned above cover all possibilities. Finally it is possible to change the
equilibrium points by using the transformations x = xk, y = yk, b = b
k
, d = d
k
and dividing both equations by k
˙x = ax(1 − (kx)2
) − bxy,
˙y = −cy(1 + (ky)2
) + dxy. (29)
Hence the unstable systems can be transformed into stable systems via variable transformations [13,14,21].
Nevertheless we will investigate each system in its own physical region and give the associated bifurcation schemes.
5.1. The case of two supercritical pitchfork bifurcations
We use the system given by
˙x = ax(1 − x2
) − bxy,
˙y = −cy(1 + y2
) + dxy. (30)
The equilibrium point (0, 0) has the eigenvalues {a, −c} and again it is a saddle point. The equilibrium points (0, −i) and
(0, i) have the eigenvalues {a + bi, 2c} and {a − bi, 2c} respectively and they are both unstable equilibrium points. The
equilibrium points (1, 0) and (−1, 0) have the eigenvalues {−2a, −c − d} and {−2a, −c + d} respectively where one of
them is a stable equilibrium point and the other one is a stable equilibrium point under the condition d < c which falls into
case II. Considering these equilibrium points we can conclude that for positive parameter values a resonant normal form
does not exist around the stable equilibrium points. With the help of the Lyapunov function V(x, y) = x2
+ y2
we have
˙V = 2ax2
− 2ax4
− 2bx2
y − 2cy2
− 2cy4
+ 2dxy2
so the system is shown to be stable and structurally stable for a > 0 and
c > 0. This is obvious if one considers that the dominant terms are
˙x = −ax3
,
˙y = −cy3
(31)
with the obvious solution x2
(t) = x2(t=0)
1+2atx2(t=0)
and y2
(t) = y2(t=0)
1+2cy2(t=0)
.
5.2. The case of two subcritical pitchfork bifurcations
We use the system given by
˙x = ax(1 + x2
) − bxy,
˙y = −cy(1 − y2
) + dxy. (32)
The equilibrium points are (0, 0), (0, 1), (0, −1) and six nontrivial points. The origin is hyperbolic saddle. (0, 1) may be either
saddle or unstable point. (0, −1) is an unstable point.
5.3. The case of a subcritical and a supercritical pitchfork bifurcation
We use the system given by
˙x = ax(1 + x2
) − bxy,
˙y = −cy(1 + y2
) + dxy. (33)
The equilibrium points are (0, 0), (0, ±i), (±i, 0) and four nontrivial points. In this case locally stable equilibrium in the
parameter plane is guaranteed.
I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387 385
Fig. 3. Bifurcation diagram showing the family of limit cycles around the Hopf point of system (35).
5.4. The case of a supercritical and a subcritical pitchfork bifurcation
The generic form describing the case of supercritical and subcritical bifurcation is given by the following system
˙x = ax(1 − x2
) − bxy,
˙y = −cy(1 − y2
) + dxy. (34)
For convenience we study the following special case of the system
˙x = x(1 − x2
) − 2xy,
˙y = −y(1 − y2
) + 2xy. (35)
The system has the equilibrium points (0, 0), (0, 1), (1, 0) and (−1 −
√
2, −1 −
√
2) as hyperbolic saddles, (0, −1) as an
unstable point with the eigenvalues {2, 3}, two complex conjugate points with cross conjugate eigenvalues, a stable point
at (−1, 0) with the eigenvalues {−2, −3} and the point (
√
2 − 1,
√
2 − 1) with a pair of purely imaginary eigenvalues as a
Hopf point. The third order normal form about this equilibrium point is, as expected
˙u = iu(−27.8611918118930381uv + 313.2622478846395944)
˙v = iv(27.8611918118930381uv − 313.2622478846395944)
higher order terms are omitted for brevity. The reason for the coefficients with large absolute values is the fact that a small
denominator 5
√
2 −7 appears. Since the exact expressions are very complicated the expressions are given to sixteen places
of decimals. This gives supercritical Hopf bifurcation. There is also observed a subcritical Hopf bifurcation at (0.5, 0) and two
symmetric pitchfork bifurcation points at (0, 1) and (0, −1). The bifurcation diagram is shown in Fig. 3. The Hopf point leads
to intermittency [8,15–17,22,23].
6. Cubic model with x2
y interaction
The equilibrium points and bifurcation scenarios for this interaction are summarized below.
6.1. The case of two supercritical pitchfork bifurcations
We use the system given by
˙x = ax(1 − x2
) − bx2
y,
˙y = −cy(1 + y2
) + dxy2
. (36)
The equilibrium points are (0, 0), (0, ±i), (±1, 0) and four nontrivial points.
6.2. The case of two subcritical pitchfork bifurcations
We use the system given by
˙x = ax(1 + x2
) − bx2
y,
˙y = −cy(1 − y2
) + dxy2
. (37)
The equilibrium points are (0, 0), (0, ±1), (±i, 0) and four nontrivial points.
386 I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387
6.3. The case of a subcritical and a supercritical pitchfork bifurcation
We use the system given by
˙x = ax(1 + x2
) − bx2
y,
˙y = −cy(1 + y2
) + dxy2
. (38)
The equilibrium points are (0, 0), (0, ±i), (±i, 0) and four nontrivial points.
6.4. The case of a supercritical and a subcritical pitchfork bifurcation
We use the system given by
˙x = ax(1 − x2
) − bx2
y,
˙y = −cy(1 − y2
) + dxy2
. (39)
The equilibrium points are (0, 0), (0, ±1), (±1, 0) and four nontrivial points.
In all cases the eigenvalues are {a, 2c} for simple equilibrium points except the origin which is saddle [22,23].
Due to the problems of stability we do not go further than the interaction x2
y.
7. Conclusion
In this work, the stability and bifurcation properties of two species Lotka–Volterra systems were studied. Both the regular
Lotka–Volterra model and its generalization via a quadratic self-coupling show transcritical and Hopf bifurcation. In both
cases, bounded orbits exist under appropriate conditions but the system is not stable in the Lyapunov sense [8,15,17,22].
Although approaches such as non-polynomial or time delay couplings have been used, an alternative way to overcome
the inherent instability of the two dimensional Lotka–Volterra system with a quadratic self-interaction is simply replacing
the quadratic self-coupling by a cubic one as suggested by Nutku. A further generalization involving couplings of the form xk
y
in the prey and xyk
has also been considered. Normal form expansions near equilibrium points leading to Hopf bifurcation
were given. In spite of its simplicity [8,15–17], the Nutku generalization provides a rich spectrum of equilibrium points
leading to Hopf, pitchfork, saddle node and cusp bifurcations. These systems are Lyapunov stable for proper choice of the
sign of the cubic term. Thus, the simplicity of the Nutku suggestion that keeps the interaction at a simple yet stable form is
demonstrated [20,22,24,23].
Acknowledgement
The authors thank Yavuz Nutku for suggesting this problem.
References
[1] I.M. Gleria, A. Figueiredo, T.M. Rocha Filho, Stability properties of a general class of nonlinear dynamical systems, Journal of Physics A: Mathematical
and General 34 (2001) 3561–3575.
[2] D. Ghosh, A.R. Chowdhury, On the bifurcation pattern and normal form in a modified predator–prey nonlinear system, Journal of Computational and
Nonlinear Dynamics 2 (2007) 267–273.
[3] P. Yu, G. Chen, The simplest parametrized normal forms of Hopf and generalized Hopf bifurcations, Nonlinear Dynamics 50 (2007) 297–313.
[4] H. Zhu, S.A. Campbell, G.S.K. Wolkowicz, Bifurcation analysis of a predator–prey system with nonmonotonic functional response, SIAM Journal on
Applied Mathematics 2 (63) (2002) 636–682.
[5] H.W. Broer, V. Naudot, R. Roussarie, K. Saleh, Bifurcations of a predator–prey model with nonmonotonic response function, Comptes Rendus
Mathematique 341 (2005) 601–604.
[6] D. Xiao, Multiple bifurcations in a delayed predator–prey system with nonmonotonic functional response, Journal of Differential Equations 176 (2001)
494–510.
[7] X.P. Yan, W.T. Li, Hopf bifurcation and global periodic solutions in a delayed predator–prey system, Applied Mathematics and Computation 177 (2006)
427–445.
[8] A. Lotka, Elements of Physical Biology, 1st edition, Williams and Wilkins Company, Baltimore, 1925.
[9] A. Hacinliyan, I. Kusbeyzi, O.O. Aybar, Approximate solutions of Maxwell Bloch equations and possible Lotka Volterra type behavior, Nonlinear
Dynamics (2010) doi:10.1007/s11071-010-9695-5.
[10] V. Pankovic, D. Banjac, R. Glavatovic, M. Predojevic, A simple solution of the Lotka–Volterra equations, preprint.
[11] H.T. Davis, Introduction to Nonlinear Differential Equations and Integral Equations, 1st edition, Dover, New York, 1960.
[12] Hernández-Bermejo, V. Fairén, Lotka–Volterra representation of general nonlinear systems, Mathematical Biosciences 1 (140) (1997) 1–32.
[13] X. Huang, Y. Wang, A. Cheng, Limit cycles in a cubic predator–prey differential system, Journal of the Korean Mathematical Society 43 (4) (2006)
829–843.
[14] K. Yan, X. Huang, Hopf bifurcation in a generalized predator–prey model, International Journal of Applied Science and Engineering 4 (3) (2006)
253–258.
[15] C. Christopher, Normalizable, integrable and linearizable saddle points in the Lotka–Volterra system, Qualitative Theory Of Dynamical Systems 5
(2004) 11–61.
[16] R.H. Enns, G.C. McGuire, Nonlinear Physics with Mathematica for Scientists and Engineers, 1st edition, Birkhauser, Boston, 2001.
[17] E.C. Zeeman, M.L. Zeeman, An n-dimensional competitive Lotka–Volterra system is generically determined by the edges of its carrying simplex,
Nonlinearity 15 (2002) 2019–2032.
I. Kusbeyzi et al. / Nonlinear Analysis: Real World Applications 12 (2011) 377–387 387
[18] X.K. Sun, H.F. Huo, H. Xiang, Bifurcation and stability analysis in predator–prey model with a stage-structure for predator, Nonlinear Dynamics 58 (3)
(2009) 497–513.
[19] C. Liu, Q. Zhang, X. Zhang, X. Duan, Dynamical behavior in a stage-structured differential algebraic prey–predator model with discrete time delay and
harvesting, Journal of Computational and Applied Mathematics 231 (2009) 612–625.
[20] Y. Nutku, Hamiltonian structure of the Lotka–Volterra equations, Physics Letters A I 145 (1990) 27–28.
[21] X. Huang, Y. Wang, L. Zhu, One and three limit cycles in a cubic predator–prey system, Mathematical Methods in the Applied Sciences 30 (2007)
501–511.
[22] I. Kusbeyzi, A. Hacinliyan, Bifurcation scenarios of some modified predator–prey nonlinear systems, Journal of Applied Functional Analysis 4 (3) (2009)
519–527.
[23] L. Gardini, R. Lupini, M.G. Messia, Hopf bifurcation and transition to chaos in Lotka–Volterra equation, Journal of Mathematical Biology 27 (3) (1989)
259–272.
[24] J.A. Vano, J.C. Wildenberg, M.B. Anderson, J.K. Noel, J.C. Sprott, Chaos in low-dimensional Lotka–Volterra models of competition, Nonlinearity 19 (2006)
2391–2404.