A new nonlinear triadic model of predator–prey based on derivative with non-local and non-singular kernel

Abstract

The model of predator–prey has been used by many researchers to predict the animal growth population in many countries in the world. However, the system of equations used in these models assumes a density of prey and predator, respectively, but when looking at the real-world situation, most of the time we have some prey that act at the same time like a predator, for instance, hyena, and also some predators that act as a prey, for instance, lions. In this research, we proposed a new model of triadic prey–prey–predator. The new model was constructed using the new fractional differentiation based on the generalized Mittag-Leffler function due to the non-locality of the dynamical system of the three species. We presented the existence of a positive set of the solutions for the new model. The uniqueness of the positive set of the solutions was presented in detail. The new model was solved numerically using the Crank–Nicolson numerical scheme.

Introduction

The description of the dynamics of biological systems, in which only two different species interrelate, where one species is considered as a predator and the other a prey, has been a concern for many reputed researchers in the past decade. In the literature, one of the most used nonlinear partial differential equations to describe this dynamical interaction is called the Lotka–Volterra equations. Many researches using different concepts of differentiations have been done and some good predictions were obtained.^1–12 However, the real observation is more complex than the suggested model. For instance, it has been observed by many field worker biologists, and also many videos played by national geography channel show a group of hyena killing and eating lions. In this situation, the lions become prey and hyena predator. In another situation, lions kill and eat hyena; in this one, the hyena are considered as prey and the lions predators. Let us now add species to this situation; for instance, a zebra that is a prey for both lions and hyena; in this situation, lions have two prey although the effort to kill both is not the same. Also, the hyena has the possibility of eating lions and zebra. The description of the dynamics of biological system of these three species cannot be achieved using the predator–prey model. Therefore, a new model is needed to portray this situation. In this work, the triadic model of prey–predator is suggested using a new concept of fractional differentiation. The choice of fractional differentiation is motivated by the fact that the interaction is not local but global and also the trend observed in the field does not follow the power law with $x^{-\alpha }$ ; therefore, classical Caputo and Riemann–Liouville differentiation could be used in this case. Also, the classical differentiation is not suitable in this case due to non-locality of the interaction. Recently, a new concept of differentiation was suggested to handle the limitation of power, and a new fractional derivative was suggested with non-local and non-singular kernel.^13,14 Many studies have been done around this new suggested derivative with great results.^15–23

Atangana–Baleanu derivative in Caputo sense

In this section, we present the definitions of the new fractional derivative with no singular and non-local kernel.^15–19

Definition 1

Let $f\in H^1(a,b)$ , $b>a$ , $\alpha \in [0,1]$ , then the definition of the new fractional derivative (Atangana–Baleanu derivative in Caputo sense) is given as

$${{}_a{}^{\mathrm{ABC}}D}_t^{\alpha }(f(t))=\frac{B(\alpha )}{1-\alpha }\underset{a}{\overset{t}{\int }}f\prime (x)E_{\alpha }[-\alpha \frac{{(t-x)}^{\alpha }}{1-\alpha }]\mathrm{dx}$$

(1)

Definition 2

Let $f\in H^1(a,b)$ , $b>a,$ $\alpha \in [0,1]$ and not necessary differentiable, then the definition of the new fractional derivative (Atangana–Baleanu fractional derivative in Riemann–Liouville sense) is given as

$${{}_a{}^{\mathrm{ABR}}D}_t^{\alpha }(f(t))=\frac{B(\alpha )}{1-\alpha }\frac{d}{\mathrm{dt}}\underset{a}{\overset{t}{\int }}f(x)E_{\alpha }[-\alpha \frac{{(t-x)}^{\alpha }}{1-\alpha }]\mathrm{dx}$$

(2)

Definition 3

The fractional integral associated with the new fractional derivative with non-local kernel is defined as

$${{}_a{}^{\mathrm{AB}}I}_t^{\alpha }\{f(t)\}=\frac{1-\alpha }{B(\alpha )}f(t)+\frac{\alpha }{B(\alpha )\Gamma (\alpha )}\underset{a}{\overset{t}{\int }}f(y){(t-y)}^{\alpha -1}\mathrm{dy}$$

(3)

The initial function is recovered when the fractional order turns to 0. Also, when the order turns to 1, we have the classical integral.

New nonlinear triadic model with Atangana–Baleanu derivative in Caputo sense

In the past year, modelers were interested in modeling the dynamic interaction of the predator and the prey. The model considered only two species, predator and prey; however, when looking at the real-world situation, we have a situation where the prey is at the same time a predator and prey. Let us consider the population of lions, hyena, and zebra; the hyena is a prey and a predator at the same time; they can kill lions with a certain probability, and they can kill zebra and can also be killed by lions. If we consider $u(x,t)$ , $v(x,t)$ , and $y(x,t)$ the density at time t and space x of lions, hyena, and zebra, respectively, then it is clear that the interaction is not local but rather global. To model this dynamic behavior, we need to use a concept of differentiation that takes into account the non-locality of the physical problem. The prey–predator model does not always follow the power law based on the function $x^{-\alpha }$ ; therefore, the model cannot be built on the fractional derivative based on power law. Nevertheless, the generalized Mittag-Leffler function can be used where the power law failed, and therefore in this work we make use of the differentiation based on the generalized Mittag-Leffler function. Taking into account the observed facts, the nonlinear triadic predator–prey model with ABC sense is given as below

