Results in Control and Optimization 8 (2022) 100153

Contents lists available at ScienceDirect

Results in Control and Optimization

journal homepage: [Link]/locate/rico

Complex dynamics of a three species predator–prey model with two

nonlinearly competing species
Prabir Panja ∗, Sailen Gayen, Tridib Kar, Dipak Kumar Jana
Department of Applied Science, Haldia Institute of Technology, Purba Midnapore W.B., 721657, India

ARTICLE INFO ABSTRACT

Keywords: In this paper, a three species predator–prey model has been proposed incorporating nonlinear
Prey Lotka–Volterra competition between two competing species 𝑥 and 𝑦. It is assumed that the
Predator growth rate of 𝑥 increases due to its commensalism on 𝑧, but predator 𝑧 species consumes
Nonlinear competition
𝑦 species as prey following Holling type II functional response. Moreover the prey refuge
Commensalism
behavior has been considered on 𝑦 species. Positivity and boundedness of solutions are shown
Hopf bifurcation
analytically. Different equilibrium points are determined and the stability of the system has been
checked around these equilibrium points. Hopf bifurcation analysis has been done with respect
to the nonlinear relationship competition parameters 𝛾12 and 𝛾21 and other parameters 𝑝, 𝑣1 , 𝑣2
and 𝑣3 respectively. Also, the Transcritical bifurcation analysis of the system with respect to 𝛾12
and 𝛾21 have been done. From the analysis, it can be concluded that the nonlinear relationship
coefficients between two species 𝑥 and 𝑦 may control the instability of the proposed system.
Finally, some numerical simulation results have been presented to study the actual dynamics
of the proposed model.

1. Introduction

The study of ecological system using mathematical models becomes one of the most important research topics among applied
mathematicians, ecologists and biologists since the pioneering work of Lotka [1] and Volterra [2]. It is a well known fact that
mathematical models can be very useful to get new insights about the complex dynamics of ecological systems. Most of the
earlier classical ecological models of interacting populations typically have focused on two species (e.g., Rosenzweig-MacArthur [3],
May [4], Hassell [5]). After that three species ecological models have been proposed and their complex dynamics have been
reported [6–11]. Most of the authors considered only one of the following interactions at a time such as prey-predation, competition
or mutualism. A simple mathematical model among the interaction of phytoplankton, zooplankton and fish has been investigated
by Panja and Mondal [12]. Again, the impact of release of toxic substance by human in the aquatic environment on the dynamics of
phytoplankton, zooplankton and fish has been studied by Panja et al. [13]. This study suggests that due to the release of toxicants,
the harvesting of fish gradually decreases. Again, Panja et al. [14] studied the impacts of supply of additional food on three species
food chain model. This study suggests that the supply of additional food to the top predator may decrease the extinction of that
species and increase the harvesting rate. Also, the impacts of supplying additional food to predator with prey refuge in a predator–
prey model has been investigated by Ghosh et al. [15]. These are some research papers on predator–prey dynamics with linear
competition among the species.
The refuge is a behavior used by several species in ecological systems to decrease the predation risk. Khajanchi et al. [16]
investigated the impacts of constant prey refuge in a stage structured predator–prey model with ratio dependent functional response.

∗ Corresponding author.
E-mail address: prabirpanja@[Link] (P. Panja).

