SHADOW SYSTEM APPROACH TO A PLANKTON MODEL GENERATING HARMFUL ALGAL BLOOM

. Spatially localized blooms of toxic plankton species have negative impacts on other organisms via the production of toxins, mechanical damage, or by other means. Such blooms are nowadays a worldwide spread environmental issue. To understand the mechanism behind this phenomenon, a two- prey (toxic and nontoxic phytoplankton)-one-predator (zooplankton) Lotka-Volterra system with diﬀusion has been considered in a previous paper. Nu- merical results suggest the occurrence of stable non-constant equilibrium solutions, that is, spatially localized blooms of the toxic prey. Such blooms appear for intermediate values of the rate of toxicity µ when the ratio D of the dif-fusion rates of the predator and the two prey is rather large. In this paper, we consider a one-dimensional limiting system (we call it a shadow system) in (0 ,L ) as D → ∞ and discuss the existence and stability of non-constant equi- librium solutions with large amplitude when µ is globally varied. We also show that the structure of non-constant equilibrium solutions sensitively depends on L as well as µ .

1. Introduction. The term harmful algal bloom (in short, HAB) indicates an algal bloom that has negative impacts on other organisms via the production of toxins, mechanical damage, or by other means. HAB includes different types of taxa such as dinoagellates, diatoms, and cyanobacteria (commonly known as blue-green algae). The latter are of special importance because of their potential impact on drinking or recreational waters. In fact, they can produce a variety of potent toxins called cyanotoxins (e.g. Falconer and Humpage [3]). These compounds have been found to be hepatotoxic or neurotoxic for a wide range of organisms, including humans, and several intoxication cases have been reported worldwide (Jochimsen et al. [7]). Therefore, in the recent years, the formation of toxic blooms of cyanobacteria in lakes and rivers has been causing more and more concern. Ecological evidence suggest that toxic and nontoxic species of freshwater phytoplankton hardly coexist in absence of other species. In particular, competition experiments have shown that 830 HIDEO IKEDA, MASAYASU MIMURA AND TOMMASO SCOTTI the toxic strain of Microcystis is a very poor competitor for light (Kardinaal et al. [9]). In these experiments the toxic strain always lost the competition against the nontoxic one, even when given a strong initial advantage. Then a natural question is: how can these species survive and actually bloom? Regarding this question, it is ecologically discussed that toxin-producing Microcystis has overall an inhibitory effect on the growth of most herbivore taxa. Nevertheless, zooplankton usually grazes on both toxic and nontoxic species (Fulton and Pearl [5]). This is interesting, since the toxic or noxious chemicals produced by blue green algae may inhibit feeding and, over long term, cause mortality of zooplankton (Porter and Orcutt [13], Lampert [10], [11], Fulton and Paerl [4]). In particular, while a few species like the rotifer Brachionus calyciorus and the cladoceran Bosmina longirostris apparently make no great distinction between toxic and non toxic prey, the feeding rates of other small-bodied cladocerans, rotifers, and copepods seem to be strongly related to the toxicity of Microcystis (Fulton and Paerl [4]). These observations suggest that predator and toxic prey have an inhibitory effect on each other. Then another naive question is: can the existence of such interaction promote the spatial pattern formation and local algal blooms? This ecological question motivates us to theoretically understand the mechanism behind the formation of spatial blooms of toxic plankton. For this purpose, we proposed a two-prey (toxic and nontoxic phytoplankton)-one-predator (zooplankton) Lotka-Volterra system with diffusion in a previous paper Scotti et al. [14], under the assumptions that (i) in absence of the predator, the toxic prey is a weaker competitor for common resources than the nontoxic one, and (ii) depending on the toxicity and another parameter, the predator is more or less inclined to avoid the toxic prey. The dynamics discussed above are described by the following reaction-diffusion system: where P (t, x), X(t, x) and Z(t, x) are respectively the population densities of the nontoxic and toxic phytoplankton and of the zooplankton, and the parameters r i , D i (i = 1, 2, 3), K 1 , K 2 , a, b, d and µ are all positive constants. Ecologically, an important parameter in (1) is µ, which is the rate of toxicity. We assume that the predation rate of the toxic prey d decreases in µ. Here we simply specify d = d(µ) = 1/(1 + µ). The ecological explanation of (1) is discussed in [14]. By using suitable transformations of time and space, (1) can be rewritten as where r = r2/r 1 , R = r 3 /r 1 , σ = D 2 /D 1 and D = D 3 /D 1 are respectively the ratios of the growth rates and the diffusion rates of the nontoxic phytoplankton and the toxic phytoplankton and/or the zooplankton.

