Potential effects of invasive Pterois volitans in coral reefs

The invasion of predatory lionfish (Pterois volitans) represents a major threat to the western Atlantic coral reef ecosystems. The proliferation of venomous, fast reproducing and aggressive P. volitans in coral reefs causes severe declines in the abundance and diversity of reef herbivores. There is also widespread cannibalism amongst P. volitans populations. A mathematical model is proposed to study the effects of predation on the biomass of herbivorous reef fishes by considering two life stages and intraguild predation of P. volitans population with harvesting of adult P. volitans. The system undergoes a supercritical Hopf bifurcation when the invasiveness of P. volitans crosses a certain critical value. It is observed that cannibalism of P. volitans induces stability in the system even with high invasiveness of adult P. volitans. The dynamic instability of the system due to higher invasiveness of P. volitans can be controlled by increasing the rate of harvesting of P. volitans. It is also proven that P. volitans goes extinct when the harvest rate is greater than some critical threshold value. These results indicate that the dynamical behaviour of the model is very sensitive to the harvesting of P. volitans, which in turn is useful in the conservation of reef herbivores.


Introduction
The invasion of Pterois volitans in the western Atlantic has brought a major change to the biodiversity of coral reefs . P. volitans are voracious in nature, spreading rapidly to new marine environments and driving down the populations of reef herbivores drastically (Benkwitt, 2015). The loss of herbivores results in the proliferation of algae, especially the brown algae Lobophora variegata, Dictyota spp. and Sargassum spp., which prevent the growth of corals (Acropora spp. and Montastraea spp.) on the seabed (Bhattacharyya & Pal, 2015). With venomous spines, P. volitans are the top-notch predators in the Atlantic and Caribbean regions. Since top predators like sharks and groupers typically avoid P. volitans and thus fail to keep the species population in check, commercial harvesting of the adult P. volitans species seems to be the only way to mitigate their impact on coral reef ecosystems (Morris & Whitfield, 2009).
Parrotfish (Sparisoma spp.) play an important role in maintaining western Atlantic coral reef ecosystems by consuming algae to control its growth and promote coral recruitment (Fishelson, 1997;Rotjan & Lewis, 2006). P. volitans affect corals by overconsuming Parrotfish and other herbivorous fish that keep algae from overgrowing corals. As observed by Goreau and Hayes (1994), Hare and Whitfield (2003), and Albins and Hixon (2008), in the presence of predatory P. volitans, there is a rapid loss of herbivorous Parrotfish and subsequent loss of corals. Apart from preying on reef herbivores, adult P. volitans exhibit cannibalism, eating juveniles of their own species ). Cannibalism has a strong impact on population dynamics because it reduces the predation pressure on reef herbivores (Rudolf, 2008).
Almost all organisms have a life history that takes them through multiple stages from juvenile to adult (Zhang, Chen, & Neumann, 2000) and cannibalistic interactions are very common in stage-structured populations (Rudolf, 2007). To model the effect of invasive P. volitans, we have considered a two-stage-structured system (Bhattacharyya & Pal, 2014), with Parrotfish and P. volitans following a Holling type III functional response (Luwig, Jones, & Holling, 1978). This response function is sigmoid, rising slowly when resources are rare, accelerating when resources become more abundant and finally reaching a saturated upper limit (Edwards & Brindley, 1999). The rate of loss from the prey population due to predation is defined as the uptake rate of the predator. In this formulation, the per capita uptake rate by the predator is given by where x is the prey abundance. This functional response is parameterized by the constants m and a, where m is the maximal prey uptake rate by the predator, and a is the value of the prey population level when the uptake rate per unit prey is half their maximum value, i.e. f (x) | x=a = m/2.
In our model, we have considered the stage structure of the juvenile and adult P. volitans under the assumption that the adult P. volitans prey on Parrotfish and have reproductive ability. For effective, widespread control of P. volitans, a non-constant harvesting policy, introduced by Lenzini and Rebaza (2010), is used in our model. We will examine the interactions of algae, Parrotfish and P. volitans to determine effective strategies for controlling the growth of P. volitans in coral reefs. We have studied the model analytically as well as numerically, with all proofs relegated to the Appendix 1.

