Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 03 November 2022
Sec. Mathematical Biology
This article is part of the Research Topic Management and Control of Infectious Disease Outbreak: From Models to Policies View all 8 articles

Prey fear of a specialist predator in a tri-trophic food web can eliminate the superpredator

\nNabaa Hassain FakhryNabaa Hassain Fakhry1Raid Kamel NajiRaid Kamel Naji1Stacey R. Smith?Stacey R. Smith?2Mainul Haque
Mainul Haque3*
  • 1Department of Mathematics, College of Science, University of Baghdad, Baghdad, Iraq
  • 2Department of Mathematics and Faculty of Medicine, The University of Ottawa, Ottawa, ON, Canada
  • 3School of Mathematical Science, The University of Nottingham Ningbo China, Ningbo, China

We propose an intraguild predation ecological system consisting of a tri-trophic food web with a fear response for the basal prey and a Lotka–Volterra functional response for predation by both a specialist predator (intraguild prey) and a generalist predator (intraguild predator), which we call the superpredator. We prove the positivity, existence, uniqueness, and boundedness of solutions, determine all equilibrium points, prove global stability, determine local bifurcations, and illustrate our results with numerical simulations. An unexpected outcome of the prey's fear of its specialist predator is the potential eradication of the superpredator.

Introduction

A functioning ecosystem depends on the framework of food webs that it supports [1]. Food webs involve many types of predator–prey interactions, which are fundamental interactions for sustaining species [2]. Interaction between the prey and the predator can be affected by factors such as refuge, disease, stage structure, competition, and fear [37]. Intraguild predation (IGP), on the other hand, is described as predator–prey interactions among consumers who may be fighting for limited resources. In natural communities, there is a growing body of literature emphasizing the relevance of intraguild predation [8, 9]. Three species are involved in the simplest intraguild predation model: a superpredator (IG predator), a specialist predator (IG prey), and a basal prey. Bai et al. recently suggested a three-species IGP food web model with the IG predator, IG prey, and basal prey, in which the basal prey grows logistically with a large Allee effect [10]. They looked into the model's local and global dynamics, focusing on the impact of the Allee effect and discovered that the intraguild predation food web model has rich and complicated dynamic behavior and that a large Allee effect in the basal prey raises the danger of extinction for not only the basal prey but also the IG prey or/and IG predator.

Many predators in food webs are superpredators, who may not restrict their diets to a specific prey species but feed also on other predators [11]. Therefore, superpredators are expected to compete not only with other predators for food and space but in many cases also through intraguild predation [5, 12, 13]. Predators induce indirect effects such as fear in prey that can change the prey's behavior [14]. Fear takes the form of sustained psychological stress on the prey, as prey species are always wary of possible attack [15]. Suraci et al. experimentally showed that fear of large carnivores reduces mesocarnivore foraging, which benefits the mesocarnivore's prey [16].

A handful of mathematical models have studied the effect of fear on food webs. Panday et al. investigated the impact of fear in a tri-trophic food chain model, with prey fear in response to both predators [15], from the middle predator to the generalist predator [17] and with delays [18]. Cong et al. introduced a fear-adjusted birth rate to a three-species food web [19]. Hossain et al. limited the growth of the prey due to fear in an intraguild predation model [20]. Mukerjee incorporated interspecific competition and fear affecting the death rate of the prey [21]. Ibrahim et al. limited prey growth due to fear of the generalist predator and showed that fear could have a stabilizing effect on the system [22]. Mondal et al. showed that the prey's fear of predators was responsible for the increase in intraspecific competition among the prey species [23]. Roy et al. showed that fear could play a destablizing role if it caused a reduction in the birth rate of susceptible prey, whereas the levels of fear responsible for the increase in the intraspecies competition of susceptible prey and eradication of the disease prevalence could stabilize an otherwise unstable system [24]. Maity et al. considered time-varying fear effects, showing that periodic solutions could arise [25]. Hossain et al. showed that perceived fear of predators could reduce the prey birth rate [26]. Tiwari et al. showed that seasonal variations in the level of prey fear generated higher-order periodic solutions [27].

Here, we examine the impact of fear on the dynamical behavior of an IGP food web system in which the prey responds with fear to the specialist predator but not the superpredator; because the superpredator has alternative food sources, it will be less dangerous to the basal prey compared to the IG prey. To the best of our knowledge, this is the first model to examine this effect.

Food-web model formulation

Our food web consists of a prey at the first level, a specialist predator at the second level, and a superpredator at the third level. Let x(t), y(t), and z(t) be the population densities at time t for the prey, specialist predator, and superpredator, respectively. The prey grows logistically in the absence of predators, while it has a fear property of predation in the presence of the specialist predator. Hence, the intrinsic growth rate of prey becomes r1+ky, which is a monotonic decreasing function of both k and y, where k represents the fear rate [28]. The food transport attack rates are given by the parameters a1, a2, and a3, with conversion rates e1, e2, and e3. Finally, the predators face natural death rates d1 and d2 for the specialist predator and superpredator, respectively. The dynamics of the food web with fear can be represented mathematically by the following set of differential equations:

dxdt=rx(1-x)(11+ky)-a1xy-a2xz=xf1(x,y,z),dydt=e1a1xy-a3yz-d1y=yf2(x,y,z),                     dzdt=e2a2xz+e3a3yz-d2z=zf3(x,y,z).                      (1)

Initial conditions satisfy x(t) ≥ 0, y(t) ≥ 0, and z(t) ≥ 0, and all parameters are assumed to be positive.

The interaction functions fi (i = 1, 2, 3) are continuous and have continuous partial derivatives and are thus Lipschitzian, so system (1) has a unique solution. Furthermore, for any initial condition in +3, the solution of system (1) is positive and uniformly bounded as shown in the following theorem. Hence, system (1) will be a dissipative system.

Theorem (1): The domain of system (1), +3, is positively invariant, and all solutions of system (1) starting in +3 are uniformly bounded.

Proof. Let x(t), y(t), and z(t) be any solution of system (1). Since the solution (x(t), y(t), z(t)) of the system (1) with initial condition in +3 exists and is unique on [0 , δ), where 0 < δ ≤ +∞, we have:

