Dynamics of asymmetric intraguild predation with time lags in reproduction and maturation

A three dimensional (3D) stage-structured predator–prey model is proposed and analyzed to study the effect of intraguild predation with harvesting of the adult species. Time lags in reproduction and maturation of the organism are introduced in the system and conditions for local asymptotic stability of steady states of delay differential forms of the ODE model are derived. The length of the delay preserving the stability is also estimated. Moreover, it is shown that the system undergoes a Hopf bifurcation when the time lags cross certain critical values. The stability and direction of the Hopf bifurcations are determined by applying the normal form method and the center manifold theory. Computer simulations have been carried out to illustrate various analytical results. Subjects: Mathematics & Statistics; Applied Mathematics; Dynamical Systems; Mathematical Biology; Mathematical Modeling; Science


Introduction
Many organisms go through multiple life stages as they proceed from birth to death. A stage-structured model of population growth consists of juvenile and adult organisms where the juveniles have no reproductive ability (Wang & Chen, 1997). Field observations by Rudolf (2007) suggest that cannibalism is a prevalent feature of stage-structured populations in both aquatic and terrestrial food webs. In this type of ecological interaction, two species of the same trophic level interact as predator and prey, known as intraguild predation (IGP). As stated by Polis, Myers, and Holt (1989), IGP can be classified as

ABOUT THE AUTHORS
Our research group is headed by Samares Pal and includes Joydeb Bhattacharyya as a post-doctoral researcher. Research in the group addresses a wide range of questions broadly concerning the dynamics of coral-reef ecosystems under natural and anthropogenic stress. The main emphasis of our work includes mathematical modeling of coral reef ecosystems affected by seaweed allelopathyinduced coral diseases, overfishing of herbivores and invasion of predatory Pterois Volitans. The research reported in this paper will help in studying the catastrophic regime shifts in coralreefs under the proliferation of predatory Pterois Volitans and the subsequent loss in herbivory.

PUBLIC INTEREST STATEMENT
Organisms go through multiple life stages as they proceed from birth to death. Cannibalism has become a prevalent feature of stage-structured populations in both aquatic and terrestrial food webs. Two fishes, Parrotfish and Pterois Volitans in coral reef ecosystem, are considered for the study of interactions within and among the stagestructured populations. With the proliferation of Pterois Volitans, there is a reduction in the density of herbivorous reef-fish. The reduction in herbivory in coral reefs can lead to a change of community structure of coral reefs for which corals decline with the proliferation of seaweeds. In this paper it is observed that cannibalism of Pterois Volitans population ensures the release of predation pressure on Parrotfish population. With the proliferation of Pterois Volitans, the system becomes oscillatory. Harvesting of adult Pterois Volitans helps in stabilizing the system. symmetrical or asymmetrical. Symmetric IGP occurs when there is mutual predation between two species, whereas in asymmetrical interaction, one species consistently preys upon the other (Zhang, Chen & Neumann, 2000). Stage structure with cannibalism is an asymmetric IGP (Rudolf, 2008).
Two fishes, Parrotfish and Pterois Volitans are considered for the study of interactions within and among the stage-structured populations. As observed by , Pterois Volitans are one of the top levels of the food web in many coral reef environments. They are known to feed mostly on small fishes, which include juveniles of their own species. Apart from preying on Parrotfish, adult Pterois Volitans exhibits a distinct cannibalistic attitude towards its juvenile species. As observed by Goreau and Hayes (1994), Hare and Whitfield (2003), and Albins and Hixon (2008), the proliferation of predatory Pterois Volitans reduces the population density of herbivorous Parrotfish by changing the community structure of coral reefs for which corals decline with an increase in abundance of seaweeds. According to NOAA (National Oceanic and Atmospheric Administration), commercial harvesting of adult Pterois Volitans is required to reduce the numbers of Pterois Volitans to mitigate their impact on coral reef ecosystems (Morris & Whitfield, 2009).
The role of time delay on ecosystem models has been investigated by Jiao, Yang, Cai and Chen (2009), and Dhar and Jatav (2013). Motivated by their works, we have studied a stage-structured system (Bhattacharyya & Pal, 2013) with Parrotfish at the first trophic level, and juvenile and adult

