An optimal adaptive time-stepping scheme for solving reaction-diffusion-chemotaxis systems

1.   Introduction

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 $\mathbb{R}^3$. 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:

where $x$, $y$, $z$ represent the density of the prey, mesopredator, and superpredator, respectively; ${D}_{t_{0}}^{\alpha}$ is the Caputo fractional derivative, $\alpha$ $\in (0, 1]$, and all parameters are nonnegative. The parameters $c_1$ and $c_3$ represent the benefits of food consumption, $k_1$ and $k_2$ are the predator rates of the mesopredator, and superpredator, respectively, $c_2$ and $d$ are the mortality rate of the corresponding predators, and $\rho$ 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 $\Omega = \{(x, y, z)\in \mathbb{R}^3 | x > 0, y > 0, z > 0 \}.$ Also, $f_{2}$, $g_{2}$ and $f_{3}$, $g_{3}$ are the types Ⅱ and Ⅲ functional responses, respectively, defined as follows:

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.

2.   Preliminaries

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 $s_{i}$ of $p$ lie in the open left half plane $\mathbb{C}_{-}$.

We assume that $p(s)$ is the following polynomial:

where $A_{i}\in\mathbb{R}$, $i = 0, 1, 2, 3$ and $A_{3} > 0$, and $A_{0}\neq 0$. Using the coefficients of $p(s)$, we construct the Table 1 (see ^[15]).

Table 1.  The Routh table.

$A_{3}$	$A_{1}$	0
$A_{2}$	$A_{0}$	0
$A_{1}-\dfrac{A_{0}A_{3}}{A_{2}}$	0	0
$A_{0}$	0	0

---

 | Show Table

DownLoad: CSV

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., $A_{3} > 0$, $A_{2} > 0$, $A_{1}-\dfrac{A_{0}A_{3}}{A_{2}} > 0$, and $A_{0} > 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)\in L_{1}(\mathbb{R}^{+})$. The Riemann-Liouville fractional integral of order $\alpha\in (0, 1]$ is defined by

where $\Gamma$ is the Gamma function.

Definition 2.3. (See ^[11]) The Caputo fractional derivative of order $\alpha\in (0, 1]$ and $x\in AC(\mathbb{R}^{+})$, the space of absolutely continuous functions on $\mathbb{R}^{+}$, is defined by

where $D = d/dt$.

The autonomous nonlinear differential system in the sense of Caputo is given as follows:

where $X(t)\in \mathbb{R}^3$, $X(t_{0}) = X_{t_{0}}$, $F(X)$ is continuous.

