370 

© 2024 The Author(s). Published by College of Education for Pure Science (Ibn Al-Haitham), 

University of Baghdad. This is an open-access article distributed under the terms of the Creative Commons 

Attribution 4.0 International License 

 

 

Ibn Al-Haitham Journal for Pure and Applied Sciences 

Journal homepage: jih.uobaghdad.edu.iq 

PISSN: 1609-4042, EISSN: 2521-3407 

 IHJPAS. 2024, 37(4) 

 

Stability and Bifurcation Analysis of the Impact of Refuge on the Dissolved 

Oxygen-Plankton Interaction 

Ahmed Ali1,* , Shireen Jawad2  and Fengde Chen3  

 
1,2Department of Mathematics, College of Science, University of Baghdad, Baghdad, Iraq. 

3College of Mathematics and Statistics, Fuzhou University, Fuzhou 350108, P.R. China. 

*Corresponding Author.  

 

Received: 21 May 2023                          Accepted: 9 July 2023                    Published: 20 October 2024  

doi.org/10.30526/37.4.3505 

 
Abstract  

The suggested mathematical model for studying the effect of refuge on the dissolved oxygen in 

the plankton ecosystem is based on measurements of dissolved oxygen, phytoplankton, and 

zooplankton populations. The aim of this work is to find out the potential equilibrium and to 

investigate their behaviour. The study shows that there are three points of equilibrium. The 

feasibility requirements and stability conditions for all steady states are determined. Using the 

consumption of oxygen by zooplankton as a bifurcation parameter, we test for the presence of 

Hopf- bifurcation for the interior equilibrium. It is shown what conditions must be met for stable 

limit cycles. Finally, a numerical simulation is conducted to back up the analytical findings. It 

shows when the stability criteria are met, the solution of the proposed system constantly oscillates 

around the positive stable state. In addition, the solution exhibits limit cycle behaviour for small 

changes in certain parameters. 

Keywords: Plankton interaction, Stability analysis, Bifurcation, Dissolved oxygen model. 
 

1. Introduction 

          Much effort has been put into understanding dissolved oxygen dynamics better since they 

are such a crucial indicator of the health of marine ecosystems [1-3]. Most of the oxygen in the 

oceans comes from phytoplankton which also serves as the foundation of the marine food chain 

through their photosynthesis [4]. It’s well knowledge that salinity, temperature, and nutrient 

availability all play significant roles in determining how much oxygen phytoplankton can create. 

Additionally, there is a considerable diurnal variation in oxygen generation by phytoplankton. 

Therefore, the link between phytoplankton and dissolved oxygen is crucial to the existence of 

organisms. Oxygen production fluctuations can have severe consequences for marine life. [5]. 

Since oxygen is used by living things for photosynthesis (during the day) and respiration (at night), 

the amount of oxygen in the water changes day and night. As a result, phytoplankton colonies can 

serve as reliable markers of ecological status. [6], [7]. Mondal,  and his colleague, for instance, 

have studied how the coupled plankton-oxygen dynamics in the ocean are affected by a low oxygen 

production rate, which can result in oxygen depletion and species extinction [8]. 

https://creativecommons.org/licenses/by/4.0/
https://creativecommons.org/licenses/by/4.0/
https://orcid.org/0009-0004-7408-6238
mailto:ahmedali.970q6@gmail.com
https://orcid.org/0000-0002-3090-8357
mailto:shireen.jawad@sc.uobaghdad.edu.iq
https://orcid.org/0000-0003-3617-5550
mailto:fdchen@263.ne


IHJPAS. 2024, 37(4) 

371 

 

Furthermore, the primary objective of the study of theoretical ecology is to identify the various 

dynamical mechanisms underlying interactions between prey and predator [9-11]. The relationship 

between phytoplankton and zooplankton is an example of a predator-prey interaction that reveals 

numerous aspects of marine ecology. Phytoplankton contributes substantially to aquatic 

ecosystems, including producing vast quantities of oxygen, managing natural resources and water 

quality, and establishing multiple food webs [12-14]. Plankton dynamics research is a fascinating 

field of study. Plankton constitutes the building elements of all aquatic food chains, with 

phytoplankton occupying the first trophic level [15]. Bagheli and Dhar examine the effect of 

dissolved oxygen on the presence of a planktonic population that interacts. They conclude that 

Hopf-bifurcation in the interior equilibrium is possible if the phytoplankton growth rate is selected 

as the bifurcation parameter [16]. 

The purpose of this research is to examine how the phytoplankton refuge affects the dynamics of 

the oxygen-plankton model. Some phytoplankton may evade their zooplankton prey by relocating 

to deeper water column layers. These sediments offer shelter from predators to the prey. Here is 

how the article is laid out: in Sec. 2, we built the structure of the proposed model. Sec. 3 explains 

the feasibility requirements and stability conditions for all steady states. The prevalence of Hopf 

bifurcations is also illustrated in Sec. 4. In Sec. 5. we undertake MATLAB-based numerical 

simulations to validate the analytical results. 

 

2. Construction of the Model  

           Include Let 𝑢(𝑡) and 𝑣(𝑡) indicate the phytoplankton and zooplankton populations at the 

time 𝑡, respectively. We assume some phytoplankton populations are safe from zooplankton 

predation because they can conceal themselves in the ocean’s floor-diverse sediments. These 

sediments offer refuge from predators to the hunted. 𝑤(𝑡) represents the dissolved oxygen 

concentration in the marine. Since phytoplankton do photosynthesis throughout the day, they also 