The basic model
We take a stage-structured model where Parrotfish is growing on the system with concentration x(t) at time t. The organisms, juvenile and adult Pterois Volitans, are introduced in the system with concentrations y J (t) and y A (t), respectively, at time t. It is assumed that adult Pterois Volitans is cannibalistic and the juveniles have no reproductive ability-an exhibition of asymmetric IGP. The parameter 1 represents the delay in reproduction and 2 represents the delay in maturation of Pterois Volitans. Adult Pterois Volitans is harvested to control its growth.
We make the following assumptions in formulating the mathematical model: (H1) In absence of Pterois Volitans, the growth equation of Parrotfish follows logistic growth with intrinsic growth rate r and carrying capacity K.
(H2) The per capita feeding rate of juvenile Pterois Volitans is assumed to follow Holling II functional response. Adult Pterois Volitans feed on Parrotfish and juvenile Pterois Volitans; accordingly the per capita feeding rate of adult Pterois Volitans follows a generalized form of Holling II functional response as a function of the concentrations of Parrotfish and juvenile Pterois Volitans.
(H3) The mortality rates of juvenile and adult Pterois Volitans are d 1 and d 2 , respectively. The rate of maturation of juvenile Pterois Volitans is proportional to the concentration of the juvenile population with proportionality constant .
(H4) Adult Pterois Volitans (y A ) prey on Parrotfish (x) with uptake rate m 1 xy A a 1 +x+b 1 y J which leads to the delayed growth of new juveniles (y J ) with growth rate e 1 m 1 x(t− 1 )y A (t− 1 ) a 1 +x(t− 1 )+b 1 y J (t− 1 ) .
(H5) Also, adult Pterois Volitans (y A ) exhibit cannibalism on juveniles (y J ) with m 2 y J y A a 2 +x+b 2 y J as uptake rate leading to the delayed growth of new juveniles (y J ) with growth rate .
Thus, m 2 represents the reduction in growth rate of juvenile The basic equations with all of the parameters are: with the initial conditions Here 1 represents the total time spent by Pterois Volitans in its juvenile stage and h is the harvesting rate of adult Pterois Volitans. Also m, m i are the maximum growth rates, a, a i are the half saturation constants, which are the concentrations of resource at which the growth rate of the organism is half maximal, x + b i y J are the interference of Parrotfish and juvenile Pterois Volitans on the per capita growth rate of adult Pterois Volitans, and e i are growth efficiencies (i = 1, 2); all of these are positive quantities.

Invariance and boundedness of the System
Obviously, the right-hand sides of the equations in system (1) are continuous nonnegative smooth functions on R 3 + = (x, y J , y A ):x, y J , y A ≥ 0 . Indeed they are Lipschitzian on R 3 + and so the solution of the system (1) exists and is unique. Therefore, it is possible of Theorem 3.1. to prove that the interior of the positive octant of R 3 is an invariant region.