x(t)=x(0)e0t[r(1-x(s))(11+ky(s))-a1y(s)-a2z(s)]ds0
y(t)=y(0)e0t[e1a1x(s)-a3z(s)-d1]ds0
z(t)=z(0)e0t[e2a2x(s)+e3a3y(s)-d2]ds 0.

From the first equation of system (1):

dxdtrx(1-x ).

Then, it is easy to verify that x(t)r4 for all t.

Let = x + y + z. Then, since ei ∈ (0, 1); i = 1, 2, 3, we have

dWdtr4-M[x+y+z ],

where M = min{1, d1, d2}. Solving the following linear differential inequality

dWdt+MW r4,

we obtain that, as t → ∞,

W(t)r4M.

Existence of equilibrium points

System (1) has at most five nonnegative biologically feasible equilibrium points. The trivial equilibrium E0 = (0, 0, 0) always exists. The axial equilibrium E1 = (1, 0, 0) always exists on the boundary of the first octant. The specialist-free equilibrium point, E2=(x̄,0,z̄), where

x̄=d2e2a2 and z̄=r(1-x̄)a2,    (2a)

exists provided that

d2<e2a2.    (2b)

The superpredator-free equilibrium point, E3=(x^,y^,0), where

x^=d1e1a1 and y^=-12k+12a1ka12+4kr(e1a1-d1)e1,    (3a)

exists, provided that

d1<e1a1.    (3b)

There is also an interior equilibrium point, E4=(x*,y*,z*), where

y*=d2-e2a2x*e3a3 and z*=e1a1x*-d1a3,    (4a)

with x* a positive root of the following second-order polynomial equation:

β1x2+β2x+β3=0,    (4b)

where

β1=e2a1a22k[e1e3-e2], 
β2=e2e3a1a2a3-re32a32+2e2a1a2kd2-e1e3a1a2a32-e1e3a1a2kd2-e2e3a22kd1, 
β3=re32a32-e3a1a3d2-a1kd22+e3a2a32d1+e3a2kd1d2.

Therefore, there is a unique interior equilibrium point in the interior of +3 provided that the following conditions hold:

d1e1a1<x*<d2e2a2    (4c)
 β1>0;β3<0ORβ1<0;β3>0}.    (4d)

Persistence

Next, we determine the requirements that ensure persistence in the system (1). Because the types of attractors available in the boundary planes affect the creation of our IGP food web model's persistence conditions, an analysis of the dynamics in the boundary planes of the IGP food web system is conducted. System (1) has two subsystems: the first occurs in the absence of the superpredator (IG predator) and the second occurs in the absence of the specialist predator (IG prey). The first subsystem is

dxdt=x[r(1-x)1+ky-a1y]=xf1(x,y)dydt=y[e1a1x-d1]=yg1(x,y).    (5)

The second subsystem can be written as follows:

dxdt=x[r(1-x)-a2z]=xf2(x,z)dzdt=z[e2a2x-d2]=zg2(x,z).          (6)

Subsystem (5) has a unique positive equilibrium point in the interior of xy-plane given by (x^,y^)E3, while subsystem (6) has a unique positive equilibrium point in the interior of xz-plane given by (x̄,z̄)E2. The interior equilibrium point of subsystem (5) corresponds to system (1)'s superpredator-free equilibrium point; similarly, the interior equilibrium point of subsystem (6) corresponds to system (1)'s specialist-predator-free equilibrium point. The dynamics around these equilibrium points can be described in the following theorem.

Theorem (2): There are no periodic dynamics in the interior of xy-plane or the xz-plane.

Proof. Define a continuously differential function B(x,y)=1xy. Then, we have:

Δ=(Bf1)x+(Bg1)y=-ry(1+ky). 

Clearly, Δ has the same sign and does not equal zero almost everywhere in a simply connected region of the xy-plane. By the Dulac–Bendixson criterion, system (1) has no periodic solutions lying entirely in the interior of xy−plane. The second part follows by using the Dulac function C(x,z)= 1xz.

Theorem (3): System (1) is uniformly persistent in the interior of +3, provided that

e1a1>d1,    (7a)
e2a2>d2,    (7b)
e1a1x̄>a3z̄+d1,    (7c)
e2a2x^+e3a3y^>d2.    (7d)

Proof. Consider the function (x,y,z)=xb1yb2zb3, where bi, i = 1, 2, 3, are positive constants and U(x, y, z) is a C1 nonnegative function in the interior of +3. Hence, we have

dUdt=Uxdxdt+Uydydt+Uzdzdt,

with

Ux=b1x(b1-1 )yb2zb3,
Uy=b2xb1y(b2-1 )zb3,
Uz=b3xb1yb2z(b3-1 ).

Therefore,

Ψ(x,y,z)=(dUdt)U(x,y,z)=b1xdxdt+b2ydydt+b3zdzdt,                    =b1[r(1-x)(11+ky)-a1y-a2z]                  +b2[e1a1x-a3z-d1]+b3[e2a2x+e3a3y-d2].

Since there are no periodic solutions in the boundary planes of the system (1), it follows that system (1) is uniformly persistent provided that Ψ(Ei) > 0 for each i = 1, 2, 3. We have

Ψ(E1)=b2[e1a1-d1]+b3[e2a2-d2 ],
Ψ(E2)=b2[e1a1x̄-a3z̄-d1 ],
Ψ(E3)=b3[e2a2x^+e3a3y^-d2 ].

Consequently, conditions (7a)–(7d) satisfy Ψ(Ei) > 0 for each i = 1, 2, 3.

Global stability analysis

Here, we use Lyapunov functions to investigate the global stability of equilibria.

Theorem (4): The axial equilibrium point E1 = (1, 0, 0) is globally asymptotically stable provided that the following sufficient condition holds:

e2<e1e3<d2a2.    (8)

Proof. Consider the following positive-definite, real-valued function around E1:

V1=c1[x-1-ln x]+c2y+c3z .

Then, we have:

dV1dt=-c1r1+ky(x-1)2-(c1-c2e1)a1xy-(c1-c3e2)a2xz      -(c2-c3e3)a3yz-(c2d1-c1a1)y-(c3d2-c1a2)z.

Then, by choosing c1 = e1, c2 = 1 and c3=1e3, we obtain

