Existence of a Limit Cycle in an Intraguild Food Web Model with Holling Type II and Logistic Growth for the Common Prey

Francisco Eduardo Castillo-Santos, Miguel Angel Dela Rosa, Iván Loreto-Hernández

División Académica de Ciencias Básicas, CONACyT-UJAT, Cunduacán, México.

DOI: 10.4236/am.2017.83030                

Abstract

In this paper, we prove the existence of a limit cycle for a given system of differential equations corresponding to an asymmetrical intraguild food web model with functional responses Holling type II for the middle and top predators and logistic grow for the (common) prey. The existence of such limit cycle is guaranteed, via the first Lyapunov coefficient and the Andronov-Hopf bifurcation theorem, under certain conditions for the parameters involved in the system.

Keywords

Hopf’s Bifurcation, Limit Cycle, Intraguild Model

1. Introduction

It is well known that interaction between three species, in which predation and competition occurs, is called intraguild predation (see [1] ). This kind of in- teraction can take place in a group of species that exploit the same resources in a similar way (see [2] ). This kind of interaction among the species in an intraguild model is of particular interest. One of the main questions when looking at inter- action of species is whether or not there will be coexistence among them. This is of importance from the ecological point of view. In the intraguild predation model, one can consider two cases, the symmetric which occurs when there is a mutual predation between two species, and the asymmetric that occurs when one species, usually called intraguild predator, always predate the middle species, which is called the intraguild prey. In both cases it is assumed that the corres- ponding species use common foods (see [1] ).

The criterion to have coexistence in the asymmetric intraguild predation sys- tem seems to be, on one hand, to impose conditions on the intraguild prey, that is, it should be superior at the competition for the resources in comparison with the intraguild predator, and on the other hand, that the intraguild predator should be substantially benefit from the consumption to the intraguild prey in the sense that its most important food source is intermediate species (see [3] ).

There are some recent papers where food chain models between three species have been studied in which the authors have obtained results about the co- existence of the species by looking at the existence of limit cycles for the corresponding model systems, for instance tritrophic models with linear growth prey (see [4] [5] [6] ) and logistic growth prey (see [7] ). These models do not consider predation of the top predator to the resource (the prey). Hence one can see that intraguild predation is a more complex interaction between species that the tritrophic model.

If the growth rate for the resource is linear, we are assuming that the density of the resource is growing exponentially. When it is assumed logistic growth rate for the resource, the corresponding carry capacity implies that the resource density is bounded, which has Ecological sense but it seems to be more difficult to have a coexistence between the species.

In this paper, we are interested in guaranteeing the coexistence of three species forming an intraguild food web model, which is an asymmetrical intraguild predation model with functional response $f_1\left(x\right)$ for the middle predator species, and functional responses $f_2\left(x\right)\mathrm{,}{f}_3\left(y\right)$ for the top predator, and logistic grow for the prey. More precisely models with the form:

$\begin{array}{l}\stackrel{\dot{}}{x}=\rho x\left(1-\frac{x}{R}\right)-y{f}_1\left(x\right)-z{f}_2\left(x\right),\\ \stackrel{\dot{}}{y}=c_{1}y{f}_1\left(x\right)-c_{2}y-z{f}_3\left(y\right),\\ \stackrel{\dot{}}{z}=z\left(c_3{f}_3\left(y\right)+c_4{f}_2\left(x\right)-\mu \right),\end{array}$ (1)

where $x$ represents the density of a prey that gets eaten by a species of density $y$ (mesopredator) and a species of density $z$ (super-predator), and the species $y$ feeds the species $z$ . Moreover $R$ represents the carry capacity of the prey and $\rho$ represents the growth rate of the prey. The parameters $c_1\mathrm{,}{c}_3$ and $c_4$ are positive constants which represent the benefit from the consumption of food and the parameters $c_2$ and $\mu$ represent the mortality rate of the correspon- ding predators. We will consider that the functions $f_1\mathrm{,}{f}_2$ and $f_3$ given in (1) are Holling type II, that is

$\begin{array}{l}{f}_1\left(x\right)=\frac{x}{a_{1}x+b_1},\\ f_2\left(x\right)=\frac{x}{a_{2}x+b_2}\text{and}{f}_3\left(y\right)=\frac{y}{a_{3}y+b_3}\mathrm{,}\end{array}$

where $a_1,b_1,a_2,b_2,a_3$ and $b_3$ are positive constants.

Consequently, the intraguild predation model that we will study is

$\begin{array}{l}\stackrel{\dot{}}{x}=\rho x\left(1-x/R\right)-y\frac{x}{a_{1}x+b_1}-z\frac{x}{a_{2}x+b_2},\\ \stackrel{\dot{}}{y}=c_{1}y\frac{x}{a_{1}x+b_1}-c_{2}y-z\frac{y}{a_{3}y+b_3},\\ \stackrel{\dot{}}{z}=z\left(c_3\frac{y}{a_{3}y+b_3}+c_4\frac{x}{a_{2}x+b_2}-\mu \right).\end{array}$ (2)

For ecological considerations the domain of interest $\Omega$ is the positive octant of ${ℝ}^3,$ that is $\Omega :=\left\{\left(x,y,z\right)\in {ℝ}^3:x>0,y>0,z>0\right\}.$

We now state our main result. We establish the existence of a unique equi- librium point $p_0$ for the system (2) in $\Omega$ , at this point, $p_0$ , we show that the system exhibits a Hopf’s Bifurcation and the limit cycle given by the bifurcation is stable. All of this is obtained under certain restrictions on the parameters involved in the system.

Theorem 1 (Main result). If the positive parameters involved in system (2) satisfy the conditions

$\begin{array}{l}{a}_1=\frac{154k_1}{75{\rho }^2{x}_0^2},a_2=\frac{b_3}{x_0},a_3=\frac{b_3\rho x_0}{2k_1},b_1=\frac{2k_1}{{\rho }^2{x}_0},b_2=b_3,\\ c_1=\frac{74k_1}{\rho x_0^2},c_2=\frac{73\rho }{4},c_3=\frac{b_3\mu \rho x_0}{2k_1},c_4=\frac{b_3\mu }{x_0},R=2x_0,\end{array}$

where $k_1\mathrm{<mo>\; >\; </mo>0}$ and $x_0\mathrm{<mo>\; >\; </mo>0,}$ then the point $p_0=\left(x_0,\frac{2k_1}{x_0\rho },\frac{b_3\rho }{76}\right)$ is the unique

equilibrium point of system (2) in $\Omega$ moreover, we have a Hopf bifurcation in $p_0$ and the limit cycle that bifurcates from the equilibrium $p_0$ of system (2) as

$\mu$ increases from the critical value ${\mu }_0=\frac{47,775,075\rho }{600,704}$ is stable.

This article is organized as follows.

In Section 2 we provide the reader with the results that allowed us to study the system. In particular we present the version of the well known Hopf’s Bifurcation Theorem.

Section 3 is devoted to study the equilibrium points for our system in the positive octant with the aim of guaranteeing the hypothesis of Hopf’s Bifurcation Theorem. For this, we consider two subsections, the subsection 3.1 in which we show, under certain conditions on the parameters, the existence of an equi- librium point $p_0$ in the positive octant of ${ℝ}^3$ ; and subsection 3.2 where we show that under certain conditions on the parameters, the eigenvalues for the linear system at the equilibrium point $p_0$ associated to the system given in (2) are $\alpha$ which is real and $\pm i\omega$ the conjugated pure imaginary, and also the Lyapunov coefficient is computed.

In Section 4 we provide the proof of our main result in this paper. Further- more in Section 5 we provide the reader with a numerical result showing the stable limit cycle of the system.

In order to obtain all the calculations and simulations in this paper, we made use of a routine in the program Mathematica. This allowed us to simplify most of the process needed to obtain our result.

2. Preliminaries

2.1. Lyapunov Coefficient and Hopf Bifurcation