The model
We consider a mathematical model consisting of algae at the first trophic level with concentration P(t) at time t, and Parrotfish at the second trophic level with concentration x(t), feeding on the algae. We also consider a two-stage structure for the top predator P. volitans, with y(t) and z(t) as the concentrations of juvenile and adult P. volitans, respectively. In our proposed model, it is assumed that adult P. volitans prey both on Parrotfish and juvenile P. volitans, whereas juvenile P. volitans do not attack Parrotfish, and have no reproductive ability (Wang & Chen, 1997). Adult P. volitans are harvested with a non-constant harvesting policy that provides diminishing marginal returns of the harvesting organization (Leard, Lewis, & Rebaza, 2008). We make the following assumptions in formulating the mathematical model: (H1) In the absence of Parrotfish, macroalgae have only intraspecic competition, and grow according to the logistic equation with intrinsic growth rate r and carrying capacity K. (H2) The death rate of Parrotfish is proportional to the existing Parrotfish population with a proportionality constant D 1 . (H3) The death rate of juvenile P. volitans and the transformation rate from juvenile to adult P. volitans are proportional to the existing juvenile population with proportionality constants D 2 and μ, respectively. (H4) The death rate of adult P. volitans is proportional to the existing adult population with a proportionality constant D 3 . (H5) m 3 y 2 z a 2 3 + y 2 represents the rate of cannibalism of adult P. volitans (z) by consuming juvenile P. volitans (y) leading to the growth of new juveniles. The growth rate of new juveniles is αm 3 y 2 z a 2 3 + y 2 , 0 < α < 1. Thus, (1 − α)m 3 y 2 z a 2 3 + y 2 represents the reduction in growth rate of juvenile P. volitans due to cannibalism.
The basic equations with all the parameters are where P(0) ≥ 0, x(0) ≥ 0, y(0) ≥ 0 and z(0) ≥ 0. Here 1/μ represents the total time spent by P. volitans in its juvenile stage, h is the maximum harvesting rate of adult P. volitans, c is the concentration of adult P. volitans for which the rate of harvesting is exactly half the maximal harvesting rate, m 1 is the maximal uptake rate of algae by Parrotfish, m 2 is the maximal uptake rate of Parrotfish by the adult P. volitans, m 3 is the maximal uptake rate of juvenile P. volitans by the adult P. volitans and a i are the corresponding half saturation constants (i = 1, 2, 3). The parameters α and α i represent the growth efficiency (0 < α, α i < 1, i = 1, 2) of the organisms; all of these are positive quantities. The parameters K, h, c are environmental variables, while r, μ, m i , a i , D i , α and α i are properties of the organisms.

Non-dimensionalization of the system
Let us change the variables of the system (1) to non-dimensional ones by substitutinḡ and defining non-dimensional parameters After we make the substitutions above and drop the bars for simplicity, the system (1) is reduced to where We assume that all initial values are non-negative. The right-hand sides of the equations in the system (2) are smooth functions of the variables P, x, y, z and the parameters. The following lemma gives the condition for which the solutions of the system (2) are positive. Lemma 3.1: If y(t) and z(t) are always positive, then all possible solutions of the system (2) are positive. Therefore, as long as y(t) > 0 and z(t) > 0 for all t, local existence and uniqueness properties hold in the region = {(P, x, y, z) : P > 0, x > 0, y > 0, z > 0}.

Boundedness and permanence of the system
We first prove that the solutions of Equation (2) with initial values in are bounded, so that Equation (2) represents a biologically meaningful system. The proofs of all the lemmas are given in the Appendix 1.
Lemma 4.1: For all > 0, there exists t > 0 such that all the solutions of (2) with positive initial values fall into the set (P, x, y, z) ∈ : Since P(t) + x(t) + y(t) + z(t) < 1/D as t → ∞ it follows that there exist positive numbers M 1 , M 2 , M 3 with M 1 + M 2 + M 3 < 1/D such that x(t) ≤ M 1 , y(t) ≤ M 2 and z(t) ≤ M 3 for large values of t.
Let us define Then λ denotes the break-even concentration of algae for which the Parrotfish population is constant in the absence of P. volitans. The following lemma states the condition under which neither Parrotfish nor P. volitans can survive in the system: According to Lemma 4.2, we have (i) If the maximal uptake rate of Parrotfish is less than or equal to its death rate, then Parrotfish and P. volitans will not survive in the system. (ii) If the maximal uptake rate of Parrotfish is greater than its death rate and the breakeven concentration λ is greater than 1/D, then Parrotfish and P. volitans will not survive in the system.
The system (2) will be permanent (Ruan, 1993) for each organism u i (t) in the system (i = 1, . . . , 4). Permanence represents convergence on an interior attractor from any positive initial conditions, and hence, can be regarded as a strong form of coexistence. From a biological point of view, the permanence of a system ensures the survival of all the organisms in the long run. Without any loss of generality, we assume that m 1 > D 1 and m 2 > D 2 . The condition given in the following Lemma rules out the possibility of extinction of any organism in the system. Lemma 4.3: If there exists p 1 , 0 < p 1 < λ, then for large t, there exists such that each solution of the system (2) with positive initial values falls into the compact set and stays there. System (2) is called competitive (Smith, 1995) if there exists a diagonal matrix H = diag( 1 , . . . , 4 ) such that HJ(X)H has non-positive off-diagonal elements, where J(X) is the Jacobian of the system (2) and i is either 1 or -1 (i = 1, . . . , 4). By choosing H = diag(1, −1, −1, 1), we see that the off-diagonal elements of HJ(X)H are non-positive if m 3 y 2 (t) a 2 3 + y 2 (t) < m 2 x 2 (t) a 2 2 + x 2 (t) for all t > 0. This leads to the following result: Lemma 4.4: The system (2) is competitive if for large t, where m 2 > m 3 .