dV1dt=-e1r1+ky(x-1)2-(e1-e2e3)a2xz-(d1-e1a1)y               -(d2e3-e1a2 )z.

Clearly, under condition (8), dV1dt is negative definite. Moreover, since V1 is radially unbounded, the axial equilibrium point E1 = (1, 0, 0) is globally asymptotically stable.

Theorem (5): The specialist-free equilibrium point E2=(x̄,0,z̄) is globally asymptotically stable provided that the following sufficient conditions hold:

a1x̄<d1e1+e3a3e2z̄,    (9a)
e1e3<e2.    (9b)

Proof. Consider the following positive-definite, real-valued function around E2:

V2=c1[x-x̄-x̄lnxx̄]+c2y+c3[z-z̄-z̄lnzz̄ ].

Then, we have:

dV2dt=-c1r(x-x̄)2-a2[c1-c3e2](x-x̄)(z-z̄)                              -(c1-c2e1)a1xy-(c2d1-c1a1x̄+c3e3a3z̄)y-(c2-c3e3)a3yz. 

Then, by choosing c1 = 1, c2=1e1 and c3=1e2, we have:

dV2dt=-r(x-x̄)2-(d1e1+e3a3e2z̄-a1x̄)y-(1e1-e3e2 )a3yz.

Obviously, under conditions (9a)–(9b), dV2dt is negative definite. Moreover, since V2 is radially unbounded, the specialist-free equilibrium point E2=(x̄,0,z̄) is globally asymptotically stable.

Theorem (6): The superpredator-free equilibrium point E3=(x^,y^,0)is globally asymptotically stable, provided that, in addition to condition (9b), the following sufficient conditions hold:

a2x^+a3e1y^<d2e2,    (10a)
(x-x^)(y-y^)>0.    (10b)

Proof. Consider the following positive-definite, real-valued function around E3:

V3=c1[x-x^-x^lnxx^]+c2[y-y^-y^lnyy^]+c3z.

Then, we have

dV3dt=-c1r(x-x^)2R-c1rk(1-x^)(x-x^)(y-y^)RR^               -(c1-c2e1)a1(x-x^)(y-y^)-(c1-c3e2)a2xz                   -(c2-c3e3)a3yz-(c3d2-c1a2x^-c2a3y^)z, 

where R = (1 + ky) and R^=(1+ky^). Then, by choosing c1 = 1, c2=1e1 and c3=1e2, we obtain:

dV3dt-r(x-x^)2R-rk(1-x^)(x-x^)(y-y^)RR^          -(1e1-e3e2)a3yz-(d2e2-a2x^-a3e1y^ )z.

Under conditions (10a)–(10b), dV3dt is negative definite. Since V3 is radially unbounded, the superpredator-free equilibrium point E3=(x^,y^,0) is globally asymptotically stable.

Theorem (7): The coexistence equilibrium point E4=(x*,y*,z*) is globally asymptotically stable provided that, in addition to condition (9b), the following sufficient condition holds:

(x-x*)(y-y*)>0.    (11)

Proof. Consider the following positive-definite, real-valued function around E4:

V4=c1[x-x*-x*lnxx*]+c2[y-y*-y*lnyy*]                                                     + c3[z-z*-z*lnzz* ].

Then, we have:

dV4dt=-c1r(x-x*)2R-c1rk(1-x*)(x-x*)(y-y*)RR*       -(c1-c2e1)a1(x-x*)(y-y*)          -(c2-c3e3)a3(y-y*)(z-z*)                     -(c1-c3e2)a2(x-x*)(z-z*), 

where R* = 1 + ky*. Then, by choosing c1 = e2, c2 = e3, and c3 = 1, we have

dV4dt=-e2r(x-x*)2R-e2rk(1-x*)(x-x*)(y-y*)RR*                                              -(e2-e1e3)a1(x-x*)(y-y* ).

Then, under condition (11) with (9b), dV4dt is negative definite. Since V4 is radially unbounded, the coexistence equilibrium point E4(x*,y*,z*) is globally asymptotically stable.

Bifurcation analysis

In this section, we investigate the effect of varying parameter values on the dynamics of the system (1) using local bifurcation analysis and the Sotomayor theorem [29]. First, we rewrite system (1) in the vector form as follows:

dXdt=F(X),X=(x,y,z)Tand F=(xf1,yf2,zf3 )T.

The second derivative of F with respect to X can be written as:

D2F(U,U)=(-2ru121+ky-2kr(1-2x)u1u2(1+ky)2-2a1u1u2-2a2u1u3+2k2rx(1-x)u22(1+ky)32e1a1u1u2-2a3u2u32e2a2u1u3+2e3a3u2u3),    (12)

where U=(u1,u2,u3)T is a general vector.

Theorem (8): Assume that d2=e2a2 (d2*). Then system (1) undergoes a transcritical bifurcation at the axial equilibrium point E1, but neither a saddle node nor a pitchfork bifurcation can occur if

e1a1<d1.    (13)

Proof. The Jacobian matrix of system (1) at the axial equilibrium point E1 with d2=e2a2(d2*) can be written as:

J1=J(E1,d2*)=[-r      -a1      -a20      e1a1-d1      00            0               0 ].

J1 has eigenvalues λ11*=-r<0, λ12*=e1a1-d1<0 under condition (13), and λ13*=0. Hence, the necessary but not sufficient condition for a bifurcation is satisfied, and E1 is a non-hyperbolic point.

Let Φ1=(v11,v12,v13)T be the eigenvector of J1 corresponding to the eigenvalue λ13*=0. Straightforward computation gives Φ1=(β1v13,0,v13)T, where v13 represents any nonzero real number and β1=-a2r< 0.

Let Ψ1=(ψ11,ψ12,ψ13)T be the eigenvector of J1T corresponding to the eigenvalue λ13*=0. Direct calculation shows that Ψ1=(0,0,ψ13)T, where ψ13 is any nonzero real number. Because Fd2=Fd2=(0,0,-z)T, we obtain that Fd2(E1,d2*)=(0,0,0)T, which yields:

Ψ1T[Fd2(E1,d2*)]= 0.

By Sotomayor's theorem, system (1) at E1 with d2=d2* does not experience a saddle-node bifurcation. Moreover, we have