2.
A mathematical model. In this paper, simply assuming that r = 1, K 1 = K 2 = K and σ = 1 in (2), we consider the one dimensional system of (2) in a finite interval 0 < x < L, that is, where a, b, K, µ, R and D are positive constants. We treat (3) under the Neumann boundary conditions and the initial conditions For the system (3), we impose the following assumptions: which implies that predator and nontoxic prey coexist in the absence of toxic prey.
(A2) a < 1 < b, which implies that in the absence of predator, the nontoxic prey is a superior competitor for resources with respect to the toxic one. Assumption (A1) corresponds to the coexistence of two species of (nontoxic) phytoplankton and zooplankton (e.g. Hutchinson [6]). Assumption (A2) is based on the idea that toxic strains are eventually outcompeted by nontoxic ones (e.g. Kardinaal et al. [9], Lampert [10]). For (3) with (4), we easily find that the spatially constant equilibrium solution E 3 = (1, 0, (K − 1)/K) exists for any µ > 0. Instead of (A1), we assume K to satisfy Putting , we know that E 3 is stable for 0 < µ < µ c , while it is unstable for µ c < µ. Furthermore, when µ increases, a spatially positive constant equilibrium solution of (3) with (4), say, E 4 = (P µ ,X µ ,Z µ ) super-critically bifurcates from E 3 at µ = µ c , that is, E 4 exists for µ > µ c . Here we assume R suitably large to require that E 4 is stable for µ > µ c in the ODEs corresponding to (3) in the absence of diffusion. In addition to (A2) -(A3), we assume µ to satisfy Here we note that E 4 is not necessarily stable in (3) with (4), that is, the stability of E 4 depends on µ, R, D and L. Figures 1(a) and (b) show respectively the linearized stability of E 4 in the (D, µ) plane (L = 30) and (L, µ) plane (D = 2500) for suitably fixed a, b, K and R ( [14]). When µ(> µ c ) is small, E 4 is stable for any fixed D (or L), while when µ is suitably large, E 4 loses its stability as D (or L) increases. This   Figure 1. Solid (resp. dashed) lines represent stable (resp. unstable) equilibrium solutions of (3) with (4). The right figure is a magnification of the left one where µ is close to µ c1 .
destabilization is called diffusion-induced instability as stated by Turing ([15]). In fact, suppose that D is suitably large (for instance, D = 2500). By suing AUTO, we can show the global structure of equilibrium solutions of (3) with (4) when µ is varied (see Figure 2). This figure shows the occurrence of two bifurcation values of The other parameters are the same as the ones in Figure 1 and D = 2500. HereP + 1 ,X + 1 andZ + 1 are drawn in blue, green and red colors, respectively. The other parameters are the same as the ones in Figure 1. Solid (resp. dashed) lines represent stable (resp. unstable) equilibrium solutions of (3) with (4).
We now address the following question: can we show the existence and stability of such non-constant equilibrium solutions of (3) and (4) analytically? One of the ways is to assume that D is rather large in (3) so that the problem (3) with (4) and (5) becomes simpler, because one can expect Z(t, x) to be spatially homogeneous. In Figure 4, we show the dependency of D on the structures of equilibrium solutions of (3) with (4) when µ is globally varied. We notice that they do not qualitatively change if D becomes larger and larger, as shown in Figures 4(c) -(e). This result motivates us to study the limiting system of (3) with (4) as D → ∞ in order to discuss the existence and stability of non-constant equilibrium solutions.
3. The shadow system. We formally derive the limiting system from (3) with (4) as D → +∞. We first start by dividing the third equation of (3) by D so that we obtain 1 If we assume that Z t , P , X and Z remain bounded as D → +∞, then (6) implies that the limit of Z, say, ξ satisfies It follows from (7) and the boundary conditions (4) that ξ must be spatially constant. On the other hand, integrating the equation of Z in (3) on the interval [0, L] with respect to x, we get by the boundary conditions (4) and then that is, (7) and (8). Consequently when D → +∞, we formally obtain the following limiting system, which is called a shadow system, of (3) with (4) for the unknowns (P (t, x), X(t, x), ξ(t)) : with (P x , X x )(t, 0) = (0, 0) = (P x , X x )(t, L), t > 0. (10) In order to obtain equilibrium solutions of (9) with (10), we first assume ξ to be suitably fixed and consider the first two equations of (9) with (10), which are the well known competition-diffusion system for two species. If we could find an equilibrium solution (P (x; ξ),X(x; ξ)) of the system and substitute it into the third equation of (9), we obtain the equation of ξ only If we can findξ to satisfy (11), we get an equilibrium solution (P (x;ξ),X(x; ξ),ξ) of the shadow system (9) with (10). In order to obtain equilibrium solutions of the competition-diffusion system for two species, we can apply the results by Kan-on [8]. For this purpose, we introduce a new parameter ε = 1/L 2 into (9) with (10) to normalize the interval [0, L] to the unit one [0, 1]. Then, the shadow system (9) with (10) is rewritten to the rescaled normal system with a parameter ε as follows: In this section, we fix µ arbitrarily to satisfy µ > µ c and ignore the µ-dependency of equilibrium solutions of the shadow system (9) with (10).
3.1. 2-component competition system of (12) with (13) for fixed ξ. In this subsection, we look for equilibrium solutions of the following 2-component competition system of (12) with (13) for a fixed ξ when ε is a free parameter: Here we note that, depending on the parameter values a, b, K, µ and ξ, the dynamics of solutions of (14) possesses three different cases, (a) mono stability, (b) bistability and (c) coexistence, as shown in Figure 5. Since our concern is to find non-constant equilibrium solutions of (14), we restrict values of the parameters a, b, K, µ and ξ to satisfy the case (b) in Figure 5, because this is the only case for which such solutions exist (Zhou and Pao [17] , Kan-On [8]). Figure 5 We define the interval I by  Figure 5. Three different structures of the nullclines of (14) with d = d(µ) = 1/(1 + µ). (a-1) P -monostability, (a-2) Xmonostability, (b) bistability and (c) coexistence. The red and white circles stand for stable and unstable equilibrium solutions of (14), respectively. and assume ξ to satisfy ξ ∈ I. Then the system (14) has a positive constant equilibrium solution (P (ξ),X(ξ)) which is explicitly represented as Remark 1. We easily find that Taking ε as a bifurcation parameter, we look for non-constant equilibrium solutions of (14), which are bifurcated from the constant solution (P (ξ),X(ξ)). By simple calculation, it can be easily checked that for arbitrarily fixed ξ ∈ I, the linearized operator of (14) around (P (ξ),X(ξ)) has the zero eigenvalue at ε = ε n 0 (ξ) = ε n=1 n=2  Figure 6. Schematic global structure of the constant and non-constant equilibrium solutions (P (ξ),X(ξ)) and (P ± n (x; ξ, ε),X ± n (x; ξ, ε)) (n = 1, 2, · · · ) of (14) when ε is varied.