$$\begin{array}{c}{{}_0{}^{\mathrm{ABC}}D}_t^{\alpha }u(x,t)+(1-{\rho }_1)v_1(x,t){\partial }_{x}u(x,t)+{\rho }_1{\varphi }_1{\partial }_{t}u(x,t)=-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi hu(x,t)}+r_{1}u(x,t)\left(1-\frac{u(x,t)}{k}\right)\\ -au(x,t)v(x,t)-bu(x,t)\\ {{}_0{}^{\mathrm{ABC}}D}_t^{\mu }v(x,t)+(1-{\rho }_2)v_2(x,t){\partial }_{x}v(x,t)+{\rho }_2{\varphi }_2{\partial }_{t}v(x,t)=-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi hv(x,t)}+r_{2}v(x,t)\left(1-\frac{v(x,t)}{k}\right)\\ -cu(x,t)v(x,t)-dv(x,t)\\ {{}_0{}^{\mathrm{ABC}}D}_t^{\theta }y(x,t)+(1-{\rho }_3)v_3(x,t){\partial }_{x}y(x,t)+{\rho }_3{\varphi }_3{\partial }_{t}y(x,t)=\frac{\beta u(x,t)v(x,t)y(x,t)}{1+\chi hu(x,t)v(x,t)}\\ -fy(x,t)-gy(x,t)v(x,t)-ey(x,t)u(x,t)\end{array}$$

(4)

where $\chi$ is the parameter of interaction of growth together of the pieces for the predator and $\beta$ is the parameter of interaction of growth together of the pieces for the prey. In the model above, $v_1$ , $v_2$ , and $v_3$ are the average velocity of lions, hyena, and zebra, respectively, and ${\rho }_1$ , ${\rho }_2$ , and ${\rho }_3$ are the factors of reduction of running of lions, hyena, and zebra, respectively, due to diseases, age, and starvation. ${\rho }_1{\varphi }_1{\partial }_{t}u(x,t)$ , ${\rho }_2{\varphi }_2{\partial }_{t}v(x,t)$ , and ${\rho }_3{\varphi }_3{\partial }_{t}y(x,t)$ terms contribute to the natural deposition of a given set of prey and predators to defend and hunt, respectively. $r_{1}u(x,t)(1-((u(x,t))/(k)))$ and $r_{2}u(x,t)(1-((u(x,t))/(k)))$ represent logistic growth where $r_1$ and $r_2$ are the intrinsic growth rates and k is the carrying capacity term. Finally, b, d, and f are the terms for natural death rate of the prey and predator and a, c, g, and e are the terms for per capita prey reduction due to consumption by the predator.

Existence and uniqueness of solutions

In this section, we will research the existence and uniqueness of the solutions. Also, we will present the proof of uniqueness of positive solutions. Let us consider $X=C[a,b]$ the Banach space of every continuous real function defined in the closed set $[a,b]$ , which contains the subnorm, and Z be the shaft defined as $Z=\{u,v\in X,u(x,t)\ge 0\mathrm{and}v(x,t)\ge 0,a\le t\le b\}$ . Now we present the following Banach fixed-point theorem that will be used for the existence of the solutions.

Definition 4

Let E be a real Banach space with a cone H. H initiates a restricted order ≤ in E in the succeeding approach ²⁴

$$x\le y\Rightarrow y-x\in H$$

For every $x,y\in E$ , the order interval is defined as $〈a,b〉=\{f\in E:a\le f\le b\}$ . A cone K is denoted normal if one can find a positive constant j such that $h,d\in K$ , $\Phi <h<d\Rightarrow \Vert h\Vert \le j\Vert d\Vert$ , where $\Phi$ denotes the zero element of K.

Theorem 1

Let H be a closed set subspace of a Banach space of D. ²⁴ Let T be a contraction mapping with Lipschitz constant $g<1$ from H to H. Thus, T possesses a fixed-point $t^{*}$ in H. In addition, if $t_0$ is a random point in H and $\{t_n\}$ is a sequence defined by $t_{n+1}=T{t}_n(n=0,1,2,\dots )$ , then for a large number n, t_n tends to $t^{*}$ in H and $d(t_n,t^{*})\le ((g^n/(1-g))d(t_1,t_0))$ .

To investigate the existence of the solutions, let us consider triadic model with ABC derivative as below

$$\begin{array}{c}{{}_0{}^{\mathrm{ABC}}D}_t^{\alpha }u(x,t)=-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi hu(x,t)}+r_{1}u(x,t)\left(1-\frac{u(x,t)}{k}\right)\\ -(1-{\rho }_1)v_1(x,t){\partial }_{x}u(x,t)-{\rho }_1{\varphi }_1{\partial }_{t}u(x,t)\\ -au(x,t)v(x,t)-bu(x,t)\\ {{}_0{}^{\mathrm{ABC}}D}_t^{\gamma }v(x,t)=-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi hv(x,t)}+r_{2}v(x,t)\left(1-\frac{v(x,t)}{k}\right)\\ -(1-{\rho }_2)v_2(x,t){\partial }_{x}v(x,t)-{\rho }_2{\varphi }_2{\partial }_{t}v(x,t)\\ -cu(x,t)v(x,t)-dv(x,t)\\ {{}_0{}^{\mathrm{ABC}}D}_t^{\theta }y(x,t)=\frac{\beta u(x,t)v(x,t)y(x,t)}{1+\chi hu(x,t)v(x,t)}-fy(x,t)-gy(x,t)v(x,t)\\ -(1-{\rho }_3)v_3(x,t){\partial }_{x}y(x,t)-{\rho }_3{\varphi }_3{\partial }_{t}y(x,t)\\ -ey(x,t)u(x,t)\end{array}$$

(5)

Now applying the AB fractional integral in equation (5), we obtain the following

