Adv Syst Sci Appl 2022; 03:1–17
Published online at https://ijassa.ipu.ru.

Stability and Bifurcation of a Delay Cancer Model
in the Polluted Environment

Ahmed A. Mohsen1*, Raid K. Naji2

1Department of Mathematics, College of Education for pure science (Ibn-Al-Haitham), University of Baghdad.
2Department of mathematics, College of science, University of Baghdad.

Abstract: It is well known that the spread of cancer or tumor growth increases in polluted
environments. In this paper, the dynamic behavior of the cancer model in the polluted
environment is studied taking into consideration the delay in clearance of the environment
from their contamination. The set of differential equations that simulates this epidemic model
is formulated. The existence, uniqueness, and the bound of the solution are discussed. The local
and global stability conditions of disease-free and endemic equilibrium points are investigated.
The occurrence of the Hopf bifurcation around the endemic equilibrium point is proved. The
stability and direction of the periodic dynamics are studied. Finally, the paper is ended with a
numerical simulation in order to validate the analytical results.

Keywords: Cancer, Environment Pollution, Time Delay, Stability, The Hopf Bifurcation.

1. INTRODUCTION

Cancer is an abnormal growth of cells that tend to proliferate in an uncontrolled way and in
some cases to spread. In fact, cancer or tumor growth can be considered as not one disease but
it is a group of more than one hundred different and distinctive diseases. Cancer can involve
any tissue of the body and have many different forms in all body areas [2].

According to the World Health Organization, the cancer is one of the global diseases
that affect the human as much as 18 million new cases annually, and more than half of the
mortalities in 2018, and the numbers are expected to nearly double in 2040. Although there
is many advanced experimental in developing interventional therapies for cancer such as
immunotherapy, virotherapy, targeted drug therapies, and chemotherapy, in addition to the
surgical resection, treatment options are still limited and the disease is considered fatal [5].
Therefore, we can consider the mathematical models as one of the ways for studying this
disease. For a brief review of mathematical models describing tumor-immune dynamics see
Tsygvintsev et al. [12]. Lestari et al [6] proposed a mathematical model of the spread of
cancer with chemotherapy and then studied their dynamical behavior. Weerasinghe et al [14]
discussed some models of plasticity, tumor progression, and metastasis using three broadly
conceived mathematical modeling techniques: discrete, continuum, and hybrid, each with
advantages and disadvantages. Simmons et al. [9] presented a brief overview of breast cancer,
focusing on its heterogeneity, and then explained the role of mathematical modeling and
simulation in teasing apart the underlying biophysical processes. In [8], Pillis and Radunskaya
proposed and investigated the immune response to tumor invasion. Trisilowati et al. [11]
proposed and analyzed the optimal control model of dendritic cell treatment of growing

∗Corresponding author: aamuhseen@gmail.com



2 A.A. MOHSEN, R.K. NAJI

tumors. Moreover, Wilson and Levy [15] formulated a model to investigate the effect of
immunotherapy on tumor growth. Recently, Abernathy et al [16], studied the long-term
dynamics of a system of nonlinear differential equations that describes the role of virotherapy
on tumors and the impact of immune response specific to fighting cancer. Al-Tuwairqi et al [1]
presented a realistic mathematical model that describes the interaction between the innate
immune system and uninfected tumor cells. This is based on the fact that both tumor and
virus-infected cells are recognized by natural killer cells that are part of the innate immune
system

On the other hand, there have been several papers proposed to show the effect of the time
delay of some epidemic disease models as seen in Zuo et al. [17], which studied the impact of
the media and delay on the spread of an epidemic. In [3], Cooke et al. proposed an epidemic
model with two delays. Wang and Wu [13], investigated an SEIR model with delay. Recently,
Naji and Majeed [7] studied the effect of delay on a stage structure prey-predator model.

In this article, it is assumed that cancer start outbreak due to environmental pollution,
and then the effect of delay in cleaning up this environment on the dynamic behavior of
the proposed system is discussed. So that, the next section deals with the formulation of the
model and presents its main properties. In section 3, the stability analysis of the disease-
free equilibrium point and endemic equilibrium point is discussed and the possibility of
occurrence of the Hopf bifurcation is also studied. Section 4, discusses the stability and
direction of the periodic dynamics that resulting from the Hopf bifurcation as the delay
starting increases. Numerical simulations are used in section 5 to further understand the
dynamics of the model. Finally, the paper is ended with a discussion section as given in
section 6.

2. THE MATHEMATICAL FORMULATION OF CANCER MODEL

It is well known that one of the main reasons for the spread of cancer is the polluted
environment. Hence, in case of existence of such disease in a population of size N(t) at time
t, then the population is divided into two compartments, namely susceptible compartment,
which is denoted to their population’s size at time t by S(t) and cancer compartment, which
contains all the individuals that infected by cancer and denoted to their population’s size at
time t by C(t), such that N (t) = S (t) + C(t). Moreover, it is assumed that the pollution
level at time t in the environment is given by E(t). On the other hand, since the delay in
clearance of the environment from their contamination has a vital effect on the spread of
cancer. Therefore, the effect of delay, with the amount of delay given by τ > 0, is considered
in the formulation of the model. Accordingly, the dynamics of the spread of cancer within
such an environment can be represented in the following block diagram.

Fig. 2.1. The block diagram of Cancer model.

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 3

And the corresponding dynamic model has been formulated by the following system of
nonlinear delay differential equations