contribute to atmospheric oxygen production. Several additional factors, such as the respiration of 

marine organisms, the consumption of oxygen by phytoplankton at night, and the gradual drop in 

oxygen concentration brought about by chemical reactions in the water, all affect the rate at which 

oxygen is depleted. The following set of ordinary differential equations serves as the governing 

structure for the dynamical system of the system (1): 

 𝑑𝑢

𝑑𝑡
=

𝑟𝑢

(𝑎1 + 𝑤0 −𝑤)
− 𝛼1𝑢(1 − 𝑚)𝑣 − 𝛿1𝑢, 

𝑑𝑣

𝑑𝑡
=
𝛼2𝑢(1 −𝑚)𝑣

(𝑎2 + 𝑤0 − 𝑤)
− 𝛿2𝑣, 

𝑑𝑤

𝑑𝑡
= 𝑠(𝑤0 − 𝑤) + 𝑑𝑢 − 𝛾𝑤 − 𝛾1𝑢𝑤 − 𝛾2𝑣𝑤, 

 

 

                    (1) 

 

with the initial conditions 𝑢(0) ≥ 0, 𝑣(0) ≥ 0 and 𝑤(0) ≥ 0. In the first equation of the system 

(1), 
𝑟𝑢

(𝑎1+𝑤0−𝑤)
  represents the absorption of dissolved oxygen from phytoplankton with the growth 

rate 𝑟. The maximum growth rate of the phytoplankton population is 𝑟/𝑎1 at 𝑤 = 𝑤0. 𝛼1 is the 

phytoplankton’s capture rate by zooplankton. 𝑚 ∈ (0,1) is the proportion of protected 

phytoplankton. (1 − 𝑚) is the ratio of unprotected phytoplankton devoured by different 

zooplankton groups. 𝛼2 is the conversion rate from phytoplankton to zooplankton. 𝛿1 and 𝛿2 are 

the phytoplankton and zooplankton’s natural death rates. 𝑎1 is the phytoplankton saturation 

constant. 𝑎2 is the zooplankton saturation constant. 𝑤0 is the constant concentration of dissolved 



IHJPAS. 2024, 37(4) 

372 

 

oxygen that comes from external sources. 𝑠 is the replenishment rate of oxygen in marine. 𝑑 is the 

amount of oxygen produced as a result of the process of photosynthesis carried out by 

phytoplankton. 𝛾 is the natural depletion rate of oxygen. 𝛾1 is the consumption of oxygen by 

phytoplankton during the night. 𝛾2 is the consumption of oxygen by zooplankton.  

The equations of the system (1) are 𝐶1(𝑅+
3), where 𝑅+

3 = {(𝑢, 𝑣, 𝑤), 𝑢 ≥ 0, 𝑣 ≥ 0, 𝑤 ≥ 0}. 

Consequently, they are Lipschitz Ian [17]. Therefore, the system’s (1) solution exists and is unique.  

 

2. Existence of equilibria 

System (1) has three non-negative steady states, namely 

1. 𝑧1 = (0,0, �̂�), where �̂� =
𝑠𝑤0

𝑠+𝛾
.  

2. 𝑧2 = (�̅�, 0, �̅�), where �̅� =
(𝑠+𝛾)(𝑎1𝛿1−𝑟)

𝑑𝛿1−𝛾1(𝛿1(𝑎1+𝑤0)−𝑟)
 and �̅� = 𝑎1 + 𝑤0 −

𝑟

𝛿1
, and. Clearly, �̅� and �̅� 

are positive if the following condition is satisfied: 

 𝑟 < 𝑚𝑖𝑛. {𝑎1𝛿1, 𝛿1(𝑑 − 𝛾1(𝑎1 + 𝑤0))}. (2) 

3. 𝑧3 = (𝑢
∗, 𝑣∗, 𝑤∗), where 𝑢 =

𝛿2(𝑎2+𝑤0−𝑤)

(1−𝑚)[𝛼2−𝑎(𝑎2+𝑤0−𝑤)]
, 𝑣 =

𝑟

𝛼1(1−𝑚)(𝑎1+𝑤0−𝑤)
−

𝛿1

𝛼1(1−𝑚)
, and 

𝑤 is the root of the following equation: 

𝐵0𝑤
3 + 𝐵1𝑤

2 + 𝐵2𝑤 + 𝐵3 = 0,                                                                                                  (3) 

where, 

𝐵0 = 𝑎𝛼1(1 − 𝑚)(𝑠 + 𝛾) > 0, 

𝐵1 = 𝛼1𝑎𝑠𝑤0[2 − (1 − 𝑚)] + 𝑑𝛿2𝛼1 − 𝛼1(𝑠 + 𝛾)(𝛼2 − 𝑎𝑎2) − 𝛼1𝑎𝛾[𝑎1(1 − 𝑚) − 2𝑤0] 

𝐵2 = 𝛼1𝑎𝑠𝑤0(1 − 𝑚)(2𝑎1 + 3𝑤0) − 𝛼1𝑤0(1 − 𝑚)(𝑠 + 𝛾)[𝛼2 − 𝑎𝑎2] − 𝑑𝛿2𝛼1𝑎1 + 𝛼1(1

− 𝑚)(𝑠 + 𝛾)[𝑎𝑎1𝑎2 − 𝛼2𝑎1 + 𝑎𝑤0
2] + 𝛼1𝑎𝛾𝑤0(1 − 𝑚)(𝑎1 + 𝑤0) 