$$\begin{array}{cc}\multicolumn{1}{c}{u(x,t)-u(x,0)}& =\frac{1-\alpha }{\mathrm{AB}(\alpha )}\left\{\begin{array}{c}-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi \mathrm{hu}(x,t)}+r_{1}u(x,t)(1-\frac{u(x,t)}{k})\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_1)v_1(x,t){\partial }_{x}u(x,t)}\\ \multicolumn{1}{c}{-{\rho }_1{\varphi }_1{\partial }_{t}u(x,t)}\\ \multicolumn{1}{c}{-\mathrm{au}(x,t)v(x,t)-\mathrm{bu}(x,t)}\end{array}\right\}\hfill \\ \multicolumn{1}{c}{}& +\frac{\alpha }{B(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\alpha -1}.\left\{\begin{array}{c}-\frac{\chi u(x,\tau )v(x,\tau )y(x,\tau )}{1+\chi \mathrm{hu}(x,\tau )}+r_{1}u(x,\tau )(1-\frac{u(x,\tau )}{k})\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_1)v_1(x,\tau ){\partial }_{x}u(x,\tau )-{\rho }_1{\varphi }_1{\partial }_{t}u(x,\tau )}\\ \multicolumn{1}{c}{-\mathrm{au}(x,\tau )v(x,\tau )-\mathrm{bu}(x,\tau )}\end{array}\right\}d\tau \hfill \end{array}$$

(6)

$$\begin{array}{cc}\multicolumn{1}{c}{v(x,t)-v(x,0)}& =\frac{1-\gamma }{\mathrm{AB}(\gamma )}\left\{\begin{array}{c}-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi \mathrm{hv}(x,t)}+r_{2}v(x,t)(1-\frac{v(x,t)}{k})\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_2)v_2(x,t){\partial }_{x}u(x,t)-{\rho }_2{\varphi }_2{\partial }_{t}v(x,t)}\\ \multicolumn{1}{c}{-\mathrm{cu}(x,t)v(x,t)-\mathrm{dv}(x,t)}\end{array}\right\}\hfill \\ \multicolumn{1}{c}{}& +\frac{\gamma }{B(\gamma )\Gamma (\gamma )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\gamma -1}.\left\{\begin{array}{c}-\frac{\chi u(x,\tau )v(x,\tau )y(x,\tau )}{1+\chi \mathrm{hv}(x,\tau )}+r_{2}v(x,\tau )(1-\frac{v(x,\tau )}{k})\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_2)v_2(x,\tau ){\partial }_{x}u(x,\tau )}\\ \multicolumn{1}{c}{-{\rho }_2{\varphi }_2{\partial }_{t}v(x,\tau )}\\ \multicolumn{1}{c}{-\mathrm{cu}(x,\tau )v(x,\tau )-\mathrm{dv}(x,\tau )}\end{array}\right\}d\tau \hfill \end{array}$$

(7)

$$\begin{array}{cc}\multicolumn{1}{c}{y(x,t)-y(x,0)}& =\frac{1-\theta }{\mathrm{AB}(\theta )}\left\{\begin{array}{c}\frac{\beta u(x,t)v(x,t)y(x,t)}{1+\chi \mathrm{hu}(x,t)v(x,t)}-\mathrm{fy}(x,t)-\mathrm{gy}(x,t)v(x,t)\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_3)v_3(x,t){\partial }_{x}y(x,t)}\\ \multicolumn{1}{c}{-{\rho }_3{\varphi }_3{\partial }_{t}y(x,t)}\\ \multicolumn{1}{c}{-\mathrm{ey}(x,t)u(x,t)}\end{array}\right\}\hfill \\ \multicolumn{1}{c}{}& +\frac{\theta }{B(\theta )\Gamma (\theta )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\theta -1}.\left\{\begin{array}{c}\frac{\beta u(x,\tau )v(x,\tau )y(x,\tau )}{1+\chi \mathrm{hu}(x,\tau )v(x,\tau )}-\mathrm{fy}(x,\tau )-\mathrm{gy}(x,\tau )v(x,\tau )\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_3)v_3(x,\tau ){\partial }_{x}u(x,\tau )-{\rho }_3{\varphi }_3{\partial }_{t}u(x,\tau )}\\ \multicolumn{1}{c}{-\mathrm{ey}(x,\tau )u(x,\tau )}\end{array}\right\}d\tau \hfill \end{array}$$

(8)

Now we can use equations (6)–(8) to show the existence of equation (5). Necessary lemmas for the existence of the solutions are given as Lemma 1 and Lemma 2 below.

Lemma 1

The mapping $T:H\to H$ is defined as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{Tu}(x,t)}& =\frac{1-\alpha }{\mathrm{AB}(\alpha )}s(x,t,u(x,t))+\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\hfill \\ \multicolumn{1}{c}{}& \underset{0}{\overset{t}{\int }}{(t-\tau )}^{\alpha -1}s(x,\tau ,u(x,\tau ))d\tau \hfill \\ \multicolumn{1}{c}{\mathrm{Tv}(x,t)}& =\frac{1-\gamma }{\mathrm{AB}(\gamma )}s(x,t,v(x,t))\hfill \\ \multicolumn{1}{c}{}& +\frac{\gamma }{\mathrm{AB}(\gamma )\Gamma (\gamma )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\gamma -1}s(x,\tau ,v(x,\tau ))d\tau \hfill \\ \multicolumn{1}{c}{\mathrm{Ty}(x,t)}& =\frac{1-\theta }{\mathrm{AB}(\theta )}s(x,t,y(x,t))+\frac{\theta }{\mathrm{AB}(\theta )\Gamma (\theta )}\hfill \\ \multicolumn{1}{c}{}& \underset{0}{\overset{t}{\int }}{(t-\tau )}^{\theta -1}s(x,\tau ,y(x,\tau ))d\tau \hfill \end{array}$$

