# Predator-prey dynamics pertaining to structuralizing predator species into
three stages coupled with maturation delay owing to juvenile hunting

Debasish Bhattacharjee1 Tapasvini Roy2 Santanu Acharjee3 Tarini Kumar Dutta4
1,2,3Department of Mathematics  Gauhati University  Assam  India
4Department of Mathematics  Assam Don Bosco University  Assam  India
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The predator-prey dynamic appertaining to two species is explored, wherein the
predator species is structured into different stages. As evidenced from
natural documentation, the immature predators possess the potential to predate
albeit not as competently as the adults. Nevertheless, this potentiality is
not acquired immediately after their incipience of life, hence, the immature
stage is branched off into the infant stage, the stage with extensive reliance
on the adults, and the juvenile stage, the stage with the potential to predate
but not to procreate. In this paper, this inaugural concept is coupled with
injuries in the juvenile stage as the repercussion of their incompetency in
predating, thereby ensuing a delay in their maturation. With the incentive to
investigate the ascendancy of these refinements over the whole system,
stability analyses along with various bifurcation analyses around the
equilibrium points of the system are corroborated. In addition to Hopf,
transcritical, and saddle node bifurcations, the existence of Bogdanov-Takens
point, cusp point, Bautin bifurcation point, ‘bloom’ phenomenon, twice
occurring Hopf bifurcation, and bi-stability phenomenon make the paper
appreciably more rich and efficacious.

2020 AMS Classifications: 34C23,92D25,92D40.
Keywords:  Population dynamics; dangerous prey; limit cycle.

##  1 Introduction

After the inaugural work of Lotka and Volterra [1, 2], a wealth of theoretical
and pragmatical research has been deliberated upon the dynamical synergy among
the interacting species regarding predation, competition, parasitism, and
other such analogous symbiosis, owing to its ubiquitous influence on other
copious aspects of the natural world. For representing these population
dynamics of the interactions, mathematical models are contrived, based upon
the experiments and surveillance, so as to prognose the distribution of
population and structure of communities, and also how miniscule alteration in
biological or environmental parameters effectuates the entire bio-system.

The standard Lotka-Volterra mathematical model has been found to be lacking in
several aspects. With the exertion of many researchers, the models epitomizing
the predator-prey relationship have become a lot more realistic and precise by
assimilating instances such as, time-delay, functional response, stage-
structure, and others into the models. Stage structured mathematical models
paint more realistic and precise scenarios as compared to that of unstructured
ones; after the assimilation, the models are known to show substantial changes
in their behaviors. Owing to the improvisation on the assumption that each and
every individual of the predator species has the same ability for predation
and survivability, the stage-structured predator-prey models have received
much attention [26, 27, 28, 29, 30, 31]. The rudimentary works of constructing
the mathematical models with stage structures are done by Gurney [3], Nisbet
et al [4, 5], and their works have been further progressed by Aiello et al [6,
7].

The life-history of predators is always divided into two stages: the immature
stage and the matured stage. The immature predators are predominantly assumed
to not being able to hunt or reproduce and are exhaustively dependent on the
adults for their sustainability. However, there are several documentations of
predators hunting as immatures, albeit not as efficiently as the matured ones.
Caimans[10], pumas[8], snakes[9], snails[12], sharks[15], Clinus
superciliosus[13], spiders[11], and raptors[14] are only a handful among the
myriad of such predators that hunt during their developing period of lives. At
the same time, it is to be noted that the immature predators do not
instantaneously start hunting from the incipience of their life. Puma (Puma
concolor) kittens are dependent on their mothers upto two years, before
dispersing away as juveniles or sub-adults [8], similarly snakes, spiders and
the likes begin their life as an egg. The same fact of there being a period of
time before the juveniles of all other predator species start hunting, can be
unearthed. Thereupon, while considering predation by immature predators, it is
quite logical to further partition the immature stage into two stages: the
infancy stage and the juvenile stage. Eliciting from all these, we have
reckoned with the idea of structuring the predator species into three stages:
a stage with its individuals unmitigatedly reliant on the adults (infant
stage), a stage whose individuals are equipped for only predation (juvenile
stage), and a stage with its individuals equipped with the ability to both
predate and proliferate (adult stage); which (to the best of our knowledge)
has not yet been done.

Predatism, being a highly complex and dynamic interaction, is not always a
simple case scenario of predator successfully killing their prey. Prey
exhausts all of their morphological and behavioural traits that enhance their
survivability during their encounter with the predators. There are countless
instances where failure on the part of predators can cost more than just a
missed opportunity to feed, specifically yielding injury and sometimes even
death. Caribou and other adult ungulates, pose an injury risk to coyotes,
despite this, coyotes single out these large prey owing to high energetic
gains [16]. The same scenario is witnessed between the wolves and their prey,
elk and bison [17]. Mountain lions may be pierced by porcupine quills or left
blind in one eye from the blow of deer hoof [18]. Jaguars are frequently seen
with fresh wounds owing to the defensive attack by peccary or giant anteater
(Myrmecophaga tridactyla) [19]. Praying Mantids vomited and showed signs of
poisoning after feeding on $O.fasciatus$ [20].

The predators that include dangerous prey in their diet evolve their hunting
strategy accordingly so as to mitigate the risk of injury [21]. Interestingly,
it is found that the younger the predator is, the greater would be the
prospect of catastrophe of predation, owing to their relatively less body size
and inexperience when compared to that of the adult individuals of their
population [22]. The ingenuousness of the juveniles towards the specific
approach orientations required while handling dangerous prey leads to them
being injured, and sometimes even killed. Thereon, the predatory behavior and
its effectiveness varies with ontogeny.

Ants, considered as dangerous prey for spiderlings because of their mandibles
and the secretion of noxious substances by stinging, lacerates the
spiderlings, which has serious implications on autonomy. Spiderlings change
their attacking tactic through experience to that of the adults, thereby
capturing their prey faster and more efficiently [23, 24]. Centipedes, capable
of delivering extremely damaging and painful strikes using their forcipules,
are known to be fed upon by snakes, which hunt as juveniles [9]. Juvenile
migratory birds, being inadequate, require stop-over sites to replenish their
fuel reserves. The unfamiliar terrain as well as handling and capturing
unfamiliar prey poses a huge risk to them [25]. Serious injuries associated
with prey capture have been witnessed among the likes of big cats such as
lions, tigers, jaguars, cougars, etc, especially in the juveniles of the
populations [22]. The injuries suffered by the juvenile predators as a
consequence of their incompetency may lead to them being unwilling/unable to
hunt while nursing the wound or some other foraging costs similar to those
arising from the risk of predation, and with similar consequences, which
inevitably would lead to a delay in their transition into the adult stage.

In [26, 27, 28, 29], the authors have studied the mathematical models with
general stage structure in the predator species, where the immature predators
are plagued with mortal perilousness while preying on dangerous prey. But in
this paper, rather than death, the concentration is solely kept on injury as
being the negative repercussion of juvenile hunting (as the frequency of
injuries is much higher than that of death), which is yet to be probed into
mathematically. We have endeavoured towards assembling a system of equations
where only the adult predators are fully competent in predating and can
reproduce; the incompetency of the juvenile predator, while predating, in
handling of the dangerous prey causes a delay in their transition rate; during
the infancy period of life, the predators are completely dependent on the
adults. The structuralization of predator species and the inclusion of injury
induced delay import the significance of each stage, on the grounds of their
individually disparate impact on the whole ecosphere.

In this paper, the construction of the equations of bio-system is done in
section 2 and its enclosedness is shown in section 3. The discussion of the
equilibrium points in conjunction with their local and global stabilities is
done in section 4; the same is done for Hopf, saddle node, and transcritical
bifurcations in section 5. Numerical simulation is done in section 6, and
finally the conclusion is given in section 7.

##  2 Construction of the equations of bio-system

To begin with, we have considered the prey species to have a logistic growth
rate, a biologically feasible scenario owing to the constraints on natural
resources, along with the existence of the predator species in the same
habitat as that of the former. The functional response, which symbolises the
interaction between both the species is taken to be Holling-type 1 functional
response. Incorporating stage structure of the predator species into the
mathematical model in the generally used manner [36, 37, 38], and taking
$X(T)$ as the biomass density of the prey and $Y_{i}(T)$, and $Y_{m}(T)$ as
the immature and matured predators’ biomass densities respectively, we have,

$\begin{split}\frac{dX}{dT}&=rX(1-\frac{X}{K})-AXY_{m},\\\
\frac{dY_{i}}{dT}&=U_{1}AXY_{m}-B_{1}Y_{i}-D_{1}Y_{i},\\\
\frac{dY_{m}}{dT}&=B_{1}Y_{i}-D_{2}Y_{m},\\\
&X(0)>0,Y_{i}(0)>0,Y_{m}(0)>0.\end{split}$ (2.1)

Here, matured predators can kill their prey hence, their rate of predation is
taken as $A$, also they possess the reproductive ability which can be
translated as the transformation of consumed prey into immature predators,
$U_{1}$ denotes this rate of transformation. $D_{1},\;\text{and }D_{2}$ are
the death rates of immature and matured predators, respectively, and $B_{1}$
is the transition rate of immature predators into the matured stage, which is
taken as constant when the nutrient required by the immature predators is not
dependent on the prey. Even among the specialist predator, diet is known to
vary during ontogeny, and the immatures can comfortably survive for a time
being on alternate prey with low energetic returns. Conjointly, the immature
predator is completely dependent on their adults.

Eliciting from all the theoretical results as given in the introduction, we
have partitioned the immature stage of predator into infancy and juvenile
stages. The infancy stage is similar to that of the immature stage of 2.1, as
the individuals of this stage are exhaustively reliant on the adults. While
the juvenile predators can hunt, they are yet to obtain the ability to
reproduce and hence, their predation do not correspond to an increase in their
total numbers. Thus, the mathematical model 2.1 is remodelled into the
following system of equations:

$\begin{split}\frac{dX}{dT}&=rX(1-\frac{X}{K})-A_{1}XY_{2}-A_{2}XY_{3},\\\
\frac{dY_{1}}{dT}&=uA_{2}XY_{3}-BY_{1}-D_{1}Y_{1},\\\
\frac{dY_{2}}{dT}&=BY_{1}-CY_{2}-D_{2}Y_{2},\\\
\frac{dY_{3}}{dT}&=CY_{2}-D_{3}Y_{3},\\\
&X(0)>0,Y_{1}(0)>0,\,Y_{2}(0)>0,\,Y_{3}(0)>0.\end{split}$ (2.2)

Here, $Y_{1}=$ infant predators, $Y_{2}=$ juvenile predators, $Y_{3}=$ adult
predators, $u=$ transformation rate of consumed prey into infant predators,
$A_{1}$= predation rate of juvenile predators, $A_{2}=$ predation rate of
adult predator. $D_{1}=$ natural death rates of infant predators, $D_{2}=$
natural death rates of juvenile predators, $D_{3}=$ natural death rates of
adult predators, $B=$ transition rate from infancy to juvenile, $C=$
transition rate from juvenile to adulthood.

Due to lack of efficiency in predating by the juveniles, consequential to lack
of robustness and experience, they are faced with injuries. The cost of being
injured leads to them taking a reprieve from predation or from hunting their
required amount of nutritive prey, resulting in causing a delay in their
transition rate. This delay has been incorporated into the model 2.2, by
remoulding $C$ into $(C-A_{3}X)$, where $A_{3}$ is the rate of injury caused
by the prey.
Therefore, our equations of bio-system is:

$\begin{split}\frac{dX}{dT}&=rX(1-\frac{X}{K})-A_{1}XY_{2}-A_{2}XY_{3},\\\
\frac{dY_{1}}{dT}&=uA_{2}XY_{3}-BY_{1}-D_{1}Y_{1},\\\
\frac{dY_{2}}{dT}&=BY_{1}-(C-A_{3}X)Y_{2}-D_{2}Y_{2},\\\
\frac{dY_{3}}{dT}&=(C-A_{3}X)Y_{2}-D_{3}Y_{3},\\\ &\text{With initial
conditions: }X(0)>0,\hskip 3.0ptY_{1}(0)>0,\hskip 3.0ptY_{2}(0)>0,\hskip
3.0ptY_{3}(0)>0.\end{split}$ (2.3)

To make the model simpler for working out various calculations, the following
steps are taken:
$t=rT,\hskip 3.0ptx=X/K,\hskip 3.0pty_{1}=Y_{1}/K,\hskip
3.0pty_{2}=Y_{2}/K,\hskip 3.0pty_{3}=Y_{3}/K$ along with
$a_{1}=A_{1}\frac{K}{r},\hskip 3.0pta_{2}=A_{2}\frac{K}{r},\hskip
3.0ptb=B/r,\hskip 3.0ptd_{1}=D_{1}/r,\hskip 3.0ptc=C/r,\hskip
3.0pta_{3}=A_{3}\frac{K}{r},\hskip 3.0ptd_{2}=D_{2}/r,\hskip
3.0ptd_{3}=D_{3}/r$ and then using them in 2.3 we have,

$\begin{split}\frac{dx}{dt}&=x(1-x)-a_{1}xy_{2}-a_{2}xy_{3},\\\
\frac{dy_{1}}{dt}&=ua_{2}xy_{3}-by_{1}-d_{1}y_{1},\\\
\frac{dy_{2}}{dt}&=by_{1}-(c-a_{3}x)y_{2}-d_{2}y_{2},\\\
\frac{dy_{3}}{dt}&=(c-a_{3}x)y_{2}-d_{3}y_{3},\\\ &\text{With initial
conditions: }x(0)>0,\hskip 3.0pty_{1}(0)>0,\hskip 3.0pty_{2}(0)>0,\hskip
3.0pty_{3}(0)>0.\end{split}$ (2.4)

