Analysis of food chain mathematical model under fractal fractional Caputo derivative

1.   Introduction

Food chains are essential phenomena for the environment in various fields such as ecological science, applied mathematics, engineering, and economics. In a food chain model species, energy and resources follow one track, while food webs have complexity because they are attached to numerous food chains. Many trophic stages have been seen in a food chain. There are many groups of creatures inside the stimulating stages, such as producers, consumers and decomposers. A formation-wise lattice arrangement is used for a food web ^[1]. Using the techniques of mathematical analysis and modeling, we can model the food chain as a differential equation. In ecology, food chains are a chain of creatures or organisms serving the organisms next to them, whereas a collection of food chains joined together forms a food web ^[2,3]. A flexible food chain theory, which shows the formation and operational characteristics of low-entity-like food webs, aims to define how to build and interact with ecosystem stability ^[4,5]. The life cycle of numerous kinds of species in nature is classified into at least two classes: mature and immature, with their behavior. The extensive study of food web models is presented here ^[6,7]. The influence of cannibalism on the environmental approach has been deliberated widely for many decades. Terrestrial and aquatic food webs have cannibalistic populations ^[8,9,10]. The stage-structured individuals commonly involve in cannibalism, whether in the inhabitants or in the aquatic food chain. Diekmann investigated and examined the cannibalism mathematical model ^[11]. An equator food web in which predator cannibalism was studied in ^[12]. Consequently, cannibalism has a very big impact on the system's dynamics. Many creatures including fish, birds, mammals, and others exhibit cannibalistic tendencies.

Fractional calculus has been developed over more than 300 years, and it is still a key idea for understanding real-world problems ^[13,14]. Numerous fractional derivatives, notably Caputo's derivative, have been presented in the literature on fractional calculus. The fractal-fractional, Atangana–Baleanu, and Caputo–Fabrizio are the most commonly used derivatives ^{[15,16,17,18,19,20,21]}. The fractal-fractional derivatives are a newly developed form of derivative that results from the recent combination of the fractal and fractional derivative concepts (FFD). Normalization of the issues of fractal-fractional orders is deliberated in ^[22,23]. The cited literature shows that the concerned models with fractal-fractional derivatives are relatively better than the integer order, which shows that these derivatives are relatively acceptable for physical and real-world problems ^[24,25]. The researchers have also revealed that FFD gives outstanding results in the development of physical modeling. The related numerical analysis and applications of FFD are given in ^{[26,27,28,29,30,31]}. The authors in ^[32], studied a fractional prey-predator with respect to harvesting rate. Bonyah et al. ^[33], proposed a listeriosis disease model, which is investigated under fractal-fractional in the sense of Caputo and Atangana-Baleanu-Caputo operators. The authors in ^[34], used a novel numerical technique for the Halvorsen system to analyze it fractionally, and discuss the chaotic behavior of the proposed system. Din and Abidin investigated a vaccinated hepatitis B model with non-singular and non-local kernels in ^[35].

The food web mathematical model is investigated through different techniques with integer and non-integer orders ^[36,37]. Motivated by the above literature, we study the food chain model ^[38] via the FFC operator, which gives a better agreement than the integer order derivative ^[39].

with initial values of $X(0), Y(0), Z(0), U(0)\geq0$. Where $0 < \delta, \beta\leq1$ and $^\mathrm{FFC}\mathrm D$ represents the FFC derivative. Motivated by the above literature, we investigate the model (1.1) via the FFC operator for the dynamics to obtain the results more precisely for parametric accuracy than integer order. Here $\delta$ represent non-integer order and $\beta$ represent a non-integer fractal dimension in $0 < \delta, \; \beta < 1$.

