www.biomathforum.org/biomath/index.php/biomath

ORIGINAL ARTICLE

Effects of Discrete Time Delays and
Parameters Variation on Dynamical Systems

Ibrahim Diakite∗, Benito M. Chen-Charpentier†
∗Department of Global Health and Social Medicine

Harvard Medical School, Boston, USA
ibrahim diakite@hms.harvard.edu
†Department of Mathematics

University of Texas Arlington, Arlington, USA
bmchen@uta.edu

Received: 29 August 2014, accepted: 20 May 2015, published: 8 June 2015

Abstract—Delay Differential Equations (DDE’s)
have received considerable attention in recent years.
While most of these articles focused on the effects
of the time delays on the stability of the equilib-
rium points and on the bifurcation that they may
raised, very few papers address the key roles that
system parameters play on if and how the discrete
delays induce stability changes of the equilibria
and produce bifurcations near such equilibria. In
this article we focus on that question in a general
setting, that is, if one has a system of DDE’s with
one or multiple discrete time delays, what are the
results of changing the system parameters values
on the effects of the discrete time delays on the
dynamic of the system. We present general results
for one equation with one and two delays and study
a specific example of one equation with one delay.
We then establish the procedure for n equations with
multiple delays and do a specific example for two
equations with two delays. We compute the steady
states and analyze their stability as both chosen
bifurcation parameters, the discrete time delay τ
and a local equation parameter µ, cross critical
values. Our analysis shows that while changes in
both parameters can destabilize the steady state, the
discrete time delay can only cause stability switches
of the steady state for certain values of µ, while

the effects of the local equation parameter on the
steady state do not necessarily depend on the value
of τ. While µ may cause the system to go through
different type of bifurcations, the discrete time delay
can only introduce a Hopf bifurcation for certain
values of µ.

Keywords-delay differential equations; bifurca-
tion; predator-prey.

I. INTRODUCTION

It is well known, that the values of the
parameters play a crucial role in the behav-
ior of dynamical systems and that changes in
the values can change the behavior significantly.
It has also been shown by many researchers
(Perelson[1],Allen[2],Bellen[3]) that there is a
need to incorporate discrete time delays in dy-
namical systems (biological systems, physical sys-
tems,...) as studied.

Models that incorporate such delays are referred
to as delay differential equations (DDE’s). DDE’s
have been extensively studied by many researchers
including pioneers Bellman[4], Driver [5], and
in more recent years by Culshaw[6], Gakkhar[7],
Bellen[3], and a superb monograph on the subject

Citation: Ibrahim Diakite, Benito M. Chen-Charpentier, Effects of Discrete Time Delays and
Parameters Variation on Dynamical Systems, Biomath 1 (2015), 1505201,
http://dx.doi.org/10.11145/j.biomath.2015.05.201

Page 1 of 17

http://www.biomathforum.org/biomath/index.php/biomath
http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

by Gopalsamy [8]. While most of these research
papers focus on issue of the stability changes
caused by the delay(s), the main motivation of this
paper is to study how a local bifurcation parameter
of the system may affect the changes in stability
caused by the delay (s).

Published papers have shown that the incorpora-
tion of discrete time delays can highly impact the
dynamics of the system, since they can switch the
stability of a steady state point, and can also cause
the system to go through a Hopf bifurcation near
that steady state point (Culshaw[6], Gakkhar[7],
Bellen[3]). In this paper we consider a system of
n delay differential equations (DDE’s) with one
parameter µ as the bifurcation parameter and also
with one or more discrete time delays, τ, which
can also behave as bifurcation parameters. We are
interested in investigating how the parameters µ
and τ affect the stability of the steady state points
of the system, and, more important, how their
effects on the system are correlated to each other.
We present general results in the one dimensional
case (propositions 1 to 3) for necessary and suffi-
cient conditions for a stability switch and present a
specific example to illustrate these conditions. For
the n dimensional case (n ≥ 2) we establish the
main ideas, but since there are multiple possible
cases, we consider only a specific example. We
present a non-Kolmogorov type of predator-prey
model similar to the model presented by Ruan [9].
In this model we introduce two delays, τ1 > 0 and
τ2 > 0, to represent the time lag in the growth
to maturity of the prey, and the time lag in the
growth to maturity of the predator, respectively.
We show how the dynamics of the system change
depending on certain conditions on τ1 and on
another bifurcation parameter R. We also point out
conditions for the system to go through stability
changes when both delays τ1 and τ2 are non-zero.
We present necessary conditions for the system
to go through a Hopf bifurcation for τ1 > 0
and τ2 = 0. Finally we show numerical results
illustrating the theoretical results.

II. ONE DIMENSIONAL FIELD

A. One Equation with One Delay

Consider the one dimensional delay differential
equation with the time delay τ, and the parameter
µ as bifurcation parameters:

dX

dt
= f(X(t),X(t− τ),µ), (1)

where f is assumed to be smooth enough to guar-
antee the existence and uniqueness of solutions to
(1) under the initial condition (R. Bellman and K.
L. Cooke [4])

X(θ) = φ(θ), θ ∈ [−τ, 0].

Unfortunately equation (1) is too general to ana-
lyze. Therefore we will consider a more special
form:

dX

dt
= f1(X(t),µ) + f2(X(t− τ),µ). (2)

This form has the advantage that it simplifies the
analytical work and also it is the form present
in many population dynamical models involving
delays [6], [7], [9], [10]. The DDE (2) may or
may not have equilibrium points (or steady states)
and these will depend on the values of µ. Let
µ∗ ∈ Dµ = {µ ∈ R : f(X∗,X∗,µ) = 0 exists}
,that is µ∗ is in the range of values of µ for
which the DDE has an equilibrium point X∗,
i.e., f(X∗,X∗,µ∗) = 0. We are interested in
studying the stability of such equilibrium point. In
particular, in studying the effect of the parameter
µ and of the discrete time τ on its stability. To do
this we linearize the DDE around the equilibrium
point. The characteristic equation is :

λ−
df1
dX
|(X∗,µ∗) −

