Modeling and verification of an intelligent tutoring system based on Petri net theory

1.   Introduction

The Lotka-Volterra equations that represent the system of continuous-time non-linear differential equations whose solution is deterministic, evolved to determine the dynamics of biological systems. This system is a more generalized version of the Kolmogorov model because it only focuses on the predator-prey interactions but ignores competition, disease and mutualism, which the Kolmogorov model includes. As compared with other mathematical model, following assumptions are made in order to formulate the Lotka-Volterra equations:

$\rm (i)$ For the prey population there is no shortage of food.

$\rm (ii)$ The amount of food supplied to the prey is directly related to the size of the prey population.

$\rm(iii)$ The rate of change of population is directly proportional to its size.

$\rm (iv)$ The environment is constant and genetic adaptation is not assumed to be negligible.

$\rm(v)$ Predators will never stop eating.

Based on said assumptions a continuous-time, two-species Lotka-Volterra system takes the following form ^[1,2,3]:

where $x$ and $y$ respectively denotes the number of prey and predator. The quantities: $\frac{dx}{dt}$ and $\frac{dy}{dt}$ respectively denotes the instantaneous rates of prey and predator, and $a$, $b$, $c$ and $d$ are positive real parameters. In $1^{st}$ equation of system (1.1), $ax$ represent that prey are reproduced exponentially whereas $bxy$ shows that the rate at which the predator kills prey. So, if there is no population of prey or no population of predator, no decrease in the population of prey can occur. The equation for prey can be summed up as the rate at which new prey is born, minus the rate at which prey is killed off. Now in $2^{nd}$ equation of system (1.1), $cxy$ represent the rate at which predators kill prey except that a different constant is used to describe this relationship since the rate at which predators kill and the rate at which they reproduce are not identical. Also, $dy$ is assumed to be an exponential decay, which mean that prey cannot kill the predators and that the decrease in predator population is due to natural death or emigration.

It is anticipated that as compared to continuous-time system, discrete-time system designated by difference equations are widely explored because these models are more reasonable, if populations have non-overlapping generations, as compared to differential type models, and also one can obtain more accurate numerical simulations ^[4,5,6]. During the last few years, discrete-time predator-prey models have been extensively investigated by studying fixed points, local and global dynamical properties, existence of possible bifurcation and chaos ^[7,8,9] and references cited therein. So, in this work, we explore the local dynamics, bifurcations and hybrid control in a $3\times3$ discrete-time predator-prey model, i.e., two prey and one predator, representing by the following system of difference equations ^[10]:

where $x_n$, $y_n$ and $z_n$ respectively denotes population densities of prey and predator respectively, and $a, b, c, d, e, f$ and $r$ are real positive parameters. More precisely, $r$ is the net reproductive rate of the prey population, $a$ denotes the intrinsic rate of growth the prey, $b$ is the rate of conversion of a consumed prey to a predator, $c$ is the per capita rate of predation of the predator, $d$ is death rate of the predator, $e$ is the rate of conversion of a consumed prey to a predator and $f$ is the per-capita rate of predation of the predator.

It is important here to mention that in 2013 Sagayaraj et al. ^[10] studied the local dynamics about three boundary fixed points only: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ and $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, out of five fixed points including the unique positive fixed point. So for the interested readers, these results regarding the local dynamics about fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ and $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ are summarized in Table 1.

From the Table 1, it is clear that for the said fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ and $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, Sagayaraj et al. ^[10] cannot give complete local dynamical classifications to interested readers, which is still further consideration and improvements towards local dynamical properties along with topological classifications and bifurcation analysis for the model under consideration. So, in this paper we will study the local dynamics about fixed points, bifurcations like fold bifurcation, Neimark-Sacker or Hopf bifurcation, supercritical or subcritical Neimark-Sacker bifurcation, period-doubling or flip bifurcation, supercritical or subcritical period-doubling bifurcation, and hybrid control in the model (1.2). More precisely, our key contributions in this paper are as follows:

● To explore the local dynamical properties along with different topological classifications about fixed points of the under consideration model.

● To explore the existence of prime period and periodic points of the model (1.2).

● To explore the existence of possible bifurcations about fixed points of $3\times3$ model.

● Presentation of numerical simulation to verify theoretical results.

● To control the N-S and P-D bifurcations about the positive fixed point by Hybrid control strategy.

This paper is structured as follows: In Section 2, we explored the existence of fixed points of the model (1.2) algebraically. The linearized form of the model (1.2) is presented in Section 3. In Section 4, we explored the local dynamics about boundary fixed points and unique positive fixed point of the model. Existence of prime period and periodic points of the model (1.2) are explored in Section 5. In Section 6, we explored the existence of bifurcations about boundary fixed points of the model (1.2) whereas detailed bifurcation analysis about boundary fixed points and unique positive fixed point of the model (1.2) are given in Section 7. Theoretical results are verified numerically in Section 8. In Section 9, the hybrid control strategy is presented to control the N-S and P-D bifurcations about the unique fixed point of the model while a brief summary of the paper is presented in Section 10.

2.   Existence of fixed points of the model (1.2)

The existence of fixed points in the interior of $\mathbb{R}_+^3$ is explored in this Section. The result regarding the existence of fixed points for the $3\times 3$ discrete-time model (1.2) can be summarized as a following Lemma:

Lemma 2.1. In the interior of $\mathbb{R}_+^3$, model (1.2) has four boundary fixed points and a unique positive fixed point. Precisely

(i) $\forall \ a, b, c, d, e, f, r$, $P_{000}(0, 0, 0)$ is a boundary fixed point of model (1.2);

(ii) $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ is a boundary fixed point of the model (1.2) if $r > 1$;

(iii) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is a boundary fixed point of the model (1.2) if $f > \frac{dr}{r-1}$ where $r > 1$;

(iv) $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is also boundary fixed point of the model (1.2), $\forall \ a, b, c, d > 0$;

(v) $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is the unique positive fixed point of (1.2) if $d > \frac{bf(r-1)-afc}{br}$ and $b > \frac{ac}{r-1}$ where $r > 1$.

Proof. If $P_{xyz}(x, y, z)$ is a fixed point of model (1.2) then

Since for $P_{000}(0, 0, 0), \ P_{0y0}\left(0, \frac{r-1}{r}, 0\right), \ P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ and $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$, Eq (2.1) satisfied identically, and hence $P_{000}(0, 0, 0), \ P_{0y0}\left(0, \frac{r-1}{r}, 0\right), \ P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ and $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ are boundary fixed points of the model (1.2). In order for the existence of the unique positive fixed point one need to solve the following equation simultaneously:

From first equation of system (2.2) one gets:

From second equation of system (2.2) and Eq (2.3), one gets:

Finally, from third equation of system (2.2) and Eq (2.4), one gets:

In view of Eqs (2.3)–(2.5) one can conclude that $P_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is the interior fixed point of model (1.2). Additionally, $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is the unique positive fixed point of the model (1.2) if $d > \frac{bf(r-1)-afc}{br}$ and $b > \frac{ac}{r-1}$ where $r > 1$.

3.   Linearized form of the model (1.2)

The linearized form of $3\times3$ model (1.2) about fixed point $P_{xyz}(x, y, z)$ is explored in this Section. So we have the following map to create the corresponding linearized form of model (1.2):

where

Finally, $J|_{P_{xyz}(x, y, z)}$ about $P_{xyz}(x, y, z)$ under the map, which is depicted in Eq (3.1), is

4.   Local dynamical properties about fixed points of the model (1.2)

Local dynamical properties along with topological classifications about fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of $3\times3$ model (1.2) are explored in this section.

4.1.   Local dynamical properties about $P_{000}(0, 0, 0)$