The explanation of the parameters used in the given model (1.1) is as, $X(t)$ represent prey(lower level species) density at time $t$; $Y(t)$ represent intermediate predator density at time $t$; top predator density (mature and immature of top-level kinds) at time $t$ is symbolized by $Z(t), U(t)$ respectively. $\eta$ represents the inherent growth rate and $\beta$ represents the rate of transferring capability with the logical growth of the prey. The intermediate predator eats the prey at the low level, with an occurrence rate of $\mu_1$ and an adaptation rate of $\epsilon_1$, according to the Lotka-Volterra functional retaliate. It continues to decay exponentially with a natural death rate $\rho_{1}$ because they lack a food supply. Mature and immature are the two classes of upper predators in the proposed model. The growth rate of Immature inhabitants is considered rapidly along with their maternities denoted by $\mu_{2}$, while the mature inhabitants with growth rate $\nu$, while a part of the population grows up to become a mature population with an expansion rate $\tau$. Furthermore, $\mu_{2}$ and $\nu$ face usual death with mortality rates of $\rho_{2}$ and $\rho_{3}$. With the concentrated dose rate $\mu_{3}$ and adaptation rate $\epsilon_{3}$, the mature top predators spells the middle predators by the Lotka–Volterra functional retaliation. The lake of the accessibility of their favorite foodstuff, they disassemble the upper immature predators centered on the Lotka–Volterra functional retaliation with risky violence rate $\mu_{3}$ and adaptation rate $\epsilon_{3}$.

This paper is organized as follows: In Section 2, we present the basic definition of fractional operators from the literature. In Section 3, by using the approach of fixed point theory, we find the existence and uniqueness of the solution along with Ulam-Hyers stability for the considered model. The approximate solution is obtained for the aforementioned model with the help of the fractional Adams-Bashforth technique in Section 4. The numerical findings of the considered model have been plotted graphically and discussed their dynamical behavior in Section 5. Finally, we conclude our work in Section 6.

2.   Basic definitions

Definition 1: ^[40] Suppose $\mathcal{P}(t)$ on $t\in(a, b)$ is continuous as well as differentiable function, the fractal-fractional operator of order $\delta$ and dimension $\beta$ is defined as

where $p-1 < \delta$, $\beta\leq p$, for $p\in \mathrm{N}$ and $\frac{d\mathcal{P}(y)}{dy^\beta} = \lim_{t\rightarrow0}\frac{\mathcal{P}(t)-\mathcal{P}(y)}{t^\beta-y^\beta}$.

Definition 2: ^[40] Suppose $\mathcal{P}(t)$ is a continuous function and $t\in(a, b)$, the fractal-fractional integral of order $\delta$ is given as

Definition 3: The model (1.1) shows U-H stability if there is a real number $\mathcal{G}_{\delta, \beta}\geq0$ such that $\forall$ $\vartheta > 0$ and all the roots $\Upsilon\in C^1(\mathcal{X}, \mathbb{R})$, the inequality as

$\mathcal{Y}\in C^1(\mathcal{X}, \mathbb{R})$ is the only one root of model (1.1), $\ni$

Definition 4: Consider a differentiable function $f(t)\in H^{1}$ in interval $(a, b)$, where $a < b$, and $\delta\in[0, 1]$, the Caputo derivative is define as

$\Gamma(.)$ represent the gamma function, and define as

Note: For the qualitative analysis, consider a Banach space $\mathrm{U} = \mathrm{X}\times\mathrm{X}\times\mathrm{X}\times\mathrm{X}$ where $\mathrm{X} = \mathcal{G}(\mathcal{X})$ with norm: $\|\Upsilon\| = \lVert {X(t)}, {Y(t)}, {Z(t)}, {U(t)} \rVert = \max_{t\in[0, T]}\bigg\{|{X}|+|{Y}|+|{Z}|+|{U}|\bigg\}$.

3.   Qualitative analysis

Here, we examine the existence and uniqueness of the solution of the given system (1.1).

3.1.   Existence

Nonlinear and non-local behavior characterized the suggested system $(1.1)$. There are no particular methods for figuring out the nonlinear system's exact roots. In rare circumstances, nevertheless, it could have an exact result. Here, we apply the functional analysis rule to determine if the system under consideration has a solution. The right side of the system $(1.1)$ as a result of the suggested integral being differentiable as

For $t\in\Lambda$ the proposed system maybe written as

Using the Riemann-Liouville integral and replace $^\mathcal{RL}\mathrm{D}^{\delta, \beta}$ by $^\mathcal{C}\mathrm{D}^{\delta, \beta}$ the solution of $(3.2)$ will be obtain as

for relation

Now, we are in a position to convert system $(1.1)$ to a fixed-point phenomenon with operator $\mathbf{{\digamma}}:{\mathbf{S}}\rightarrow {\mathbf{S}}$ defined as

The following theorem is used to analyze the existence results for the suggested model $(1.1$).