df2
dX
|(X∗,µ∗)e

−λτ = 0, (3)

and the stability of the equilibrium point (X∗,µ∗)
is determined by the sign of the real part of the
eigenvalues λ of equation (3).

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 2 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

1) Stability of the Steady State: If τ = 0 then
the characteristic equation (3) becomes

λ−
df1
dX
|(X∗,µ∗) −

df2
dX
|(X∗,µ∗)=0.

The stability of the steady state then depends only
on values of µ∗ within Dµ. We have two cases:
(a) The steady state (X∗,µ∗) is stable if

df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗) < 0.

(b) The steady state (X∗,µ∗) is unstable if
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗) > 0.

Assume that condition (a) holds, namely the
steady state (X∗,µ∗) is stable when there is no
delay (τ = 0). We want to know if there exists
τ > 0 for which the steady state will lose stability.
So for τ ≥ 0, let λ(τ) = α(τ) + iω(τ). The
characteristic equation (3) becomes:

α + iω =
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗)e

−ατcosωτ+

i
df2
dX
|(X∗,µ∗)e

−ατsinωτ,

(4)

where, for clarity in the notation, we have not
explicitly shown the dependence on τ. Separating
the real and imaginary parts, we have:

α =
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗)e

−ατcosωτ, (5)

ω =
df2
dX
|(X∗,µ∗)e

−ατsinωτ. (6)

The steady state will lose stability when the real
part of the eigenvalue λ crosses the zero axis from
negative to positive as τ passes a critical value.
By Rouche’s Theorem (Dieudonne[11], Theorem
9.17.4) and by the continuity in τ, the transcenden-
tal equation (3) has roots with positive real parts if
and only if it has pure imaginary roots. Therefore,
we look at when the real part of the eigenvalue λ
becomes zero. In other words, we want to find if
there exists a τc > 0 such that α(τc) = 0. Since

α(0) =
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗),

and α(0) < 0 by assumption (a) , therefore if τc >
0 exists such that α(τc) = 0 then by the continuity

(Michael Y. Li and Hogying Shu [10]) of α we
have:
• α(τ) < 0 for any 0 ≤ τ < τc,
• α(τ) > 0 for any τ > τc.

Namely the steady state (X∗,µ∗) will lose stability
as the delay parameter τ crosses a critical value
τc. Such τc exists if and only if α(τc) = 0 and
ω(τc) = ωc satisfies :

df1
dX
|(X∗,µ∗) = −

df2
dX
|(X∗,µ∗)cosωcτc (7)

ωc =
df2
dX
|(X∗,µ∗)sinωcτc. (8)

Squaring equations (7) and (8), and adding them
up, we obtain:

ω2c = [
df2
dX
|(X∗,µ∗)]

2 − [
df1
dX
|(X∗,µ∗)]

2. (9)

If equation (9) has at least a positive root ωc, then
there exists a τc > 0 such that α(τ) > 0 whenever
τ > τc (see proof in Appendix A). An important
question we want to address is, since equation (9)
depends on the bifurcation parameter µ∗, can one
chose µ∗ within Dµ so that equation (9) does not
have a positive root ωc? That is, are there values
of µ∗ within Dµ such that the delay does not
have any effect on the stability of the steady state
(X∗,µ∗)? This question motivates the following
propositions (see Appendix B for the proof).

Proposition 1: Consider the one dimensional
delay differential equation

dX

dt
= f1(X(t),µ) + f2(X(t− τ),µ).

And assume that the steady state (X∗,µ∗) is
stable for τ = 0 then we have
(i) If df1

dX
|(X∗,µ∗) < 0 and

df2
dX
|(X∗,µ∗) > 0 then

the steady state (X∗,µ∗) remains stable for
all τ ≥ 0.

(ii) If df1
dX
|(X∗,µ∗) > 0 and

df2
dX
|(X∗,µ∗) < 0 then

there exists a critical value of the delay such
that the steady state loses stability as the delay
crosses its critical value.

(iii) If df1
dX
|(X∗,µ∗) < 0 and

df2
dX
|(X∗,µ∗) < 0 then:

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 3 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

(a) the steady state remains stable for all
τ ≥ 0 if

|
df2
dX
|(X∗,µ∗)| < |

df1
dX
|(X∗,µ∗)|,

(b) there exists a τc > 0 such that the steady
state becomes unstable for all τ > τc if

|
df2
dX
|(X∗,µ∗)| > |

df1
dX
|(X∗,µ∗)|.

Proposition 2: Consider the one dimensional
delay differential equation

dX

dt
= f1(X(t),µ) + f2(X(t− τ),µ).

And assume that the steady state (X∗,µ∗) is stable
for τ = 0, that conditions of proposition 1(iii)
hold, and that further more for some µ∗ within
Dµ we have:

df1
dX
|(X∗,µ∗) = g(X

∗)O(µ∗),
df2
dX
|(X∗,µ∗) = h(X

∗)O(
1

µ∗
),

then there exists a critical value for µ∗ within Dµ
such that the steady state (X∗,µ∗) will stay stable
for all τ ≥ 0 when µ∗ > µc.

2) Example: Consider the one dimensional
DDE{

dY
dt

= µ
Y (t)
Y (t)+1

− 1
µ
Y (t− τ)2, if µ 6= 0

Y (t) = 0, if µ = 0

where µ is a bifurcation parameter and τ ≥ 0
is a discrete time delay. For µ ∈ Dµ = R, the
equation has two non-negative equilibrium points:
the trivial one Y ∗0 = 0, and the positive equilib-
rium point Y ∗1 =

−1+
√
1+4µ2

2
. The characteristic

equation is given as

λ−µ
1

(Y ∗ + 1)2
−

2

µ
Y ∗e−λτ = 0. (10)

• For the trivial equilibrium point Y ∗ = 0, its
stability only depends on µ since equation
(10) evaluated at Y ∗ = 0 becomes λ = µ.
The trivial equilibrium is unstable for µ > 0
and all τ ≥ 0.

The trivial equilibrium is stable for µ < 0
and all τ ≥ 0.

• At Y ∗1 =
−1+

√
1+4µ2

2
, equation (10) becomes:

λ− 4µ
(1+
√
1+4µ2)2

+
√
1+4µ2−1

µ
e−λτ = 0,

(11)

then the stability of Y ∗1 depends on both µ
and τ.

1) If τ = 0 then equation (11)
becomes