Ψ1T[DFd2(E1,d2*)Φ1]=-v13Ψ13 0,

where DFd2 represents the derivative of Fd2 with respect to X. Applying equation (12) at (E1,d2*) with the eigenvector Φ1, we obtain that

Ψ1T[D2F(E1,d2*)(Φ1,Φ1)]=2e2a2β1v132ψ13 0.

By Sotomayor's theorem, system (1) near the equilibrium point E1 with d2=d2* undergoes a transcritical bifurcation, but a pitchfork bifurcation cannot occur.

Theorem (9): System (1) at the specialist-free equilibrium point E2 undergoes a transcritical bifurcation when d1=e1a1x̄+a3z̄ (d1*), but neither a saddle-node nor a pitchfork bifurcation can occur, provided

e1e3a1a2+re3a3e2a2 [kr(1-x̄)+a1].    (14)

Proof. The Jacobian matrix of system (1) at the specialist-free equilibrium point E2 with d1=d1* takes the form:

J2=J(E2,d1*)=[-rx̄-krx̄(1-x̄)-a1x̄-a2x̄000e2a2z̄e3a3z̄0 ].

J2 has eigenvalues

λ21* =-rx̄2+12(rx̄)2-4e2a22x̄z ̄,λ23* =-rx̄2-12(rx̄)2-4e2a22x̄z̄,

while λ22*=0; hence, the necessary but not sufficient condition for a bifurcation is satisfied, and E2 is a non-hyperbolic point.

Let Φ2=(v21,v22,v23)T be the eigenvector of J2 corresponding to the zero eigenvalue. Straightforward computation gives Φ2=(α1v22,v22,α2v22)T, where v22 represents any nonzero real number, with

α1=-e3a3e2a2 and α2=re3a3-[kr(1-x̄)+a1]e2a2e2a22.

Let Ψ2=(ψ21,ψ22,ψ23)T be the eigenvector of J2T corresponding to the zero eigenvalue. Direct calculation shows that Ψ2=(0,ψ22,0)T, where ψ22 is any nonzero real number. Because Fd1=Fd1=(0,-y,0)T, we obtain that Fd1(E2,d1*)=(0,0,0)T, which yields

Ψ2T[Fd1(E2,d1*)]= 0.

Thus, by Sotomayor's theorem, system (1) at E2 with d1=d1* does not experience a saddle-node bifurcation. Moreover, we have:

Ψ2T[DFd1(E2,d1*)Φ2]=-v22ψ22 0,

where DFd1 represents the derivative of Fd1 with respect to X. Using equation (12) at (E2,d1*) with the eigenvector Φ2, we obtain that

Ψ2T[D2F(E2,d1*)(Φ2,Φ2)]=(2e1a1α1-2a3α2 )v222ψ22.

It is easy to verify that Ψ2T[D2F(E2,d1*)(Φ2,Φ2)]0 due to condition (14). Hence, by Sotomayor's theorem, system (1) near the equilibrium point E2 with d1=d1* undergoes a transcritical bifurcation but a pitchfork bifurcation cannot occur.

Theorem (10): At the superpredator-free equilibrium point E3, system (1) undergoes a transcritical bifurcation when d2=e2a2x^+e3a3y^ (d2*), but neither a saddle-node nor a pitchfork bifurcation can occur, provided

e2a2[kr(1-x^)+a1(1+ky^)2]e3(1+ky^)[a3r+e1a1a2(1+ky^)].    (15)

Proof. The Jacobian matrix of system (1) at the superpredator-free equilibrium point E3 with d2=d2* is

J3=J(E3,d2*)=[-rx^1+ky^-krx^(1-x^)(1+ky^)2-a1x^-a2x^e1a1y^0-a3y^000 ].

J3 has eigenvalues

λ31* =-rx^2(1+ky^)+ 12(rx^1+ky^)2-4e1a1y^(krx^(1-x^)(1+ky^)2+a1x^ ),
λ32*=-rx^2(1+ky^)- 12(rx^1+ky^)2-4e1a1y^(krx^(1-x^)(1+ky^)2+a1x^ ),

while λ33*=0; hence, the necessary but not sufficient condition for bifurcation is satisfied, and E3 is a non-hyperbolic point.

Let Φ3=(v31,v32,v33)T be the eigenvector of J3 corresponding to the zero eigenvalue. Straightforward computation gives Φ3=(s1v33,s2v33,v33)T, where v33 represents any nonzero real number, while

s1=a3e1a1 and s2=-[ra3(1+ky^)+e1a1a2(1+ky^)2e1a1[kr(1-x^)+a1(1+ky^)2] ].

Let Ψ3=(ψ31,ψ32,ψ33)T be the eigenvector of J3T corresponding to the zero eigenvalue. Direct calculation shows that Ψ3=(0,0,ψ33)T, where ψ33 is any nonzero real number. Because Fd2=Fd2=(0,0,-z)T, we obtain Fd2(E3,d2*)=(0,0,0)T, which yields

Ψ3T[Fd2(E3,d2*)]= 0.

Thus, by Sotomayor's theorem, system (1) at E3 with d2=d2* does not experience a saddle-node bifurcation. Moreover, we have:

Ψ3T[DFd2(E3,d2)Φ3]=-v33ψ33 0,

where DFd2 represents the derivative of Fd2 with respect to X. By using equation (12) at (E3,d2*) with the eigenvector Φ3, we obtain that:

Ψ3T[D2F(E3,d2*)(Φ3,Φ3)]=(2e2a2s1+2e3a3s2 )v332ψ33.

Straightforward computation shows that Ψ3T[D2F(E3,d2*)(Φ3,Φ3)]0 due to condition (15). Thus, by Sotomayor's theorem, system (1) near the equilibrium point E3 with d2=d2* undergoes a transcritical bifurcation, but a pitchfork bifurcation cannot occur.

Theorem (11): System (1) undergoes a saddle-node bifurcation at the coexistence equilibrium point E4, but neither a transcritical nor a pitchfork bifurcation can occur when e3 passes through the value e3*=b12b23b31(b11b23-b13b21)a3z*, provided the following conditions hold:

(krx*(1-x*)(1+ky*)2+a1x*)e2a2>(rx*1+ky*)e3*a3,    (16a)
a1a2a3(e2e3*-e1)-r(1+ky*)×[1-k(1-2x*)e2a2(1+ky*)e3*a3-k2x*(1-x*)(e1a1)2(1+ky*)2a32]0.    (16b)

Proof. Direct computation shows that system (1) at the coexistence equilibrium point and e3=e3* has Jacobian matrix

J4= J(E4,e3*)=[-rx*1+ky*-krx*(1-x*)(1+ky*)2-a1x*-a2x*e1a1y*0-a3y*e2a2z*e3*a3z*0]                                                             =[bij]3 ×3.

Straightforward computation shows that |J4| = 0, under condition (16a). Hence, J4 has two eigenvalues with negative real parts, and the third one is λ* = 0. It follows that E4 becomes a non-hyperbolic equilibrium point.

Let Φ4=(v41,v42,v43)T be the eigenvector of J4 corresponding to the eigenvalue λ* = 0. Straightforward computation gives Φ4=(v41,δ1v41,δ2v41)T, where v41 represents any nonzero real number, δ1=-b31b32*<0 and δ2=-b21b23> 0.

Let Ψ4=(ψ41,ψ42,ψ43)T be the eigenvector of J4T corresponding to the eigenvalue λ* = 0. Direct calculation shows that Ψ4=(ψ41,μ1ψ41,μ2ψ41)T, where ψ33 is any nonzero real number, with μ1=-b13b23<0 and μ2=-b12b32> 0.

Since Fe3=Fe3=(0,0,a3yz)T, we have Fe3(E4,e3*)=(0,0,a3y*z*)T, which yields

Ψ4T[Fe3(E4,e3*)]=a3μ2y*z*ψ41 0.

By Sotomayor's theorem, transcritical and pitchfork bifurcations cannot occur while the first condition of the saddle-node bifurcation is satisfied. Moreover, from equation (12) with E4,e3* and Φ4, we obtain that:

D2F(E4,e3*)(Φ4,Φ4)=2v412×(-r1+ky*-kr(1-2x*)δ1(1+ky*)2-a1δ1-a2δ2+k2rx*(1-x*)δ22(1+ky*)3e1a1δ1-a3δ1δ2e2a2δ2+e3*a3δ1δ2). 

Hence,

Ψ4T[D2F(E4,e3*)(Φ4,Φ4)]=2v412ψ41×[-r1+ky* -kr(1-2x*)δ1(1+ky*)2-a1δ1-a2δ2 +k2rx*(1-x*)δ22(1+ky*)3+e1a1δ1μ1-a3δ1δ2μ1+e2a2δ2μ2+ e3*a3δ1δ2μ2].

Further computation shows that

Ψ4T[D2F(E4,e3*)(Φ4,Φ4)]=2v412ψ41[a1a2a3(e2e3*-e1)        -r(1+ky*)[1-k(1-2x*)e2a2(1+ky*)e3*a3-k2x*(1-x*)(e1a1)2(1+ky*)2a32]]. 

Using condition (16b), we have Ψ4T[D2F(E4,e3*)(Φ4,Φ4)]0. Hence, system (1) has a saddle-node bifurcation at E4 when e3=e3*.

Theorem (12): System (2) undergoes a Hopf bifurcation around E4 when e3=e3**, if and only if the following conditions hold:

(krx*(1-x*)(1+ky*)2+a1x*)e2a2<(rx*1+ky*)e3a3,    (17a)
[A(e3** )B(e3** )]< C(e3** ),    (17b)

where

A=-b11> 0,
B=-b12b21-b13b31-b23b32*> 0,
C=b23(b11b32*-b12b31)-b13b21b32*, 

with b32*=e3**a3z* and,

e3**=1e1a1a2a3y*z*[[kr(1-x*)(1+ky*)2-a1]×(e1a1y*(rx*1+ky*)+e2a2a3y*z*)+e2a22z*(rx*1+ky*)].

Proof. Consider the Jacobian matrix J4 given in Theorem (11), where e3=e3**. It is simple to determine that the characteristic equation is

λ3+Aλ2+Bλ+C=0.    (18)

Obviously, condition (17a) guarantees that C > 0. A straightforward computation shows that AB = C when e3=e3**. Hence, the characteristic equation (18) becomes

(λ2+B)(λ+A)=0.    (19)

Consequently, we obtain that λ1 = −A and λ2,3=±iB=±iδ2(e3**). Therefore, the Jacobian matrix has one negative real eigenvalue and two pure imaginary complex conjugates when e3=e3**. As a result, the first criterion for having a Hopf bifurcation is met.

Moreover, where e3 belongs to the neighborhood of e3**, the eigenvalues become λ2,3 = δ1(e3) ± 2(e3).

Next, we substitute δ1(e3) + 2(e3) into equation (18) and take the derivative of the resulting equation with respect to e3. Equating their real and imaginary parts, we obtain that

H1(e3)δ1(e3)-H2(e3)δ2(e3)=-H3(e3),H2(e3)δ1(e3)+H1(e3)δ2(e3)=-H4(e3),    (20)

where

H1(e3)=3[δ1(e3)]2+2A(e3)δ1(e3)-3[δ2(e3)]2+B(e3).
H2(e3)=6δ1(e3)δ2(e3)+2A(e3)δ2(e3 ).
H3(e3)=A(e3)[δ1(e3)]2-A(e3)[δ2(e3)]2+B(e3)δ1(e3)+C(e3 ).
H4(e3)=2A(e3)δ1(e3)δ2(e3)+B(e3)δ2(e3 ).

Solving system (20) for δ1(e3) gives

δ1(e3)=-H1(e3)H3(e3)+H2(e3)H4(e3)[H2(e3)]2+(e3)[H2(e3) ]2.

The result follows if and only if δ1(e3)0 or, equivalently, H1(e3)H3(e3) + H2(e3)H4(e3) ≠ 0 when e3=e3**. Note that we have the following:

δ1(e3**)=0 and δ2(e3**)=B(e3** ).
H1(e3**)=-2B(e3** ).
H2(e3**)=2A(e3**)B(e3** ).
H3(e3**)=-A(e3**)[B(e3**)]+C(e3** ).
H4(e3**)=B(e3**)B(e3** ).