Theorem 1. ^[41] Let the operator ${{\digamma}}:{\mathbf{S}}\rightarrow {\mathbf{S}}$, is said to be a completely-continuous mapping if

is bounded, then $\mathrm{{\digamma}}$ has at least one fixed-point in $\mathbf{S}.$

Theorem 2. Suppose $\Upsilon:\Xi\times\mathbf{S}\rightarrow \mathrm{R}$ be a continuous operator, then $\mathrm{{\digamma}}$ shows the compactness.

Proof. To prove this theorem, first we have show that ${{\digamma}}:{\mathbf{S}}\rightarrow {\mathbf{S}}$ in Eq $(3.3)$ is continuous. Consider a bounded subset $\mathrm{E}$ of $\mathbf{S}$, then there exists $\mathcal{H}_{\Upsilon} > 0$ with $|\Upsilon(t, \Psi(t))|\leq\mathcal{H}_{\Upsilon}, \forall \, \, \, \Psi\in\mathrm{E}$. And any $\Psi\in\mathrm{E}$ we have

Hence, Eq $(3.6)$, shows that the operator $\mathrm{\digamma}$ is uniformly bounded, where $\mathrm{E}(\delta, \beta)$ represent Beta function. Furthermore, we have to prove equi-continuity of the operator $\mathrm{\digamma}, \forall \, \, \, t_{1}, t_{2}\in \Xi$ and $\Psi\in \mathrm{E}$, we get

Which shows that $\digamma$ is equi-continuous, hence the operator is bounded and as well as continuous, therefore by "Arzel$\acute{a}$-Ascoli" theorem, $\mathcal{\digamma}$ is relatively-compact and so completely continuous. Further, we use the hypothesis

(a) There exists a constant $M_{\Upsilon} > 0$ such that for every $\Psi, {\bar{\Psi}}\in \mathbb{f}$ we have

3.2.   Existence of unique solution

Here, we are going to study the uniqueness of the solution for model $(1.1)$ with the aid of fixed-point theory ^[41].

Theorem 3. The suggested model $(1.1)$ has unique solution by using hypothesis $(a)$ and for $\beta < 1$ as

Proof. Suppose that, $\max_{0 < t < T}|\Upsilon(t, 0)| = \mathcal{V}_\Upsilon < \infty$, such that

here, we investigate that $\digamma(E_{r})$ is a subset of $E_{r}$ and $E_{r} = \big\{ \Psi\in\mathcal{V}:||\Psi||\leq r\big\}$ where $\Psi\in E_{r}$, so we have

According to Eq $(3.5)$, the operator $\digamma$ is defined and by hypothesis $(a)$ for all $t\in \Xi, \Psi, \bar{\mathcal{W}\in \Xi}$, we get

Hence, the operator $\digamma$ has contraction by Eq $(3.9)$. Therefore the equation Eq $(3.3)$ has unique solution and so model $(1.1)$ has a unique solution.

3.3.   Ulam-Hyers stability

Here, we have to study the stability analysis for the proposed problem $(1.1)$ for this we use the well-known theorem of functional analysis the "Hyers-Ulam" type stability analysis, holding that $\vartheta$ is independent, i.e., $\vartheta(0) = 0$ and $\vartheta\in C(\Xi)$, then

● $|\vartheta(t)|\leq \alpha, \ \text{for}\ \alpha > 0$;

● $^\mathit{FF}\mathrm{D}_t^{\delta, \beta} \Psi(t) = \vartheta(t, \Psi(t))+\vartheta(t).$

Lemma 1. The solution of a perturb equation is

satisfying

Theorem 4. From Eq $(3.11)$ and supposition $({a})$, the solution of the suggested model (1.1) is U-H stable, so the obtain result of the given system is U-H stable if $\beta < 1$, where $\beta$ is in Eq $(3.7)$.

Proof. Let $\mathcal{G}\in \mathbf{S}$ has a unique solution and $\Psi\in \mathbf{S}$ is the solutions of Eq $(3.3)$, further we use fractal-fractional integral Definition 2, we obtain

hence

Equation $(3.12)$, can be written as

Hence, from Eq $(3.13)$ satisfying all the conditions of Ulam-Hyers stability so we claimed that Eq $(3.3)$ shows that the solution of the proposed system is stable.