[Link]
Received 7 March 2022; Received in revised form 11 June 2022; Accepted 15 June 2022
Available online 20 June 2022
2666-7207/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license
([Link]
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Jana et al. [17] studied the effects of prey refuge and foraging ability of predator in a predator–prey system. The impacts of functional
response dependent predator–prey model has been investigated by Haque et al. [18]. There exist also many research articles where
prey refuge has been considered. Commensalism is a form of mutualism between individuals of two species in which one is benefited
and the other is unaffected. The commensal may obtain food, shelter or support from the host species. But there exist few research
articles on mutualism. Goh [19] studied the stability of some mathematical models on mutualism. After that Dhakne et al. [20]
investigated the stability of mutualistic interaction among three species. Kumar and Pattabhiramacharyulu [21] studied the impacts
of commensal in a three species predator–prey model. Gakkhar et al. [22] investigated the role of use of refuge as a protection for
prey species in a three species predator–prey model.
In most of the studies the linear relationship has been considered between two competing species. Taylor and Crizer [23]
introduced the concept of the nonlinear relationship between two competing species. In fact the nonlinear relationship between two
competing species is a better representation of reality than linear relationship. There exist some research articles [24–26] which
supports the consideration of the nonlinear competing relationship between species. In their model, they add nonlinear competition
terms to prevent the population of species 2 to have a smaller effect on the population of species 1 when the population density
of species 1 is very small compare to the population density of species 2 and vice versa. These facts motivated us to formulate and
analyze a mathematical model considering the nonlinear relationship between two competing species 𝑥 and 𝑦. In particular, we have
modified the model of Gakkhar and Gupta [27] by considering the nonlinear relationship between two competing species 𝑥 and 𝑦.
The paper is organized as follows: The mathematical model is formulated in Section 2. Positivity and boundedness of all solutions
of the proposed system are analyzed in Section 3. The existence of different equilibrium points and their stability have been analyzed
in Section 4. The Hopf bifurcation point is calculated in Section 5 with respect to the parameters 𝛾12 . We have presented numerical
simulation results in Section 6. Finally, in the last section, we give some of the main outcomes of our work.

2. Model formulation

In this section, we have considered a three species food chain model. Let 𝑋, 𝑌 and 𝑍 are the densities of two competing species
and the densities of predator of 𝑌 and a host of species 𝑋 at time 𝑇 respectively. We impose the following assumptions to formulate
the differential equations which describe the model system.

Assumption 1. The species 𝑋 grows logistically in the absence of species 𝑌 and in the absence of commensalism. Also, the density of
species 𝑋 decreases due to the competition with 𝑌 species. Again, the species 𝑋 is benefiting from 𝑍 species due to commensalism.
Hence, from the assumption, we have
[ ]
𝑑𝑋 𝑋 𝛼 𝑌
= 𝑟1 𝑋 1 − − 12 + 𝛿13 𝑋𝑍
𝑑𝑇 𝐾1 𝐾1

Assumption 2. The species 𝑌 grows logistically in the absence of the species 𝑋 and 𝑍. The density of 𝑌 decreases due to the
competition with 𝑋 species. Also, the density of 𝑌 decrease due to the predation of 𝑌 by 𝑍 species through Holling type II functional
response. Also, 𝑌 species show refuge behavior for the predator species 𝑍 during predation. Then we have
[ ]
𝑑𝑌 𝑌 𝛼 𝑋 𝑎(1 − 𝑝)𝑌 𝑍
= 𝑟2 𝑌 1 − − 21 −
𝑑𝑇 𝐾2 𝐾2 𝑏 + (1 − 𝑝)𝑌

Assumption 3. The density of the species 𝑍 increase due to the predation of 𝑌 species following Holling type II functional form.
Also, the density of 𝑍 decreases due to the death of that species. Then, we have
[ ]
𝑑𝑍 𝑎𝑐(1 − 𝑝)𝑌
= 𝑍 −𝑒 +
𝑑𝑇 𝑏 + (1 − 𝑝)𝑌
Hence, combining the above three equations, we can write
𝑑𝑋
= 𝑟1 𝑋[1 − 𝑋
− 𝐾
𝛼12 𝑌
] + 𝛿13 𝑋𝑍 ⎫
𝑑𝑇 𝐾1 1 ⎪
𝛼21 𝑋 ⎪
𝑑𝑌
𝑑𝑇
= 𝑟2 𝑌 [1 − 𝑌
𝐾2
− 𝐾 ] − 𝑎(1−𝑝)𝑌 𝑍
𝑏+(1−𝑝)𝑌 ⎬ (1)
2
𝑑𝑍 𝑎𝑐(1−𝑝)𝑌
⎪
= 𝑍[−𝑒 + ] ⎪
𝑑𝑇 𝑏+(1−𝑝)𝑌 ⎭
which satisfies the initial conditions 𝑋(0) ≥ 0, 𝑌 (0) ≥ 0 and 𝑍(0) ≥ 0. The biological meaning of the model parameters are given in
Table 1.
The interaction among three species such as Fox, Deer and Tiger can be considered as a real life example of our proposed model.
It is observed that Fox and Deer consume grasses and different parts of plants and they are competing for the same food source. The
intra-species competition among Fox and Deer Species has been found. It is seen that Tiger hunting Deer for food but tiger may not
consume the body parts of Deer completely. It is observed that Fox are consuming the remaining body parts of dead Deer hunted
by Tiger. So, it can be said that Fox is benefited from the Tiger species.
Model (2.1) does not investigate the effects on the ecological system when the population sizes of 𝑋 and 𝑌 are very small. But
different literature survey [23–26] shows that sizes of 𝑋 and 𝑌 population has great effects in their interactions. The model (1)
𝛼12 𝑌 𝛼21 𝑋
assumes that 𝑑𝑋𝑑𝑇
is negatively affected by the linear function − 𝐾 and 𝑑𝑌
𝑑𝑇
is negatively affected by the linear function − 𝐾 , that is,
1 2

2
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Table 1
Descriptions of the parametric values.
Name Biological meaning
𝑋, 𝑌 Densities of two logistically growing competing species.
𝑍 Density of a species which is a predator of 𝑌 .
𝑟1 , 𝑟2 Intrinsic growth rate of species 𝑋 and 𝑌 respectively.
𝐾1 , 𝐾2 Carrying capacity of species 𝑋 and 𝑌 respectively.
𝛼12 , 𝛼21 Inter-species competition between two species 𝑋 and 𝑌 respectively.
𝛿13 Commensal coefficient of 𝑋 over 𝑍 species.
𝑎 Attack rate of 𝑌 species by 𝑍 species.
𝑝 Refuge rate of 𝑌 species.
𝑐 Conservation rate of 𝑌 species.
𝑏 Half saturation constant for Holling type II function.
𝑒 Death rate of 𝑍 species.

population 𝑌 interferes with population 𝑋 in a linear fashion and vice versa. Biologically, it makes sense that this relationship could
be non-linear, and we suppose that as the population size of 𝑋 grows the population is more efficient at gathering the resource that
is shared with population 𝑌 thus resulting in population 𝑋 having a larger effect on population 𝑌 . On the other hand, if population
𝑋 is very small it may be very inefficient at gathering the shared resource and hence population 𝑋 has a smaller effect on population
𝑌 , when 0 < 𝑋 < 1, 𝑋 2 < 𝑋 and 𝑋 > 1, 𝑋 2 > 𝑋. The equation below is a modification of the model (1) which suggests a non-linear
effect of population 𝑌 on population 𝑋 and vice versa. So, we have introduced the concept of Nonlinear-relationship in the model
system (1), then it converted in the following form
𝑑𝑋
= 𝑟1 𝑋[1 − 𝑋 𝛼
− 𝐾12 𝑌 2 ] + 𝛿13 𝑋𝑍 ⎫
𝑑𝑇 𝐾1 1 ⎪
𝛼 ⎪
𝑑𝑌
𝑑𝑇
= 𝑟2 𝑌 [1 − 𝑌
𝐾2
− 𝐾21 𝑋 2 ] − 𝑎(1−𝑝)𝑌
𝑏+(1−𝑝)𝑌
𝑍
⎬ (2)
2
𝑑𝑍 𝑎𝑐(1−𝑝)𝑌
⎪
= 𝑍[−𝑒 + ] ⎪
𝑑𝑇 𝑏+(1−𝑝)𝑌 ⎭
The modified model (2) is more realistic than the model (1) because in the modified model, if population 𝑋 is very small (less than
1) then it has a smaller effect on population 𝑌 and if population 𝑋 is greater than 1 then it has a larger effect on population 𝑌 and
vice versa, where 1 is regarded as some ‘‘threshold’’ in a certain biology application.
Now to non-dimensionless the proposed mathematical model, we have introduced some variables and constants as 𝑡 = 𝑟1 𝑇 , 𝑥 =
𝑋 𝛼12 𝐾 2 𝛿 𝐾 𝑟 𝛼21 𝐾12
𝐾1
,𝑦 = 𝐾𝑌 , 𝑧 = 𝑟𝑎𝑍 , 𝛾12 = 𝐾 2 , 𝛾 = 13𝑎 2 , 𝑟 = 𝑟2 , 𝛾21 = , 𝑣1 = 𝑏
, 𝑣2 = 𝑒
and 𝑣3 = 𝑎𝑐
.
2 1 𝐾2 1 1 𝐾2 𝐾2 𝑟1 𝑟1
Then model (2.2) reduced to the following form
𝑑𝑥
𝑑𝑡
= 𝑥[1 − 𝑥 − 𝛾12 𝑦2 ] + 𝛾𝑥𝑧 ⎫
𝑑𝑦
= 𝑟𝑦[1 − 𝑦 − 𝛾21 𝑥2 ] − 𝑣 (1−𝑝)𝑦𝑧 ⎪
𝑑𝑡 1 +(1−𝑝)𝑦 ⎬
(3)
𝑑𝑧
=
𝑣 (1−𝑝)𝑦
𝑧[−𝑣2 + 𝑣 3+(1−𝑝)𝑦 ] ⎪
𝑑𝑡 1
⎭
with initial conditions 𝑥(0) ≥ 0, 𝑦(0) ≥ 0 and 𝑧(0) ≥ 0.

3. Positivity and boundedness of solutions

In this section, positivity and boundedness of solutions of the proposed system (3) have been investigated.

Theorem 1. All solutions of system (3) are non-negative.

Proof. Since 𝑥(0) ≥ 0, 𝑦(0) ≥ 0 and 𝑧(0) ≥ 0 then from system (3), we have
( 𝑡[ ] )
𝑥(𝑡) = 𝑥(0) exp 1 − 𝑥(𝑠) − 𝛾12 𝑦2 (𝑠) + 𝛾𝑧(𝑠) 𝑑𝑠 ≥ 0.
∫0
( 𝑡[ ] )
(1 − 𝑝)𝑧(𝑠)
𝑦(𝑡) = 𝑦(0) exp 𝑟(1 − 𝑦(𝑠) − 𝛾21 𝑥2 (𝑠)) − 𝑑𝑠 ≥ 0.
∫0 𝑣1 + (1 − 𝑝)𝑦(𝑠)
( 𝑡[ ] )
𝑣 (1 − 𝑝)𝑦(𝑠)
𝑧(𝑡) = 𝑧(0) exp −𝑣2 + 3 𝑑𝑠 ≥ 0.
∫0 𝑣1 + (1 − 𝑝)𝑦(𝑠)
Hence all solutions of system (3) are non-negative. □

Theorem 2. All solutions of system (3) are bounded.

Proof. From second equation of system (3), we have

𝑑𝑦 (1 − 𝑝)𝑦𝑧
= 𝑟𝑦[1 − 𝑦 − 𝛾21 𝑥2 ] −
𝑑𝑡 𝑣1 + (1 − 𝑝)𝑦
𝑑𝑦
≤ 𝑟𝑦(1 − 𝑦).
𝑑𝑡

3
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Then solving above equation, it is obtained that

𝑐1
𝑦(𝑡) ≤ , where 𝑐1 is a constant.
𝑐1 + 𝑒−𝑟𝑡
Then using the concept of differential inequality [28], as 𝑡 → ∞ we have

0 ≤ 𝑦(𝑡) ≤ 1.

Now, we define a function

1
𝑋1 = 𝑦 + 𝑧. (4)
𝑣3
Then differentiating the above Eq. (4) with respect to 𝑡, we have
𝑑𝑋1 𝑑𝑦 1 𝑑𝑧
= + .
𝑑𝑡 𝑑𝑡 𝑣3 𝑑𝑡
(1 − 𝑝)𝑦𝑧 (1 − 𝑝)𝑦𝑧 𝑣
= 𝑟𝑦[1 − 𝑦 − 𝛾21 𝑥2 ] − + − 2𝑧
𝑣1 + (1 − 𝑝)𝑦 𝑣1 + (1 − 𝑝)𝑦 𝑣3
𝑣2
≤ 𝑟𝑦(1 − 𝑦) − 𝑧.
𝑣3
𝑑𝑋1
+ 𝑣2 𝑋1 ≤ 𝑟𝑦(1 − 𝑦) + 𝑣2 𝑦. (5)
𝑑𝑡
𝑟+𝑣2 (𝑟+𝑣2 )2
Let us consider 𝑓1 (𝑦) = 𝑟𝑦(1 − 𝑦) + 𝑣2 𝑦, then the maximum value of 𝑓1 (𝑦) at 𝑦 = 2𝑟
is 4𝑟
. Then from Eq. (5), we have

𝑑𝑋1 (𝑟 + 𝑣2 )2
+ 𝑣2 𝑋 1 ≤ .
𝑑𝑡 4𝑟
Then solving above equation, we have
(𝑟 + 𝑣2 )2
𝑋1 ≤ + 𝑐2 𝑒−𝑣2 𝑡 , where 𝑐2 is a constant. (6)
4𝑟𝑣2
Now using differential inequality [28] and 𝑡 → ∞ in (6), we have
(𝑟 + 𝑣2 )2
0 < 𝑋1 ≤ .
4𝑟𝑣2
Next we define another function 𝑋2 as follows
1
𝑋2 = 𝑥 + 𝑦 + 𝑧.
𝑣3
Then differentiating above equation with respect to 𝑡, we have
𝑑𝑋2 𝑑𝑥 𝑑𝑦 1 𝑑𝑧
= + + .
𝑑𝑡 𝑑𝑡 𝑑𝑡 𝑣3 𝑑𝑡
𝑑𝑋2 (1 − 𝑝)𝑦𝑧 (1 − 𝑝)𝑦𝑧
i.e., = 𝑥[1 − 𝑥 − 𝛾12 𝑦2 ] + 𝛾𝑥𝑧 + 𝑟𝑦[1 − 𝑦 − 𝛾21 𝑥2 ] − +
𝑑𝑡 𝑣1 + (1 − 𝑝)𝑦 𝑣1 + (1 − 𝑝)𝑦
𝑣
− 2 𝑧.
𝑣3
𝑑𝑋2 𝑣
i.e., ≤ 𝑥 − 𝑥2 + 𝛾𝑥𝑧 + 𝑟𝑦 − 𝑟𝑦2 − 2 𝑧.
𝑑𝑡 𝑣3
𝑑𝑋2
i.e., + 𝑣2 𝑋2 ≤ 𝑥 + 𝑣2 𝑥 − 𝑥 + 𝛾𝑥𝑧 + 𝑟𝑦 − 𝑟𝑦2 + 𝑣2 𝑦.
2
𝑑𝑡
𝑑𝑋2 (𝑟 + 𝑣2 )2 (𝑟 + 𝑣2 )2
i.e., + 𝑣2 𝑋2 ≤ 𝑥 + 𝑣2 𝑥 − 𝑥2 + 𝛾𝑥𝑧 + , since max{𝑟𝑦 − 𝑟𝑦2 + 𝑣2 𝑦} =
𝑑𝑡 4𝑟 4𝑟
𝑑𝑋2 𝛾𝑣 (𝑟 + 𝑣 )2 (𝑟 + 𝑣 )2
3 2 2
i.e., + 𝑣2 𝑋 2 ≤ 𝑥 + 𝑣2 𝑥 − 𝑥 2 + 𝑥 + . (7)
𝑑𝑡 4𝑟𝑣2 4𝑟
𝛾𝑣3 (𝑟+𝑣2 )2 𝑎 𝑎 2 𝑎𝛾𝑣3 (𝑟+𝑣2 )2
Let us assume that 𝑓2 (𝑥) = 𝑥 + 𝑣2 𝑥 − 𝑥2 + 𝑥 4𝑟𝑣2
, then the maximum value of 𝑓2 (𝑥) at 𝑥 = 2
is 𝑏 where 𝑏 = 2
+ 𝑣2 𝑎2 − 𝑎4 + 8𝑟𝑣2
𝛾𝑣3 (𝑟+𝑣3 )2
and 𝑎 = 1 + 𝑣2 + 4𝑟𝑣2
. Then from Eq. (7), we have

𝑑𝑋2 (𝑟 + 𝑣2 )2
+ 𝑣2 𝑋 2 ≤ 𝑏 + .
𝑑𝑡 4𝑟
Then solving above equation, we have

1 (𝑟 + 𝑣2 )2
𝑋2 ≤ (𝑏 + ) + 𝑐3 𝑒−𝑣2 𝑡 .
𝑣2 4𝑟

4
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Using the concept of differential inequality [28] as 𝑡 → ∞, we have

1 (𝑟 + 𝑣2 )2
0 < 𝑋2 ≤ (𝑏 + ).
𝑣2 4𝑟
Hence all the solutions of system (3) are bounded. □

4. Equilibrium points and stability analysis

In this section, the possible equilibrium points are determined and then the stability of the system around these equilibrium
points has been discussed.

4.1. Equilibrium points

The system (3) have following possible equilibrium points

(i) The trivial equilibrium point 𝐸0 = (0, 0, 0).

(ii) The axial equilibrium point 𝐸1 = (1, 0, 0).
(iii) Another axial equilibrium point 𝐸2 = (0, 1, 0).
(iv) The boundary equilibrium point 𝐸3 = (𝑥, ̄ 0) where 𝑦̄ = 1 − 𝛾21 𝑥̄ 2 and 𝑥̄ is a positive root of the equation 𝛾12 𝛾21
̄ 𝑦, 2 𝑥
̄ 4 − 2𝛾12 𝛾21 𝑥̄ 2 +
𝑥̄ + 𝛾12 − 1 = 0. If 𝛾12 < 1, the biquadratic equation have either three or one positive root as per Descartes’ rule of sign. Thus
at least one positive root 𝑥∗ of the equation exists if 𝛾12 < 1.
𝑣1 𝑣2 ̂ 1 +(1−𝑝)𝑦)
𝑟(1−𝑦)(𝑣 ̂
(v) Another boundary equilibrium point 𝐸4 = (0, 𝑦, ̂ where 𝑦̂ = (1−𝑝)(𝑣
̂ 𝑧) and 𝑧̂ = , 𝐸4 exists if 𝑣3 > 𝑣2 and 𝑦̂ < 1.
3 −𝑣2 ) (1−𝑝) [
𝑣 𝑣 1
(vi) The interior equilibrium point 𝐸 ∗ = (𝑥∗ , 𝑦∗ , 𝑧∗ ) where 𝑦∗ = (1−𝑝)(𝑣 1 2
, exists if 𝑣 3 > 𝑣 2 , 𝑧∗ =
𝛾(1−𝑝)2 (𝑣3 −𝑣2 )2
(𝑥∗ − 1){𝑣23 + 𝑣22 −
3 −𝑣2 )
] 2 2
𝛾12 𝑣 𝑣 (1−𝑝)(𝑣 −𝑣 )
2𝑣2 𝑣3 } + 𝛾12 𝑣21 𝑣22 , exists if 𝑥∗ > 1 − 2 2 1 2 and 𝑥∗ is the root of the equation 𝑥∗ 2 + 𝐴𝑥∗ + 𝐵 = 0 where 𝐴 = 𝑟𝛾𝑣 3𝑣 2 and
𝑣2 +𝑣3 −2𝑣2 𝑣3
√1 3
𝐵 = (1−𝑝)𝑟𝛾𝑣 1𝑣 (𝑣 −𝑣 ) [(𝑟𝛾 + 𝛾12 )𝑣21 𝑣22 + 𝑣1 𝑣2 (𝑣3 − 𝑣2 )[𝑟𝛾𝑣1 − 𝑟𝛾(1 − 𝑝)] − (1 − 𝑝)(𝑣3 − 𝑣2 )2 [(1 − 𝑝) + 𝑟𝛾𝑣1 ]]. Then 𝑥∗ = 21 (−𝐴 ± 𝐴2 − 4𝐵)
1 3 3 2
and 𝑥∗ admits at least one positive root under the conditions (1) 𝐴 < 0 and 𝐵 < 0, (2) 𝐴 < 0, 𝐵 > 0, 𝐴2 −4𝐵 > 0 (3) 𝐴 > 0, 𝐵 < 0.

4.2. Local stability

To discuss the stability of the proposed system, we first determine the Jacobian matrix as follows
⎛ 1 − 2𝑥 − 𝛾12 𝑦2 + 𝛾𝑧 −2𝑥𝑦𝛾12 𝛾𝑥 ⎞
⎜ −2𝑟𝛾21 𝑥𝑦
(1−𝑝)𝑣1 𝑧
𝑟 − 2𝑟𝑦 − 𝑟𝛾21 𝑥2 − (𝑣 +(1−𝑝)𝑦) − 𝑣 (1−𝑝)𝑦 ⎟
𝐽 =⎜ 1
2
1+(1−𝑝)𝑦 ⎟
⎜ 0
𝑣1 𝑣3 (1−𝑝)𝑧
−𝑣2 +
𝑣3 (1−𝑝)𝑦 ⎟
⎝ (𝑣1 +(1−𝑝)𝑦)2 𝑣1 +(1−𝑝)𝑦 ⎠

Theorem 3. The equilibrium point 𝐸0 = (0, 0, 0) is a saddle point.

Proof. The characteristic equation of the Jacobian matrix 𝐽 at the point 𝐸0 is given by

(𝜆 − 1)(𝜆 − 𝑟)(𝜆 + 𝑣2 ) = 0.

So, the roots of the above equation are 1, 𝑟 and −𝑣2 . Since all the parameters related to the proposed mathematical model are
positive, hence the trivial equilibrium point is a saddle point. □

Theorem 4. The axial equilibrium point 𝐸1 = (1, 0, 0) is locally asymptotically stable if 𝛾21 > 1.

Proof. The characteristic equation of the Jacobian matrix 𝐽 at the point 𝐸1 is given by

(𝜆 + 1)(𝜆 − 𝑟(1 − 𝛾21 ))(𝜆 + 𝑣2 ) = 0.

Then the roots of the above equation are −1, 𝑟(1 − 𝛾21 ) and −𝑣2 . So, the axial equilibrium point is locally asymptotically stable if
𝛾21 > 1. □

Note. The proposed system undergoes transcritical bifurcation if 𝛾21 = 1.

𝑣1 𝑣2
Theorem 5. Another axial equilibrium point 𝐸2 = (0, 1, 0) is locally asymptotically stable if 𝛾12 > 1 and 𝑝 > 1 − 𝑣3 −𝑣2
.

Proof. The characteristic equation of the Jacobian matrix 𝐽 at 𝐸2 is given by

𝑣3 (1 − 𝑝)
(𝜆 − 1 + 𝛾12 )(𝜆 + 𝑟)(𝜆 + 𝑣2 − ) = 0.
𝑣1 + (1 − 𝑝)

5
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

𝑣 (1−𝑝)
The eigenvalues of the jacobian matrix 𝐽 at 𝐸2 are 1−𝛾12 , −𝑟 and −𝑣2 + 𝑣 3+(1−𝑝) . So, the equilibrium point 𝐸2 is locally asymptotically
1
stable if 𝛾12 > 1 and
𝑣3 (1 − 𝑝)
−𝑣2 + < 0.
𝑣1 + (1 − 𝑝)
i.e., (1 − 𝑝)(𝑣3 − 𝑣2 ) < 𝑣1 𝑣2 .
𝑣 𝑣
i.e., 𝑝 > 1 − 1 2 .
𝑣3 − 𝑣2
Hence the theorem. □
𝑣3 (1−𝑝) 𝑣1 𝑣2 +(1−𝑝)𝑣2
Note. The proposed system undergoes transcritical bifurcation about 𝐸2 = (0, 1, 0) if 𝛾12 = 1, 𝑣2 = 𝑣1 +(1−𝑝)
, 𝑣3 = 1−𝑝
and
𝑣1 𝑣2
𝑝=1− 𝑣3 −𝑣2
.

𝑣3 (1−𝑝)𝑦̄ ̄ 12 𝑦̄2
1−2𝑥−𝛾 1−2𝑦−2̄ 𝑥+4 ̄ 𝑥̄ 𝑦+2𝛾
̄ 12 𝑦̄
3
Theorem 6. Boundary equilibrium point 𝐸3 is locally asymptotically stable if 𝑣1 +(1−𝑝)𝑦̄
< 𝑣2 , 𝑟 > ̄ 12 𝑥̄ 2 −1
2𝑦+𝛾
and 𝛾21 < 𝑥̄ 2 −2𝑥̄ 3 +𝑦̄2 +3𝛾12 𝑥̄ 2 𝑦̄2
.

Proof. The characteristic equation of the Jacobian matrix at 𝐸3 = (𝑥, ̄ 0) is given by

̄ 𝑦,
( )
𝑣3 (1 − 𝑝)𝑦̄
− 𝑣2 − 𝜆 (𝜆2 + 𝐶1 𝜆 + 𝐶2 ) = 0.
𝑣1 + (1 − 𝑝)𝑦̄

where 𝐶1 = 2𝑥̄ + 𝛾12 𝑦̄2 + 2𝑟𝑦̄ + 𝑟𝛾21 𝑥̄ 2 − 𝑟 − 1 and 𝐶2 = 𝑟 − 2𝑟𝑦̄ − 𝑟𝛾21 𝑥̄ 2 − 2𝑟𝑥̄ + 4𝑟𝑥̄ 𝑦̄ + 2𝑟𝛾21 𝑥̄ 3 − 𝑟𝛾21 𝑦̄2 + 2𝑟𝛾12 𝑦̄3 − 3𝑟𝛾12 𝛾21 𝑥̄ 2 𝑦̄2 . One of the
𝑣 (1−𝑝)𝑦̄
root of the above characteristic equation is 𝑣 3+(1−𝑝)𝑦̄ − 𝑣2 and other two roots can be evaluated from the equation 𝜆2 + 𝐶1 𝜆 + 𝐶2 = 0.
1
According to Routh–Hurwitz criteria, the roots of the equation 𝜆2 + 𝐶1 𝜆 + 𝐶2 = 0 are negative real number or complex roots with
̄ 12 𝑦̄2
1−2𝑥−𝛾 1−2𝑦−2̄ 𝑥+4 ̄ 𝑥̄ 𝑦+2𝛾
̄ 𝑦̄3
negative real part if 𝐶1 > 0 and 𝐶2 > 0 i.e., 𝑟 > 2𝑦+𝛾
̄ 𝑥̄ 2 −1
and 𝛾21 < 𝑥̄ 2 −2𝑥̄ 3 +𝑦̄2 +3𝛾 𝑥̄122 𝑦̄2 . So, the boundary equilibrium point is locally
12 12
𝑣3 (1−𝑝)𝑦̄ ̄ 12 𝑦̄2
1−2𝑥−𝛾 1−2𝑦−2̄ 𝑥+4 ̄ 𝑥̄ 𝑦+2𝛾
̄ 12 𝑦̄
3
asymptotically stable if 𝑣1 +(1−𝑝)𝑦̄
< 𝑣2 , 𝑟 > ̄ 12 𝑥̄ 2 −1
2𝑦+𝛾
and 𝛾21 < 𝑥̄ 2 −2𝑥̄ 3 +𝑦̄2 +3𝛾12 𝑥̄ 2 𝑦̄2
. □

(1−𝑝)𝑣1 𝑧̂ 𝑣3 (1−𝑝)𝑦̂
Theorem 7. Another boundary equilibrium point 𝐸4 is locally asymptotically stable if 1+𝛾 𝑧̂ < 𝛾12 𝑦̂2 , 𝑟 < 2𝑟𝑦+
̂ (𝑣 2 and 𝑣1 +(1−𝑝)𝑦̂
< 𝑣2 .
1 +(1−𝑝)𝑦)
̂

Proof. The characteristic equation of the Jacobian matrix at 𝐸4 = (0, 𝑦, ̂ is given by

̂ 𝑧)
(1 − 𝑝)𝑣1 𝑧̂ 𝑣3 (1 − 𝑝)𝑦̂
(𝜆 − 1 + 𝛾12 𝑦̂2 − 𝛾 𝑧)(𝜆
̂ − 𝑟 + 2𝑟𝑦̂ + )(𝜆 − + 𝑣2 ) = 0.
̂2
(𝑣1 + (1 − 𝑝)𝑦) 𝑣1 + (1 − 𝑝)𝑦̂
(1−𝑝)𝑣1 𝑧̂ 𝑣3 (1−𝑝)𝑦̂
The roots of the above equation are 1 − 𝛾12 𝑦̂2 + 𝛾 𝑧,
̂ 𝑟 − 2𝑟𝑦̂ − (𝑣 2 and 𝑣1 +(1−𝑝)𝑦̂
− 𝑣2 . According to stability theory, the equilibrium
1 +(1−𝑝)𝑦)
̂
(1−𝑝)𝑣1 𝑧̂ 𝑣3 (1−𝑝)𝑦̂
point 𝐸4 is locally asymptotically stable if 1 + 𝛾 𝑧̂ < 𝛾12 𝑦̂2 , 𝑟 < 2𝑟𝑦̂ + ̂2
(𝑣1 +(1−𝑝)𝑦)
and 𝑣1 +(1−𝑝)𝑦̂
< 𝑣2 . □

Theorem 8. The interior equilibrium point 𝐸 ∗ is conditionally locally asymptotically stable.

Proof. The characteristic equation of the Jacobian matrix 𝐽 at 𝐸 ∗ is given by

𝜆3 + 𝜙1 𝜆2 + 𝜙2 𝜆 + 𝜙3 = 0. (8)
(1−𝑝)𝑣1 𝑧∗
where 𝑀11 = 1 − 2𝑥∗− 𝛾12 𝑦∗ 2
+ 𝛾𝑧∗ , 𝑀12 = −2𝛾12 𝑥∗ 𝑦∗ , 𝑀13 = 𝛾𝑥∗ , 𝑀21 = 2𝑟𝛾21 𝑥∗ 𝑦∗ , 𝑀22 = 𝑟 − 2𝑟𝑦∗ − 𝑟𝛾21 𝑥∗ 2 − (𝑣 +(1−𝑝)𝑦 ∗ )2 , 𝑀23 =
∗ ∗ 1
(1−𝑝)𝑦 ∗ 𝑣1 𝑣3 (1−𝑝)𝑧 𝑣3 (1−𝑝)𝑦
𝑣1 +(1−𝑝)𝑦∗
, 𝑀32 = (𝑣 +(1−𝑝)𝑦∗ )2 and 𝑀33 = 𝑣 +(1−𝑝)𝑦∗ −𝑣2 , 𝜙1 = −𝑀11 −𝑀22 −𝑀33 , 𝜙2 = 𝑀11 𝑀22 +𝑀11 𝑀33 +𝑀22 𝑀33 +𝑀23 𝑀32 +𝑀12 𝑀21
1 1
and 𝜙3 = 𝑀13 𝑀21 𝑀32 − 𝑀12 𝑀21 𝑀33 − 𝑀11 𝑀23 𝑀32 − 𝑀11 𝑀22 𝑀33 . According to Routh–Hurwitz criteria, the above characteristic
equation have negative real root or pair of imaginary roots with negative real part if 𝜙1 > 0, 𝜙3 > 0 and 𝜙1 𝜙2 > 𝜙3 holds. Now,
(1 − 𝑝)𝑣1 𝑧∗
𝜙1 = 2𝑥∗ + 𝛾12 𝑦∗ 2 + 2𝑟𝑦∗ + 𝑟𝛾21 𝑥∗ 2 + − 1 − 𝛾𝑧∗ − 𝑟.
(𝑣1 + (1 − 𝑝)𝑦∗ )2
2𝑟𝛾𝛾21 𝑣1 𝑣3 (1 − 𝑝)𝑥∗ 2 𝑦∗ 𝑧∗ 4𝑟𝛾12 𝛾21 𝑣3 (1 − 𝑝)𝑥∗ 2 𝑦∗ 3
𝜙3 = + − 4𝑟𝛾12 𝛾21 𝑣2 𝑥∗ 2 𝑦∗ 2 .
(𝑣1 + (1 − 𝑝)𝑦∗ )2 𝑣1 + (1 − 𝑝)𝑦∗
𝑟𝑣1 𝑣3 (1 − 𝑝)2 𝑦∗ 2 𝑧∗ (1 − 𝑝)4 𝑥∗ 𝑦∗ 2 𝑧∗ 2
𝜙1 𝜙2 − 𝜙3 = 𝑟𝑥∗ 2 𝑦∗ + 𝑟2 𝑥∗ 𝑦∗ 2 + +
(𝑣1 + (1 − 𝑝)𝑦∗ )3 (𝑣1 + (1 − 𝑝)𝑦∗ )4
(1 − 𝑝)𝑥∗ 𝑦∗ 𝑧∗
+ {4𝑟(1 − 𝑝)𝛾12 𝛾21 𝑥∗ 𝑦∗ − 3𝑟(1 − 𝑝)𝑦∗ − (1 − 𝑝 + 2𝑟𝛾𝛾21 𝑣1 𝑣3 )𝑥∗ }
(𝑣1 + (1 − 𝑝)𝑦∗ )2
𝑣1 𝑣3 (1 − 𝑝)4 𝑦∗ 2 𝑧∗ 2
− 4𝑟2 𝛾12 𝛾21 𝑥∗ 2 𝑦∗ 3 − 4𝑟𝛾21 𝛾12 𝑥∗ 3 𝑦∗ 2 − .
(𝑣1 + (1 − 𝑝)𝑦∗ )5
The interior equilibrium point 𝐸 ∗ is locally asymptotically stable if 𝜙1 > 0, 𝜙3 > 0 and 𝜙1 𝜙2 > 𝜙3 . □

6
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

5. Hopf bifurcation

For the variation of the model parameters, different dynamical behavior may occur in the mathematical model. The critical
parameter value at which qualitative change of dynamics occur is called bifurcation point [29]. The objective of this study is to
determine the stability of the system with the variation of different parameters of the system. We have considered here the parameter
𝛾12 associated to our mathematical model as a bifurcation parameter.

∗ are stated as follows:

Theorem 9. The necessary and sufficient conditions for the occurrence of Hopf bifurcation at 𝛾12 = 𝛾12
∗ ) > 0, 𝑖 = 1, 2, 3.
(𝑖) 𝜙𝑖 (𝛾12
∗ )𝜙 (𝛾 ∗ ) − 𝜙 (𝛾 ∗ ) = 0.
(𝑖𝑖) 𝜙1 (𝛾(12 2) 12 3 12
𝑑𝜆𝑖
(𝑖𝑖𝑖) 𝑅𝑒 𝑑𝛾12 ∗
≠ 0, 𝑖 = 1, 2, 3,
𝛾12 =𝛾12
where 𝜆𝑖 are the roots of the characteristic equation stated in Theorem 8.

Proof . The characteristic equation of system (8) around interior equilibrium 𝐸 ∗ is given by

𝜆3 + 𝜙1 𝜆2 + 𝜙2 𝜆 + 𝜙3 = 0. (9)
∗ then the above characteristic equation becomes
When the bifurcation parameter reaches its critical value i.e., when 𝛾12 = 𝛾12

(𝜆2 + 𝜙2 )(𝜆 + 𝜙1 ) = 0. (10)

since at the critical value of the parameter 𝛾12 = 𝛾12 ∗ , 𝜙 (𝛾 ∗ )𝜙 (𝛾 ∗ ) − 𝜙 (𝛾 ∗ ) = 0.

√ 1 12 2√12 3 12
So, the Eq. (9) have three roots such as 𝜆1 = 𝑖 𝜙2 , 𝜆2 = −𝑖 𝜙2 and 𝜆3 = −𝜙1 .
Let us assume that for 𝛾12 ∈ (𝛾12∗ − 𝜖, 𝛾 ∗ + 𝜖), the roots of the characteristic equation are
12

𝜆1 (𝛾12 ) = 𝜉1 (𝛾12 ) + 𝑖𝜉2 (𝛾12 ),

𝜆2 (𝛾12 ) = 𝜉1 (𝛾12 ) − 𝑖𝜉2 (𝛾12 ),
𝜆3 (𝛾12 ) = −𝜙1 .
( )
𝑑𝜆𝑖
Next, we want to verify the transversality condition 𝑅𝑒 𝑑𝛾12 ∗
≠ 0, 𝑖 = 1, 2, 3.
𝛾12 =𝛾12
Substituting 𝜆1 (𝛾12 ) = 𝜉1 (𝛾12 ) + 𝑖𝜉2 (𝛾12 ) in Eq. (8) and calculating the derivatives we have

𝑃 (𝛾12 )𝜉1′ (𝛾12 ) − 𝑄(𝛾12 )𝜉2 (𝛾12 ) + 𝑈 (𝛾12 ) = 0,

𝑄(𝛾12 )𝜉1′ (𝛾12 ) + 𝑃 (𝛾12 )𝜉2 (𝛾12 ) + 𝑉 (𝛾12 ) = 0,

where

𝑃 (𝛾12 ) = 3𝜉12 (𝛾12 ) + 2𝜙1 (𝛾12 )𝜉1 (𝛾12 ) + 𝜙2 (𝛾12 ) − 3𝜉22 (𝛾12 ),
𝑄(𝛾12 ) = 6𝜉1 (𝛾12 )𝜉2 (𝛾12 ) + 2𝜙1 (𝛾12 )𝜉2 (𝛾12 ),
𝑈 (𝛾12 ) = 𝜉12 (𝛾12 )𝜙′1 (𝛾12 ) + 𝜙′2 (𝛾12 )𝜉1 (𝛾12 ) + 𝜙′3 (𝛾12 ) − 𝜙′1 (𝛾12 )𝜉22 (𝛾12 ),
𝑉 (𝛾12 ) = 2𝜉1 (𝛾12 )𝜉2 (𝛾12 )𝜙′1 (𝛾12 ) + 𝜙′2 (𝛾12 )𝜉2 (𝛾12 ).
√ √
∗ ) = −2𝜙 (𝛾 ∗ ), 𝑄(𝛾 ∗ ) = 2𝜙 (𝛾 ∗ ) 𝜙 (𝛾 ∗ ), 𝑈 (𝛾 ∗ ) = 𝜙′ (𝛾 ∗ ) − 𝜙′ (𝛾 ∗ )𝜙 (𝛾 ∗ ), 𝑉 (𝛾 ∗ ) = 𝜙′ (𝛾 ∗ ) 𝜙 (𝛾 ∗ ). Therefore,
Now, 𝑃 (𝛾12 2 12 12 1 12 2 12 12 3 12 1 12 2 12 12 2 12 2 12

∗ ∗ ∗
𝑄(𝛾12 )𝑉 (𝛾12 ) + 𝑃 (𝛾12 )𝑈 (𝛾12 ) ∗
𝑑
(𝑅𝑒(𝜆𝑖 (𝛾12 )))|𝛾12 =𝛾 ∗ = − ∗ ) + 𝑄2 (𝛾 ∗ )
≠ 0, for 𝑖 = 1, 2
𝑑𝛾12 12 𝑃 2 (𝛾12 12
∗ )𝑉 (𝛾 ∗ ) + 𝑃 (𝛾 ∗ )𝑈 (𝛾 ∗ ) ≠ 0 and 𝜆 (𝛾 ∗ ) = −𝜙 ≠ 0.
if 𝑄(𝛾12 12 12 12 3 12 1
∗ . Hence the theorem.
Therefore, transversality conditions hold. This indicates that Hopf bifurcation occurs at 𝛾12 = 𝛾12

5.1. Stability and direction of Hopf bifurcation

In this section, we have determined the direction and stability criteria of Hopf bifurcating periodic solution by reducing the
system of differential Eq. (3) into normal form following the procedure stated by Hassard et al. [29]. For this reason, we have
introduced new variables 𝑥 = 𝑥∗ + 𝑤1 , 𝑦 = 𝑦∗ + 𝑤2 and 𝑧 = 𝑧∗ + 𝑤3 then the system (3) reduces to the form

𝑋̇ = 𝐴𝑋 + 𝐵. (11)

where 𝑋̇ denotes the time derivative of 𝑋. Also, 𝐴𝑋 and 𝐵 are the linear part and the nonlinear part of the system respectively.
Moreover,
⎡ 𝑤1 ⎤ ⎡ 𝑎11 𝑎12 𝑎13 ⎤ ⎡ 𝑏11 ⎤
𝑋 = ⎢ 𝑤2 ⎥,𝐴 = ⎢ 𝑎 𝑎22 𝑎23 ⎥,𝐵 = ⎢ 𝑏 ⎥.
⎢ ⎥ ⎢ 21 ⎥ ⎢ 12 ⎥
⎣ 𝑤3 ⎦ ⎣ 0 𝑎32 𝑎33 ⎦ ⎣ 𝑏13 ⎦

7
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

where the expression of 𝑎11 , 𝑎12 , 𝑎13 , 𝑎21 , 𝑎22 , 𝑎23 , 𝑎32 , 𝑎33 and 𝑏11 , 𝑏12 , 𝑏13 , 𝐻1 are given as follows:

𝑎11 = 1 − 2𝑥∗ − 𝛾12 𝑦∗ 2 + 𝛾𝑧∗ , 𝑎12 = −2𝑥∗ 𝑦∗ 𝛾12 , 𝑎13 = 𝛾𝑥∗ , 𝑎21 = −2𝑟𝛾21 𝑥∗ 𝑦∗ ,
(1 − 𝑝)𝑣1 𝑧∗ (1 − 𝑝)𝑦∗ 𝑣1 𝑣3 (1 − 𝑝)𝑧∗
𝑎22 = 𝑟 − 2𝑟𝑦∗ − 𝑟𝛾21 𝑥∗ 2 − , 𝑎23 = − ∗
, 𝑎32 = ,
(𝑣1 + (1 − 𝑝)𝑦∗ )2 𝑣1 + (1 − 𝑝)𝑦 (𝑣1 + (1 − 𝑝)𝑦∗ )2
𝑣 (1 − 𝑝)𝑦 ∗
1
𝑎33 = −𝑤2 + 3 , 𝐻1 = ,
𝑣1 + (1 − 𝑝)𝑦∗ 𝑣1 + (1 − 𝑝)𝑦∗
𝑏11 = −𝑤21 − 2𝑤1 𝑤2 𝛾12 𝑦∗ − 𝑤22 𝛾12 𝑥∗ − 𝑤1 𝑤22 𝛾12 + 𝑤1 𝑤3 𝛾,
𝑏12 = −𝑟𝑤22 − 𝑟𝛾21 𝑤21 𝑦∗ − 2𝑤1 𝑤2 𝑟𝛾21 𝑥∗ − 𝑤21 𝑤2 𝑟𝛾21 − 𝑤1 𝑤2 (1 − 𝑝)𝑧∗ 𝐻12 ,
− 𝑤1 𝑤2 𝑤3 (1 − 𝑝)𝐻12 + 𝑤1 𝑤22 (1 − 𝑝)2 𝑧∗ 𝐻13 − 𝑤1 𝑤52 (1 − 𝑝)3 𝑧∗ 𝐻14 + 𝑤1 𝑤22 𝑤3 (1 − 𝑝)2 𝐻13 ,
𝑏13 = 𝑤1 𝑤2 𝑤3 (1 − 𝑝)𝑧∗ 𝐻12 + 𝑤23 (1 − 𝑝)𝑦∗ 𝐻1 + 𝑤1 𝑤2 𝑤23 (1 − 𝑝)𝐻12 − 𝑤1 𝑤22 𝑤3 (1 − 𝑝)2 𝑧∗ 𝐻13 ,
+ 𝑤1 𝑤52 𝑤3 (1 − 𝑝)3 𝑧∗ 𝐻14 − 𝑤1 𝑤22 𝑤3 (1 − 𝑝)2 𝐻13 .
√
Again, from Theorem 9 we have the equation (𝜆2 + 𝜙2 )(𝜆 + 𝜙1 ) = 0. It is clear that 𝜆1,2 = ±𝑖𝜔0 where 𝜔0 = 𝜙2 and one eigenvalues
have negative real part −𝜙1 .
Now, we attempt to find a transformation matrix 𝑃 which reduces the matrix 𝐴 in the form

⎡ 0 −𝜔0 0 ⎤
𝑃 −1 𝐴𝑃 = ⎢ 𝜔0 0 0 ⎥.
⎢ ⎥
⎣ 0 0 −𝜙1 ⎦
and also the corresponding non-singular matrix 𝑃 is given by

⎡ 𝑃11 𝑃12 𝑃13 ⎤

𝑃 = ⎢ 𝑃21 𝑃22 𝑃23 ⎥.
⎢ ⎥
⎣ 𝑃31 𝑃32 𝑃33 ⎦

where 𝑃11 = 𝑎12 𝑎23 − 𝑎13 𝑎22 , 𝑃12 = 𝜔0 𝑎13 , 𝑃21 = 𝑎13 𝑎21 − 𝑎11 𝑎23 , 𝑃22 = 𝜔0 𝑎23 , 𝑃31 = 𝑎11 𝑎22 − 𝑎12 𝑎21 − 𝜔20 , 𝑃32 = −(𝜔0 𝑎11 + 𝜔0 𝑎22 ), 𝑃13 =
𝑎12 𝑎23 − 𝑎13 (𝑎22 + 𝜙1 ), 𝑃23 = 𝑎13 𝑎21 − 𝑎23 (𝑎11 + 𝜙1 ), 𝑃33 = (𝑎11 + 𝜙1 )(𝑎22 + 𝜙1 ).
To fulfill the normal form of the Eq. (11), we use another change of variable i.e., 𝑋 = 𝑃 𝑌 , where 𝑌 = (𝑦1 , 𝑦2 , 𝑦3 )𝑇 .
After some algebraic calculations the Eq. (11) takes the following form

𝑌̇ = 𝛴𝑌 + 𝐹 . (12)
𝐹 1 (𝑦1 , 𝑦2 , 𝑦3 )
⎡ ⎤ ⎡ 𝑓 1 (𝑦1 , 𝑦2 , 𝑦3 ) ⎤
where 𝛴 = 𝑃 −1 𝐴𝑃 and 𝐹 = 𝑃 −1 𝑓 = ⎢ ⎥ and 𝑓 is given by 𝑓 = ⎢
𝐹 2 (𝑦1 , 𝑦2 , 𝑦3 ) 𝑓 2 (𝑦1 , 𝑦2 , 𝑦3 ) ⎥ where the expression of
⎢ ⎥ ⎢ ⎥
⎣ ⎦
𝐹 3 (𝑦1 , 𝑦2 , 𝑦3 ) ⎣ 𝑓 3 (𝑦1 , 𝑦2 , 𝑦3 ) ⎦
1 2 3
𝑓 (𝑦1 , 𝑦2 , 𝑦3 ), 𝑓 (𝑦1 , 𝑦2 , 𝑦3 ) and 𝑓 (𝑦1 , 𝑦2 , 𝑦3 ) are given below.

𝑓 1 (𝑦1 , 𝑦2 , 𝑦3 ) = 𝛾(𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 ) − (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )2
− 2𝛾12 𝑦∗ (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 ) − 𝛾12 𝑥∗ (𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2
− 𝛾12 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2 .
𝑓 2 (𝑦1 , 𝑦2 , 𝑦3 ) = −𝑟(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2 − 𝑟𝛾21 𝑦∗ (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )2
− 2𝑟𝛾21 𝑥∗ (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )
− 𝑟𝛾21 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )2 (𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )
− (1 − 𝑝)𝑧∗ 𝐻12 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )2 (𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )
− (1 − 𝑝)𝐻12 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )(𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 )
+ (1 − 𝑝)2 𝑧∗ 𝐻13 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2
− (1 − 𝑝)3 𝑧∗ 𝐻14 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )5
+ (1 − 𝑝)2 𝐻13 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2 (𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 ).
𝑓 3 (𝑦1 , 𝑦2 , 𝑦3 ) = (1 − 𝑝)𝑧∗ 𝐻12 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )(𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 )
+ (1 − 𝑝)𝑦∗ 𝐻1 (𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 )2
+ (1 − 𝑝)𝐻12 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )(𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 )2
− (1 − 𝑝)2 𝑧∗ 𝐻13 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2 (𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 )
+ (1 − 𝑝)3 𝑧∗ 𝐻14 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )5 (𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 )
− (1 − 𝑝)2 𝐻13 (𝑃11 𝑦1 + 𝑃12 𝑦2 + 𝑃13 𝑦3 )(𝑃21 𝑦1 + 𝑃22 𝑦2 + 𝑃23 𝑦3 )2 (𝑃31 𝑦1 + 𝑃32 𝑦2 + 𝑃33 𝑦3 ).