Consequently,

H1(e3**)H3(e3**)+H2(e3**)H4(e3**)                                               =2A(e3**)B2(e3**)-2B(e3**)                                               C(e3**)+2A(e3**)B(e3**)B(e3**)                                               =-2B(e3**)[C(e3**)-(A(e3**)                                               B(e3**)+A(e3**)B(e3**)].      

Condition (17b) ensures that H1(e3**)H3(e3**)+H2(e3**)H4(e3**) 0.

Hence, the system has a Hopf bifurcation because δ1(e3)>0 under the conditions (17a)–(17b).

Numerical simulations

To illustrate the global dynamics of the system and confirm our analytical findings, we numerically simulated a hypothetical set of parameter values. Consider the following set of parameters:

r=1,k=1,a1=1,a2=0.5,a3=0.75,e1=0.75        e2=0.25,e3=0.25,d1=0.05,d2=0.1.    (21)

Using these data, system (1) asymptotically approaches the coexistence equilibrium point E4 = (0.5, 0.2, 0.43), as shown in Figure 1.

FIGURE 1
www.frontiersin.org

Figure 1. Globally asymptotically stable coexistence equilibrium point using data (17) and multiple initial conditions. (A) 3D-Phase portrait of the system (1). (B) Time series for trajectories of x. (C) Time series for trajectories of y. (D) Time series for trajectories of z.

Next, we varied specific parameters to understand the effect of that parameter on the dynamical behavior of the system. For the values of the parameter r satisfying r ≥ 1.84 and other parameters as in equation (17), system (1) approaches the specialist-free equilibrium point, as shown in Figure 2A for the typical value r = 2. For r ≤ 0.76 and other parameters as in equation (17), system (1) asymptotically approaches the superpredator-free equilibrium point, as shown in Figure 2B.

FIGURE 2
www.frontiersin.org

Figure 2. Trajectories of system (1) for different values of r, with other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for r = 2, which approach E2 = (0.8, 0, 0.8). (B) Time series for trajectories of x, y, and z for r = 0.3, which approach E3 = (0.06, 0.22, 0).

For different values of the fear rate k with the rest of the data as given in equation (17), system (1) is solved numerically and illustrated in Figure 3.

FIGURE 3
www.frontiersin.org

Figure 3. Trajectories of system (1) using data given in equation (17) with different values of k. (A) 3D-Phase portrait of system (1) for different values of k. (B) Time series for the trajectory of x, in which x decreases but remains positive as k increases. (C) Time series for the trajectory of y, in which y increases as k increases. (D) Time series for the trajectory of z, in which z decreases as k increases and approaches zero for k > 2.

As shown in Figure 2, the superpredator decreases as k increases, with extinction for k > 2. We also varied the parameter a1 with other parameters still fixed as in equation (17). For the range a1 ≥ 1.3, system (1) asymptotically approaches the superpredator-free equilibrium, as shown in Figure 4A for the typical value a1 = 1.75. Conversely, system (1) approaches asymptotically to the specialist-free equilibrium point, as shown in Figure 4B for the typical value a1 = 0.4.

FIGURE 4
www.frontiersin.org

Figure 4. Trajectories of system (1) for different values of a1 and other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for a1 = 1.75, which approach E3 = (0.03, 0.39, 0). (B) Time series for trajectories of x, y, and z for a1 = 0.4, which approach E2 = (0.8, 0, 0.4).

System (1) still asymptotically approaches the coexistence equilibrium point for the data given in equation (17), when e1 or e3 increases or when e2 decreases. However, for e1 ≤ 0.44, system (1) asymptotically approaches the specialist-free equilibrium point E2 = (0.8, 0, 0.4), as shown in Figure 4B. For e2 ≥ 0.3 or e3 ≤ 0.2 and the other parameters as in equation (17), the trajectories of system (1) asymptotically approach the specialist-free equilibrium point or the superpredator-free equilibrium point as shown in Figure 5.

FIGURE 5
www.frontiersin.org

Figure 5. Trajectories of system (1) for different values of e2 and e3, with other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for e2 = 0.4, which approach E2 = (0.5, 0, 1). (B) Time series for trajectories of x, y, and z for e3 = 0.1, which approach E3 = (0.06, 0.58, 0).

Further investigation for the effect of varying parameters a2 and a3 while keeping the rest of the parameters as in equation (17) shows that the parameter a2 has similar effects as the parameter e2 with a bifurcation point at a2 = 0.63. However, parameter a3 has similar effects as parameter r, with two bifurcation points: a3 = 1.38 and a3 = 0.62.

Finally, for d1 ≥ 0.31 and the rest of the parameters as in equation (17), the trajectories of system (1) asymptotically approach the specialist-free equilibrium point, as shown in Figure 6A for the typical value d1 = 0.4. However, system (1) approaches the coexistence equilibrium point otherwise. On the other hand, for d2 ≥ 0.12 or d2 ≤ 0.08 and the other parameters as in equation (17), trajectories of system (1) asymptotically approach either the superpredator-free equilibrium point or the specialist-free equilibrium point, as shown in Figures 6B,C for the typical values d2 = 0.15 and d2 = 0.08, respectively.

FIGURE 6
www.frontiersin.org

Figure 6. Trajectories of system (1) for different values of d1 and d2, with other parameters as in equation (17). (A) Time series for trajectories of x, y, and z for d1 = 0.4, which approach E2 = (0.8, 0, 0.4). (B) Time series for trajectories of x, y, and z for d2 = 0.15, which approach E3 = (0.06, 0.58, 0). (C) Time series for trajectories of x, y, and z for d2 = 0.08, which approach E2 = (0.64, 0, 0.72).

Varying parameters d1 and d2 together so that they satisfy conditions (7a) and (7b) while keeping other parameters as in equation (17) makes the trajectories of system (1) asymptotically approach the axial equilibrium point E1 = (1, 0, 0) as shown in Figure 7.

FIGURE 7
www.frontiersin.org

Figure 7. Globally asymptotically stable axial equilibrium point for d1 = 0.8, d2 = 0.15 and other parameters as in equation (17). (A) 3D-Phase portrait. (B) Time series for trajectories of x. (C) Time series for trajectories of y. (D) Time series for trajectories of z.

Discussion

Classical intraguild predation food web models have been extensively studied, but relatively little work has been done on the effect of fear on the prey. To the best of our knowledge, this is the first model considering prey fear of only the specialist predator. We determined the equilibria and global stability properties of the model, found bifurcations, and illustrated the theoretical behavior with numerical simulations.

Many previous studies have found that three-species intraguild predation models, in which a superpredator (IG predator) both attacks and competes with a specialist predator (IG prey), are often unstable, either because one consumer is excluded or because long feedback loops produce sustained oscillations [30]. Despite this, many natural IGP systems continue to thrive. Many empirical intraguild predation systems are entrenched in communities with alternative prey species, and standard models of intraguild predation simplify actual systems in significant ways that could affect persistence. Holt and Huxel [30] presented results of theoretical explorations of how alternative prey can influence the persistence and stability of a focal intraguild predation interaction. They reviewed the key conclusions of standard three-species IGP theory and then presented results of theoretical explorations of how alternative prey can influence the persistence and stability of a focal intraguild predation interaction. Conversely, Bai et al. [10] indicated that the intraguild predation food web model has rich and complicated dynamic behavior, and a large Allee effect in the basal prey raises the extinction risk of not only the basal prey but also the IG prey or/and IG predator.

Hossain et al. investigated fear in an intraguild predation model, in which the growth rate of a specialist predator (IG prey) is reduced due to the cost of fear of a superpredator (IG predator), and the growth rate of a basal prey is suppressed due to the cost of fear of both the IG prey and the IG predator [20]. In a three-species food chain system, they found that omnivory can cause chaos in the absence of fear. Fear, on the other hand, can help to keep the chaos under control. They also discovered that the system exhibits bistability between the IG prey-free and IG predator-free equilibrium, as well as bistability between IG prey-free and interior equilibria. Furthermore, they showed that the system can display numerous stable limit cycles for a given set of parameter values.

However, in our study, instilling the fear of a specialist predator (IG prey) in the presence of a superpredator (IG predator) had the unintended consequence of eradicating the superpredator. This demonstrates the often-surprising complexity of food-chain dynamics and has not been observed in previous articles dedicated to the study of fear [2327]. In contrast, if the prey's intrinsic growth rate is sufficiently high, the IG prey can be eradicated, while the IG predator can be removed if the prey's intrinsic growth rate is sufficiently low. Coexistence or the eradication of both predators are other possibilities. As a result, system (1) has rich dynamical behavior, is sensitive to parameter changes, and at least one bifurcation point exists for each of the parameters. This result is comparable to decreasing the productivity of a resource (comparable to a higher fear effect in the prey in this manuscript) resulting in a reduction in the density of the highest trophic level [31].

The most important parameters that influenced the outcome were, in order, the growth rate r, the fear effect k, predation upon the prey by the specialist predator a1, the conversion factor between the prey and the superpredator e2, the conversion factor between the specialist predator and the superpredator e3, the death rate of the specialist predator d1, and the natural death rate of the superpredator d2.

Our model contains several limitations, which should be acknowledged. We only looked at the prey's fear of the specialized predator, not of the superpredator, and we ignored the specialist predator's fear of the superpredator. Our interpretation of such a case is that the pressure of the superpredator on the basal prey is much lower than that of a specialist predator due to the existence of alternative food sources. Mass-action kinetics and a conversion factor are used to model each predator's attack rates, which is a simplification of the genuine dynamics. Finally, the predators were completely reliant on a single food source, which is not always the case.

As a result, the dynamics of a food web in which the prey is afraid of one predator but not the other can produce unexpected results. In the absence of other fear effects, future research will consider the prey's fear of the superpredator but not the specialized predator; the specialist predator's fear of the superpredator could also be incorporated.

Data availability statement

The original contributions presented in the study are included in the article/supplementary files, further inquiries can be directed to the corresponding author.

Author contributions

NF developed the model, wrote the initial draft, conducted analysis, and performed numerical simulations. RN conducted analysis and performed numerical simulations. MH designed the study. SS? rewrote and edited and manuscript. All authors contributed to the article and approved the submitted version.

Funding

SS? is supported by an NSERC Discovery Grant. For citation purposes, please note that the question mark is part of her name.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Thompson RM, Brose U, Dunne JA, Hall RO, Hladyz S, Kitching RL, et al. Food webs: reconciling the structure and function of biodiversity. Trends Ecol Evolut. (2012) 27:689–697. doi: 10.1016/j.tree.2012.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Allesina S, Tang S. Stability criteria for complex ecosystems. Nature. (2012) 483:205–08. doi: 10.1038/nature10832

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Aziz-Alaoui MA. Study of a leslie gower type tri-trophic population model. Chaos Solitons Fractals. (2002) 14:1275–1293. doi: 10.1016/S0960-0779(02)00079-6

CrossRef Full Text | Google Scholar

4. Gakkhar S, Singh B, Naji RK. Dynamical behavior of two predators competing over a single prey. BioSystem. (2007) 90:808–817. doi: 10.1016/j.biosystems.2007.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gakkhar S, Naji RK. On a food web consisting of a specialist and generalist predator. J Biol Syst. (2003) 11:365–376. doi: 10.1142/S0218339003000956

CrossRef Full Text | Google Scholar

6. Gakkhar S, Naji RK. Order and chaos in a food web consisting of a predator and two independent preys. Commun Nonlinear Sci Numer Simulat. (2005) 10:105–120. doi: 10.1016/S1007-5704(03)00120-5

CrossRef Full Text | Google Scholar

7. Klebanoff A, Hastings A. Chaos in one predator two prey model: general results from bifurcation theory. Math. Biosciences. (1994) 122:221–223. doi: 10.1016/0025-5564(94)90059-0

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Holt RD, Polis GA. A theoretical framework for intraguild predation. Am Natur. (1997) 149:745–764. doi: 10.1086/286018

CrossRef Full Text | Google Scholar

9. Tanabe K, Namba T. Omnivory creates chaos in simple food web models. Ecology. (2005) 86:3411–3414. doi: 10.1890/05-0720

CrossRef Full Text | Google Scholar

10. Bai YK, Ruan S, Wang L. Dynamics of an intraguild predation food web model with strong Allee effect in the basal prey. Nonlinear Anal Real World Appl. (2021) 58:103206. doi: 10.1016/j.nonrwa.2020.103206

CrossRef Full Text | Google Scholar

11. Sabelis MW. Predatory arthropods. In: Natural enemies, the population biology of predators, parasites and diseases, eds M. J. Crawley (Oxford: Blackwell Scientific) (1992). p. 225–264. doi: 10.1002/9781444314076.ch10

CrossRef Full Text | Google Scholar

12. Crowder DW, Snyder WE. Eating their way to the top? Mechanisms underlying the 456 success of invasive generalist predators. Biol Invas. (2010) 12:2857–2876. doi: 10.1007/s10530-010-9733-8

CrossRef Full Text | Google Scholar

13. Schreiber SJ. Generalist and specialist predators that mediate permanence in ecological communities. J. Math. Biol. (1997) 36:133–148 doi: 10.1007/s002850050094

CrossRef Full Text | Google Scholar

14. Lima S, Dill LM. Behavioral decisions made under the risk of predation: a review and prospectus. Can. J. Zool. (1990) 68:619–640. doi: 10.1139/z90-092

CrossRef Full Text | Google Scholar

15. Panday P, Pal N, Samanta S, Chattopadhyay J. Stability and bifurcation analysis of a three-species food chain model with fear. Int J Bifurcation Chaos. (2018) 28:1850009. doi: 10.1142/S0218127418500098

CrossRef Full Text | Google Scholar

16. Clinchy LMDi, Roberts D, Zanette LY. Fear of large carnivores causes a trophic cascade. Nat Commun. (2016) 7:10698. doi: 10.1038/ncomms10698

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Panday NP, Samanta S, Chattopadhyay J. A three species food chain model with fear induced trophic cascade. Int J Appl Comput Math. (2019) 5:1–26. doi: 10.1007/s40819-019-0688-x

CrossRef Full Text | Google Scholar

18. Panday P, Samanta S, Pal N, Chattopadhyay J. Delay induced multiple stability switch and chaos in a predator–prey model with fear effect. Math Comput Simulat. (2020) 172:134–58. doi: 10.1016/j.matcom.2019.12.015

CrossRef Full Text | Google Scholar

19. Cong MF, Zou X. Dynamics of a three-species food chain model with fear effect. Commun Nonlinear Sci Numer Simulat. (2021) 99:105809. doi: 10.1016/j.cnsns.2021.105809

CrossRef Full Text | Google Scholar

20. Hossain NP, Samanta S, Chattopadhyay J. Fear induced stabilization in an intraguild predation model international. J Bifurcation Chaos. (2020) 30:2050053. doi: 10.1142/S0218127420500534

CrossRef Full Text | Google Scholar

21. Mukherjee D. Role of fear in predator–prey system with intraspecific competition. Math Comput Simulat. (2020) 177: 263–275. doi: 10.1016/j.matcom.2020.04.025

CrossRef Full Text | Google Scholar

22. Ibrahim HA, Naji RK. Chaos in beddington–deangelis food chain model with fear. J Phys Conf Ser. (2020) 1591:012082. doi: 10.1088/1742-6596/1591/1/012082

CrossRef Full Text | Google Scholar

23. Mondal B, Roy S, Ghosh U, Tiwari PK. A systematic study of autonomous and nonautonomous predator–prey models for the combined effects of fear, refuge, cooperation and harvesting. Eur Phys J Plus. (2022) 137:724. doi: 10.1140/epjp/s13360-022-02915-0

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Roy S, Tiwari PK, Nayak H, Martcheva M. Effects of fear, refuge and hunting cooperation in a seasonally forced eco-epidemic model with selective predation. Eur Phys J Plus. (2022) 137:1–31. doi: 10.1140/epjp/s13360-022-02751-2

CrossRef Full Text | Google Scholar

25. Maity PKTi, Pal S. An ecoepidemic seasonally forced model for the combined effects of fear, additional foods and selective predation. J Biol Syst. (2022) 12:1–37. doi: 10.1142/S0218339022500103

CrossRef Full Text | Google Scholar

26. Pal PKTi, Pal N. Bifurcations, chaos, and multistability in a nonautonomous predator–prey model with fear Chaos: an interdisciplinary. J Nonlinear Sci. (2021) 31:123–34. doi: 10.1063/5.0067046

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Tiwari PK, Verma M, Pal S, Kang Y, Misra AK. A delay nonautonomous predator–prey model for the effects of fear, refuge and hunting cooperation. J Biol Syst. (2021) 29:927–69. doi: 10.1142/S0218339021500236

CrossRef Full Text | Google Scholar

28. Wang X, Zanette L, Zou X. Modelling the fear effect in predator–prey interactions. J Math Biol. (2016) 73:1179–1204. doi: 10.1007/s00285-016-0989-1

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Perko L. Differential Equations and Dynamical Systems, Volume 7. New York: Springer (1996). doi: 10.1007/978-1-4684-0249-0

CrossRef Full Text | Google Scholar

30. Holt RD, Huxel GR. Alternative prey and the dynamics of intraguild predation: theoretical perspectives. Ecology. (2007) 88:2706–2712. doi: 10.1890/06-1525.1

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Gurney W, Nisbet RM. Ecological Dynamics. Oxford: Oxford University Press. (1998).

Google Scholar

Keywords: intraguild predation, food web, fear effect, stability, specialist predator

Citation: Fakhry NH, Naji RK, Smith? SR and Haque M (2022) Prey fear of a specialist predator in a tri-trophic food web can eliminate the superpredator. Front. Appl. Math. Stat. 8:963991. doi: 10.3389/fams.2022.963991

Received: 08 June 2022; Accepted: 05 October 2022;
Published: 03 November 2022.

Edited by:

Aurelio A. de los Reyes V, University of the Philippines Diliman, Philippines

Reviewed by:

Pankaj Tiwari, University of Kalyani, India
Abdelalim A. Elsadany, Suez Canal University, Egypt

Copyright © 2022 Fakhry, Naji, Smith? and Haque. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Mainul Haque, Mainul.Haque@nottingham.edu.cn

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.