A3 | A1 | 0 |
A2 | A0 | 0 |
A1−A0A3A2 | 0 | 0 |
A0 | 0 | 0 |
Citation: Chichia Chiu, Jui-Ling Yu. An optimal adaptive time-stepping scheme for solving reaction-diffusion-chemotaxis systems[J]. Mathematical Biosciences and Engineering, 2007, 4(2): 187-203. doi: 10.3934/mbe.2007.4.187
[1] | Ihtisham Ul Haq, Shabir Ahmad, Sayed Saifullah, Kamsing Nonlaopon, Ali Akgül . Analysis of fractal fractional Lorenz type and financial chaotic systems with exponential decay kernels. AIMS Mathematics, 2022, 7(10): 18809-18823. doi: 10.3934/math.20221035 |
[2] | M. Adel, M. M. Khader, Hijaz Ahmad, T. A. Assiri . Approximate analytical solutions for the blood ethanol concentration system and predator-prey equations by using variational iteration method. AIMS Mathematics, 2023, 8(8): 19083-19096. doi: 10.3934/math.2023974 |
[3] | Mohammed A. Almalahi, Satish K. Panchal, Fahd Jarad, Mohammed S. Abdo, Kamal Shah, Thabet Abdeljawad . Qualitative analysis of a fuzzy Volterra-Fredholm integrodifferential equation with an Atangana-Baleanu fractional derivative. AIMS Mathematics, 2022, 7(9): 15994-16016. doi: 10.3934/math.2022876 |
[4] | Nauman Ahmed, Ali Raza, Ali Akgül, Zafar Iqbal, Muhammad Rafiq, Muhammad Ozair Ahmad, Fahd Jarad . New applications related to hepatitis C model. AIMS Mathematics, 2022, 7(6): 11362-11381. doi: 10.3934/math.2022634 |
[5] | Amjad E. Hamza, Arshad Ali, Khaled Aldwoah, Hicham Saber, Ria Egami, Amel Touati, Amal F. Alharbi . Mathematical analysis of tri-trophic food webs with carrying capacity and Holling-type predation using fractal-fractional Caputo derivatives. AIMS Mathematics, 2025, 10(6): 13130-13150. doi: 10.3934/math.2025589 |
[6] | Sumbal Ahsan, Rashid Nawaz, Muhammad Akbar, Saleem Abdullah, Kottakkaran Sooppy Nisar, Velusamy Vijayakumar . Numerical solution of system of fuzzy fractional order Volterra integro-differential equation using optimal homotopy asymptotic method. AIMS Mathematics, 2022, 7(7): 13169-13191. doi: 10.3934/math.2022726 |
[7] | Guillaume Cantin, Cristiana J. Silva . Complex network near-synchronization for non-identical predator-prey systems. AIMS Mathematics, 2022, 7(11): 19975-19997. doi: 10.3934/math.20221093 |
[8] | Jagdev Singh, Behzad Ghanbari, Ved Prakash Dubey, Devendra Kumar, Kottakkaran Sooppy Nisar . Fractional dynamics and computational analysis of food chain model with disease in intermediate predator. AIMS Mathematics, 2024, 9(7): 17089-17121. doi: 10.3934/math.2024830 |
[9] | Yasir Nawaz, Muhammad Shoaib Arif, Kamaleldin Abodayeh, Mairaj Bibi . Finite difference schemes for time-dependent convection q-diffusion problem. AIMS Mathematics, 2022, 7(9): 16407-16421. doi: 10.3934/math.2022897 |
[10] | Rasool Shah, Abd-Allah Hyder, Naveed Iqbal, Thongchai Botmart . Fractional view evaluation system of Schrödinger-KdV equation by a comparative analysis. AIMS Mathematics, 2022, 7(11): 19846-19864. doi: 10.3934/math.20221087 |
In a tritrophic food chain, successive predation occurs between three species in the order prey-mesopredator-superpredator, in which there is a transfer of energy and nutrients when one organism eats another. The importance of the coexistence of species in the food chain, lies in avoiding overpopulation and extinction of species. There are several interesting results on the study of these systems involving Holling-type functional responses, and among them we mention: Rao and Narayan [1] studied the stability of the interior equilibrium point using the Routh-Hurwitz criterion and the global stability of a three-species food chain model with harvesting by constructing an appropriate Lyapunov function. Mammat et al. [2] studied an ecological model with a trophic chain with a classical Lotka-Volterra functional response and found a parameter space where a Hopf bifurcation occurs. Moreover, it is possible to find bifurcation points analytically and to show that the system has periodic solutions around these points. Rihand et al. [3] studied the dynamics of a two-prey one-predator system, where the growth of both prey populations is subject to Allee effects, and there is a direct competition between them. Castillo et al. [4] demonstrated the existence of a limit cycle, via the first Lyapunov coefficient and the Andronov-Hopf bifurcation theorem, for an asymmetric intragremial food web model with Holling type Ⅱ functional response for intermediate and top predators and logistic growth for (common) prey. Cheng et al. [5] reviewed the existence and local stabilities of all equilibria of the classical Holling type Ⅱ model of the three-species food chain. They obtained the existence of a single limit cycle when the limit equilibrium loses its stability, and they also demonstrated the global stability of the limit cycle in R3. Alsakaji et al. [6] studied the dynamics of a delay differential model of predator-prey system involving teams of two-prey and one-predator, with Monod-Haldane and Holling type Ⅱ functional responses, and a cooperation between the two-teams of preys against predation. Blé et al. [7] demonstrated that a tritrophic chain model with Holling type Ⅲ functional response has an equilibrium point where it presents a supercritical Hopf bifurcation independently of the prey growth rate. In the logistic case, they demonstrate the existence of at least three equilibrium points in the positive octant and one of them presents a supercritical Hopf bifurcation.
On the other hand, the theory of fractional order differential equations has gained great popularity in several scientific areas; such as biomathematics, control theory, and financial mathematics, among others (see [8,9,10,11]). In particular, in the case of biomathematics, tritrophic models with fractional derivative have been studied; for example, Maria et al. [12] calculated the equilibrium points and analyzed their stability to exhibit the dynamic behavior of a fractional order prey-predator model (3-Species). Mondal et al. [13] studied the dynamics of a three-dimensional discrete fractional-order eco-epidemiological model with Holling type Ⅱ functional response. They determined analytical conditions for the local stability of different fixed points using the Jury criterion and showed that the stability of the fractional-order discrete system depends strongly on the step size and the fractional order. More specifically, the critical value of the step size, at which stability switching occurs, decreases as the order of the fractional derivative decreases.
In the literature, we did not find results concerning the study of fractional order tritrophic systems with Holling type Ⅲ functional response.
In this work, we study the stability for the positive solution of the fractional tritrophic food chain models:
Dαt0x=ρx−k1fi(x)y,Dαt0y=c1fi(x)y−k2gi(y)z−c2y,Dαt0z=c3gi(y)z−dz, | (1.1) |
where x, y, z represent the density of the prey, mesopredator, and superpredator, respectively; Dαt0 is the Caputo fractional derivative, α ∈(0,1], and all parameters are nonnegative. The parameters c1 and c3 represent the benefits of food consumption, k1 and k2 are the predator rates of the mesopredator, and superpredator, respectively, c2 and d are the mortality rate of the corresponding predators, and ρ is the growth rate of the prey in the absence of predators. Since system (1.1) is an ecological model, our region of interest is the positive octant Ω={(x,y,z)∈R3|x>0,y>0,z>0}. Also, f2, g2 and f3, g3 are the types Ⅱ and Ⅲ functional responses, respectively, defined as follows:
fi(x)=xi−1xi−1+a1, gi(y)=yi−1yi−1+a2, i=2,3. |
The type Ⅱ functional response is used to describe predator organisms which take some time capture and ingest their prey. On the other hand, type Ⅲ functional response is characteristic of predator organisms which do not capture prey intensively below a certain level of threshold density; however, above that density level, predator organisms increase their feeding rates until some saturation level is reached.
We begin this section by defining the Hurwitz polynomial as follows.
Definition 2.1. A polynomial p is said to be a Hurwitz polynomial or Hurwitz stable if all the roots si of p lie in the open left half plane C−.
We assume that p(s) is the following polynomial:
p(s)=A0+A1s+A2s2+A3s3, | (2.1) |
where Ai∈R, i=0,1,2,3 and A3>0, and A0≠0. Using the coefficients of p(s), we construct the Table 1 (see [15]).
A3 | A1 | 0 |
A2 | A0 | 0 |
A1−A0A3A2 | 0 | 0 |
A0 | 0 | 0 |
Theorem 2.1. Consider p(s) given in (2.1). Then, p(s) is Hurwitz if, and only if, each element of the first column of the Routh table is positive, i.e., A3>0, A2>0, A1−A0A3A2>0, and A0>0.
Proof. The proof can be found in [15].
Theorem 2.2. Consider p(s) given in (2.1). Suppose when calculating the Routh table that no element in the first column is zero. Then, the number of sign changes in the first column of the Routh table is the number of open right half-plane zeros of p(s).
Proof. The proof can be found in [15].
Now, we give some definitions and properties of fractional calculus theory.
Definition 2.2. (See [11]) Let x(t)∈L1(R+). The Riemann-Liouville fractional integral of order α∈(0,1] is defined by
(Iαt0x)(t)=1Γ(α)∫tt0x(τ)(t−τ)1−α>dτ, t>t0≥0, |
where Γ is the Gamma function.
Definition 2.3. (See [11]) The Caputo fractional derivative of order α∈(0,1] and x∈AC(R+), the space of absolutely continuous functions on R+, is defined by
(Dαt0x)(t)=(I1−αt0Dx)(t), |
where D=d/dt.
The autonomous nonlinear differential system in the sense of Caputo is given as follows:
Dαt0X(t)=F(X(t)), | (2.2) |
where X(t)∈R3, X(t0)=Xt0, F(X) is continuous.
Definition 2.4. Suppose that P is an equilibrium point of system (2.2) and that all the eigenvalues λ of the linearized matrix J(P)=∂F/∂X|X=P evaluated at P satisfy: |λ|≠0 and |arg(λ)|≠πα2, then we call P a hyperbolic equilibrium point.
Theorem 2.3. If P is a hyperbolic equilibrium point of (2.2), then vector field F(x) is topologically equivalent with its linearization vector field J(P)x in the neighborhood δ(P).
Proof. The proof can be found in [16].
Theorem 2.4. We consider the fractional-order system
DαX(t)=F(X(t)), | (2.3) |
where 0<α<1. System (2.3) is asymptotically stable at the equilibrium point P if, and only if,
|arg(λ)|>πα2 |
for all roots λ of the following equation det(λI−J(P))=0.
Proof. The proof can be found in [17].
In this section, we present two results obtained by analyzing the stability around the coexistence equilibrium points of three species of the tritrophic system (1.1). First, we consider the system (1.1) with Holling type Ⅱ functional response, which has only one unstable equilibrium point in the first octant.
Theorem 3.1. We assume α=1, i=2, c3>d, a2dk1>a1ρ(c3−d), (c1−c2)a2dk1>a1c1ρ(c3−d). Then, there is only one equilibrium point P in the first octant for the system (1.1),
P=(k1a2d−ρa1(c3−d)ρ(c3−d),a2dc3−d,c3((c1−c2)k1a2d−ρa1c1(c3−d))dk1k2(c3−d)). |
Also, the system (1.1) is unstable around the equilibrium point P.
Proof. Let
F1(x,y)=ρx−k1xyx+a1,F2(x,y,z)=c1xyx+a1−k2yzy+a2−c2y,F3(y,z)=c3yzy+a2−dz. | (3.1) |
The equilibrium point P∈Ω for system (1.1) when i=2 is obtained by solving system F1(x,y)=F2(x,y,z)=F3(y,z)=0.
On the other hand, we obtain the Jacobian matrix J of (3.1)
J=(ρ−a1k1y(x+a1)2−k1xx+a10a1c1y(x+a1)2c1xx+a1−a2k2z(y+a2)2−c2−k2yy+a20a2c3z(y+a2)2c3yy+a2−d). |
Now, evaluate P in J
J(P) = (ρ(a2dk1−a1ρ(c3−d))a2dk1ρa1(c3−d)−k1a2da2d0a1c1ρ2(c3−d)a2dk21a2dk1(c1−c2)−a1c1ρ(c3−d)a2c3k1−dk2c30(c3−d)((c1−c2)a2dk1−a1c1ρ(c3−d))a2dk1k20). |
We denote by J(P)=(mjk), j,k=1,2,3. The characteristic polynomial is
p(λ)=det[J(P)−λI], | (3.2) |
where I is the identity matrix 3x3. That is,
p(λ)=A3λ3+A2λ2+A1λ+A0, | (3.3) |
with
A3=1, A2=−m11−m22, A1=m11m22−m12m21−m23m32, A0=m11m23m32. |
Let P=(x0,y0,z0) be the equilibrium point in the first octant. We reduce m11 and m32:
m11=ρ2x0k1y0, m32=a2d2z0c3y20. |
Then
A0=−a2d3k2ρ2x0z0c23k1y30<0. |
Therefore, by Theorem 2.2, the system is unstable.
Now, we consider the system (1.1) with Holling type Ⅲ functional response, which has two equilibrium points in the first octant, one stable and the other unstable.
Theorem 3.2. Assuming i=3, k1=2ρ√r, c3=rd, c1=rc2, a1=a2, and r>4, the system (1.1) has points of equilibrium in the first octant:
P1=(√a1(√r+1)√r−1,√a1r−1,√a1c2r(r+√r−2)2k2√r−1),P2=(√a1(√r−1)√r−1,√a1r−1,√a1c2r(r−√r−2)2k2√r−1). |
Then, the system (1.1) is unstable around the equilibrium point P1 and asymptotically stable around P2.
Proof. Let
G1(x,y)=ρx−2ρ√rx2yx2+a1,G2(x,y,z)=c2rx2yx2+a1−k2y2zy2+a1−c2y,G3(y,z)=dry2zy2+a1−dz. | (3.4) |
Then, the Jacobian matrix J of system (3.4) is
J = (ρ−4a1ρ√rxy(x2+a1)2−2ρ√rx2x2+a102a1c2xy(x2+a1)2c2rx2x2+a1−2a1k2yz(y2+a1)2−c2−k2y2y2+a102a1dryz(y2+a1)2dry2y2+a1−d). |
We evaluate J in P1,
J(P1) = (ρ√r−ρ(√r+1)0c2(√r−1)2c2(−r2−√r2+1√r−2r+2)−k2r0c2dk2(r2+r3/2−3r−√r+2)0). |
Let J(P1)=(mjk), j,k=1,2,3. We calculate the characteristic polynomial:
p(λ)=det[J(P1)−λI]. |
That is,
p(λ)=A3λ3+A2λ2+A1λ+A0, | (3.5) |
where
A0=−c2dρr3/2(r−1)(r+√r−2),A1=c2ρ√r(r3/22−r2−√r+1√r−2r+2)+c2dr(r2+r3/2−3r−√r+2),A2=c2(r2+√r2−2−1√r+2r)−ρ√r,A3=1. |
We note that A0<0, when r>4. Then, by Theorem 2.1, p(λ) is not a Hurwitz polynomial.
On the other hand,
J(P2) = (−ρ√r−ρ(√r−1)0c2(√r+1)2c2(−r2+√r2+2−1√r−2r)−k2r0c2dk2(r2−r3/2−3r+√r+2)0). |
The characteristic polynomial has the form (3.5), where Ai are
A0=c2dρr3/2(r−1)(r−√r−2),A1=c2ρ√r(r3/22+r2−√r−2+1√r+2r)+c2dr(r−1)(r−√r−2),A2=c2(r2−√r2−2+1√r+2r)+ρ√r,A3=1. |
Also, r>4, from which we obtain
r−√r−2>0,r3/22+r2−√r−2+1√r+2r>r2−√r2−2+1√r+2r>(√r2−2√r)2+r4−√r2>0. | (3.6) |
For (3.6), we have Ai>0, for each i=0,1,2. Now, we calculate A1A2−A0.
A1A2−A0=c2ρA1√r(r3/22+r2−√r−2+1√r+2r)+c22dr(r−1)(r−√r−2)(r2−√r2−2+1√r+2r). |
By (3.6) and A1>0, we have A1A2−A0>0. Then, by Theorem 2.1, p(λ) is a Hurwitz polynomial. Therefore, the system (1.1), when α=1, is stable around P2.
Now, let us consider α∈(0,1). By Theorem 2.3, the system (1.1) and its linearization vector field J(P2)X are topologically equivalent. By Theorem 2.4, DαX=J(P2)X is stable around P2. Therefore, the system (1.1) is stable around P2.
In this section, we present the construction of the analytical solution for the system (1.1), with Holling type Ⅲ functional response, k1=2ρ√r, c3=rd, c1=rc2, and a1=a2, which are obtained by applying the multistage homotopy perturbation method (see [18]):
First, let's define a regular partition of the interval [0,T] by t0=0<t1<t2<⋯<tm=T and the following family of fractional-order systems
Dαtk−1xk=ρxk−2ρ√rx2kykx2k+a1,Dαtk−1yk=c2rx2kykx2k+a1−k2y2kzky2k+a1−c2yk,Dαtk−1zk=dry2kzky2k+a1−dzk, | (4.1) |
with t∈[tk−1,tk] and initial conditions (x0,y0,z0)=(x1(t0),y1(t0),y1(t0)) and (xk(tk−1),yk(tk−1),zk(tk−1))=(xk−1(tk−1),yk−1(tk−1),zk−1(tk−1)) for k=1,…,m.
Now, we define the homotopy for each k:
(1−p)(Dαtk−1uk−Dαtk−1xk(tk−1))+p(Dαtk−1uk−ρuk+2ρ√ru2kvku2k+a1)=0,(1−p)(Dαtk−1vk−Dαtk−1yk(tk−1))+p(Dαtk−1vk−c2ru2kvkx2k+a1+k2v2kwkv2k+a1+c2vk)=0,(1−p)(Dαtk−1wk−Dαtk−1zk(tk−1))+p(Dαtk−1wk+cwk−dv2kwkv2k+e)=0, | (4.2) |
where p∈[0,1]. Let's suppose that the solution for (4.2) is given by
uk=u0k+pu1k+p2u2k+p3u3k+⋯,vk=v0k+pv1k+p2v2k+p3v3k+⋯,wk=w0k+pw1k+p2w2k+p3w3k+⋯, | (4.3) |
where uik, vik, wik, i=1,2,…, are functions to be determined. Then, the solution for t∈[tk−1,tk] is
xk=limp→1uk=u0k+u1k+u2k+u3k+⋯,yk=limp→1vk=v0k+v1k+v2k+v3k+⋯,zk=limp→1wk=w0k+w1k+w2k+w3k+⋯. | (4.4) |
Therefore, the analytic solution for t∈[0,T] is given by
x(t)=m∑k=1I[tk−1,tk]∞∑j=1ujk(t),y(t)=m∑k=1I[tk−1,tk]∞∑j=1vjk(t),z(t)=m∑k=1I[tk−1,tk]∞∑j=1wjk(t), | (4.5) |
where I[tk−1,tk] is the characteristic function, tk=kTm, and m⋃k=1[tk−1,tk]=[0,T]. We show the first addends of the series (4.4),
u1k(t)=ρu0k(a1+u0k(u0k−2√rv0k))(a1+u20k)Γ(1+α)(t−tk−1)α,v1k(t)=−v0k(c2(a1−(r−1)u20k)(a1+v20k)+k2(a1+u20k)v0kw0k)α(a1+u20k)(a1+v20k)Γ(α)(t−tk−1)α,w1k(t)=−dw0k(a1−(r−1)v20k)(a1+v20k)Γ(1+α)(t−tk−1)α,u2k(t)=ρu0k(a1+u20k)3Γ2(1+α)(a1ρ(a1+u20k)(a1+u0k(u0k−2√rv0k))+3ρu20k(a1+u20k)(a1+u0k×(u0k−2√rv0k))−4ρ√ru0kv0k(a1+u20k)(a1+u0k(u0k−2√rv0k))−2ρu20k(a1+u0k×(u0k−2√rv0k))2+2√ru0kv0ka1+v20k(a1+u20k)(c2(a1−(r−1)u20k)(a1+v20k)+k2(a1+u20k)×v0kw0k))(t−tk−1)2α,v2k(t)=v0k(a1+u20k)3(a1+v20k)3Γ2(1+α)((a1+u20k)((a1+v20k)c2(a1−(r−1)u20k)+k2(a1+u20k)×v0kw0k)(c2(a1−(r−1)u20k)(a1+v20k)2+2a1k2(a1+u20k)v0kw0k)+(a1+v20k)(2a1c2ρ×u20kr(a1+v20k)2(a1+u0k(u0k−2√rv0k))+dk2(a1+u20k)3v0k(a1−(r−1)v20k)w0k))×(t−tk−1)2α,w2k(t)=dw0k(a1+u20k)(a1+v20k)3αΓ2(1+α)(d(a1+u20k)(a1+v20k)(a1−(r−1)v20k)2α−2a1rv20k(c2×(a1−(r−1)u20k)(a1+v20k)+k2(a1+u20k)vk0w0k)α)(t−tk−1)2α. |
In this section, we illustrate the results obtained in Theorems 3.1 and 3.2, through some particular examples, where its graphs show the dynamics of the solutions around the equilibrium points. Furthermore, we can observe that as the order of the derivative α moves away from 1, the length of the trajectories increases. For this purpose, we use the analytical solution given in (4.5).
Example 5.1. (See Figure 1) For i=2, fix a1=c2=ρ=r=1, c1=d=k1=k2=2, c3=2.5, a2=5, m=50, t∈[0,0.5], (16,8,5) is the initial condition, and P=(39,20,11.875) is an unstable equilibrium point.
Example 5.2. (See Figure 2) For i=3:
● Let a1=c2=d=k2=ρ=r=1, r=5, m=100, t∈[0,1], (2,1,7) is the initial condition, and P1=(1.61803,0.5,6.54508) is an unstable equilibrium point.
● We take a1=c2=d=k2=ρ=r=5, m=100, t∈[0,1], (1,1,2) is the initial condition, and P2=(1.38197,1.11803,2.13525) is a locally asymptotically stable equilibrium point.
In this work, the local stability around the coexistence equilibrium points of two tritrophic fractional models were studied, with Holling type Ⅱ and Ⅲ functional responses, respectively. Under certain conditions of the parameters, it was found that only the second model has a stable equilibrium point P2, when α∈(0,1]. The multistage homotopic perturbation method was used to obtain the solution of the tritrophic fractional model with Holling type Ⅲ function response, and it was observed that the trajectories tend faster toward the equilibrium point when α<1. These results reveal the possibility of tritrophic coexistence when the interaction is Holling type Ⅲ. A future work is to use the analytical solution, obtained via the homotopic perturbation method, to solve the inverse problem associated to system (1.1); that is, to estimate the parameters involved in the system, from a Bayesian analysis perspective, which is important to model experimental problems in ecology.
Anel Esquivel-Navarrete: investigation, formal analysis, writing – review & editing; Jorge Sanchez-Ortiz: investigation, supervision, writing – review & editing; Gabriel Catalan-Angeles: formal analysis, writing- original draft preparation, visualization; Martin P. Arciga-Alejandre: investigation, methodology, writing – review & editing.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
The authors declare there is no conflict of interest.
A3 | A1 | 0 |
A2 | A0 | 0 |
A1−A0A3A2 | 0 | 0 |
A0 | 0 | 0 |