λ = −
4µ
√

1 + 4µ2

(1 +
√

1 + 4µ2)2

then
λ < 0 if µ > 0, therefore the equilibrium Y ∗1
is stable (Fig 2)
λ > 0 if µ < 0, therefore the equilibrium Y ∗1
is unstable (Fig 2).

Remark: To better understand the situation,
the stability of both equilibria when there is
no delay is shown in the following table:

TABLE I
STABILITY REGIONS

Case Y ∗0 = 0 Y
∗
1 =

−1+
√

1+4µ2

2

µ < 0 stable unstable
µ = 0 stable stable
µ > 0 unstable stable

At the equilibrium (Y,µ) = (0, 0), there
is an exchange of stability. This is a
transcritical bifurcation (Guckenheimer[12]).
Geometrically, there are two curves of
equilibria which intersect at the origin and
lie on both sides of µ = 0. Stability of the
equilibrium changes along either curve on
passing through µ = 0.

2) If τ > 0 and λ(τ) = α(τ) + iω(τ),
there exists a critical τc such that α(τc) = 0
and λ(τc) = ±iω(τc) = ±iωc (a pair of pure

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 4 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

Fig. 1. Transcritical bifurcation around µ = 0. Unstable
equilibrium (red), and Stable equilibrium (blue).

imaginary eigenvalues) is solution of (8) if
and only if

ω2c =
16[µ2(1+

√
1+4µ2)2−1]

(1+
√
1+4µ2)4

has a positive root
ωc,

and that is the case if and only if µ2(1 +√
1 + 4µ2)2 − 1 > 0.

(2a) If µ ≥ 1
2

then µ2(1+
√

1 + 4µ2)2−1 > 0
therefore there exists τc > 0 such that
the equilibrium loses stability whenever
τ > τc (Fig 3 left).

(2b) If µ ≤−1
2

then µ2(1+
√

1 + 4µ2)2−1 >
0 therefore there exists τc > 0 such that
the equilibrium gains stability whenever
τ > τc .

(2c) If −1
2

< µ < 1
2

then µ2(1 +√
1 + 4µ2)2 − 1 < 0 therefore the de-

lay has no effect on the stability of the
equilibrium.

For µ ≥ 1
2
, the equilibrium is unstable for all τ >

0.55, and for 0 < µ < 1
2

the equilibrium remains
stable for all τ.

B. One Equation with Multiple Delays

Consider the one dimensional delay differential
equation with the time lags τk, k = 1, 2, ..., and
µ as bifurcation parameters:
dX

dt
= f1(X(t),µ) + f2(X(t− τ1), ...,X(t− τk),µ)

(12)

Fig. 2. The positive equilibrium is stable for τ = 0 and
µ = 2, top graph. The equilibrium still remains stable for
τ = 0.4 (τ < τc = 0.55) and µ = 2, bottom graph

Let (X∗,µ∗) = (X∗,X∗, ...,X∗,µ∗) be the
steady state of equation (12), i.e., f1(X∗,µ∗) +
f2(X

∗,X∗, ...,X∗,µ∗) = 0. To study the stability
of the steady state we compute the characteristic
equation:

λ−
df1
dX
|(X∗,µ∗) −

k∑
j=1

df2
dX
|(X∗,µ∗)e

−λτj = 0.

(13)

For clarity of the presentation we consider the case
of only two delays. Therefore the characteristic
equation is written as

λ−
df1
dX
|(X∗,µ∗) −

df2
dX
|(X∗,µ∗)(e

−λτ1 + e−λτ2 ) = 0

(14)

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 5 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

Note that if τ1 = τ2 = τ or τ1 = 0 or τ2 =
0 then we are back to the previous case of one
equation with one delay. We will assume that τ1
is in its stable domain, i.e., 0 < τ1 < τ1c and
τ2 > 0. We now examine how variation of τ2 and
µ affects the stability of the steady state. Consider
λ(τ2) = α(τ2) + iω(τ2) as solution of equation
(14). We look for a critical value τ2c of τ2 such
that α(τ2c) = 0 and λ(τ2c) = iω(τ2c) = iω2c is
solution of equation (14). Such τ2c exists if and
only if:

iω2c −
df2
dX
|(X∗,µ∗)(cos ω2cτ1 − i sin ω2cτ1)

−
df2
dX
|(X∗,µ∗)(cos ω2cτ2c − i sin ω2cτ2c)

−
df1
dX
|(X∗,µ∗) (15)

Separate real and imaginary parts:

−
df2
dX
|(X∗,µ∗) cos ω2cτ2c =

df2
dX
|(X∗,µ∗) cos ω2cτ1, +

df1
dX
|(X∗,µ∗)

df2
dX
|(X∗,µ∗) sin ω2cτ2c =

−
df2
dX
|(X∗,µ∗) sin ω2cτ1 −ω2c.

(16)

Adding the square of (16) and (II-B) we have

[
df2
dX
|(X∗,µ∗)]

2 = (ω2c +
df2
dX
|(X∗,µ∗) sin ω2cτ1)

2

+(
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗) sin ω2cτ1)

2

(17)

Clearly τ2c exists if and only the function:

H(ω2c) = (ω2c +
df2
dX
|(X∗,µ∗) sin ω2cτ1)

2

+ (
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗) cos ω2cτ1)

2

− [
df2
dX
|(X∗,µ∗)]

2 (18)

has at least a positive root.

Fig. 3. For µ = 2 and τ = 0.9 the equilibrium is unstable,
top graph. For µ = 0.2 and τ = 1 the equilibrium is stable,
bottom graph

Proposition 3: Consider the one dimensional
delay differential equation with the time lag τ1,
τ2, and µ as bifurcation parameters:

dX

