DYNAMICS OF HARMFUL ALGAE WITH SEASONAL TEMPERATURE VARIATIONS IN THE COVE-MAIN LAKE

. In this paper, we investigate two-vessel gradostat models describ-ing the dynamics of harmful algae with seasonal temperature variations, in which one vessel represents a small cove connected to a larger lake. We ﬁrst deﬁne the basic reproduction number for the model system, and then show that the trivial periodic state is globally asymptotically stable, and algae is washed out eventually if the basic reproduction number is less than unity, while there exists at least one positive periodic state and algal blooms occur when it is greater than unity. There are several types of productions for dissolved tox- ins, related to the algal growth rate, and nutrient limitation, respectively. For the system with a speciﬁc toxin production, the global attractivity of positive periodic steady-state solution can be established. Numerical simulations from the basic reproduction number show that the factor of seasonality plays an important role in the persistence of harmful algae.


1.
Introduction. Harmful algal blooms emerged as an important water quality issue in recent decades, and have apparently increased in frequency and intensity worldwide, in both coastal and inland waters. For example, blooms of the haptophyte algae Prymnesium parvum have become common in western Texas and other parts of the American Southwest, where it is referred to as golden algae [21,24]. Blooms of this species were documented to cause large fish mortalities. However, it has been known that water flow in reservoirs can wash out the populations of planktonic algae, and their toxins [9]. This raises many ecological paradoxes [19]. Indeed, the shorelines and the bed can retard flow, producing slow-flowing regions 314 FENG-BIN WANG, SZE-BI HSU AND WENDI WANG which provide a hydraulic storage zone, so that the persistence of algae is enhanced (e.g., [7,Section 3.3

]).
A recent study suggests a potential technique that is possible for managing and mitigating harmful algal blooms through flow manipulations in some river systems [17,18,20]. This fact motivates the theoretical exploration of harmful algal dynamics in riverine reservoirs considered in [9]. To understand longitudinal patterns arising along the axis of flow, advection-dispersion-reaction systems were proposed to incorporate the effects of spatial variations of harmful algae and its toxin production and decay, in riverine reservoirs [9]. The models are one-dimensional systems with simple habitat geometry and transport processes [9], and they were analyzed in [13].
In order to study differences between a single fringing cove and a main lake, the authors in [9] also proposed two-vessel gradostat models with constant volumes to represent the dynamics of harmful algae that can excrete a dissolved toxin. Indeed, fringing coves along the shoreline of riverine reservoirs also provide a storage zone that promotes the persistence of both algae and their toxins. In this paper, we shall give rigorous proofs in the two compartment gradostat models with the influences of seasonal temperature variations, that is, we incorporate the periodic time dependence in the parameters of the model system.
Following [9], we first describe model settings and assumptions. The limiting nutrient for algal growth enters the main lake and the single cove at a constant concentration R in 1 (t) and R in 2 (t), respectively. To simplify the model, constant volume is assumed and the system dilution rate is denoted by D(t) and the system exchange rate between the main lake and the cove is denoted by E(t). We assume that ψ represents the fraction of total volume in cove and φ represents the fraction of total inflow entering cove. In the following, we set   i (t), α i (t), and β i (t), i = 1, 2, is to save notations when we describe our governing systems (1.3) and (1.5) later.
The nutrient uptake function is assumed to be a Monod function of the limiting nutrient concentration (R) at a given location: Here, µ max (t) represents the maximal growth rate; K represents the half saturation constant; the mortality of algal is assumed to be a constant rate m. To simplify our model, we suppose that the nutrient content of algae that die is instantaneously and locally recycled [9]. Further, q N represents nutrient quota of algae. To reflect temporal variations, we assume that there exists a τ > 0 such that R 2 (t), α i (t), β i (t) and µ max (t) satisfy 1 (t + τ ) = R   2 (t), α i (t + τ ) = α i (t), β i (t + τ ) = β i (t), µ max (t + τ ) = µ max (t), i = 1, 2.
The following governing system represents the general form of two models, which is true for many flagellate toxins [9]: where R 1 (t), N 1 (t) and C 1 (t) denote dissolved nutrient concentration, algal abundance and dissolved toxin concentration at time t in the main lake, respectively; R 2 (t), N 2 (t) and C 2 (t) denote dissolved nutrient concentration, algal abundance and dissolved toxin concentration at time t in the cove, respectively. The terms kC 1 and kC 2 stand for toxin degradation. There are two types of productions for dissolved toxins. The first assumes that the algae produce toxin more rapidly when there is little nutrient in the system [9], where is a constant coefficient. It has been observed that toxins produced by Prymnesium parvum (toxic flagellates) are proportional to the degree of algal nutrient limitation ( [3,8,15,16]). The second type of toxin production assumes that the toxin is produced proportional to the algal productivity, This case assumes that toxin is produced in proportion to other cellular products and released into the water at a constant rate [9]. We refer to this as the case of cylindrospermopsin [6,12], which is a cyanotoxin produced by a variety of freshwater cyanobacteria. It was also known that cyanobacteria excrete some toxins that contain nitrogen, a potential limiting nutrient for algae. Hence, chemical decomposition of the toxin results in nutrient recycling [9]. We assume that represents a dimensionless coefficient that specifies the allocation to toxin production [9]. The nutrient content of the toxin is denoted by q C and then the governing equations take the form [9]: (1.5) The terms kq C C 1 and kq C C 2 in the first two equations of (1.5) reveal that the toxin can get recycled back into the system as available nutrient. From the third and fourth equations of (1.5), we realize that only a part, (1 − ), of the nutrient consumed is used for algal growth, which is discounted by the cost of toxin production.