𝐵3 = 𝛼1𝑠𝑤0(1 − 𝑚)[𝛼2𝑎1 − 𝑎𝑎1𝑎2 − 𝑎𝑎1𝑤0 + 𝛼2𝑤0 − 𝑎𝛼2𝑤0 − 𝑎𝑤0
2] + 𝑑𝛿2𝛼1[𝑎1𝑎2 +

𝑎1𝑤0 + 𝑎2𝑤0 + 𝑤0
2]. 

Using Descartes’s rule of sign [15], equation (3) has a unique positive root, say 𝑤 = 𝑤∗, if one of 

the following sets conditions hold: 

 𝐵1 > 0 and 𝐵3 < 0, 

𝐵2 < 0 and 𝐵3 < 0. 

(4) 

 

For 𝑢∗and 𝑣∗to be positive, the following two conditions must be satisfied: 

 𝛼2 > 𝑎(𝑎2 +𝑤0 − 𝑤), 

𝑟 > 𝛿1(𝑎1 +𝑤0 − 𝑤). 

(5) 

 

3. Stability Analysis  

The feature of the eigenvalues of the Jacobian matrix 𝐽(𝑢, 𝑣, 𝑤) at an equilibrium point is directly 

related to the behaviour of the system (1) near an equilibrium. The 𝐽(𝑢, 𝑣, 𝑤) at any point, say 

(𝑢, 𝑣, 𝑤), can be written as: 

𝐽 =

[
 
 
 
 

𝑟

(𝑎1 +𝑤0 − 𝑤)
− 𝛼1𝑣(1 −𝑚) − 𝛿1 −𝛼1𝑢(1 − 𝑚)

𝑟𝑢

(𝑎1 + 𝑤0 − 𝑤)
2

𝛼2𝑣(1 − 𝑚)

(𝑎2 + 𝑤0 − 𝑤)

𝛼2𝑢(1 − 𝑚)

(𝑎2 + 𝑤0 − 𝑤)
− 𝛿2

𝛼2𝑢(1 − 𝑚)𝑣

(𝑎2 +𝑤0 − 𝑤)
2

𝑑 − 𝛾1𝑤 −𝛾2𝑤 −(𝑠 + 𝛾 + 𝛾1𝑢 + 𝛾2𝑣)]
 
 
 
 

 

The local stability of system (1) around each equilibrium is: 

 



IHJPAS. 2024, 37(4) 

373 

 

 

1. The Jacobian matrix at 𝑧1 = (0,0, �̂�) is given as: 

 

𝐽(𝑧1) = [

𝑟

(𝑎1 + 𝑤0 − �̂�)
− 𝛿1 0 0

0 −𝛿2 0
𝑑 − 𝛾1�̂� −𝛾2�̂� −𝑠 − 𝛾

] 

 

 

Then,  𝐽(𝑧1) has the eigenvalues 𝜆11 =
𝑟

(𝑎1+𝑤0−�̂�)
− 𝛿1,  𝜆12 = −𝛿2 < 0, and 𝜆13 = −𝑠 − 𝛾. 𝑧1is a 

locally asymptotically stable point if and only if 

 𝑟 < 𝛿1(𝑎1 +𝑤0 − �̂�) (6) 

 

2. The Jacobian matrix at 𝑧2 = (�̅�, 0, �̅�) is given as: 

𝐽(𝑧2) =

[
 
 
 
 
 

𝑟

(𝑎1 +𝑤0 − �̅�)
− 𝛿1 −𝛼1�̅�(1 − 𝑚)

𝑟�̅�

(𝑎1 + 𝑤0 − �̅�)
2

0
𝛼2�̅�(1 − 𝑚)

(𝑎2 +𝑤0 − �̅�)
− 𝛿2 0

𝑑 − 𝛾1�̅� −𝛾2�̅� −𝑠 − 𝛾 − 𝛾1�̅� ]
 
 
 
 
 

 

Then, |𝐽(𝑧2) − 𝐼𝜆| = 0 gives: 

(
𝛼2�̅�(1 − 𝑚)

(𝑎2 + 𝑤0 − �̅�)
− 𝛿2 − 𝜆) [𝜆

2 − 𝑇𝑟(𝐽(𝑧2))𝜆 + 𝐷𝑒𝑡(𝐽(𝑧2))] 

 The eigenvalues of the above equation can be written as follows 

 𝜆21 =
𝛼2�̅�(1−𝑚)

(𝑎2+𝑤0−�̅�)
− 𝛿2, 

𝑇𝑟(𝐽(𝑧2)) =
𝑟

(𝑎1+𝑤0−�̅�)
− 𝛿1 − (𝑠 + 𝛾 + 𝛾1�̅�),  

 𝐷𝑒𝑡(𝐽(𝑧2)) = 𝛿1(𝑠 + 𝛾 + 𝛾1�̅�) +
𝑟𝑢 ̅ (𝛾1�̅�−𝑑)

(𝑎1+𝑤0−�̅�)
2 −

𝑟(𝑠+𝛾+𝛾1�̅�)

(𝑎1+𝑤0−�̅�)
. 

Clearly, 𝑧2 is a locally asymptotical stable point if and only if the following conditions are satisfied: 

𝛿2 >
𝛼2�̅�(1 − 𝑚)

(𝑎2 + 𝑤0 − �̅�)
,

𝑟 < [𝛿2 + (𝑠 + 𝛾 + 𝛾1�̅�)](𝑎1 + 𝑤0 − �̅�),

𝛿1(𝑠 + 𝛾 + 𝛾1�̅�) +
𝑟𝑢 ̅ (𝛾1�̅� − 𝑑)