Equilibria and their stability
The system (2) possesses the following feasible equilibria: We see that E 0 and E 1 always exist, E 2 exists if λ < 1, and E * exists if P * > λ. We analyse the local stability of system (2) using eigenvalue analysis of the Jacobian matrix, J(X), evaluated at the appropriate equilibrium. The eigenvalues of J(X) at E 0 are 1, −D 1 , −μ − D 2 and −D 3 − h/c. This gives the following result: Lemma 5.1: The organism-free equilibrium E 0 of the system (2) is always a saddle point. Therefore, the system (2) is very unlikely to collapse. Lemma 5.2: The critical point E 1 of the system (2) is locally asymptotically stable if D 1 > m 1 /(a 2 1 + 1). Therefore, with a high mortality rate of Parrotfish, the system (2) stabilizes at the Parrotfish-and P. volitans-free equilibrium E 1 . The decrease in the mortality rate of Parrotfish changes the stability of the system (2) from an algae-dominated state in the absence of Parrotfish to an algae-Parrotfish coexistence state. The following Lemma gives the condition for which coexistence of Parrotfish and algae is possible in the absence of P. volitans.
The system is persistent at E * if all the boundary equilibria repel interior trajectories (Ruan, 1993). The following Lemma gives the condition of persistence of the system (2) at E * : Lemma 5.4: All the organisms will persist if λ < 1 and Therefore, with a low mortality rate of Parrotfish, all the organisms in the system (2) coexist. Lemma 5.5: The system (2) has no periodic solution around the positive equilibrium E * if where L is the minimum of the following six quantities: (a 2 3 +y * 2 ) 2 . Corollary 5.1: If the conditions stated in Lemmas 5.4 and 5.5 both hold, then the positive equilibrium is locally asymptotically stable. Now we use the Routh-Hurwitz criterion to find the necessary and sufficient conditions for stability of the system (2) at E * . Lemma 5.6: The positive equilibrium E * of the system (2) is locally asymptotically stable if Due to the complexity in the algebraic expressions involved, it is difficult to interpret the results in ecological terms; however, numerical simulations are used to illustrate the dynamical behaviour of the system about E * .

Hopf bifurcation
We will study the Hopf bifurcation of the system (2) at E * , taking m 2 as a bifurcation parameter.
The characteristic equation of the Jacobian matrix at β 1 (m 2 ) and β 2 (m 2 ) are real and imaginary parts, respectively, of a pair of eigenvalues for all m 2 ∈ (m 2 cr − , m 2 cr + ).
The condition (ii) is equivalent to dg(m 2 )/dm 2 | (m 2 =m 2cr ) = 0. Thus, using numerical methods, condition (i) can be verified by showing that the curves y = f 1 (m 2 ) and y = f 2 (m 2 ) intersect at m 2 = m 2 cr , whereas the condition (ii) can be verified by showing that the tangent to the curve y = g(m 2 ) at m 2 = m 2 cr is not parallel to the m 2 axis (Siekmann, Malchow, & Venturino, 2008). Corollary 6.1: The period τ of the bifurcating periodic orbits close to m 2 = m 2 cr is given by τ (m 2 cr ) = 2π Q 1 (m 2 cr ) Q 3 (m 2 cr ) .