One of the main tools to determine the existence of a stable or unstable limit cycle is the first Lyapunov coefficient. This, in general, is not easy to calculate. To compute the first Lyapunov coefficient $l_1\left(p_0\right)$ of a differential system at an equilibrium point $p_0\mathrm{,}$ we make use of result by Kuznetsov (see [ [8] , p. 175]) whose statement is given in the following Theorem (cf. [4] ).

Theorem 2. Let $F\mathrm{<mo>\; :\; </mo>}\Omega \to {ℝ}^n$ be a differentiable map of class $C^4$ in an open subset $\Omega$ of ${ℝ}^n$ whose third order Taylor approximation of $F$ around $p_0=0$ is

$F\left(x\right)=Ax+\frac{1}{2!}B\left(x,x\right)+\frac{1}{3!}C\left(x,x,x\right)+O\left({\left|x\right|}^4\right).$

being $B$ and $C$ bilinear and trilinear forms, respectively. More over, assume that $A$ has a pair of purely imaginary eigenvalues $\pm \omega i$ . Let $q$ be the eigenvector of $A$ corresponding to the eigenvalue $\omega i$ , normalized so that the hermitian product satisfies $q\cdot \bar{q}=1$ being $\bar{q}$ the conjugate vector of $q$ . Let $p$ be the adjoint eigenvector such that $A^{\text{T}}p\mathrm{<mo>\; =\; </mo>}-\omega ip$ and $\bar{p}\cdot q=1$ . If $I$ denotes the identity matrix, then the first Lyapunov constant ${\mathcal{l}}_1\left(p_0\right)$ of the system of Ordinary Differential Equations $\stackrel{\dot{}}{x}=F\left(x\right)$ with an equilibrium point at $p_0=0$ is