dt
= f1(X(t),µ) + f2(X(t− τ1),X(t− τ2),µ).

(19)

Assume that the steady state (X∗,µ∗) =
(X∗,X∗,µ∗) of (19) is stable for 0 < τ1 < τ1c.
If df2

dX
|(X∗,µ∗) > 0 and

df1
dX
|(X∗,µ∗) < 0, then there

exists a critical value τ2c > 0 for τ2 such that
(X∗,µ∗) losses stability as τ2 crosses τ2c.
Proof: Such τ2c exists if and only if equation
H(ω2c) = 0 has at least a positive equation. Or

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 6 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

If df2
dX
|(X∗,µ∗) > 0 and

df1
dX
|(X∗,µ∗) < 0 then

H(0) = (
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗))

2

− (
df2
dX
|(X∗,µ∗))

2 < 0 (20)

And also H(ω2c) → ∞ as ω2c → ∞. Then the
intermediate value theorem assures that equation
H(ω2c) = 0 has at least a positive root. �

We now extend our analysis to a system of n-
delay differential equations with multiple discrete
time delays τ1,τ2, ...,τk, and a local bifurcation
parameter µ.

III. N DIMENSIONAL FIELD

Consider the following system non-linear delay
differential equations:

dx

dt
= f(x(t),x(t− τ1), ...,x(t− τk),µ), (21)

where x ∈ Rn, τj ≥ 0, 1 ≤ j ≤ k are constant
discrete times,
f : Rn+1 × Ck → Rn is assumed to be smooth
enough to guarantee existence and uniqueness of
solutions to (21) under the initial value condition
(R. Bellman and K. L. Cooke [4] and J. K. Hale
and S. M. Verduyn Lunel [13])

x(θ) = φ(θ), θ ∈ [−τ, 0],

where C = C([−τ, 0],Rn), τ = max
1≤j≤k

τj.

Suppose f(x∗,x∗, ...,x∗,µ∗) = 0, that is, (x∗,µ∗)
is a steady state of system (21). We are interested
in studying the stability of such equilibrium point.
In particular studying the effect of the parameter
µ and the discrete time delays τ1,τ2, ...,τk on its
stability. The linearization of (21) at (x∗,µ∗) has
the form (Ruan [9]):

dX

dt
= A0(µ

∗)X(t) +

k∑
j=1

Aj(µ
∗)X(t− τj),

(22)

where X ∈ Rn, each Aj(µ∗) (0 ≤ j ≤ k) is an
n×n constant matrix that depends on values of µ∗

within Dµ. The transcendental equation associated
with system (21) is given as :

det

[
λI −A0(µ∗) −

k∑
j=1

Aj(µ
∗)e−λτj

]
= 0

(23)

Equation (23) has been studied by many re-
searchers
(Ruan [9], R. Bellman and K. L. Cooke [4] and
J. K. Hale and S. M. Verduyn Lunel [13]). The
following result, which was proved by Chin [14]
for k = 1 and by Datko [15] and Hale et al.
[13] for k ≥ 1, gives a necessary and sufficient
condition for the absolute stability of system (22).

Lemma 1: System (22) is stable for all delays
τj(1 ≤ j ≤ k) if and only if
(i) Reλ(

∑k
j=0 Aj(µ

∗)) < 0;
(ii) det[iωI−A0(µ∗)−

∑k
j=1 Aj(µ

∗)e−iωτj ] 6= 0
for all ω > 0

Clearly, the stability of the steady state (x∗,µ∗)
and the effects of the discrete times τj on its
stability depend on values of µ∗ within Dµ. To
further investigate the effects of µ, and the discrete
time delays τj on the stability of (x∗,µ∗), the
exact entries of the matrices Aj(µ∗) are needed to
avoid doing a large number of cases. Note that the
difficulty of the analysis is not due to the number
of delays but to the number of equations. Even
in the case of two equations with one delay, one
needs to consider:

det[λI −A0(µ∗) −A1(µ∗)eλτ = 0],

where

Ai(µ
∗) =

∂f

∂Xi
|(X∗,µ∗), i = 0, 1.

So the stability depends on all the entries of the
Ai, i = 0, 1, we have many different cases.

Therefore to present the ideas we consider a
specific example with n = 2, k = 2, that is a
two dimensional delay differential equations with
two discrete time delays, and a local bifurcation
parameter.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 7 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

A. Two Dimensional Field Example

Consider the non-Kolmogorov type (Holling)
predator-prey model

dx

dt
= r1x(t− τ1) −a1

x(t)y(t)

x(t) + 1
, (24)

dy

dt
= −r2y(t) + a2

x(t− τ2)y(t− τ2)
x(t− τ2) + 1

(25)

where the parameters are described in the fol-
lowing table:

TABLE II
PARAMETER VALUES

Parameters Description Values
x(t) the prey population
y(t) the predator population
r1 the growth rate of the prey in the absence of predators 0.5
r2 the death rate of predators in the absence of the prey 0.5
a1 the predation rate of the prey by the predators 0.5
a2 the conversion rate for the predators 5
τ1 the time lag in the growth to maturity of the prey varies
τ2 the time lag in the growth to maturity of the predators varies

Note that r1 > 0, r2 > 0, a1 > 0, a2 > 0, τ1 > 0,
τ2 > 0.

Proposition 4: If the basic reproductive ratio
(Ameh[16]) R > 1, the system has two non-
negative steady states:

(x∗0,y
∗
0) = (0, 0), and (x

∗
1,y
∗
1) = (

1
R−1,

RR′

R−1 ),

where R = a2
r2

R′ = a1
r1
.

We consider R, τ1 and τ2 as the bifurcation
parameters for the system (24-25) since changes
of them may affect the existence and stability of
the equilibrium points.

B. Stability Analysis

Proposition 5: There exists a critical value for
τ1 such that
(i) The steady state (x∗0,y

∗
0) is unstable for τ1 =

0, and all τ2 ≥ 0.
(ii) The steady state (x∗0,y

∗
0) is stable for τ1 ≥