Stability of bifurcating periodic solution
We investigate the orbital stability of the Hopf-bifurcating periodic solution using Poore's sufficient condition (Poore, 1976). The supercritical/subcritical nature of Hopf-bifurcating periodic solution is determined by the positive/negative sign of the real part of , respectively, where the repeated indices within each term imply a sum from 1 to 4, all the derivatives of F l are evaluated at the equilibrium E * with u 1 = P, u 2 = x, u 3 = y, u 4 = z, and J E * is the Jacobian matrix of (2) calculated at E * . (J E * ) −1 mr denotes the element in row m, column r of (J E * ) −1 .
The left and right normalized eigenvectors of J E * with respect to the eigenvalues ±iω 0 at m 2 = m 2 cr are given by where ξ 1 and ξ 2 are complex numbers, Using a.b = 1, we obtain ξ 1 ξ 2 . If ( ) m 2 =m 2cr > 0, then the system (2) undergoes a supercritical Hopf bifurcation as m 2 is increased through m 2 cr , so that the bifurcating periodic orbit is asymptotically orbitally stable.

Numerical simulations
In this section, we investigate the effects of various parameters on the qualitative behaviour of our system using the numerical approach of Bhattacharyya and Pal (2013) using MATLAB. The default set of parameter values, mostly taken from Bhattacharyya and Pal (2013), is given in Table 1. Under this set of parameter values, it is observed that the system becomes locally asymptotically stable at E * (cf. Figure 1).
We will now verify the feasibility of the criteria of stability in Section 5.   Table 1, the system has a stable focus at E * .
Example 1: If the maximal harvesting rate of adult P. volitans is increased (viz. h = 2.5), and all other parameters are as given in Table 1, then we obtain λ < 1 and 0.0503 satisfying the condition of stability at E 2 = (0.2828, 2.7046, 0, 0) as given in Lemma 5.3. E 2 is a stable focus with eigenvalues −0.0477, −1.9023 and −0.1195 ± i0.2396 (cf. Figure 2).  Table 1, the system has a stable focus at E 2 .  Table 1, the system has a stable node at E 1 .

Example 2:
If the death rate of Parrotfish is increased (viz. D 1 = 0.85), leaving all other parameters unaltered, then the system approaches a stable node at E 1 = (1, 0, 0, 0), with eigenvalues −1, −0.0185, −0.525 and −1.3. In this case, we obtain D 1 > m 1 /(a 2 1 + 1) = 0.194, satisfying Lemma 5.2 (cf. Figure 3).  Table 1 with initial value I 1 . The system is oscillatory around E * (in black). For m 1 = 1.2 and other parameters as given in Table 1, the system is locally asymptotically stable at E 2 (in blue). For m 1 = 0.5 and other parameters as given in Table 1, the system is locally asymptotically stable at E 1 (in red).  Table 1 with initial value I 1 = (0.1, 1, 0.01). The system is oscillatory around E * (in blue). For m 2 = 7 and h = 3.2, other parameters as given in Table 1, the system is locally asymptotically stable at E 2 (in black).
−0.3302 and −1.0076 ± i0.1422. In this case we obtain satisfying the analytical conditions of persistence as given in Lemma 5.4. We also obtain  Table 1, the system is oscillatory around E * (solid). For m 2 = 7, a 3 = 0.5 and other parameter values as given in Table 1, the system is LAS at E * (dotted). satisfying the analytical conditions of stability at E * as given in Lemma 5.6 (cf. Figure 1). Example 4: For m 1 = 2.5, the system is oscillatory around E * (cf. Figure 4). When we decrease the value of m 1 (viz. m 1 = 1.2), the system becomes stable at E 2 . On further lowering the value of m 1 (viz. m 1 = 0.5), we find that the system stabilizes at E 1 .  Table 1.

Example 5:
We observe that the system is oscillatory around E * for m 2 = 7 (cf. Figure 5). In this case by increasing the maximal harvesting rate of adult P. volitans (viz. h = 3.2), the system becomes locally asymptotically stable at E 2 . Example 6: For high invasiveness of adult P. volitans, the system becomes oscillatory around E * (cf. Figure 5). When the value of a 3 is lowered (viz. a 3 = 0.5), the system stabilizes at E * (cf. Figure 6).