$\frac{1}{2\omega }\mathrm{Re}\left(\bar{p}\cdot C\left(q,q,\bar{q}\right)-2\bar{p}\cdot B\left(q,A^{-1}B\left(q,\bar{q}\right)\right)+\bar{p}\cdot B\left(\bar{q},{\left(2\omega iI-A\right)}^{-1}B\left((q,q\right)\right)\right).$ (3)

The next theorem was proved by E. Hopf in 1942 (see [9] and for a proof in the bidimensional case see [10] and the general case see [ [11] , Section 5], and [ [8] , Section 5.4]). This theorem guarantees the existence of a Hopf’s bifurcation at an equilibrium point of a system of ordinary differential equations $\stackrel{\dot{}}{x}=F\left(x,\mu \right)$ whenever $\mu$ reaches a critical value ${\mu }_0$ .

Theorem 3 (Hopf's Theorem.). Suppose that the $C^4$ -system

$\stackrel{\dot{}}{x}=f\left(x,\mu \right),$ (4)

with $x\in {ℝ}^n$ and $\mu \in ℝ\mathrm{,}$ has a critical point $p_0$ for $\mu ={\mu }_0$ and that $Df\left(p_0\mathrm{,}{\mu }_0\right)$ has a simple pair of pure imaginary eigenvalues and no other eigenvalues with zero real part. Then there is a smooth curve of equilibrium points $p\left(\mu \right)$ with $p\left({\mu }_0\right)=p_0$ and the eigenvalues, $\lambda \left(\mu \right)$ and $\bar{\lambda }\left(\mu \right)$ of $Df\left(p\left(\mu \right)\mathrm{,}\mu \right)$ , which are pure imaginary at $\mu ={\mu }_0$ vary smoothly with $\mu$ . Furthermore, if

${\xi }_0:=\frac{\text{d}\mathrm{Re}\left(\lambda \left(\mu \right)\right)}{\text{d}\mu }\left({\mu }_0\right)\ne 0,$

then there is a unique two-dimensional center manifold passing through the point $\left(x_0\mathrm{,}{\mu }_0\right)$ and a smooth transformation of coordinates such that the sys- tem (4) on the center manifold is transformed into the normal form

$\begin{array}{l}\stackrel{\dot{}}{x}=-y+ax\left(x^2+y^2\right)-by\left(x^2+y^2\right)+O\left({\left|x\right|}^4\right)\\ \stackrel{\dot{}}{y}=x+bx\left(x^2+y^2\right)+ay\left(x^2+y^2\right)+O\left({\left|x\right|}^4\right)\end{array}$

in a neighborhood of the origin which, for $a\ne 0$ , has a weak focus of multiplicity one at the origin and

$\begin{array}{l}\stackrel{\dot{}}{x}=\mu x-y+ax\left(x^2+y^2\right)-by\left(x^2+y^2\right)\\ \stackrel{\dot{}}{y}=x+\mu y+bx\left(x^2+y^2\right)+ay\left(x^2+y^2\right)\end{array}$

is a universal unfolding of this normal form in a neighborhood of the origin on the center manifold. Moreover a periodic solution bifurcates from the point $p_0$ for $\mu >{\mu }_0$ if ${\mathcal{l}}_1\left(p_0\right){\xi }_0<0$ or for $\mu <{\mu }_0$ if ${\mathcal{l}}_1\left(p_0\right){\xi }_0>0$ . This periodic solution is stable if ${\mathcal{l}}_1\left(p_0\right)<0$ and unstable if ${\mathcal{l}}_1\left(p_0\right)>0$ . For ${\xi }_0<\mathrm{0,}$ the equilibrium point $p_0$ is a locally stable point for $\mu >{\mu }_0$ and locally unstable point for $\mu <{\mu }_0$ . For ${\xi }_0>\mathrm{0,}$ the equilibrium point $p_0$ is locally unstable point for $\mu >{\mu }_0$ and locally stable point for $\mu <{\mu }_0$ .

2.2. A Little of Linear Algebra

In this subsection we show a few results from Linear Algebra that allowed us to simplify our calculations in the next sections. This will provide us with a dif- ferent technique to find the eigenvalues of a given matrix.

If $M$ is a $3\times 3$ matrix with $m_{33}=0$ , its characteristic polynomial is de- termined by the entries of $M$ as a classical computation shows. In fact,

$\begin{array}{c}{P}_M\left(\lambda \right)=\det \left(M-\lambda Id\right)\\ =-{\lambda }^3+{\lambda }^2\left(m_{11}+m_{22}\right)+\lambda \left(-m_{11}{m}_{22}+m_{12}{m}_{21}+m_{13}{m}_{31}+m_{23}{m}_{32}\right)\\ -m_{11}{m}_{23}{m}_{32}+m_{12}{m}_{23}{m}_{31}+m_{13}{m}_{21}{m}_{32}-m_{13}{m}_{22}{m}_{31}.\end{array}$ (5)

Lemma 1. Let $M=\left(m_{ij}\right)$ be a $3\times 3$ matrix with $m_{33}=0$ . Then $M$ has eigenvalues $\alpha \in ℝ$ and $\pm i\omega$ with $\omega >0$ if and only if:

$\begin{array}{l}\alpha =m_{11}+m_{22}\\ {\omega }^2=m_{11}{m}_{22}-m_{12}{m}_{21}-m_{13}{m}_{31}-m_{23}{m}_{32}\\ \alpha {\omega }^2=-m_{11}{m}_{23}{m}_{32}+m_{12}{m}_{23}{m}_{31}+m_{13}{m}_{21}{m}_{32}-m_{13}{m}_{22}{m}_{31}\end{array}$ (6)

Proof. $M$ has the given eigenvalues of $\alpha \in ℝ$ and $\pm i\omega$ with $\omega >0$ , if and only if its characteristic polynomial takes the form:

$p\left(\lambda \right)=-\left(\lambda -\alpha \right)\left({\lambda }^2+{\omega }^2\right)=-{\lambda }^3+\alpha {\lambda }^2-\lambda {\omega }^2+\alpha {\omega }^2$

Comparing to (5) we obtain the result. $\square$

Corollary 1. If $M=\left(m_{ij}\right)$ is a $3\times 3$ matrix with $m_{33}=0$ . Then $M$ has eigenvalues $\alpha \in ℝ$ and $\pm i\omega$ with $\omega >0$ if and only if $D\left(M\right)=0$ , where

$\begin{array}{c}D\left(M\right)=-m_{11}{m}_{22}^2+m_{11}{m}_{12}{m}_{21}+m_{13}{m}_{31}{m}_{11}-m_{11}^2{m}_{22}\\ +m_{12}{m}_{21}{m}_{22}+m_{12}{m}_{23}{m}_{31}+m_{13}{m}_{21}{m}_{32}+m_{22}{m}_{23}{m}_{32}\mathrm{.}\end{array}$

Proof. Use that the system (6) is satisfied. $\square$

3. Equilibrium Points in the Positive Octant

3.1. Existence of an Equilibrium Point $p_0$ in the Positive Octant of ${ℝ}^3$

In order to find the equilibrium points and the restrictions in the parameters involved in the system (2) we use a different approach. We think of the equi- librium point $p_0=\left(x_0,y_0,z_0\right)$ as a new three parameters of the system. In this way our system (2) will have as new parameters the values of $x_0,y_0,z_0$ which we are setting to be positive. This approach differs from the usual method that is applied to find the equilibrium points (See [4] ).

In the next lemma we proceed to show the existence of an equilibrium point given conditions on the parameters involved in the system of differential equations. Moreover we can guarantee that the equilibrium point will be in $\Omega$ .

Lemma 2. Assume that the parameters in the system (2) are given by

$\begin{array}{l}{a}_1=\frac{b_3{k}_1\rho \left(a_2{x}_0+b_2\right)}{2\left(a_2{b}_1{\rho }^2{x}_0^2+b_3\left(2b_1{c}_2\rho x_0+k_1\right)\right)},\\ a_3=\frac{b_2{b}_3}{a_2{b}_1{x}_0\rho },c_3=\frac{b_2{b}_3\mu }{a_2{b}_1{x}_0\rho },c_4=a_2\mu ,R=2x_0,\end{array}$

where, $k_1>0$ and $x_0>0$ then

$p_0=\left(x_0,b_1\rho ,\frac{b_3{k}_1\rho \left(a_2{x}_0+b_2\right)}{2\left(a_2{b}_1{\rho }^2{x}_0{}^2+b_3\left(2b_1{c}_2\rho x_0+k_1\right)\right)}\right),$

is an equilibrium point of the system (2) in the region $\Omega \mathrm{.}$

Proof. The equilibrium points of the system are solutions of the following equations.

$\begin{array}{l}\rho x\left(1-\frac{x}{R}\right)-y\frac{x}{a_{1}x+b_1}-z\frac{x}{a_{2}x+b_2}=0,\\ c_{1}y\frac{x}{a_{1}x+b_1}-c_{2}y-z\frac{y}{a_{3}y+b_3}=0,\\ z\left(c_3\frac{y}{a_{3}y+b_3}+c_4\frac{x}{a_{2}x+b_2}-\mu \right)=0,\end{array}$ (7)

By multiplying the above equations by the denumerators (which are always non zero), involved in each corresponding equation we obtain that the equi- librium point must satisfy (8). Correspondingly each solution of (8) must also be an equilibrium point of the system (2)

$\begin{array}{l}R\left(b_{2}y+a_{2}xy+\left(b_1+a_{1}x\right)z\right)-\left(R-x\right)\left(b_1+a_{1}x\right)\left(b_2+a_{2}x\right)\rho =0,\\ \left(b_1{c}_2-c_{1}x+a_1{c}_{2}x\right)\left(b_3+a_{3}y\right)+\left(b_1+a_{1}x\right)z=0,\\ b_2\left(b_3\mu -c_{3}y+a_3\mu y\right)-x\left(b_3\left(c_4-a_2\mu \right)\right)+\left(a_2{c}_3+a_3{c}_4-a_2{a}_3\mu \right)y=0,\end{array}$ (8)

By taking $c_3=a_3\mu$ and $c_4=a_2\mu$ (8) reduces to

$\begin{array}{l}R\left(b_{2}y+a_{2}xy+\left(b_1+a_{1}x\right)z\right)-\left(R-x\right)\left(b_1+a_{1}x\right)\left(b_2+a_{2}x\right)\rho =0,\\ \left(b_1{c}_2-c_{1}x+a_1{c}_{2}x\right)\left(b_3+a_{3}y\right)+\left(b_1+a_{1}x\right)z=0,\\ b_2{b}_3\mu -x{a}_2{a}_3\mu y=0.\end{array}$ (9)

Notice that the last equation in the system above, is linear with respect to the variable $a_3$ . Solving this equation and substituting this value on the second one we obtain a system of two equations in $a_1$ and $z$ where the exponents of $a_1$ and $z$ in each equation is 1. From there we can obtain that the solutions of the system (9) are the following.

$\begin{array}{l}{a}_1=-\frac{b_1}{x}+\frac{R\left(b_3{c}_1+a_{2}y\right)}{b_3{c}_{2}R+a_2\left(R-x\right)x\rho }\\ a_3=\frac{b_2{b}_3}{a_{2}xy}\\ z=-\frac{b_3\left(b_2+a_{2}x\right)\left(c_{2}Ry+c_{1}x\left(-R+x\right)\rho \right)}{Rx\left(b_3{c}_1+a_{2}y\right)},\end{array}$

Taking $R=2x$ and $c_1=\frac{k_1+c_{2}Ry}{\left(R-x\right)x\rho },$ where $k_1>0,$ we have:

$\begin{array}{l}{a}_1=\left(a_2{x}^2\rho \left(2y-b_1\rho \right)+2b_3\left(k_1+c_{2}x\left(2y-b_1\rho \right)\right)\right)/\left(x^2\rho \left(2b_3{c}_2+a_{2}x\rho \right)\right),\\ a_3=\left(b_2{b}_3\right)/\left(a_{2}xy\right),\\ z=\left(b_3{k}_1\left(b_2+a_{2}x\right)\rho \right)/\left(2\left(b_3\left(k_1+2c_{2}xy\right)+a_2{x}^{2}y\rho \right)\right).\end{array}$

From there if $y=b_1\rho$ then all parameters involved in (2) becomes positive and

$p_0=\left(x_0,b_1\rho ,\frac{b_3{k}_1\rho \left(a_2{x}_0+b_2\right)}{2\left(a_2{b}_1{\rho }^2{x}_0^2+b_3\left(2b_1{c}_2\rho x_0+k_1\right)\right)}\right)$

is a solution of (7). Thus proving the lemma. $\square$

3.2. A Pair of Pure Imaginary Eigenvalues and the First Lyapunov Coefficient

Now our goal is to determine when the equilibrium point $p_0$ exhibits a Hopf's bifurcation. In order to show this, we show the existence of parameters where the equilibrium point has a pair of pure imaginary eigenvalues and a negative real eigenvalue. Making use of Hopf’s Theorem, we shall prove the existence of a Hopf bifurcation.

Theorem 4. If the parameters involved in system (2) satisfy the conditions of Lemma 2 and additionally

$c_2=\frac{73\rho }{4},b_1=\frac{2k_1}{{\rho }^2{x}_0},b_2=b_3,a_2=\frac{b_3}{x_0}\text{and}\mu =\frac{47,775,075\rho }{600,704},$

then the equilibrium point $p_0$ is

$p_0=\left(x_0,\frac{2k_1}{x_0\rho },\frac{b_3\rho }{76}\right)$

and the eigenvalues of the linear approximation of system (2) at $p_0$ are

$\alpha =-\frac{5,625\rho }{23,104}\text{and}\pm i\omega \mathrm{,}$ (10)

where

${\omega }^2=\frac{\rho \left(23,104\mu +31,213,049\rho \right)}{7,023,616}>0.$ (11)

Proof. The Jacobian matrix $M\left(p\right)=\left(m_{ij}\right)$ of the system (2) at $p_0$ is

$M=\left(\begin{array}{ccc}{m}_{11}& m_{12}& m_{13}\\ m_{21}& m_{22}& m_{23}\\ m_{31}& m_{32}& m_{33}\end{array}\right)$

where

$\begin{array}{l}{m}_{11}=\frac{1}{4}\rho \left(\frac{b_3{k}_1\left(a_2{x}_0\left(2a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+4b_1{c}_2{x}_0\rho \right)\right)-b_2{b}_3{k}_1\right)}{\left(b_2+a_2{x}_0\right){\left(a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+2b_1{c}_2{x}_0\rho \right)\right)}^2}-1\right)\\ m_{12}=-\frac{x_0^2\rho \left(2b_3{c}_2+a_2{x}_0\rho \right)}{2\left(a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+2b_1{c}_2{x}_0\rho \right)\right)}\\ m_{13}=-\frac{x_0}{b_2+a_2{x}_0}\\ m_{21}=\frac{b_1^2{\rho }^2{\left(2b_3{c}_2+a_2{x}_0\rho \right)}^2\left(k_1+2b_1{c}_2{x}_0\rho \right)}{4{\left(a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+2b_1{c}_2{x}_0\rho \right)\right)}^2}\\ m_{22}=\frac{a_2{b}_2{k}_1{x}_0\rho }{2\left(b_2+a_2{x}_0\right)\left(a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+2b_1{c}_2{x}_0\rho \right)\right)}\\ m_{23}=-\frac{a_2{b}_1{x}_0\rho }{b_2{b}_3+a_2{x}_0{b}_3}\\ m_{31}=\frac{a_2{b}_2{b}_3\mu k_1\rho }{2\left(b_2+a_2{x}_0\right)\left(a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+2b_1{c}_2{x}_0\rho \right)\right)}\\ m_{32}=\frac{a_2{b}_2{b}_3\mu k_1{x}_0}{2b_1\left(b_2+a_2{x}_0\right)\left(a_2{b}_1{x}_0^2{\rho }^2+b_3\left(k_1+2b_1{c}_2{x}_0\rho \right)\right)}\\ m_{33}=0.\end{array}$