##  3 Boundedness of the equations of bio-system

From the first equation of (2.4), we can observe $\dot{x}\leq x(1-x)$, i.e.,
$\text{Hence }\lim_{t\to\infty}\,sup\,x(t)\leq 1$.
Next, we consider a variable $\beta$ in the following manner,

$\displaystyle\beta=$ $\displaystyle ux+y_{1}+y_{2}+y_{3},$ differentiating
$\beta$ with respect to time t and then adding $\zeta\beta$, where
$\zeta=min\,\\{d_{1},d_{2},d_{3}\\}$, we have,
$\displaystyle\frac{d\beta}{dt}+\zeta\beta=$ $\displaystyle
ux(1-x)-ua_{1}xy_{2}+\zeta
ux-(d_{1}-\zeta)y_{1}-(d_{2}-\zeta)y_{2}-(d_{3}-\zeta)y_{3}\leq ux(1-x+\zeta)$
$\displaystyle i.e.\,\frac{d\beta}{dt}+\zeta\beta\leq$ $\displaystyle
u(\frac{1+\zeta}{2})^{2}$ $\displaystyle i.e.,\beta\leq$ $\displaystyle
u\frac{(1+\zeta)^{2}}{4\zeta}(1-e^{-\zeta
t})+\beta(x(0),y_{1}(0),y_{2}(0),y_{3}(0))e^{-\zeta t}$
$\displaystyle\therefore\lim_{t\to\infty}\,sup\,\beta(t)$
$\displaystyle\leq\frac{u(1+\zeta)^{2}}{4\zeta}$
$\displaystyle\,i.e.\,\lim_{t\to\infty}\,sup\,(y_{1}(t)+y_{2}(t)+y_{3}(t))\leq\frac{u(1+\zeta)^{2}}{4\zeta}-u\leq\frac{u(1-\zeta)^{2}}{4\zeta}.$

Therefore, we arrive at the following theorem.

######  Theorem 1.

Every solution of the equations of bio-system (2.4) originating from the
region $R_{+}^{4}$ is confined in
$\\{(x,y_{1},y_{2},y_{3}):0\leq x\leq
1,0\leq(y_{1}+y_{2}+y_{3})\leq\frac{u(1-\zeta)^{2}}{4\zeta}\\}$ where
$\zeta=min\\{d_{1},d_{2},d_{3}\\}$.

##  4 Equilibrium points and their stabilities

Equilibrium points are the stationary points experienced by the system of
equations where the augmentation and reduction velocities of the species
vanish. Our given equations of bio-system (2.4) is encountered with four types
of such situations. First of them is the extinction equilibrium point
$E_{1}(0,0,0,0)$, and it always exists. In the absenteeism of predators, the
prey can attain the highest ecological load of its environment, resulting in
the axial equilibrium $E_{2}(1,0,0,0)$, the prey-only equilibrium point, which
too always exists, also with the predator being a specialist type, other axial
equilibriums are not possible. The most advantageous and desirable case
scenario in the ecological system is the co-existence of all the species. The
equations of bio-system (2.4) manifests two compresent equilibrium points
$E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$ and
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$. Stability analysis is
requisite for the alteration of behavior of the equilibrium points into a
desirable one. With this inducement, the theorems in this section have been
deliberated upon.

###  4.1 Extinction equilibrium point $E_{1}(0,0,0,0)$

The equilibrium point $E_{1}(0,0,0,0)$ is a saddle point because:

$\displaystyle
J=\left(\begin{array}[]{cccc}-2x-a_{1}y_{2}-a_{2}y_{3}+1&0&-a_{1}x&-a_{2}x\\\
a_{2}uy_{3}&-b-d_{1}&0&a_{2}ux\\\ a_{3}y_{2}&b&-c-d_{2}+a_{3}x&0\\\
-a_{3}y_{2}&0&c-a_{3}x&-d_{3}\\\ \end{array}\right)$ (4.5)

is the Jacobian matrix of the equations of bio-system (2.4), the
characteristic polynomial of which at $(0,0,0,0)$ is
$(1-\omega)(\omega+b+d_{1})(cd_{3}+d_{2}d_{3}+c\omega+d_{2}\omega+d_{3}\omega+\omega^{2})=0,\;\;i.e.,\;\;\omega=1,-b-d_{1},-c-d_{2},-d_{3}$.
$\therefore$ All the eigen values are negative except one.

###  4.2 Prey-only equilibrium point $E_{2}(1,0,0,0)$

######  Theorem 2.

The equilibrium point $E_{2}(1,0,0,0)$ can be either a center, stable point,
or an unstable point. All the conditions required for those related to
stability are given within the following proof.

Proof: Putting $(x,y_{1},y_{2},y_{3})=(1,0,0,0)$ at the Jacobian matrix (4.5),
we arrive at the required matrix $J_{1}$, the determinant and trace of whose
are:
$Det\,[J_{1}]=d_{1}(c+d_{2})d_{3}+b(cd_{3}+d_{2}d_{3}-a_{2}cu)-a_{3}(d_{1}d_{3}+b(d_{3}-a_{2}u))$,
and
$Trace\,[J_{1}]=a_{3}-(1+b+d_{1}+c+d_{2}+d_{3})$.
For (1,0,0,0) to be a saddle point, $Det\,[J_{1}]<0$ is required,

$i.e.,\hskip 3.0pteither.\hskip
5.0pta_{3}>c+d_{2},a_{2}<\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{b(a_{3}-c)u}\hskip
5.0ptor\hskip
4.0pta_{3}<c,a_{2}>\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{b(a_{3}-c)u}.$

For (1,0,0,0) to be a unstable point, $Det\,[J_{1}]>0$, and $Trace\,[J_{1}]>0$
are required,

$i.e.,\hskip 3.0pta_{3}>1+b+c+d_{1}+d_{2}+d_{3},\text{ and
}a_{2}>\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{bu(a_{3}-c)}.$

For (1,0,0,0) to be a centre, $Det\,[J_{1}]>0$, and $Trace\,[J_{1}]=0$ are
required,

$i.e.,\hskip
3.0pta_{2}>\frac{(b+d_{1})(1+b+d_{1}+d_{3})d_{3}}{b(1+b+d_{1}+d_{2}+d_{3})u},\text{
and }a_{3}=1+b+c+d_{1}+d_{2}+d_{3}.$

For (1,0,0,0) to be a stable point, $Det\,[J_{1}]>0$, and $Trace\,[J_{1}]<0$
are required, resulting the following cases:

$\displaystyle Case1:$ $\displaystyle\hskip 5.0pta_{3}<c,\text{ and
}a_{2}<\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{bu(a_{3}-c)},$ $\displaystyle
Case2:$ $\displaystyle\hskip 5.0pta_{3}=c,$ $\displaystyle Case3:$
$\displaystyle\hskip 5.0ptc<a_{3}<1+c,\text{ along with, either
}\left(a_{2}>\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{bu(a_{3}-c)},a_{3}>c+d_{2}\right)\text{
or }\left(a_{3}\leq c+d_{2}\right),$ $\displaystyle Case4:$
$\displaystyle\hskip
5.0pta_{3}>1+c,a_{2}>\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{bu(a_{3}-c)},\text{
and }a_{3}<1+c+b+d_{1}+d_{2}+d_{3}.$

(a)

(b)

Figure 1: Local and global stability of prey-only equilibrium point
$(1,0,0,0)$: (a)depicts extinction of juvenile and adult predators, (b)depicts
extinction of infant and juvenile predators.

######  Theorem 3.

The prey-only equilibrium point $E_{2}(1,0,0,0)$ is globally asymptotically
stable when $\;(1-x)(a_{1}y_{2}+x-1)+a_{2}y_{3}<0$.

Proof: Considering the Lyapunov function
$L=\tau_{1}(x-x_{1}-x_{1}ln\frac{x}{x_{1}})+\tau_{2}y_{1}+\tau_{3}y_{2}+\tau_{4}y_{3}$
so as to investigate the global stability of prey-only equilibrium point
$E_{1}(1,0,0,0)$, where $\tau_{i}(i=1,2,3,4)$ are positive constants, which
are to be taken accordingly.
Differentiating the Lyapunov function with respect to time t, we have,

$\displaystyle\frac{dL}{dt}=$
$\displaystyle\tau_{1}\frac{x-x_{1}}{x}\dot{x}+\tau_{2}\dot{y_{1}}+\tau_{3}\dot{y_{2}}+\tau_{4}\dot{y_{3}}$
Now, putting $x_{1}=1$ and $\tau_{1}=u\tau_{2},\tau_{2}=\tau_{3}=\tau_{4}$, we
have, $\displaystyle\frac{dL}{dt}=$
$\displaystyle\tau_{2}(-(1-x)^{2}-a_{1}xy_{2}+a_{1}y_{2}+a_{2}y_{3}-d_{1}y_{1}-d_{2}y_{2}-d_{3}y_{3})$
$\displaystyle\leq$
$\displaystyle\tau_{2}\\{(1-x)(a_{1}y_{2}-1+x)+a_{2}y_{3}\\}.$

Therefore, the sufficient condition for the equilibrium point to be globally
asymptotically stable would be:

$(1-x)(a_{1}y_{2}+x-1)+a_{2}y_{3}<0.$

###  4.3 Compresent equilibrium points
$E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$ and
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$

The two compresent equilibrium points manifested by the equations of bio-
system (2.4) are:

$\displaystyle
E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})=\left(\frac{\delta_{3}+\mu}{2a_{2}a_{3}bu},-\frac{(\delta_{4}-2a_{2}bd_{2}u+\mu)\mu_{1}}{4a_{2}^{2}a_{3}b^{2}u\delta_{1}+\delta_{2})},\frac{\mu_{1}}{2a_{2}a_{3}\delta_{1}+\delta_{2})},-\frac{(\delta_{4}+\mu)\mu_{1}}{4a_{2}^{2}a_{3}bd_{3}u(\delta_{1}+\delta_{2})}\right),$
$\displaystyle
E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})=\left(\frac{\delta_{3}-\mu}{2a_{2}a_{3}bu},\frac{(-\delta_{4}+2a2bd_{2}u+\mu)\mu_{2}}{4a_{2}^{2}a_{3}b^{2}u(\delta_{1}+\delta_{2})},\frac{\mu_{2}}{2a_{2}a_{3}(\delta_{1}+\delta_{2})},\frac{(-\delta_{4}+\mu)\mu_{2}}{4a_{2}^{2}a_{3}bd_{3}u(\delta_{1}+\delta_{2})}\right),$

where,
$\mu=\sqrt{a_{3}^{2}(b+d_{1})^{2}d_{3}^{2}-2a_{2}a_{3}b(b+d_{1})(c+2d_{2})d_{3}u+a_{2}^{2}b^{2}c^{2}u^{2}}$,
$\mu_{1}=a_{2}^{2}b(a_{3}-c)cu–a_{1}d_{3}(a_{3}(b+d_{1})d_{3}+\mu)+a_{2}(-a_{3}^{2}(b+d_{1})d_{3}-c(a_{1}bd_{3}u+\mu)+a_{3}(cd_{1}d_{3}+2d_{1}d_{2}d_{3}+bd_{3}(c+2(d_{2}+a_{1}u))+\mu))$,
$\mu_{2}=a_{2}^{2}b(a_{3}-c)cu+a_{1}d_{3}(-a_{3}(b+d_{1})d_{3}+\mu)+a_{2}(-a_{3}^{2}(b+d_{1})d_{3}+a_{3}(cd_{1}d_{3}+2d_{1}d_{2}d_{3}+bd_{3}(c+2(d_{2}+a_{1}u))-\mu)+c(-a_{1}bd_{3}u+\mu)),$
$\delta_{1}=(a_{1}d_{3}(-a_{3}(b+d_{1})+a_{1}bu),\;\delta_{2}=a_{2}(a_{3}(b+d_{1})d_{2}+a_{1}bcu),\;\delta_{3}=a_{3}(b+d_{1})d_{3}+a_{2}bcu,\;\text{and
}\delta_{4}=a_{3}(b+d_{1})d_{3}–a_{2}bcu$.

Existence
For both the compresent equilibrium points $E_{3}$, and $E_{4}$ to exist, the
following necessary and sufficient conditions should be satisfied:

$\displaystyle\hskip
28.45274ptd_{2}<\frac{(a_{3}(b+d_{1})d_{3}-a_{2}bcu)^{2}}{4a_{2}a_{3}b(b+d_{1})d_{3}u}\text{
along with}$ $\displaystyle Case1:\hskip 5.0pt$ $\displaystyle a_{3}=c,\text{
and }a_{2}>\frac{a_{3}bd_{3}+a_{3}d_{1}d_{3}}{2a_{3}bu-bcu},$ $\displaystyle
Case2:\hskip 5.0pt$ $\displaystyle a_{3}>c,\text{ and
}a_{2}bcu>a_{3}(b+d_{1})d_{3},$ $\displaystyle Case3:\hskip 5.0pt$
$\displaystyle a_{3}<c,\;2a_{3}>c,\;a_{2}>\frac{a_{3}d_{3}(b+d_{1})}{2a_{3}bu-
bcu},\text{ and
}\frac{(a_{3}-c)(d_{1}d_{3}+b(d_{3}-a_{2}u))}{(b+d_{1})d_{3}}<d_{2}.$

######  Theorem 4.

The compresent equilibrium points $E_{3}$ and $E_{4}$ are locally
asymptotically stable iff
$\epsilon_{1}>0,\;\epsilon_{4}>0,\;\epsilon_{1}\epsilon_{2}-\epsilon_{3}>0,\;\text{
and
}\;\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}>0.$
The values of $\epsilon_{1},\,\epsilon_{2},\,\epsilon_{3},\,\epsilon_{4}$ are
explained within the proof.