Eq. (12) is the normal form of Eq. (11) from which the stability and direction of the Hopf bifurcation can be computed. First term
on the right hand side of Eq. (11) is linear and the second term is nonlinear in y’s. For evaluating the direction of periodic solution,

8
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

we can evaluated the following quantities at 𝛾12 = 𝛾12 ∗ .

[ 2 1 ( )]
1 𝜕 𝐹 𝜕2 𝐹 1 𝜕2 𝐹 2 𝜕2 𝐹 2
𝑔11 = + +𝑖 + .
4 𝜕𝑦2 𝜕𝑦 2 𝜕𝑦21 𝜕𝑦22
1 2
[ ( 2 2 )]
1 𝜕2 𝐹 1 𝜕2 𝐹 1 𝜕2 𝐹 2 𝜕 𝐹 𝜕2 𝐹 2 𝜕2 𝐹 1
𝑔02 = − −2 +𝑖 − +2 .
4 𝜕𝑦2 𝜕𝑦22 𝜕𝑦1 𝜕𝑦2 𝜕𝑦21 𝜕𝑦22 𝜕𝑦1 𝜕𝑦2
[ 2 11 ( )]
1 𝜕 𝐹 𝜕2 𝐹 1 𝜕2 𝐹 2 𝜕2 𝐹 2 𝜕2 𝐹 2 𝜕2 𝐹 1
𝑔20 = − +2 +𝑖 − −2 .
4 𝜕𝑦2 𝜕𝑦22 𝜕𝑦 𝜕𝑦
1 2 𝜕𝑦12 𝜕𝑦22 𝜕𝑦 1 𝜕𝑦2
1