τ1c, and all τ2 ≥ 0.
Proposition 6: If [(b−d)2−r21−2a1f] < 0 and

∆ = [(b−d)2 −r21 − 2a1f]
2 − 4a21f

2 ≥ 0 then
there exists a critical τ ′1c such that

(i) The steady state (x∗1,y
∗
1) = (

1
R−1,

RR′

R−1 ) is
unstable for 0 ≤ τ1 < τ ′1c and τ2 = 0.

(ii) The steady state (x∗1,y
∗
1) is stable for τ1 > τ

′
1c

and τ2 = 0.
Note that τ1 affects the stability of the positive
equilibrium only for values of R such that condi-
tions C(0) are satisfied.
Remark: For our parameter values, we have

[(b−r2)2 −r21 − 2a1f] = −0.9475 < 0 and

∆ = [(b−r2)2−r21−2a1f]
2−4a21f

2 = 0.0878 > 0

Proposition 7: Consider system (24-25) with τ1
in its unstable interval (0 ≤ τ1 < τ ′1c). If a1 ≥ 2,
then there exists a critical τ2 > 0, such that the
positive equilibrium becomes stable for τ2 > τ2c.
Note that the effect of τ2 on the stability of the
positive equilibrium does not depend on the values
of R.

C. Hopf Bifurcation Analysis

According to the Hopf Bifurcation Theorem
(Culshaw [6]), the discrete time delay τ1 will
cause the system to go through a Hopf bifurcation
near the steady state (x∗1,y

∗
1), if the following

transversality condition is satisfied:

dα(τ1)

dτ1
|τ1=τ′1c 6= 0. (26)

To check this condition we recall that the char-
acteristic equation of the system at (x∗1,y

∗
1) when

τ2 = 0 is given as :

λ2 + (b−r2)λ−r1e−λτ1λ + a1f = 0. (27)

Substituting λ(τ1) = α(τ1) + iω(τ1) in equation
(27), we have :

α2 −ω2 + 2αωi + (b−r2)α + (b−r2)ωi
−r1eατ1 (cos ωτ1 − i sin ωτ1)(α + iω) + a1f = 0.

(28)

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 8 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

We equate the real and the imaginary parts to zero,
and we have :

α2 −ω2 + a1f + r1eατ1 (α cos ωτ1 + ω sin ωτ1)
+(b−r2)α = 0,

(29)

2αω −r1eατ1 (ω cos ωτ1 −α sin ωτ1)
+(b−r2)ω = 0.

(30)

We differentiate equations (29) and (30) with
respect to τ1 and evaluate at τ1 = τ ′1c for which
α(τ ′1c) = 0 and ω(τ

′
1c) = ω

′
c. We obtain

A
dω(τ1)

dτ1
|τ1=τ′1c −B

dα(τ1)

dτ1
|τ1=τ′1c =

C cos ω′cτ
′
1c + D sin ω

′
cτ
′
1c

(31)

B
dω(τ1)

dτ1
|τ1=τ′1c + A

dα(τ1)

dτ1
|τ1=τ′1c =

C sin ω′cτ
′
1c −D cos ω

′
cτ
′
1c

(32)

where

A := 2ω′c −r1τ
′
1c sin ω

′
cτ
′
1c,

B := (b−r2) + r1τ ′1c cos ω
′
cτ
′
1c

C := r1τ
′
1c(ω

′2
c + ω

′
cτ
′
1c),

D := r1τ
′
1cω
′
c sin ω

′
cτ
′
1c.

By solving equations (31) and (32) we have:

dα(τ1)

dτ1
|τ1=τ′1c =

(AC −BD) sin ω′cτ ′1c − (AD + BC) cos ω
′
cτ
′
1c

A2 + B2
.

(33)

The system undergoes through a Hopf bifurcation
near (x∗1,y

∗
1) if:

(AC−BD) sin ω′cτ
′
1c−(AD+BC) cos ω

′
cτ
′
1c 6= 0.

D. Numerical Results

To illustrate the effect of the parameter R
and the discrete time delay on the stability
of the steady state (x∗,y∗), and to support
the theoretical predictions discussed above,
we conducted numerical simulations for the
system (24-25). We used DDE-BIFTOOL
(Engelborghs[17]) for the stability and bifurcation
analysis and also used the Matlab solvers ode23
and dde23 (Shampine[18],Shampine[19]) to see
the behavior of the predator and prey populations
through time. All the parameter values are given
in Table II.

Fig. 4. The positive equilibrium is unstable for τ1 = τ2 = 0
and R = 10 > 1.The system exhibits a spiral out from the
equilibrium (x∗1,y

∗
1) = (0.111,1.111).

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 9 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

For the given parameters values we have R =
10 > 1, and a positive equilibrium exists and is
given as (x∗1,y

∗
1) = (0.111, 1.111). When there

is no delay the prey and predator populations
variation through time is shown on Figure 4.

For our parameter values we have:

[(b−d)2 −r21 − 4f] = −0.9475 < 0 and

∆ = [(b−d)2 −r21 − 4f]
2 − 16f2 = 0.0878 > 0.

Then there exists a τ1c = 6 such that the steady
state remains unstable for 0 ≤ τ1 < τ1c and τ2 = 0
(see Figure 5), it becomes stable as τ1 crosses τ1c
and τ2 = 0 as shown on Figure 6.

Fig. 5. The positive equilibrium remains unstable for τ1 =
1 < τ1c = 6 and τ2 = 0.

Fig. 6. The positive equilibrium is stable for τ1 = 7 >
τ1c = 6 and τ2 = 0.

We examine closely the stability switch
introduces by τ1. We use DDE-BIFTOOL to
compute the eigenvalues of the characteristic
equation (38) for τ2 = 0 and 0 ≤ τ1 ≤ 10.
In Figure 7 we plot the real parts versus the
imaginary parts of these eigenvalues.

We see that the equilibrium (x∗1,y
∗
1) stabilizes