Qualitative analysis of the system without delay
In this section, we restrict ourselves to analyze the model in the absence of delay. Thus, only the interaction parts of the model system are taken into account. In the absence of delay, the system (1) reduces to where x(0) ≥ 0, y J (0) ≥ 0, and y A (0) ≥ 0.
The system will be permanent if there exist (Hofbauer & Sigmund, 1989). Permanence represents convergence on an interior attractor from any positive initial conditions and so it can be regarded as a strong form of coexistence (Ruan, 1993). From a biological point of view, permanence of a system ensures the survival of all the organisms in the long run.
The following lemma rules out the possibility of extinction of any organism in the system under suitable conditions.
, then for large values of t each solution of (2) enters in the compact set ≤ M 2 and remains in it finally.
We analyze the stability of system (2) using eigenvalue analysis of the Jacobian matrix evaluated at the appropriate equilibrium. The detailed calculations are given in appendix.
Lemma 4.2 The critical point E 0 of the system (2) is always unstable. (3) Thus, under no circumstances the system (2) collapses.
The characteristic equation of the system (2) at E 1 = (K, 0, 0) is = 0, without one eigenvalue −r and the remaining two eigenvalues will have negative real parts if ( Thus, with high harvesting rate of adult Pterois Volitans, the system stabilizes at E 1 . The characteristic equation of the system (2) at E * is 3 Lemma 4.4 The critical point E * of the system (2) is locally asymptotically stable if h > h * and > 0, Thus, harvesting of adult pterois volitans plays a critical role on the stability of the system (2) at the interior equilibrium.
Lemma 4.5 The system (2) undergoes a Hopf bifurcation at h = h cr iff (h); 1 (h) and 2 (h) are real and imaginary parts, Respectively, of a pair of eigenvalues in h ∈ (h cr − , h cr + ).
The condition (iii) is equivalent to Thus, using numerical methods, condition (ii) can be verified by showing that the curves y = f 1 (h) and y = f 2 (h) intersect at h = h cr , whereas the condition (iii) can be verified by showing that the tangent to the curve y = g(h) at h = h cr is not parallel to the h axis (Siekmann, Malchow & Venturino, 2008).

Qualitative analysis of the system with delays
In this section, we analyze the stability of the delayed systems at different equilibria and existence of local Hopf bifurcations.
Since r > 0, it follows that the critical point E 0 of the systen (1) is always unstable.
The characteristic equation of the system (1) at E 1 is The Jacobian of the system (1) at E * is In order to investigate the distribution of roots of the characteristic equations, we use the following lemma by Ruan and Wei (2003).
Lemma 5.1 For the transcendental equation

System with reproduction time delay only
( 1 > 0, 2 = 0) We now consider the case in which reproduction of the adult species is not instantaneous, but mediated by some discrete time lag 1 .
In the absence of maturation time delay, the system (1) reduces to Lemma 5.1.1 If 1 is small, stability of the system (2) at E 1 implies the stability of system (3) at E 1 . Proof Let the system (2) is locally asymptotically stable at E 1 . Then h > e 1 m 1 K (a 1 +K)(d 1 + ) − d 2 holds. If 1 is small, we can write e − 1 = 1 − 1 .
In this case, the characteristic equation of the system (3) at E 1 becomes , all the roots of this equation have negative real parts, and so the stability of the system without delay at E 1 implies the stability of system (3) at E 1 .
Lemma 5.1.2 A sufficient condition for the system (3) to be locally asymptotically stable at E 1 is Proof The system (3) is locally asymptotically stable at E 1 for all 1 ≥ 0 if the following conditions given by Gopalsamy (1992) and Beretta and Kuang (2002) hold: (a) the real parts of all the roots of D E 1 ( , 0, 0) = 0 are negative, (b) for all real and any Sufficient conditions for the nonexistence of ∈ R satisfying D E 1 (i , 1 , 0) = 0 can be written as: Thus, for all real and for any Therefore, the statement of lemma (5.1.2) holds.
The characteristic equation of the system (3) at E * is Lemma 5.1.3 If 1 is small, stability of the system (2) at E * implies the stability of system (3) at E * . Proof Let the system (2) is locally asymptotically stable at E * . Then we have A > 0, C > 0 and AB > C.
For small 1 , we can write e − 1 = 1 − 1 and so the characteristic equation of (3) at E * becomes and so for Therefore, for small 1 > 0, the stability of the system without delay at E * implies the stability of system (3) at E * .
Proof Let u(t), v J (t), andv A (t) be the respective linearized variables of the model.
Taking the Laplace transformation of system (4), we have The inverse Laplace transform of ū(s),v J (s), and v A (s) will have terms which exponentially increase with time if ū(s),v J (s), and v A (s) have a pole with positive real parts. Since E * needs to be locally asymptotically stable, it is necessary and sufficient that all poles of ū(s),v J (s), and v A (s) have negative real parts. We shall employ the Nyquist criterion, which states that if s is the arc length of a curve encircling the right half-plane, the curve , and will encircle the origin a number of times equal to the difference between the number of poles and the number of zeroes of , and in the right half-plane.

Let
. Also, let be the smallest positive root of Then is locally asymptotically stable if and .
Also, Therefore, the positive solution of is always greater than or equal to .

We obtain
where is independent of .
where and Therefore, is an estimate for the length of delay for which the stability of the system at is preserved.
We The smallest is given by and we take it as . Let .
Thus, sign d(Re ) , =i k ≠ 0 and so the system (3) undergoes a Hopf bifurcation at

System with only maturation time delay
Now we consider the case in which the maturity time from juvenile to adult is mediated by a discrete time lag 2 , while the reproduction process is instantaneous.
In the absence of reproduction time delay, the system (1) reduces to The characteristic equation of the system (5) at E 0 is In this case, the characteristic equation of the system (5) at E 0 becomes Since 0 < 2 < 1 and r > 0, it follows that at least one eigenvalue of the characteristic equation will always have positive real part and consequently, E 0 is always a saddle point of the system (5).
The characteristic equation of the system (5) at E 1 is Lemma 5.2.2 If 2 is small and 0 < 2 < 1 , the stability of the system (2) at E 1 implies the stability of system (5) at E 1 . Proof Let the system (2) be locally asymptotically stable at E 1 . Then we have h > If 1 is small, we can write e − 2 = 1 − 2 . Then, the characteristic equation of the system (5) at E 1 > 0 and so, the stability of the system (2) at E 1 implies the stability of system (5) at E 1 .
Lemma 5.2.3 A sufficient condition for the system (5) to be locally asymptotically stable at E 1 is Sufficient conditions for the nonexistence of ∈ R satisfying D E 1 (i , 0, 2 ) = 0 can be written as: Thus, for all real and for any 2 ≥ 0, Therefore, the statement of lemma (5.2.3) holds.
For small 2 , we can write e − 2 = 1 − 2 and so the characteristic equation of (5) at E * becomes Thus, for small 2 > 0, the stability of the system (2) at E * implies the stability of system (5) at E * .
Taking the Laplace transformation of system (6), we have We shall employ the Nyquist criterion, which states that if s is the arc length of a curve encircling the right half-plane, the curves ū(s),v J (s), and v A (s) will encircle the origin a number of times equal to the difference between the number of poles and the number of zeroes of ū(s),v J (s), and v A (s) in the right half-plane.

System with both delays ( 1 , 2 > 0)
The characteristic equation of the system (1) at E 0 is independent of 1 and is identical to the characteristic equation of the system (5) at E 0 . Therefore, for all 1 > 0 and for small 2 satisfying 0 < 2 < 1 , the critical point E 0 is always a saddle point of the system (1).
Lemma 5.3.1 For small 1 , 2 satisfying 0 < 2 < 1 , the stability of the system (2) at E 1 implies the stability of system (1) at E 1 . Proof Let the system (2) be locally asymptotically stable at E 1 .
The characteristic equation of the system (1) at E 1 becomes ( 1 + 2 ) e 1 m 1 K a 1 +K > 0 and so, the stability of the system without delay at E 1 implies the stability of system (1) at E 1 .
Thus, for small 1 , 2 , the stability of the system (2) at E * implies the stability of system (1) at E * .

Lemma 5.3.3 If the system (2) is stable at E * and if there exists
Proof Let u(t), v J (t), andv A (t) be the respective linearized variables of the model. Then system (1) can be expressed as Taking the Laplace transformation of system (7), we have We shall employ the Nyquist criterion, which states that if s is the arc length of a curve encircling the right half-plane, the curves ū(s),v J (s), and v A (s) will encircle the origin a number of times equal to the difference between the number of poles and the number of zeroes of ū(s),v J (s), and v A (s) in the right half-plane.
Therefore, the positive solution v * of Ā * v 2 − vB * −C * = 0 is always greater than or equal to v * , are estimates for the length of delays for which the stability of the system at E * is preserved (i = 1, 2). Now, we consider D E * ( , 1 , 2 ) = 0 in its stable interval and regard 1 as a parameter.
Without any loss of generality, we assume that Lemma 4.4 and 5.2.6(ii) hold.

Direction and stability of the Hopf bifurcation
In this section, we shall study the bifurcation properties using techniques from normal form and center manifold theory by Hassard, Kazarinoff, and Wan (1981). Throughout this section, we assume that system (1) undergoes Hopf bifurcation at the positive equilibrium E * for 1 = * 1 , and ±i * are the corresponding roots of D( , 1 , 2 ) = 0 at E * .
The sign of 2 indicates the direction of bifurcation, while that 2 determines the stability of z(t, 1 ( )). z(t, 1 ( )) is stable if 2 < 0 and unstable if 2 > 0. To derive the coefficients in these expansions, we first construct the coordinates to describe a center manifold Ω 0 near 1 = 0, which is a local invariant, attracting a two dimensional (2D) manifold.
On the manifold Ω 0 we have W(t, In fact, z and z are local coordinates of center manifold Ω 0 in the direction of q and q * , respectively.

Numerical simulations
In this section, we investigate numerically the effect of the various parameters on the qualitative behavior of the system using parameter values given in Table 1 throughout, unless otherwise stated.
We carried out numerical simulations and interpret the resulting figures by varying the parameter(s) under investigation, keeping all other parameters fixed. (12)

The effect of varying the carrying capacity (K)
With low carrying capacity of Parrotfish (viz., K = 0.2), the juvenile and adult Pterois Volitans cannot survive and all the systems become stable at E 1 (cf. Figure 2). By increasing the value of K (viz., K = 2.7), the system without delay and the system with maturation time delay only are stable at E * , whereas, the systems with reproduction delay only and with both delays are oscillatory around E * (cf. Figure 3 (solid)). With high carrying capacity of Parrotfish (viz., K = 3), all the systems are oscillatory around E * (cf. Figure 4 (solid)).

The effect of varying the intrinsic growth rate (r)
With low intrinsic growth rate of Parrotfish (viz., r = 1.5), all the systems are oscillatory around E * (cf. Figure 5 (solid)).
By increasing the value of r (viz., r = 3.5), the system without delay and the system with maturation time delay only are stable at E * , whereas, the systems with reproduction delay only and with both delays are oscillatory around E * (cf. Figure 6 (solid)).  Table  1, all the systems are LAS at E * (dotted).

The effects of invasion
With high invasiveness of adult Pterois Volitans, the maximal growth rate of adult Pterois Volitans on Parrotfish becomes high. It is observed that, for m 1 = 3.2 the system without delay and the system with maturation time delay only are stable at E * , whereas, the systems with reproduction delay only and with both delays are oscillatory around E * (cf. Figure 7 (solid)).

The effects of cannibalism
The rate of cannibalism is dependent on the maximal growth rate of adult maximal growth rate of on its juveniles. In absence of cannibalism (viz., m 2 = 0), all the systems become oscillatory around E * (cf. Figure 9 (solid)).  Table  1, all the systems are LAS at E * (dotted).

The effects of harvesting
In absence of harvesting of adult Pterois Volitans the systems are oscillatory around E * (cf. Figure 11 (solid)). By increasing the harvesting rate (viz., h = 0.15), the system without delay becomes stable at E * and the systems with delay are oscillatory around E * (cf. Figure 11 (dotted)). Further increase of harvesting (viz. h = 0.36) stabilizes all the systems at E * . With high rate of harvesting, juvenile as well as adult Pterois Volitans are eliminated from the system. Further, we observe the following effects:   Figure 4 (dotted)).
(b) With low intrinsic growth rate of Parrotfish (viz., r = 1.5), all the systems become oscillatory around E * . By increasing the rate of harvesting of adult Pterois Volitans (viz., h = 0.5), the systems stabilize at E * (cf. Figure 5 (dotted)).

Hopf bifurcation
We observe that in absence of harvesting of adult Pterois Volitans, other parameter values as in Table 1, the system (2) is oscillatory around E * (Figure 11 (solid)). Increasing the maximal rate of harvesting, the system (2) becomes locally asymptotically stable at E * (Figure 11 (dotted)). We therefore consider h as a bifurcation parameter. From Figure 12(a) it is observed that f 1 (h) and f 2 (h) intersect at h = 0.268, indicating that the system (2) changes its stability when the parameter h crosses the threshold h cr = 0.268. More specifically, for h > h cr we see that f 1 (h) > f 2 (h), satisfying the condition of stability at E * as as given in Lemma 4.4. Also, for h < h cr we see that f 1 (h) < f 2 (h) and so the system (2) is unstable at E * . Moreover, we observe that the tangent to g(h) at h cr = 0.268 is not parallel to the h axis ( Figure 12(b)), satisfying the condition dg dh | (h=0.268) ≠ 0 (Lemma 4.5(iii)). Thus the system (2) undergoes a subcritical Hopf bifurcation when the parameter h is increased through h cr = 0.268.

Discussion
We have considered an intraguild system in which adult Pterois Volitans exhibit cannibalism toward juveniles of its species and are subjected to harvesting. The main objective of this paper is to study the dynamic behavior of the system in the presence of discrete time lags in reproduction and    maturation of Pterois Volitans. Experimental observations by Murray et al. (2013) reveal that Holling II functional response is quite accurate in predicting the observed functional response of fishes. This prompts us to use Holling II response function in our model. Based upon the observations of Castillo-Chavez et al. (2002) it is reasonable to assume that the mortality and maturity rates of fishes are proportional to the number of fishes present in our system. We have shown that solutions of the system are bounded in the long run. We have obtained conditions for permanence and the stability of the system at the coexistence equilibrium. It is observed that if the maximal harvesting rate of adult Pterois Volitans exceeds a certain threshold, the system stabilizes at the coexistence equilibrium through a subcritical Hopf bifurcation. Also, we have established that the systems remain locally asymptotically stable and all solutions approach E * whenever the magnitude of the delays lies below some threshold values. In order to find the expression for these threshold values, we have obtained the estimated length of stability preserving delays. The stability as well as the direction of bifurcation is obtained by applying the algorithm due to Hassard, Kazarinoff, and Wan (1981) that depends on the centre manifold theorem. From numerical simulation it is observed that in the absence of cannibalism and with low harvesting rate of adult Pterois Volitans, the IG-systems become oscillatory around the interior equilibrium. In this case, increase of cannibalism up to a certain threshold stabilizes the system at the interior equilibrium, justifying the observations of Bosch, Roos, and Gabriel (1988) that decrease of predator density due to cannibalism releases prey from predation pressure by inducing stability. Also, with high rate of cannibalism of adult Pterois Volitans, the coexisting populations show oscillatory dynamics. This supports the observations from previous modeling analyses by Diekmann et al. (1986), Hastings (1987, and Magnússon (1999) that high cannibalism level can have destabilizing effects leading to oscillations. The IG-systems become oscillatory around the interior equilibrium with low intrinsic growth rate of Parrotfish in absence of harvesting of adult Pterois volitans, a phenomenon which has not been observed in the previous studies by Dhar and Jatav (2013), Bhattacharyya and Pal (2013), and deserve further investigation. If the reproduction time delay is high, the system with both time lags becomes oscillatory around the positive equilibrium. Also, high rate of invasion of adult Pterois Volitans on Parrotfish induces oscillations around the positive equilibrium, leading to dynamic instability. This represents the phenomenon of ecological imbalance due to the presence of the invasive Pterois Volitans in a coral reef ecosystem, justifying the observations of Albins and Hixon (2008) that high predation rates of adult Pterois Volitans are detrimental to coral reef ecosystems. The dynamic instability can be controlled by increasing the maximal harvesting rate of adult Pterois Volitans. This too justifies that harvesting of these species is beneficial for coral reef ecosystem as observed by Morris, Shertzer and Rice (2011).
Throughout the article an attempt is made to search for a suitable way to control the growth of Parrotfish and Pterois Volitans and maintain stable coexistence of all the species. From analytical and numerical observations, it is seen that increase of the harvesting of adult Pterois Volitans induces stability of the system. Moreover, we observe that harvesting at higher rates are necessary to obtain stability of the systems with delay than that of the system without delay. Also, x ≤ K as t → ∞. Therefore, there exists T 2 > 0 such that x(t) ≤ K for all t ≥ T 2 .
Again, for all t ≥ max{T 1 , T 2 }, Then dy J dt > 0 for y A (t) ≥ y A 1 > 0 and t > max{T 1 , T 2 } and so in this case there exists T 3 > 0 and 0 < y J 1 < M 1 such that y J (t) ≥ y J 1 for all t ≥ T 3 .
Therefore, for all t ≥ T 3 , if y A 1 ≤ y A (t) ≤ M 2 , then y J 1 ≤ y J (t) ≤ M 1 .
Therefore, by the Routh-Hurwitz criterion of stability, the system (2) is stable at E * if > 0 and h > h * .

Proof of Lemma (4.5)
The characteristic equation of the variational matrix at E * is 3 + A 2 + B + C = 0, where The necessary and sufficient conditions for Hopf bifurcation to occur at h = h cr are that (i)h cr > h * and C(h cr ) > 0, (ii)f 1 (h cr ) = f 2 (h cr ),