[ ( 3 2 )]
1 𝜕3 𝐹 1 𝜕3 𝐹 1 𝜕3 𝐹 2 𝜕3 𝐹 2 𝜕 𝐹 𝜕3 𝐹 2 𝜕3 𝐹 1 𝜕3 𝐹 1
𝐺21 = + + + +𝑖 + − − .
8 𝜕𝑦 3 2
𝜕𝑦1 𝜕𝑦2 𝜕𝑦1 𝜕𝑦2 2 𝜕𝑦23 3
𝜕𝑦1 2 2
𝜕𝑦1 𝜕𝑦2 𝜕𝑦1 𝜕𝑦2 𝜕𝑦23
1
[ ( 2 2 )]
𝑗 1 𝜕2 𝐹 1 𝜕2 𝐹 2 𝜕 𝐹 𝜕2 𝐹 1
𝐺110 = + +𝑖 − .
2 𝜕𝑦1 𝜕𝑦𝑗 𝜕𝑦2 𝜕𝑦𝑗 𝜕𝑦1 𝜕𝑦𝑗 𝜕𝑦2 𝜕𝑦𝑗
[ 2 1 ( )]
𝑗 1 𝜕 𝐹 𝜕2 𝐹 2 𝜕2 𝐹 2 𝜕2 𝐹 1
𝐺101 = − +𝑖 + .
2 𝜕𝑦1 𝜕𝑦𝑗 𝜕𝑦2 𝜕𝑦𝑗 𝜕𝑦1 𝜕𝑦𝑗 𝜕𝑦2 𝜕𝑦𝑗
[ 2 𝑗 ] [ ]
1 𝜕 𝐹 𝜕2 𝐹 𝑗 1 𝜕2 𝐹 𝑗 𝜕2 𝐹 𝑗 𝜕2 𝐹 𝑗
ℎ𝑗11 = + , ℎ𝑗20 = − − 2𝑖 .
4 𝜕𝑦 2 𝜕𝑦22 4 𝜕𝑦 2 𝜕𝑦22 𝜕𝑦1 𝜕𝑦2
1 1
ℎ𝑗11 ℎ𝑗20
𝑤𝑗11 = , 𝑤𝑗20 = , 𝑗 = 1, 2.
𝑞𝑗 𝑞𝑗 + 2𝑖𝜔0
1
𝑔21 = 𝐺21 + 2(𝐺110 𝑤111 + 𝐺110
1
𝑤211 ) + (𝐺101
1
𝑤120 + 𝐺101
2
𝑤220 ).