Proof: The characteristic equation of the Jacobian (4.5) is:

$\displaystyle a_{2}bux$
$\displaystyle[a_{3}x(-1+\omega+2x)-c(-1+\omega+2x+a_{1}y_{2})+a_{1}(d_{3}+\omega)y_{3}]+(-b-d_{1}-\omega)[a_{2}a_{3}(d_{2}+\omega)xy_{2}$
$\displaystyle+(-d_{3}-\omega)(a_{1}a_{3}xy_{2}+(c+d_{2}+\omega-
a_{3}x)(-1+\omega+2x+a_{1}y_{2}+a_{2}y_{3}))]=0.$

This is rewritten as:

$\displaystyle\omega^{4}+\epsilon_{1}\omega^{3}+\epsilon_{2}\omega^{2}+\epsilon_{3}\omega+\epsilon_{4}=0.$
(4.6)

Using the Routh-Hurwitz criteria, and putting either
$(x,y_{1},y_{2},y_{3})=(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$, or
$(x,y_{1},y_{2},y_{3})=(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$ according to
the requirement, we would have $E_{3}$ or $E_{4}$ to be locally asymptotically
stable iff $\epsilon_{1}>0,\hskip 7.0pt\epsilon_{4}>0,\hskip
7.0pt\epsilon_{1}\epsilon_{2}-\epsilon_{3}>0,\hskip 7.0pt\text{and
}\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}>0$.

######  Theorem 5.

A sufficient condition for the equations of bio-system (2.4) to be globally
asymptotically stable around the compresent equilibrium point
$E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*}),$ or
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$ is:
$(x-X^{*})(1-x-a_{1}y_{2})+a_{2}X^{*}y_{3}-a_{2}xy_{3}\frac{Y_{1}^{*}}{y_{1}}<0,\;\;\;Y_{1}^{*}(d_{1}+b)<d_{1}y_{1}+d_{2}y_{2},\;\;\;c+d_{2}<a_{3}x+by_{1}/y_{2},\;\;\text{and
}\;Y_{3}^{*}(d_{3}-c\frac{y_{2}}{y_{3}}+a_{3}x\frac{y_{2}}{y_{3}})<d_{3}y_{3}$.

Proof:  The positive definite Lyapunov function to be considered to
investigate the global stability of the system of equations around the
equilibrium point $E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$ is taken as:

$L_{1}=\eta_{1}(x-x^{*}-x^{*}ln\frac{x}{x^{*}})+\eta_{2}(y_{1}-y_{1}^{*}-y_{1}^{*}ln\frac{y_{1}}{y_{1}^{*}})+\eta_{3}(y_{2}-y_{2}^{*}-y_{2}^{*}ln\frac{y_{2}}{y_{2}^{*}})+\eta_{4}(y_{3}-y_{3}^{*}-y_{3}^{*}ln\frac{y_{3}}{y_{3}^{*}}).$

Here $\eta_{1},\eta_{2},\eta_{3},\text{ and }\eta_{4}$ are all positive
constants that are to be determined. Similarly, for the equilibrium point
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$, the Lyapunov function would
be

$L_{2}=\eta_{1}(x-x^{**}-x^{**}ln\frac{x}{x^{**}})+\eta_{2}(y_{1}-y_{1}^{**}-y_{1}^{**}ln\frac{y_{1}}{y_{1}^{**}})+\eta_{3}(y_{2}-y_{2}^{**}-y_{2}^{**}ln\frac{y_{2}}{y_{2}^{**}})+\eta_{4}(y_{3}-y_{3}^{**}-y_{3}^{**}ln\frac{y_{3}}{y_{3}^{**}}).$

Because of their similarities, both the Lyapunov functions can be written as
follows:

$L^{*}=\eta_{1}(x-X^{*}-X^{*}ln\frac{x}{X^{*}})+\eta_{2}(y_{1}-Y_{1}^{*}-Y_{1}^{*}ln\frac{y_{1}}{Y_{1}^{*}})+\eta_{3}(y_{2}-Y_{2}^{*}-Y_{2}^{*}ln\frac{y_{2}}{Y_{2}^{*}})+\eta_{4}(y_{3}-Y_{3}^{*}-Y_{3}^{*}ln\frac{Y_{3}}{Y_{3}^{*}}),$

where, $(X^{*},Y^{*},Z^{*})$ represent either of the equilibrium points
$E_{3}$, $E_{4}$.
Therefore, differentiating $L^{*}$ with respect to time t,

$\displaystyle\frac{dL^{*}}{dt}=$
$\displaystyle\eta_{1}\frac{x-X^{*}}{x}\dot{x}+\eta_{2}\frac{y_{1}-Y_{1}^{*}}{y_{1}}\dot{y_{1}}+\eta_{3}\frac{y_{2}-Y_{2}^{*}}{y_{2}}\dot{y_{2}}+\eta_{4}\frac{y_{3}-Y_{3}^{*}}{y_{3}}\dot{y_{3}}$
$\displaystyle\text{Taking }\eta_{1}=u\eta_{2},\text{ and
}\eta_{2}=\eta_{3}=\eta_{4},\text{ we get,}$ $\displaystyle\frac{dL^{*}}{dt}=$
$\displaystyle\eta_{2}[u(x-x^{2}-X^{*}+xX^{*}-a_{1}xy_{2}+a_{1}X^{*}y_{2}-a_{2}xy_{3}+a_{2}X^{*}y_{3})+ua_{2}xy_{3}-ua_{2}xy_{3}\frac{Y_{1}^{*}}{y_{1}}-by_{1}$
$\displaystyle+bY_{1}^{*}-d_{1}y_{1}+d_{1}Y_{1}^{*}+by_{1}-by_{1}\frac{Y_{2}^{*}}{y_{2}}-(c-a_{3}x)y_{2}+(c-a_{3}x)Y_{2}^{*}-d_{2}y_{2}+d_{2}Y_{2}^{*}+(c-a_{3}x)y_{2}$
$\displaystyle-(c-a_{3}x)y_{2}\frac{Y_{3}^{*}}{y_{3}}-d_{3}y_{3}+d_{3}Y_{3}^{*}]$
$\displaystyle=$
$\displaystyle\eta_{2}[\\{u(x-x^{2}-X^{*}+xX^{*}-a_{1}xy_{2}+a_{1}X^{*}y_{2}+a_{2}X^{*}y_{3}-a_{2}xy_{3}\frac{Y_{1}^{*}}{y_{1}})\\}+\\{bY_{1}^{*}-d_{1}y_{1}+d_{1}Y_{1}^{*}-d_{2}y_{2}\\}$
$\displaystyle+\\{-by_{1}\frac{Y_{2}^{*}}{y_{2}}+(c-a_{3}x)Y_{2}^{*}+d_{2}Y_{2}^{*}\\}+\\{-(c-a_{3}x)y_{2}\frac{Y_{3}^{*}}{y_{3}}-d_{3}y_{3}+d_{3}Y_{3}^{*}\\}].$

Therefore, the sufficient condition requisite for $\dot{L^{*}}$ to be negative
would be

$(x-X^{*})(1-x-a_{1}y_{2})+a_{2}X^{*}y_{3}-a_{2}xy_{3}\frac{Y_{1}^{*}}{y_{1}}<0,\;\;Y_{1}^{*}(d_{1}+b)<d_{1}y_{1}+d_{2}y_{2},\;\;c+d_{2}<a_{3}x+by_{1}/y_{2},\;\text{and
}\;Y_{3}^{*}(d_{3}-c\frac{y_{2}}{y_{3}}+a_{3}x\frac{y_{2}}{y_{3}})<d_{3}y_{3}.$