(𝑎1 + 𝑤0 − �̅�)
2
>
𝑟(𝑠 + 𝛾 + 𝛾1�̅�)

(𝑎1 + 𝑤0 − �̅�)
.
}
 
 

 
 

 

 

 

 

   (7) 

3. The Jacobian matrix at 𝑧3 = (𝑢
∗, 𝑣∗, 𝑤∗) is given as: 

𝐽(𝑧3) =

[
 
 
 
 

𝑟

(𝑎1+𝑤0−𝑤
∗)
− 𝛿1 − 𝛼1𝑣

∗(1 − 𝑚) −𝛼1𝑢
∗(1 − 𝑚)

𝑟𝑢∗

(𝑎1+𝑤0−𝑤
∗)2

𝛼2𝑣
∗(1−𝑚)

(𝑎2+𝑤0−𝑤
∗)

0
𝛼2𝑢

∗(1−𝑚)𝑣∗

(𝑎2+𝑤0−𝑤
∗)2

𝑑−𝛾1𝑤
∗ −𝛾2𝑤

∗ −𝑠 − 𝛾 − 𝛾1𝑢
∗ − 𝛾2𝑣

∗]
 
 
 
 

=

(𝑎𝑖𝑗)3×3                   (8) 

So, the characteristic equation of 𝐽(𝑧3) can be written as: 

 𝜆3 + 𝐴1𝜆
2 + 𝐴2𝜆 + 𝐴3 = 0, (9) 

where,   

𝐴1 = −(𝑎11 + 𝑎33), 𝐴2 = −(𝑎13𝑎31 + 𝑎23𝑎32 + 𝑎12𝑎21 − 𝑎11𝑎33), 

𝐴3 = 𝑎11𝑎23𝑎32 + 𝑎12𝑎21𝑎33 − 𝑎13𝑎21𝑎32 − 𝑎12𝑎23𝑎31, 

∆= 𝐴1𝐴2 − 𝐴3 = (𝑎11 + 𝑎33)(𝑎13𝑎31 − 𝑎11𝑎33) + 𝑎11𝑎12𝑎21 + 𝑎23𝑎32𝑎33 + 𝑎12𝑎23𝑎31 +



IHJPAS. 2024, 37(4) 

374 

 

𝑎13𝑎21𝑎32. 

Now, from the Routh-Hurwitz criteria [18], 𝑧3 is a LAS  point, under the condition that 𝐴1 > 0, 𝐴3 >

0 and ∆> 0.  

 

4. Hop Bifurcation    

           From Theorem 2, the steady state 𝑧3 changes as the parameter 𝛾2 crosses the threshold value 

𝛾2
∗, which implies that 𝑧3 may become unstable due to Hopf bifurcation when forced to operate within 

particular restrictions on its parameters [19-25]. In the case where we use 𝛾2
∗ as the bifurcation 

parameter, the Hopf bifurcation threshold and its conditions are clearly clarified in the following 

theorem []. 

Theorem 2. Under the following assumptions 

𝐴𝑖 > 0, 𝑖 = 1,2
𝑤∗ > 𝑣∗

𝛾2
∗ > 0

 

(11) 

(12) 

(13) 

 

where Ai’s are the coefficients of the characteristic equation given in equation (9) with 𝛾2 = 𝛾2
∗ and 

the formula for 𝛾2
∗ is shown in the following proof. Then, there exists a Hopf bifurcation for 𝑧3 at 𝛾2 =

𝛾2
∗. 

Proof: - The value of the bifurcation parameter can be found if we set 𝐴1(𝛾2
∗)𝐴2(𝛾2

∗) − 𝐴3(𝛾2
∗) = 0 

in equation (9). This gives: 

𝛾2
∗ =

( 𝑎11+𝑎33)(𝑎13𝑎31−𝑎11𝑎33)+𝑎11𝑎12𝑎21+𝑎12𝑎23𝑎31)

(𝑎23𝑎33+𝑎13𝑎21)𝑤
∗ . 

Clearly, 𝛾2
∗ > 0 if condition (13) holds. Now, at 𝛾2 = 𝛾2

∗, equation (9) can be written as 

(𝜆 + 𝐴1)(𝜆
2 + 𝐴2) = 0.  

According to condition (11), the above equation has three roots, a negative root 𝜆1 = −𝐴1 and two 

purely imaginary roots  𝜆2,3 = ±𝑖√𝐴2. In a neighbourhood of 𝛾2
∗, the roots have the following forms: 

𝜆1 = −𝐴1, 𝜆2,3 = 𝜌1(𝛾2) ± 𝑖𝜌2(𝛾2). 

Clearly, 𝑅𝑒( 𝜆2,3)|𝛾2=𝛾2∗ = 𝜌1(𝛾2
∗) = 0 indicates that the first condition for Hopf bifurcation has been 

met at 𝛾2 = 𝛾2
∗. Now to confirm the transversality condition, we substitute 𝜌1(𝛾2) ± 𝑖𝜌2(𝛾2) into 

equation (9) and then compute its derivative with respect to 𝑑∗, 𝛩(𝛾2
∗)𝜓(𝛾2

∗) + 𝛤(𝛾2
∗)𝜙(𝛾2

∗) ≠ 0, 

where the form of 𝛩(𝛾2
∗), 𝜓(𝛾2

∗), 𝛤(𝛾2
∗) and 𝜙(𝛾2

∗) are 

𝜓(𝛾2) = 3𝜌1
2(𝛾2) + 2𝐴1(𝛾2)𝜌1(𝛾2) + 𝐴2(𝛾2) − 3𝜌2