Hopf bifurcation
We observe that the system becomes oscillatory when the maximal uptake rate of adult P. volitans is high. We therefore consider m 2 as a bifurcation parameter. For m 2 < 5.2 we see that f 1 (m 2 ) > f 2 (m 2 ), satisfying the Routh-Hurwitz condition, so that the system is locally asymptotically stable at E * . For m 2 > 5.2 we see that f 1 (m 2 ) < f 2 (m 2 ) and so the system is unstable at E * (cf. Figure 7(a)). Moreover, we observe that the tangent to the graph of g(m 2 ) at m 2 = 5.2 is not parallel to the m 2 axis (cf. Figure 7(b)), satisfying the condition dg/dm 2 | (m 2 =5.2) = 0. In Figure 8, we observe that a supercritical Hopf bifurcation occurs when the parameter m 2 is increased through the critical value m 2 cr = 5.2.

Discussion
We have considered a tri-trophic food chain model consisting of algae and Parrotfish in the first two trophic levels, respectively, while juvenile and adult P. volitans reside in the third trophic level. We analyse the effect of predation with stage-structured cannibalism and study the effect of harvesting of the adult P. volitans on the dynamics of the system. The threshold values for the existence and stability of various steady states of the system are worked out. In order to keep sustainable development of the coral reef ecosystem, it is desirable to have a positive equilibrium which is asymptotically stable. Keeping this view in mind, we have established some strong criteria for existence of the positive equilibrium. We studied bifurcation with respect to the parameter representing the invasiveness of adult P. volitans in the system. The critical parameter value at which bifurcation occurs is determined to preserve the system under consideration in its natural state. We observe that when the maximal uptake rate of adult P. volitans on Parrotfish crosses a certain critical value, the system enters into Hopf bifurcation that induces oscillation around the positive equilibrium. The stability as well as the direction of Hopf bifurcation near the interior equilibrium is obtained by applying the algorithm due to Poore (1976). We have also provided numerical simulations to substantiate our analytic results. From analytical and numerical observations we obtain the following conclusions: (i) If the growth rate of Parrotfish is low, then the Parrotfish would become extinct.
(ii) Higher mortality rate of Parrotfish can lead to the extinction of both Parrotfish and P. volitans from the system. This represents the fact that rapid elimination of herbivorous fish can be fatal for the coral reef ecosystem. (iii) High rate of predation of adult P. volitans on Parrotfish induces oscillation around the positive equilibrium leading to dynamic instability, representing the phenomenon of ecological imbalance due to high invasiveness of P. volitans in the coral reef ecosystem. This dynamic instability can be controlled by increasing the rate of harvesting of adult P. volitans. Moreover, a high harvesting rate of adult P. volitans can eliminate P. volitans from the system. (iv) The increase of cannibalism of P. volitans stabilizes the system even with high invasiveness of adult P. volitans.
Throughout the article we focus on searching for a suitable way to control the growth of algae, Parrotfish and P. volitans, and maintain a stable coexistence of all the species. Our numerical simulations suggest that the maximal harvesting rate of adult P. volitans can be used as a control parameter to maintain the stability of the system at the coexistence steady state.
Our results are based on a model that has no growth equation for the corals. It would be interesting to incorporate corals in our system to study the dynamics of coral reefs in the presence of invasive P. volitans.

Disclosure statement
No potential conflict of interest was reported by the authors.

A.5. Proof of Lemma 4.4
Proof: The Jacobian of the system (2) at (P, x, y, z) is given by Let us consider Then we have All the off-diagonal elements of the matrix HJH are negative if m 3 y 2 /(a 2 3 + y 2 ) < m 2 x 2 /(a 2 2 + x 2 ). Therefore, if m 2 > m 3 and all the off-diagonal elements of HJH are negative, and consequently, the system (2) is competitive.

A.6. Proof of Lemma 5.1
Proof: The Jacobian matrix at E 0 is At E 0 , the eigenvalues of the Jacobian matrix are 1, −D 1 , −μ − D 2 and −D 3 − h/c. Therefore, the system (2) is always unstable at E 0 .

A.9. Proof of Lemma 5.4
Proof: In order to prove the persistence of the system, we shall show that all the boundary equilibria of the system are repellers. It is observed that the system is always unstable at E 0 . If D 1 ≤ m 1 /(a 2 1 +1), then the system is unstable at E 1 . The system is unstable at E 2 if D 1 ≤ λ(1 − λ) a 2 cμm 2 (μ + D 2 )(cD 3 + h) − 1 or D 1 ≥ m 1 2(1 − λ) .
Since m 1 /(a 2 1 + 1) < m 1 /{2(1 − λ)}, it follows that all the boundary equilibria are repellers if We have also proved that the system is bounded. Therefore, the system is persistent under the aforesaid conditions.