4.
Dependency of equilibrium solutions of the shadow system on the toxicity µ. In this section, by suitably fixing the interval length L and using µ as a free parameter, we consider (9) with (10). Since Figure 4(e) suggests that when D is sufficiently large, 1-mode equilibrium solutions of (9) with (10) are exist and are stable, we restrict the dependency of µ to 1-mode equilibrium solutions Π + 1 (x; 1/L 2 , µ) of (9) with (10) for different values of L.
For the case of L = 35, AUTO exhibits the structure of 1-mode equilibrium solutions which is almost similar to the one for the case of L = 30 in a sense that the 1-mode equilibrium solutions super-critically bifurcate, as shown in Figure  14. This result is verified exactly in Section 4.2. However, it exhibits an S-shaped structure with respect to µ so that there exist three 1-mode equilibrium solutions for suitable range of µ. The difference on global structures of equilibrium solutions for L = 30 and L = 35 is well-consistent with the ones obtained from the information on H + 1 (ξ; 1/L 2 , µ) in Figures 10 and 11, respectively. Let us show one more important example. It is the case L = 60 where H + 1 (ξ; 1/L 2 , µ) is shown in Figure 15. We can observe one hump in H + 1 (ξ; 1/L 2 , µ) so that there are two solutions ξ of H + 1 (ξ; 1/L 2 , µ) = 0 in a suitable range of µ. AUTO exhibits the corresponding global structure of equilibrium solutions in Figure 16. It is emphasized that 1-mode equilibrium solutions sub-critically bifurcate from the constant one E 4 at µ = µ c1 so that there are two 1-mode equilibrium solutions for suitable range for µ.
(3) The number and stability of 1-mode equilibrium solutions also depends on the interval length L.
The assertion (1) is proved in Theorem 3.1. In the next subsection, the assertion (2) will be discussed by using the local bifurcation theory ([1]).