Using Corollary 1 the characteristic polynomial of $M$ has roots $\alpha =m_{11}+m_{22}$ and $\pm i\omega$ where

${\omega }^2=m_{11}{m}_{22}-m_{12}{m}_{21}-m_{13}{m}_{31}-m_{23}{m}_{32}$

if and only if

$\begin{array}{c}0=-m_{11}{m}_{22}^2+m_{11}{m}_{12}{m}_{21}+m_{13}{m}_{31}{m}_{11}-m_{11}^2{m}_{22}+m_{12}{m}_{21}{m}_{22}\\ +m_{12}{m}_{23}{m}_{31}+m_{13}{m}_{21}{m}_{32}+m_{22}{m}_{23}{m}_{32}\mathrm{.}\end{array}$ (12)

In this case the value of $\alpha$ is given by:

$\begin{array}{c}\alpha =-[\rho (2b_2{b}_3^2{k}_1^2-2a_2{b}_2{b}_3{k}_1^2{x}_0+4b_1{b}_2{b}_3^2{c}_2{k}_1{x}_0\rho -4a_2{b}_1{b}_2{b}_3{c}_2{k}_1{x}_0^2\rho \\ +4b_1^2{b}_2{b}_3^2{c}_2^2{x}_0^2{\rho }^2+2a_2{b}_1{b}_2{b}_3{k}_1{x}_0^2{\rho }^2+4a_2{b}_1^2{b}_3^2{c}_2^2{x}_0^3{\rho }^2-2a_2^2{b}_1{b}_2{k}_1{x}_0^3{\rho }^2\\ +4a_2{b}_1^2{b}_2{b}_3{c}_2{x}_0^3{\rho }^3+4a_2^2{b}_1^2{b}_3{c}_2{x}_0^4{\rho }^3+a_2^2{b}_1^2{b}_2{x}_0^4{\rho }^4\\ +a_2^3{b}_1^2{x}_0^5{\rho }^4)]/\left[4\left(b_2+a_2{x}_0\right){\left(b_3{k}_1+2b_1{b}_3{c}_2{x}_0\rho +a_2{b}_1{x}_0^2{\rho }^2\right)}^2\right]\end{array}$ (13)

Choose $k_2>0$ and

$c_2=-\frac{k_2+2a_2{b}_2{b}_3{k}_1^2{x}_0+2a_2^2{b}_1{b}_2{k}_1{x}_0^3{\rho }^2}{4a_2{b}_1{b}_2{b}_3{x}_0^2\rho \left(k_1-b_1{x}_0{\rho }^2\right)},$

Now taking $b_1=\left(2k_1\right)/\left(x_0{\rho }^2\right)$ we have $c_2>0$ and with this choices $\alpha <0$ .

Now taking $b_2=b_3$ , the expression for ${\omega }^2$ simplifies to:

$\begin{array}{c}{\omega }^2=[\rho [32a_2^3{b}_3^4\mu k_1^4{x}_0^3\left(b_3+a_2{x}_0\right){\left(k_2+4a_2{b}_3{k}_1^2{x}_0\left(b_3+2a_2{x}_0\right)\right)}^2\\ +(b_3^2{k}_2^4+2a_2{b}_3{k}_2^3\left(5b_3^3{k}_1^2+k_2\right)x_0+a_2^2{k}_2^2\left(36b_3^6{k}_1^4+48b_3^3{k}_1^2{k}_2+k_2^2\right)x_0^2\\ +2a_2^3{b}_3^2{k}_1^2{k}_2\left(28b_3^6{k}_1^4+140b_3^3{k}_1^2{k}_2+33k_2^2\right)x_0^3\\ +4a_2^4{b}_3{k}_1^2\left(8b_3^9{k}_1^6+144b_3^6{k}_1^4{k}_2+187b_3^3{k}_1^2{k}_2^2+7k_2^3\right)x_0^4\\ +8a_2^5{b}_3^3{k}_1^4\left(40b_3^6{k}_1^4+311b_3^3{k}_1^2{k}_2+99k_2^2\right)x_0^5\\ +16a_2^6{b}_3^2{k}_1^4\left(140b_3^6{k}_1^4+297b_3^3{k}_1^2{k}_2+18k_2^2\right)x_0^6+32a_2^7{b}_3^4{k}_1^6\left(221b_3^3{k}_1^2+128k_2\right)x_0^7\\ +128a_2^8{b}_3^3{k}_1^6\left(83b_3^3{k}_1^2+10k_2\right)x_0^8+7680a_2^9{b}_3^5{k}_1^8{x}_0^9\\ +2048a_2^{10}{b}_3^4{k}_1^8{x}_0^{10})\rho ]]/\left[32a_2{b}_3^2{k}_1^2{x}_0{\left(b_3+a_2{x}_0\right)}^2{\left(k_2+4a_2{b}_3{k}_1^2{x}_0\left(b_3+2a_2{x}_0\right)\right)}^3\right].\end{array}$ (14)

Now solving Equation (12) for the parameter $\mu$ in terms of $a_2,b_3,k_1,x_0,k_2$ and $\rho \mathrm{,}$ we obtain:

$\begin{array}{c}\mu =-[(b_3{k}_2^2+8a_2{b}_3^3{k}_1^2{k}_2{x}_0+a_2{k}_2^2{x}_0+20a_2^2{b}_3^5{k}_1^4{x}_0^2+16a_2^2{b}_3^2{k}_1^2{k}_2{x}_0^2\\ +52a_2^3{b}_3^4{k}_1^4{x}_0^3+16a_2^3{b}_3{k}_1^2{k}_2{x}_0^3+64a_2^4{b}_3^3{k}_1^4{x}_0^4+64a_2^5{b}_3^2{k}_1^4{x}_0^5)\\ (b_3^2{k}_2^4+10a_2{b}_3^4{k}_1^2{k}_2^3{x}_0+2a_2{b}_3{k}_2^4{x}_0+36a_2^2{b}_3^6{k}_1^4{k}_2^2{x}_0^2+48a_2^2{b}_3^3{k}_1^2{k}_2^3{x}_0^2\\ +a_2^2{k}_2^4{x}_0^2+56a_2^3{b}_3^8{k}_1^6{k}_2{x}_0^3+280a_2^3{b}_3^5{k}_1^4{k}_2^2{x}_0^3+66a_2^3{b}_3^2{k}_1^2{k}_2^3{x}_0^3\\ +32a_2^4{b}_3^{10}{k}_1^8{x}_0^4+576a_2^4{b}_3^7{k}_1^6{k}_2{x}_0^4+748a_2^4{b}_3^4{k}_1^4{k}_2^2{x}_0^4+28a_2^4{b}_3{k}_1^2{k}_2^3{x}_0^4\\ +320a_2^5{b}_3^9{k}_1^8{x}_0^5+2488a_2^5{b}_3^6{k}_1^6{k}_2{x}_0^5+792a_2^5{b}_3^3{k}_1^4{k}_2^2{x}_0^5+2240a_2^6{b}_3^8{k}_1^8{x}_0^6\\ +4752a_2^6{b}_3^5{k}_1^6{k}_2{x}_0^6+288a_2^6{b}_3^2{k}_1^4{k}_2^2{x}_0^6+7072a_2^7{b}_3^7{k}_1^8{x}_0^7+4096a_2^7{b}_3^4{k}_1^6{k}_2{x}_0^7\\ +10624a_2^8{b}_3^6{k}_1^8{x}_0^8+1280a_2^8{b}_3^3{k}_1^6{k}_2{x}_0^8+7680a_2^9{b}_3^5{k}_1^8{x}_0^9\\ +2048a_2^{10}{b}_3^4{k}_1^8{x}_0^{10})\rho ]/[8a_2^2{b}_3^3{k}_1^2{x}_0^2\left(b_3+a_2{x}_0\right){\left(k_2+4a_2{b}_3^2{k}_1^2{x}_0+8a_2^2{b}_3{k}_1^2{x}_0^2\right)}^2\\ (-k_2^3-4a_2{b}_3^2{k}_1^2{k}_2^2{x}_0+12a_2^2{b}_3^4{k}_1^4{k}_2{x}_0^2-12a_2^2{b}_3{k}_1^2{k}_2^2{x}_0^2+64a_2^3{b}_3^6{k}_1^6{x}_0^3\\ -16a_2^3{b}_3^3{k}_1^4{k}_2{x}_0^3+112a_2^4{b}_3^5{k}_1^6{x}_0^4+128a_2^5{b}_3^4{k}_1^6{x}_0^5+256a_2^6{b}_3^3{k}_1^6{x}_0^6)]\end{array}$ (15)

Choosing

$k_2=\frac{4a_2{k}_1^2{x}_0\left(4b_3^3+7a_2{b}_3^2{x}_0+8a_2^2{b}_3{x}_0^2+16a_2^3{x}_0^3\right)}{-3b_3+4a_2{x}_0},$

$\mu >0,$ and if $a_2=\left(3b_3+k_5\right)/\left(4x_0\right),$ with $k_5>0,$ we also obtain that, $k_2>0$ . If we take $k_5=b_3$ then $k_2=140b_3^3{k}_1^2,$ hence all expressions of the as- signed parameters of system (2) are simplified:

$\begin{array}{c}{a}_1=\frac{154k_1}{75{\rho }^2{x}_0{}^2},a_2=\frac{b_3}{x_0},a_3=\frac{b_3\rho x_0}{2k_1},b_1=\frac{2k_1}{{\rho }^2{x}_0},b_2=b_3,\\ c_1=\frac{74k_1}{\rho x_0{}^2},c_2=\frac{73\rho }{4},c_3=\frac{b_3\mu \rho x_0}{2k_1},c_4=\frac{b_3\mu }{x_0},R=2x_0,\end{array}$

and the expression for $\mu$ given by (15) simplifies to $\mu =\frac{47,775,075\rho }{600,704}.$ Thus

the equilibrium point is

$p_0=\left(x_0,\frac{2k_1}{x_0\rho },\frac{b_3\rho }{76}\right),$

and from (13) and (14) the eigenvalues of the linear approximation of system (2) at $p_0$ are given by Equations (10) and (11), which proves the theorem. □

Remark 5. Notice that by Theorem 4 and Subsection 2.2, the characteristic polynomial of the linear approximation of system (2) at the equilibrium point $p_0$ has the form $P_M\left(\lambda \right)=-\left(\lambda -\alpha \right)\left(\lambda +{\omega }^2\right)$ if and only if the fundamental Equation (12) is satisfied, which, in this case, Equation (12) reduces to

$D\left(M\right)=\frac{3,675{\rho }^2\left(47,775,075\rho -600,704\mu \right)}{162,273,624,064}=\mathrm{0,}$

thus, the linear approximation of system (2) at the equilibrium point $p_0$ has a

pair of pure imaginary eigenvalues if and only if $\mu =\frac{47,775,075\rho }{600,704}$ .

Applying the Theorem 2 to system (2) at the equilibrium point $p_0$ we get the following result.

Theorem 6. If the parameters involved in system (2) satisfy the hypothesis of Lemma 2 and Theorem 4 then the eigenvalues of the linear approximation of system (2) at the equilibrium point

$p_0=\left(x_0,\frac{2k_1}{x_0\rho },\frac{b_3\rho }{76}\right)$ (16)

are $\alpha =-\frac{5,625\rho }{23,104}$ and $\pm i\omega ,$ where ${\omega }^2=\frac{45,227,071{\rho }^2}{9,611,264},$ and the first

Lyapunov coefficient ${\mathcal{l}}_1\left(p_0\right)$ of the differential system (2) at the equilibrium point $p_0$ is given by

${\mathcal{l}}_1\left(p_0\right)=-\frac{s_1{x}_0^2{\rho }^2}{s_2\left(s_3{k}_1^2+s_4{x}_0^4{\rho }^2+s_5{b}_3^2{x}_0^2{\rho }^4\right)},$

where $s_1,s_2,s_3,s_4$ y $s_5$ are the positive constants defined by

$\begin{array}{c}{s}_1=2,247,304,865,615,302,723,837,224,164,005,159,392,947,216\sqrt{\frac{26}{45,227,071}},\\ s_2=1,538,018,750,999,306,870,027,605,\\ s_3=664,087,191,252,992,\\ s_4=9,099,199,882,240,\\ s_5=2,552,956,682,775.\end{array}$

Proof. Let $p_0$ as in (16) and consider a linear change of variables to translate $p_0$ to the origin of coordinates, after that change, we obtain a differential sys- tem

$\stackrel{\dot{}}{X}=F\left(X\right)$ (17)

with $X=\left(x,y,z\right)$ and $F\left(0\right)=0.$ Denote the vector field $F\left(x,y,z\right)=\left(F_1\left(x,y,z\right),F_2\left(x,y,z\right),F_3\left(x,y,z\right)\right),$ associated to this differential system. Now, we compute the linear part $A\mathrm{,}$ the bilinear $B\mathrm{,}$ and trilinear $C$ forms of the Taylor expansion of the function $F\mathrm{.}$

The linear part of system (17) at 0 is

$A=\left(\begin{array}{ccc}-\frac{5,701\rho }{23,104}& -\frac{75x_0^2{\rho }^2}{304k_1}& -\frac{x_0}{2b_3}\\ \frac{208,125k_1}{11,552x_0^2}& \frac{\rho }{304}& -\frac{k_1}{b_3{x}_0\rho }\\ \frac{47,775,075b_3{\rho }^2}{182,614,016x_0}& \frac{47,775,075b_3{x}_0{\rho }^3}{365,228,032k_1}& 0\end{array}\right)$

It follows immediately from Theorem 4 that the eigenvalues of $A$ are

$\alpha =-\frac{5,625\rho }{23,104}$ and $\pm i\omega \mathrm{,}$ where ${\omega }^2=\frac{45,227,071{\rho }^2}{9,611,264}.$

The bilinear function $B$ at vectors $\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\left(u\mathrm{,}v\mathrm{,}w\right)$ is given by

$B\left(\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\left(u\mathrm{,}v\mathrm{,}w\right)\right)=\left(B_1\mathrm{,}{B}_2\mathrm{,}{B}_3\right)$

where