2(𝛾2), 

𝜙(𝛾2) = 6𝜌1(𝛾2)𝜌2(𝛾2) + 2𝐴1(𝛾2)𝜌2(𝛾2), 

𝛩(𝛾2) = 𝜌1
2(𝛾2)𝐴1 

′(𝛾2) + 𝐴2 
′(𝛾2)𝜌1(𝛾2) + 𝐴3 

′(𝛾2) − 𝐴1 
′(𝛾2)𝜌2

2(𝛾2), 

𝛤(𝛾2) = 2𝜌1(𝛾2)𝜌2(𝛾2)𝐴1
′ (𝛾2) + 𝐴2

′ (𝛾2)𝜌2(𝛾2). 

Now at𝛾2 = 𝛾2
∗, substitution 𝜌1 = 0 and 𝜌2 = √𝐴2, into equation (9), the following is obtained: 

 

𝜓(𝛾2
∗) = −2𝐴2(𝛾2

∗),

𝜙(𝛾2
∗) = 2𝐴1(𝛾2

∗)√𝐴2(𝛾2
∗),

𝛩(𝛾2
∗) = 𝐴3

′ (𝛾2
∗) − 𝐴1

′ (𝛾2
∗)𝐴2(𝛾2

∗),

𝛤(𝛾2
∗) = 𝐴2

′ (𝛾2
∗)√𝐴2(𝛾2

∗),

 

where 

 𝐴1
′ (𝛾2

∗) = 𝑣∗,  



IHJPAS. 2024, 37(4) 

375 

 

𝐴2
′ (𝛾2

∗) = 𝑤∗ − 𝑣∗, 

𝐴3
′ (𝛾2

∗) = −2𝑤∗ − 𝑣∗. 
 

Hence, condition (12) gives   

𝛩(𝛾2
∗)𝜓(𝛾2

∗) + 𝛤(𝛾2
∗)𝜙(𝛾2

∗) = 2𝐴2(𝛾2
∗)[𝑣∗ + 2𝑤∗ + 2𝑣∗𝐴2(𝛾2

∗) + (𝑤∗ − 𝑣∗)2𝐴1(𝛾2
∗)𝐴2(𝛾2

∗)]

≠ 0. 

That means the Hop bifurcation has occurred at 𝛾2
∗. 

From Theorem 3, the stability condition of the stable limit cycle in 𝑅(𝑢,𝑣,𝑤)
3  is presented using the 

coefficient of curvature of the limit cycle. For a detailed discussion, we refer to  [18].  

 

Theorem 3 The system (1) has a stable limit cycle in 𝑅(𝑢,𝑣,𝑤)
3 , if the following conditions are true: 

𝑟

(𝑎1+𝑤0−𝑢3−𝑤
∗)2
≠