4.1.
Direction of the pitchfork bifurcation. In the above, we have numerically shown the existence of the n-mode equilibrium solutions of the shadow system (12) with (13) when L is varied. Here we are interested in showing the local stability of the 1-mode equilibrium solutionsΠ ± 1 (x; 1/L 2 , µ) of (12) with (13). In order to do it, we can apply the local bifurcation theory [1] to the stationary problem of Figure 14. Bifurcation diagram of the shadow system (9) with (10) when µ is a free parameter. Other parameters are fixed at a = 0.95, b = 1.2, K = 2.9, R = 0.43 and L = 35. Solid (resp. dashed) lines represent stable (resp. unstable) equilibrium solutions of (9) with (10). In the right corner, it is shown a magnification around the primary bifurcation point (n = 1).
Though Theorem 3.1 is identical to this proposition limited to a neighborhood of the bifurcation points, the representation of the solution is different. For example, (P + n (x; ε),X + n (x; ε),ξ + n (ε)) in Theorem 3.1 corresponds to (P n (x; s),X n (x; s),ξ n (s)) for s ∈ (0, δ) in Proposition 2. For the proof of Proposition 2, we refer to the paper by Q. Wang, C. Gai and J. Yan [16], so we omit it.   (9) with (10). In the right corner, it is shown a magnification around the primary bifurcation point (n = 1).
Remark 2. If K 2 (n) < 0, the pitchfork bifurcation of the n-mode solutions is supercritical, on the other hand, if K 2 (n) > 0, it is subcritical.

4.2.
Bifurcation of the n-mode equilibrium solutions from (P * ,X * ,ξ * ). In this subsection, we study the relation between the sign of the principal eigenvalue of the linearized problem (see (29)) and that of K 2 (n) of the n-mode equilibrium solutions (P n (x; s),X n (x; s),ξ n (s)) of the following time dependent system: (P x , X x )(t, 0) = (0, 0) = (P x , X x )(t, 1), t > 0.
For this purpose, we apply the classical results by Crandall and Rabinowitz [2] on the linearized stability to (28). Then we consider the following linearized eigenvalue problem around (P n (x; s),X n (x; s),ξ n (s)): where λ(s; n) is an eigenvalue and (P (x; s; n), X(x; s; n), ξ(s; n)) is the corresponding eigenfunction of the linearized operator.
The proof will be stated in Section 5.3. Noting the order relation of the bifurcation values ε = ε n 0 (ξ * ) for each n ∈ N, we easily find that the 1-mode equilibrium solutions (P 1 (x; s),X 1 (x; s),ξ 1 (s)) of (28) bifurcate primarily from the constant equilibrium solution (P * ,X * ,ξ * ). Then we have Corollary 1. Assume (A2) -(A5). In a neighborhood of the 1-mode bifurcation value ε = ε 1 0 (ξ * ), the 1-mode equilibrium solutions (P 1 (x; s),X 1 (x; s),ξ 1 (s)) of (28) are asymptotically stable (resp. unstable) if their bifurcation is supercritical (resp. subcritical). On the other hand, for n 2 the n-mode equilibrium solutions (P n (x; s),X n (x; s),ξ n (s)) of (28) are unstable in a neighborhood of the n-mode bifurcation value ε = ε n 0 (ξ * ). Finally, we consider the type of bifurcation and the stability of the 1-mode equilibrium solutions in a neighborhood of the bifurcation points. For these purpose, we have to calculate the sign of K 2 (1) (see Theorem 4.1). In Figure 17, the solid curves represent the relation between L and µ on which 1-mode equilibrium solutions bifurcate from the constant equilibrium solution E 4 . The dashed curves represent K 2 = K 2 (1) on which 1-mode equilibrium solutions bifurcate depending on µ. First, for a given L we can get the bifurcation values of µ by using the solid curves and, next, we find the sign of K 2 (1) corresponding to the bifurcation values µ by using the dashed curves. Figure 17 shows the following: when L = 30, 1-mode equilibrium solutions bifurcate at µ = µ c1 and µ = µ c2 . Both of them super-critically bifurcate and are stable near the bifurcation points. When L = 35, a similar situation occurs. On the other hand, when L = 60, 1-mode equilibrium solutions bifurcate at µ = µ * c1 and µ = µ * c2 . Figure 17(c) says that K 2 (1) > 0 at µ = µ * c1 , which means that it sub-critically bifurcates and is unstable near the bifurcation point, while Figure 17(a) says that K 2 (1) < 0 at µ = µ * c2 , which means that it super-critically bifurcates and is stable near the bifurcation point. Of course, these results agree with the ones obtained by AUTO (see Figures 12, 14 and 16).