The $J|_{P_{000}(0, 0, 0)}$ about $P_{000}(0, 0, 0)$ is

The eigenvalues of $J|_{P_{000}(0, 0, 0)}$ about $P_{000}(0, 0, 0)$ are

So by stability theory ^[11], the local dynamics of $3\times3$ model (1.2) about $P_{000}(0, 0, 0)$ can be conclude as follows.

Lemma 4.1. For $P_{000}(0, 0, 0)$ following conditions hold:

(i) $P_{000}(0, 0, 0)$ is never sink;

(ii) $P_{000}(0, 0, 0)$ is a source if

(iii) $P_{000}(0, 0, 0)$ is a saddle if

(iv) $P_{000}(0, 0, 0)$ is non-hyperbolic if

or

4.2.   Local dynamical properties about $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$

The $J|_{P_{0y0}\left(0, \frac{r-1}{r}, 0\right)}$ about $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ is

The eigenvalues of $J|_{P_{0y0}\left(0, \frac{r-1}{r}, 0\right)}$ about $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ are

We also summarized the local dynamics of $3\times3$ model (1.2) about $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ as follows:

Lemma 4.2. For $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, following statements hold:

(i) $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ is never sink;

(ii) $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ is a source if

(iii) $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ is a saddle if

(iv) $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$ is non-hyperbolic if

or

4.3.   Local dynamical properties about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$

The $J|_{P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)}$ about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is

The eigenvalues of $J|_{P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)}$ about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ are

where

In the following Lemmas, we give the topological classifications about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ if $\Delta\ge0\ (\mathrm{resp.}\ \Delta < 0)$ by existing literature ^[9,12].

Lemma 4.3. If $\Delta = \frac{4d^2fr+d^2r^2-4drf^2+4d^2fr}{f^2}\ge0$ then for $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ of model (1.2), following statements hold:

(i) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is a stable node if

(ii) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is an unstable node if

(iii) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is saddle node if

or

(iv) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is non-hyperbolic if

or

Lemma 4.4. If $\Delta = \frac{4d^2fr+d^2r^2-4drf^2+4d^2fr}{f^2} < 0$ then for $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ of model (1.2), following statements hold:

(i) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is stable focus-node if

(ii) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is an unstable focus-node if

(iii) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is saddle focus-node if

or

(iv) $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is non-hyperbolic if

or

4.4.   Local dynamical properties about $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$

The $J|_{P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)}$ about $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is

The eigenvalues of $J|_{P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)}$ about $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ are

Lemma 4.5. For $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$, following statements hold:

(i) $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is never stable focus-node;

(ii) $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is an unstable focus-node if

(iii) $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is a saddle focus-node if

(iv) $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is non-hyperbolic if

4.5.   Local dynamical properties about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$

The $J|_{P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)}$ about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is

The characteristic equation of $J|_{P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)}$ about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is

where

Now in the following lemma, we explore the necessary and sufficient condition under which the roots of (4.34) lies on open unit disk by utilizing Theorem 1.4 of ^[11].

Lemma 4.6. $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of $3\times3$ model (1.2) is a sink iff the following conditions hold:

where $\delta_1, \delta_2$ and $\delta_3$ are depicted in (4.35).

5.   Existence of prime period and periodic points of the model (1.2)

We explore that the fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of the model (1.2) are periodic points of prime period- $1$. Additionally, we also explore that $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of the model (1.2) are periodic points having period- $2, 3, \cdots, n$, by existing Literature ^[13,14,15].

Theorem 5.1. Fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of model (1.2) are periodic points of prime period- $1$.

Proof. From model (1.2), one has

where $f, g$ and $h$ are represented in system (3.2). After straightforward computation, from Eq (5.1) one gets:

From Eq (5.2), one can conclude that $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of model (1.2) is a periodic point of prime period- $1$.

