On some models for cancer cell migration through tissue networks

  • We propose some models allowing to account for relevant processes at the various scales of cancer cell migration through tissue, ranging from the receptor dynamics on the cell surface over degradation of tissue fibers by protease and soluble ligand production towards the behavior of the entire cell population.
       For a genuinely mesoscopic version of these models we also provide a result on the local existence and uniqueness of a solution for all biologically relevant space dimensions.

    Citation: Jan Kelkel, Christina Surulescu. On some models for cancer cell migration throughtissue networks[J]. Mathematical Biosciences and Engineering, 2011, 8(2): 575-589. doi: 10.3934/mbe.2011.8.575

    Related Papers:

    [1] Ya-Lei Li, Da-Bin Wang, Jin-Long Zhang . Sign-changing solutions for a class of p-Laplacian Kirchhoff-type problem with logarithmic nonlinearity. AIMS Mathematics, 2020, 5(3): 2100-2112. doi: 10.3934/math.2020139
    [2] Xia Li, Wen Guan, Da-Bin Wang . Least energy sign-changing solutions of Kirchhoff equation on bounded domains. AIMS Mathematics, 2022, 7(5): 8879-8890. doi: 10.3934/math.2022495
    [3] Jinfu Yang, Wenmin Li, Wei Guo, Jiafeng Zhang . Existence of infinitely many normalized radial solutions for a class of quasilinear Schrödinger-Poisson equations in $ \mathbb{R}^3 $. AIMS Mathematics, 2022, 7(10): 19292-19305. doi: 10.3934/math.20221059
    [4] Qing Yang, Chuanzhi Bai . Sign-changing solutions for a class of fractional Kirchhoff-type problem with logarithmic nonlinearity. AIMS Mathematics, 2021, 6(1): 868-881. doi: 10.3934/math.2021051
    [5] Ziqing Yuan, Jing Zhao . Solutions for gauged nonlinear Schrödinger equations on $ {\mathbb R}^2 $ involving sign-changing potentials. AIMS Mathematics, 2024, 9(8): 21337-21355. doi: 10.3934/math.20241036
    [6] Ye Xue, Zhiqing Han . Existence and multiplicity of solutions for Schrödinger equations with sublinear nonlinearities. AIMS Mathematics, 2021, 6(6): 5479-5492. doi: 10.3934/math.2021324
    [7] Jie Yang, Lintao Liu, Haibo Chen . Multiple positive solutions to the fractional Kirchhoff-type problems involving sign-changing weight functions. AIMS Mathematics, 2024, 9(4): 8353-8370. doi: 10.3934/math.2024406
    [8] Zhongxiang Wang, Gao Jia . Existence of solutions for modified Kirchhoff-type equation without the Ambrosetti-Rabinowitz condition. AIMS Mathematics, 2021, 6(5): 4614-4637. doi: 10.3934/math.2021272
    [9] Zhongxiang Wang . Existence and asymptotic behavior of normalized solutions for the modified Kirchhoff equations in $ \mathbb{R}^3 $. AIMS Mathematics, 2022, 7(5): 8774-8801. doi: 10.3934/math.2022490
    [10] Yuan Shan, Baoqing Liu . Existence and multiplicity of solutions for generalized asymptotically linear Schrödinger-Kirchhoff equations. AIMS Mathematics, 2021, 6(6): 6160-6170. doi: 10.3934/math.2021361
  • We propose some models allowing to account for relevant processes at the various scales of cancer cell migration through tissue, ranging from the receptor dynamics on the cell surface over degradation of tissue fibers by protease and soluble ligand production towards the behavior of the entire cell population.
       For a genuinely mesoscopic version of these models we also provide a result on the local existence and uniqueness of a solution for all biologically relevant space dimensions.


    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]:

    $ dxdt=axbxy, dydt=cxydy,
    $
    (1.1)

    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]:

    $ xn+1=(1+a)xnbxnzn, yn+1=ryn(1yn)cynzn, zn+1=(1d)zn+exnzn+fynzn,
    $
    (1.2)

    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.

    Table 1.  Summarized results of Sagayaraj et al. [10] 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)$ of the model (1.2).
    Fixed pointsCorresponding behavior for the model (1.2)
    $P_{000}$sink if $r < 1, -2 < a < 0\ \mathrm{and}\ 0 < d < 2$; source if $a > 0, r > 1 \ \mathrm{and}\ d > 2$;
    saddle if $a > 0, r < 1 \ \mathrm{and}\ 0 < d < 2$; non-hyperbolic if $r = 1$ or $d = 2$.
    $P_{0y0}$sink if $1 < r < 3$, $-2 < a < 0 \ \mathrm{and}\ r < \frac{f}{f-d}$; source if $r > 3, a > 0 \ \mathrm{and}\ r > \frac{f}{f-d}$;
    saddle if $1 < r < 3, a > 0 \ \mathrm{and}\ r < \frac{f}{f-d}$.
    $P_{0yz}$ sink if $r > \frac{f(b+ac)}{b(f-d)}\ \rm{and} \ r > \frac{f}{f-d}$; source if $r < \frac{f(b+ac)}{b(f-d)}\ \rm{and} \ r < \frac{f}{f-d}$;
    saddle if $\frac{f}{f-d} < r < \frac{f(2c+ac+b)}{b(f-d)}$.

     | Show Table
    DownLoad: CSV

    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.

    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

    $ x=(1+a)xbxz, y=ry(1y)cyz, z=(1d)z+exz+fyz.
    $
    (2.1)

    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:

    $ abz=0, r(1y)cz=1, d+ex+fy=0.
    $
    (2.2)

    From first equation of system (2.2) one gets:

    $ z=ab.
    $
    (2.3)

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

    $ y=brbacbr.
    $
    (2.4)

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

    $ x=br(df)+f(b+ac)ber.
    $
    (2.5)

    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 $.

    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):

    $ (f,g,h)(xn+1,yn+1,zn+1),
    $
    (3.1)

    where

    $ f=(1+a)xbxz, g=ry(1y)cyz, h=(1d)z+exz+fyz.
    $
    (3.2)

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

    $ J|Pxyz(x,y,z)=(1+abz0bx0r2ryczcyezfz1d+ex+fy).
    $
    (3.3)

    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.

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

    $ J|P000(0,0,0)=(1+a000r0001d).
    $
    (4.1)

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

    $ λ1=1+a, λ2=r, λ3=1d.
    $
    (4.2)

    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

    $ a>0, r>1 and d>2;
    $
    (4.3)

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

    $ a>0, 0<r<1 and 0<d<2;
    $
    (4.4)

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

    $ r=1,
    $
    (4.5)

    or

    $ d=2.
    $
    (4.6)

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

    $ J|P0y0(0,r1r,0)=(1+a0002rc(1r)r001d+f(r1)r).
    $
    (4.7)

    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

    $ λ1=1+a, λ2=2r, λ3=1d+f(r1)r.
    $
    (4.8)

    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

    $ a>0, and 3<r<ffd+2;
    $
    (4.9)

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

    $ a>0, and ffd+2<r<3;
    $
    (4.10)

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

    $ r=3,
    $
    (4.11)

    or

    $ r=ffd+2.
    $
    (4.12)

    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

    $ J|P0yz(0,df,rffdrcf)=(1+a+bc(1r+drf)0001drfdcfec(r1drf)fc(r1drf)1).
    $
    (4.13)

    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

    $ λ1=1+a+bc(1r+drf),λ2,3=2fdrf±Δ2,
    $
    (4.14)

    where

    $ Δ=(dr2ff)24(fdr+drffdd2rf),=4df2+d2r24drf2+4d2frf2.
    $
    (4.15)

    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

    $ f(ac+b)b(fd)<r<f(2c+ac+b)b(fd)  and  ffd<r<f(4d)d(d+2f);
    $
    (4.16)

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

    $ r>max{f(4d)d(d+2f),f(2c+ac+b)b(fd)};
    $
    (4.17)

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

    $ f(ac+b)b(fd)<r<f(2c+ac+b)b(fd)  and  r>f(4d)d(d+2f),
    $
    (4.18)

    or

    $ r>f(2c+ac+b)b(fd)  and  ffd<r<f(4d)d(d+2f),
    $
    (4.19)

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

    $ r=f(4d)d(d+2f),
    $
    (4.20)

    or

    $ r=f(2c+ac+b)b(fd).
    $
    (4.21)

    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

    $ f(ac+b)b(fd)<r<f(2c+ac+b)b(fd)  and  r<ff1d;
    $
    (4.22)

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

    $ r>max{f(2c+ac+b)b(fd),ff1d};
    $
    (4.23)

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

    $ f(ac+b)b(fd)<r<f(2c+ac+b)b(fd)  and  r>ff1d,
    $
    (4.24)

    or

    $ r>f(2c+ac+b)b(fd)  and  r<ff1d;
    $
    (4.25)

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

    $ r=ff1d,
    $
    (4.26)

    or

    $ r=f(2c+ac+b)b(fd).
    $
    (4.27)

    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

    $ J|Px0z(de,0,ab)=(10bde0racb0acbafb1).
    $
    (4.28)

    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

    $ λ1=racb,λ2,3=1±ιab.
    $
    (4.29)

    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

    $ r<acb1;
    $
    (4.30)

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

    $ r>acb1;
    $
    (4.31)

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

    $ r=acb1.
    $
    (4.32)

    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

    $ J|P+xyz(br(df)+f(b+ac)ber,brbacbr,ab)=(10br(df)+f(b+ac)re02r+acbc(b+ac)brcaebafb1).
    $
    (4.33)

    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

    $ P(λ)=λ3+δ1λ2+δ2λ+δ3,
    $
    (4.34)

    where

    $ δ1=r4acb,δ2=52r+adaf+2acbacfbra2c2fb2r+acfb+afr+a2cfbr,δ3=r2+adrafr2ad+3af+acb+acfbr+a2c2fb2racfb3a2cfbr+2a2cfba2cdb2afra3c2fb2r.
    $
    (4.35)

    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:

    $ |δ1+δ3|<1+δ2,  |δ13δ3|<3δ2,  δ23+δ2δ3δ1<1,
    $
    (4.36)

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

    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

    $ F(x,y,z):=(f(x,y,z),g(x,y,z),h(x,y,z)),
    $
    (5.1)

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

    $ F|P000(0,0,0)=P000(0,0,0),F|P0y0(0,r1r,0)=P0y0(0,r1r,0),F|P0yz(0,df,rffdrcf)=P0yz(0,df,rffdrcf),F|Px0z(de,0,ab)=Px0z(de,0,ab),F|P+xyz(br(df)+f(b+ac)ber,brbacbr,ab)=P+xyz(br(df)+f(b+ac)ber,brbacbr,ab).
    $
    (5.2)

    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

    $ F2(x,y,z)=((1+a)f(x,y,z)bf(x,y,z)h(x,y,z), rg(x,y,z)(1g(z,y,z))cg(x,y,z)h(x,y,z), (1d)h(x,y,z)+ef(x,y,z)h(x,y,z)+ fg(x,y,z)h(x,y,z))F2|P000(0,0,0)=P000(0,0,0),F3(x,y,z)=((1+a)f2(x,y,z)bf2(x,y,z)h2(x,y,z), rg2(x,y,z)(1g2(z,y,z))cg2(x,y,z)h2(x,y,z),(1d)h2(x,y,z)+ef2(x,y,z)h2(x,y,z)+ fg2(x,y,z)h2(x,y,z))F2|P000(0,0,0)=P000(0,0,0), Fi(x,y,z)=((1+a)fi1(x,y,z)bfi1(x,y,z)hi1(x,y,z), rgi1(x,y,z)(1gi1(z,y,z))cgi1(x,y,z)hi1(x,y,z),(1d)hi1(x,y,z)+efi1(x,y,z)hi1(x,y,z)+ fgi1(x,y,z)hi1(x,y,z))Fi|P000(0,0,0)=P000(0,0,0).
    $
    (5.3)

    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.

    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:

    $ FB|P000(0,0,0)={Λ:r=1}.
    $
    (6.1)

    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:

    $ PDB|P000(0,0,0)={Λ:d=2}.
    $
    (6.2)

    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:

    $ PDB|P0y0(0,r1r,0)={Λ:r=3}.
    $
    (6.3)

    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:

    $ PDB|P0y0(0,r1r,0)={Λ:r=ffd+2}.
    $
    (6.4)

    $ \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:

    $ NSB|P0yz(0,df,rffdrcf)={4d2fr+d2r24drf2+4d2frf2<0, r=ff1d}.
    $
    (6.5)

    $ \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:

    $ PDB|Px0z(de,0,ab)={Λ:r=acb1}.
    $
    (6.6)

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

    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).

    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:

    $ yn+1=ryn(1yn).
    $
    (7.1)

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

    $ f(y):=ryry2.
    $
    (7.2)

    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.

    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

    $ un=xn, vn=yndf, wn=znrffdrcf,
    $
    (7.3)

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

    $ (un+1vn+1vn+1)=(1+a+bc(1r+drf)0001drfdcfec(r1drf)fc(r1drf)1)(unvnvn)+(bununrv2ncvnwneunwn+fvnwn),
    $
    (7.4)

    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:

    $ |λ2,3|(4.26)=1 and d|λ2,3|dr|(4.26)=fddd22f0,
    $
    (7.5)

    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:

    $ (unvnwn)=T(xtytzt),
    $
    (7.6)

    Equation (7.4) takes the following form:

    $ (xt+1yt+1zt+1)=(2fdr2f4drf24df24d2frr2d22f04drf24df24d2frr2d22f2fdr2f0001+a+bc(1r+drf))×(xtytzt)+(Γ(xt,yt,zt)Φ(xt,yt,zt)Ψ(xt,yt,zt)),
    $
    (7.7)

    where

    $ T=(00c2dfe(cdr+acfbfr+bf+bdr)bf+acf+bdrbfre(f+drfr)cdr2f(frfdr)c4drf24df24d2frd2r22f(frfdr)c2dbf+acf+bdr+cdrbfr101),
    $
    (7.8)

    and

    $ Γ(xt,yt,zt)=bxtzt+bz2t+(xtzt+z2t)(c2dfbf+acf+bdr+cdrbfrbf+acf+bdrbfrf+drfr)cdr2(frfdr)x2t+c4drf24df24d2frr2d22(frfdr)xtytc2fdbf+acf+bdr+cdrbfrztxtcdr2(frfdr)xtzt+c4drf24df24d2frr2d22(frfdr)ytztc2fdbf+acf+bdr+cdrbfrz2t,Φ(xt,yt,zt)=bcdr2bc2df(frfdr)bf+acf+bdr+cdrbfrc4drf24df24d2frr2d2(xtzt+z2t)+2f(frfdr)c4drf24df24d2frr2d2×[r(cdr2f(frfdr)xt+c4drf24df24d2frr2d22f(frfdr)ytc2dbf+acf+bdr+cdrbfrzt)2+(xt+zt)(c2dr2f(frfdr)xtc24drf24df24d2frr2d22f(frfdr)yt+c3dbf+acf+bdr+cdrbfrzt)]+dr4drf24df24d2frr2d2×[zt(xt+zt)(c2dfbf+acf+bdr+cdrbfrbf+acf+bdrbfr(f+drfr))+(xt+zt)(cdr2(frfdr)xt+c4drf24df24d2frr2d22(frfdr)ytc2dfbf+acf+bdr+cdrbfrzt)],Ψ(xt,yt,zt)=bxtztbz2t.
    $
    (7.9)

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

    $ Wc(0)={(xt,yt,zt)zt=h(xt,yt),h(0,0)=0,Dh(0,0)=0},
    $
    (7.10)

    where

    $ h(xt,yt)=α1x2t+α2xtyt+α3y2t+O((|xt|+|yt|)3).
    $
    (7.11)

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

    $ h(2fdr2fxt+4drf24df24d2frr2d22fyt+Γ(xt,yt,h),4drf24df24d2frr2d22fxt+2fdr2fyt+Φ(xt,yt,h))=(1+abc(r1drf))h(xt,yt)+Ψ(xt,yt,h).
    $
    (7.12)

    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

    $ (xt+1yt+1)=(2fdr2f4drf24df24d2frr2d22f4drf24df24d2frr2d22f2fdr2f)(xtyt)+(ζ1(xt,yt)ζ2(xt,yt)),
    $
    (7.13)

    where

    $ ζ1(xt,yt)=cdr2(frfdr)x2t+c4drf24df24d2frr2d22(frfdr)xtyt,ζ2(xt,yt)=2f(frfdr)c4drf24df24d2frr2d2[r(cdr2f(frfdr)xt+c4drf24df24d2frr2d22f(frfdr)yt)2+(c2dr2f(frfdr))x2tc24drf24df24d2frr2d22f(frfdr)xtyt]+dr4drf24df24d2frr2d2(cdr2(frfdr)x2t+c4drf24df24d2frr2d22(frfdr)xtyt).
    $
    (7.14)

    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]:

    $ σ=Re((12λ)ˉλ21λL11L12)12|L11|2|L21|2+Re(ˉλL22),
    $
    (7.15)

    where

    $ L11=14[ˉζ1xtxt+ˉζ1ytyt+ι(ˉζ2xtxt+ˉζ2ytyt)],L12=18[ˉζ1xtxtˉζ1ytyt+2ˉζ2xtyt+ι(ˉζ2xtxtˉζ2ytyt2ˉζ1xtyt)],L21=18[ˉζ1xtxtˉζ1ytyt2ˉζ2xtyt+ι(ˉζ2xtxtˉζ2ytyt+2ˉζ1xtyt)],L22=116[ˉζ1xtxtxt+ˉζ1xtytyt+ˉζ2xtxtyt+ˉζ2ytytyt+ι(ˉζ2xtxtxt+ˉζ2xtytytˉζ1xtxtytˉζ1ytytyt)].
    $
    (7.16)

    After straightforward calculation, one gets:

    $ L11=14[cdrfrfdr+ι(14drf24df24d2frr2d2(2c2drc2d2r3f(frfdr))cd2r2(frfdr)4drf24df24d2frr2d2rc4drf24df24d2frr2d2f(frfdr))],L12=18[cdrf+2cdr22cf2r+2cf2+cdrf(frfdr)+ι(1c4drf24df24d2frr2d2×(2c2drc2d2r3f(frfdr))cd2r2(frfdr)4drf24df24d2frr2d2+c4drf24df24d2frr2d2f(frfdr)(rf))],L21=18[2cf2r3cdrf2cdr22cf2cdrf(frfdr)+ι(1c4drf24df24d2frr2d2(2c2drc2d2r3f(frfdr))cd2r2(frfdr)4drf24df24d2frr2d2+c4drf24df24d2frr2d2f(frfdr)(r+f))],L22=0.
    $
    (7.17)

    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].

    Theorem 7.2. If

    $ 1δ2+δ3(δ1δ3)=0,1+δ2δ3(δ1+δ3)>0,1+δ1+δ2+δ3>0,1δ1+δ2δ3>0,ddr(1δ2+δ3(δ1δ3))|r=r00,cos2πl11+δ1+δ2+δ32(1+δ2), l=3,4,,
    $
    (7.18)

    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:

    $ Δ2(r)=1δ2+δ3(δ1δ3)=0,Δ+2(r)=1+δ2δ3(δ1+δ3)>0,Pr(1)=1+δ1+δ2+δ3>0,(1)3Pr(1)=1δ1+δ2δ3>0,ddr(Δ2(r))|r=r0=ddr(1δ2+δ3(δ1δ3))|r=r00.
    $
    (7.19)

    Finally

    $ 10.5Pr(1)Δ0(r)Δ+1(r)=11+δ1+δ2+δ32(1+δ2).
    $

    Theorem 7.3. If

    $ 1δ2+δ3(δ1δ3)>0,1+δ2δ3(δ1+δ3)>0,1±δ2>0,1+δ1+δ2+δ3>0,1+δ1δ2+δ3=0,δ1δ2+δ332δ1+δ20,
    $
    (7.20)

    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:

    $ Δ2(r)=1δ2+δ3(δ1δ3)>0,Δ+2(r)=1+δ2δ3(δ1+δ3)>0,Δ±1(r)=1±δ2>0,Pr(1)=1+δ1+δ2+δ3>0,Pr(1)=1+δ1δ2+δ3=0.
    $
    (7.21)

    Finally

    $ ni=1(1)niδini=1(1)ni(ni+1)δi1=δ1δ2+δ332δ1+δ20.
    $

    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.

    Figure 1.  P-D bifurcation diagram about $ P_{0y0}\left(0, \frac{r-1}{r}, 0\right) $ of model (1.2) by varying $ r\in [1.2, 10.4] $.

    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.

    Figure 2.  Topological classification of model (1.2) about $ P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right) $ if $ r < 2.2772277227722775 $.
    Figure 3.  Supercritical N-S bifurcation of model (1.2) about $ P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right) $ if $ r > 2.2772277227722775 $.
    Figure 4.  Bifurcation diagram and their corresponding Maximum Lyapunov Exponent of the model (1.2) about $ P_{0yz}\left(0, \frac{d}{f}, \frac{rf-f-dr}{cf}\right) $. (4a, 4b) Bifurcation diagram of the model if $ r\in[1.2, 10.4] $ and initial condition $ (0, 0.17, 0.54) $. (4c) Maximum Lyapunov Exponent corresponding to Figure 4a, 4b.

    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:

    $ P(λ)=λ31.5239267741935483λ2+0.15468797005336343λ+0.4321468661912955=0.
    $
    (8.1)
    Figure 5.  Topological classification of model (1.2) about $ P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right) $ if $ r < 2.86317 $.

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

    $ \lambda_1 = -0.43214706529074226, \ \lambda_{2, 3} = 0.9780369197421452\pm0.208430618911403\iota $

    with

    $ |\lambda_{2, 3}| = |0.9780369197421452\pm0.208430618911403\iota| = 1. $

    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.,

    $ 1δ2+δ3(δ1δ3)=0,1+δ2δ3(δ1+δ3)=1.6264972358671579>0,1+δ1+δ2+δ3=0.06290806205111066>0,1δ1+δ2δ3=2.246467878055616>0,ddr(1δ2+δ3(δ1δ3))|r=2.86317=1.94299184399591770,cos2πl11+δ1+δ2+δ32(1+δ2)=0.9727597135838336,  l=3,4,5,,
    $

    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 6a6c). 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.

    Figure 6.  Supercritical N-S bifurcation of model (1.2) about $ P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right) $ if $ r > 2.86317 $.
    Figure 7.  Bifurcation diagram and their corresponding Maximum Lyapunov Exponent of the model (1.2) about $ P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right) $. (7a, 7b) Bifurcation diagram of the model if $ r\in[1.1, 10.4] $ and initial condition $ (0.2, 0.2, 0.2) $. (7c) Maximum Lyapunov Exponent corresponding to Figure 7a, 7b.

    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:

    $ P(λ)=λ30.6049781838999324λ20.8791823027251384λ+0.7257958811747927=0.
    $
    (8.2)

    From Eq (8.2) one gets:

    $ \lambda_1 = -1, \ \lambda_{2, 3} = 0.802489091949966\pm 0.28601947219745705\iota, $

    with

    $ |\lambda_{2, 3}| = |0.802489091949966\pm 0.28601947219745705\iota| = 0.8519365476224112\not = 1. $

    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.,

    $ 1δ2+δ3(δ1δ3)=0.9133119675196675>0,1+δ2δ3(δ1+δ3)=0.03312871021974505>0,1+δ2=0.12081769727486158>0,1δ2=1.8791823027251384>0,1+δ1+δ2+δ3=0.24163539454972183>0,1+δ1δ2+δ3=0,δ1δ2+δ332δ1+δ2=0.078240040995743830,
    $

    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.

    Figure 8.  Bifurcation diagram and their corresponding Maximum Lyapunov Exponent of the model (1.2) about $ P^+_{xyz}\left(\frac{br(d-f)+f(b+ac)}{ber}, \frac{br-b-ac}{br}, \frac{a}{b}\right) $. (8a, 8b) Bifurcation diagram of the model if $ r\in[1.5, 5.2] $ and initial condition $ (0.5, 0.15, 0.15) $. (8c) Maximum Lyapunov Exponent corresponding to Figure 8a, 8b.
    Figure 9.  Complex dynamics of model (1.2) in the Case Ⅳ.

    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.

    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:

    $ xn+1=α((1+a)xnbxnzn)+(1α)xn,yn+1=α(ryn(1yn)cynzn)+(1α)yn,zn+1=α((1d)zn+exnzn+fynzn)+(1α)zn,
    $
    (9.1)

    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

    $ J|P+xyz(19.044546850998465,0.20737327188940094,0.48387096774193544)=(106.63919α011.47607α0.41243α0.0725806α1.42742α1).
    $
    (9.2)

    The characteristic polynomial of Eq (9.2) becomes

    $ F(λ)=λ3+(1.47607α3)λ2+(0.106834α22.95215α+3)λ0.711285α30.106834α2+1.47607α1=0,
    $
    (9.3)

    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.

    Figure 10.  Phase portraits for controlled system (9.1).

    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

    $ J|P+xyz(13.30555251494148,0.604750275583457,0.1507537688442211)=(1010.4777α012.39502α2.26781α0.0354271α0.565327α1).
    $
    (9.4)

    The characteristic polynomial of Eq (9.4) becomes

    $ F(λ)=λ3+(2.39502α3)λ2+(0.910861α24.79004α+3)λ0.889018α30.9108661α2+2.39502α1=0,
    $
    (9.5)

    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.

    Figure 11.  Phase portraits for controlled system (9.1).

    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.

    Table 2.  Fixed points with corresponding behavior of the $3\times 3$ model (1.2).
    Fixed pointsCorresponding behavior
    $P_{000}$never sink; source if $r > 1\ \mathrm{and}\ \ d > 2$; saddle if $0 < r < 1\ \mathrm{and}\ 0 < d < 2$;
    non-hyperbolic if $r = 1 \ \ {\rm{or}} \ \ d = 2$.
    $P_{0y0}$never sink; source if $3 < r < \frac{f}{f-d+2}$; saddle if $\frac{f}{f-d+2} < r < 3$;
    non-hyperbolic if $r = 3 \ {\rm{or}}\ r = \frac{f}{f-d+2}$.
    $P_{0yz}$stable node if
    $\frac{f(ac+b)}{b(f-d)} < r < \frac{f(2c+ac+b)}{b(f-d)}\ \ \mathrm{and}\ \ \frac{f}{f-d} < r < \frac{f(4-d)}{d(d+2-f)}$;
    unstable node if
    $ r > \max\left\{\frac{f(4-d)}{d(d+2-f)}, \frac{f(2c+ac+b)}{b(f-d)}\right\}; $
    saddle node if
    $\frac{f(ac+b)}{b(f-d)} < r < \frac{f(2c+ac+b)}{b(f-d)}\ \ \mathrm{and}\ \ r > \frac{f(4-d)}{d(d+2-f)}$;
    or
    $r > \frac{f(2c+ac+b)}{b(f-d)}\ \ \mathrm{and}\ \ \frac{f}{f-d} < r < \frac{f(4-d)}{d(d+2-f)}$;
    non-hyperbolic if
    $r = \frac{f(4-d)}{d(d+2-f)} \ {\rm{or}} \ \ r = \frac{f(2c+ac+b)}{b(f-d)}; $
    stable focus-node if
    $\frac{f(ac+b)}{b(f-d)} < r < \frac{f(2c+ac+b)}{b(f-d)}\ \ \mathrm{and}\ \ r < \frac{f}{f-1-d}; $
    unstable focus-node if
    $r > \max\left\{\frac{f(2c+ac+b)}{b(f-d)}, \frac{f}{f-1-d}\right\}; $
    saddle focus-node if
    $\frac{f(ac+b)}{b(f-d)} < r < \frac{f(2c+ac+b)}{b(f-d)}\ \ \mathrm{and}\ \ r > \frac{f}{f-1-d}$;
    or
    $r > \frac{f(2c+ac+b)}{b(f-d)}\ \ \mathrm{and}\ \ r < \frac{f}{f-1-d}$;
    non-hyperbolic if
    $r = \frac{f}{f-1-d}\ {\rm{or}} \ \ r = \frac{f(2c+ac+b)}{b(f-d)}.$
    $P_{x0z}$never stable focus-node; unstable focus-node if $r < \frac{ac}{b}-1$;
    saddle focus-node if $r > \frac{ac}{b}-1$; non-hyperbolic if $r = \frac{ac}{b}-1$.
    $P_{xyz}$sink if Eq (4.36) along with Eq (4.35) hold.

     | Show Table
    DownLoad: CSV

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

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

  • This article has been cited by:

    1. Sabrina H. Streipert, Gail S. K. Wolkowicz, 2023, Chapter 22, 978-3-031-25224-2, 473, 10.1007/978-3-031-25225-9_22
    2. Subhadipa Das, Amar Nath Chatterjee, Sourav Rana, Complex Dynamics of a Discrete Time One Prey Two Predator System with Prey Refuge, 2025, 11, 2349-5103, 10.1007/s40819-025-01896-9
  • Reader Comments
  • © 2011 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3130) PDF downloads(466) Cited by(15)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog