Modeling and bifurcation analysis of an intraguild predation system

This paper proposes a mathematical model of an asymmetric intraguild (IG) predation system with an exclusive alternative resource. In particular, this paper analyzes the effects that the exclusive alternative resource has on the consumption/predation behaviors of both the IG predator and IG prey species in the system. The results presented on this paper show that, if the IG predator is less competitive in consumption and less efficient in conversion of the shared resource than that of the IG prey, then there exists a lower bound on the value of the predation rate parameter that should be maintained by IG predator species to ensure its survival and co-existence in the system. This is an open access article under the CC BY-SA license.

INTRODUCTION An intraguild predation (IGP) system refers to simultaneous predation and competition interactions where both IG predator and IG prey species utilize a common resource while also took benefit from preying upon one another [1]. The IGP phenomenon is common in nature and its important implications have been discussed to a large extent from both theoretical and practical perspectives [2]- [4]. IGP systems can exhibit complex and shifting dynamics between different contrasting behaviors [5] which include among others the existence of multiple stable states (e.g. species extinction and coexistence) as well as random shifts among such different dynamics [6], [7]. As the number of species in an IGP system increases, the overall dynamical complexity and observed biodiversity also increase and eventually shape the functioning of an ecosystem [8].
To date, studies on the dynamics of IGP systems have attracted significant interests from researchers across various disciplines. Since its introduction in [9], various modeling and analyses of complex dynamics arising from IGP systems have been investigated in several studies, see for examples [10]- [14]. The work in [16] demonstrates how IGP systems are found to be widspread in various ecosystems and significantly shape the trophic structure/biodiversity and the stability/persistence of an ecosystem. The work in [17] reports that the resource preference of IGP predators play important roles in determining which predators occupy a particular ecosystem. More recently, an investigation of the impact that prey abundance and diversity have on the strength of predator-prey interaction in IGP systems are reported in [18]. Despite these developments, studies about IGP dynamics in the presence of alternative resources remain limited.
This paper proposes a model of an IGP system with exclusive alternative resource and evaluates its qualitative dynamical properties using concepts from stability and bifurcation theories [19], [20]. Specifically, this paper identifies the parameter values that ensure the survival/coexistence of both IG predator prey in the  ❒  ISSN: 2252-8814 presence of alternative resources. The main finding in this paper is that if the IG predator is less efficient than the IG prey in converting the consumed shared resource into its growth utilizations, then the IG predator should satisfies a minimum predation rate parameter value ensure its survival and co-existence in the ecosystem.

AN IGP SYSTEM MODEL
Consider an asymmetric IGP system as shown in Figure 1 which consists of one predator and one prey which consume three resources. Let N 1 and N 2 , respectively, be the IG prey and IG predator biomasses. Assume both species compete for a common resource N 0 which is essential for their growth/reproduction. Each species N 1 and N 2 also has an exclusive alternative resource N 3 and N 4 , respectively, that is essential for survival. All resources satisfy logistic growth functions and their consumptions by the predator/prey species are modeled with type 2 response function [7], [21]. The consumption of the IG prey by IG predator is modeled as type 3 response function, while the death rates of both species are linearly proportional to the number/density of each species. The IGP system model considered in this paper is therefore defined as (1).
In (1), parameters r i , K i , and h i denote, respectively, the growth rate, carrying capacity, and half saturation constant of the ith resource (for i = 0, 3, 4), while a ji and e ji denote the grazing rate and grazing efficiency of species i on resource j (for j = 1, 2), respectively. Furthermore, d j is the death rate of species j, and c k is the attack rate of species j on species k = 1, 2 with k ̸ = j. The time variable is denoted with t. Consider a non-dimensional parameterization [21], [22] of state variables dan parameters of (1) as (2).
(2) Defineẋ i = dx i /dt. Then a parameterization of the state equation model (1) can be written as (3).
All parameters in (3) are chosen to satisfy predation rate efficiency constraints. For instance, |ϵ i | < |α ij | must hold to account for conversion efficiency of species i in consuming resource j, or |α i0 | < 1 to account for shared consumption of resource 0 by species i. Assumptions below are also used in model (3). -Assumption 1: the IG predator is less competitive than the IG prey in consuming the shared resource -Assumption 2: the IG predator is also less efficient than that of the IG prey in converting the consumed shared resource into their energy use or reproduction -Assumption 3: the IG predator and prey have similar consumption rate on their exclusive resource Assumptions 1 and 2 (as defined by setting ϵ 2 < ϵ 1 ) suggest that IG prey may survives if its growth from consuming the shared resource is sufficient enough to overcome the predation pressure from IG predator. The IG predator, on the other hand, should maintains certain level of predation rate on IG prey to balance its loss from being less competitive in the consumption and less efficient in the conversion of the shared resource.