as τ1 crosses the critical value τ ′1c = 6. We also
plot in Figure 7 the eigenvalues of equation (38)
for τ1 = τ ′1c = 6 and observe a pair of two pure
imaginary eigenvalues. The system undergoes
through a Hopf bifurcation as τ1 crosses τ ′1c. We
compute the Hopf bifurcations branches using
Matlab and show them in Figure 8.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 10 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

Fig. 7. The eigenvalues of the characteristic equation (38) for
τ1 = 3 (left) and τ1 = 8(center) with τ2 = 0.At τ1 = τ′1c = 6
we can clearly observe a pair of 2 pure imaginary eigenvalues
(right).

Note τ ′1c = 6 and τ2c = 2.5 For τ1 = 2 and τ2 =
0.5 the equilibrium is unstable as shown in Figure
9.

Fig. 8. Global Hopf bifurcations branches as we vary τ1 and
a1 (same as varying R).

Fig. 9. The positive equilibrium is unstable for τ1 = 2 and
τ2 = 0.5.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 11 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

Fig. 10. The positive equilibrium is stable for τ1 = 7
and τ2 = 1.2.The system exhibits a spiral in toward the
equilibrium (x∗1,y

∗
1) = (0.111,1.111).

For τ1 = 7 and τ2 = 1.2 the equilibrium becomes
stable as shown in Figure 10.

For the case of two non-zero delays, we use Matlab
to compute numerical simulations illustrating the
effects of the two delays. The analysis is summa-
rized in Table III

TABLE III
STABILITY REGIONS IN CASE OF TWO NON-ZERO DELAYS