Hereafter we proved that fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of $3\times3$ model (1.2) are periodic points having period- $2, 3, \cdots, n$.

Theorem 5.2. The fixed point: $P_{000}(0, 0, 0)$ of $3\times3$ model (1.2) is a periodic point of period- $2, 3, \cdots, n$.

Proof. From Eq (5.1), one has

Thus Eq (5.3) indicates that $P_{000}(0, 0, 0)$ of model (1.2) is a periodic point of period-$2, 3, \cdots, n$.

Theorem 5.3. The fixed points: $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of $3\times3$ model (1.2) are periodic points of period- $2, 3, \cdots, n$.

Proof. Same as the proof of Theorem 5.2.

6.   Existence of bifurcations about fixed points of the $3\times3$ model (1.2)

Existence of bifurcations about boundary fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ and $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ are explored in this Section. Moreover occurrence of bifurcations about boundary fixed points can be summarized as follows:

$\rm(i)$ About $P_{000}(0, 0, 0)$, one of the eigenvalues $\lambda_2|_{(4.5)} = 1$ but $\lambda_{1, 3} = 1+a, 1-d \ne 1\ \mathrm {or }\ -1$. So $3\times3$ model (1.2) may undergoes a fold bifurcation if parameters $\Lambda = (a, b, c, d, e, f, r)$ goes through the following curve:

For the interested readers, it is important here to mention that we have introduced new notation like $\lambda_2|_{(4.5)}$ which stands that we have evaluated the eigenvalue $\lambda_2$ at non-hyperbolic condition, which is depicted in Eq (4.5). Similarly we use such notation for same purpose in the rest of the paper.

Moreover, about $P_{000}(0, 0, 0)$, one of the eigenvalues $\lambda_3|_{(4.6)} = -1$ but $\lambda_{1, 2} = 1+a, r \ne 1\ \mathrm {or }\ -1$. So $3\times3$ model (1.2) may undergoes a P-D bifurcation if parameters $\Lambda$ goes through the following curve:

It is noted that for $d = 2$ the predator leaves a negative number of offspring, i.e., in the absence of prey and therefore its no make sense biology. Thus we will exclude this case from further study.

$\rm (ii)$ About $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, one of the eigenvalues $\lambda_2|_{(4.11)} = -1$ but $\lambda_{1, 3} = 1+a, 1-d+\frac{f(r-1)}{r} \ne 1\ \mathrm {or}\ -1$. So $3\times3$ model (1.2) undergoes a P-D bifurcation if parameters $\Lambda$ goes through the following curve:

Also $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, one of the eigenvalues $\lambda_3|_{(4.12)} = -1$ but $\lambda_{1, 2} = 1+a, 2-r \ne 1\ \mathrm {or}\ -1$. So $3\times3$ model (1.2) undergoes a P-D bifurcation if $\Lambda$ goes through the following curve:

$\rm(iii)$ About $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, pair of complex conjugate eigenvalues $|\lambda_{2, 3}|_{(4.26)} = 1$ but $\lambda_1 = 1+a+\frac{b}{c}\left(1-r+\frac{dr}{f}\right) \ne 1\ \mathrm {or}\ -1$. So $3\times3$ model (1.2) undergoes a N-S bifurcation if parameters $\Lambda$ goes through the following curve:

$\rm(iv)$ About $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$, $\lambda_1|_{(4.32)} = -1$ but $|\lambda_{2, 3}| = 1+ad > 1$. So $3\times3$ model (1.2) may undergoes P-D bifurcation if $\Lambda$ goes through the following curve:

7.   Detailed bifurcation analysis about fixed points of the model (1.2)

For the completeness of this Section, the detailed bifurcation analysis about fixed points of the model (1.2) is presented as follows:

7.1.   Fold bifurcation about $P_{000}(0, 0, 0)$