Putting
$(X^{*},Y_{1}^{*},Y_{2}^{*},Y_{3}^{*})=(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$
or
$(X^{*},Y_{1}^{*},Y_{2}^{*},Y_{3}^{*})=(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$
accordingly, we would have
the sufficient requirements for global stability of equilibrium points
$E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$ and
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$.

(a)

(b)

Figure 2: Local and global stability of compresent equilibrium point:
(a)depicts co-existence of prey, infant predators, and juvenile predators,
(b)depicts co-existence of prey, juvenile predators, and adult predators.

##  5 Bifurcations

The focal point of this section is to discuss the existence of various
transition points where the behavior of equations of bio-system (2.4)
modulates. The following transition points are investigated:

###  5.1 Hopf bifurcation

The transition point whose existence causes the system’s stability to switch
and the birth or death of a periodic solution around an equilibrium point is
the Hopf bifurcation point.

First we strive towards the acquisition of the requisites for the
manifestation of Hopf bifurcation around either of the compresent equilibrium
points $E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$ and
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$ with respect to any one of
the parameters.
The characteristic equation of the variational matrix (4.5) of the equations
of bio-system (2.4), as given in (4.6) is

$\omega^{4}+\epsilon_{1}\omega^{3}+\epsilon_{2}\omega^{2}+\epsilon_{3}\omega+\epsilon_{4}=0.$

If the point $c_{h}$ of parameter c (transition rate of juvenile predator into
adult) is presumed to be a bifurcation point, then we make the following
assumptions:

$(a)\;\;\epsilon_{1}>0,\epsilon_{4}>0,\epsilon_{1}\epsilon_{2}-\epsilon_{3}>0,\hskip
14.22636pt(b)\;\;\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0,\hskip
14.22636pt(c)\;\;[\frac{d}{dc}(\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{3}-\epsilon_{1}^{2}\epsilon_{4})]_{c=c_{h}}\neq
0.$

Taking the roots of the characteristic equation (4.5) as: $-j,-k,-l,-m$ ; we
would have

$\epsilon_{1}=j+k+l+m,\hskip 19.91684pt\epsilon_{2}=jk+lm+(j+k)(l+m),\hskip
19.91684pt\epsilon_{3}=(j+k)lm+(l+m)jk,\hskip 19.91684pt\epsilon_{4}=jklm.$
$\text{From (b): }\hskip
28.45274pt\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0,\hskip
19.91684pti.e.,\;(j+k)(j+l)(j+m)(k+l)(k+m)(l+m)=0.$

Therefore, at least one of the pair of roots should be additionally reciprocal
to each other. Say $j+k=0$.

$\text{Hence,}\hskip
14.22636pt\epsilon_{1}\epsilon_{2}-\epsilon_{3}=lm(l+m),\hskip
19.91684pt\epsilon_{1}=l+m,\hskip 19.91684pt\epsilon_{4}=-j^{2}lm.$

Also, from $(a)$, we would have, $\hskip 11.38092ptlm(l+m)>0,\;\;l+m>0,\text{
and }-j^{2}lm>0$.
Therefore, no other pair of roots of the characteristic equation would be
additionally reciprocal to each other, and $Re(l)>0,\;Re(m)>0$ and $j,k$ are
the purely complex conjugate numbers, i.e., we would have two eigen values
whose real part are negative and one pair of purely complex conjugates as the
roots of the characteristic equation (4.6).
Next, because of the property of continuity of roots, $\exists$ interval
$(c_{h}-\varsigma,c_{h}+\varsigma)$ for some $\varsigma>0$, such that for $c$
belonging to this interval we would have,

$j(c)=\psi_{1}(c)+\psi_{2}(c)i,\;\text{ and }\;k(c)=\psi_{1}(c)-\psi_{2}(c)i.$

For the transversality condition, we need to prove
$[\frac{dRe(j)}{dc}]_{c=c_{h}}\neq 0$ and $[\frac{dRe(k)}{dc}]_{c=c_{h}}\neq
0$.
For this we substitute $j(c)=\psi_{1}(c)+\psi_{2}(c)i$ in the characteristic
equation (4.6) and we get the following:

$\displaystyle(\psi_{1}+\psi_{2}i)^{4}+\epsilon_{1}(\psi_{1}+\psi_{2}i)^{3}+\epsilon_{2}(\psi_{1}+\psi_{2}i)^{2}+\epsilon_{3}(\psi_{1}+\psi_{2}i)+\epsilon_{4}=0$
$\displaystyle i.e.$
$\displaystyle\;[4(\psi_{1}+\psi_{2}i)^{3}+3\epsilon_{1}(\psi_{1}+\psi_{2}i)^{2}+2\epsilon_{2}(\psi_{1}+\psi_{2}i)+\epsilon_{3}](\psi_{1}’+\psi_{2}’i)+[\epsilon_{1}’(\psi_{1}+\psi_{2}i)^{3}$
$\displaystyle\hskip
14.22636pt+\epsilon_{2}^{\prime}(\psi_{1}+\psi_{2}i)^{2}+\epsilon_{3}’(\psi_{1}+\psi_{2}i)+\epsilon_{4}’]=0\hskip
14.22636pt[\text{where
}\epsilon^{\prime}_{i}=\frac{d\epsilon_{i}}{dt},\;\psi_{i}^{\prime}=\frac{d\psi_{i}}{dt}].$
(5.1)

Also using
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0$
in the characteristic equation at $c=c_{h}$ we get,

$\omega^{4}+\epsilon_{1}\omega^{3}+(\frac{\epsilon_{3}}{\epsilon_{1}}+\frac{\epsilon_{1}\epsilon_{4}}{\epsilon_{3}})\omega^{2}+\epsilon_{3}\omega+\epsilon_{4}=0,\hskip
14.22636pti.e.\;(\omega^{2}+\frac{\epsilon_{3}}{\epsilon_{1}})(\omega^{2}+\epsilon_{1}\omega+\frac{\epsilon_{1}\epsilon_{4}}{\epsilon_{3}})=0.$

Utilizing $\psi_{2}=(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}$ and
$\psi_{1}=0$ at $c=c_{h}$ in 5.1, we have

$\displaystyle[-4i(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{3}{2}}-3\epsilon_{1}(\frac{\epsilon_{3}}{\epsilon_{1}})+2\epsilon_{2}(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}i+\epsilon_{3}](\psi_{1}’+\psi_{2}’i)+[-\epsilon_{1}’(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{3}{2}}i-\epsilon_{2}’(\frac{\epsilon_{3}}{\epsilon_{1}})+\epsilon_{3}’(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}i+\epsilon_{4}’]=0$
Now, separating real and imaginary parts, we get the following: Real part
$\displaystyle
4(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{3}{2}}\psi_{2}’-3\epsilon_{3}\psi_{1}’-2\epsilon_{2}(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}\psi_{2}’+\epsilon_{3}\psi_{1}’-(\frac{\epsilon_{3}}{\epsilon_{1}})\epsilon_{2}’+\epsilon_{4}=0$
$\displaystyle\implies$
$\displaystyle\psi_{2}’=\frac{2\epsilon_{3}\psi_{1}’+\epsilon_{3}\epsilon_{2}’/\epsilon_{1}-\epsilon_{4}’}{2(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}(2\frac{\epsilon_{3}}{\epsilon_{1}}-\epsilon_{2})}.$
and Imaginary part
$\displaystyle-4(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{3}{2}}\psi_{1}’–3\epsilon_{3}\psi_{2}’+2\epsilon_{2}(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}\psi_{1}’+\epsilon_{3}\psi_{2}’-\epsilon_{1}’(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{3}{2}}+(\frac{\epsilon_{3}}{\epsilon_{1}})^{\frac{1}{2}}\epsilon_{3}’=0$
$\displaystyle\implies$ $\displaystyle
2(\epsilon_{2}-2\frac{\epsilon_{3}}{\epsilon_{1}})\psi_{1}+\frac{-2\epsilon_{1}\epsilon_{3}\psi_{1}’-\epsilon_{3}\epsilon_{2}’+\epsilon_{1}\epsilon_{4}’}{2\frac{\epsilon_{3}}{\epsilon_{1}}-\epsilon_{2}}+(\epsilon_{3}’-\frac{\epsilon_{1}’\epsilon_{3}}{\epsilon_{1}})=0\hskip
10.0pt\text{[putting the value of $\psi_{2}’$]}$ $\displaystyle\implies$
$\displaystyle\psi_{1}’=-\epsilon_{1}\frac{\epsilon_{1}’(\epsilon_{2}\epsilon_{3}-2\epsilon_{1}\epsilon_{4})+\epsilon_{2}’\epsilon_{1}\epsilon_{3}+\epsilon_{3}’(-2\epsilon_{3}+\epsilon_{1}\epsilon_{2})-\epsilon_{4}’\epsilon_{1}^{2}}{2[(\epsilon_{3}-\epsilon_{2}\epsilon_{2})^{2}+\epsilon_{1}^{3}\epsilon_{3}]}.$

Therefore, transversality condition holds when

$\epsilon_{1}’(\epsilon_{2}\epsilon_{3}-2\epsilon_{1}\epsilon_{4})+\epsilon_{2}’\epsilon_{1}\epsilon_{3}+\epsilon_{3}’(-2\epsilon_{3}+\epsilon_{1}\epsilon_{2})-\epsilon_{4}’\epsilon_{1}^{2}\neq
0.$

i.e. when our third assumption holds. Hence, we arrive at the following
theorem:

######  Theorem 6.

The sufficient condition for the occurence of Hopf bifurcation in the given
equations of bio-system (2.4) with respect to any parameter, say $c$ is:

$(a)\;\;\epsilon_{1}>0,\epsilon_{4}>0,\epsilon_{1}\epsilon_{2}-\epsilon_{3}>0,\hskip
19.91684pt(b)\;\;\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0,\text{
and}\hskip
14.22636pt(c)\;\;[\frac{d}{dc}(\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{3}-\epsilon_{1}^{2}\epsilon_{4})]_{c=c_{h}}\neq
0.$

Now, we discuss the stability and orientation of the bifurcating periodic
trajectory

######  Theorem 7.

The orientation of the trajectory of Hopf bifurcation around a compresent
equilibrium point can be ascertained by the sign of $\vartheta$, positivity
would imply supercritical Hopf bifurcation while negativity would imply
subcritical Hopf bifurcation; also positive $\beta_{2}$ indicates unstability
while negative $\beta_{2}$ implies stability. The values of
$\vartheta,\beta_{2}$ are given within the proof.

Proof:  So as to arbitrate the stability and orientation of bifurcating
periodic solution occurring because of Hopf bifurcation, we use a procedure
homologous to the methodology given by Hassard et al [34]. At Hopf bifurcation
point, we must have a pair of conjugate imaginary eigen values, say $\pm\alpha
i$ and the other two eigen values are $v_{1},$ and $v_{2}$.
For the purpose of reducing the given equations of bio-system into normal
form, new variables $(\bar{x},\bar{y_{1}},\bar{y_{2}},\bar{y_{3}})$ are
introduced as: $x=\bar{x}+X^{*}$, $y_{1}=\bar{y_{1}}+Y_{1}^{*}$,
$y_{2}=\bar{y_{2}}+Y_{2}^{*}$, $y_{3}=\bar{y_{3}}+Y_{3}^{*}$. For sake of
simplicity while writing we take,
$x\to\bar{x},y_{1}\to\bar{y_{1}},y_{2}\to\bar{y_{2}},y_{3}\to\bar{y_{3}}$.
Also $(X^{*},Y_{1}^{*},Y_{2}^{*},Y_{3}^{*})$ is used in lieu of either
$E_{3}(x^{*},y_{1}^{*},y_{2}^{*},y_{3}^{*})$ or
$E_{4}(x^{**},y_{1}^{**},y_{2}^{**},y_{3}^{**})$.
After the introduction of variables, the system of equations (2.4) are:

$\displaystyle\dot{x}=$
$\displaystyle(x+X^{*})(1-x-X^{*})-a_{1}(x+X^{*})(y_{2}+Y_{2}^{*})-a2(x+X^{*})(y_{3}+Y_{3}^{*}),$
$\displaystyle\dot{y_{1}}=$ $\displaystyle
ua_{2}(x+X^{*})(y_{3}+Y_{3}^{*})-b(y_{1}+Y_{1}^{*})-d_{1}(y_{1}+Y_{1}^{*}),$
$\displaystyle\dot{y_{2}}=$ $\displaystyle
b(y_{1}+Y_{1}^{*})-(c-a_{3}(x+X^{*}))(y_{2}+Y_{2}^{*})-d_{2}(y_{2}+Y_{2}^{*}),$
$\displaystyle\text{and }\dot{y_{3}}=$
$\displaystyle(c-a_{3}(x+X^{*}))(y_{2}+Y_{2}^{*})-d_{3}(y_{3}+Y_{3}^{*}).$

The above system of equations can be written as:

$\displaystyle\dot{X}=UX+V,$ (5.2)

where, $X=\left(\begin{array}[]{cccc}x&y_{1}&y_{2}&y_{3}\\\
\end{array}\right)’,\;\;$ with $U$ being the linear part and $V$ the non-
linear part.
Next, a non-singular matrix $A=\left(\begin{array}[]{cccc}1&0&1&1\\\
a_{21}&a_{22}&a_{23}&a_{24}\\\ a_{31}&a_{32}&a_{33}&a_{34}\\\
a_{41}&a_{42}&a_{43}&a_{44}\\\ \end{array}\right)$ is procured so as to get
$A^{-1}UA=\left(\begin{array}[]{cccc}0&\alpha&0&0\\\ -\alpha&0&0&0\\\
0&0&v_{1}&0\\\ 0&0&0&v_{2}\\\ \end{array}\right).$
For this, we would need

$\displaystyle a_{21}=$
$\displaystyle[-u_{13}u_{24}\left(u_{24}u_{32}\left(\alpha^{2}-u_{11}u_{33}\right)+u_{14}\left(u_{21}u_{32}u_{33}+u_{22}u_{31}u_{33}+\alpha^{2}(-u_{31})\right)\right)+u_{14}\left(\alpha^{2}+u_{33}^{2}\right)(u_{24}\left(\alpha^{2}+u_{11}u_{22}\right)$
$\displaystyle-
u_{14}u_{21}u_{22})-u_{13}^{2}u_{24}^{2}u_{31}u_{32}]/[u_{13}^{2}u_{24}^{2}u_{32}^{2}+2u_{13}u_{14}u_{24}u_{32}\left(u_{22}u_{33}-\alpha^{2}\right)+u_{14}^{2}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right)],$
$\displaystyle a_{22}=$
$\displaystyle-\frac{\alpha\left(u_{13}u_{24}(u_{24}u_{32}(u_{11}+u_{33})+u_{14}(-u_{21}u_{32}+u_{22}u_{31}+u_{31}u_{33}))+u_{14}\left(\alpha^{2}+u_{33}^{2}\right)(u_{24}(u_{22}-u_{11})+u_{14}u_{21})\right)}{u_{13}^{2}u_{24}^{2}u_{32}^{2}+2u_{13}u_{14}u_{24}u_{32}\left(u_{22}u_{33}-\alpha^{2}\right)+u_{14}^{2}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right)},$
$\displaystyle a_{23}=$
$\displaystyle-\frac{(u_{33}-v_{1})(u_{24}(v_{1}-u_{11})+u_{14}u_{21})+u_{13}u_{24}u_{31}}{u_{13}u_{24}u_{32}+u_{14}(u_{22}-v_{1})(u_{33}-v_{1})},$
$\displaystyle a_{24}=$
$\displaystyle-\frac{(u_{33}-v_{2})(u_{24}(v_{2}-u_{11})+u_{14}u_{21})+u_{13}u_{24}u_{31}}{u_{13}u_{24}u_{32}+u_{14}(u_{22}-v_{2})(u_{33}-v_{2})},$
$\displaystyle a_{31}=$
$\displaystyle-[u_{13}u_{24}u_{32}(u_{11}u_{24}u_{32}-u_{14}u_{21}u_{32}+u_{14}u_{22}u_{31})+u_{14}(u_{24}u_{32}\left(u_{11}\left(u_{22}u_{33}-\alpha^{2}\right)+\alpha^{2}(u_{22}+u_{33})\right)+u_{14}(-u_{21}u_{22}u_{32}u_{33}$
$\displaystyle+\alpha^{2}(u_{21}u_{32}+u_{31}u_{33})+u_{22}^{2}u_{31}u_{33}))]/[u_{13}^{2}u_{24}^{2}u_{32}^{2}+2u_{13}u_{14}u_{24}u_{32}\left(u_{22}u_{33}-\alpha^{2}\right)+u_{14}^{2}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right)],$
$\displaystyle a_{32}=$
$\displaystyle-\frac{\alpha\left(u_{14}u_{24}u_{32}\left(\alpha^{2}+u_{11}(u_{22}+u_{33})-u_{13}u_{31}-u_{22}u_{33}\right)-u_{13}u_{24}^{2}u_{32}^{2}+u_{14}^{2}\left(-u_{21}u_{22}u_{32}-u_{21}u_{32}u_{33}+u_{22}^{2}u_{31}+\alpha^{2}u_{31}\right)\right)}{u_{13}^{2}u_{24}^{2}u_{32}^{2}+2u_{13}u_{14}u_{24}u_{32}\left(u_{22}u_{33}-\alpha^{2}\right)+u_{14}^{2}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right)},$
$\displaystyle a_{33}=$
$\displaystyle\frac{u_{24}u_{32}(v_{1}-u_{11})+u_{14}(u_{21}u_{32}-u_{22}u_{31}+u_{31}v_{1})}{u_{13}u_{24}u_{32}+u_{14}(u_{22}-v_{1})(u_{33}-v_{1})},$
$\displaystyle a_{34}=$
$\displaystyle\frac{u_{24}u_{32}(v_{2}-u_{11})+u_{14}(u_{21}u_{32}-u_{22}u_{31}+u_{31}v_{2})}{u_{13}u_{24}u_{32}+u_{14}(u_{22}-v_{2})(u_{33}-v_{2})},$
$\displaystyle a_{41}=$
$\displaystyle[u_{13}\left(u_{24}u_{32}\left(u_{11}\left(\alpha^{2}-u_{22}u_{33}\right)+\alpha^{2}(u_{22}+u_{33})\right)+u_{14}\left(-u_{21}u_{22}u_{32}u_{33}+\alpha^{2}(u_{21}u_{32}+u_{31}u_{33})+u_{22}^{2}u_{31}u_{33}\right)\right)-u_{11}u_{14}(\alpha^{2}$
$\displaystyle+u_{22}^{2})\left(\alpha^{2}+u_{33}^{2}\right)+u_{13}^{2}u_{24}u_{32}(u_{22}u_{31}-u_{21}u_{32})]/[u_{13}^{2}u_{24}^{2}u_{32}^{2}+2u_{13}u_{14}u_{24}u_{32}\left(u_{22}u_{33}-\alpha^{2}\right)+u_{14}^{2}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right)],$
$\displaystyle a_{42}=$
$\displaystyle[\alpha(u_{13}\left(u_{24}u_{32}\left(-\alpha^{2}+u_{11}(u_{22}+u_{33})+u_{22}u_{33}\right)+u_{14}\left(-u_{21}u_{22}u_{32}-u_{21}u_{32}u_{33}+u_{22}^{2}u_{31}+\alpha^{2}u_{31}\right)\right)-u_{13}^{2}u_{24}u_{31}u_{32}$
$\displaystyle+u_{14}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right))]/[u_{13}^{2}u_{24}^{2}u_{32}^{2}+2u_{13}u_{14}u_{24}u_{32}\left(u_{22}u_{33}-\alpha^{2}\right)+u_{14}^{2}\left(\alpha^{2}+u_{22}^{2}\right)\left(\alpha^{2}+u_{33}^{2}\right)],$
$\displaystyle a_{43}=$
$\displaystyle\frac{u_{13}(-u_{21}u_{32}+u_{22}u_{31}-u_{31}v_{1})-(u_{11}-v_{1})(v_{1}-u_{22})(v_{1}-u_{33})}{u_{13}u_{24}u_{32}+u_{14}(u_{22}-v_{1})(u_{33}-v_{1})},$
$\displaystyle a_{44}=$
$\displaystyle\frac{u_{13}(-u_{21}u_{32}+u_{22}u_{31}-u_{31}v_{2})-(u_{11}-v_{2})(v_{2}-u_{22})(v_{2}-u_{33})}{u_{13}u_{24}u_{32}+u_{14}(u_{22}-v_{2})(u_{33}-v_{2})}.$