Unstable Stable Stable Unstable{
0 ≤ τ1 < τ′1c,
0 ≤ τ2 < τ2c

{
τ1 > τ

′
1c,

0 ≤ τ2 < τ2c

{
0 ≤ τ1 < τ′1c,
τ2 > τ2c

{
τ1 > τ

′
1c,

τ2 > τ2c
see Fig 4 and Fig9 see Fig6 and Fig10 see Fig11 see Fig12

Fig. 11. The positive equilibrium is stable for τ1 = 0.7 and
τ2 = 8.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 12 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

Fig. 12. The positive equilibrium is unstable for τ1 = 7 and
τ2 = 3.1.The system exhibits an unstable periodic solutions.

For τ1 = 7 and τ2 = 3.1 the equilibrium becomes
unstable again as shown in Figure 12.

IV. CONCLUSIONS AND DISCUSSION

It is well known that changes in the parameters
play a crucial role in understanding dynamical sys-
tems. There is a need to incorporate discrete time
delays in dynamical systems (Biological systems,
physical systems,...) as has been shown and stud-
ied by many researchers (Perelson[1],Bellen[3],..).

Published papers have shown that the incorpo-
ration of discrete time delays can highly impact
the dynamics of the system, since they can cause
stability switches of a steady state point, and can
also cause the system to go through a Hopf bi-
furcation near that steady state point (Culshaw[6],
Bellen[3],...). The highlight of this paper is on
how a local bifurcation parameter of the system
may modify the stability changes caused by the
delay(s).To understand the effects of discrete time
delays and parameter variations on certain biolog-
ical system models, we carried out a bifurcation
analysis of a system of delay differential equations
in detail for n=1 with specific examples, gave the
procedure for higher n, and did a concrete example
for n=2. We investigated the stability of the steady
states as both bifurcation parameters, the discrete
time delay τ and a local bifurcation parameter µ,
cross critical values. Our analysis shows that while
both parameters can destabilize the steady state,
the discrete time delay can cause stability switches
of the steady state only upon certain values of
µ. The local bifurcation parameter effects on the
stability of the steady state do not depend on the
value of τ. We also showed that both parameters
act differently in term of bifurcation. While the
discrete time delay may only introduce a Hopf
bifurcation, the parameter µ can introduce other
type of bifurcations.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 13 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

V. APPENDIX A

Theorem 1: Consider the transcendental equa-
tion

λn +

n∑
i=1

an−iλ
n−i +

n∑
i=1

bn−ie
−λτλn−i = 0,

(34)

if there exists a τc > 0 such that λ(τc) is a purely
imaginary eigenvalue of (34), then for τ > τc
the transcendental equation (34) has at least one
eigenvalue with a strictly positive real part.

Before we prove the above theorem , let just
first consider a much simpler case. Consider the
analytic function

h(λ,a) = λ + e−λτ + a,

with τ ≥ 0, and a ∈ R.
Then h(λ, 0) = 0 if and only if

λ = −e−λτ. (35)

Equation (35) has purely imaginary roots if and
only if τ = τc = 2jπ + π2 , j = 0, 1, 2, . . .

The proof of the following lemma can be found
in Cooke and Van den Driessche [20]; see also
Bellman and Cooke [4].

Lemma 2: If τ ∈ [0, π
2
), then all roots of

equation (35) have strictly negative real parts. If
τ ∈

(
2jπ + π

2
, (2j + 1)π + π

2

]
, then equation (35)

has exactly 2j + 1 roots with strictly positive real
parts.

We have h(λ,a) is an analystic function in
λ, a. When τ 6= 2jπ + π

2
, the function h(λ, 0)

has no zeros on the boundary of Ω, where Ω =
{λ, |Re(λ) ≥ 0, |λ| ≤ ρ}. Thus, Rouche’s theo-
rem (Dieudonne[11], Theorem 9.17.4) implies that
there exists a δ > 0 such that :
(1) for any a < δ, h(λ,a) has no zero on the

boundary of Ω
(2) for any a < δ, h(λ,a) and h(λ, 0) have the

same sum of the orders of zeros belonging to
Ω.

It follows from lemma 2 that when τ > π
2

, the
sum of the orders of the zeros of h(λ, 0) belonging
to Ω is at least 1. Thus when τ > π

2
, τ 6= 2jπ+ π

2
,

and a < δ then h(λ,a) has at least a root with
strictly positive real part.

Now we can prove the more general form which
is theorem 1

Proof : Consider the analytic function in λ, A

h(λ,A) =λn+

n∑
i=1

an−iλ
n−i+

n∑
i=1

bn−ie
−λτλn−i,

(36)

where λ ∈ C, and
A = (an−1, ...,a1,a0,bn−1, ...,b1) ∈ Rn×(n−1).
Then

h(λ,A0) = λ
n + b0e

−λτ

where A0 = (0, ..., 0) is the null vector. h(λ,A0)
has purely imaginary roots if and only if

τ = τjc =
2jπ

b
1/n
0

j = 1, 2, ... when n is even,

or

τ = τjc =
(4j + 1)π

2b
1/n
0

j = 0, 1, 2, ... when n is odd,

and here we assume that b0 > 0, otherwise we
multiple by a − sign. When τ 6= τjc the function
h(λ,A0) has no zero on the boundary of Ω, where
Ω = {λ, |Re(λ) ≥ 0, |λ| ≤ ρ}. Thus, Rouche’s
theorem implies that there exists a δ > 0 such
that :
(1) when ‖A‖∞ < δ, h(λ,A) has no zero on the

boundary of Ω
(2) when ‖A‖∞ < δ, h(λ,A) and h(λ,A0) have

the same sum of the orders of zeros belonging
to Ω.

It follows from lemma 2 that when τ > τc =
2π
b
1/n
0

and τ 6= 2jπ
b
1/n
0

(or τ > τc = 1π
2b

1/n
0

and τ 6=
(4j+1)π

2b
1/n
0

), the sum of the orders of the zeros of
h(λ,A0) belonging to Ω is at least 1. Thus when
τ > τc , τ 6= τ

j
c and ‖A‖∞ < δ then h(λ,A) has

at least a root with strictly positive real part. �

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 14 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

VI. APPENDIX B

A. Proof of Proposition 1

The characteristic equation of the one
dimensional DDE is given by (3) and because the
steady state is assumed to be stable at τ = 0 then

α(0) =
df1
dX
|(X∗,µ∗) +

df2
dX
|(X∗,µ∗) < 0. (37)

If (i) holds then equation (37) implies
|df2
dX
|(X∗,µ∗)| < |

df1
dX
|(X∗,µ∗)|

therefore equation (9) has no positive root
meaning the steady state remains stable for all
τ ≥ 0.

If (ii) holds then equation (37) implies
|df2
dX
|(X∗,µ∗)| > |

df1
dX
|(X∗,µ∗)|

therefore equation (9) has a positive root then
there exists a τc > 0 such that α(τ) > 0 whenever
τ > τc.

If (iii)(a) holds then again equation (9) has no
solution therefore α(τ) < 0 for all τ ≥ 0 meaning
the steady state remains stable.

If (iii)(b) holds then equation (9) has a positive
root then there exists a τc > 0 such that α(τ) > 0
whenever τ > τc. �

B. Proof of Proposition 2

If conditions of proposition 1(iii)(a) hold then
there is nothing to prove. Assume that conditions
of proposition 1(iii)(b) hold then equation (9) has
a positive solution, therefore the delay can affect
the stability of the equilibrium point. But if for
some µ∗ in Dµ we have the extra condition

df1
dX
|(X∗,µ∗) = g(X

∗)O(µ∗)

and
df2
dX
|(X∗,µ∗) = h(X

∗)O(
1

µ∗
),

then one can rewrite equation (9) as

ω2c = [
h(X∗)

µ∗
]2 − [g(X∗)µ∗]2.

Then there exists a critical value µc ∈ Dµ of µ
such that

h(X∗)

µ∗
≈ 0 as µ∗ → µc.

Therefore equation (9) becomes

ω2c = −[g(X
∗)µc]

2 < 0,

which has no real positive root ωc, therefore
α(τ) < 0 for all τ ≥ 0. This implies the delay
does not have any effect on the stability of the
equilibrium point when µ∗ > µc. �

C. Proof of Proposition 5

The Jacobian matrix of the system (24-25) is
given by :

J =

[
r1e
−λτ1 − a1y

∗

(x∗+1)2
− a1x

∗

x∗+1
a2y
∗

(x∗+1)2
e−λτ2 −r2 + a2x

∗

x∗+1
e−λτ2

]
.

Evaluating at (x∗,y∗) = (0, 0), the characteristic
equation is given as

(λ−r1e−λτ1 )(λ + r2) = 0.

We note that the stability of (x∗,y∗) = (0, 0)
depends only on τ1.
• If τ1 = 0 then the eigenvalues are :
λ = r1 > 0 and λ = −r2 < 0. Therefore the
(0, 0) is unstable.

• If τ1 > 0, we have λ = r1e−λτ1 , let λ(τ) =
α(τ) + iω(τ) then we have
λ = r1e

−ατ1 (cos ωτ1 − i sin ωτ1).
One can choose ωcτ1c =

π(2n+1)
2

(n=0,1,2,...)
or τ1c =

π(2n+1)
2ωc

such that the real part of
λ(τ) = α(τ) + iω(τ) at τ1c is zero (α(τ1c) =
0) and the imaginary part ω(τ1c) = ωc is a
solution of the characteristic equation.
Then by the continuity of α we have :

– α(τ) > 0 for τ1 < τ1c,
– α(τ) < 0 for τ1 > τ1c. �

D. Proof of Proposition 6

The characteristic equation of the system eval-
uating at (x∗1,y

∗
1) is given by

λ2 + (b−r1e−λτ1 −r2e−λτ2 )λ + (f −f1e−λτ1 )+
(a1 − 1)fe−λτ2 + f1e−λ(τ1+τ2) = 0,

(38)

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 15 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

where

b = r2 +
r1

x∗1 + 1
, f =

r1r2
x∗1 + 1

, f1 = r1r2.

• If τ1 = τ2 = 0 we have:

λ2 − (b−r1 −r2)λ + a1f = 0

with

b−r1 −r2 = −
r1
R
< 0, and a1f > 0.

Then the characteristic equation has at least
a positive eigenvalue
(if the eigenvalues are real) λ = r1

2R
+
√
δ
2

where

δ = (b−r1 −d)2 − 4a1f,

or, all its eigenvalues (if complex) have a
positive real part ( r1

2R
). Therefore the steady

state (x∗1,y
∗
1) is unstable.

• If τ1 > 0 and τ2 = 0 Then the characteristic
equation becomes:

λ2 + (b−r2)λ−r1e−λτ1λ + a1f = 0.

Since we know that the steady state is unsta-
ble when τ1 = τ2 = 0, the question becomes:
does there exist a τ ′1c such that the steady state
stabilizes as τ1 crosses τ ′1c? In other words if
λ(τ1) = α(τ1) + iω(τ1), does there exist τ ′1c
such that α(τ ′1c) = 0 and ω(τ

′
1c) = ω

′
c which

satisfies

−ω′2c + i(b−r2)ω
′
c−

ir1ω
′
c(cos ω

′
cτ
′
1c − i sin ω

′
cτ
′
1c) + a1f = 0.

(39)

Setting the real and imaginary parts equal
zero, we obtain:

−ω′2c + a1f = r1ω
′
c sin ω

′
cτ
′
1c (40)

(b−d)ω′c = r1ω
′
c cos ω

′
cτ
′
1c. (41)

Adding the square of both equations, we
obtain:

ω′4c +[(b−r2)
2 −r21 − 2a1f]ω

′2
c + a

2
1f

2 = 0
(42)

Such τ ′1c exists if and only if the above
equation has at least a positive root ω′c. Let
M = ω′2c , then we have the quadratic equa-
tion:

M2 + [(b−r2)2 −r21 − 2a1f]M + a
2
1f

2 = 0
(43)

which has at least a positive root if:
C(0) : [(b − r2)2 − r21 − 2a1f] < 0, and
∆ = [(b − r2)2 − r21 − 2a1f]

2 − 4a21f
2 ≥

0, consequently, equation (42) has at least a
positive root ω′c. Which implies there exist a
τ ′1c > 0 such that the steady state changes
stability as τ1 crosses τ ′1c for τ2 = 0. In fact
τ ′1c is the smallest of :

τ
′j
1c =

1

ω′c
arccos

b−r2
r1

+
2πj

ω′c
, j = 1, 2, ...�

REFERENCES

[1] A. S. Perelson, D. E. Kirschner, R. DE Boer, Dynamics
of HIV infection of CD4+ T cells,Math Biosci. 114
(1993) 81-125

[2] L. Allen, An Introduction to Mathematical Biology,
Pearson-Prentice Hall, Upper Saddle River, NJ, 2007.

[3] A. Bellen, M. Zennaro, Numerical Methods for Delay
Differential equations, Oxford University Press, NY
2005.

[4] R. Bellman, K. L. Cooke, Differential-Difference Equa-
tions, Academic Press, New York, 1963.

[5] Driver, Rodney D. Ordinary and Delay Differential
Equations, New York: Springer Verlag (1977). ISBN 0-
387-90231-7.

[6] R. V. Culshaw, S. Ruan, A delay-differential equation
model of HIV infection of CD4+ T-cells, Mathematical
Biosciences 165 (2000) 27–39.

[7] S. Gakkhar, A. Singh, Complex dynamics in a prey
predator system with multiple delays, Commun Non-
linear SCi Numer Simulat. 17 (2012) 914-929.

[8] K. Gopalsamy Stability and Oscillations in Delay Dif-
ferential Equations of Population Dynamics, Kluwer
Academics Publishers, 1992.

[9] S. Ruan, On Nonlinear Dynamics of Predator-Prey
Models with discrete Delay, Math. Model. Nat. Phenom.
4 (2009) 140-188.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 16 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201


I. Diakite et al., Effects of Discrete Time Delays and Parameters ...

[10] M. Y. Li, H. Shu, Joint effects of mitosis and intracellu-
lar delay on viral dynamics:two-parameter bifurcation
analysis, Math. Biol. 64 (2011) 1005-20.

[11] J. Dieudonne, Foundations of Modern Analysis, Aca-
demic Press, New York. 1960

[12] J. Guckenheimer, P. Holmes, Nonlinear Oscillations,
Dynamical Systems, and Bifurcations of Vector Fields,
Springer-Verlag, New York, 1997.

[13] J. K. Hale, S. M. V. Lunel, Introduction to Functional
Differential Equations, Springer Verlag, New York,
1993.

[14] Y. S. Chin, Unconditional stability of systems with time-
lags., Acta Math. Sinica, 1 (1960) 138-155.

[15] R. Datko, A procedure for determination of the exponen-
tial stability of certain differential difference equations.
Quart. Appl. Math. 36 (1978) 279-292.

[16] J. E. Ameh and R. Ouifki, The Basic Reproductive Num-

ber: Bifurcation and Stability, Thesis, African Institute
for Mathematical Sciences. 2009.

[17] K. Engelborghs, T. Luzyanina, G. Samaey, DDE-
BIFTOOL v. 2.00: a Matlab package for bifurcation
analysis of delay differential equations, Report TW 330,
2001.

[18] L. F. Shampine, I. Gladwell, S. Thompson,Solving
ODEs with MATLAB, Cambridge University Press,
Cambridge 2003.

[19] L. F. Shampine, S. Thompson, Solving DDEs in MAT-
LAB, Applied Numerical Mathematics 37 (2001) 441-
458.

[20] K.L. Cooke, P.vanden Driessche, and X. Zou, Interac-
tion of maturation delay and nonlinear birth in popu-
lation and epidemic models, J. Math. Biol, 39 (1999)
332-352.

Biomath 1 (2015), 1505201, http://dx.doi.org/10.11145/j.biomath.2015.05.201 Page 17 of 17

http://dx.doi.org/10.11145/j.biomath.2015.05.201

	Introduction
	One Dimensional Field
	One Equation with One Delay
	Stability of the Steady State
	Example

	One Equation with Multiple Delays

	n Dimensional Field
	Two Dimensional Field Example
	Stability Analysis
	Hopf Bifurcation Analysis
	Numerical Results

	Conclusions and Discussion
	Appendix A
	Appendix B
	Proof of Proposition 1
	Proof of Proposition 2
	Proof of Proposition 5
	Proof of Proposition 6

	References