$\begin{array}{c}{B}_1=-\frac{uz+wx}{4b_3}-\frac{3\rho \left(439,001k_{1}ux+71,250\rho x_0^2\left(uy+vx\right)\right)}{1,755,904k_1{x}_0},\\ B_2=\frac{-\frac{219,488\left(vz+wy\right)}{b_3}-\frac{16,025,625k_{1}ux}{x_0{}^3}+\frac{1,444{\rho }^{2}v{x}_{0}y}{k_1}+\frac{7,908,750\rho \left(uy+vx\right)}{x_0}}{877,952}\mathrm{,}\\ B_3=47,775,075\rho [4k_1{}^2\left(76x_0\left(uz+wx\right)-b_3\rho ux\right)-b_3{\rho }^{3}v{x}_0{}^{4}y\\ +152k_1\rho x_0{}^3\left(vz+wy\right)]/\left(730,456,064k_1^2{x}_0^2\right)\mathrm{.}\end{array}$

The trilinear function $C$ at vectors $\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\left(u\mathrm{,}v\mathrm{,}w\right)\mathrm{,}\left(r\mathrm{,}s\mathrm{,}t\right)$ is given by

$C\left(\left(x\mathrm{,}y\mathrm{,}z\right)\mathrm{,}\left(u\mathrm{,}v\mathrm{,}w\right)\mathrm{,}\left(r\mathrm{,}s\mathrm{,}t\right)\right)=\left(C_1\mathrm{,}{C}_2\mathrm{,}{C}_3\right)$

where

$\begin{array}{c}{C}_1=\frac{ruz+rwx+tux}{4b_3{x}_0}+\frac{3\rho \left(10,972,500\rho x_0^2\left(ruy+rvx+sux\right)-3,378,9601k_{1}rux\right)}{266,897,408k_1{x}_0^2},\\ C_2=\frac{\rho x_0\left(svz+swy+tvy\right)}{8b_3{k}_1}-\frac{3{\rho }^{3}sv{x}_0^{2}y}{2,432k_1^2}+\frac{3,701,919,375k_{1}rux}{133,448,704x_0^4}\\ -\frac{16,025,625\rho \left(ruy+rvx+sux\right)}{1,755,904x_0^2},\\ C_3=(47,775,075\rho (8k_1^3\left(3b_3\rho rux-152x_0\left(ruz+rwx+tux\right)\right)+3b_3{\rho }^{4}sv{x}_0^{6}y\\ -304k_1{\rho }^2{x}_0^5\left(svz+swy+tvy\right)))/\left[2,921,824,256k_1^3{x}_0^3\right]\end{array}$

The normalized eigenvector $q=\left(q_1,q_2,q_3\right)$ of $A$ corresponding to eigen- value $i\omega$ has coordinates

$\begin{array}{c}{q}_1=\frac{\left(-\frac{7,358,624}{4,007,775}+\frac{32,703,712i\sqrt{\frac{1,846}{637,001}}}{4,007,775}\right)x_0}{b_3\rho },\\ q_2=\frac{23,104\left(405,769,637+50,606i\sqrt{1,175,903,846}\right)k_1}{2,552,956,682,775b_3{\rho }^2{x}_0},\\ q_3=1.\end{array}$

The adjoint eigenvector $p=\left(p_1,p_2,p_3\right)$ of the transpose matrix of $A$ corresponding to the eigenvalue $-i\omega$ has coordinates

$\begin{array}{c}{p}_1=\frac{\left(-33,124,052+4,249i\sqrt{1,175,903,846}\right)b_3\rho }{36,342,592x_0},\\ p_2=\frac{\left(33,124,052+349i\sqrt{1,175,903,846}\right)b_3{\rho }^2{x}_0}{72,685,184k_1},\\ p_3=1.\end{array}$

Taking into account the formula of the first Lyapunov constant ${\mathcal{l}}_1\left(p_0\right)$ of Theorem 2, the values of $A\mathrm{,}B\mathrm{,}C$ and $q,\bar{q},\bar{p},$ we have that the expression of the first Lyapunov coefficient at the equilibrium point $p_0$ is

${\mathcal{l}}_1\left(p_0\right)=-\frac{s_1{x}_0^2{\rho }^2}{s_2\left(s_3{k}_1^2+s_4{x}_0^4{\rho }^2+s_5{b}_3^2{x}_0^2{\rho }^4\right)},$

where $s_1\mathrm{,}{s}_2\mathrm{,}{s}_3\mathrm{,}{s}_4$ y $s_5$ are the positive constants defined by

$\begin{array}{c}{s}_1=224,730,486,561,530,272,383,722,416,400,515,939,2947,216\sqrt{\frac{26}{45,227,071}},\\ s_2=1,538,018,750,999,306,870,027,605,\\ s_3=664,087,191,252,992,\\ s_4=9,099,199,882,240,\\ s_5=2,552,956,682,775.\end{array}$

$\square$

Remark 7. Notice that with the parameters as in Theroem 4 and Theorem 6 and according with the above result the first Lyapunov coefficient of the system at the given equilibrium point is always negative.

4. Proof of the Main Result

In this section, using the results given in Section 2 and results obtained in Subsection 3.2, we give a proof of our main result given by Theorem 1.

Proof of Theorem 1. If

$a_1=\frac{154k_1}{75{\rho }^2{x}_0{}^2},a_3=\frac{b_3\rho x_0}{2k_1},c_3=\frac{b_3\mu \rho x_0}{2k_1},c_4=\frac{b_3\mu }{x_0},R=2x_0,$

it follows immediately from Lemma (2), that $p_0=\left(x_0,\frac{2k_1}{x_0\rho },\frac{b_3\rho }{76}\right)$ is an equi-

librium point of system (2) in the positive octant of ${ℝ}^3\mathrm{.}$

If

$a_2=\frac{b_3}{x_0},b_1=\frac{2k_1}{{\rho }^2{x}_0},b_2=b_3,c_2=\frac{73\rho }{4},$

and taking into account the above assignments of $a_1\mathrm{,}{a}_3\mathrm{,}{c}_3\mathrm{,}{c}_4$ and $R\mathrm{,}$ it

follows from Theorem 4 and Remark 5 that if $\mu ={\mu }_0=\frac{47,775,075\rho }{600,704}$ the eigen-

values of the linear approximation of system (2) at the equilibrium point $p_0\mathrm{,}$

are $\alpha \left({\mu }_0\right)=-\frac{5,625\rho }{23,104}$ and the pure imaginary complex numbers $\pm i\omega \left({\mu }_0\right)\mathrm{,}$ where $\omega \left({\mu }_0\right)=\frac{1}{608}\sqrt{\frac{45,227,071}{26}}\rho \mathrm{.}$ For $\mu \in ℝ\mathrm{,}$ with the above assignments

of $a_1,a_2,a_3,b_1,b_2,c_3,c_4$ and $R\mathrm{,}$ the real part of the complex eigenvalues $\lambda \left(\mu \right)$ and $\bar{\lambda }\left(\mu \right)$ of the linear approximation of system (2) at the equilibrium point $p_0\mathrm{,}$ (the equilibrium point $p_0$ does not depend of the parameter $\mu$ ), is

$\mathrm{Re}\left(\lambda \left(\mu \right)\right)=\frac{\rho \left(-1,755,904\sqrt{3}\mu -2,361,644,849\sqrt{3}\rho +\sqrt{3}{r}_1^{1/3}-11,250r_2^{1/6}\right)}{138,624r_3^{1/6}},$

where,

$\begin{array}{c}{r}_1=5,413,801,206,181,003,264{\mu }^3+432,269,257,343,222,802,481,152{\mu }^2\rho \\ +13,437,254,340,777,791,771,358,085,799{\rho }^3-505,553,051,317,500\sqrt{3}\\ [{\rho }^3(3,749,169,810,374,656{\mu }^3+157,241,523,352,487,313,408{\mu }^2\rho \\ +13,117,228,409,565,237,258,048\mu {\rho }^2\\ {+9,213,647,144,507,986,638,030,671{\rho }^3)]}^{1/2}\\ +5,267,712\mu (1,614,096,551,032,357,801{\rho }^2+3,773,400\sqrt{3}\\ [\rho (374,916,981,037,4656{\mu }^3+157,241,523,352,487,313,408{\mu }^2\rho \\ +13,117,228,409,565,237,258,048\mu {\rho }^2{+9,213,647,144,507,986,638,030,671{\rho }^3)]}^{1/2})\mathrm{,}\end{array}$