$u_{ij}$ $(i,j=1,2,3,4)$ are the elements of the matrix U,
i.e.$U=[u_{ij}]_{4\times 4}$.
To obtain the normal form, another transformation $X=AZ$ is done, where
$Z=\left(\begin{array}[]{cccc}z_{1}&z_{2}&z_{3}&z_{4}\\\
\end{array}\right)^{\prime}$. Therefore, (5.2) becomes
$\dot{X}=UX+V\;\;\;\;\implies\dot{Z}=A^{-1}U(AZ)+A^{-1}V\;\;\;\;\implies\dot{Z}=A^{-1}UAZ+W$.
The matrix $W=\left(\begin{array}[]{cccc}w_{1}&w_{2}&w_{3}&w_{4}\\\
\end{array}\right)’$ is obtainable by transformation of V done by the
following
$x=z_{1}+z_{3}+z_{4},\;\;$
$y_{1}=a_{21}z_{1}+a_{22}z_{2}+a_{23}z_{3}+a_{24}z_{4},\;\;$
$y_{2}=a_{31}z_{1}+a_{32}z_{2}+a_{33}z_{3}+a_{34}z_{4},$ and
$y_{3}=a_{41}z_{1}+a_{42}z_{2}+a_{43}z_{3}+a_{44}z_{4}.$
Now, we arrive at the stage where we can derive the stability and orientation
of the bifurcating periodic solution with the help of the expressions,
$g_{11},g_{02},g_{20},G_{21},h^{1}_{11},h^{1}_{20},w_{11},w_{20},G_{110}^{1},G_{101}^{1},g_{21}$,
all of which can be obtained as is given by Hassard et al [34] at
$(z_{1},z_{2},z_{3},z_{4})=(0,0,0,0)$. Finally, we arrive at the required
expression for determining the direction of Hopf bifurcation,

$C_{1}(0)=\frac{i}{2\alpha}\left(g_{20}g_{11}-2|g_{11}|^{2}-\frac{1}{3}|g_{02}|^{2}\right)+\frac{g_{21}}{2}.$
$\therefore\;\;\vartheta=-\frac{ReC_{1}(0)}{Re\omega’(c_{h})},\text{
and}\hskip 5.69046pt\beta_{2}=2ReC_{1}(0).$

$\vartheta>0$ implies supercritical Hopf bifurcation and $\vartheta<0$ implies
subcritical Hopf bifurcation. Moreover, $\beta_{2}<0$ and $\beta_{2}>0$ would
mean stability and unstablity for the bifurcating periodic solution,
respectively.

(a) at $c=0.037$

(b) at $c=0.03598345$

(c) at $c=0.033$

Figure 3: Time series depicting equilibrium states in the neighbourhood of
Hopf bifurcation point of $c$, the transition rate of juvenile predators into
adult predators: $(a)c>c_{h}$, $(b)c=c_{h}$, $(c)c<c_{h}$, $c_{h}$ is the Hopf
bifurcation point.

###  5.2 Transcritical bifurcation

The transition point where two equilibria convene and the system experiences
changes in stability in regards to its equilibrium points is the transcritical
bifurcation point.

######  Theorem 8.

The equations of bio-system (2.4) experiences interchangeability of stability
between the compresent equilibrium point and prey-only equilibrium when
$a_{2}=[(b+d¬_{1})(a_{3}-c-d_{2})d_{3}]/[b(a_{3}-c)u]$.

Proof: The variational matrix of the equations of bio-system (2.4) at
$E_{2}(1,0,0,0)$ is:

$\left(\begin{array}[]{cccc}-1&0&-a_{1}&-a_{2}\\\ 0&-b-d_{1}&0&a_{2}u\\\
0&b&a_{3}-c-d_{2}&0\\\ 0&0&c-a_{3}&-d_{3}\\\ \end{array}\right)$ (5.3)

For $a_{2}=[(b+d_{1})(a_{3}-c-d_{2})d_{3}]/[b(a_{3}-c)u](=a_{2t})$, the
determinant of the above matrix would be zero, i.e., zero would be the value
of one of the eigen values. Corresponding to this zero eigen value for the
matrices (5.3) and its transpose at $a_{2t}$, the respective eigen vectors
are:
$l=\left(-a_{1}l_{3}+\frac{(b+d_{1})(a_{3}-c–d_{2})l_{3}}{bu},\frac{(-a_{3}+c+d_{2})l_{3}}{b},l_{3},\frac{(-a_{3}+c)l_{3}}{d3}\right)$,
and
$m=\left(0,m_{2},((b+d_{1})m_{2})/b),\frac{(b+d_{1})(a_{3}-c–d_{2})m_{2}}{b(a_{3}–c)}\right)$
respectively, where $l_{3}$ and $m_{2}$ are arbitrary
Next, the given model (2.4) is rewritten in the form

$\dot{Y}=BY,$

where $Y=\left(\begin{array}[]{c}x\\\ y_{1}\\\ y_{2}\\\ y_{3}\\\
\end{array}\right)$, and
$B=\left(\begin{array}[]{c}(1-x)x-a_{1}y_{2}x-a_{2}y_{3}x\\\
-by_{1}-d_{1}y_{1}+a_{2}uxy_{3}\\\ by_{1}-d_{2}y_{2}-(c-a_{3}x)y_{2}\\\
(c-a_{3}x)y_{2}-d_{3}y_{3}\\\ \end{array}\right)$.
According to the Sotomayor theorem [33], the instance of transcritical
bifurcation can be proven by

$mB_{a_{2}}(E_{2},a_{2t})=0$, $m(DB_{a_{2}}(E_{2},a_{2t}))\hskip
3.0ptl^{t}\neq 0$, $m[D^{2}(B_{a_{2}}(E_{2},a_{2t}))]\hskip
3.0pt(l^{t},l^{t})\neq 0$.

In our case, we have

$\displaystyle
B_{a_{2}}(E_{2},a_{2t})=\left(\begin{array}[]{cccc}\frac{\partial
B_{1}}{\partial a_{2}}&\frac{\partial B_{2}}{\partial a_{2}}&\frac{\partial
B_{3}}{\partial a_{2}}&\frac{\partial B_{4}}{\partial a_{2}}\\\
\end{array}\right)_{(1,0,0,0);a_{2t}}$ $\displaystyle\text{Thus,
}\;mB_{a_{2}}(E_{2},a_{2t})=0,$ $\displaystyle
DB_{a_{2}}(E_{2},a_{2t})=\left(\begin{array}[]{cccc}\frac{\partial^{2}B_{1}}{\partial
x\,\partial a_{2}}&\frac{\partial^{2}B_{1}}{\partial y_{1}\,\partial
a_{2}}&\frac{\partial^{2}B_{1}}{\partial y_{2}\,\partial
a_{2}}&\frac{\partial^{2}B_{1}}{\partial y_{3}\,\partial a_{2}}\\\
\frac{\partial^{2}B_{2}}{\partial x\,\partial
a_{2}}&\frac{\partial^{2}B_{2}}{\partial y_{1}\,\partial
a_{2}}&\frac{\partial^{2}B_{2}}{\partial y_{2}\,\partial
a_{2}}&\frac{\partial^{2}B_{2}}{\partial y_{3}\,\partial a_{2}}\\\
\frac{\partial^{2}B_{3}}{\partial x\,\partial
a_{2}}&\frac{\partial^{2}B_{3}}{\partial y_{1}\,\partial
a_{2}}&\frac{\partial^{2}B_{3}}{\partial y_{2}\,\partial
a_{2}}&\frac{\partial^{2}B_{3}}{\partial y_{3}\,\partial a_{2}}\\\
\frac{\partial^{2}B_{4}}{\partial x\,\partial
a_{2}}&\frac{\partial^{2}B_{4}}{\partial y_{1}\,\partial
a_{2}}&\frac{\partial^{2}B_{4}}{\partial y_{2}\,\partial
a_{2}}&\frac{\partial^{2}B_{4}}{\partial y_{3}\,\partial a_{2}}\\\
\end{array}\right)_{(1,0,0,0);a_{2t}}$ $\displaystyle\text{Thus,
}\;m(DB_{a_{2}}(E_{2},a_{2t}))l^{t}=m_{2}u\frac{(-a_{3}+c)l_{3}}{d3}\neq 0,$
$\displaystyle
D^{2}B(E_{2},a_{2t})(l^{t},l^{t})=\left(\begin{array}[]{c}\sum_{i=1}^{4}\sum_{j=1}^{4}[\frac{\partial^{2}B_{1}}{\partial
p_{i}\,\partial p_{j}}l_{i}l_{j}]\\\
\sum_{i=1}^{4}\sum_{j=1}^{4}[\frac{\partial^{2}B_{2}}{\partial p_{i}\,\partial
p_{j}}l_{i}l_{j}]\\\
\sum_{i=1}^{4}\sum_{j=1}^{4}[\frac{\partial^{2}B_{3}}{\partial p_{i}\,\partial
p_{j}}l_{i}l_{j}]\\\
\sum_{i=1}^{4}\sum_{j=1}^{4}[\frac{\partial^{2}B_{4}}{\partial p_{i}\,\partial
p_{j}}l_{i}l_{j}]\\\
\end{array}\right)^{\prime}_{(1,0,0,0);\;\text{{\small$p_{i}=(x,y_{1},y_{2},y_{3})$}}}$
Thus, $\displaystyle
mD^{2}B_{a_{2}}(E_{2},a_{2t})(l^{t},l^{t})=-\frac{2l_{3}^{2}m_{2}((b+d_{1})(a_{3}-c-d_{2})-a_{1}bu)(-a_{3}(b+d_{1})d_{2}d_{3}+a_{2}b(a_{3}-c)^{2}u)}{b^{2}(a_{3}-c)d_{3}u}\;\;\neq
0.$

Hence, the theorem is proved.

Note: From direct observation, it is to be noted that for the occurance of
trancritical bifurcation, we must have $a_{3}<c$ or $c+d_{2}<a_{3}$. Also in
lieu of $a_{2}$ other parametric values can be used as bifurcating parameter.

###  5.3 Saddle node bifurcation

The point of collision and disappearance of two branches of equilibria is the
saddle node bifurcation point (also known as limit point bifurcation).

######  Theorem 9.

The equations of bio-system (2.4) experiences saddle node bifurcation when the
parameters satisfy the expression
$a_{2}bux(a_{3}x(-1+2x)-c(-1+2x+a_{1}y_{2})+a_{1}d_{3}y_{3})+(-b-d_{1})[a_{2}a_{3}d_{2}xy_{2}-d_{3}(a_{1}a_{3}xy_{2}+(c+d_{2}-a_{3}x)(-1+2x+a_{1}y_{2}+a_{2}y_{3}))]=0$
around the equilibrium points, along with
$\chi_{2}^{T}[B_{b}(E_{4},b_{s})]\neq 0,\text{ and
}\chi_{2}^{T}[D^{2}B_{b}(E_{4},b_{s})(\chi_{1},\chi_{1})]\neq 0.$

Proof: For saddle node bifurcation to occur at an equilibrium point, one of
the eigen values of the variational matrix (4.5) of our model needs to be
zero, i.e., the determinant of (4.5) would be zero. Hence, we would need-

$a_{2}bux(a_{3}x(-1+2x)-c(-1+2x+a_{1}y_{2})+a_{1}d_{3}y_{3})+(-b-d_{1})[a_{2}a_{3}d_{2}xy_{2}-d_{3}(a_{1}a_{3}xy_{2}+(c+d_{2}-a_{3}x)(-1+2x+a_{1}y_{2}+a_{2}y_{3}))]=0,$

where, $x,y_{1},y_{2},y_{3}$ are to be replaced by the points of equilibrium
where the bifurcation is supposed to occur. Any one parameter from the above
expression can be used as the control parameter, say $b=b_{s}$ is the point of
saddle node bifurcation.
Suppose $E_{4}$ is the equilibrium point, then we have

$\displaystyle\frac{1}{2a_{2}^{2}a_{3}^{2}b^{2}u^{2}}(a_{3}^{3}(b+d_{1})^{2}d_{3}^{2}((b+d_{1})d_{3}–a_{2}bu)+a_{2}^{2}b^{2}c^{2}u^{2}(-R+a_{2}bcu)+a_{2}a_{3}bu(2(b+d_{1})d_{2}d_{3}R+a_{2}bc(-(b+d_{1})(c+4d_{2})d_{3}$
$\displaystyle+R)u–a_{2}^{2}b^{2}c^{2}u^{2})–a_{3}^{2}(b+d_{1})d_{3}((b+d_{1})d_{3}R+a_{2}b((b+d_{1})(c+4d_{2})d_{3}-R)u-2a_{2}^{2}b^{2}(c+2d_{2})u^{2}))=0,$
$\displaystyle\hskip 28.45274pt\text{where,
}R=\sqrt{a_{2}^{2}b^{2}c^{2}u^{2}-2a_{2}a_{3}bd_{3}u(b+d_{1})(c+2d_{2})+a_{3}^{2}d_{3}^{2}(b+d_{1})^{2}}.$

To simplify the calculations, the matrix (4.5) at compresent equilibrium point
is be rewritten as

$\displaystyle\left(\begin{array}[]{cccc}e_{1}&0&e_{2}&e_{3}\\\
e_{4}&e_{5}&0&e_{6}\\\ e_{7}&e_{8}&e_{9}&0\\\ e_{10}&0&e_{11}&e_{12}\\\
\end{array}\right).$ (5.8)

Hence, the determinant of the above matrix would be zero, when
$e_{5}=\frac{(e_{12}e_{2}e_{4}–e_{11}e_{3}e_{4}+e_{1}e_{11}e_{6}–e_{10}e_{2}e_{6})e_{8}}{e_{12}e_{2}e_{7}–e_{11}e_{3}e_{7}–e_{1}e_{12}e_{9}+e_{10}e_{3}e_{9}}$.
Therefore, the eigen vectors corresponding to zero eigen value of matrix (5.8)
with the value of
$e_{5}=\frac{(e_{12}e_{2}e_{4}–e_{11}e_{3}e_{4}+e_{1}e_{11}e_{6}–e_{10}e_{2}e_{6})e_{8}}{e_{12}e_{2}e_{7}–e_{11}e_{3}e_{7}–e_{1}e_{12}e_{9}+e_{10}e_{3}e_{9}}$
and its transpose are:

$\displaystyle\chi_{1}=$
$\displaystyle\left(\phi,-\frac{(e_{12}e_{2}e_{7}–e_{11}e_{3}e_{7}–e_{1}e_{12}e_{9}+e_{10}e_{3}e_{9})\phi}{(e_{12}e_{2}–e_{11}e_{3})e_{8}},-\frac{(e_{1}e_{12}–e_{10}e_{3})\phi}{e_{12}e_{2}–e_{11}e_{3}},-\frac{(e_{1}e_{11}–e_{10}e_{2})\phi}{-e_{12}e_{2}+e_{11}e_{3}}\right)^{t},$
$\displaystyle\text{ and }\chi_{2}=$
$\displaystyle\left(\frac{-(e_{11}e_{6}e_{7}+e_{12}e_{4}e_{9}–e_{10}e_{6}e_{9})\phi_{1}}{-e_{12}e_{2}e_{7}+e_{11}e_{3}e_{7}+e_{1}e_{12}e_{9}–e_{10}e_{3}e_{9}},\phi_{1},\frac{-(e_{12}e_{2}e_{4}–e_{11}e_{3}e_{4}+e_{1}e_{11}e_{6}–e_{10}e_{2}e_{6})\phi_{1}}{e_{12}e_{2}e_{7}–e_{11}e_{3}e_{7}–e_{1}e_{12}e_{9}+e_{10}e_{3}e_{9}},\frac{-(e_{2}e_{6}e_{7}+e_{3}e_{4}e_{9}–e_{1}e_{6}e_{9})\phi_{1}}{e_{12}e_{2}e_{7}–e_{11}e_{3}e_{7}–e_{1}e_{12}e_{9}+e_{10}e_{3}e_{9}}\right)^{t},$
$\displaystyle\text{where }\phi,\phi_{1}$ $\displaystyle\text{ are arbitrary
real numbers}.$

According to Sotomayor theorem, saddle node bifurcation would occur when
$\chi_{2}^{T}[B_{b}(E_{4},b_{s})]\neq 0,\text{ and
}\chi_{2}^{T}[D^{2}B_{b}(E_{4},b_{s})(\chi_{1},\chi_{1})]\neq 0.$ This can be
proven in a similar manner as done in Theorem 8. This is verified numerically
in section 6 for figure 4.

##  6 Numerical simulations

In this section, the aim is to validate the analytical findings, and also to
explicate the dynamics of the equations of bio-system using ode45, Matcont
package present in Matlab software [32]. For the illustrations, biologically
feasible, and hypothetical values of the parameters are considered as given in
tables 1 and 2. The values may vary within a suitable range.

Biological interpretation | Parameter | Value
---|---|---
Predation rate by juvenile predators | $a_{1}$ | 0.46
Predation rate by matured predators | $a_{2}$ | 0.625
Transformation rate of consumed prey into infant predator | $u$ | 0.8
Transition rate of infant predators into juvenile stage | $b$ | 0.112
Transition rate of juvenile predators into adult stage | $c$ | 0.09
Rate of maturation delay in juvenile predators | $a_{3}$ | 0.06
Natural death rate of infant predators | $d_{1}$ | 0.15
Natural death rate of juvenile predators | $d_{2}$ | 0.1
Natural death rate of adult predators | $d_{3}$ | 0.05
Table 1: Parametric values. Parameter | Value
---|---
$a_{1}$ | 0.6
$a_{2}$ | 0.8
$u$ | 0.82
$b$ | 0.031
$c$ | 0.035
$a_{3}$ | 0.075
$d_{1}$ | 0.026
$d_{2}$ | 0.023
$d_{3}$ | 0.013
Table 2: Parametric values.

For the data set as given in table 1 except for $u$, which is taken $u=0.7$,
it is observed that the equations of bio-system admits a locally
asymptotically stable prey-only equilibrium $E_{2}(1,0,0,0)$. Analytically, as
given in the proof of Theorem 2, either one of the three cases requisite for
the stability of $E_{2}(1,0,0,0)$ should be satisfied. And since we have
$(a_{3}=0.06)<(c=0.09)$ and
$(a_{2}=0.625)<(\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{bu(a_{3}-c)}=0.6336)$,
the equilibrium point is locally stable, which is illustrated in figure 1.
Also, the trajectories of widely varying co-existing initial populations can
be seen leading to the prey-only equilibrium, even if not all of them satisfy
the global stability condition 3, because of its sufficient nature.

For the parametric values from table 1, both of the compresent equilibrium
points exist, but only one of them is espied to be stable. For the equilibrium
point $E_{3}(0.97,0.04,0.03,0.02)$, from Theorem 4 we have
$\epsilon_{1}=1.41,\,\epsilon_{4}=-0.000019,\,\epsilon_{1}\epsilon_{2}-\epsilon_{3}=0.6329,\,\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0.0335$;
and for equilibrium point $E_{4}(0.762,0.30159,0.23419,0.207)$,
$\epsilon_{1}=1.219,\,\epsilon_{4}=0.000126,\,\epsilon_{1}\epsilon_{2}-\epsilon_{3}=0.44539,\,\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0.02096$.
Hence, the system is locally asymptotically stable only around the equilibrium
point $E_{4}$, due to satisfying the Routh-Hurwitz criteria. Figure 2
illustrates this scenario, also at the same time, the widely varying initial
populations approaching equilibrium $E_{4}$ depicts the global stability of
the system.

To investigate the bifurcation behavior of the system, we begin by varying the
transition rate of infant predators into the juvenile stage $b$ (the rest of
the parametric values are taken from table 1). Graphically, a saddle node
bifurcation (also known as fold bifurcation) is espied at $b=0.108186$, the
point where a stable compresent equilibrium $E_{4}$ and an unstable compresent
equilibrium $E_{3}$ encounter and annihilate each other. In addition,
transcritical bifurcation is espied at $b=0.114706$, where the unstable
compresent equilibrium encounters the prey-only equilibrium. Figure 4
illustrates all this, where the red line depicts unstable equilibrium and
green line depicts stable equilibrium. Numerically, for $b=0.114706$, we get
$\frac{(b+d_{1})(a_{3}-c-d_{2})d_{3}}{b(a_{3}-c)u}=0.625\;(=a_{2})$, which is
the condition of transcritical bifurcation of Theorem 8, and for $b=0.108186$,
we have
$a_{2}bux(a_{3}x(-1+2x)-c(-1+2x+a_{1}y_{2})+a_{1}d_{3}y_{3})+(-b-d_{1})[a_{2}a_{3}d_{2}xy_{2}-d_{3}(a_{1}a_{3}xy_{2}+(c+d_{2}-a_{3}x)(-1+2x+a_{1}y_{2}+a_{2}y_{3}))]\approx
0$, which is the condition of saddle node bifurcation of Theorem 9 along with
$\chi_{2}^{T}[B_{b}(E_{4},b_{s})]=0.02745\neq 0$ and
$\chi_{2}^{T}[D^{2}B_{b}(E_{4},b_{s})(\chi_{1},\chi_{2})]=-0.00421\neq 0$. On
that account, it can be said that the theoretical part agrees with the
numerical simulations.

Figure 5 depicts the way the equations of bio-system (2.4) behaves in the
vicinity of the saddle node bifurcation point. Varying the parametric values
of the transition rate of infant predators, $b$ in the neighborhood of
$0.108186$, it is seen that for $b<0.108186$, equilibrium $E_{2}$ is stable,
but the closer it is to the bifurcation point, the more time would the
compresent initial population take to reach the prey-only equilibrium point.
Contrariwise, for $b>0.108186$, the compresent equilibrium is espied to be
stable.

Figure 4: Equilibrium curve in context of $b$, the transition rate of infant
predators into the juvenile stage. Here, $LP$ denotes the saddle node
bifurcation point, and $BP$ denotes the transcritical bifurcation point.

(a)

(b)

(c)

(d)

Figure 5: Equilibrium states of the system 2.4 at different parametric values
of $b$ in the neighbourhood of its saddle node bifurcation point in context of
(a)prey population, (b)infant predator population, (c)juvenile predator
population, (d)adult predator population.

(a)

(b)

(c)

Figure 6: Depicting bi-stability of the equations of bio-system 2.4 for the
same parametric values with different initial population values, i.e., (a) at
initial value $(0.2,0.1,0.1,0.01)$, (b) at initial value
$(0.2,0.1,0.01,0.01)$, (c) at initial value $(0.01,0.3,0.2,0.3)$.

(a)

(b)

Figure 7: Bi-parametric bifurcations: (a) transition rate of juvenile predator
i.e., $c$, and transition rate of infant predator i.e., $b$ (b) rate of
maturation delay in juvenile predators i.e, $a_{3}$ and transition rate of
juvenile predator i.e., $c$. Here, $CP$ denotes cusp point

Bi-stability is a fascinating circumstance whereby there is an existence of
multiple attractors for the same set of parametric values, i.e., the solution
trajectories converge to different attractors, which is extensively
conditional to the initial populations. In our system (2.4), there is a bi-
stability between the compresent equilibrium point $E_{4}$ and the prey-only
equilibrium point $E_{2}$. As can be evidenced from figure 4, for
$b<0.114706$, the prey-only equilibrium point is stable, and for $b>0.108186$,
the compresent equilibrium is stable; which is to say that for
$b\in(0.108186,0.114706)$, both the equilibrium points $E_{2}$ and $E_{4}$ are
stable. Using the set of parametric values from table 1 (i.e., for
$b=0.112\;\;\in(0.108186,0.114706)$ ) at different initial populations, figure
6 is visualized. In figure 6(a), the initial population of the predator
species is sparse, the time evolution of the prey population can be espied to
rapidly increase towards the environment’s carrying capacity and then
gradually decrease to the optimum population size for the co-existence of both
species, while at the same time, all the stages of the predator population
steadily increase to their respective optimum population sizes. But when the
dimensions of the population of predator species is incommensurate, as is in
figure 6(b), then unlike the previous scenario, the predator species are not
in the position to revive again, and thus the equations of bio-system settles
into the prey-only equilibrium point. Furthermore, when the initial biomass of
prey population is low, as is in figure 6(c), the total biomass of the
predator species decreases initially due to scarcity of food resources, giving
the prey an opportunity to increase its biomass, after which the predator
population size too increases, and thereby settling into the co-existing
optimum population sizes, i.e., $E_{4}(0.762,0.30159,0.23419,0.207)$. This can
be espied in the magnified window of the figure 6(c). Interestingly, this same
scenario, identified as ‘bloom’ phenomenon [35, 36], is witnessed even when
the initial prey population is taken incredibly low (for example the initial
prey population is $0.1*10^{-9}$). Therefore, it can be concluded that the
prey-only equilibrium is attainable in this region of bi-stability only when
the total predator population size is sufficiently low.
At this particular set of values as given in table 1, the similar scenario as
witnessed in figures 4,5,6 for the transition rate of infant predators, also
occurs for transition rate of juvenile predators($c$), transformation
rate($u$) and, predation rate of adult predators($a_{2}$) , and the scenario
corresponding to the rate of maturity is also akin to this scenario, albeit in
the reverse orientation, i.e., increasing the parametric value of rate of
maturation delay of juvenile predators ($a_{3}$), spawns the stability of
equilibrium $E_{2}(1,0,0,0)$. While the variations in the parametric values of
predation rate by juvenile predator ($a_{1}$) causes no changes in the system
of equations.

So as to investigate the joint impact of the transition rate of infant
predator $b$, and the transition rate of juvenile predator $c$, the bi-
parametric bifurcation illustration is given in figure 7(a). The entire bi-
parametric region is segregated into regions of stability and bi-stability for
different equilibrium points by the curves of transcritical bifurcation (solid
violet colored lines) and saddle node bifurcation (solid blue colored lines).
Region $I$ and region $III$ are the areas of mono-stability for the prey-only
equilibrium $E_{2}$ and compresent equilibrium $E_{4}$ respectively, while
region $II$ is the area of bi-stability for both the equilibrium points
$E_{2}$ and $E_{4}$. For the parametric values from the left and lower side of
transcritical bifurcation curve, equilibrium point $E_{2}(1,0,0,0)$ is stable
else it is unstable. Similarly, for the values from only the right and upper
side of the saddle node bifurcation curve, the equilibrium point $E_{4}$ is
stable, region $II$ which is common to both the cases becomes a bi-stable
area. The same bi-stable scenario should have been for region $IV$, but the
presence of cusp bifurcation point at $(c,b)=(0.1022,0.07622)$ re-configures
that. The point where transcritical bifurcation curve meets saddle node
bifurcation curve is known as cusp bifurcation point. Because of cusp
bifurcation, there is an existence of a phenomenon called ‘hysterics’.
‘Hysterics’, spawned because of the disappearance of the traced stable
compresent equilibrium $E_{4}$ via the saddle node bifurcation, is a
catastrophic ‘jump’ to a different stable equilibrium, which in our case is
the prey-only equilibrium $E_{2}$. Thus, in region $IV$ only $E_{2}$ is
stable.
In similar manner, figure 7(b) illustrates the dynamical relation between the
parameters, rate of maturation delay $a_{3}$ and the transition rate $c$ of
the juvenile predator. The cusp bifurcation point is at
$(a_{3},c)=(0.04002,0.7054)$. (Stability regions are same as that of the
previous figure 7(a)).