𝛼2(𝑢1+𝑢
∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2

. (14) 

 

Proof: - by shifting the 𝑧3 = (𝑢
∗, 𝑣∗, 𝑤∗) to (0, 0, 0) by using the following transformations 𝑢 = 𝑢1 +

𝑢∗, 𝑣 = 𝑢2 + 𝑣
∗, 𝑤 = 𝑢3 +𝑤

∗. Then the system (1) becomes: 

𝑑𝑢1
𝑑𝑡

=
𝑟(𝑢1 + 𝑢

∗)

(𝑎1 + 𝑤0 − 𝑢3 − 𝑤
∗)
− 𝛼1(𝑢1 + 𝑢

∗)(𝑢2 + 𝑣
∗)(1 − 𝑚) − (𝑢1 + 𝑢

∗)𝛿1 

𝑑𝑢2
𝑑𝑡

=
𝛼2(𝑢1 + 𝑢

∗)(𝑢2 + 𝑣
∗)(1 − 𝑚)

(𝑎2 + 𝑤0 − 𝑢3 − 𝑤
∗)

− 𝛿2(𝑢2 + 𝑣
∗) 

𝑑𝑢3
𝑑𝑡

= 𝑠[𝑤0 − (𝑢3 + 𝑤
∗)] + 𝑑(𝑢1 + 𝑢

∗) − 𝛾(𝑢3 + 𝑤
∗) − 𝛾1(𝑢1 + 𝑢

∗)(𝑢3 + 𝑤
∗)

− 𝛾2(𝑢2 + 𝑣
∗)(𝑢3 +𝑤

∗), 

where the nonlinear part of the above system is presented in the following matrix is 

℧ = (
℧1
℧2
℧3

) =

(

  
 

𝑟(𝑢1 + 𝑢
∗)

(𝑎1 + 𝑤0 − 𝑢3 − 𝑤
∗)
− 𝛼1(1 − 𝑚)𝑢1𝑢2

𝛼2(𝑢1 + 𝑢
∗)(𝑢2 + 𝑣

∗)(1 − 𝑚)

(𝑎2 + 𝑤0 − 𝑢3 − 𝑤
∗)

−𝛾1𝑢1𝑢3 − 𝛾2𝑢2𝑢3 )

  
 

 

We derive the following characteristic quantities from the nonlinear part: 

𝑔20
0 =

1

4
{
𝜕2℧1

𝜕𝑢1
2 −

𝜕2℧1

𝜕𝑢2
2 + 2

𝜕2℧2

𝜕𝑢1𝜕𝑢2
+ 𝑖 (

𝜕2℧2

𝜕𝑢1
2 −

𝜕2℧2

𝜕𝑢2
2 − 2

𝜕2℧1

𝜕𝑢1𝜕𝑢2
)} =

1

2
{

𝛼2(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)
− 𝛼1(1 −

𝑚)𝑖}, 

𝑔11
0 =

1

4
{
𝜕2℧1

𝜕𝑢1
2 +

𝜕2℧1

𝜕𝑢2
2 + 𝑖 (

𝜕2℧2

𝜕𝑢1
2 +

𝜕2℧2

𝜕𝑢2
2 )} = 0, 

𝐺110
0 =

1

2
{
𝜕2℧1

𝜕𝑢1𝜕𝑢3
+

𝜕2℧2

𝜕𝑢2𝜕𝑢3
+ 𝑖 (

𝜕2℧2

𝜕𝑢1𝜕𝑢3
−

𝜕2℧1

𝜕𝑢2𝜕𝑢3
)} =

1

2
{

𝑟

(𝑎1+𝑤0−𝑢3−𝑤
∗)2
+

𝛼2(𝑢1+𝑢
∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2
+

𝑖 (
𝛼2(𝑢2+𝑣

∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2
)}, 

𝐺101
0 =

1

2
{
𝜕2℧1

𝜕𝑢1𝜕𝑢3
−

𝜕2℧2

𝜕𝑢2𝜕𝑢3
+ 𝑖 (

𝜕2℧2

𝜕𝑢1𝜕𝑢3
+

𝜕2℧1

𝜕𝑢2𝜕𝑢3
)} =

1

2
{

𝑟

(𝑎1+𝑤0−𝑢3−𝑤
∗)2
−

𝛼2(𝑢1+𝑢
∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2
+

𝑖 (
𝛼2(𝑢2+𝑣

∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2
)}, 

𝑊11
0 = −

1

4𝜆3(𝑎1(𝑘
∗)
(
𝜕2℧3

𝜕𝑢1
2 +

𝜕2℧3

𝜕𝑢2
2 ) = 0, 

𝑊20
0 = −

1

4𝜆3(𝑎1(𝑘
∗)
(
𝜕2℧3

𝜕𝑢1
2 +

𝜕2℧3

𝜕𝑢2
2 − 2𝑖

𝜕2℧3

𝜕𝑢1𝜕𝑢2
) = 0, 



IHJPAS. 2024, 37(4) 

376 

 

𝐺21
0 =

1

8
{
𝜕3℧1

𝜕𝑢1
3 +

𝜕3℧1

𝜕𝑢1𝜕𝑢2
2 +

𝜕3℧2

𝜕𝑢2
3 +

𝜕3℧2

𝜕𝑢1
2𝜕𝑢2

+ 𝑖 (
𝜕3℧2

𝜕𝑢1
3 +

𝜕3℧2

𝜕𝑢1𝜕𝑢2
2 −

𝜕3℧1

𝜕𝑢2
3 −

𝜕3℧1

𝜕𝑢1
2𝜕𝑢2

)} = 0, 

Thus, the coefficient of the curvature of the limit cycle of the DOPZ system (1) is given by 

𝜎1
0 = 𝑅𝑒 {

𝑔20
0 𝑔11

0

4
𝑖 + 𝐺110

0 𝑊11
0 +

𝐺21
0 +𝐺101

0 𝑊20
0

2
}, 

𝜎1
0 = 𝑅𝑒

1

4
{

𝑟

(𝑎1+𝑤0−𝑢3−𝑤
∗)2
−

𝛼2(𝑢1+𝑢
∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2
+ 𝑖 (

𝛼2(𝑢2+𝑣
∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2
)} =

𝑟

(𝑎1+𝑤0−𝑢3−𝑤
∗)2
−

𝛼2(𝑢1+𝑢
∗)(1−𝑚)

(𝑎2+𝑤0−𝑢3−𝑤
∗)2

. 

Thus, Condition (14) guarantees that system (1) has a stable limit cycle. 

 

5. Numerical Simulations and Discussion  

          Numerical simulations support our theoretical predictions and reveal the system’s numerous 

dynamics (1). The ode45 solver was used to find the numerical solution to our system, and all figures 

were made in MATLAB 2019b. We aim to study the kinetics of dissolved oxygen depletion for the 

phytoplankton-zooplankton interaction with the following data: 

 

𝑟 = 0.35, 𝛼1 = 0.26, 𝛼2 = 0.17, 𝑎1 = 0.21, 𝑎2 = 0.21, 𝛿1 = 0.11, 𝛿2 = 0.11,𝑚 =

0.25, 𝛾 = 0.21, 𝛾1 = 0.19, 𝛾2 = 0.41,𝑤0 = 3, 𝑠 = 2.86, 𝑑 = 0.41, 

   (15) 

 

To examine the effect of varying 𝛾2 (the consumption of oxygen by zooplankton), system (1) has been 

numerically solved for the data in (15) with different values. It is clear from Figure 1 the solution 

converges to 𝑧3 for  𝛾2 > 0.22. Further, the solution approaches a periodic behaviour for 𝛾2 ≤ 0.22. 

The latter result confirms the one obtained in Theorems 2, which establishes the existence of Hopf 

bifurcation at 𝛾2 = 0.22. 
 
 

 
Figure 1.  Dynamics of the system (1) (a) time series with  γ2 = 0.41; (b) phase portrait of (a); (c) time series with γ2 =

0.22; (d) phase portrait of (c). 

 

Further, Figure 2 investigates the effect of change in the proportion of protected phytoplankton (𝑚) 

on the stability properties of the system (1). It shows for 𝑚 < 0.65, and the solution settles to the 



IHJPAS. 2024, 37(4) 

377 

 

positive equilibrium point. While for 𝑚 ≥ 0.65, the solution delivers a periodic attractor behaviour.  

 
Figure 2. Dynamics of the system (1) (a) time series with  𝑚 = 0.15; (b) phase portrait of (a); (c) time series with m =

0.65; (d) phase portrait of (c). 

 

Further, Figure 3 shows for different values of 𝛾 (the natural depletion rate of oxygen), the solution 

stabilizes at  𝑧3 for  𝛾 < 0.62. While for 𝛾 ≥ 0.62, the solution shows a periodic behaviour.  
 
 

 
Figure 3. Dynamics of the system (1) (a) time series with  γ = 0.21; (b) phase portrait of (a); (c) time series with γ =

0.62; (d) phase portrait of (c). 

 

Further, Figure 4 investigates the effect of change in the replenishment rate of oxygen in the marine 

(𝑠) on the stability properties of the system (1). It shows for 𝑠 > 1.85; the solution settles down to 𝑧3. 

While for 𝑠 ≤ 1.85., the solution shows a periodic behaviour. 
 



IHJPAS. 2024, 37(4) 

378 

 

 
Figure 4. Dynamics of system (1) (a) time series with s=2.86; (b) phase portrait of (a); (c) time series with s=1.85; (d) 

phase portrait of (c); (e) time series with s=0.44; (f) phase portrait of (e). 

 

 Now the effect of changing the concentration of dissolved oxygen that comes from several sources 

(𝑤0) is explored in Figure 5. It shows that the solution settles asymptotically to the, 𝑧3, for 𝑤0 >

2.48. Further, the solution approaches a periodic attractor for 𝑤0 ≤ 2.48. 
 

 
Figure 5. Dynamics of system (1) (a) time series with w0 =3; (b) phase portrait of (a); (c) time series with  w0 = 2.48; 

(d) phase portrait of (c). 

 

6. Conclusion 

This study modified the dissolved oxygen-plankton model by considering that zooplankton feeds 

only on the available phytoplankton. The objective is to determine how this type of interaction 



IHJPAS. 2024, 37(4) 

379 

 

impacts the dynamics of an aquatic ecosystem. The system was analyzed theoretically and 

numerically. The results of the theoretical analysis revealed three stable states. Depending on the 

conditions, the behaviour of the three constant states was either stable or unstable. The conditions 

necessary for a Hopf bifurcation around the positive stable state have been identified. 

Nonetheless, the numerical simulation deduced that system (1) always sways about the positive 

steady state when the stability criteria are met. Further, for small changing in some parameters, 

such that 𝑤0, 𝑠 and 𝛾, the system (1) shows limit cycle behaviour. For future work, we suggest 

considering climate change’s impact on the ocean’s oxygen-plankton dynamics. 

 

Acknowledgment 

         The authors are very grateful to the executive manager and editorial board members of the 

Ibn AL-Haitham Journal of Pure and Applied Sciences. 

 

Conflict of Interest  

         The authors declare that they have no conflicts of interest. 

 

Funding 

        There is no funding for the article. 
 

References 

1. Hull, V.; Parrella, L.; Falcucci, M. Modelling dissolved oxygen dynamics in coastal lagoons. 

Ecological Modelling, 2008, 211(3–4), 468–480. https://doi.org/10.1016/j.ecolmodel.2007.09.023 

2. Misra, A. K. Modeling the depletion of dissolved oxygen in a lake due to submerged macrophytes. 

Nonlinear Analysis: Modelling and Control, 2010, 15(2), 185–198. 

https://doi.org/10.15388/NA.2010.15.2.14353 

3. Misra, A. K.; Chandra, P.; Raghavendra, V. Modeling the depletion of dissolved oxygen in a lake due 

to algal bloom: Effect of time delay. Advances in Water Resources, 2011, 34(10), 1232–1238. 

https://doi.org/10.1016/j.advwatres.2011.05.010 

4. Gökçe, A. A mathematical study for chaotic dynamics of dissolved oxygen-phytoplankton interactions 

under environmental driving factors and time lag. Chaos, Solitons & Fractals, 2021, 151, 111268. 

https://doi.org/10.1016/j.chaos.2021.111268 

5. Sekerci, Y.; Petrovskii, S. Mathematical modelling of plankton–oxygen dynamics under the climate 

change. Bulletin of Mathematical Biology, 2015, 77(12), 2325–2353. https://doi.org/10.1007/s11538-

015-0126-0 

6. Hancke, K.; Glud, R. N. Temperature effects on respiration and photosynthesis in three diatom-

dominated benthic communities. Aquatic Microbial Ecology, 2004, 37(3), 265–281. 

7. Mandal, S.; Ray, S.; Ghosh, P. B. Modeling nutrient (dissolved inorganic nitrogen) and plankton 

dynamics at Sagar island of Hooghly–Matla estuarine system, West Bengal, India. Natural Resource 

Modeling, 2012, 25(4), 629–652. https://doi.org/10.1111/j.1939-7445.2011.00116.x 

8. Mondal, S.; Samanta, G.; De la Sen, M. Dynamics of Oxygen-Plankton Model with Variable 

Zooplankton Search Rate in Deterministic and Fluctuating Environments. Mathematics, 2022, 10(10), 

1641. https://doi.org/10.3390/math10101641 

9. Turner, J.; Vollrath, F.; Hesselberg, T. Wind speed affects prey-catching behaviour in an orb web 

spider. Naturwissenschaften, 2011, 98, 1063–1067. https://doi.org/10.1007/s00114-011-0854-4 

10. Das, A.; Samanta, G. P. Modeling the fear effect on a stochastic prey–predator system with additional 

food for the predator. Journal of Physics A: Mathematical and Theoretical, 2018, 51(46), 465601. 

https://doi.org/10.1088/1751-8121/aae4c6 

11. Arditi, R.; Ginzburg, L. R. Coupling in predator-prey dynamics: ratio-dependence. Journal of 

https://doi.org/10.1016/j.ecolmodel.2007.09.023
https://doi.org/10.15388/NA.2010.15.2.14353
https://doi.org/10.1016/j.advwatres.2011.05.010
https://doi.org/10.1016/j.chaos.2021.111268
https://doi.org/10.1007/s11538-015-0126-0
https://doi.org/10.1007/s11538-015-0126-0
https://doi.org/10.1111/j.1939-7445.2011.00116.x
https://doi.org/10.3390/math10101641
https://doi.org/10.1007/s00114-011-0854-4
https://doi.org/10.1088/1751-8121/aae4c6


IHJPAS. 2024, 37(4) 

380 

 

Theoretical Biology, 1989, 139(3), 311–326. https://doi.org/10.1016/S0022-5193(89)80211-5 

12. Sajan; Sasmal, S. K.; Dubey, B. A phytoplankton–zooplankton–fish model with chaos control: In the 

presence of fear effect and an additional food. Chaos: An Interdisciplinary Journal of Nonlinear 

Science, 2022, 32(1), 13114. https://doi.org/10.1063/5.0069474 

13. Gard, T. C.; Hallam, T. G. Persistence in food webs—I Lotka-Volterra food chains. Bulletin of 

Mathematical Biology, 1979, 41(6), 877–891. https://doi.org/10.1016/S0092-8240(79)80024-5 

14. Paine, R. T. Food web complexity and species diversity. The American Naturalist, 1966, 100(910), 65–

75. 

15. Meng, X.-Y.; Xiao, L. Stability and Bifurcation for a Delayed Diffusive Two-Zooplankton One-

Phytoplankton Model with Two Different Functions. Complexity, 2021, 2021, 5560157. 

https://doi.org/10.1155/2021/5560157 

16. Dhar, J.; Baghel, R. S. Role of dissolved oxygen on the plankton dynamics in spatio-temporal domain. 

Modeling Earth Systems and Environment, 2016, 2(1), 1–15. https://doi.org/10.1007/s40808-015-0061-

y 

17. Hirsch, M. W.; Smale, S.; Devaney, R. L. Differential equations, dynamical systems, and an 

introduction to chaos. 2012, Academic Press. 

18. Perko, L. Differential equations and dynamical systems (Vol. 7). 2013, Springer Science & Business 

Media. 

19. LaSalle, J. P. Stability theory and invariance principles. In Dynamical Systems, 1976, 211–222. 

Elsevier. https://doi.org/10.1016/B978-0-12-164901-2.50021-0 

20. Liu, Y.; Zhao, L.; Huang, X.; Deng, H. Stability and bifurcation analysis of two species amensalism 

model with Michaelis–Menten type harvesting and a cover for the first species. Advances in Difference 

Equations, 2018, 2018(1), 1–19. https://doi.org/10.1186/s13662-018-1752-2 

21. Kuznetsov, V. A.; Makalkin, I. A.; Taylor, M. A.; Perelson, A. S. Nonlinear dynamics of immunogenic 

tumors: parameter estimation and global bifurcation analysis. Bulletin of Mathematical Biology, 1994, 

56(2), 295–321. https://doi.org/10.1016/S0092-8240(05)80260-5 

22. Collings, J. B. Bifurcation and stability analysis of a temperature-dependent mite predator-prey 

interaction model incorporating a prey refuge. Bulletin of Mathematical Biology, 1995, 57(1), 63–76. 

https://doi.org/10.1007/BF02458316 

23. Yu, X.; Zhu, Z.; Li, Z. Stability and bifurcation analysis of two-species competitive model with 

Michaelis–Menten type harvesting in the first species. Advances in Difference Equations, 2020, 

2020(1), 1–25. https://doi.org/10.1186/s13662-020-02817-4 

24. Ali, N. Stability and bifurcation of a prey predator model with Qiwu’s growth rate for prey. 

International Journal of Mathematics and Computation, 2016, 27(2), 30–39. 

25. Jawad, S. R.; Al Nuaimi, M. Persistence and bifurcation analysis among four species interactions with 

the influence of competition, predation and harvesting. Iraqi Journal of Science, 2023, 64(3), 1369–

1390. https://doi.org/10.24996/ijs.2023.64.3.30 

26. Mukherjee, D. Study of fear mechanism in predator-prey system in the presence of competitor for the 

prey. Ecology and Genetics and Genomics, 2020, 15, 100052. 

https://doi.org/10.1016/j.egg.2020.100052 

https://doi.org/10.1016/S0022-5193(89)80211-5
https://doi.org/10.1063/5.0069474
https://doi.org/10.1016/S0092-8240(79)80024-5
https://doi.org/10.1155/2021/5560157
https://doi.org/10.1007/s40808-015-0061-y
https://doi.org/10.1007/s40808-015-0061-y
https://doi.org/10.1016/B978-0-12-164901-2.50021-0
https://doi.org/10.1186/s13662-018-1752-2
https://doi.org/10.1016/S0092-8240(05)80260-5
https://doi.org/10.1007/BF02458316
https://doi.org/10.1186/s13662-020-02817-4
https://doi.org/10.24996/ijs.2023.64.3.30
https://doi.org/10.1016/j.egg.2020.100052