dS
dt

= ψ − β (1−m)S (t)E (t)− µS (t)

dC
dt

= β (1−m)S (t)E (t)− (µ+ α)C (t)

dE
dt

= P − θE (t− τ)

 , (2.1)

where S (r) = S (0) > 0, C (r) = C(0) ≥ 0 and E (r) = E (0) ≥ 0 for r ∈ [−τ, 0). It is
assumed that all parameters in system (2.1) have the following description: ψ > 0 is the
natural birth of susceptible; µ > 0 represents the death rate of all populations; β > 0
represents the contact rate of the susceptible with the contaminated environment. However,
the constant m represents the body resistance rate due to chromosomal fixed against
environment effect such that 0 ≤ m ≤ 1. While P ≥ 0 is the pollution rate in the environment
and θ > 0 is the clearance rate in the environment.

Note that, all the interaction functions on the right-hand side of the system (2.1) are
continuous and have continuous partial derivatives and hence the system has a unique solution
in R3

+. Furthermore, the boundedness of the solution is shown in the following theorem.

Theorem 2.1:
All the solutions of system (2.1) with initial conditions belong to R3

+, are bounded.

Proof
Let (S (t) , C (t) , E (t)) be any solution of system (2.1) with initial conditions belongs to R3

+.
Hence, by adding the first two equations in system (2.1) to each other we get:

dS

dt
+
dC

dt
= ψ − µN − αC.

Therefore, it is obtained that:
dN

dt
+ µN ≤ ψ.

Then, solving the above inequality by using the Gronwall lemma we obtain for t→ ∞
that:

N ≤ ψ

µ
.

Similarly, from the 3rd equation in the system (2.1) we have:

dE

dt
= P − θE.

This implies that:
dE

dt
+ θE = P.

Hence for t→ ∞, we have E(t) ≤ P
θ

. This completes the proof.

Now, since the variable C in the second equation of system (2.1) dose not appear in the
other two equations, then the following system, say system (2.3), can be solved independently.
Then by substituting the obtained solution, say (Sυ, Eυ) in second equation of system (2.1),
we obtain the following solution:

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



4 A.A. MOHSEN, R.K. NAJI

C (t) =
β (1−m)SυEυ

µ+ α
, (2.2)

where (Sυ, Eυ) is a solution of the system:

dS
dt

= ψ − β (1−m)S (t)E (t)− µS (t)

dE
dt

= P − θE (t− τ)

 . (2.3)

3. ANALYSIS OF SYSTEM 2.3

In this section, the existence and their stability analysis are carried out. It is easy to verify that
system (2.3) has two equilibrium points and can be described as follows:

1. The first point is called the disease-free equilibrium point DFEP, which is denoted by
e0 = (S0, 0), where:

S0 =
ψ

µ
. (3.4)

Clearly, e0 exists provided that:
P = 0. (3.5)

2. The other point is called the endemic equilibrium point EEP and that denoted by
e1 = (S1, E1), where:

S1 =
ψθ

β(1−m)P+µθ

E1 =
P
θ

 . (3.6)

Obviously, e1 exists uniquely in the interior of the positive quadrant of the SE−plane
provided that:

P > 0. (3.7)

Now, to analyze the stability of the above equilibrium points the general Jacobian matrix
for the system (2.3) at any point, say (S,E), can be written as:

J (S,E) =

[ − [β (1−m)E + µ] −β (1−m)S

0 −θe−λτ

]
. (3.8)

Accordingly, the Jacobian matrix of system (2.3) at the DFEP that given in Eq. (3.4),
becomes:

J (e0) =

[ −µ −β (1−m)S0

0 −θe−λτ

]
. (3.9)

Then, the characteristic equation of the matrix (3.9) can be written as:

λ2 + µλ+ (θλ+ µθ) e−λτ = 0. (3.10)

Obviously, when τ = 0, Eq. (3.10) becomes:

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 5

λ2 + (µ+ θ)λ+ µθ = 0. (3.11)

Clearly, the eigenvalues in such a case can be written as:

λ
[τ=0]
1 = −µ, λ[τ=0]

2 = −θ.

Hence, e0 is always locally asymptotically stable for τ = 0.
Furthermore, by using the function h(S,E) =

(
1
SE

)
as the Dulac function, which

is a continuously differentiable unction in a simply connected Ω ⊂ R2
+ such that