(9)

To be dealt with more easily, let us consider below

$$\begin{array}{c}s(x,t,u(x,t))=\left\{\begin{array}{c}-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi \mathrm{hu}(x,t)}+r_{1}u(x,t)(1-\frac{u(x,t)}{k})\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_1)v_1(x,t){\partial }_{x}u(x,t)-{\rho }_1{\varphi }_1{\partial }_{t}u(x,t)}\\ \multicolumn{1}{c}{-\mathrm{au}(x,t)v(x,t)-\mathrm{bu}(x,t)}\end{array}\right\}\hfill \\ s(x,t,v(x,t))=\left\{\begin{array}{c}-\frac{\chi u(x,t)v(x,t)y(x,t)}{1+\chi \mathrm{hv}(x,t)}+r_{2}v(x,t)(1-\frac{v(x,t)}{k})\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_2)v_2(x,t){\partial }_{x}u(x,t)-{\rho }_2{\varphi }_2{\partial }_{t}v(x,t)}\\ \multicolumn{1}{c}{-\mathrm{cu}(x,t)v(x,t)-\mathrm{dv}(x,t)}\end{array}\right\}\hfill \\ s(x,t,y(x,t))=\left\{\begin{array}{c}\frac{\beta u(x,t)v(x,t)y(x,t)}{1+\chi \mathrm{hu}(x,t)v(x,t)}-\mathrm{fy}(x,t)-\mathrm{gy}(x,t)v(x,t)\hfill \\ \multicolumn{1}{c}{-(1-{\rho }_3)v_3(x,t){\partial }_{x}y(x,t)-{\rho }_3{\varphi }_3{\partial }_{t}y(x,t)}\\ \multicolumn{1}{c}{-\mathrm{ey}(x,t)u(x,t)}\end{array}\right\}\hfill \end{array}$$

(10)

Lemma 2

Let $M\subset H$ be bounded implying, we can find $n,r,s>0$ for system such that

$$\begin{array}{c}\Vert u(x,t_2)-u(x,t_1)\Vert \le n\Vert t_2-t_1\Vert \hfill \\ \Vert v(x,t_2)-v(x,t_1)\Vert \le r\Vert t_2-t_1\Vert \hfill \\ \Vert y(x,t_2)-y(x,t_1)\Vert \le s\Vert t_2-t_1\Vert \forall u,v,y\in M\hfill \end{array}$$

(11)

Then $\overline{T(M)}$ is compact.

Proof

Let $N=max\{((1-\alpha )/(\mathrm{AB}(\alpha )))+s(x,t,u(x,t))\}$ , $0\le u(x,t)\le P$ . For $u(x,t)\in M$ , then we have the following

$$\begin{array}{cc}\multicolumn{1}{c}{\Vert \mathrm{Tu}(x,t)\Vert }& =\frac{1-\alpha }{\mathrm{AB}(\alpha )}\Vert s(x,t,u(x,t))\Vert \hfill \\ \multicolumn{1}{c}{}& +\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\alpha -1}\Vert s(x,\tau ,u(x,\tau ))\Vert d\tau \hfill \\ \multicolumn{1}{c}{}& \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}·N+\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}·N\frac{t^{\alpha }}{\alpha }\hfill \\ \multicolumn{1}{c}{}& \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}·N+\frac{\alpha t^{\alpha }N}{\mathrm{AB}(\alpha )\Gamma (\alpha +1)}\hfill \end{array}$$

(12)

If $R=max\{((1-\gamma )/(\mathrm{AB}(\gamma )))+s(x,t,v(x,t))\}$ , $0\le v(x,t)\le Q$ . For $v(x,t)\in M$ , then

$$\begin{array}{cc}\multicolumn{1}{c}{\Vert \mathrm{Tv}(x,t)\Vert }& \le \frac{1-\gamma }{\mathrm{AB}(\gamma )}\Vert s(x,t,v(x,t))\Vert \hfill \\ \multicolumn{1}{c}{}& +\frac{\gamma }{\mathrm{AB}(\gamma )\Gamma (\gamma )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\gamma -1}\Vert s(x,\tau ,v(x,\tau ))\Vert d\tau \hfill \\ \multicolumn{1}{c}{}& \le \frac{1-\gamma }{\mathrm{AB}(\gamma )}·R+\frac{\gamma R{t}^{\gamma }}{\mathrm{AB}(\gamma )\Gamma (\gamma +1)}\hfill \end{array}$$

(13)

And if $S=max\{((1-\theta )/(\mathrm{AB}(\theta )))+s(x,t,y(x,t))\}$ , $0\le y(x,t)\le V$ . For $y(x,t)\in M$ , then we have the following

$$\begin{array}{cc}\multicolumn{1}{c}{\Vert \mathrm{Gv}(x,t)\Vert }& \le \frac{1-\theta }{\mathrm{AB}(\theta )}\Vert s(x,t,y(x,t))\Vert \hfill \\ \multicolumn{1}{c}{}& +\frac{\theta }{\mathrm{AB}(\theta )\Gamma (\theta )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\theta -1}\Vert s(x,\tau ,y(x,\tau ))\Vert d\tau \hfill \\ \multicolumn{1}{c}{}& \le \frac{1-\theta }{\mathrm{AB}(\theta )}·S+\frac{\theta S{t}^{\theta }}{\mathrm{AB}(\theta )\Gamma (\theta +1)}\hfill \end{array}$$

(14)

So equations (12)–(14) imply that function T is bounded.

Now in the following part, we will consider $u(x,t)\in M$ , $t_1,t_2$ , and $t_1<t_2$ ; then for a given $\epsilon >0$ , if $|t_2-t_1|<\delta$ . Then

