BIFURCATION AND STABILITY ANALYSIS FOR A NUTRIENT-PHYTOPLANKTON MODEL WITH TOXIC EFFECTS

. In this paper, we analyze a nutrient-phytoplankton model with toxic eﬀects governed by a Holling-type III functional. We show the model can undergo two saddle-node bifurcations and a Hopf bifurcation. This results in very interesting dynamics: the model can have at most three positive equilibria and can exhibit relaxation oscillations. Our results provide some insights on understanding the occurrence and control of phytoplankton blooms.


Introduction.
Recently, there has been a great concern on the outbreak of phytoplankton blooms, which are a sudden increase in the population of phytoplankton in water systems. Many factors such as nutrient enrichment and the increasing temperatures due to global warming contribute to triggering the phytoplankton blooms. However, the underlying mechanism is still not well understood. There has been evidence that many phytoplankton species produce toxin [4,7]. This, according to the work of [9], might benefit the ecosystems with diversity and coexistence of competing phytoplankton species. It is also of great importance to understand if toxin producing phytoplankton (TPP) plays a role in causing the outbreak of phytoplankton blooms. In this regard, studies in [1,3,8] have shown that toxic chemicals produced by TPP may contribute to the outbreak of phytoplankton blooms.
The model proposed by Chakraborty et al. [1] reads as follows where u and v represent the nutrient concentration and the phytoplankton abundance, respectively. The constant b > 0 is the rate of constant external nutrient input into the system, µ is the per capita-loss rate of nutrient, a is the maximal nutrient uptake rate of phytoplankton, ε is the conversion factor, θ is the rate of release of toxic chemicals by the TPP population, δ is the per capita-mortality rate of phytoplankton, and ξ is the half-saturation constant.
Note that for Model (1), only simple analysis was carried out in [1]. In this paper, we will give detailed stability and bifurcation analysis. More precisely, we will show that the model undergoes saddle-node bifurcations and a Hopf bifurcation. Moreover, we (numerically) show that the model exhibits relaxation oscillations. Our analysis and numerical result will provide some insights on understanding the occurrence of phytoplankton blooms.
We organize the rest of the paper as follows. We first discuss the existence of equilibria and locate the saddle-node bifurcation curve numerically in Section 2, then we carry out stability analysis in Section 3. In Section 4, we provide some numerical simulations and we briefly summarize our results in the last section.
2. Existence of equilibria. It has been shown in [1] that Model (1) is well-posed in the sense that for any given nonnegative initial condition, Model (1) admits a unique solution, which remains nonnegative for t ≥ 0. Moreover, there exist constantsC > 0 and T > 0 such that for t ≥ T . We now discuss the existence of equilibria. It it straightforward to see that Model (1) always has a trivial equilibrium E 0 = ( b µ , 0). Next we seek positive equilibria.
However, conclusion (iii) is incomplete. In fact, under the conditions p 3 < 0 and p 1 < 0, it is possible for (1) to have exactly two positive equilibria, which corresponds to the occurrence of saddle-node bifurcations. Applying [10, Lemma 12], we can give more precise conditions to specify the existence of equilibria for Model (1).
is nondecreasing. Consequently, Model (1) has a unique positive equilibrium. We next assume that p 2 Summarizing the above analysis, we have the following result on the existence of equilibria for Model (1). Theorem 2.1. Consider Model (1) with p 1 , p 2 and p 3 defined in (5). There is always a trivial equilibrium E 0 = ( b µ , 0). Furthermore, the following conclusions hold: Case i: p 3 ≥ 0.: Model (1) has no positive equilibrium; Case ii: p 3 < 0, p 1 ≥ 0.: Model (1) has a unique positive equilibrium; Case iii: It is seen from the above analysis that Model (1) undergoes a saddle-node bifurcation at h(v − )h(v + ) = 0. However, it is very difficult to analytically find the bifurcation curves in two parameter space. Next we geometrically calculate the bifurcation curves in (θ, b) space. It follows from (3) that a positive equilibrium of Model (1) corresponds to an intersection point of the two curves defined by and The condition for the occurrence of a saddle-node bifurcation is that the two curves intersect tangentially. Consequently, we require and Note that (9) reduces to As it is impossible to explicitly express θ in terms of b by solving (8) and (10), we locate the bifurcation curves in the parametric form by expressing both θ and b in terms of v. From (8) and (10), we obtain where g 1 (v) = µv 2 − 2aξ 2 v − µξ 2 . The positivity of θ requires that v must be restricted to the range v > v * , where v * is the unique positive root of g 1 (v) = 0. That is v * = aξ 2 + ξ a 2 ξ 2 + µ 2 µ .
If we fix parameter values for a, µ, ε, δ, then (11) defines the bifurcation curves. Figure 1 gives the resulting saddle-node bifurcation curves. Model (1) has exactly two positive equilibria when θ and b are located on the bifurcation curves, there are three positive equilibria when (θ, b) is between the bifurcation curves. If (θ, b) is in Region II, then besides the trivial equilibrium, there is a unique positive equilibrium, while if (θ, b) is in Region I, then the trivial equilibrium is the only equilibrium of Model (1). Region I is obtained from p 3 ≥ 0 ⇐⇒ b ≤ δµ aε . In view of b(v) and θ(v), a direct calculation yields that and where g 2 (v) = µv 3 −3aξ 2 v 2 −3µξ 2 v +aξ 4 . Note that g 2 (v) = 3g 1 (v) and g 1 (v * ) = 0. Thus g 2 (v) > 0 for v > v * . This implies that g 2 (v) is increasing for v > v * . A calculation shows that g 2 (v * ) = −2 µ + a 2 ξ 2 µ v * < 0. Thus g 2 (v) has a unique zero v =v ∈ (v * , ∞). This implies that b(v) and θ(v) are decreasing for v ∈ (v * ,v) and increasing for v >v and at v =v, the derivatives of b(v) and θ(v) both vanish. This implies that the two bifurcation curves meet tangentially at the point C = (b(v), θ(v)), which gives a cusp point, and at C a codimension-2 bifurcation occurs.
3. Stability of equilibria. In this section, we discuss the stability of equilibria for Model (1). At any equilibrium (u, v), the Jacobian matrix reads as About the trivial equilibrium E 0 , we have the following stability result.
Note that if p 3 ≥ 0, then the trivial equilibrium E 0 is the only equilibrium. By the index theory (see [11]), if there is a closed trajectory, then it must enclose the trivial equilibrium. As a result, the closed trajectory must intersect the u-axis. However, for any initial condition (u(0), v(0)) = (u(0), 0), the resulting solution (u(t), v(t)) satisfies u(t) → b/µ and v(t) ≡ 0. That is, the u-axis consists of solution trajectories. Thus, there is no closed trajectory. Therefore E 0 attracts all solutions. Remark 1. If abε < δµ (i.e., p 3 > 0), then the local stability and global attractivity together imply that the trivial equilibrium E 0 is indeed globally asymptotically stable.
At any positive equilibrium (û,v), the Jacobian matrix can be written as: The corresponding characteristic equation is given by Clearly ifv ≤ ξ, then T < 0 and D > 0 and the positive equilibrium (û,v) is locally asymptotically stable. Standard stability analysis also gives the following result on the stability of a positive equilibrium (û,v) withv > ξ.
Since nutrient concentrations play an important role in the phytoplanktonnutrient interactions, we first vary the nutrient inflow rate b to explore how the nutrient inflow rate affects the dynamics. The bifurcation diagram shown in Figure 2 suggests that Model (1) undergoes two saddle-node bifurcations at b = b 1 ≈ 2.25 and b = b 2 ≈ 2.37, respectively. When b ∈ (b 1 , b 2 ), Model (1) admits three positive equilibria. Such saddle-node bifurcations provide a trigger for the outbreak of phytoplankton blooms as the equilibrium level of phytoplankton can suddenly change from a small value to a much larger value. As discussed above, Model (1) can admit three positive equilibria under certain conditions. The phase portrait shown in Figure 3 illustrates the dynamics of Model (1) when three positive equilibria (E 1 , E 2 and E 3 ) and a trivial equilibrium E 0 coexist: E 1 is a stable equilibrium, E 2 and E 0 are saddle points, E 3 is an unstable equilibrium.
Moreover, a Hopf bifurcation occurs at b H ≈ 1.94, periodic solutions emerge for b ∈ (b H , b 1 ). A numerical solution is shown in Figure 4, where the solution of Model (1) approaches a stable periodic solution. In addition, the solution exhibits relaxation oscillations [5,6]: we observe that the nutrient concentration and the phytoplankton abundance exhibit a sudden change in a very short time. This also leads to the outbreak of phytoplankton blooms.
To demonstrate the impact of the toxin effect from TPP, we take θ as the bifurcation parameter. Fixing all other parameters and varying the value of θ, we obtain a bifurcation diagram in Figure 5. As can be seen from Figure 5, increasing θ results in the decrease of the phytoplankton abundance. When θ 1 ≈ 2.11 and θ 2 ≈ 2.24, Model (1) undergoes two saddle-node bifurcations. A further increase in  θ results in the occurrence of a Hopf-bifurcation at θ H ≈ 2.60. This indicates that toxins produced by TPP have a substantial influence on phytoplankton-nutrient interactions. Control of the rate of toxin production may also play an important role in terminating phytoplankton blooms.

5.
Summary. In this work, we have carried out detailed bifurcation and stability analysis for a phytoplankton-nutrient model (Model (1)). We have derived explicit conditions for the existence of positive equilibria and we found that, besides the trivial equilibrium, Model (1) can have exactly 1, 2 and 3 positive equilibria. Further, we provided a practical approach to calculate the saddle-node bifurcation curves. We also discussed the stability of the resulting equilibria and showed that Model (1) can also have a Hopf bifurcation. The occurrence of saddle-node bifurcations and a Hopf bifurcation indicates that both parameters b and θ play a tremendous role in influencing the dynamics of phytoplankton-nutrient interactions. These bifurcations, together with relaxation oscillations observed from numerical simulations provide some insights on understanding the occurrence and control of phytoplankton blooms.