Now, we can compute the following quantities:

( )
𝑖 |𝑔 |2 𝑔
𝐶1 (0) = 𝑔11 𝑔20 − 2|𝑔11 |2 − 02 + 21 .
2𝜔0 3 2
𝑅𝑒{𝐶1 (0)}
𝜇2 = − ∗ )} .
𝑅𝑒{𝜆′ (𝛾12
∗ )}
𝐼𝑚{𝐶1 (0)} + 𝜇2 𝐼𝑚{𝜆′ (𝛾12
𝜏2 = − .
𝜔0
𝛽2 = 2𝑅𝑒{𝐶1 (0)}.

Hence, the Hopf bifurcation of the system (3) at 𝐸 ∗ is non-degenerate and supercritical provided the sign of 𝜇2 is positive and
∗ (𝛾 < 𝛾 ∗ ); 𝛽 determines the stability of bifurcating
subcritical if 𝜇2 is negative. The bifurcating periodic solutions exist for 𝛾12 > 𝛾12 12 12 2
periodic solutions. The bifurcating periodic solutions are orbitally asymptotically stable (unstable) if 𝛽2 < 0(𝛽2 > 0) and 𝜏2 determines
the period of the bifurcating periodic solutions, the period increases (decreases) if 𝜏2 > 0(𝜏2 < 0).

6. Numerical simulation

We have solved the system (3), numerically using MATLAB to get a better insight of the proposed model. All the parametric
values have been considered from the research paper published by Gakkhar and Gupta [27]. We have taken the parametric values
for Fig. 1 are 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 0.4, 𝛾21 = 0.1, 𝑝 = 0.4, 𝑣1 = 0.14, 𝑣2 = 0.32009, 𝑣3 = 0.5. The Fig. 1(a) shows that the system (3)
is locally asymptotically stable around the interior equilibrium point (0.9376, 0.4152, 0.1611). Also, the eigenvalues of the Jacobian
matrix at the interior equilibrium point (0.9376, 0.4152, 0.1611) are −0.9508, −0.0177−0.1672𝑖 and −0.0177+0.1672𝑖 i.e., the conditions
of local stability are verified. We plot the phase portrait diagram of the system (3) in Fig. 1(b) by taking the same set of parametric
values used in Fig. 1(a). Hence it can be concluded that the interior equilibrium point (0.9376, 0.4152, 0.1611) is locally asymptotically
stable for this set of parameters.
Time evolution of the solution and bifurcation diagram of system (3) with respect to 𝛾12 have been drawn in Fig. 2 by taking the
parametric values 𝑟 = 0.5, 𝛾 = 0.04, 𝛾21 = 0.1, 𝑝 = 0.4, 𝑣1 = 0.14, 𝑣2 = 0.32009, 𝑣3 = 0.5. Fig. 2(a) shows existence of oscillatory behavior
of system (3) for this set of parametric values for 𝛾12 = 3.0. From Fig. 2(b–d), it is seen that when we change the interspecies
coefficient 𝛾12 of 𝑦 to 𝑥 from 1.0 to 5.0, then system (3) undergoes Hopf bifurcations. For this set of parametric values, we have
calculated the critical value of the interspecies competition coefficient 𝛾12 for Hopf bifurcation using Theorem 9 and obtain the
critical value of 𝛾12 as 1.64. Also, the system (3) shows the oscillatory behavior in (1.64 < 𝛾12 < 4.6) which is clear from Fig. 2(b–d).
The system (3) shows a stable steady state in (1 ≤ 𝛾12 < 1.64) and (4.6 < 𝛾12 ≤ 5).
Using the parametric values 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 1.5, 𝑝 = 0.4, 𝑣1 = 0.2001, 𝑣2 = 0.32009, 𝑣3 = 0.5, the solution of a system (3) with
respect to time and bifurcation diagram of system (3) with respect to 𝛾21 have been drawn in Fig. 3. The existence of oscillatory
behavior of the system (3) for this set of parametric values for 𝛾21 = 0.5 has been shown in Fig. 3(a). From Fig. 3(b–d), it is seen that
when we change the interspecies coefficient 𝛾21 of 𝑥 to 𝑦 from 0.3 to 1.2, then system (3) undergoes Hopf bifurcations. For this set
of parametric values, we have calculated the critical value of the interspecies competition coefficient 𝛾21 for Hopf bifurcation using
Theorem 9 and obtain the critical value of 𝛾21 as 0.408. Also, the system (3) shows the oscillatory behavior in (0.408 < 𝛾21 < 0.903)
which is clear from Fig. 2(b–d). The system (3) shows stable steady state for 𝑥 species in (0.3 ≤ 𝛾21 < 0.408) and (0.903 < 𝛾21 ≤ 1.2)

9
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Fig. 1. (a) Solving system (3) time evolution of 𝑥, 𝑦 and 𝑧 are plotted by taking 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 0.4, 𝛾21 = 0.1, 𝑝 = 0.4, 𝑣1 = 0.14, 𝑣2 = 0.32009, 𝑣3 = 0.5. (b)
Phase portrait diagram has been shown for the same parametric values used in (a) with initial population (0.7, 0.4, 0.2).

Fig. 2. (a) Solving system (3) time evolution of 𝑥, 𝑦 and 𝑧 are plotted by considering 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 3.0, 𝛾21 = 0.1, 𝑝 = 0.4, 𝑣1 = 0.14, 𝑣2 = 0.32009, 𝑣3 = 0.5.
Bifurcation diagram of system (3) with respect to 𝛾12 of (b) 𝑥, (c) 𝑦 and (d) 𝑧.

and also the extinction of 𝑥 species for (0.579 < 𝛾21 < 0.903) is observed in Fig. 3(b). Also, stable steady state and extinction of 𝑦
species have been seen in Fig. 3(c) within (0.3 ≤ 𝛾21 < 0.408 and 0.903 < 𝛾21 < 1.0) and (1.0 < 𝛾21 ≤ 1.2) respectively. Again, the stable
steady state behavior and the extinction of 𝑧 species have been observed in Fig. 3(d) within (0.3 ≤ 𝛾21 < 0.408) and (0.786 < 𝛾21 ≤ 1.2)
respectively.
For 𝑝 = 0.1 and the set of other parametric values 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 1.5, 𝛾21 = 0.3, 𝑣1 = 0.2, 𝑣2 = 0.35, 𝑣3 = 0.5, the interior
equilibrium point 𝐸 ∗ will be unstable and oscillating behavior is shown in Fig. 4(a). Again, for 𝑣1 = 0.1 and the set of other parametric
values 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 0.4, 𝛾21 = 0.1, 𝑝 = 0.4, 𝑣2 = 0.32009, 𝑣3 = 0.5, the existence of oscillatory behavior has been observed in
Fig. 4(b). Also, for 𝑣2 = 0.2 and the set of other parametric values 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 1.5, 𝛾21 = 0.3, 𝑝 = 0.4, 𝑣1 = 0.2, 𝑣3 = 0.5 the
interior equilibrium point 𝐸 ∗ will be unstable and oscillatory dynamics of the system (3) has been observed in Fig. 4(c). Again, the
unstable solution of a system (3) around interior equilibrium point 𝐸 ∗ has been shown in Fig. 4(d) for 𝑣3 = 1.0 and other parametric

10
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Fig. 3. (a) Solving system (3) time evolution of 𝑥, 𝑦 and 𝑧 are plotted by choosing 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 1.5, 𝛾21 = 0.5, 𝑝 = 0.4, 𝑣1 = 0.2001, 𝑣2 = 0.32009, 𝑣3 = 0.5.
Bifurcation diagram of system (3) with respect to 𝛾21 of (b) 𝑥, (c) 𝑦 and (d) 𝑧.

Fig. 4. Time evolution of solutions of system (3) are plotted by taking parametric values 𝑟 = 0.5, 𝛾 = 0.04 ∶ (a) for 𝛾12 = 1.5, 𝛾21 = 0.3, 𝑝 = 0.1, 𝑣1 = 0.2, 𝑣2 =
0.35, 𝑣3 = 0.5, (b) for 𝛾12 = 0.4, 𝛾21 = 0.1, 𝑝 = 0.4, 𝑣1 = 0.1, 𝑣2 = 0.32009, 𝑣3 = 0.5, (c) for 𝛾12 = 1.5, 𝛾21 = 0.3, 𝑝 = 0.4, 𝑣1 = 0.2, 𝑣2 = 0.2, 𝑣3 = 0.5 and (d) for
𝛾12 = 1.5, 𝛾21 = 0.3, 𝑝 = 0.4, 𝑣1 = 0.2, 𝑣2 = 0.35, 𝑣3 = 1.0.

11
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Fig. 5. Bifurcation diagrams of 𝑥 with respect to parameters 𝑝, 𝑣1 , 𝑣2 and 𝑣3 .

Fig. 6. Bifurcation diagrams of 𝑦 with respect to parameters 𝑝, 𝑣1 , 𝑣2 and 𝑣3 .

values 𝑟 = 0.5, 𝛾 = 0.04, 𝛾12 = 1.5, 𝛾21 = 0.3, 𝑝 = 0.4, 𝑣1 = 0.2, 𝑣2 = 0.35. Then the Hopf bifurcation diagram of the system (3) for 𝑥, 𝑦
and 𝑧 species with respect to 𝑝, 𝑣1 , 𝑣2 and 𝑣3 for the above mentioned set of parametric values have been presented in Figs. 5–7
respectively.
From Fig. 5(a–d), it is seen that when we change the prey refuge rate 𝑝, half saturation constant 𝑣1 , the death rate of predator
𝑣2 and conservation rate of prey 𝑣3 from 0.3 to 1.2, 0 to 0.4, 0 to 0.4 and 0 to 30 respectively then the 𝑥 species of system (3)

12
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

Fig. 7. Bifurcation diagrams of 𝑧 with respect to parameters 𝑝, 𝑣1 , 𝑣2 and 𝑣3 .

undergoes Hopf bifurcations. For the same set of parametric values used in Fig. 4(a–d), we have calculated the critical value of the
prey refuge rate 𝑝, half saturation constant 𝑣1 , the death rate of predator 𝑣2 and conservation rate of prey 𝑣3 for Hopf bifurcation
using Theorem 9 and obtain the critical value of 𝑝 ≅ 0.15, 𝑣1 ≅ 0.1287, 𝑣2 ≅ 0.304 and 𝑣3 ≅ 0.696. From Fig. 5(a), it is observed
that 𝑥 species show oscillatory behavior in (0 ≤ 𝑝 < 0.15). Also, 𝑥 species show stable steady state solution in (0.15 < 𝑝 < 0.47).
The extinction of 𝑥 species is observed in (0.47 < 𝑝 ≤ 1). Again, from Fig. 5(b), it is seen that 𝑥 species show oscillatory behavior
in (0 ≤ 𝑣1 < 0.1287) and also stable steady state behavior is observed in (0.1287 < 𝑣1 ≤ 0.4). Also, from Fig. 5(c) it is observed that
𝑥 species show oscillatory behavior in (0 ≤ 𝑣2 < 0.304). Also, 𝑥 species shows stable steady state behavior in (0.304 < 𝑣2 ≤ 0.364).
But the extinction of 𝑥 species is observed in (0.364 < 𝑣2 < 0.4). Again, Fig. 5(d) shows the oscillatory behavior of 𝑥 species in
0.3654 ≤ 𝑣3 ≤ 30.
Fig. 6(a–d) represents a bifurcation diagram of 𝑦 species of system (3) with respect to prey refuge rate 𝑝, half saturation constant
𝑣1 , the death rate of predator 𝑣2 and conservation rate of prey 𝑣3 respectively. The oscillatory or periodic dynamics of 𝑦 species
have been seen in (0 ≤ 𝑝 < 0.14) in Fig. 6(a). Also stable steady state has been observed for 𝑦 species in (0.14 < 𝑝 ≤ 1.0). Again, from
Fig. 6(b) it is seen that the 𝑦 species undergoes high amplitude oscillation in (0 ≤ 𝑣1 < 0.1287) and stable steady state is observed in
(0.1287 < 𝑣1 < 0.4). Also, Fig. 6(c) shows periodic or oscillatory dynamics of 𝑦 species in (0 ≤ 𝑣2 ≤ 0.304) and the stable steady state
behavior is observed in 0.304 < 𝑣2 ≤ 0.4. Again, Fig. 6(d) shows periodic behavior in (0.3654 < 𝑣3 ≤ 30) and amplitude of oscillation
gradually decreases.
Fig. 7(a–d) represents a bifurcation diagram of 𝑧 species of system (3) with respect to prey refuge rate 𝑝, half saturation constant
𝑣1 , the death rate of predator 𝑣2 and conservation rate of prey 𝑣3 respectively. The oscillatory or periodic dynamics of 𝑧 species
have been seen in (0 ≤ 𝑝 < 0.14) and stable steady state has been observed for 𝑧 species within (0.14 < 𝑝 < 0.67) in Fig. 7(a). Also,
the extinction of 𝑧 species is observed in (0.67 < 𝑝 ≤ 1). Again, from Fig. 7(b) it is seen that the 𝑧 species undergoes high amplitude
oscillation in (0 ≤ 𝑣1 < 0.1287) and stable steady state is observed in (0.1287 < 𝑣1 < 0.4). But the extinction of 𝑧 species may occur if
𝑣1 ≥ 0.4. Also, Fig. 7(c) shows periodic or oscillatory dynamics of 𝑧 species in (0 ≤ 𝑣2 ≤ 0.304) and the stable steady state behavior
is observed in 0.304 < 𝑣2 ≤ 0.4. Again, Fig. 7(d) shows high amplitude oscillation in (0.3654 < 𝑣3 ≤ 30).

7. Conclusion

In this paper, a three species predator–prey mathematical model has been formulated considering interspecies competition
coefficient’s as nonlinear type. Conditions for positivity and boundedness of solutions have been derived. Different equilibrium
points of the proposed mathematical model are determined and their stability conditions are determined. Theoretical analysis of
Hopf bifurcation of our proposed model has been done. Also, transcritical bifurcation of our proposed system has been investigated
with respect to the inter-species competition coefficient (𝛾12 ) of 𝑦 on 𝑥 species, inter-species competition coefficient (𝛾21 ) of 𝑥
on 𝑦, prey refuge rate (𝑝), half saturation constant (𝑣1 ), the natural death rate of predator species (𝑣2 ) and conservation rate
of prey species (𝑣3 ) respectively. The persistence condition of all the three species is reported. Conditions for orbital stability of
periodic solution has been determined analytically. After analyzing the proposed model, it is observed that the nonlinear interspecies
competition coefficients of the species have a significant role in the stability of an ecological system. Depending on the competition
coefficients species extinction is also possible. The importance of prey refuge rate and half saturation constant are also reported in

13
P. Panja, S. Gayen, T. Kar et al. Results in Control and Optimization 8 (2022) 100153

this study. Complex dynamical behavior can be observed in the system for suitable biologically meaningful parameter values. We
can conclude form our theoretical and numerical investigations that the interspecies competition coefficients, prey refuge rate, half
saturation constant, the death rate of predator species and conservation rate of prey species have a significant role in the stability
of an ecological system.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared
to influence the work reported in this paper.

References

[1] Lotka AJ. Elements of physical biology. Baltimore: Williams and Wilkins; 1925.
[2] Volterra V. Variazioni e fluttuazioni del numero dindividui in specie animali conviventi. Mem Acad Lincei Roma 1926;2:31–113.
[3] Rosenzweig M, MacArthur R. Graphical representation and stability conditions of predator–prey interaction. Amer Nat 1963;97:209–23.
[4] May RM. On relationships among various types of population models. Amer Nat 1973;107:46–57.
[5] Hassell MP. The dynamics of arthropod predator–prey systems. Monogr Popul Biol 1978;13:1–237.
[6] Freedman HI, Waltman P. Mathematical analysis of some three-species food-chain models. Math Biosci 1977;33:257–76.
[7] Hastings A, Powell T. Chaos in a three-species food chain. Ecology 1991;72:896–903.
[8] Chattopadhyay J, Sarkar RR. Chaos to order: Preliminary experiments with a population dynamics models of three trophic levels. Ecol Model
2003;163:45–50.
[9] Jana A, Roy SK. Behavioural analysis of two prey-two predator model. Ecol Complex 2021;47:100942.
[10] Roy SK, Roy B. Analysis of prey-predator three species Fishery model with harvesting including prey refuge and migration. Int J Bifurc Chaos
2016;26(02):1650022.
[11] Jana A, Roy SK. Holling-tanner prey-predator model with Beddington–DeAngelis functional response including delay. Int J Modelling Simul 2022;42:86–100.
[12] Panja P, Mondal SK. Stability analysis of coexistence of three species prey-predator model. Nonlinear Dynam 2015;81:373–82.
[13] Panja P, Mondal SK, Jana DK. Effects of toxicants on Phytoplankton-Zooplankton-Fish dynamics and harvesting. Chaos Solitons Fractals 2017;104:389–99.
[14] Panja P, Poria S, Mondal SK. Analysis of a harvested tritrophic food chain model in the presence of additional food for top predator. Int J Biomath
2018;11. [Link]
[15] Ghosh J, Sahoo B, Poria S. Prey-predator dynamics with prey refuge providing additional food to predator. Chaos Solitons Fractals 2017;96:110–9.
[16] Khajanchi S, Banerjee S. Role of constant prey refuge on stage structure predator–prey model with ratio dependent functional response. Appl Math Comput
2017;314:193–8.
[17] Jana D, Banerjee A, Samanta GP. Degree of prey refuges: Control the competition among prey and foraging ability of predator. Chaos Solitons Fractals
2017;104:350–62.
[18] Haque M, Rahman S, Venturino E, Li B. Effect of a functional response-dependent prey refuge in a predator–prey model. Ecol Complex 2014;20:248–56.
[19] Goh B. Stability in models of mutualism. Amer Nat 1979;113:261–75.
[20] Dhakne M, Munde A. Stability analysis of mutualistic interactions among three species with limited resources for first species and unlimited resources for
second and third species. Differ Equ Dyn Syst 2012;20:405–14.
[21] Kumar NP, Pattabhiramacharyulu NC. A three species ecosystem consisting of a prey, predator and a host commensal to the prey. Int J Open Probl Comput
Math 2010;3:92–113.
[22] Gakkhar S, Priyadarshi A, Banerjee S. Role of protection in a tri-trophic food chain dynamics. J Biol Systems 2012;20:155–75.
[23] Taylor A, Crizer A. A modified Lotka–Volterra competition model with a non-linear relationship between species. Rose-Hulman Undergrad Math J
2005;6:1–14.
[24] Gavina MKA, et al. Multi-species coexistence in Lotka–Volterra competitive systems with crowding effects. Sci Rep 2018;8. [Link]
s41598-017-19044-9.
[25] Neuhauser BC, Pacala SW. An explicitly spatial version of the Lotka–Volterra model with interspecific competition. Ann Appl Probab 1999;9:1226–59.
[26] Mohammed-Awel J, Lazari A. Investigation of the qualitative behavior of the equilibrium points for a modified Lotka–Volterra model. Georgian J Sci
2009;67:46–60.
[27] Gakkhar S, Gupta K. A three species dynamical system involving prey-predation, competition and commensalism. Appl Math Comput 2016;273:54–67.
[28] Birkhoff G, Rota GC. Ordinray differential equations. Ginn Boston; 1982.
[29] Hassard B, Kazarinoff N, Wan Y. Theory and application of hopf bifurcation. London mathematical society lecture note series, Cambridge: Cambridge
University Press; 1981.

14