$$\begin{array}{c}\Vert \mathrm{Tu}(x,t_2)-\mathrm{Tu}(x,t_1)\Vert \hfill \\ \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}\Vert s(x,t_2,u(x,t_2))-s(x,t_1,u(x,t_1))\Vert \hfill \\ +\Vert \begin{array}{c}\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t_2}{\int }}{(t_2-\tau )}^{\alpha -1}\Vert s(x,\tau ,u(x,\tau ))\Vert d\tau \hfill \\ \multicolumn{1}{c}{-\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t_1}{\int }}{(t_1-\tau )}^{\alpha -1}\Vert s(x,\tau ,u(x,\tau ))\Vert d\tau }\end{array}\Vert \hfill \\ \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}\Vert s(x,t_2,u(x,t_2))-s(x,t_1,u(x,t_1))\Vert \hfill \\ +\frac{\alpha P}{\mathrm{AB}(\alpha )\Gamma (\alpha )}\{\underset{0}{\overset{t_2}{\int }}{(t_2-y)}^{\alpha -1}\mathrm{dy}-\underset{0}{\overset{t_1}{\int }}{(t_1-y)}^{\alpha -1}\mathrm{dy}\}\hfill \end{array}$$

(15)

We first start with the integral part

$$\begin{array}{c}\underset{0}{\overset{t_2}{\int }}{(t_2-y)}^{\alpha -1}\mathrm{dy}-\underset{0}{\overset{t_1}{\int }}{(t_1-y)}^{\alpha -1}\mathrm{dy}\hfill \\ =\underset{0}{\overset{t_1}{\int }}\{{(t_1-y)}^{\alpha -1}-{(t_2-y)}^{\alpha -1}\}\mathrm{dy}\hfill \\ +\underset{t_1}{\overset{t_2}{\int }}{(t_2-y)}^{\alpha -1}\mathrm{dy}=2{\frac{(t_2-t_1)}{\alpha }}^{\alpha }\hfill \end{array}$$

(16)

Now we will investigate the following

$$\Vert s\left(x,t_2,u(x,t_2)\right)-s\left(x,t_1,u(x,t_1)\right)\Vert =\Vert \begin{array}{c}u(x,t_2)\left[-\frac{\chi v(x,t_2)y(x,t_2)}{1+\chi hu(x,t_2)}\right]-u(x,t_1)\left[-\frac{\chi v(x,t_1)y(x,t_1)}{1+\chi hu(x,t_1)}\right]\\ +\frac{r_1}{k}\left(u(x,t_2)-u(x,t_1)\right)\left[1-u(x,t_2)-u(x,t_1)\right]\\ +({\rho }_1-1)\left[v_1(x,t_2){\partial }_{x}u(x,t_2)-v_1(x,t_1){\partial }_{x}u(x,t_1)\right]\\ -{\rho }_1{\varphi }_1\left[{\partial }_{t}u(x,t_2)-{\partial }_{t}u(x,t_1)\right]\\ -a\left[u(x,t_2)v(x,t_2)-u(x,t_1)v(x,t_1)\right]-b\left[u(x,t_2)-u(x,t_1)\right]\end{array}\Vert$$

(17)

Because all the solutions are bounded, let us find appropriate different positive constants, a₁, b₁, c₁, d₁, and e₁ for all t. Also, we can consider $[-(\chi v(x,t_2))/(1+\chi \mathrm{hu}(x,t_2))]<1$ and $[-(\chi v(x,t_1))/(1+\chi \mathrm{hu}(x,t_1))]$ . Finally, if we use Lipschitz condition of the derivative, equality (17) can be reconsidered as below

$$\begin{array}{c}\Vert s(x,t_2,u(x,t_2))-s(x,t_1,u(x,t_1))\Vert \hfill \\ \le \Vert u(x,t_2)-u(x,t_1)\Vert a_1+\frac{r_1}{k}{b}_1\Vert u(x,t_2)-u(x,t_1)\Vert \hfill \\ +({\rho }_1-1)c_1\Vert u(x,t_2)-u(x,t_1)\Vert -{\rho }_1{\varphi }_1{d}_1\Vert u(x,t_2)-u(x,t_1)\Vert \hfill \\ -a{e}_1\Vert u(x,t_2)-u(x,t_1)\Vert -b\Vert u(x,t_2)-u(x,t_1)\Vert \hfill \\ \le \underset{F_1}{\underbrace{\{a_1+\frac{r_1}{k}{b}_1+({\rho }_1-1)c_1-{\rho }_1{\varphi }_1{d}_1-a{e}_1-b\}}}\Vert u(x,t_2)-u(x,t_1)\Vert \hfill \\ \le F_1·n\Vert t_2-t_1\Vert \hfill \\ \le A\Vert t_2-t_1\Vert \hfill \end{array}$$

(18)

Now putting equations (16)–(18) in equation (15), we obtain

$$\begin{array}{c}\Vert \mathrm{Tu}(x,t_2)-\mathrm{Tu}(x,t_1)\Vert \hfill \\ \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}A\Vert t_2-t_1\Vert +\frac{\alpha P}{\mathrm{AB}(\alpha )\Gamma (\alpha )}2{\frac{\Vert t_2-t_1\Vert }{\alpha }}^{\alpha }\hfill \\ \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}A\Vert t_2-t_1\Vert +\frac{2\alpha P}{\mathrm{AB}(\alpha )\Gamma (\alpha +1)}\Vert t_2-t_1{\Vert }^{\alpha }\hfill \end{array}$$

(19)

$$\delta =\frac{\epsilon }{\frac{1-\alpha }{\mathrm{AB}(\alpha )}A+\frac{2\alpha P}{\mathrm{AB}(\alpha )\Gamma (\alpha +1)}}$$