FENG-BIN WANG, SZE-BI HSU AND WENDI WANG
The organization of the paper is as follows. The model analysis of systems (1.3) and (1.5) are presented in sections 2 and 3, respectively. The simulation results are presented in section 4. In section 5, we give some biological interpretations on the analysis of the models.
2. Mathematical analysis of system (1.3). This section is devoted to the study of the global attractivity of system (1.3). We first show that R 6 + is positively invariant for (1.3). For any (2.1) Then W 1 (t) and W 2 (t) satisfy the following coupled differential equations We have the following results concerned with the global dynamics of (2.2): Lemma 2.1. The system (2.2) admits a unique positive τ -periodic solution (W * 1 (t), W * 2 (t)) which is globally attractive in R 2 , that is, for any (W 1 (0), W 2 (0)) ∈ R 2 , we have Proof. We first consider the following system: Let From (1.1), it is easy to see that H(t) is irreducible, cooperative, and the sum of each row in H(t) is negative, which imply that the Floquet multipliers of H(t) are both less than 1 (see, e.g., [ The solution of (2.2) can be expressed as where T (t) is the solution semigroup generated by (2.3). It is easy to see that (W 1 (t), W 2 (t)) is a τ -periodic solution of (2.2) if and only if (2.4) From above discussions, it follows that r(T (τ )) < 1. This implies that I − T (τ ) is invertible, and hence, (2.2) admits a unique τ -periodic solution (W * . Then (Ŵ 1 (t),Ŵ 2 (t)) satisfies (2.3), and hence, lim t→∞ (Ŵ 1 (t),Ŵ 2 (t)) = (0, 0). Thus, we complete the proof.
It follows from Lemma 2.1 that R i (t) and N i (t) are ultimately bounded, i = 1, 2. From the last two equations of (1.3), it is easy to see C 1 (t) and C 2 (t) are also ultimately bounded. Thus, the solutions of system (1.3) exist globally on the interval [0, ∞). Let P : R 6 + → R 6 + be the Poincaré map associated with system (1.3), that is, is the unique solution of system (1.3). Then P is point dissipative (i.e., ultimately bounded) in R 6 + , and P is compact. By [11,Theorem 3.4.8], it follows that P : R 6 + → R 6 + has a global compact attractor in R 6 + . We summarize our discussions above: Lemma 2.2. R 6 + is positively invariant for (1.3) and system (1.3) has a unique and bounded solution with the initial value in R 6 + . Further, the Poincaré map associated with system (1.3) admits a connected global attractor on R 6 + which attracts all positive orbits in R 6 + .
2.1. The basic reproduction number of system (1.3). In order to find the trivial solution of (1.3), we assume N 1 = N 2 = 0 in (1.3). Then R 1 and R 2 satisfy the system (2.2) and it follows from Lemma 2.1 that Furthermore, C 1 and C 2 satisfy the following system (2.5) By the same arguments as those in the dynamics of the system (2.3), we are able to show that lim t→∞ (C 1 (t), C 2 (t)) = (0, 0).
Bacaër and Guernaoui [2] investigated a vector-borne disease model with seasonality, and proposed a general definition of the basic reproduction number in periodic habitats. Wang and Zhao [27] further gave a computational formula for periodic compartmental epidemic models and showed that it is a threshold parameter for the local stability of the disease-free periodic solution. Here, we introduce the basic reproduction number R 0 for the periodic compartmental ecological model (1.3). Linearizing system (1.3) at the trivial solution , 0, 0, 0, 0), we get the following cooperative system: Next, we review some basic results related to a monodromy matrix that is needed for our subsequent discussions. Let A(t) be a continuous, cooperative, irreducible, and τ -periodic k × k matrix function. Suppose Φ A(·) (t) is the monodromy matrix of the linear ordinary differential system and is a matrix with all entries positive for each t > 0. By the Perron-Frobenius theorem, r(Φ A(·) (τ )) is the principal eigenvalue of Φ A(·) (τ ) in the sense that it is simple and admits a positive eigenvector. We further have the following results: . Then there exists a positive, τ -periodic function v(t) such that e µt v(t) is a solution of (2.7).
Motivated by (2.6), we define , is the monodromy matrix of the linear τ -periodic differential sys- We assume that φ(s), τ -periodic in s, is the initial distribution of algae individuals. Then F(s)φ(s) is the rate of new population generated by initial fertile algae individuals who were introduced at time s. Given t ≥ s, then Y (t, s)F(s)φ(s) gives the distribution of those fertile individuals who were newly produced at time s and remain in the fertile compartments at time t. It follows that is the distribution of accumulative new individuals at time t produced by all those fertile individuals φ(s) introduced at time previous to t.
Let C τ be the ordered Banach space of all τ -periodic functions from R to R 2 , which is equipped with the maximum norm · and the positive cone (2.10) Then we call L the next generation operator [27], and define the basic reproduction number as the spectral radius of L.
We present the following result, which will be used in the proof of our main result.

It follows from Lemma 2.4 that the trivial solution
In the following, we shall discuss a special case where µ max (t) ≡ µ max and all the parameters in (1.1) are positive constants. Then it is easy to see that is the unique equilibrium for (2.2). Hence, where we may choose f (R) = µmaxR K+R as a classical example. Then the basic reproduction number corresponds to the spectral radius of FV −1 (see [4,5]), that is, By computations, it follows that where det V := (α 1 + m)(β 2 + m) − α 2 β 1 > 0 ( using (1.1) ).
2.2. Global attractivity of system (1.3). We first study the following systems: If we linearize the system (2.15) at the trivial solution (0, 0), we get the linear system (2.6). This means the trivial solution (0, 0) is locally asymptotically stable for (2.15) if R 0 < 1, and unstable if R 0 > 1. From the biological view of point, the feasible domain Λ(t) for (2.15) should be It is not hard to see that for any (N 0 1 , N 0 2 ) ∈ Λ(0), system (2.15) has a unique solu- is the unique positive τ -periodic solution of (2.15). Proof. Let P : X → X be the Poincaré map associated with system (2.15), that is, is the unique solution of system (2.15). It is easy to see that P : X → X is strongly monotone, and strongly subhomogeneous in the sense that where DP (0, 0) denotes the Fréchet derivative of P at (0, 0), and Φ F(·)−V(·) (t) denotes the monodromy matrix of the linear ordinary differential system (2.6) (see [30,  Since the equations of R i and N i , i = 1, 2 in (1.3) are independent of the equations of C 1 and C 2 , we first study the following subsystem: (2.17) We are now in the position to address the dynamics of system (2.17) in the sense of the following theorem.
) be a unique positive τ -periodic solution of the system (2.2). Then the following statements hold. ( is the unique positive τ -periodic solution of (2.15).
Proof. We rewrite the system (2.17) as follows: where W 1 and W 2 are defined in (2.1). Let Thus, our claim is true.

FENG-BIN WANG, SZE-BI HSU AND WENDI WANG
From (3.2), we conclude that the limiting systems of (1.5) take the forms It is easy to see that for any (N 0 We will always restrict the initial values of (3.3) in this region Σ(0). The solutions of the system (3.3) are in the following forms: (i) Trivial solution0 := (0, 0, 0, 0) always exists; (ii) There may be additional solutions as well and these must be positive.
Now we are ready to state the main result of this section, which indicates that R 0 is a threshold index for the persistence of algae.
Theorem 3.1. The following statements hold.
By Theorem 3.1, [30, Theorem 1.3.1] and the similar arguments as in Theorem 2.1, we can obtain the dynamics of the full system (1.5).
Theorem 3.2. The following statements hold.
Moreover, system (1.5) admits at least one positive τ -periodic solution where (Ñ 1 (t),Ñ 2 (t),C 1 (t),C 2 (t)) is a positive τ -periodic solution of system  averages of R in 1 and R in 2 in (4.2) are given at u = 0, it suggests that the more time heterogeneity of resource inputs, the lower of the basic reproduction number. One more interesting thing from the simulation is that suitable heterogenous resource inputs can reduce the algae growth/bloom threshold below unity under the condition that the time averages are fixed. Now, we select K = 0.2, take and fix the other parameters as in (4.1). Through the synthetic effects of heterogenous resource inputs and heterogenous nutrient conversions, the algae growth/bloom threshold becomes an increasing function of u (see Figure 4.2). This means that the higher heterogenous conversion rate, the more risk algae blooming. In contrast, and fix R in i (t) by (4.3). Simulations indicate that R 0 is below 1 when u > 0.5 (see Figure 4.4). In the rest of this section, we focus on the numerical simulations of system (1.5) to demonstrate influences of allocation coefficient and recycling coefficient q C on the dynamics of (1.5). For illustration purpose, we take K = 0.2, q N = 0.0139, k = 0.1, and use the parameter values in (4.1) and the time-dependent coefficients: R in 1 (t) = 0.0056(1.01 + cos(2πt/365)), R in 2 (t) = 0.005(1.01 + cos(2πt/365)), µ max (t) = 0.15(1.01 + cos(2πt/365)). (4.5) In the first step, we fix k = 0.1 and then take = 0.012, 0.022, 0.032 respectively. Numerical simulations indicate that algae concentrations decreases as increases, which results in the decrease of toxin (see Figure 4.5). Further simulations indicate that C 2 increases as increases in the left of = 0.012 and approximates the maximum at = 0.012 (see Figure 4.6). It should be noted that the enhance of allocation coefficient adds more production coefficient of toxin, but lowers the production of algae. Since algae are the main pool to produce toxin in the former case, the concentration of toxin is reduced. However, the production coefficient of toxin plays the major role in the latter, thus giving more toxin concentrations when increases in the left of = 0.012.
In the second step, we fix = 0.01 and take q C = 0.01, 0.02, 0.03 respectively. Numerical calculations show that toxin concentrations are reduced due to larger inhibition coefficient q C in the production rates of toxins, but this coefficient gives little effect on the concentrations of resources for large time (see Figure 4.7).  5. Discussion. In this paper, we analyze two mathematical models which were established to represent the growth of populations of toxic flagellates or cyanobacteria in the cove-main lake. Large lakes often have a number of smaller adjoining coves, and that previous models either treat these as an ensemble, or focus on a single cove [9]. For the approach where an ensemble of fringing coves is represented as a hydraulic storage zone, the authors in [9] proposed an advection-dispersion-reaction system to model the interactions of nutrient, algae and their toxins in a riverine reservoir. Those mathematical analyses of this model are given in [13]. Here, we concentrate on the study of another approach, that is, we focus on the investigation of the dynamics of nutrient, algae and their toxins with exchange between a single cove and a large lake. To incorporate the influence of seasonal temperature variations on non-steady dynamics, the coefficients in the model systems are time-dependent.
For many flagellate toxins, the toxin contains little or none of the limiting nutrient and the model (1.3) is more appropriate. Since the first four equations of (1.3) are independent of the equations of C 1 and C 2 , this makes the analyses easier and the results more complete. We first define the basic reproduction number, R 0 , for the system (1.3) in (2.11). Then, we are able to prove that the washout periodic state is globally stable if R 0 < 1, and there exists a globally stable coexistence    periodic state if R 0 > 1 (see Theorem 2.1 and Theorem 2.2) by appealing to the theory of monotone dynamical systems and chain transitive sets (see, e.g., [23,26]). Although R 0 is abstractly given in (2.11), we can calculate it numerically (see our numerical simulations). For the special where the coefficients of (1.3) are all positive constants ( i.e., the system (1.3) is temporally homogenous ), R 0 is reduced into the one defined in (2.13). For this homogenous case, it follows that R 0 is the largest real eigenvalue of FV −1 by (2.13) and the Perron-Frobenius Theorem (see, e.g., [23,Theorem 4.3.1]), where FV −1 is given in (2.14).
The toxin production of cyanobacteria (cylindrospermopsin case) is assumed to be proportional to the product of growth rate and abundance, and the governing system is (1.5). In (3.7), we also define the basic reproduction number, R 0 , for the system (1.5). Then we determine the threshold dynamics of system (1.5), namely, the trivial periodic state is globally asymptotically stable, and the algae will be washed out eventually if R 0 < 1, while there exists at least one positive periodic state and the algae can persist continuously when R 0 > 1 (see Theorem 3.1 and Theorem 3.2). In general, the uniqueness and stability of this positive periodic state remain open. However, we point out that we are able to prove the positive periodic state is unique and globally stable for system (1.5) under the additional assumption: m (algal mortality) = k (toxin decay).
Mathematically, (1.5) is more tractable if we further impose this extra assumption (mathematical proofs not shown). When the coefficients of (1.5) are all positive constants, we can also show that R 0 is the largest real eigenvalue of F V −1 := (1 − )FV −1 , where FV −1 is given in (2.14).
There might be opportunities to mitigate algal blooms in some reservoirs where toxic blooms have occurred. It was recognized that mathematical modeling can provide assistance in management of harmful algae within coves through flow manipulations. Here, we rigorously determine a threshold index, the basic reproduction number for the algae, and show that it can predict the algal persistence or extinction for the two-compartment model of algal dynamics with temporal variations proposed by the authors in [9]. The basic reproduction number for the algae involves: the main-cove characteristics (i.e. D(t) and E(t)); the inflowing nutrient concentration ( R in 1 (t) and R in 2 (t) ); a value of the nutrient uptake rate (or growth rate) f (t, R); and the value of period τ . We numerically compute the basic reproduction number and show the influences of seasonality on algae growth in Section 4. This study provides a computable formula that predicts the algal dynamics with temporal variations.