4.   Approximate solution

In this part of the manuscript, we find the approximate solution of the proposed system $(1.1)$ under numerical scheme of fractional Adams Bashforth iterative technique ^[42]. The given system may be written as

and

to calculate the approximate solution for Eq $(4.2)$ using the advance iterative technique $t_{k+1}$. The approximate solution for the first class of the proposed system as

we get the approximate integral as in the form

For an infinite values of $[t_\alpha, t_{\alpha+1}]$ in the form of Lagrange interpolation polynomial with function $\mathcal{N}_1({X}, y)$ along with $\hbar = [t_\alpha-t_{\alpha-1}]$ such that

put Eq $(4.8)$ in Eq $(4.7)$ we get

In right hand side of the integral of Eq $(4.9)$ gives an approximate solution for the class $X(t)$ in proposed system by using FF derivative operator with Caputo derivative operator.

similarly, the remaining terms might be like as

5.   Numerical simulation

In this portion, the desired analytical results are simulated via MATLAB-18. The validity and efficiency are verified via numerical simulation of the numerical results of the food web under the FFC operator. The evolution of all the classes of the proposed food web model is provided for a few sets of fractional order $\delta$ and fractal dimension $\beta$. For the required simulation, we choose the parameter values that are given in Table 1 from ^[38].

In the numerical simulation, we have provided the graphical representation of all four compartments of the proposed food-web model on different fractional orders $\delta$ and fractal dimensions $\beta$ in Figure 1(a)–(d). All the quantities show chaotic behavior as the food web depends on all the compartmental values. The Figure 1(a) shows the dynamics of the prey population fluctuating and then becoming stable on different fractional orders. The Figure 1(b) shows the dynamics of an intermediate predator class, which also fluctuates and then becomes convergent. Figure 1(c) is for a mature predator population, which shows chaotic behavior along with an increase in its numbers and then becomes stable. Figure 1(d) is for the immature predator population, which also shows the oscillation and increase moving towards stability.

Next, Figure 2(a)–(d) shows the chaotic behaviors in the 3D figures related to each other. In Figure 2(a), the three quantities of ${X}, {Y}$ and ${Z}$ are presented showing their relation depending on each other and converges with the passage of time. Figure 2(b) shows the dynamics of ${X}$ and ${Y}$ with circulating motion and then converges to a point with zero radius. Figure 2(c) represents the dynamics of $\mathcal{BC}$ showing chaotic behaviors and dependence of ${Y}$ and ${Z}$ on each other.

In the rest of the figures, Figure 3(a)–(d) shows the dynamics of four food-web agents with the same data but only changes the initial data of first agent prey populations ${X}$ in 2D representation.

Figure 4(a)–(d) represents the dynamics of four food-web agents with the same data but only changes the initial data of first agent prey populations ${X}$ in 3D representation.

6.   Conclusions

We have investigated a three-species food chain mathematical model under fractal-fractional derivative in the sense of Caputo, where top predators are stage-structured with a mature predator having a cannibalistic feature. At the first level, the prey grows logistically in the absence of the predators. The existence and uniqueness of the solution are studied using fixed point theory. The stability analysis of the proposed model is studied with the help of the Ulam-Hyers technique. The fractional Adams-Bashforth iterative scheme is applied for the numerical calculations. The results are studied for various fractional orders $\delta$ and fractal dimension $\beta$. This analysis gives us the results for the food chain of prey and predator in the ecosystem. It also provides us a stable situation for both the species on different fractal dimension and fractional orders in the complex geometrical analysis. In the graphical representation, we also provide the dependence of each species in the environment which shows the chaotic behavior. Each and every quantity has been shown in a spectral format which shows the total density of all compartments which will be effective for checking inside behavior lying between 0 and 1. We may also study the proposed system by global piecewise derivative for the crossover dynamical behavior along with the existence and uniqueness of the solution and numerical solutions.

Acknowledgements

The fourth author would like to thank to his mentor Professor Dr. Poom Kuman from King Mongkut's University of Technology Thonburi, Thailand for his advice and comments to improve quality the results of this paper. This work (Grant No. RGNS 65-168) was financially supported by Office of the Permanent Secretary, Ministry of Higher Education, Science, Research and Innovation.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have influenced or appeared to influence the work reported in this paper.