(20)

Such that $\Vert \mathrm{Tu}(x,t_2)-\mathrm{Tu}(x,t_1)\Vert \le \epsilon$ is satisfied.

With same rules step by step, we can obtain the following for other two equations

$$\begin{array}{l}\Vert Tv(x,t_2)-Tv(x,t_1)\Vert \le \frac{1-\mu }{AB(\mu )}B\Vert t_2-t_1\Vert \\ +\frac{2\mu Q}{AB(\mu )\Gamma (\mu +1)}{\Vert t_2-t_1\Vert }^{\alpha }\end{array}$$

(21)

For each $\epsilon >0$ , we can find

$$\delta =\frac{\epsilon }{\frac{1-\mu }{\mathrm{AB}(\mu )}B+\frac{2\mu Q}{\mathrm{AB}(\mu )\Gamma (\mu +1)}}$$

(22)

and

$$\begin{array}{c}\Vert \mathrm{Ty}(x,t_2)-\mathrm{Ty}(x,t_1)\Vert \hfill \\ \le \frac{1-\theta }{\mathrm{AB}(\theta )}C\Vert t_2-t_1\Vert +\frac{2\mu V}{\mathrm{AB}(\theta )\Gamma (\theta +1)}\Vert t_2-t_1{\Vert }^{\alpha }\hfill \\ \multicolumn{1}{c}{\delta =\frac{\epsilon }{\frac{1-\theta }{\mathrm{AB}(\theta )}C+\frac{2\mu V}{\mathrm{AB}(\theta )\Gamma (\theta +1)}}}\end{array}$$

(23)

Consequently

$$\Vert \mathrm{Tv}(x,t_2)-\mathrm{Tv}(x,t_1)\Vert \le \epsilon$$

(24)

and

$$\Vert \mathrm{Ty}(x,t_2)-\mathrm{Ty}(x,t_1)\Vert \le \epsilon$$

(25)

are satisfied. So $T(M)$ is equi-continuous and with the help of Arzela–Ascoli theorem, $\overline{T(M)}$ is compact.

Theorem 2

$S:[u_1,u_2]\times [0,\infty )\to [0,\infty )$ be a continuous function and $S(x,t,.)$ increasing for each t in $[u_1,u_2]$ . Let us assume that one can find w₁ and w₁ satisfying $V(s)w_1\le S(x,t,w_1)$ , $V(s)w_2\ge S(x,t,w_2)$ , $0\le w_1(x,t)\le w_2(x,t),u_1\le t\le u_2$ . Then our new equation has a positive solution.

Proof

To investigate the presence of positive solution, we need fixed-point operator of T. With framework of Lemma 1, the considered operator $T:H\to H$ is completely continuous. Let us take two arbitrary densities of population of predator $v_1$ and $v_2$ in H satisfying $v_1\le v_2$ and also densities of population of prey $s_1$ and $s_2$ in H satisfying $s_1\le s_2$ , then we have that S is a positive function, then the following are satisfied

$$\begin{array}{cc}\multicolumn{1}{c}{T{v}_1(x,t)}& \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}\Vert s(x,t,v_1(x,t))\Vert \hfill \\ \multicolumn{1}{c}{}& +\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\alpha -1}\Vert s(x,\tau ,v_1(x,\tau ))\Vert d\tau \hfill \\ \multicolumn{1}{c}{}& \le T{v}_2(x,t)\hfill \end{array}$$

(26)

and

$$\begin{array}{cc}\multicolumn{1}{c}{\Vert T{s}_1(x,t)\Vert }& \le \frac{1-\mu }{\mathrm{AB}(\mu )}\Vert s(x,t,s_1(x,t))\Vert \hfill \\ \multicolumn{1}{c}{}& +\frac{\mu }{\mathrm{AB}(\mu )\Gamma (\mu )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\mu -1}\Vert s(x,\tau ,s_1(x,\tau ))\Vert d\tau \hfill \\ \multicolumn{1}{c}{}& \le T{s}_2(x,t)\hfill \end{array}$$

(27)

So the mapping T is increasing. With the help of conjecture, we have $T{w}_2\ge w_2$ and $T{w}_1\le w_1$ . So the operator $T:〈w_1,w_2〉\to 〈w_1,w_2〉$ is compact via the framework of Lemma 2 and continuous in view of Lemma 1. Since H is a normal cone of T.

Uniqueness of solution

To work on the uniqueness of the solution is a good idea for modeling which will make our problem more clear. So the uniqueness of the solution is presented as below

$$\begin{array}{c}\Vert T{v}_1(x,t)-T{v}_2(x,t)\Vert \hfill \\ \le \Vert \begin{array}{c}\frac{1-\alpha }{\mathrm{AB}(\alpha )}(s(x,t,v_1(x,t))-s(x,t,v_2(x,t)))\hfill \\ \multicolumn{1}{c}{\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\alpha -1}(s(x,\tau ,v_1(x,\tau ))-s(x,\tau ,v_2(x,\tau )))d\tau }\end{array}\Vert \hfill \\ \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}\Vert s(x,t,v_1(x,t))-s(x,t,v_2(x,t))\Vert \hfill \\ +\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}\underset{0}{\overset{t}{\int }}{(t-\tau )}^{\alpha -1}\Vert s(x,\tau ,v_1(x,\tau ))-s(x,\tau ,v_2(x,\tau ))\Vert d\tau \hfill \\ \le \frac{1-\alpha }{\mathrm{AB}(\alpha )}{F}_1\Vert v_1(x,t)-v_2(x,t)\Vert \hfill \\ +\frac{\alpha }{\mathrm{AB}(\alpha )\Gamma (\alpha )}{F}_1\underset{0}{\overset{t}{\int }}(t-\tau {)}^{\alpha -1}\Vert v_1(x,\tau )-v_2(x,\tau )\Vert d\tau \hfill \end{array}$$