Local stability analysis
Based on the equilibriums data obtained in section 3.1, the system's local stability property at each equilibrium can be analyzed using linearization method. The main objective in this analysis is to examine which equilibrium points will ensure stable species co-existence of both IG predator and IG prey species.
Consider nonlinear system model in (3) which satisfies the following general form.
The elements of the linearized system matrix A in (8) are detailed in [23]. The local stability property of the IGP system at each equilibrium point may then be evaluated using (8), and the results are summarized in Table 1. It can be seen in this table that most of the equilibrium points of model (3)

RESULTS AND DISCUSSION
Since the local stability analysis for a particular parameter set suggests most equilibrium points to be unstable, bifurcation analysis is further conducted to search for possible existence of multiple stable equilibrium points for different combinations of parameter values. The existence of such multiple (nonzero) stable equilibriums essentially indicates the existence several states where all species can survive and co-exist.
The computation of stable co-existence equilibriums of (3) is challenging as it requires the solution of polynomial function equations of order 6 with 18 parameters. To compute possible co-existence equilibriums, bifurcation analysis which starts at a nominal, stable species co-existence equilibrium was conducted. The following initial parameter values are chosen in the bifurcation analysis: k 0 = 20, k 3 = k 4 = 15, α 1,0 = 0.6, α 2,0 = 0.4, α 1,3 = 0.1, α 2,4 = 0.15, c = 0.01, q = 2, δ 1 = δ 2 = 0.01, γ 3 = γ 4 = 1, ϵ 1 = 0.5, ϵ 2 = ϵ 4 = 0.1, and ϵ 3 = 0.01. Figure 2 shows the resulting trajectories and phase portraits, especially in Figures 2(a) and 2(b), Figure 2 show a stable co-existence equilibrium at [x 0 , x 1 , x 2 , x 3 , x 4 ] T = [0.8545, 2.1211, 4.4054, 0.9867, 0.9724] T is reached. This stable point is therefore used as the starting point for the one-parameter bifurcation analysis of model (3). The bifurcation analysis reported in this paper were obtained using MATCONT toolbox [24] in MATLAB [25]. Note that a bifurcation point is essentially a point where the number and qualitative behaviors of an equilibrium change as the parameters are varied. A Hopf (H) bifurcation is a point where an equilibrium changes from stable to unstable, or vice versa. A Limit Point (LP ) bifurcation is a threshold point in which an equilibrium point jumps between two stable states with contrasting qualitative behavior. Finally, a Branching Point (BP ) is the point from which periodic doubling of equilibriums (such as in logistic map) starts to occur.

Bifurcation analysis to identify species coexistence equilibriums
The parameters that are examined in this analysis case are ε i1 (i 1 = 1, . . . , 4), δ i2 (i 2 = 1, 2), and c. Other parameters are kept constant. The results of one-parameter bifucation analysis are plotted in Figure 3. The indicated INIT point in Figure 3 is the initial equilibrium point with response shown in Figure 2. Each red-star mark is a bifurcation point whose corresponding symbol θ i j indicates the i-th bifurcation point of type θ ∈ {H, BP , LP } under variation of one bifurcation parameter j ∈ {ε j , δ j , c}. Each bifurcation point is obtained by computing the forward (increasing) and backward (decreasing) variations of one bifurcation parameter while keeping the other fixed. For example, Figure 3 shows that forward and backward variations of parameter δ i2 produce three bifurcation points which include two Hopf (H 1 δ , H 2 δ ) and one LimitPoint (LP 3 δ ) bifurcation points. Moreover, there exist multiple species co-existence attractors among which the IGP system may jump from one to another when one bifurcation parameter is varied. Table 2 summarizes all bifurcation points and their equilibrium points x * i which occur in the oneparameter bifurcation analysis. Figure 4 illustrates the response (4(a)) and phase portrait (4(b)) of model (3) at the bifurcation parameter ε 1 = 6.0507 which differ significantly from the nominal response shown in Figure 2. The biological interpretation of the results in Figure 3 and Table 2 can be examined based on the corresponding bifurcation parameter. For instance, consider the impacts of variation on IGP prey's consumption rate ε 1 from its nominal value of ε 1 = 0.5 at the indicated INIT point in Figure 3. It can be seen that forward computation of ε 1 drives the IGP system towards bifurcation points H 1 ε2 and H 3 ε1 which correspond to ε 1 values of ε 1 = 1.199 and ε 1 = 6.0507, respectively. However, from biological system modeling perspectives, these values are not reasonable as the consumption rate parameters should always satisfy a constraint of the form |ε 1 | < 1. Thus, one may concludes that consumption rate parameter ε i that is higher than the nominal value will not change the stability of IGP system (3). Similar interpretation may also be inferred for other parameters.