$\begin{array}{c}{r}_2={\rho }^3(5,413,801,206,181,003,264{\mu }^3+432,269,257,343,222,802,481,152{\mu }^2\rho \\ +13,437,254,340,777,791,771,358,085,799{\rho }^3-505,553,051,317,500\sqrt{3}\\ [{\rho }^3(3,749,169,810,374,656{\mu }^3+157,241,523,352,487,313,408{\mu }^2\rho \\ +13,117,228,409,565,237,258,048\mu {\rho }^2\\ {+9,213,647,144,507,986,638,030,671{\rho }^3)]}^{1/2}\\ +5,267,712\mu (1,614,096,551,032,357,801{\rho }^2+3,773,400\sqrt{3}\\ [\rho (3,749,169,810,374,656{\mu }^3+157,241,523,352,487,313,408{\mu }^2\rho \\ \text{+}13,117,228,409,565,237,258,048\mu {\rho }^2\\ {+9,213,647,144,507,986,638,030,671{\rho }^3)]}^{1/2}))\mathrm{,}\end{array}$

and

$\begin{array}{c}{r}_3={\rho }^3(5,413,801,206,181,003,264{\mu }^3+432,269,257,343,222,802,481,152{\mu }^2\rho \\ +13,437,254,340,777,791,771,358,085,799{\rho }^3-50,555,3051,317,500\sqrt{3}\\ [{\rho }^3(3,749,169,810,374,656{\mu }^3+157,241,523,352,487,313,408{\mu }^2\rho \\ +13,117,228,409,565,237,258,048\mu {\rho }^2\\ {+9,213,647,144,507,986,638,030,671{\rho }^3)]}^{1/2}\\ +5,267,712\mu (1,614,096,551,032,357,801{\rho }^2+3,773,400\sqrt{3}\\ [\rho (3,749,169,810,374,656{\mu }^3+157,241,523,352,487,313,408{\mu }^2\rho \\ \text{+}13,117,228,409,565,237,258,048\mu {\rho }^2\\ {+9,213,647,144,507,986,638,030,671{\rho }^3)]}^{1/2}))\mathrm{.}\end{array}$

Hence,

${\xi }_0=\frac{\text{d}\mathrm{Re}\left(\lambda \left(\mu \right)\right)}{\text{d}\mu }\left({\mu }_0\right)=\frac{47,201,700}{33,065,273,387}>0.$

Moreover, by Theorem 6 the first Lyapunov coefficient of the differential sys- tem (2) at the equilibrium point $p_0\mathrm{,}$ is $l_1\left(p_0\right)<0,$ then applying Hopf’s Theorem, (Theorem 3), we have a Hopf’s bifurcation at $p_0$ and that the limit cycle that bifurcates from the equilibrium $p_0$ of system (2) as $\mu$ increases

from the critical value ${\mu }_0=\frac{47,775,075\rho }{600,704}$ is stable.

Now, taking into account the assignments for parameters given above, the system (2) has the form:

$\begin{array}{c}\stackrel{\dot{}}{x}=x\left(-\frac{x_{0}z}{b_{3}x+b_3{x}_0}-\frac{75{\rho }^2{x}_0{}^{2}y}{154k_{1}x+150k_1{x}_0}+\rho -\frac{\rho x}{2x_0}\right),\\ \stackrel{\dot{}}{y}=y\left(\rho \left(\frac{2775x}{77x+75x_0}-\frac{73}{4}\right)-\frac{2k_{1}z}{2b_3{k}_1+b_3\rho x_{0}y}\right),\\ \stackrel{\dot{}}{z}=\mu z\left(\frac{x}{x+x_0}-\frac{2k_1}{2k_1+\rho x_{0}y}\right),\end{array}$ (18)

which is in terms of the free positive parameters $x_0,b_3,k_1,\rho$ and $\mu \mathrm{.}$ We have that the rest of equilibrium points of system (18) are

$p_1=\left(\frac{5,475x_0}{5479},\frac{60,861,300k_1}{30,019,441\rho x_0},0\right),p_2=\left(2x_0,0,0\right),\text{and}{p}_3=\left(0,0,0\right),$ (19)

thus, the unique equilibrium point of system (18) in the positive octant of ${ℝ}^3$ is $p_0$ and the theorem is proved. $\square$

Remark 8. Notice that system (18) has, additionally to $p_0$ , the equilibrium points given by (19). For $\mu ={\mu }_0$ the eigenvalues of the linear approximation of the system are:

For $p_1$

$\frac{101,822,611,784,888,325\rho }{397,768,345,977,357,056}\mathrm{,}-\frac{\left(29,952,557\pm \sqrt{265,848,086,191,023,551}i\right)\rho }{243,267,600}\mathrm{.}$

For $p_2$

$-\frac{15,925,025\rho }{600,704},-\rho ,\frac{5,483\rho }{916}.$

For $p_3$

$-\frac{47,775,075\rho }{600,704}\mathrm{,}-\frac{1}{4}\left(73\rho \right)\mathrm{,}\rho \mathrm{.}$

As a consequence these equilibrium points are hyperbolic, moreover they are saddle points.

5. Numerical Result

Theorem 1 guarantees the existence of a Hopf’s bifurcation if we have the fol- lowing assignments for the parameters of system (2):

$\begin{array}{c}{a}_1=\frac{154k_1}{75{\rho }^2{x}_0{}^2},a_2=\frac{b_3}{x_0},a_3=\frac{b_3\rho x_0}{2k_1},b_1=\frac{2k_1}{{\rho }^2{x}_0},b_2=b_3,\\ c_1=\frac{74k_1}{\rho x_0{}^2},c_2=\frac{73\rho }{4},c_3=\frac{b_3\mu \rho x_0}{2k_1},c_4=\frac{b_3\mu }{x_0},R=2x_0.\end{array}$

With these assignments of the parameters the system (2) is in terms of the free positive parameters $x_0\mathrm{,}{b}_3\mathrm{,}{k}_1\mathrm{,}\rho$ and $\mu \mathrm{,}$ the unique equilibrium point of

system (2) in the positive octant of ${ℝ}^3$ is $p_0=\left(x_0,\frac{2k_1}{x_0\rho },\frac{b_3\rho }{76}\right).$ By Theorem 1, for $\mu$ close enough to ${\mu }_0=\frac{47,775,075\rho }{600,704}$ and $\mu >{\mu }_0$ then a limit cycle

bifurcates from the equilibrium $p_0$ of system (2).

For example, if we consider the parameters values

$\begin{array}{l}{x}_0=1,b_3=76,\\ k_1=1,\rho =1,\end{array}$ (20)

then the linear approximation of system (2) at $p_0=\left(\mathrm{1,2,1}\right)$ is

$DF\left(p_0\right)=\left(\begin{array}{ccc}-\frac{5,701}{23,104}& -\frac{75}{304}& -\frac{1}{152}\\ \frac{208,125}{11,552}& \frac{1}{304}& -\frac{1}{76}\\ \frac{\mu }{4}& \frac{\mu }{8}& 0\end{array}\right).$

The real part of the complex eigenvalues is

$\mathrm{Re}\left(\lambda \left(\mu \right)\right)=\frac{\sqrt{3}\sqrt[3]{S_2}-11,250\sqrt[6]{S_2}-1,755,904\sqrt{3}\mu -2,361,644,849\sqrt{3}}{138,624\sqrt[6]{S_2}},$

where