(28)

$$\begin{array}{l}\Vert T{v}_1(x,t)-T{v}_2(x,t)\Vert \\ \le \left\{\frac{1-\alpha }{AB(\alpha )}{F}_1+\frac{\alpha F_1{b}^{\alpha }}{AB(\alpha )\Gamma (\alpha +1)}\right\}\Vert v_1(x,t)-v_2(x,t)\Vert \end{array}$$

(29)

And similarly

$$\begin{array}{l}\Vert T{s}_1(x,t)-T{s}_2(x,t)\Vert \\ \le \left\{\frac{1-\mu }{AB(\mu )}{F}_2+\frac{\mu F_2{b}^{\mu }}{AB(\mu )\Gamma (\mu +1)}\right\}\Vert s_1(x,t)-s_2(x,t)\Vert \end{array}$$

(30)

Therefore, if the following conditions hold

$$\{\frac{1-\alpha }{\mathrm{AB}(\alpha )}{F}_1+\frac{\alpha F_1{b}^{\alpha }}{\mathrm{AB}(\alpha )\Gamma (\alpha +1)}\}<1$$

(31)

and

$$\{\frac{1-\mu }{\mathrm{AB}(\mu )}{F}_2+\frac{\mu F_2{b}^{\mu }}{\mathrm{AB}(\mu )\Gamma (\mu +1)}\}<1$$

(32)

then mapping T is a contraction, which implies fixed point, and thus the model has a unique positive solution.

Numerical analysis

In this section, we will introduce the numerical solution of the nonlinear triadic model. First, we will obtain numerical approximation of model’s derivative, that is, Atangana–Baleanu derivative in Caputo sense. For some positive integer N, the grid sizes in time for finite difference technique are given by

$$k=\frac{1}{N}$$

(33)

The grid points in the time interval $[0,T]$ are marked with $t_n=\mathrm{nj},n=0,1,2,\dots ,\mathrm{TN}$ .

The value of function g at the grid point is $g_i=g(t_i)$ . Now we start finding discrete approximation of the Atangana–Baleanu derivative ( $t=t_n$ , $n>1$ )

$${{}_0{}^{\mathrm{ABC}}D}_t^{\alpha }(g(t_n))=\frac{B(\alpha )}{1-\alpha }\underset{0}{\overset{t_n}{\int }}\frac{\mathrm{dg}(y)}{\mathrm{dy}}{E}_{\alpha }\left[\frac{-\alpha }{1-\alpha }{(t_n-y)}^{\alpha }\right]\mathrm{dy}$$

(34)

Putting the first-order approximation equality on equation (34), the definition of ABC derivative can be rewritten as below

$$\begin{array}{c}{{}_0{}^{\mathrm{ABC}}D}_t^{\alpha }\left(g(t_n)\right)=\frac{B(\alpha )}{1-\alpha }{\displaystyle \sum _{j=0}^{n-1}}{\displaystyle \underset{t_j}{\overset{t_{j+1}}{\int }}\left(\frac{g^{j+1}-g^j}{\Delta t}+O(\Delta t)\right)}{E}_{\alpha }\left[\frac{-\alpha }{1-\alpha }{(t_n-y)}^{\alpha }\right]dy,{}_0^{ABC}{D}_t^{\alpha }\left(g\left(t_n\right)\right)\\ =\frac{B(\alpha )}{1-\alpha }{\displaystyle \sum _{j=0}^{n-1}\left(\frac{g^{j+1}-g^j}{\Delta t}+O(\Delta t)\right)}{\displaystyle \underset{t_j}{\overset{t_{j+1}}{\int }}{E}_{\alpha }\left[\frac{-\alpha }{1-\alpha }{(t_n-y)}^{\alpha }\right]dy}\end{array}$$

(35)

If we compute the integral expression of equation (36)

$$\underset{t_j}{\overset{t_{j+1}}{\int }}{E}_{\alpha }\left[\frac{-\alpha }{1-\alpha }{(t_n-y)}^{\alpha }\right]\mathrm{dy}$$

(36)

then equation (37) will be obtained via Mittag-Leffler definition as below

$$\begin{array}{c}\underset{t_j}{\overset{t_{j+1}}{\int }}{E}_{\alpha }\left[\frac{-\alpha }{1-\alpha }{(t_n-y)}^{\alpha }\right]=\underset{t_j}{\overset{t_{j+1}}{\int }}\sum _{p=0}^{\infty }\frac{{\left(\frac{\alpha }{1-\alpha }\right)}^p(-{(t_n-y)}^{\alpha p})}{\Gamma (\alpha p+1)}37\hfill \\ =\sum _{p=0}^{\infty }\underset{t_j}{\overset{t_{j+1}}{\int }}\frac{{\left(\frac{\alpha }{1-\alpha }\right)}^p(-{(t_n-y)}^{\alpha p})}{\Gamma (\alpha p+1)}\hfill \\ =\sum _{p=0}^{\infty }\frac{{\left(\frac{\alpha }{1-\alpha }\right)}^p}{(\alpha p+1)!}\underset{{\delta }_{\alpha \mathrm{pnj}}}{\underbrace{({(t_n-t_{j+1})}^{\alpha p+1}-{(t_n-t_j)}^{\alpha p+1})}}\hfill \\ =\sum _{p=0}^{\infty }\frac{{\left(\frac{\alpha }{1-\alpha }\right)}^p}{(\alpha p+1)!}{\delta }_{\alpha \mathrm{pnj}}\hfill \end{array}$$