∂(hf(S,E)/∂S + ∂(hg(S,E)/∂E dose not change the sign in Ω and vanishes at most on
a set of measure zero then the system (2.3) not have periodic orbits in Ω. Therefore, it is
easy to verify that system (2.3) with τ = 0 has no periodic dynamics and hence e0 is always
globally asymptotically stable. On the other hand, since the DFEP exists under the condition
(3.5), hence system (2.3) can be reduce to a single equation due to the fact É < 0. Therefore,
the system (2.3) has no periodic dynamics when τ > 0 and P = 0, and hence e0 is always
globally absolutely stable [10] for all τ ≥ 0.

Now the stability of the EEP of system 2.3 that given by (3.6), is investigated using the
linearization technique. Accordingly, the Jacobian matrix that given in (3.8) of system (2.3)
at the EEP can be written as:

J (e1) =

[ − [β (1−m)E1 + µ] −β (1−m)S1

0 −θe−λτ

]
. (3.12)

Hence, the characteristic equation of J(e1) can be written:

λ2 +B1λ+ (B2λ+B3) e
−λτ = 0, (3.13)

where B1 = β (1−m)E1 + µ, B2 = θ, and B3 = θ [β (1−m)E1 + µ] .
So for τ = 0, then Eq. (3.13) becomes:

λ2 + (B1 +B2)λ+B3 = 0. (3.14)

It easy to see that the roots of the above equation have negative real parts. Consequently,
whenever the EEP exists, it is always locally asymptotically stable (LAS). Indeed, it is
globally asymptotically stable as mentioned above with the help of the Dulac criterion.

Now, for τ > 0, assume that Eq. (3.13) has a pair of purely imaginary roots, say λ = ±iω
with ω > 0, cross the imaginary axis. Hence, by substituting λ = iω in Eq. (3.13) we get:

−ω2 + iB1ω + (iB2ω +B3) (cosωτ − isin ωτ ) = 0.

Therefore, by separating the real and imaginary parts we obtain:

B1ω = B3sinωτ −B2ω cosωτ
−ω2 = −B2ωsinωτ −B3cosωτ

}
. (3.15)

Clearly, squaring the equations in Eq. (3.15) and adding to each other gives that:

ω4 +
(
B2

1 −B2
2

)
ω2 −B2

3 = 0, (3.16)

set K = ω2, then Eq. (3.14) becomes:

K2 +
(
B2

1 −B2
2

)
K −B2

3 = 0. (3.17)

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



6 A.A. MOHSEN, R.K. NAJI

Clearly, by using Descartes’ rule of signs, Eq. (3.17) has a unique positive root, namely
K = ω2

0 . Hence, ω0 is a positive root of Eq. (3.16) too. Therefore, there is at least a pair of
purely imaginary roots ±iω0 satisfying Eq.(3.13).

Substituting ω0 in Eq. (3.15) and then solving the resulting equation with respect to τ we
get that:

τj =
1

ω0

sin−1ω0 (B1B3 +B2ω0
2)(

B2
2ω0

2 +B3
2
) +

2jπ

ω0

; j = 0, 1, 2, · · · (3.18)

Now, define that τ0 = minj≥0 τj , then λ (τ) = γ (τ) + iω(τ) be a root of Eq. 3.13, such
that γ (τ0) = 0 and ω (τ0) = ω0. Then we have the following theorem.

Theorem 3.1:
Assume that condition (3.7) holds. Then the EEP is conditionally stable.

Proof
It is well known that the equilibrium point e1 is conditionally stable if it is asymptotically
stable for τ ∈ [0, τ0). Moreover, as shown above, the equilibrium point e1 is globally
asymptotically stable for τ = 0 and the transcendental characteristic equation given by Eq.
(3.13) has at least a pair of purely imaginary roots ±iω0 at τ = τ0.

Consequently, if we assume that λ (τ) = γ (τ) + iω(τ) is an eigenvalue of Eq. (3.13) with
γ (τ0) = 0 and ω (τ0) = ω0 > 0, where τ0 is given by Eq. (3.18) then the proof follows if we
can show that the following transversality condition hold.[

d(Reλ (τ))

dτ

]
τ=τ0

̸= 0. (3.19)

By deriving Eq. (3.13) with respect to τ we obtain that:[
2λ+B1 +B2e

−λτ − τ(B2λ+B3)e
−λτ] dλ

dτ
= λ(B2λ+B3)e

−λτ . (3.20)

Therefore, we obtain that:[
dλ

dτ

]−1

=
(2λ+B1)

−λ(λ2 +B1λ)
+

B2

λ(B2λ+B3)
− τ

λ
. (3.21)

Since λ = iω0 at τ = τ0, then Eq. (3.21) can be rewritten as:[
dλ

dτ

]−1

τ=τ0

=
2iω0 +B1

B1ω0
2 + iω0

3
+

B2

−B2ω0
2 + iω0B3

− τ0
iω0

.

Now since

sgn

[
d(Reλ)

dτ

]
τ=τ0

= sgn

[
Re

(
dλ

dτ

)−1
]
λ=iω0

. (3.22)

Accordingly, from the fact:

Re

[
2iω0 +B1

B1ω0
2 + iω0

3

]
=

2ω0
2 +B1

2

ω0
2(B1

2 + ω0
2)
,

Re

[
B2

iω0B3 −B2ω0
2

]
=

−B2
2

B2
2 ω0

2 +B3
2 ,

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 7

Re

[
τ0
iω0

]
= 0.

Hence, we obtain that:[
Re

(
dλ

dτ

)]−1

τ=τ0

=
2ω0

2 +B1
2

ω0
2(B1

2 + ω0
2)

− B2
2

B2
2 ω0

2 +B3
2 > 0.

Therefore, d(Reλ(τ))
dτ

at τ = τ0 is not equal zero and does not change sign. Clearly, the
obtained result shows that the eigenvalue of characteristic equation Eq. (3.13) crosses the
imaginary axis from left to right as τ passes through τ0. Hence system (2.3) losses its stability
and undergoes the Hopf bifurcation at τ = τ0. Thus the proof is complete.

4. STABILITY AND DIRECTION OF THE HOPF BIFURCATION

In the previous section, we have shown that system (2.3) undergoes the Hopf bifurcation near
the EEP at τ = τj . In this section, however, the direction of the Hopf bifurcation and their
stability at τ = τ0 is investigated. The normal form theory and center manifold theorem are
used in this investigation, see [4].

Let x1 = S − S1, x2 = E − E1, xi(t) = xi(τt) and τ = τ0+u, then we drop the bars for
simplification of notations. Therefore system (2.3) is transformed into functional differential
equations in C = C ([−1, 0] , R2) as:

ẋ (t) = Lu (xt) +H (u, xt) , (4.23)

where x (t) = (x1, x2)
T ∈ R2, Lu : C → R2 and H:R×C → R2 are given as:

Lu (Γ) = (τ0 + u)

[
− [β (1−m)E1 + µ] −β(1−m)S1

0 0

] [
Γ1(0)
Γ2(0)

]
+(τ0 + u)

[
0 0
0 −θ

] [
Γ1(−1)
Γ2(−1)

]
 , (4.24)

and

H (u,Γ) = (τ0 + u)

[
−β(1−m)Γ1(0)Γ2(0)

0

]
, (4.25)

here Γ(v) = (Γ1(v),Γ2(v))
T ∈ C. Now, by the Riesz representation theorem, there exists a

function of bounded variation, say η(v, u) for v ∈ [−1, 0], such that:

LuΓ =

∫ 0

−1

dη (v, u) Γ (v),Γ ∈ C. (4.26)

Here we can choose:

η (v, u) = (τ0 + u)

[
− [β (1−m)E1 + µ] −β(1−m)S1

0 0

]
δ (v)

− (τ0 + u)

[
0 0
0 −θ

]
δ (v + 1)

 , (4.27)

where δ(v) is the Dirac function, which define as follows, then Eq. (4.26) is satisfied.

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



8 A.A. MOHSEN, R.K. NAJI

δ (v) =

{
1 v = 0.
0 v ̸= 0. (4.28)

For Γ ∈ C1 ([−1, 0] , R2) , define:

A (u) Γ(v) =

{
dΓ(v)
dv

− 1 ≤ v < 0.∫ 0

−1
dη (ς, u) Γ (ς) v = 0.

(4.29)

And

R (u) Γ(v) =

{
0 − 1 ≤ v < 0.
H (u,Γ) v = 0. (4.30)

Hence, system (4.23) can be transformed into an operator differential equation of the
form:

x́t = A (u)xt +R (u)xt, (4.31)

where xt = x (t+ v), v ∈ [−1, 0].
Now, we defined the adjoint operator A∗ of A, where ψ ∈ C1

(
[−1, 0] , (R2)

∗), by:

A∗ (0)ψ(s) =

{
−dψ(s)

ds
0 < s ≤ 1,∫ 0

−1
dηT (ς, 0)ψ (−ς) s = 0,

(4.32)

with a bilinear inner product form:

⟨ψ (s) ,Γ(v)⟩ = ψ(0)Γ (0)−
∫ 0

v=−1

∫ v

ξ=0

ψ
T
(ξ − v) dη (v) Γ (ξ) dξ, (4.33)

where, η (v) = η (v, 0). Clearly, A (0) and A∗(0) are adjoint operators. Now since ±iτ0ω0 are
the eigenvalues of A (0) as we show in the previous section then they are also eigenvalues
of A∗(0). Furthermore, we need to find the eigenvector q of A corresponding to iτ0ω0 and
the eigenvector q∗ of A∗ corresponding to −iτ0ω0 so that they satisfying the normalization
conditions ⟨q∗, q⟩ = 1 and ⟨q∗, q⟩ = 0, as determined in the following theorem.
Theorem 4.1:
The eigenvectors of A (0) and A∗ corresponding to iτ0ω0 and −iτ0ω0 are given by q (v) =
(1, r1)

T eivτ0ω0 and q∗ (s) = D(1, r2)
T eisτ0ω0 respectively, where:

r1 = − [iω0 + β (1−m)E1 + µ]

β (1−m)S1

,

r2 =
β (1−m)S1

iω0 − θeiω0τ0
,

D =
[
1 + r1r2 − τ0θr1r2e

−iω0τ0
]−1

.

Proof
Since q (v) = (1, r1)

T eivτ0ω0 , is the eigenvector of A(0) corresponding to iτ0ω0, then:

A (0) q (v) = iτ0ω0 q (v) .

So, we obtain:

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 9

A (0) q (0) eivτ0ω0 = iτ0ω0 q (0) e
ivτ0ω0 .

Now according to the definition of A (0), we get

τ0

(
iω0 + β (1−m)E1 + µ β (1−m)S1

0 iω0 + θe−iω0τ0

)(
1
r1

)
=

(
0
0

)
.

Consequently, we obtain that:

r1 = − [iω0 + β (1−m)E1 + µ]

β (1−m)S1

.

Similarly, since q∗ (s) = D(1, r2)
T eisτ0ω0 is the eigenvector of A∗ corresponding to

−iτ0ω0, hence:
A∗q∗ (s) = −iτ0ω0q

∗ (s) .

From the definition of A∗ we obtain:

τ0

(
iω0 − [β (1−m)E1 + µ] 0

−β (1−m)S1 iω0 − θe−iω0τ0

)(
1
r2

)
=

(
0
0

)
.

Such that:

r2 =
β (1−m)S1

iω0 − θeiω0τ0
.

To guarantee that ⟨ q∗ (s) , q (v) ⟩ = 1, then the parameter D should be determined, so by
using Eq. (4.33), we get that:

⟨q∗ (s) , q(v)⟩ = D (1, r2)

(
1
r1

)
−
∫ 0

v=−1

∫ v
ξ=0

D (1, r2) e
−i(ξ−v)ω0τ0dη (v)

(
1
r1

)
eiξω0τ0dξ

= D

{
1 + r1r2 −

∫ 0

−1
(1, r2) ve

ivω0τ0dη (v)

(
1
r1

)}
= D

{
1 + r1r2 +

[
τ0 (1, r2)

(
0

−θr1

)]
e−iω0τ0

}
= D {1 + r1r2 − τ0θr1r2e

−iω0τ0} .

Thus we can choose D so that:

D =
[
1 + r1r2 − τ0θr1r2e

−iω0τ0
]−1

.

In addition, it can be easily verify that ⟨q∗, q⟩ = 0 by applying the adjoint property
⟨φ,AΓ⟩ = ⟨A∗φ,Γ⟩, as follows:

iτ0ω0 ⟨q∗, q⟩ = ⟨−iτ0ω0q
∗, q⟩ = ⟨A∗q∗, q⟩ = ⟨q∗, Aq⟩

= ⟨q∗, iτ0ω0q⟩ = −iτ0ω0 ⟨q∗, q⟩ .

Then, we have ⟨q∗, q⟩ = 0.
In the following, the technique given by Hassard et al. [4] to compute the coordinates

describing center manifold C0 at u = 0 is used.
Define

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



10 A.A. MOHSEN, R.K. NAJI

z (t) = ⟨q∗, xt⟩
w (t, v) = xt (v)− z (t) q (v)− z (t) q (v) = xt (v)− 2Re {z (t) q (v)}

}
, (4.34)

where xt (v) = x(t+ v) be the solution of Eq. (4.24). Also, on the center manifold C0, we
have:

w (t, v) = w (z (t) , z (t) , v) , (4.35)

here w (z (t) , z (t) , v) = w20 (v)
z2

2
+ w11 (v) zz + w02 (v)

z2

2
+ . . . , such that, z(t) and z(t)

are local coordinates of center manifold C0 in the direction of q∗ and q∗, respectively. Clearly,
we know that xt ∈ C0 of Eq. (4.31), since u = 0, we know that ⟨φ,AΓ⟩ = ⟨A∗φ,Γ⟩, for
⟨Γ, φ⟩ ∈ D(A)×D(A), then:

ź(t) = ⟨q∗, x′(t)⟩ = ⟨q∗, A(0)xt +R(0)xt⟩ = ⟨q∗, A (0)xt⟩+ ⟨q∗, R (0)xt⟩ = ⟨A∗(0)q∗, xt⟩+ ⟨q∗, R (0)xt⟩ .

Thus

ź (t) = iω0τ0z (t) + q∗
T
(0)H(0, w (t, 0) + 2Re {z (t) q(0)})

= iω0τ0z (t) + q∗
T
(0)H0 (z, z)

}
, (4.36)

here H0 (z, z) = H(0, z, z), Now rewrite Eq. (4.36) as:

ź (t) = iω0τ0z (t) + g (z, z) , (4.37)
where

g (z, z) =
1

2
g20z

2 + g11zz +
1

2
g02z

2 +
1

2
g21z

2z + . . .

So according to the Eq. (4.34) with the Eq. (4.31), we obtain:

ẃ (t, v) = x′t (v)− ź (t) q (v)− ź (t) q (v)

= A(0)xt +R(0)xt − iω0τ0z(t)q(v)− q∗
T
(0)H0(z, z)q(v)

+iω0τ0z (t) q (v)− q∗T (0)H0 (z, z) q (v)

}
= A(0)xt +R(0)xt − A(0)z(t)q (v)− A(0)z(t)q (v)

−2Re
{
q∗
T
(0)H0(z, z)q(v)

} }

= A(0)w(t, v) +R(0)xt − 2Re
{
q∗
T
(0)H0(z, z)q(v)

}
.

Therefore, it is obtain that:

ẃ (t, v) = A (0)w (t, v)− 2Re
{
q∗
T
(0)H0 (z, z) q (v)

}
for v ∈ [−1, 0) .

A (0)w (t, 0) +H0 (z, z)− 2Re
{
q∗
T
(0)H0 (z, z) q (0)

}
for v = 0.

(4.38)

Accordingly, we can rewrite Eq. (4.38) for v ∈ [−1, 0) as follows:

ẃ (t, v) = Aw (t, v) +G (z, z, v) , (4.39)

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 11

where

G (z, z, v) = G20 (v)
z2

2
+G11zz +G02

z2

2
+ . . .

So, differentiate the two sides of (4.35), gives that:

ẃ = wz ź + wz ź. (4.40)

Then, using the Eqs. (4.35), (4.36) and (4.39) give that:

(A− 2iτ0ω0)w20 (v) = −G20 (v) ,
Aw11 (v) = −G11 (v) ,

(A+ 2iτ0ω0)w02 (v) = −G02 (v) .
(4.41)

According to Eqs. (4.36) and (4.37), we have that:

g (z, z) = q∗
T
(0)H0(z, z) =

1

2
g20z

2 + g11zz +
1

2
g02z

2 +
1

2
g21z

2z + . . . , (4.42)

where q∗T (0) = D (1, r2).
Also from Eq. (4.34), it gets that:

x (t+ v) = (x1 (t+ v) , x2(t+ v))T = zq (v) + zq (v) + w (t, v) ,

here q (v) = (1, r1)
T eivτ0ω0 .

Hence by Eq. (4.35), it is obtained that:

x1 (t+ 0) = z + z + w(1) (t, 0)

= z + z + w
(1)
20 (0) z

2

2
+ w

(1)
11 (0) zz + w

(1)
02 (0) z

2

2
+ . . . ,

x2 (t+ 0) = zr1 + zr1 + w(2) (t, 0)

= zr1 + zr1 + w
(2)
20 (0) z

2

2
+ w

(2)
11 (0) zz + w

(2)
02 (0) z

2

2
+ . . . ,

x1 (t− 1) = ze−iτ0ω0 + zeiτ0ω0 + w(1) (t,−1)

= ze−iτ0ω0 + zeiτ0ω0 + w
(1)
20 (−1) z

2

2

+w
(1)
11 (−1) zz + w

(1)
02 (−1) z

2

2
+ . . . ,

x2 (t− 1) = zr1e
−iτ0ω0 + zr1e

iτ0ω0 + w2 (t,−1)

= zr1e
−iτ0ω0 + zr1e

iτ0ω0 + w
(2)
20 (−1) z

2

2

+w
(2)
11 (−1) zz + w

(2)
02 (−1) z

2

2
+ . . . .

(4.43)

Now, substituting Eq. (4.25), when u = 0 and Γ(v) ≡ x(t+ v) in Eq. (4.42), gives that:

g (z, z) = Dτ0 (1, r2)

(
−β (1−m)x1 (t+ 0)x2 (t+ 0)

0

)
.

Therefore, the following is obtained:

g (z, z) = −Dτ0β (1−m)
[
r1z

2 + (r1 + r1) zz + r1z
2

+
(
2w

(2)
11 (0) + w

(2)
20 (0) + 2r1w

(1)
11 (0) + r1w

(1)
20 (0)

)
z2z
2

+
(
2w

(2)
11 (0) + w

(2)
20 (0) + 2r1w

(1)
11 (0) + r1w

(1)
02 (0)

)
z z2

2
+ · · ·

]
 . (4.44)

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



12 A.A. MOHSEN, R.K. NAJI

So comparing the coefficients in Eq. (4.44) with those in Eq. (4.42) gives that:

g20 = −2Dτ0β (1−m) r1,
g11 = −Dτ0β (1−m) (r1 + r1) ,
g02 = −2Dτ0β (1−m) r1,

g21 = −Dτ0β (1−m)
(
2w

(2)
11 (0) + w

(2)
20 (0) + 2r1w

(1)
11 (0) + r1w

(1)
20 (0)

)
.

(4.45)

Since Eq. (4.45) contains w11 and w20, hence it needs to compute them.
Now, from Eqs. 4.38 and (4.39), for v ∈ [−1, 0), we have:

G (z, z, v) = −q∗T (0)H0 (z, z) q (v)− q∗T (0)H0 (z, z) q (v) .

Then by using Eqs. (4.36) and (4.37), it obtains that:

G (z, z, v) = −g (z, z) q (v)− g (z, z) q (v) . (4.46)

By comparing coefficients in both sides gives:

G20 (v) = −g20q (v)− g02q (v) ,
G11 (v) = −g11q (v)− g11q (v) ,

(4.47)

Now, substituting Eqs. (4.47) into Eqs. (4.41) respectively gives that:

w′
20 (v) = 2iτ0ω0w20 (v) + g20q (v) + g02q (v) ,

w′
11 (v) = g11q (v) + g11q (v) .

(4.48)

Now by solving Eqs. (4.48), it is easy to verify that the solutions are written respectively
as:

w20 (v) = −g20q(0)
iτ0ω0

eiτ0ω0v − g02q(0)
3iτ0ω0

e−iτ0ω0v +K1e
2iτ0ω0v,

w11 (v) =
g11q(0)
iτ0ω0

eiτ0ω0v − g11q(0)
iτ0ω0

e−iτ0ω0v +K2.
(4.49)

whereK1 =
(
K

(1)
1 , K2

1

)
andK2 =

(
K

(1)
2 , K2

2

)
are constants vectors to be determined in the

following.
From the definition of A when v = 0, and Eqs. (4.41), it is obtained that:∫ 0

−1

dη (v)w20 (v) = 2iτ0ω0w20 (0)−G20(0). (4.50)

∫ 0

−1

dη (v)w11 (v) = −G11(0). (4.51)

Also from Eqs. (4.38) and (4.39), it is easy to verify that:

G20 (0) = −g20q (0)− g02q (0) + 2τ0

[
r1β (1−m)

0

]
. (4.52)

G11 (0) = −g11q (0)− g11q (0) + 2τ0

[
β (1−m) (r1 + r1)

0

]
. (4.53)

Therefore by substituting Eqs. (4.49) and (4.52) into Eq. (4.50), and then using the results:

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 13

(
iτ0ω0I −

∫ 0

−1
eiτ0ω0vdη (v)

)
q (0) = 0,(

−iτ0ω0I −
∫ 0

−1
e−iτ0ω0vdη (v)

)
q (0) = 0.

(4.54)

It is obtained that:(
2iτ0ω0 I −

∫ 0

−1

e2iτ0ω0vdη (v)

)
K1 = 2τ0

[
r1β(1−m)

0

]
. (4.55)

Similarly, by substituting Eqs. (4.49) and (4.53) into Eq. (4.51), and then using Eqs.
(4.54), it is obtained that:

K2

∫ 0

−1

dη (v) = −2τ0

[
β (1−m) (r1 + r1)

0

]
. (4.56)

Now, using Eq. (4.28) with u = 0 into Eqs. (4.55) and (4.56), gives that:

[
2iω0 + β (1−m)E1 + µ β (1−m)S1

0 2iω0 + θe−2iω0τ0

]
K1 = 2

[
r1β(1−m)

0

]
. (4.57)

[
−β (1−m)E1 − µ −β (1−m)S1

0 −θ

]
K2 = −2

[
β (1−m) (r1 + r1)

0

]
. (4.58)

Accordingly, solving the linear systems given by Eqs. (4.57) and (4.58), gives that:

K
(1)
1 = 2r1β(1−m)[2iω0+θe−2iω0τ0 ]

[2iω0+β(1−m)E1+µ][2iω0+θe−2iω0τ0 ]
,

K
(2)
1 = 0.

(4.59)

And

K
(1)
2 = 2β(1−m)(r1+r1)

[β(1−m)E1+µ]
,

K
(2)
2 = 0.

(4.60)

Consequently, all the values of gij; i = 0, 1, 2; j = 0, 1, 2 given in Eqs. (4.45) can be
determined by the parameters and delay. Thus, it can calculate the following quantities:

C1 (0) =
i

2τ0ω0

(
g20g11 − 2|g11|2 − |g02|2

3

)
+ g21

2
,

µ2 =
−Re{C1(0)}
Re{λ′(τ0)} ,

β2 = 2Re {C1(0)} ,
T2 = − Im {C1(0)}+µ2Im {λ′(τ0)}

τ0ω0
.

(4.61)

These quantities are used to determine the direction of the Hopf bifurcation and stability of
bifurcated periodic solutions of system (2.3) at the critical value τ0 as shown in the following
theorem.

Theorem 4.2:
The direction of the Hopf bifurcation is determined by the sign of µ2 at the critical value τ0,
so that

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



14 A.A. MOHSEN, R.K. NAJI

1. The the Hopf bifurcation is supercritical if µ2 > 0 and the Hopf bifurcation is subcritical
if µ2 < 0.

2. The stability of bifurcated periodic solutions is determined by β2 so that the periodic
solutions are stable if β2 < 0 and unstable if β2 > 0.

3. The period of bifurcated periodic solutions is determined by β2 so that the period
increases if T2 > 0 and decreases if T2 < 0.

5. NUMERICAL SIMULATION

In this section, the obtained results are illustrated using numerical simulation. The following
set of hypothetical parameters set is adopted throughout this section.

ψ = 0.87, β = 0.010, m = 0.05, µ = 0.015,
P = 0, α = 0.3, θ = 0.2, τ = 7.7 < τ0 ∼= 7.85. (5.62)

All the obtained trajectories of the system (2.1) are drawn using the MATLAB of version
8. Starting from different sets of initial data, system (2.1) is solved numerically using the
parameters set given in Eq. (5.62) and then the obtained trajectories are drawn in Fig. (5.2).

Fig. 5.2. : The trajectories of the system (2.1) using data given by Eq. (5.62) approach to DFEP. (a) Trajectories
of Susceptible. (b) Trajectories of Cancer. (c) Trajectories of Environment. (d) 3D phase plot for globally

asymptotically stable DFEP.

Clearly, Fig. (5.2) illustrates that system (2.1) has a globally asymptotically stable DFEP
for that data (5.62). However, for the same data given by Eq. (5.62) with P = 0.6, it is
observed that, although that system (2.1) is solved with different sets of initial points, it has a
globally asymptotically stable EEP given by E1 = (17.9, 1.6, 3) as shown in Fig. (5.3) below.

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 15

Fig. 5.3. : The trajectories of the system (2.1) using data given by Eq. (5.62) with P = 0.6 approach to EEP. (a)
Trajectories of Susceptible. (b) Trajectories of Cancer. (c) Trajectories of Environment. (d) 3D phase plot for

globally asymptotically stable EEP.

Now, to show the effect of varying the body resistance rate (m) due to chromosomal
fixed against environment effect on the system behavior, the system is solved numerically for
different values, say m = 0.1, 0.7 respectively, keeping other parameters fixed as given in
Eq. (5.62) with P = 0.6, and then the solution of system (2.1) is drawn in Figs. (5.4a) and
(5.4b).

Fig. 5.4. : Time series of the trajectory of the system (2.1) using data given by Eq. (5.62) with P = 0.6 and
different values of m. (a) Trajectory approaches to EEP given by e1 = (18, 1.5, 3) when m = 0.1. (b) Trajectory

approaches to EEP given by e1 = (32, 0.92, 3) when m = 0.7.

Clearly, in Fig. (5.4) as m increases, the trajectory of the system (2.1) still approaches
asymptotically to the EEP point. In fact, it is observed that as m increases the number of
individuals with cancer decreases, and the number of susceptible increases without any effect
of environment.

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



16 A.A. MOHSEN, R.K. NAJI

Now, the dynamical behavior of the system (2.1) near the EEP point under the effect of
increasing the time delay is investigated. The system (2.1) is solved numerically for the set
of parameter values given by Eq. (5.62) with P = 0.6 and τ = 7.86 then the trajectory of the
system (2.1) is drawn in Figs. (5.5a-5.5e).

Fig. 5.5. : The existence of periodic solution near EEP of the system (2.1) for data given by equation (5.62)
with p = 0.6 and τ = 7.86. (a) Trajectories of Susceptible. (b) Trajectories of Cancer. (c) Trajectories of

Environment. (d) 3D periodic solution.

6. CONCLUSION

In this paper, a mathematical model that describes the spread of cancer in a polluted
environment incorporating delay in cleaning up the environment from the contaminated has
been proposed and studied. The properties of the solution are discussed. It is observed that the
proposed model has two equilibrium points namely DFEP and EEP. The stability analysis
of the model shows that the DFEP is globally asymptotically stable for all τ ≥ 0. While the
EEP is conditionally stable so that it’s globally asymptotically stable for all τ ∈ [0, τ0) and
there is the Hopf bifurcation at τ0.

However, it is an unstable point for τ > τ0 and periodic dynamics occurred. The stability
and direction of the periodic dynamics are also investigated analytically by finding the
normal form using the center manifold theory as well as numerically. It is observed
that for the hypothetical set of data given by Eq. (5.62) with p = 0.6 the EEP is still
globally asymptotically stable for 0 ≤ ô < τ0 ∼= 7.853. While the Hopf bifurcation occurs
and periodic solutions bifurcate near the EEP point as τ passes through the above critical
value τ0. On the other hand, all the quantities, which are given in Eqs. (4.61), are determined
for the periodic dynamics drawn in Fig. 5.5 as C1 (0) = 7.360− 77.426i, β2 = 14.72 > 0,
µ2 = −0.294 < 0 and T2 = 51.873 > 0. Therefore, according to theorem (4.2), the periodic

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)