About $P_{000}(0, 0, 0)$, $3\times3$ model (1.2) undergoes a fold bifurcation when parameters $\Lambda$ goes through the curve, which is depicted in Eq (6.1). But by computation, fold bifurcation do not occur and therefore $P_{000}(0, 0, 0)$ is degenerate with high co-dimension for the original parameters $\Lambda$ of $3\times3$ model (1.2).

7.2.   P-D bifurcation about $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$

About $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $3\times3$ model (1.2) undergoes P-D bifurcation when $\Lambda$ goes through the curve, which is depicted in Eq (6.3). Since model (1.2) is invariant $w.r.t$ $x = z = 0$ then we can restrict the model (1.2) on $x = z = 0$, to determine the bifurcation, where model (1.2) takes the form:

Also the map corresponds to Eq (7.1) takes the form:

Denote $r = r^* = 3, \ y = y^* = \frac{r-1}{r}$. By computation one gets: $f_y|_{r = r^* = 3, \ y = y^* = \frac{r-1}{r}} = -1$ which is non-hyperbolic condition for P-D bifurcation. Moreover, the non-degenerate conditions: $f_{yy}|_{r = r^* = 3, \ y = y^* = \frac{r-1}{r}} = -6\neq0$ and $f_r|_{r = r^* = 3, \ y = y^* = \frac{r-1}{r}} = \frac{2}{9}\neq0$ hold, which implies that $3\times3$ model (1.2) undergoes the P-D bifurcation when $\Lambda$ goes through the curve, which is depicted in Eq (6.3). Additionally, $\left(f_{yy}f_r+2f_{yr}\right)\left(2f_{yyy}+3f^2_{yy}\right)|_{r = r^* = 3, \ y = y^* = \frac{r-1}{r}} = -81 < 0$, and hence by existing literature ^[15,16,17] model (1.2) undergoes the supercritical P-D bifurcation.

7.3.   N-S bifurcation about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$

In this Subsection, N-S bifurcation of the $3\times3$ model (1.2) is explored about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, if $\Lambda$ goes through the curve, which is depicted in Eq (6.5) by assuming that $r\neq\frac{f\left(2c+ac+b\right)}{b\left(f-d\right)}$. Let

then $3\times3$ model (1.2) takes the form:

whose linear part is same as $J|_{P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)}$, which is depicted in Eq (4.13). Moreover at Eq (4.26), following non-degenerate conditions hold:

and additionally it is require that $\lambda^m_{2, 3}\neq 1 \ \mathrm{where}\ m = 1, 2, 3, 4$ which also true by computation. Now using the transformation:

Equation (7.4) takes the following form:

where

and

Hereafter, we have the system (7.7) on the center manifold, i.e.,

where

In view of Eqs (7.10) and (7.7) becomes:

After some manipulation, from Eq (7.12) one gets: $\alpha_1 = \alpha_2 = \alpha_3 = 0$. So, the system (7.7) restricted to $W^c(0)$ becomes

where

In order for Eq (7.13) to undergoes a N-S bifurcation, one require that the following discriminatory quantity $\sigma$ should be non-zero ^[12,18,19]:

where

After straightforward calculation, one gets:

Based on N-S bifurcation analysis about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, one has

Theorem 7.1. $3\times3$ model (1.2) undergoes a N-S bifurcation at $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ if $\sigma\neq0$. Moreover, if $\sigma < 0\ (\mathrm{resp.} \ \sigma > 0)$, then an attracting (resp. repelling) closed curve bifurcates from $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$. Finally Supercritical (resp. subcritical) N-S bifurcation take place if $\sigma < 0 \ (\mathrm{resp.} \ \sigma > 0)$.

Remark 1: Recall that about $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$, $3\times3$ model (1.2) may undergoes P-D bifurcation if $\Lambda$ goes through the curve, which is depicted in Eq (6.6). But by computation, P-D bifurcation do not occur and therefore $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ is degenerate with high co-dimension.