$\begin{array}{l}{S}_1=23,104\mu \left(92,416\mu \left(1,755,904\mu +73,643,242,047\right)+567,747,074,513,730,837\right)\hfill \\ +9,213,647,144,507,986,638,030,671,\hfill \\ S_2=-505,553,051,317,500\sqrt{3}\sqrt{S_1}+1,755,904\mu (1,755,904\mu (1,755,904\mu \hfill \\ +140,201,549,547)+11,320,200\sqrt{3}\sqrt{S_1}+4,842,289,653,097,073,403)\hfill \\ +13,437,254,340,777,791,771,358,085,799\hfill \end{array}$

and its derivative is

$\begin{array}{l}\begin{array}{l}\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}\mu }\left(\mu \right)\\ =-19(521,660,579,269,757,164,085,418,393,600{\mu }^4\end{array}\hfill \\ +284,936,905,588,473,856{\mu }^3{S}_5-162,273,624,064{\mu }^2{S}_6\hfill \\ -3,228,193,102,064,715,602\sqrt{3}\sqrt{S_3}\sqrt[3]{S_4}\hfill \\ -204,823,995,417,137,033,019,192,719,132,775\sqrt[3]{S_4}\hfill \\ +46,125,171,752,042,700,219,319,109,098\sqrt{3}\sqrt{S_3}-92,416\mu S_7\hfill \\ -539,617,382,860,506,682,442,178,492,636,645,334,335,975)/\left(6\sqrt{S_3}{S}_4^{7/6}\right)\mathrm{,}\hfill \end{array}$

where

$\begin{array}{l}{S}_3=23,104\mu (92,416\mu \left(1,755,904\mu +73,643,242,047\right)\hfill \\ +567,747,074,513,730,837)+9,213,647,144,507,986,638,030,\mathrm{671,}\hfill \\ S_4=-505,553,051,317,500\sqrt{3}\sqrt{S_3}+1,755,904\mu (1,755,904\mu (1,755,904\mu \hfill \\ +140,201,549,547)+11,320,200\sqrt{3}\sqrt{S_3}+4,842,289,653,097,073,403)\hfill \\ +13,437,254,340,777,791,771,358,085,\mathrm{799,}\hfill \\ S_5=-744,750\sqrt[3]{S_4}+38\sqrt{3}\sqrt{S_3}+85,934,584,097,414,\mathrm{325,}\hfill \\ S_6=38\sqrt{3}\sqrt{S_3}\sqrt[3]{S_4}+43,856,687,570,103,225\sqrt[3]{S_4}-7,013,802,672,786\sqrt{3}\sqrt{S_3}\hfill \\ +94,011,905,468,216,518,536,340,\mathrm{125,}\hfill \\ S_7=3,551,772,588,524\sqrt{3}\sqrt{S_3}\sqrt[3]{S_4}+3,840,506,095,057,118,582,941,800\sqrt[3]{S_4}\hfill \\ +8,081,347,093,810,953,130,686\sqrt{3}\sqrt{S_3}\hfill \\ -9,957,828,728,837,768,681,342,262,239,204,175.\hfill \end{array}$

If $\mu ={\mu }_0:=\frac{47,775,075}{600,704},$ then $DF\left(p_0\right)$ has eigenvalues

$-\frac{5,625}{23,104},\pm i\frac{1}{608}\sqrt{\frac{45,227,071}{26}},$

and $\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}\mu }\left({\mu }_0\right)=\frac{47,201,700}{33,065,273,387}.$

The Lyapunov coefficient is

${\mathcal{l}}_1\left(p_0\right)=-\frac{389,076,327,149,463,767,977,358,754,156,017,900,441\sqrt{\frac{26}{45,227,071}}}{4,105,749,627,489,518,890,132,734,733,438,379,235}<0,$

hence, we have a supercritical Hopf bifurcation, and then the periodic orbit obtained from the bifurcation is stable.

In Figure 1, we exhibit the stable limit cycle of differential system (2) with the above parameters values, that is, we show an orbit tending to the local attractor defined by a stable Hopf periodic orbit with $\mu ={\mu }_0+1/100$ and initial con-

dition $q_0=p_0+\left(\frac{1}{50},\frac{1}{50},\frac{1}{50}\right)$ .

Figure 2 shows the same behavior but with different initial condition $q_0=\left(1-4/{10}^4,2,0.025\right)$ .

Figure 1. Stable limit cycle and Time series with initial condition $q_0=p_0+\left(\frac{1}{50},\frac{1}{50},\frac{1}{50}\right)$ .

Figure 2. Stable limit cycle and Time series with initial condition $q_0=\left(1-4/{10}^4,2,0.025\right)$ .

Finally, notice that, under the assignations in (20) one has the following:

・ $f_1\left(x\right)=\frac{x}{\frac{154x}{75}+2}>\frac{x}{76x+76}=f_2\left(x\right)$ and hence in the competition for the

resource, the meso-predator is superior in comparison with the super- predator;

・ $f_3\left(y\right)=\frac{y}{38y+76}>f_2\left(x\right)$ on the corresponding domains, that is

$y\in \left[\frac{9}{10},\frac{11}{10}\right]$ and $x\in \left[\frac{19}{10}\mathrm{,}\frac{21}{10}\right]$ (see Series Time in Figure 1 and Figure 2),

therefore the super-predator is substantially benefit from the consumption to the meso-predator in the sense that its most important food source is in- termediate species.

From the above, one can conclude that our model makes ecological sense.

Acknowledgements

The first author was partially supported by CONACyT grant number CB- 2014-243722. The authors would like to thank Prof. Gamaliel Blé González and Prof. Víctor Castellanos Vargas for their helpful discussions and corrections in the preparation of this paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Polis, G.A. and Myers, C.A. (1989) The Ecology and Evolution of Intraguild Predation: Potential Competitors That Eat Each Other. Annual Review of Ecology, Evolution, and Systematics, 20, 297-330. https://doi.org/10.1146/annurev.es.20.110189.001501
[2]	Root, R. (1967) The Niche Exploitation Pattern of the Blue-Gray Gnat Catcher. Ecological Monographs, 37, 317-350. https://doi.org/10.2307/1942327
[3]	Holt, R.D. and Polis, G.A. (1996) A Theoretical Framework for Intraguild Predation. The American Naturalist, 149, 745-764. https://doi.org/10.1086/286018
[4]	Blé, G., Castellanos, V. and Llibre, J. (2016) Existence of Limit Cycles in a Tritrophic Food Chain Model with Holling Functional Responses of Type II and III. Mathematical Methods in the Applied Sciences, 39, 3996-4006. https://doi.org/10.1002/mma.3842
[5]	Castellanos, V., Llibre, J. and Quilantan, I. (2013) Simultaneous Periodic Orbits Bifurcating from Two Zero-Hopf Equilibria in a Tritrophic Food Chain Model. Journal of Applied Mathematics and Physics, 1, 31-38. https://doi.org/10.4236/jamp.2013.17005
[6]	Castellanos, V. and Chan-López, R.E. (2017) Existence of Limit Cycles in a Three Level Trophic Chain with Lotka-Volterra and Holling Type II Functional Responses. Chaos, Solitons & Fractals, 95, 157-167. https://doi.org/10.1016/j.chaos.2016.12.011
[7]	Francoise, J.P. and Llibre, J. (2011) Analytical Study of a Triple Hopf Bifurcation in a Tritrophic Food Chain Model. Applied Mathematics and Computation, 217, 7146-7154. https://doi.org/10.1016/j.amc.2011.01.109
[8]	Kuznetsov, Y.A. (2004) Elements of Applied Bifurcation Theory. 3rd Edition, Springer-Verlag. https://doi.org/10.1007/978-1-4757-3978-7
[9]	Perko, L. (2001) Differential Equations and Dynamical Systems. 3rd Edition, Springer-Verlag. https://doi.org/10.1007/978-1-4613-0003-8
[10]	Andronov, A.A., Leontovich, E.A., Gordon, I.I. and Maier, A.G. (1971) Theory of Bifurcations of Dynamic Systems on a Plane. Israel Program for Scientific Translations, Jerusalem.
[11]	Marsden, J.E. and McCracken, M. (1976) The Hopf Bifurcation and Its Applications. Springer-Verlag, New York. https://doi.org/10.1007/978-1-4612-6374-6