STABILITY AND BIFURCATION OF A DELAY CANCER MODEL ... 17

resulting from the Hopf bifurcation is subcritical, unstable, and the size of the period
increases.

ACKNOWLEDGEMENTS

The authors are thankful to the referees for their valuable suggestions.

REFERENCES

1. Al-Tuwairqi, S. M., Al-Johani N., O., & Eman A. S., (2020). Modeling dynamics of
cancer virotherapy with immune response. Advances in Difference Equations, 438.

2. Chaffer, C. L. & Weinberg, R. A., (2011). A perspective on cancer cell metastasis.
Science Journal, 331, 1559–1564.

3. Cooke, K. L., & Driessche, V. D., (1996). Analysis of an SEIRS epidemic model with
two delays, J. math. Boil., 35, 240–260.

4. Hassard, B., Kazarino D., & Wan, Y., (1981). Theory and application of the Hopf
bifurcation, Cambridge University Press.

5. International Agency for Research on Cancer. (16 October 2019). [Online] Available:
https://gco.iarc.fr/tomorrow/home.

6. Lestari, D., Sari, E. R. & Arifah, H., (2019). Dynamical of a mathematical model of
cancer cells with chemotherapy. Journal of Physics Conf. Series, 1320.

7. Naji, R. K. & Majeed, S. J., (2020). The dynamical Analysis of a Delay Pray-Predator
Model with a Refuge-Stage Structure Prey Population. Iranian Journal of Mathematical
Sciences ad Informatics, 15(1), pp 135–159.