Hereafter, we explore the N-S as well as P-D bifurcations of the $3\times3$ model (1.2) without calculating the eigenvalues of $J|_{P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)}$ about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$. Precisely, in the following two theorems, we explore N-S and P-D bifurcations by utilizing explicit criteria, i.e., the criteria which is formulated by using a set of inequalities or equalities that contain the coefficients of the characteristic equation derived from $J|_{P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)}$ about the unique positive fixed point ^[20,21,22].

7.4.   N-S and P-D bifurcations about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$

Theorem 7.2. If

then $3\times3$ model (1.2) undergoes a N-S bifurcation about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$, where $\delta_1, \delta_2$ and $\delta_3$ are depicted in Eq (4.35).

Proof. By utilizing explicit criterion, which is depicted in ^[20,21], for $n = 3$ one gets:

Finally

Theorem 7.3. If

then $3\times3$ model (1.2) undergoes a P-D bifurcation about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$.

Proof. By utilizing explicit criterion, which is depicted in ^[21,22], for $n = 3$ one gets:

Finally

8.   Numerical simulations

Theoretical results are illustrated numerically in this Section. In this regard, following cases are presented to discuss the correctness of obtained theoretical results about fixed points for the $3\times3$ model (1.2):

Case I: If $a = 0.007, b = 0.001, c = 0.1, d = 0.39, e = 0.5, f = 2.7$ and varying $r\in [1.2, 10.4]$ then Figure 1a clearly indicates that about $P_{0y0}\left(0, 0.6884735202492211, 0\right)$, model (1.2) undergoes P-D bifurcation, which correctness the theoretical discussion in Subsection 7.2. More $3D$ bifurcation diagrams are plotted and depicted in Figure 1b, 1c.

Case II: If $a = 0.45, \ b = 0.9, c = 0.99, d = 0.29, e = 0.3, f = 2.3$ then from Eq (4.26) one gets: $r = 2.2772277227722775$. By Lemma 4.4, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ is locally asymptotically focus if $r < 2.2772277227722775$, exchange of stability if $r = 2.2772277227722775$, and unstable focus if $r > 2.2772277227722775$. For instance, if $r = 1.4 < 2.2772277227722775$ then it is clear from Figure 2a that $P_{0yz}\left(0, 0.12608695652173912, 0.2257356170399649\right)$ of model (1.2) is locally asymptotically focus. Similarly others values of $r$, if $r < 2.2772277227722775$ then $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ of model (1.2) is also locally asymptotically focus (Figure 2b, 2c). But if $r > 2.2772277227722775$ then discrete-time model (1.2) loss its stability and meanwhile stable closed curve appears (Figure 3a, 3c). The appearance of these closed curve implies that model (1.2) undergoes supercritical N-S bifurcation when parameters goes through the curve, which is depicted in Eq (6.5). This correctness of theoretical results obtain in Subsection 7.3. Finally, the bifurcation diagrams along with Maximum Lyapunov Exponent in this case, are plotted and drawn in Figure 4.

Case III: If $a = 0.45, \ b = 0.93, c = 0.8, d = 0.95, e = 0.15, f = 0.45$ and parameter $r$ vary from $1.23$ to $10.4$, i.e., $r\in[1.23, 10.4]$ then at $r = 2.86317$ the $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ loss its stability. That means that $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is locally stable if $r < 2.86317$. For instance, if $r = 1.75 < 2.86317$ then it is clear from Figure 5a that unique positive fixed point $P^+_{xyz}\left(19.044546850998465, 0.2737327188940094, 0.48387096774193544\right)$ of model (1.2) is locally stable. Similarly others values of $r$, if $r < 2.86317$ then $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ of model (1.2) is also locally stable (Figure 5b, 5c). Also if $a = 0.45, b = 0.93, c = 0.8, d = 0.95, e = 0.15, f = 0.45, r = 2.86317$ then Eq (4.34) along with Eq (4.35) takes the following form:

By computer computation, from Eq (8.1) one gets:

with

This grantee that for above parametric values, discrete-time model (1.2) undergoes N-S bifurcation when parameter $r$ goes through the bifurcation values, i.e., $r = 2.86317$. Moreover the parametric conditions, which are depicted in Eq (7.18), under which the model (1.2) undergoes a N-S bifurcation hold, i.e.,

which again grantee that model (1.2) undergoes a N-S bifurcation if parameters goes through the bifurcation value. Hence if $r > 2.86317$ then discrete-time model (1.2) loss its stability and meanwhile stable closed curve appears (Figure 6a–6c). The appearance of these closed curve implies that model (1.2) undergoes supercritical N-S bifurcation when parameters goes through $r > 2.86317$. This correctness of theoretical results obtain in Subsection 7.4. Finally, the bifurcation diagrams along with Maximum Lyapunov Exponent in this case, are plotted and drawn in Figure 7.

Case IV: If $a = 0.3, \ b = 1.99, c = 3.75, d = 0.75, e = 0.235, f = 1.0305$ and parameter $r$ vary from $1.6$ to $5.2$, i.e., $r\in[1.6, 5.2]$ then at $r = 3.9603484492658967$ the $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ loss its stability. That means that $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ is stable if $r < 3.9603484492658967$, unstable if $r = 3.9603484492658967$ and P-D bifurcation appears if $r > 3.9603484492658967$. For instance, if $a = 0.3, b = 1.99, c = 3.75, d = 0.75, e = 0.235, f = 1.0305, r = 3.97$ then Eq (4.34) along with Eq (4.35) takes the following form:

From Eq (8.2) one gets:

with

This grantee that for above parametric values, discrete-time model (1.2) undergoes P-D bifurcation when parameter $r$ goes through the bifurcation values, i.e., $r = 3.9603484492658967$. Moreover the parametric conditions, which are depicted in Eq (7.20), under which the model (1.2) undergoes a P-D bifurcation hold, i.e.,

which again grantee that model (1.2) undergoes a P-D bifurcation. The bifurcation diagrams along with Maximum Lyapunov Exponent in this case, are plotted and drawn in Figure 8. More figure are plotted and drawn in Figure 9 that indicates that model (1.2) complex dynamics of period-2, 10, 11, 12 respectively. For the qualitative behavior of the continuous-time model, we refer the reader ^[23,24,25] and references cited therein.

9.   Hybrid control of bifurcations

In discrete-time models, chaos control can be obtained using various methods, some of them are the state feedback method, pole-placement technique and hybrid control method. Among all these methods, hybrid control technique is the simplest to implement. Moreover this method can be implemented for controlling chaos under the influence of both N-S and P-D bifurcations. To control the chaos in model (1.2), we will utilize Hybrid control strategy to control the N-S and P-D bifurcations about the positive fixed point: $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$. For more detailed study about the Hybrid control strategy, we refer the reader to ^[26,27,28] and reference cited therein.

9.1.   Control of N-S bifurcation

The Hybrid control strategy is applied for control of N-S bifurcation in this Subsection. Recall that if the uncontrolled discrete-time model (1.2) undergoes a N-S bifurcation about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ then the controlled discrete-time model correspond to model (1.2) becomes:

where $0 < \alpha < 1$. Here it is noted that control model (9.1) reduces to model (1.2) if $\alpha = 1$, and chaos for model (9.1) can be delayed, advance or even completely eliminate by choosing the suitable choice of the parametric value of $\alpha$. Now if one reconsider the Case Ⅲ of Section 8 then $J|_{P^+_{xyz}\left(19.044546850998465, 0.20737327188940094, 0.48387096774193544\right)}$ about $P^+_{xyz}\left(19.044546850998465, 0.20737327188940094, 0.48387096774193544\right)$ for the controlled system (9.1) becomes

The characteristic polynomial of Eq (9.2) becomes