(37)

So ABC derivative can be reconsidered like equation (37)

$${{}_0{}^{\mathrm{ABC}}D}_t^{\alpha }(g(t_n))=\frac{B(\alpha )}{1-\alpha }\sum _{j=0}^{n-1}\frac{g^{j+1}-g^j}{\Delta t}\sum _{p=0}^{\infty }\frac{{\left(\frac{\alpha }{1-\alpha }\right)}^p}{(\alpha p+1)!}{\delta }_{\alpha \mathrm{pnj}}$$

(38)

In the rest of the section, we will replace the model with equation (38) and first approximation of local derivative definition

$$\begin{array}{c}\frac{B(\alpha )}{1-\alpha }\sum _{k=0}^{j-1}\frac{u_i^{k+1}-u_i^k}{\Delta t}\sum _{p=0}^{\infty }\frac{{\left(\frac{\alpha }{1-\alpha }\right)}^p}{(\alpha p+1)!}{\delta }_{\alpha \mathrm{pnj}}\hfill \\ =\frac{-\chi v_i^j{y}_i^j\left\{\frac{u_i^{j+1}+u_i^j}{2}\right\}}{1+\chi h\left\{\frac{u_i^{j+1}+u_i^j}{2}\right\}}+r_1\left\{\frac{u_i^{j+1}+u_i^j}{2}\right\}(1-\frac{\left\{\frac{u_i^{j+1}+u_i^j}{2}\right\}}{k})\hfill \\ -(1-{\rho }_1)v_{1i}^j\left\{\frac{(u_{i+1}^{j+1}-u_{i-1}^{j+1})+(u_{i+1}^j-u_{i-1}^j)}{4\Delta x}\right\}\hfill \\ -{\rho }_1{\Phi }_1\left\{\frac{(u_{i+1}^{j+1}-u_{i-1}^{j+1})+(u_{i+1}^j-u_{i-1}^j)}{4\Delta t}\right\}\hfill \\ -{\mathrm{av}}_i^j\left\{\frac{u_i^{j+1}+u_i^j}{2}\right\}\hfill \\ -b\left\{\frac{u_i^{j+1}+u_i^j}{2}\right\}\hfill \end{array}$$

(39)

$$\begin{array}{c}\frac{B(\gamma )}{1-\gamma }\sum _{k=0}^{j-1}\frac{v_i^{k+1}-v_i^k}{\Delta t}\sum _{p=0}^{\infty }\frac{{\left(\frac{\gamma }{1-\gamma }\right)}^p}{(\gamma p+1)!}{\delta }_{\gamma \mathrm{pnj}}\hfill \\ =\frac{-\chi u_i^j{y}_i^j\left\{\frac{v_i^{j+1}+v_i^j}{2}\right\}}{1+\chi h\left\{\frac{v_i^{j+1}+v_i^j}{2}\right\}}+r_2\left\{\frac{v_i^{j+1}+v_i^j}{2}\right\}(1-\frac{\left\{\frac{v_i^{j+1}+v_i^j}{2}\right\}}{k})\hfill \\ -(1-{\rho }_2)v_{2i}^j\left\{\frac{(v_{i+1}^{j+1}-v_{i-1}^{j+1})+(v_{i+1}^j-v_{i-1}^j)}{4\Delta x}\right\}\hfill \\ -{\rho }_2{\Phi }_2\left\{\frac{(v_{i+1}^{j+1}-v_{i-1}^{j+1})+(v_{i+1}^j-v_{i-1}^j)}{4\Delta t}\right\}\hfill \\ -{\mathrm{cu}}_i^j\left\{\frac{v_i^{j+1}+v_i^j}{2}\right\}\hfill \\ -d\left\{\frac{v_i^{j+1}+v_i^j}{2}\right\}\hfill \end{array}$$

(40)

$$\begin{array}{c}\frac{B(\theta )}{1-\theta }\sum _{k=0}^{j-1}\frac{y_i^{k+1}-y_i^k}{\Delta t}\sum _{p=0}^{\infty }\frac{{\left(\frac{\theta }{1-\theta }\right)}^p}{(\theta p+1)!}{\delta }_{\theta \mathrm{pnj}}\hfill \\ =\frac{\beta v_i^j{u}_i^j\left\{\frac{y_i^{j+1}+y_i^j}{2}\right\}}{1+\chi {\mathrm{hu}}_i^j{v}_i^j}-f\left\{\frac{y_i^{j+1}+y_i^j}{2}\right\}\hfill \\ -g\left\{\frac{y_i^{j+1}+y_i^j}{2}\right\}{v}_i^j\hfill \\ -(1-{\rho }_3)v_{3i}^j\left\{\frac{(y_{i+1}^{j+1}-y_{i-1}^{j+1})+(y_{i+1}^j-y_{i-1}^j)}{4\Delta x}\right\}\hfill \\ -{\rho }_3{\Phi }_3\left\{\frac{(y_{i+1}^{j+1}-y_{i-1}^{j+1})+(y_{i+1}^j-y_{i-1}^j)}{4\Delta t}\right\}\hfill \\ -e\left\{\frac{y_i^{j+1}+y_i^j}{2}\right\}{u}_i^j\hfill \end{array}$$

(41)

Conclusion

This work addresses a new mathematical model able to describe the dynamical biology of three species where two are considered as predator and prey but is only acting as prey. The mathematical model is built using a new fractional differentiation that has derivative with the generalized Mittag-Leffler function as kernel due to the non-locality of the physical problem and also the wider applicability of the generalized Mittag-Leffler function. An analytical study was undertaken to prove the existence of a unique set of the solutions of the new model using the fixed-point theorem.