The parametric values as given in table 2 is considered, except $c$, the
transition rate of juvenile predators into adult stage, which is varied, to
investigate the occurrence of Hopf bifurcation (also known as Andronov-Hopf
bifurcation) point. Analytically, at $c=0.03598345(=c_{h})$, we have for the
equilibrium point $E_{3}(0.452951,0.609794,0.755778,0.116978)$,
$\epsilon_{1}=0.547,\,\epsilon_{4}=-0.0001,\,\epsilon_{1}\epsilon_{2}-\epsilon_{3}=0.02,\,\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0.00007$;
and for the equilibrium point $E_{4}(0.063266,0.649761,0.37137,0.89239)$, we
have
$\epsilon_{1}=0.187,\,\epsilon_{4}=0.00003,\,\epsilon_{1}\epsilon_{2}-\epsilon_{3}=0.0013,\,\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}=0$,
therefore the conditions of Hopf Bifurcation along with the transversality
condition
($\frac{\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{3}^{2}-\epsilon_{1}^{2}\epsilon_{4}}{a_{3}}=-0.000025\neq
0$) is satisfied for $E_{4}$. Therefore, there is an emergence of limit cycle
around the equilibrium point $E_{4}$, at $c_{h}$, according to Theorem 6. Next
for the examination of direction and stability of the limit cycle, as
specified in Theorem 7, we have,
$g_{11}=-0.018595+0.0328119i,\;g_{02}=0.0751727+0.397665i,\;g_{20}=-0.112363-0.332041i,\;g_{21}=-0.00766054+0.438565i$.
Hence we get, $C_{1}(0)=-0.0214226+0.134726i$, $\gamma_{2}=0.0580115$, and
$\beta_{2}=-0.0428451$. Here, $\gamma_{2}>0$ implies a supercritical Hopf
bifurcation, and $\beta_{2}<0$ implies the limit cycle being stable. Figure
3(b) depicts the emergence of limit cycles at $c_{h}$. For the values of $c$
greater than $c_{h}$, the system is unstable, and the Hopf bifurcation being a
supercritical one, stable oscillation is witnessed around the equilibrium
$E_{4}$ (figure 3(a)), the greater the distance of the value of $c$ from
$c_{h}$ the greater would be the amplitude of the oscillation. On the
contrary, for $c<c_{h}$, the equations of bio-system evidences a stable focus
point $E_{4}$ (figure 3(c)).

The Hopf bifurcation plot with respect to $a_{3}$, the rate of maturity delay
in juvenile predators is illustrated in figure 8, the bifurcating parametric
value is $a_{3}=0.0604877(=a_{3h})$, when the values of the other parameters
are taken from table 2. The first Lyapunov coefficient at $a_{3h}$ is found to
be $-1.461335e^{-02}$, its negativity signifies the supercritical nature of
Hopf bifurcation. For $a_{3}>a_{3h}$ (figure 8(c)), the trajectory with the
initial co-habitance population spirals into the compresent equilibrium point
$E_{4}$, while for $a_{3}<a_{3h}$, it spirals around the formerly stable
equilibrium point $E_{4}$ (figure 8(a)). The stable limit cycle being born
with its focus at the unstable compresent equilibrium point $E_{4}$ is seen in
figure 8(b)

Next, the bifurcations of the system with respect to a vital parameter, $u$,
the transformation rate of consumed prey into infant predators, is examined in
figure 9 (other parametric values are as is in table 2). At $u=0.281804$, the
system evinces a saddle node bifurcation where a stable compresent equilibrium
$E_{4}$ and an unstable compresent equilibrium $E_{3}$ collide, while at
$u=0.833189$ the system experiences the Hopf bifurcation, which is
supercritical as the first Lyapunov coefficient is $-1.502700e^{-02}$.
Therefore, for $0.281804<u<0.833189$, the equations of bio-system (2.4)
harbours one stable compresent equilibrium $E_{4}$, and when $u$ crosses the
threshold value at $0.833189$, it exchanges its stability with stable limit
cycles around the now unstable equilibrium $E_{4}$. An augmentation of the
modulations in the populations of infant predators, juvenile predators, adult
predators and prey can be seen with further increment of the value of $u$.

(a) at $a_{3}=0.055$

(b) at $a_{3}=0.0604877$

(c) at $a_{3}=0.084$

Figure 8: Phase portrait in a bi-dimensional space depicting predator-prey
population dynamics in the neighborhood of Hopf bifurcation point of the rate
of maturation delay in juvenile predators, $a_{3}$: $(a)a_{3}<a_{3h}$,
$(b)a_{3}=a_{3h}$, $(c)a_{3}>a_{3h}$, $a_{3h}$ is the Hopf bifurcation point.

(a)

(b)

(c)

(d)

Figure 9: Equilibrium states of the prey, infant predator, juvenile predator
and, adult predator populations with respect to the transformation rate of
consumed prey into infant predator, $u$ for (a) prey population, (b) infant
predator population, (c) juvenile predator population, (d) adult predator
population. $LP$ denotes the saddle node bifurcation point and, $H$ denotes
the Hopf bifurcation point.

(a)

(b)

Figure 10: Bi-parametric bifurcations: (a) predation rate of juveniles,
$a_{1}$, and transition rate of juvenile predator, $c$ (b) predation rate of
juveniles, $a_{1}$, and predation rate of adult predators, $a_{2}$.

(a)

(b)

(c)

(d)

Figure 11: Equilibrium states of (a) prey population, (b) infant predator
population, (c) juvenile predator population, (d) adult predator population
depicting bifurcations with respect to the transition rate of infant predator
into juvenile stage, $b$.

In figure 9, the red line denotes unstability, the green line denotes
stability, and the blue line denotes the amplitude of fluctuation at that
particular parametric value of $u$. An illustration alike figure 9 can be
espied for the parameter $a_{2}$, the predation rate of adult predators. While
varying $a_{2}$, the saddle node bifurcation is experienced by the system at
$a_{2}=0.274931$, the Hopf bifurcation point is detected at $a_{2}=0.810103$.

Figure 10(a) illustrates the change in dynamical behavior of the equations of
bio-system (2.4) with variations in the predation rate of juvenile predators
($a_{1}$) and the transition rate of juvenile predators ($c$), by continuing
the curves of saddle node bifurcation and Hopf bifurcation. These two curves
partition the whole bi-parametric region into sub-regions. The Hopf
bifurcation curve separates the unstable compresent equilibrium region
($region\,III$) from the stable one. The Hopf curve being a supercritical one,
there exists stable limit cycle in the neighborhood of the curve. The region
$IV$ bounded between the Hopf curve and the saddle node curve is the area of
stability of compresent equilibrium point $E_{4}$, after crossing the saddle
node curve and stepping into region $I$, the stability of $E_{4}$ vanishes
because of collision with the unstable equilibrium $E_{3}$. Furthermore, the
two bifurcation curves meet tangentially at $(a_{1},c)=(0.149588,0.018589)$,
called Bogdanov-Takens ($BT$) point. Because of $BT$ point, the non-saddle
equilibrium $E_{4}$ which undergoes Hopf bifurcation, loses the occurrence of
limit cycle via the saddle node bifurcation, in addition to that $BT$ point
causes the separation of the saddle node curve into two parts: stable curve
and unstable curve. Hence, in region $II$, neither the occurrence of limit
cycle nor the stability of $E_{4}$ is found. Similarly, in figure 10(b), the
Hopf curve and saddle node curve meet tangentially at Bogdanov-Takens point,
$(a_{1},a_{2})=(0.019133,0.274931)$, thereby making region $II$ synonymous
with region $II$ of figure 10(a). Also, in figure 10(b), there exists a
generalized Hopf bifurcation point (also known as Bautin bifurcation point) at
$(a_{1},a_{2})=(0.042954,0.379816)$, where the first Lyapunov coefficient
vanishes. A division of the Hopf curve into supercritical part and subcritical
part, is done by Bautin bifurcation point. For higher values of $a_{1}$ and
$a_{2}$ (i.e., at the right hand side of Bautin bifurcation point) the Hopf
curve is supercritical. Colossally, like the previous scenario, region $III$
is the area of unstability of compresent equilibrium, and regions $I$ and $IV$
are the regions of stability of compresent equilibrium $E_{4}$, while $E_{4}$
loses its stability via saddle node bifurcation curve in region $V$. The
horizontal line corresponding to saddle node bifurcation curve in figures
signifies that the values of both transition rates of juvenile predator ($c$)
or predation rate of adult predator ($a_{2}$) remain the same in
correspondence to all values of $a_{1}$, the predation rate of juvenile
predators, i.e., $a_{1}$ is not responsible for the occurrence of saddle node
bifurcation. However, it should be noted that for a completely different set
of values, $a_{1}$ may cause the said bifurcation, by reason of the presence
of the parameter in the condition given in theorem 9 for saddle node
bifurcation.

The local bifurcations while varying the transition rate of infant predators,
$b$ are visualized in figure 11. The equations of bio-system changes its
stability through Hopf bifurcation not once but twice, the bifurcating points
being $b=0.032488$ and $b=0.136940$ whose Lyapunov coefficients are
$-1.619062e^{-02}$ and $-7.177191e^{-02}$ respectively. Hence, both the points
are supercritical Hopf bifurcation points. For $b\in(0.032488,0.136940)$,
there is the existence of a stable limit cycle whose amplitude increases and
then gradually decreases as the value of $b$ travels from one bifurcation
point to another, forming a ‘balloon’ like structure as can be seen in the
figure 11. For $b>0.136940$ and $b_{s1}<b<0.032488$, compresent equilibrium
$E_{4}$ is stable, where $b_{s1}=0.005977$ is the saddle node bifurcation
point. The solid red and green lines denote unstable and stable equilibriums
respectively, while the blue lines depict the amplitudes.

##  7 Discussion and conclusion

Juveniles of the predator have been chronicled to hunt their prey. But
predation without any prior experience, i.e., without knowing the prey
specific approach would be a hellacious venture. Concomitantly, the prey oft-
times counterattacks as a desperate remedy for intercepting their demise,
thereby effectuating injuries, which are induced particularly in the
inexperienced predator. Also, as evidenced from all the documentations, there
is always a period of time in the beginning of the predator’s life where they
are completely reliant on their parents for survival. Thereupon, while
considering predation by immature predators, it would be perspicuous for the
predator species to be structured into three stages (infant, juvenile, and
adult), while simultaneously incorporating injury, as the baleful consequence
of predation by juveniles, ensuing a delay in the maturation rate for the
juvenile predator into the adult stage. In this paper, equations of bio-system
representing this very scenario has been assembled, thus providing a much more
detailed and richer dynamic of the prey-predator interaction.

The boundedness of the equations of bio-system (2.4) is authenticated, proving
its feasibility in the real-life scenario, cause none of the species should be
able to grow exponentially owing to the constraints of natural resources. The
system of equations (2.4) harbours four biologically viable equilibrium
points: one each of extinction equilibrium and prey-only equilibrium, and two
of compresent equilibrium points. With the intent of studying the behavior of
the system in the proximity of the equilibrium points, the conditions for
local and global stability are examined for the equilibrium points in Theorems
2, 3, 4 and, 5, which are then authenticated by visualizing graphically in
figures 1, and 2.

An interesting state of affairs known as ‘bi-stability’ (a situation with
multiple attractors for the system of equations) is observed between
compresent equilibrium and prey-only equilibrium points in figures 4, 6. The
latter equilibrium point is obtained only with a considerably low initial
population of one or more stages of the predator species. It is noteworthy
that even with negligible initial prey population, the system of equations
instead of settling into extinction equilibrium, goes towards compresent
equilibrium, which figure 6(c) shows happening, by virtue of rapidly
decreasing predator population (ascribable to inanition), in parallelism the
prey population rapidly increases, before both of the species converges to the
compresent equilibrium point. This circumstance is known as ‘bloom’
phenomenon.

Analysis of bifurcations is a prerequisite for uncovering the rich dynamics of
the equations of bio-system (2.4), because it enables the determination of the
transition points while tracing the paths of steady as well as unstable
equilibriums. The equations of bio-system exhibits the occurrence of Hopf
bifurcation around compresent equilibrium, i.e., the transition point where
the said equilibrium point exchanges its stability with either a stable or
unstable limit cycle (Theorem 6). The analytical methodology used to determine
the stability and direction of the thus occurred limit cycle is shown Theorem
7 and then verified numerically through figure 3. Likewise, the manifestation
of saddle node bifurcation is corroborated both analytically (Theorem 9) and
graphically (figure 4). The expression
$a_{2}bu(a_{3}-c)=(b+d_{1})(a_{3}-c-d_{2})d_{3}$ denotes the surface on which
transcritical bifurcation occurs (Theorem 8), the bifurcation responsible for
the interchangeability of stability between the equilibrium points. The
transcritical bifurcation is found to be existent between the prey-only
equilibrium and one of the compresent equilibrium. As can be witnessed from
the expression above, all the parameters are responsible for causing this
bifurcation except for $a_{1}$, the predation rate of juveniles.