8. Pillis, L. C. & Radunskaya, A., (2003). A mathematical model of immune response
to tumor invasion, Computational fluid and solid mechanics, pp. 1661–1668, Elsevier
science, Amsterdam.

9. Simmons, A., Burrage, P. M., Nicolau, D.V., Lakhani S.R. & Burrage, K., (2017).
Environmental factors in breast cancer invasion: a mathematical modeling review.
Pathology, 49(2), 172–180.

10. Tipsri, S. & Chinviriyasit W., (2015). The effect of time delay on the dynamics of an
SEIR model with nonlinear incidence. Chaos Solitons & Fractals 75, 153–172.

11. Trisilowati, T., McCue, S., & Mallet, D., (2013). Numerical solution of an optimal
control model of dendritic cell treatment of growing tumour, ANZIAM journal, 54, 664–
680.

12. Tsygvintsev, A., Marino, S. & Kirschner, D. E., (2012). A mathematical model of
gene therapy for the treatment of cancer, in Mathematical Models and Methods in
Biomedicine, 357–373. Springer-Verlag.

13. Wang, L., & Wu, X., (2018). Stability and the Hopf bifurcation for an SEIR epidemic
model with delay, Advances in the theory of nonlinear analysis and its applications,
113–127.

14. Weerasinghe, H. N., Burrage, P. M., Burrage, K. & Nicolau D.V., (2019). Mathematical
models of cancer call plasticity. Journal of Oncology, 2019, 2403483.

15. Wilson, S. & Levy, D., (2012). A mathematical model of enhancement of tumor vaccine
efficacy by immunotherapy, Bulletin of mathematical biology, 74(7), 1485–1500.

16. Zachary, A., Kristen A., & Jessica S., (2020). A mathematical model for tumor growth
and treatment using virotherapy. AIMS Mathematics, 5(5), 4136–4150.

17. Zuo, L., Liu, M. & Wang, J., (2015). The impact of awareness program with recruitment
and delay on the spread of an epidemic, Mathematical problems in engineering, 235935.

Copyright © 2022 ASSA. Adv Syst Sci Appl (2022)


	Introduction
	THE MATHEMATICAL FORMULATION OF CANCER MODEL 
	ANALYSIS OF SYSTEM 2.3
	STABILITY AND DIRECTION OF the Hopf BIFURCATION
	NUMERICAL SIMULATION
	Conclusion