Bifurcation analysis in the presence of consumer's alternative resource
One main question on the IGP model (3) is whether the IG predator can survive by only consuming the shared alternative resource and without preying upon the IG prey. We examine this question for the case when the IG predator is less efficient than the IG prey in converting the consumed shared resource into energy/reproduction uses. Formally, this case can be defined as a constraint on parameters of the form: |ε 2 | < |ε 1 |.
The backward computation on predation rate c was conducted and the bifurcation curve of the IG predator density with respect to parameter c is plotted in Figure 5. This plot shows that, depending on the values of c, two types of equilibriums for the IG predator may occur, namely a stable equlibrium (horizontal curve) and an unstable equilibrium (vertical curve). On one hand, the unstable equilibrium curve shown in Figure 5 tends to result in IG predator extinction at a BP bifurcation point for the value of c = 0.007238. This suggests that c = 0.007238 is a lower bound value of predation rate c which should be maintained by IG predator to avoid extinction. On the other hand, the LP bifurcation which occurs at c = 0.006269 corresponds to the threshold value c which will guarantee the IG predator survival for some time before eventually extincts. Based on the results shown in Figure 5, one may thus defines the following three partitioned values of parameter c which also correspond to different equilibriums/behaviors of the IGP system dynamics: -if 0 ≤ c ≤ 0.006269 (LP ), then the IG predator will extinct without any chance to grow/reproduce. -if (LP ) 0.006269 < c ≤ 0.007238 (BP ), then the IG predator will exist and survive for some period of time before eventually extincts. -if c > 0.007238 (BP ), then the IG predator will survive and co-exist with other species. Its dynamics will exhibit limit cycles below c = 0.0085 and become asymptotically stable above c = 0.0085. From biological perspective, the partitioned values of parameter c suggest that, even when alternative resources exist, if the IG predator is less efficient than the IGP prey in converting the consumed shared resource into energy/reproduction uses, there still exists a lower bound value of c that the IG predator should maintains to ensure its co-existence survival. Moreover, even when the IGP predator consumption rate on alternative resource is c = 1 (100% efficiency), it cannot survive without predating on the IGP prey. One may thus concludes that, to ensure species co-existence in IGP system, the IG predator predation rate on IG prey is a more dominant factor than its consumption rate on the shared alternative resource.
In contrast to the IG predator species, the IG prey consumption rate ε 3 on its alternative resource does not play significant roles to ensure its co-existence in the system. As shown in Figure 6, the IG prey density which starts at the INIT point will increase towards the LP bifurcation point as ε 3 is increased. However, since the constraint 0 < ε 3 < 1 must holds, the LP bifurcation point (ε 3 = 1.109 ) will never be realized. Thus, the IG predator is guaranteed to survive even without consuming the alternative resource. Specifically, the Hopf bifurcation point H 4 where IG prey extincts will never be reached since the system dynamics must follow the bifurcation curve by first visiting points LP, H 2 , H 3 , BP which, in particular, are also not reachable.

CONCLUSION
This paper has developed and analyzed the qualitative properties of an IGP system model in the presence of exclusive alternative resource. Local stability and bifurcation analyses were performed on the model to examine the values of parameter set which guarantee the existence of stable equilibriums where both IG predator and IG prey species can survive and co-exist. One main finding in this paper suggests that, in the presence of alternative exclusive resource, if the IG predator is less efficient than the IGP prey in converting the consumed shared resource into energy/reproduction utilizations, then there is a lower bound on the value of predation rate parameter which the IG predator should maintains to ensure the co-existing survival of both the IG predator and IG prey in the ecosystem. The results presented in this paper also suggest that in order to ensure species co-existence in the ecosystem, the IG predator's predation rate on IG prey is a more dominant parameter/factor than its consumption rate parameter on the alternative resource.