When the influence of transition rate of juvenile predators into the adult
stage ($c$) is investigated, an exchange of stability from co-existing
population to the stability of prey-only population is witnessed after
decreasing the value of $c$ after a threshold point because of transcritical
bifurcation for a certain parametric set of values (table 1). While for
another set of values for the parameters i.e., table 2, the compresent
equilibrium $E_{4}$ is stable, but increasing the parametric value of $c$
after the threshold value of Hopf bifurcation point, each and every stage of
the predator population and the prey population, incapable of abiding by a
specific population size due to the unstability of $E_{4}$, keeps on
oscillating (figure 3). A similar sequence of events can be espied for the
parameter $a_{3}$, the rate of maturity delay of juvenile predators, but in
the reverse orientation, for example, a decrease in parametric value after the
Hopf bifurcation point leads to the birth of the limit cycle around the
compresent equilibrium (figure 8). Also, the existence of the fold bifurcation
point where the unstable compresent equilibrium $E_{3}$ and stable compresent
equilibrium $E_{4}$ collide, can be found for both $c$ and $a_{3}$. Both of
these parameters combined together as $(c-a_{3}x)$ gives the comprehensive
maturation rate of juvenile predators into the adult stage. Hence,
biologically we can surmise that when this rate increases, the juvenile
population decreases leading to lesser juvenile predator maturing in the next
instance, correspondingly the adult predator would increase at first but their
number too would gradually decrease both due to lesser recruitment and lesser
availability of prey (decreased number of prey is attributable to more of them
being killed because of increased adult predator population). Then when the
biomass of the predator population would be less, the prey species would
increase their number, resulting in the predator population following suit
because of the now availability of food, and thus creating the oscillations.
On the contrary, with low comprehensive maturation rate, the predator
population not being able to survive goes towards extinction, i.e. the
equilibrium point $E_{2}(1,0,0,0)$.

Predation by juveniles can be appraised as ancillary, owing to the fact that
it doesn’t equate to proliferation of the species, but rather to amplification
of expertise of the juveniles while hunting for their later stage. Varying
this parameter can cause Hopf-bifurcation in the system, since in figures
7(a), 7(b), for different values of $c$ and $a_{2}$, the Hopf bifurcating
point of $a_{1}$ too differs (when the values of other parameters are taken
from table 2). For values lesser than the threshold value of $a_{1}$, the
equations of bio-system has unstable compresent equilibrium. That is for
lesser value of $a_{1}$, more prey survives, which corresponds to a decrease
in comprehensive maturation rate of juvenile predators ($c-a_{3}x$), the
biological implication of which can be comprehended as the juveniles are not
able to kill the prey with the reason of the anti-predator behavior of the
prey being too strong, which also results into more juveniles getting injured.
Also, since $a_{1}$ doesn’t possess the ascendancy in causing transcritical
bifurcation, the prey-only equilibrium remains unstable.

From a new vantage point, a stage in the predator species is introduced,
identified in the present study as infant predators, which is similar to that
of the immature predators as is considered in the conventional prey- predator
model with stage structure in predator [35, 36]. Their transition rate into
the juvenile stage is extensively investigated through graphical
visualization. A decrease in $b$ is espied to cause the extinction of the
predator species, that is the attainment of prey-only equilibrium (figure 4).
Biologically, this can be perceived as- with lesser predator individuals
attaining the juvenile stage leads to a lesser number of adult predators,
which corresponds to frivolous birth of infants, thereby the eventual
extinction. Also, as is espied in figure 11, an increment of this transition
rate beyond a threshold value results in destabilisation of compresent
equilibrium, similar to the scenario of the birth of oscillations in all the
populations with the increment of transition rate of juvenile predators, the
solution trajectory keeps on encircling the equilibrium point beyond a
threshold parametric value of the transition rate of infant predators. But the
circumference of the solution trajectory at first increases with the increment
of the parametric value, and then again decreases due to the manifestation of
another supercritical Hopf bifurcation, ensuing the stability of compresent
equilibrium point all over again. Therefore, there exists a neighbourhood in
the parametric value of the infants’ transition rate where none of the
populations can conform to its optimal number, that is, after crossing the
second threshold parametric value, both of the species can obtain their
respective optimal populations

The existence of the predators obviously relies a lot on the predation rate of
adult predators ($a_{2}$), and the transformation of consumed prey into infant
predators ($u$), which can be clearly seen with the existence of supercritical
Hopf bifurcation and saddle node bifurcation with respect to these two
parameters. Furthermore, diminishing either of them beyond a threshold value
insinuates non-viability of the subsistence of the predator species;
contrarily, an increment beyond a certain value interdicts the existence of
both the species.

Occurrence of bifurcation of equilibria in a two-parameter region such as cusp
bifurcation point, Bogdanov Takens point, and Bautin bifurcation point implies
the richer dynamics of the equations of bio-system. The ascendancy of these
bifurcation points over the system is shown in figures 7, and 10 and explained
in section 6.

With no biological data being accessible, the parametric values are based on
hypothetical data, but the analysis is hoped to be useful for researchers
doing their work on related fields with experimental data. For further
analysis, other instances such as refuge, fear effect, group hunting, and
others can be incorporated in the current equations of bio-system (2.4).

Competing interests:  The authors declare no competing interests.

## References

  * [1] Lotka, A. J. (1910). Contribution to the theory of periodic reactions. The Journal of Physical Chemistry, 14(3), 271-274.
  * [2] Volterra, V. (1927). Fluctuations in the abundance of a species considered mathematically. Nature, 119(2983), 12-13.
  * [3] Gurney, W. (1986). The systematic formulation of models of stage-structured populations.
  * [4] Nisbet, R. M., & Gurney, W. S. C. (1984). “Stage-structure” models of uniform larval competition. In Mathematical ecology (pp. 97-113). Springer, Berlin, Heidelberg.
  * [5] Nisbet, R. M., Gurney, W. S. C., & Metz, J. A. J. (1989). Stage structure models applied in evolutionary ecology. In Applied Mathematical Ecology (pp. 428-449). Springer, Berlin, Heidelberg.
  * [6] Aiello, W. G., & Freedman, H. I. (1990). A time-delay model of single-species growth with stage structure. Mathematical biosciences, 101(2), 139-153.
  * [7] Aiello, W. G., Freedman, H. I., & Wu, J. (1992). Analysis of a model representing stage-structured population growth with state-dependent time delay. SIAM Journal on Applied Mathematics, 52(3), 855-869.
  * [8] Elbroch, L. M., Feltner, J., & Quigley, H. B. (2017). Stage‐dependent puma predation on dangerous prey. Journal of Zoology, 302(3), 164-170.
  * [9] Farrell, T. M., Smiley-Walters, S. A., & McColl, D. E. (2018). Prey species influences foraging behaviors: Rattlesnake (Sistrurus miliarius) predation on Little Brown Skinks (Scincella lateralis) and Giant Centipedes (Scolopendra viridis). Journal of Herpetology, 52(2), 156-161.
  * [10] Figueiredo, A., Alves-Martins, N., & Nogueira-Costa, P. (2021). Predation attempt by the Spectacled Caiman, Caiman crocodilus (Linnaeus, 1758), on the microhylid Elachistocleis carvalhoi Caramaschi, 2010 in the southeastern Amazon of Brazil. Herpetology Notes, 14, 1227-1229.
  * [11] Gajski, D., Petráková, L., & Pekár, S. (2020). Ant-eating spider maintains specialist diet throughout its ontogeny. Journal of Zoology, 311(3), 155-163.
  * [12] Gosselin, L. A., & Chia, F. S. (1996). Prey selection by inexperienced predators: do early juvenile snails maximize net energy gains on their first attack?. Journal of experimental marine biology and ecology, 199(1), 45-58.
  * [13] Gibbons, M. J. (1988). Impact of predation by juvenile Clinus superciliosus on phytal meiofauna: are fish important as predators?. Marine ecology progress series, 45, 13-22.
  * [14] Zuluaga, S., Vargas, F. H., Aráoz, R., & Grande, J. M. (2022). Main aerial top predator of the Andean Montane Forest copes with fragmentation, but may be paying a high cost. Global Ecology and Conservation, 37, e02174.
  * [15] Dhellemmes, F., Smukall, M. J., Guttridge, T. L., Krause, J., & Hussey, N. E. (2021). Predator abundance drives the association between exploratory personality and foraging habitat risk in a wild marine meso‐predator. Functional Ecology, 35(9), 1972-1984.
  * [16] Huang, R. K., Webber, Q. M., Laforge, M. P., Robitaille, A. L., Bonar, M., Balluffi-Fry, J., Zabihi-Seissan, S., & Vander Wal, E. (2021). Coyote (Canis latrans) diet and spatial co-occurrence with woodland caribou (Rangifer tarandus caribou). Canadian Journal of Zoology, 99(5), 391-399.
  * [17] MacNulty, D. R. (2002). The predatory sequence and the influence of injury risk on hunting behavior in the wolf (Doctoral dissertation, University of Minnesota).
  * [18] Brown, J. S., Embar, K., Hancock, E., & Kotler, B. P. (2016). Predators risk injury too: the evolution of derring-do in a predator–prey foraging game. Israel Journal of Ecology and Evolution, 62(3-4), 196-204.
  * [19] Miranda, E. B., Menezes, J. F. D., & Rheingantz, M. L. (2016). Reptiles as principal prey? Adaptations for durophagy and prey selection by jaguar (Panthera onca). Journal of Natural History, 50(31-32), 2021-2035.
  * [20] Pokharel, P., Sippel, M., Vilcinskas, A., & Petschenka, G. (2020). Defense of milkweed bugs (Heteroptera: Lygaeinae) against predatory lacewing larvae depends on structural differences of sequestered cardenolides. Insects, 11(8), 485.
  * [21] Mukherjee, S., & Heithaus, M. R. (2013). Dangerous prey and daring predators: a review. Biological reviews, 88(3), 550-563.
  * [22] Fernández Moya, S., Iglesias Pastrana, C., Marín Navas, C., Ruíz Aguilera, M. J., Delgado Bermejo, J. V., & Navas González, F. J. (2021). The Winner Takes it All: Risk Factors and Bayesian Modelling of the Probability of Success in Escaping from Big Cat Predation. Animals, 12(1), 51.
  * [23] Escalante, I. (2015). Predatory behaviour is plastic according to prey difficulty in naïve spiderlings. Journal of insect behavior, 28(6), 635-650.
  * [24] García, L. F., Franco, V., Robledo-Ospina, L. E., Viera, C., Lacava, M., & Willemart, R. H. (2016). The predation strategy of the recluse spider Loxosceles rufipes (Lucas, 1834) against four prey species. Journal of Insect Behavior, 29(5), 515-526.
  * [25] Perlman, Y., & Tsurim, I. (2008). Daring, risk assessment and body condition interactions in steppe buzzards Buteo buteo vulpinus. Journal of Avian Biology, 39(2), 226-228.
  * [26] Kaushik, R., & Banerjee, S. (2022). Predator–prey system with multiple delays: prey’s countermeasures against juvenile predators in the predator–prey conflict. Journal of Applied Mathematics and Computing, 68(4), 2235-2265.
  * [27] Kaushik, R., & Banerjee, S. (2021). Predator-prey system: Prey’s counter-attack on juvenile predators shows opposite side of the same ecological coin. Applied Mathematics and Computation, 388, 125530.
  * [28] Li, J., Liu, X., & Wei, C. (2022). Stationary distribution of a stage-structure predator–prey model with prey’s counter-attack and higher-order perturbations. Applied Mathematics Letters, 129, 107921.
  * [29] Bhattacharjee, D., Roy, T., Acharjee, S., & Dutta, T. K. (2022). Stage structured prey-predator model incorporating mortal peril consequential to inefficiency and habitat complexity in juvenile hunting. Heliyon, 8, e11365.
  * [30] Yao, P., Wang, Z., & Wang, L. (2022). Stability Analysis of a Ratio-Dependent Predator-Prey Model. Journal of Mathematics, 2022, 46052.
  * [31] Jawad, S., & Naji, R. K. (2022). The Influence of Stage Structure and Prey Refuge on the Stability of the Predator-Prey Model. I. J. Engineering and Manufacturing, 3, 51-59.
  * [32] Dhooge, A., Govaerts, W., Kuznetsov, Y. A., Meijer, H. G. E., & Sautois, B. (2008). New features of the software MatCont for bifurcation analysis of dynamical systems. Mathematical and Computer Modelling of Dynamical Systems, 14(2), 147-175.
  * [33] Perko, L. (2013). Differential equations and dynamical systems (Vol. 7). Springer Science & Business Media.
  * [34] Hassard, B. D., Hassard, B. D., Kazarinoff, N. D., Wan, Y. H., & Wan, Y. W. (1981). Theory and applications of Hopf bifurcation (Vol. 41). Cambridge, New York, Cambridge University Press.
  * [35] Khajanchi, S. (2017). Modeling the dynamics of stage-structure predator-prey system with Monod–Haldane type response function. Applied Mathematics and Computation, 302, 122-143.
  * [36] Khajanchi, S. (2014). Dynamic behavior of a Beddington–DeAngelis type stage structured predator–prey model. Applied Mathematics and Computation, 244, 344-360.
  * [37] Georgescu, P., & Hsieh, Y. H. (2007). Global dynamics of a predator-prey model with stage structure for the predator. SIAM Journal on Applied Mathematics, 67(5), 1379-1395.
  * [38] Liu, Q., Jiang, D., Hayat, T., & Alsaedi, A. (2018). Dynamics of a stochastic predator–prey model with stage structure for predator and Holling type II functional response. Journal of Nonlinear Science, 28, 1151-1187.