Definition 2.4. Suppose that $P$ is an equilibrium point of system (2.2) and that all the eigenvalues $\lambda$ of the linearized matrix $J(P) = \partial F / \partial X |_{X = P}$ evaluated at $P$ satisfy: $|\lambda| \neq 0$ and $|arg(\lambda)|\neq \dfrac{\pi \alpha}{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 $\delta (P)$.

Proof. The proof can be found in ^[16].

Theorem 2.4. We consider the fractional-order system

where $0 < \alpha < 1$. System (2.3) is asymptotically stable at the equilibrium point $P$ if, and only if,

for all roots $\lambda$ of the following equation $\det(\lambda I - J(P)) = 0$.

Proof. The proof can be found in ^[17].

3.   Main result

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 $\alpha = 1$, $i = 2, \ c_{3} > d, \ a_{2}dk_{1} > a_{1}\rho(c_{3}-d), \ (c_{1}-c_{2})a_{2}dk_{1} > a_{1}c_{1}\rho (c_{3}-d)$. Then, there is only one equilibrium point $P$ in the first octant for the system (1.1),

Also, the system (1.1) is unstable around the equilibrium point $P$.

Proof. Let

The equilibrium point $P \in \Omega$ for system (1.1) when $i = 2$ is obtained by solving system $F_{1}(x, y) = F_{2}(x, y, z) = F_{3}(y, z) = 0$.

On the other hand, we obtain the Jacobian matrix $J$ of (3.1)

Now, evaluate $P$ in $J$

We denote by $J(P) = (m_{jk})$, $j, k = 1, 2, 3$. The characteristic polynomial is

where $I$ is the identity matrix 3x3. That is,

with

Let $P = (x_0, y_0, z_0)$ be the equilibrium point in the first octant. We reduce $m_{11}$ and $m_{32}$:

Then

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, \ k_{1} = 2\rho \sqrt{r}$, $c_{3} = rd$, $c_{1} = rc_{2}$, $a_{1} = a_{2}$, and $r > 4$, the system (1.1) has points of equilibrium in the first octant:

Then, the system (1.1) is unstable around the equilibrium point $P_{1}$ and asymptotically stable around $P_{2}$.

Proof. Let

Then, the Jacobian matrix $J$ of system (3.4) is

We evaluate $J$ in $P_{1}$,

Let $J(P_{1}) = (m_{jk})$, $j, k = 1, 2, 3$. We calculate the characteristic polynomial:

That is,

where

We note that $A_{0} < 0$, when $r > 4$. Then, by Theorem 2.1, $p(\lambda)$ is not a Hurwitz polynomial.

On the other hand,

The characteristic polynomial has the form (3.5), where $A_{i}$ are

Also, $r > 4$, from which we obtain

For (3.6), we have $A_{i} > 0$, for each $i = 0, 1, 2$. Now, we calculate $A_{1}A_{2}-A_{0}$.

By (3.6) and $A_{1} > 0$, we have $A_{1}A_{2}-A_{0} > 0$. Then, by Theorem 2.1, $p(\lambda)$ is a Hurwitz polynomial. Therefore, the system (1.1), when $\alpha = 1$, is stable around $P_2$.

Now, let us consider $\alpha \in (0, 1)$. By Theorem 2.3, the system (1.1) and its linearization vector field $J(P_{2})X$ are topologically equivalent. By Theorem 2.4, $D^{\alpha}X = J(P_{2})X$ is stable around $P_2$. Therefore, the system (1.1) is stable around $P_2$.

4.   Analytical solution

In this section, we present the construction of the analytical solution for the system (1.1), with Holling type Ⅲ functional response, $k_{1} = 2\rho \sqrt{r}$, $c_{3} = rd$, $c_{1} = rc_{2}$, and $a_{1} = a_{2}$, 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 $t_{0} = 0 < t_{1} < t_{2} < \dots < t_{m} = T$ and the following family of fractional-order systems

with $t \in [t_{k-1}, t_{k}]$ and initial conditions $(x_{0}, y_{0}, z_{0}) = (x_{1}(t_0), y_{1}(t_0), y_{1}(t_0))$ and $(x_k(t_{k-1}), y_k(t_{k-1}), z_k(t_{k-1})) = (x_{k-1}(t_{k-1}), y_{k-1}(t_{k-1}), z_{k-1}(t_{k-1}))$ for $k = 1, \dots, m$.

Now, we define the homotopy for each $k$:

where $p\in[0, 1]$. Let's suppose that the solution for (4.2) is given by

where $u_{ik}$, $v_{ik}$, $w_{ik}$, $i = 1, 2, \dots$, are functions to be determined. Then, the solution for $t\in[t_{k-1}, t_{k}]$ is

Therefore, the analytic solution for $t\in [0, T]$ is given by

where $I_{[t_{k-1}, t_k]}$ is the characteristic function, $t_{k} = k\dfrac{T}{m}$, and $\bigcup\limits_{k = 1}^{m} [t_{k-1}, t_k] = [0, T]$. We show the first addends of the series (4.4),

5.   Examples

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 $\alpha$ 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 $a_{1} = c_{2} = \rho = r = 1$, $c_{1} = d = k_{1} = k_{2} = 2$, $c_{3} = 2.5$, $a_{2} = 5$, $m = 50$, $t \in [0, 0.5]$, $(16, 8, 5)$ is the initial condition, and $P = (39, 20, 11.875)$ is an unstable equilibrium point.

Figure 1.  Unstable equilibrium point $P$.

DownLoad: Full-Size Img PowerPoint

Example 5.2. (See Figure 2) For $i = 3$:

Figure 2.  Unstable equilibrium point $P_{1}$ and asymptotic stability around $P_{2}$.

DownLoad: Full-Size Img PowerPoint

● Let $a_{1} = c_{2} = d = k_{2} = \rho = r = 1$, $r = 5$, $m = 100$, $t \in [0, 1]$, $(2, 1, 7)$ is the initial condition, and $P_{1} = (1.61803, 0.5, 6.54508)$ is an unstable equilibrium point.

● We take $a_{1} = c_{2} = d = k_{2} = \rho = r = 5$, $m = 100$, $t \in [0, 1]$, $(1, 1, 2)$ is the initial condition, and $P_{2} = (1.38197, 1.11803, 2.13525)$ is a locally asymptotically stable equilibrium point.

6.   Conclusions

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 $P_2$, when $\alpha \in (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 $\alpha < 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.

Author contributions

Anel Esquivel-Navarrete: investigation, formal analysis, writing – review $\And$ editing; Jorge Sanchez-Ortiz: investigation, supervision, writing – review $\And$ editing; Gabriel Catalan-Angeles: formal analysis, writing- original draft preparation, visualization; Martin P. Arciga-Alejandre: investigation, methodology, writing – review $\And$ editing.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.