whose roots lies in an open disk $|\lambda| < 1$ iff $\alpha\in (0, 0.66675)$. For instance, if $\alpha = 0.56$ the phase portraits for controlled system (9.1) are plotted and drawn in Figure 10.

9.2.   Control of P-D bifurcation

Recall that if Case Ⅳ of Section 8 hold then $J|_{P^+_{xyz}\left(13.30555251494148, 0.604750275583457, 0.1507537688442211\right)}$ about $P^+_{xyz}\left(13.30555251494148, 0.604750275583457, 0.1507537688442211\right)$ for the controlled system, which is depiction in system (9.1), becomes

The characteristic polynomial of Eq (9.4) becomes

whose roots lies in an open disk $|\lambda| < 1$ iff $\alpha\in (0, 0.543543)$. For instance, if $\alpha = 0.50$ the phase portraits for controlled system (9.1) are plotted and drawn in Figure 11.

10.   Conclusions and future work

This work deals with the study of local dynamics along with different topological classifications, bifurcations and Hybrid control of the model (1.2) in $\mathbb{R}_+^3$. In the interior of $\mathbb{R}_+^3$, it is explored that the model has boundary fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and a unique positive fixed point: $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ under certain parametric conditions. We have explored the local dynamics along with different topological classifications about fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ and $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ by applying the method of Linearization and stability theory, and conclusions are presented in Table 2. We have also investigated the local dynamics about positive fixed point: $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$, and proved that it is a sink if and only if Eq (4.36) along with Eq (4.35) hold. We have also investigated the existence of prime period and periodic points for the model (1.2). Further, about each fixed point of the under consideration model, we have discussed the possible bifurcation scenarios. In general, bifurcation theory is the mathematical study of changes in the qualitative or topological structure of a given family, such as the integral curves of a family of vector fields and the solutions of family of differential equations. It is most commonly applied to the mathematical study of dynamical systems, a bifurcation occurs when a small smooth change made to the parameter values namely the bifurcation parameters of the system causes a sudden qualitative or topological change in its behavior. So, in the current study, we proved that about $P_{000}(0, 0, 0)$ there exist a fold bifurcation as well as P-D bifurcation if $\Lambda$ respectively goes through the curves: $FB|_{P_{000}(0, 0, 0)} = \left\{\Lambda: r = 1\right\}$ and $PDB|_{P_{000}(0, 0, 0)} = \left\{\Lambda: d = 2\right\}$. But by computation it goes to degenerate to high dimension and hence there does not such bifurcations about $P_{000}(0, 0, 0)$. About $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, it is proved that there exist a P-D bifurcation when $\Lambda$ goes through the curve: $PDB|_{P_{0y0}\left(0, \frac{r-1}{r}, 0\right)} = \left\{\Lambda: r = 3\right\}$. We have also proved that about $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$ the model (1.2) undergoes a N-S bifurcation when $\Lambda$ goes through the curve: $NSB|_{P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)} = \left\{\frac{4d^2fr+d^2r^2-4drf^2+4d^2fr}{f^2} < 0, \ r = \frac{f}{f-1-d}\right\}$. By explicit criterion(without finding eigenvalues) it is proved that about $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$ the model (1.2) undergoes N-S as well as P-D bifurcations. Numerical simulations are presented to verify theoretical results. Further, we have computed maximum Lyapunov exponents numerically. At the end, hybrid control strategy is applied to stabilize chaos existing in the model (1.2). Global dynamics about fixed points: $P_{000}(0, 0, 0)$, $P_{0y0}\left(0, \frac{r-1}{r}, 0\right)$, $P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right)$, $P_{x0z}\left(\frac{d}{e}, 0, \frac{a}{b}\right)$ and $P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right)$, existence of global bifurcations and construction of forbidden set for the model (1.2) are our next aim to study.

Acknowledgments

This research is partially supported by the Higher Education Commission of Pakistan.

Conflict of interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.