Deterministic and stochastic nutrient-phytoplankton- zooplankton models with periodic toxin producing phytoplankton

Applied Mathematics and Computation 271 (2015) 52–67

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Deterministic and stochastic nutrient-phytoplanktonzooplankton models with periodic toxin producing phytoplankton Sophia R.-J. Jang∗, Edward J. Allen Department of Mathematics and Statistics, Texas Tech University, Lubbock, TX 79409-1042, USA

a r t i c l e

i n f o

MSC: 92D25 92D40 Keywords: Phytoplankton Zooplankton Uniform persistence Itô stochastic differential equations

a b s t r a c t Deterministic and stochastic models of nutrient-phytoplankton-zooplankton interaction are proposed to investigate the impact of toxin producing phytoplankton upon persistence of the populations. The toxin liberation by phytoplankton is modeled periodically in the deterministic system. We derive two thresholds in terms of the parameters for which both plankton populations can persist if these thresholds are positive. We construct stochastic models with Itô differential equations to model variability in the environment. It is concluded that the input nutrient concentration along with the toxin liberation rate play critical roles in the dynamics of the planktonic interaction. In particular, toxin producing phytoplankton can terminate harmful algal blooms and the planktonic interaction is more stable if either the input nutrient concentration is smaller or if the toxin production rate is larger. © 2015 Elsevier Inc. All rights reserved.

1. Introduction There has been a global increase in harmful algal blooms in the last few decades [1] and a considerable number of scientific studies have been carried out to investigate their causes and effects [2–7]. Harmful Algal Blooms (HABs) are broadly classified into two categories. One group is the toxin producing phytoplankton (TPP), which are the toxin producers that can contaminate seafood or kill higher trophic organisms. The other group is the high biomass producers which can cause anoxia and indiscriminately kill marine life after reaching dense concentrations. Some of the phytoplankton populations can be classified into both of these groups. There are many toxin producing phytoplankton species and we refer the reader to Hallegraeff [1] for a list of TPP species. Mathematical modeling is an important tool to study plankton dynamics [8–10] and to understand the various mechanisms involved in toxin producing phytoplankton. The work in [3] introduces an autonomous system of ordinary differential equations of two trophic levels to study the effects of toxin liberation by phytoplankton. Their simulations conclude that TPP may terminate algal blooms. In a sequel paper of [3], the authors incorporate a discrete delay in the toxin liberation and use a simple mass action for the zooplankton grazing rate [11]. The toxin liberation rate is modeled by a Holling type II function. Due to the discrete delay, the model exhibits persistent oscillations that mimic phytoplankton blooms. The simple mass action assumption in [11] is then extended to the Holling type II function in [12] and it is also assumed that the toxin liberation is modeled by a Holling type II function. In particular, the

∗

Corresponding author. Tel.: +806 834 7006; fax: +806 742 1112. E-mail address: [email protected] (S.R.-J. Jang).

http://dx.doi.org/10.1016/j.amc.2015.08.065 0096-3003/© 2015 Elsevier Inc. All rights reserved.

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

53

half saturation constants of the two rates are equal. The authors derive a threshold below which there are no positive periodic solutions. In [13], a nutrient-plankton reaction-diffusion model with toxin producing phytoplankton is studied. The toxin liberation, however, is not explicitly modeled but built into the zooplankton’s grazing rate. The zooplankton’s functional response decreases with increasing phytoplankton. It is concluded that diffusion can stabilize the interaction by eliminating periodic solutions. The system in [3] is extended to include periodic forcing of toxin liberation and passive diffusion of both populations in [14]. Several thresholds are derived to provide criteria for coexistence of both plankton populations. It is observed numerically in [14] that periodic toxin production by phytoplankton can destabilize the interactions and may induce chaos. However, diffusion of plankton populations can on the other hand stabilize the interactions and promote plankton patchiness. The effect of inflows on population dynamics and assemblage structures has long been of interest to plankton ecologists [15–17]. The magnitude and the timing of inflows influence the physicochemical environment and bring in variations of salinity and nutrient availability that affect phytoplankton growth. In particular, inflows affect plankton community composition and productivity [6,18]. By shifting nutrient availability and alternating hydrology, inflows can affect the incidence of harmful algal blooms over a wide range of taxa, including cyanobacteria, dinoflagellates, diatoms and prymnesiophytes [6]. In addition, several studies demonstrate that inflow-induced changes to the environment can also lead to bloom decline [6]. Other studies suggest that inflows of nutrients alter water chemistry and can increase toxic blooms [18]. Furthermore, populations in natural systems live in a periodically varying environment and population interactions are therefore qualitatively affected by seasonal forcing. Consequently, periodic toxin liberation by TPP can have profound effects on the excitation and succession of algal blooms as demonstrated in [14]. There is also evidence suggesting that toxin release by phytoplankton is not constant and periodicity can either terminate or initiate periodic outbreaks of phytoplankton [5,19]. Motivated by these biological findings, we propose three trophic levels of nutrient-plankton models with general uptake and grazing functions to investigate the effects of toxin production by phytoplankton. The toxin liberation in this study is incorporated directly into the interaction so that the modeling aspect is different from that of [13]. In addition, toxin liberation is modeled periodically in the deterministic system. This periodic toxin production is only considered very recently in [20] and [14], where two trophic levels of nutrient-phytoplankton and of phytoplankton-zooplankton are studied respectively. Demographic stochasticity is low for the dynamical system due to the high population levels of phytoplankton and zooplankton. However, randomly-varying input nutrient concentrations in the environment, in particular, result in variability in the plankton levels. To examine the effects of variability of the environment, we extend this periodic nutrient-plankton system to a stochastic model, where the input nutrient concentration is modeled by a continuous, mean-reverting, stochastic process satisfying an Itô stochastic differential equation. We provide analytical results for both systems and numerically simulate the models. We compare both models and make biological conclusions. In the following section, a deterministic nutrient-plankton model is proposed and analyzed. We first provide analysis for the nutrient-phytoplankton subsystem. We then derive two thresholds in terms of the model parameters and the average value of the toxin liberation rate and prove that both plankton populations can persist if these two thresholds are positive. Numerical investigation with parameter values taken from the literature is also presented. A stochastic model using Itô’s stochastic differential equations is proposed in Section 3. Simulations of the stochastic model with the same parameter values and functional forms as the deterministic system are performed. The final section provides a brief summary and conclusions. 2. The deterministic model The deterministic nutrient-phytoplankton-zooplankton system is presented and analyzed in this section. 2.1. The model and analysis Let N(t), P(t) and Z(t) denote the nutrient concentration, the concentrations of phytoplankton and zooplankton populations at time t, respectively. The two plankton levels are modeled in terms of the nutrient content and therefore their units are nitrogen or nitrate per unit volume. Let δ and μ be the per capita death rate of phytoplankton and zooplankton, respectively. In a natural nutrient-plankton system, water flowing into the system bring an input of nutrients and outflows carry out nutrients [15]. It is assumed that the input nutrient concentration N0 is constant with a constant inflow and outflow rate D. Both plankton populations may also flow out of the system and we let D1 and D2 be the washout rates for phytoplankton and zooplankton respectively, where D, D1 , and D2 are in general different to account for other physical considerations such as sinking of phytoplankton, etc. The phytoplankton’s nutrient uptake h and zooplankton’s grazing f are modeled using general functions and are assumed to satisfy the following hypotheses: (H1) h, f ∈ C1 [0, ∞), h(0) = f (0) = 0, h (x), f (x) > 0 for x ≥ 0, and limx→∞ h(x) = 1 = limx→∞ f (x). The well known Holling types II and III functions clearly satisfy (H1). Let a be the maximal nutrient uptake rate of phytoplankton and let α be the maximal zooplankton ingestion rate. Parameters b and β are the fractions of phytoplankton nutrient conversion and zooplankton grazing conversion, respectively, 0 < b, β < 1. Let  (t) denote the maximum toxin production rate satisfying (H2)  (t ) = γ (1 + Aq(t )), q(t) is τ -periodic, τ > 0, q ∈ C1 [0, ∞), |q(t)| ≤ 1 for t ≥ 0, 0 ≤ A < 1, and γ > 0. If γ = 0, then the phytoplankton is not toxin producing. If γ > 0, then the parameter A in (H2), 0 ≤ A < 1, is the magnitude of the periodic part. The maximum toxin liberation rate is a constant if A = 0.

54

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67 Table 1 Parameters and their biological meanings. Parameters

Biological meanings

D D1 D2 N0 a b

Nutrient input and washout rate Phytoplankton washout rate Zooplankton washout rate Constant input nutrient concentration Maximal nutrient uptake rate by phytoplankton Fraction of nutrient conversion for phytoplankton Maximal zooplankton ingestion rate Fraction of zooplankton conversion Phytoplankton mortality rate Zooplankton mortality rate Maximal toxin liberation rate by phytoplankton

α β δ μ  (t)

The phytoplankton’s toxin production rate g is density dependent and satisfies similar properties as those of f and h: (H3) g ∈ C1 [0, ∞), g(0) = 0, g (x) > 0 for x ≥ 0, and limx→∞ g(x) = 1. In particular, toxin release by phytoplankton reduces the growth of zooplankton with an increase in mortality. This additional mortality depends on the phytoplankton level via g(P) that satisfies (H3). Under the above biological assumptions, the nutrient-plankton interactions are described by the following first-order system of periodic ordinary differential equations

⎧  0 ⎪ ⎨N = D(N − N) − ah(N)P  P = abh(N)P − α f (P )Z − δ P − D1 P ⎪ ⎩Z  = αβ f (P)Z − μZ −  (t )g(P)Z − D Z 2

(2.1)

with nonnegative initial conditions. For example, if  (t ) = γ (1 + A sin (Tt )), 0 < A < 1, then (2.1) models annual, semi-annual, 2π 2π 2π , 180 , 30 , and 2π , respectively. monthly, and daily periodicities of the toxin liberation by phytoplankton when T are 365 We summarize the parameters and their biological interpretations in Table 1. Let

F (t, N, P, Z ) = (F1 (N, P, Z ), F2 (N, P, Z ), F3 (t, N, P, Z )) be the right hand side of (2.1), where

F1 = D(N0 − N) − ah(N)P,

F2 = abh(N)P − α f (P )Z − δ P − D1 P,

and

F3 = αβ f (P )Z − μZ −  (t )g(P )Z − D2 Z. We first provide existence, uniqueness and boundedness of the solutions of (2.1). The global asymptotic behavior of the nutrientphytoplankton subsystem can be readily analyzed. Proposition 2.1. Solutions of (2.1) exist, remain nonnegative and are bounded for t > 0. Proof. Notice F is locally Lipschitz with F1 |N=0 ≥ DN0 > 0, F2 |P=0 = 0 and F3 |Z=0 = 0. It follows from Theorem A.4 of [21] that the initial value problem (2.1) has a unique nonnegative solution defined on some interval [0, tˆ), where tˆ > 0 may depend on initial conditions. We prove that tˆ = ∞. Let S = bβ N + β P + Z. Then S = bβ DN0 − bβ DN − (δ + D1 )β P − (μ + D2 )Z −  (t )g(P )Z ≤ bβ DN0 − bβ DN − (δ + D1 )β P − (μ + D2 )Z ≤ bβ DN0 − mS, where m = min{D, δ + D1 , μ + D2 } > 0. Consider the scalar linear equation x = bβ DN0 − mx with 0 DN0 DN0 x(0) = S(0). Since x(t ) = bβm + (x(0) − bβm )e−mt and S(t) ≤ x(t) for t ≥ 0, we have lim supt→∞ S(t ) ≤ bβmDN . Therefore, so-

DN0 lutions of (2.1) are defined for all t > 0 and are uniformly bounded with lim supt→∞ (bβ N(t ) + β P (t ) + Z (t )) ≤ bβm for all solutions (N(t), P(t), Z(t)) of (2.1). 

It follows from the proof of Proposition 2.1 that system (2.1) is point dissipative. Observe that the N-component of any solution of (2.1) satisfies N ≤ D(N0 − N) and hence solutions of (2.1) satisfy

lim sup N(t ) ≤ N0 .

(2.2)

t→∞

On the other hand, Z (t ) = 0 for t > 0 if Z (0) = 0 and (2.1) reduces to the following two-dimensional autonomous system of nutrient-phytoplankton interactions



N = D(N0 − N) − ah(N)P P  = abh(N)P − δ P − D1 P.

(2.3)

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

55

The subsystem (2.3) always has a boundary steady state (N0 , 0) and its stability depends on the Jacobian matrix J1 of (2.3) evaluated at (N0 , 0) given by





−D −ah(N0 ) . 0 abh(N0 ) − δ − D1

J1 (N0 , 0) =

(2.4)

Therefore, (N0 , 0) is locally asymptotically stable if abh(N0 ) < δ + D1 and unstable if abh(N0 ) > δ + D1 . Furthermore, components ¯ P¯ ) satisfy of an interior steady state (N,

abh(N¯ ) = δ + D1 ,

(2.5)

and

P¯ =

D(N0 − N¯ ) > 0. ah(N¯ )

(2.6)

¯ P¯ ) if and only if We conclude that (2.3) has a unique interior steady state (N,

abh(N0 ) > δ + D1 .

(2.7)

¯ P¯ ) depends on the Jacobian matrix J1 evaluated at the steady state, The stability of (N,



¯ P¯ ) = J1 (N,

−D − ah (N¯ )P¯ abh (N¯ )P¯



−ah(N¯ ) . 0

(2.8)

¯ P¯ ) is −D − ah (N¯ )P¯ < 0 and the determinant of J1 (P, ¯ N¯ ) equals a2 bh(N¯ )h (N¯ )P¯ > 0. Therefore, (N, ¯ P¯ ) is locally The trace of J1 (N, asymptotically stable whenever it exists. Furthermore, we have the following global results for (2.3). Proposition 2.2. The following statements hold for system (2.3). (a) If abh(N0 ) < δ + D1 , then the unique steady state (N0 , 0) is globally asymptotically stable in R2+ . ¯ P¯ ) which is moreover (b) If abh(N0 ) > δ + D1 , then steady state (N0 , 0) is unstable and there exists a unique interior steady state (N, globally asymptotically stable in {(N, P ) ∈ R2+ : P > 0}. Proof. (a) Notice R2+ is positively invariant for (2.3) and solutions of (2.3) are bounded by Proposition 2.1. By the assumption, (N0 , 0) is the only steady state of (2.3) and is on the boundary of R2+ . It follows from the Poincaré–Benxison Theorem that (2.3) has no periodic solutions and (N0 , 0) is globally asymptotically stable. (b) By the assumption abh(N0 ) − δ − D1 > 0, (N0 , 0) is unstable with its stable manifold lying on the nonnegative N-axis. ¯ P¯ ), There exists a unique positive solution N¯ for (2.5), where N¯ < N0 and hence (2.3) has a unique interior steady state (N, where P¯ is given by (2.6). The interior steady state is always locally asymptotically stable. Applying the Dulac criterion with G(N, P ) = 1P for N ≥ 0, P > 0, we have

    D ∂ N ∂ P + = − − ah (P ) < 0 for N ≥ 0, P > 0. ∂N P ∂P P P

Therefore (2.3) has no periodic solutions in {(N, P ) ∈ R2+ : P > 0} by the Dulac criterion [21,22]. Since solutions of (2.3) are bounded by Proposition 2.1, the interior steady state is globally asymptotically stable in {(N, P ) ∈ R2+ : P > 0} by the Poincaré–Bendixson Theorem [23].  Since the explicit time dependency of system (2.1) is in the third trophic level, the analytical results obtained for the NPsubsystem can be applied to both the autonomous NPZ system A = 0 and the periodic NPZ system A > 0. Let

ˆ =

1

τ

τ

0

 (s)ds,

M = maxt∈[0,τ ]  (t ),

m = mint∈[0,τ ]  (t )

(2.9)

ˆ ≤ M and  ˆ = γ if  (t ) = γ is denote the average, maximum and minimum values of  (t), respectively. Then m ≤  a constant. Moreover, m ≥ γ (1 − A) > 0 and M ≤ γ (1 + A) by (H2). Let P : R3+ → R3+ be the Poincaré map of (2.1), i.e., P (N(0), P (0), Z (0)) = (N(τ ), P (τ ), Z (τ )), where X (t ) = (N(t ), P (t ), Z (t )) is the solution of (2.1) with initial condition X(0). Observe that a τ -periodic solution of (2.1) corresponds to a fixed point of P. The stability of any τ -periodic solution depends on the derivative of the Poincaré map P evaluated at the corresponding fixed point [23]. Clearly system (2.1) has a trivial periodic solution E0 = (N0 , 0, 0). Its stability depends on

P  (N, P, Z )|E0 = 0 (τ ),

(2.10)

where 0 (t) is the fundamental matrix of z = J(E0 )z with 0 (0) = I and



J(E0 ) =

−D 0 0

−ah(N0 ) abh(N0 ) − δ − D1 0



0 . 0 −μ − D2

(2.11)

56

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

Define

R0 = abh(N0 ) − δ − D1 .

(2.12)

A direct computation yields

⎛

ah(N0 ) R0 +D

e−Dt

0 (t ) = ⎝ 0 0

(e−Dt − eR0 t ) eR0 t 0

0 0 e−(μ+D2 )t

⎞ ⎠.

Therefore, the Floquet multipliers of E0 are e−Dτ < 1, e−(μ+D2 )τ < 1 and eR0 τ . If R0 < 0, then it can be shown that E0 is globally asymptotically stable for (2.1). Theorem 2.3. If R0 < 0, then the trivial periodic solution E0 = (N0 , 0, 0) is globally asymptotically stable in R3+ for system (2.1). Proof. Let (N(t), P(t), Z(t)) be a solution of (2.1) with P(0) > 0 and Z(0) > 0. Since lim supt→∞ N(t ) ≤ N0 by (2.2), for any  > 0, there exists t0 > 0 such that N(t ) < N0 +  . By the assumption of R0 < 0, we can choose  > 0 so that abh(N0 + ) < δ + D1 . It follows from the second equation of (2.1) that P  (t ) ≤ (abh(N0 + ) − δ − D1 )P (t ) for t ≥ t0 . Thus limt→∞ P (t ) = 0 and it can then be verified that limt→∞ Z (t ) = 0. Therefore, E0 is globally attracting in R3+ and it is globally asymptotically stable.  The threshold R0 can be viewed as the maximal net growth rate of the phytoplankton. The phytoplankton population goes extinct and so does the zooplankton population if this threshold is negative. The extinction of both populations is due to the small input nutrient concentration and/or the large death and washout rates of the phytoplankton. ¯ P, ¯ 0), where N¯ and P¯ satisfy (2.5) and (2.6), respectively. Let R0 > 0. Then (2.1) has another trivial periodic solution E1 = (N, The stability of E1 is determined by the fundamental matrix 1 (t) to the following linear τ -periodic system,

 x y z



x = J(E1 ) y , z

where

(2.13)

⎛

−D − ah (N¯ )P¯ J(E1 ) = ⎝ abh (N¯ )P¯ 0

−ah(N¯ ) 0 0

⎞

0 ⎠. −α f (P¯ ) ¯ ¯ αβ f (P) − μ − D2 −  (t )g(P)

(2.14)

Letting z = 0, we can obtain two linearly independent solutions using the fundamental matrix for the subsystem (2.3), namely ¯ ¯ eJ1 (N,P)t . To find the third linearly independent solution, we let

z˜(t ) = e

t 0

(αβ f (P¯)−μ−D2 −(s)g(P¯))ds .

(2.15)

Then (2.13) reduces to the linear nonhomogeneous system

  x y

  x y

¯ P¯ ) = J1 (N,





0 , −α f (P¯ )z˜(t )

+

which can be solved easily with



u(t ) =

u1 (t ) u2 (t )



=e

¯ P¯ )t J1 (N,





t

e 0

¯ P¯ )s −J1 (N,



0 ds. −α f (P¯ )z˜(s)

We have the following fundamental matrix for (2.13)



1 (t ) =

eJ1 (N,P)t 0 ¯ ¯



u(t ) , z˜(t )

¯ P¯ ) is always locally asymptotically stable for the subsystem (2.3), E1 is asymptotically where 0 is the 1 × 2 zero vector. Since (N, ¯ P, ¯ 0) is asymptotically stable if αβ f (P¯ ) − μ − D2 − stable for (2.1) if z˜(τ ) < 1 and unstable if z˜(τ ) > 1. It follows that E1 = (N, ˆ < 0 and unstable if αβ f (P¯ ) − μ − D2 − g(P¯ ) ˆ > 0. Denote g(P¯ )

∗ = {(N, P, Z ) ∈ R3+ : P > 0}.

(2.16)

¯ P, ¯ 0) is asymptotically stable if αβ f (P¯ ) − μ − D2 − Proposition 2.4. Let R0 > 0. Then E0 = (N0 , 0, 0) is unstable and E1 = (N, ˆ < 0 and unstable if αβ f (P¯ ) − μ − D2 − g(P¯ ) ˆ > 0. If αβ − μ − D2 ≤ 0, then the trivial periodic solution E1 is globally asympg(P¯ ) totically stable in ∗ . Proof. We only need to prove the last statement. If αβ < μ + D2 , then Z  < (αβ − μ − D2 )Z for t ≥ 0 implies limt→∞ Z (t ) = 0. If αβ − μ + D2 = 0, then Z  < (αβ − μ − D2 )Z implies limt→∞ Z (t ) = z∗ ≥ 0 exists. If z∗ > 0, then we must have limt→∞ Z  (t ) = 0, which is impossible. This shows limt→∞ Z (t ) = 0 and (2.1) is asymptotically autonomous to system (2.3). Furthermore, E1 is

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

57

ˆ g(P¯ ) − D2 < αβ − μ − D2 < 0 by the assumption. Therefore, E1 is globally locally asymptotically stable since αβ f (P¯ ) − μ −  asymptotically stable in ∗ by Proposition 2.2(b) and [24].  Applying LaSalle’s invariance principle for nonautonomous systems [25], we have the following results on the global stability of E1 . Theorem 2.5. Let R0 > 0 and αβ − μ − D2 > 0. (a) If f (P ) = (b) If f (P ) =

P ¯ P, ¯ 0) is globally asymptotically stable in ∗ if , then E1 = (N, k+P 2 P ¯ P, ¯ 0) is globally asymptotically stable in ∗ , then E1 = (N, k2 +P 2

αβ P¯ ≤ (μ + D2 )k. if αβ P¯ ≤ 2k(μ + D2 ).

Proof. By the assumption, E0 is unstable and E1 exists. Let (N(t), P(t), Z(t)) be an arbitrary solution of (2.1) with (N(0), P(0), Z(0)) ∈ ∗ . Notice that if N(0) = 0 then N(t) > 0 for t > 0 since N (0) > 0. Therefore to show E1 is globally asymptotically stable in ∗ it is enough to consider the subset = {(N, P, Z ) ∈ R3+ : N > 0, P > 0}. On , we define

V (N, P, Z ) =



N N¯

abh(x) − δ − D1 dx + ah(x)



P P¯

x − P¯ dx + CZ, x

(2.17)

where C > 0 will be determined below. Then V(N, P, Z) ≥ 0 and V (N, P, Z ) = 0 if and only if (N, P, Z ) = E1 . Since V does not depend on t explicitly, the time derivative of V along with trajectories of (2.1) is given by

abh(N) − δ − D1 P − P¯ N˙ + P˙ + C Z˙ P ah(N) P − P¯ P − P¯ abh(N) − abh(N¯ ) (abh(N) − abh(N¯ ))P − α f (P)Z = (D(N0 − N) − ah(N)P) + P P ah(N)

V˙ =

+ C (αβ f (P ) − μ −  (t )g(P ) − D2 )Z



= ab(h(N) − h(N¯ ))

D(N0 − N) − P¯ ah(N)



− α f (P )Z +

f (P ) ¯ + C αβ f (P )Z − C (μ +  (t )g(P ) + D2 )Z. α PZ P

(2.18)

N0 −N) Notice that h(N) is increasing and D(ah (N) is decreasing. Hence

  D(N0 − N) ¯ (h(N) − h(N¯ )) − P¯ < 0 for all 0 < N = N, ah(N)

and we need to determine C so that the remaining terms in V˙ become non-positive. (a) If f (P ) =

P , k+P

then

f (P ) P

=

1 k+P

≤

1 k

ˆ g(P¯ ) = αβ P¯ − μ − D2 −  ˆ g(P¯ ) ≤ αβ P¯ − μ − D2 −  ˆ g(P¯ ) ≤ and αβ f (P¯ ) − μ − D2 −  ¯ k k+P

P −m g(P¯ ) < 0 so that E1 is locally asymptotically stable. We can, from (2.18), choose C > 0 so that C ≤ 1/β and C ≥ k(μα+D 2) ¯ ¯ ˙ ˙ since by the assumption αβ P ≤ k(μ + D2 ) holds. Then V ≤ 0 and the set V = 0 on , the closure of , is M = {(N, P, Z ) ∈ ¯ Z = 0}. R3+ : N = N, 2 αβ P¯2 αβ P¯ P 1 ˆ ¯ ˆ ¯ ¯ (b) If f (P ) = k2P+P2 , then f (PP) = k2 +P 2 ≤ 2k and αβ f (P ) − μ − D2 −  g(P ) = 2 ¯ 2 − μ − D2 −  g(P ) ≤ 2k − μ − D2 − ¯

k +P

ˆ g(P¯) ≤ −m g(P¯) < 0 implying that E1 is also asymptotically stable. In (2.18) we choose C > 0 so that C ≤ 1/β and P¯ ¯ . Notice that such a C exists by the assumption given in (b). It follows that V˙ ≤ 0 and the set V˙ = 0 on C ≥ 2k(μα+D 2) 3 ¯ is M = {(N, P, Z ) ∈ R+ : N = N, Z = 0}.

¯ is M = {(N, P, Z ) ∈ R3+ : N = N, ¯ As a result, V is a Liapunov function on for cases (a)–(b). Since the set V˙ = 0 on Z = 0} and the only invariant set in M is {E1 } by Proposition 2.2(b), LaSalle’s invariance principle [23] implies that E1 is globally attracting in and hence it is globally asymptotically stable in .  Notice that in (a) and (b) f is of Holling types II and III, respectively. These functional forms are frequently adopted in modeling plankton dynamics. However, if f is neither of this form, we can apply a similar technique to derive a sufficient condition. ¯ In particular, the conditions are satisfied if P¯ is small. Furthermore, the sufficient conditions derived in Theorem 2.5 depend on P. Therefore, E1 is expected to be globally asymptotically stable when it just appears, that is, as soon as P¯ becomes positive provided that (2.7) is satisfied. We now turn to investigate interior steady states of (2.1) when A = 0. Let R0 > 0. It is clear that there is no interior steady state if αβ < (μ + D2 ) by Proposition 2.4. Define

ˆ R1 = αβ f (P¯ ) − μ − D2 − g(P¯ )

(2.19)

ˆ = γ since A = 0. The P-component of an interior steady state of (2.1) with A = 0 satisfies and assume R1 > 0. Notice 

αβ f (P) − μ − γ g(P) − D2 = 0.

(2.20)

58

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67 Table 2 Different combinations of functional forms for f and g. Cases 1 2 3 4

f(P)

g(P)

P k+P P2 k2 + P 2 P2 k2 + P 2 P k+P

P m+P P m+P P2 m2 + P 2 P2 m2 + P 2

Holling types Type II, Type II Type III, Type II Type III, Type III Type II, Type III

Let F(P) be the left hand side of (2.20). Then F (0) = −μ − D2 < 0 and F (P¯ ) > 0 since R1 > 0. Therefore, F (P ) = 0 has at least one ¯ The N-component of an interior steady state then solves positive solution P∗ , where P ∗ < P.

D(N0 − N) = ah(N)P ∗ . Clearly (2.21) always has a positive solution N∗ of N∗ and P∗ as

Z∗ =

(2.21) with N∗

< N0 . The Z-component of an interior steady state can be written in terms

(abh(N∗ ) − δ − D1 )P∗ , α f (P∗ )

(2.22)

¯ where F (P ∗ ) = 0 where Z∗ > 0 if abh(N∗ ) > δ + D1 . We conclude that an interior steady state E ∗ = (N∗ , P ∗ , Z ∗ ) exists if N∗ > N, ¯ (2.21) implies N∗ > N. ¯ Consequently, (2.1) has at least one interior steady state when A = 0 and N∗ satisfies (2.21). Since P ∗ < P, and R1 > 0. Although (2.1) with A = 0 has at least one interior steady state when R0 > 0 and R1 > 0, the number of interior steady states are not easy to determine analytically for general functional forms of f, g and h. We assume

αβ − ˆ − μ − D2 > 0

(2.23)

holds and consider the following cases for different combinations of Holling-type functions listed in Table 2. We prove that the interior steady state is unique for cases 1–3 if R0 > 0, R1 > 0 and (2.23) hold. F (P ) For case 1, F(P) can be written as F (P ) = (k+P1)(m+P) , where

F1 (P ) = (αβ − γ − μ − D2 )P 2 + (αβ m − (μ + D2 )(m + k) − γ k)P − (μ + D2 )km. Since αβ − γ − μ − D2 > 0 by (2.23), F1 (P ) = 0 has a unique positive solution P∗ with P ∗ < P¯ and hence (2.1) with A = 0 has a unique interior steady state. F (P ) For case 2, F(P) can be rewritten as F (P ) = (k2 +P22 )(m+P) , where

F2 (P ) = (αβ − γ − μ − D2 )P 3 + m(αβ − μ − D2 )P 2 − k2 (γ + μ + D2 )P − (μ + D2 )k2 m. Since there is only one sign change in the coefficients of F2 (P), F2 (P ) = 0 has a unique positive solution P∗ with P ∗ < P¯ and (2.1) with A = 0 therefore has a unique interior steady state. F (P ) For case 3, F(P) results in F (P ) = (k2 +P23)(m2 +P2 ) , where

F3 (P ) = (αβ − γ − μ − D2 )P 4 + (αβ m2 − (μ + D2 )(m2 + k2 ) − γ k2 )P 2 − (μ + D2 )k2 m2 . ¯ Hence (2.1) with A = 0 has a unique interior It is clear that F3 (P ) = 0 also has a unique positive solution which is less than P. steady state. F (P ) The above finding, however, is not conclusive for case 4. Indeed, F (P ) = 0 becomes (k+P)(4 m2 +P2 ) = 0, where

F4 (P ) = (αβ − γ − μ − D2 )P 3 − k(γ + μ + D2 )P 2 + m2 (αβ − μ − D2 )P − (μ + D2 )km2 . Since there are three sign changes in the coefficients of F4 , Descartes’ Rule of Signs [22] implies that F4 (P ) = 0 has at least one ¯ Consequently, (2.1) with positive solution that is less than P¯ but may have three positive solutions each of which is less than P. A = 0 may have three interior steady states. We now summarize. Proposition 2.6. If R0 > 0, R1 > 0 and (2.23) are satisfied, then system (2.1) with A = 0 has a unique interior steady state for cases 1–3 and has at least one interior steady state for case 4. Although it is not easy to determine the number of interior steady states analytically for case 4, one can verify that both plankton populations coexist when R0 > 0 and R1 > 0 under the general assumptions of (H1)–(H3). One of the important math= f˜(t, x), ematical tools in studying coexistence of populations is the concept of uniform persistence. An ecological system dx dt where f˜ : R+ × Rn → Rn , is said to be uniformly persistent if there exists η > 0 such that lim inft→∞ xi (t ) ≥ η for all solutions +

+

x(t ) = (x1 (t ), . . . , xn (t )) with xi (0) > 0 for 1 ≤ i ≤ n. Similar to a periodic nutrient-phytoplankton-zooplankton model studied in

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

59

[26], we can apply Theorem 3.1 of Butler and Waltman [27] to show uniform persistence of the periodic system (2.1) when the trivial τ -periodic solutions E0 and E1 are unstable. In such case, the system also has a positive τ -periodic solution by Yang and Freedman [28, Theorem 4.11]. To this end, we summarize some terminology needed for the analysis which can also be found in [27]. Let F be the flow generated by (2.1) and let ∂ F be F restricted to the boundary ∂(R3+ ). A nonempty invariant set M is called isolated invariant if M is the maximal invariant set in some neighborhood of itself. The neighborhood is then called an isolated neighborhood of M. Let W + (M) and W − (M) denote the stable manifold and the unstable manifold of M, respectively. Two isolated invariant sets M and    N are said to be chained if there exists x ∈ / M N such that either x ∈ W − (M) W + (N) or x ∈ W + (M) W − (N). For the former case, M is chained to N and we write M → N and for the latter case N is said to be chained to M and it is written as N → M. A finite sequence M1 , M2 , . . . , Mk , k > 1, of isolated invariant sets will be called a chain if M1 → M2 →  → Mk . The chain is called a cycle if Mk = M1 .  Let ω(x) denote the ω-limit set of x ∈ R3+ and (∂ F ) be the collection of all ω-limit sets of ∂(R3+ ), i.e., (∂ F ) = x∈∂(R3 ) ω(x). +

The boundary flow ∂ F is said to be isolated if there exists a covering M of (∂ F ) by pairwise disjoint, compact, isolated invariant sets M1 , M2 , . . . , Mk for ∂ F such that each Mi is also isolated for F. M is then called an isolated covering. The boundary flow ∂ F  is called acyclic if there exists some isolated covering M = ki=1 Mi of (∂ F ) such that no subset of the {Mi } forms a cycle. Applying Theorem 3.1 of Butler and Waltman [27], it can be shown that system (2.1) is uniformly persistent if R0 > 0 and R1 > 0 are satisfied. Theorem 2.7. If R0 > 0 and R1 > 0, then system (2.1) is uniformly persistent and (2.1) has a positive τ -periodic solution. Proof. Let F be the continuous flow generated by (2.1) and let ∂ F be F restricted on the boundary ∂(R3+ ). Let M1 = {E0 } and M2 = {E1 }. Note that F is point dissipative by Proposition 2.1 and (∂ F ) consists of only {E0 , E1 }. Then M = {M1 , M2 } is composed of disjoint, compact and isolated invariant sets for ∂ F. The stable manifold W0+ of E0 is two-dimensional and is the nonnegative NZ-plane. The unstable eigenvalue of J(E0 ) is R0 . A simple calculation shows that an eigenvector (x, y, z)T of J(E0 ) belonging to −D−R R0 can be chosen to be z = 0, x = 1 and y = ah(N00) < 0. Therefore the unstable manifold W0− of E0 lies outside of R3+ . Similarly, it

can be seen from J(E1 ) given by (2.14) that the stable manifold W1+ of E1 is two-dimensional and lies on the nonnegative NP-plane. As a result, M1 and M2 cannot be chained.  It remains to verify that each Mi is isolated for F, i = 1, 2. Let c0 = max[0,P] ¯ f (P ) > 0. Using R0 > 0, we can choose  > 0 so that

abh(N0 − ) − α c0  − δ − D1 > 0. We prove that the set N1 = {(N, P, Z ) ∈ R3+ : d((N, P, Z ), M1 ) < } is an isolated neighborhood of M1 in R3+ , where d is the Euclidean metric on R3 . Indeed, if N1 is not an isolated neighborhood of M1 in R3+ , then we can find an invariant set V1 such that M1 ⊂ V1 ⊂ N1 and V1 \M1 = ∅. Let (N(0), P(0), Z(0)) ∈ V1 with P(0), Z(0) > 0. Then (N(t), P(t), Z(t)) ∈ V1 for t > 0 and there exists ξ ∈ (0, P) ⊂ [0,  ] such that

P = abh(N) − α f  (ξ )Z − δ − D1 > abh(N0 − ) − α c0  − δ − D1 > 0 P for t > 0. It follows that P(t) → ∞ as t → ∞ and we obtain a contradiction by Proposition 2.1. Therefore, M1 is isolated for F. Similarly, R1 > 0 implies that we can find  0 > 0 sufficiently small such that

αβ f (P¯ − 0 ) − μ − D2 − ˆ g(P¯ + 0 ) > 0. We prove that the set N2 = {(N, P, Z ) ∈ R3+ : d((N, P, Z ), M2 ) < 0 } is an isolated neighborhood of M2 in R3+ . If not, then we can find an invariant set V2 with M2 ⊂ V2 ⊂ N2 and V2 \M2 = ∅. Let (N(0), P(0), Z(0)) ∈ V2 with Z(0) > 0. Then (N(t), P(t), Z(t)) ∈ V2 for t > 0 and in particular, P¯ − 0 < P (t ) < P¯ + 0 for t > 0. Consequently, 

t ¯ ¯ Z (t ) > Z (0)e 0 (αβ f (P−0 )−μ− (s)g(P+0 )−D2 )ds → ∞

as t → ∞ and we obtain another contradiction. This shows that M2 is also isolated for F and ∂ F is acyclic. Since Wi+ ∩ int (R3+ ) = ∅ for i = 0, 1, (2.1) is uniformly persistent by [27, Theorem 3.1]. Moreover, the system has a positive τ -periodic solution by [28, Theorem 4.11] since (2.1) is point dissipative and uniformly persistent.  The threshold R1 can be interpreted as the maximal net growth rate of the zooplankton population when the phytoplankton population is stabilized at the P¯ level. The zooplankton population can persist if this threshold is positive. Therefore, the ˆ of the phytoplankton is small. zooplankton population is more likely to survive if the average toxin production rate  2.2. Examples and numerical simulations Although both plankton populations can coexist under the assumptions of R0 > 0 and R1 > 0, their interactions may be complicated. In this subsection, we consider specific examples with numerical simulations to study the models in more detail. We simulate a particular model using Holling type II for f, g and h with  (t ) = γ (1 + A sin (Tt )):

60

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67 Table 3 Default parameter values1 . Parameters

Default values

Units

a D D1 D2 N0 m1

0.6 0.2 0.03 0.01 6 0.1 0.5 0.01 0.01 10 10 0.5 0.6

day−1 day−1 day−1 day−1 μgNl −1 μgNl −1 day−1 day−1 day−1 μgNl −1 μgNl −1

α δ μ

m k b

β

1 These parameter values are within the ranges given in reference [29].

⎧ aNP ⎪ ⎪ N = D(N0 − N) − ⎪ ⎪ m 1 +N ⎪ ⎨ abNP α PZ P = − − δ P − D1 P m1 + N k+P ⎪ ⎪ ⎪ ⎪ ⎪ ⎩Z  = αβ PZ − μZ − γ (1 + A sin (Tt )) PZ − D2 Z. k+P

(2.24)

m+P

The time unit is a day and we adopt the following default parameter values given in Table 3, where b and β have no units since they are fractions. These parameter values are within the ranges given in [29], where the data were collected from different 2π so that the toxin liberation by phytoplankton research articles using Holling Type II functional forms for f and h. Let T = 365 ˆ has an annual periodicity, τ = 365. Moreover,  = γ , m = γ (1 − A) and M = γ (1 + A). We will vary γ (day−1 ) and A in the following simulations. When the phytoplankton population is not toxin producing, γ = 0, then the system has a positive periodic solution and the two plankton populations are oscillating over time under the default parameter values given in Table 3. In the following simulations we shall assume R0 > 0 and R1 > 0 so that the trivial periodic solutions E0 and E1 are unstable. Since system (2.24) is of case 1 listed in Table 2, there exists a unique interior steady state E∗ when A = 0 and its stability depends on N0 and D. Fig. 1 plots the PZ components of trajectories with γ = 0.01 and different values of A after the transient behavior is discarded. We see that increasing amplitude A of the toxin liberation will make the interactions more complicated including chaotic behavior. Chaotic interactions are also observed in a nutrient-phytoplankton model with periodic releasing of toxin by phytoplankton [20]. In this parameter regime, decreasing nutrient inflow rate D can further destabilize the planktonic interactions. See Fig. 1(d). However, if the input nutrient concentration N0 is decreased to a smaller level, then it can make the interaction more stable even with a smaller nutrient inflow rate D as illustrated in Fig. 1(e) and (f). If the maximum toxin liberation rate γ is increased to γ = 0.02 and with a large input nutrient concentration N0 , then the earlier observation remains hold. That is, the planktonic interaction is destabilized if the amplitude of the toxin liberation is increased. However, since the maximum toxin production rate is larger, decreasing nutrient inflow rate D will not further destabilize the nutrient-plankton system as in the case of γ = 0.01. See Fig. 2. We perform similar numerical investigations with γ = 0.03, 0.05, 0.1, and 0.25 and conclude that the nutrient-plankton dynamics behave similarly as in the case of γ = 0.02. If the maximum toxin liberation rate γ is increased above 0.25, then R1 > 0 is no longer satisfied and the zooplankton population becomes extinct. We provide a particular simulation result for N0 = 6, D = 0.2, γ = 0.5 and A = 0.5 in Fig. 3(a), where the phytoplankton population is stabilized in about 120 days. Using the same parameter values of N 0 = 6 and D = 0.2, we numerically approximate the maximum Liapunov exponents [30] with γ = 0.01 in Fig. 3(b), where positive Liapunov exponent is an indication of chaos. 3. The stochastic model and analysis Natural populations can have large variations among individuals in patterns of survival and reproduction over time which can be associated with differences in state variables such as age, size or spatial location. However, a majority of the variations are random due to environmental fluctuations. This is particularly true for the plankton populations because of unpredictability of weather, temperature and many other physical factors embedded in the water environment. The deterministic model studied in the previous section does not take environmental randomness into consideration. The evolution of the population interactions is determined precisely by the differential equations. In this section we will investigate the role of unpredictable environmental factors in plankton dynamics. Specifically, we seek to understand whether the biological conclusions obtained from the deterministic systems remain valid if the environment is randomly varied.

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

61

Fig. 1. Phase plane trajectories of (2.24) are plotted with γ = 0.01 after the transient behavior is discarded.

In the stochastic model, it is hypothesized that the environment is randomly varying, in particular, the input nutrient concentration N0 is varying randomly with respect to time. Continuous, mean-reverting processes have been proposed as being biologically reasonable models for parameters undergoing environmental variability [31]. In the present investigation, it is assumed that the input nutrient concentration is given by a continuous, mean-reverting, stochastic process that satisfies an Itô stochastic differential equation of the form

dN0 (t ) = αN (Ne0 − N0 (t )) dt + βN (N0 (t ))ν dW (t )

(3.1)

where αN , βN , Ne0 and ν are positive constants and W(t) is a standard Wiener process. For 1/2 ≤ ν ≤ 1, it is known that the process N0 (t) that satisfies (3.1) is continuous and nonnegative with probability unity [32]. With this biologically-reasonable stochastic model (3.1), the input concentration N0 (t) varies randomly and continuously about an asymptotic mean value Ne0 with a variance, for example, of

βN2 Ne0 2αN when ν = 1/2.

Another possible reasonable form for N0 (t) is given by a continuous, mean-reverting, stochastic process that satisfies an Itô stochastic differential equation (SDE) of the form

d ln (N0 (t )) = αN ( ln (Ne0 ) − ln (N0 (t ))) dt + βN dW (t ).

(3.2)

In this case, N0 (t) can be shown to satisfy a log-normal distribution for any time t. Here, also, N0 (t) is continuous and nonnegative with probability unity. In the present investigation, however, N0 (t) is assumed to satisfy (3.1).

62

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

Fig. 2. Phase plane trajectories of (2.24) are plotted with γ = 0.02 after the transient behavior is discarded.

With this change, system (2.1) becomes the SDE system:

⎧ 0 ⎪ ⎪dN = (D(N (t ) − N) − ah(N)P) dt ⎨ dP = (abh(N)P − α f (P )Z − δ P − D1 P ) dt dZ = (αβ f (P )Z − μZ −  (t )g(P )Z − D2 Z ) dt ⎪ ⎪ ⎩ 0 dN = αN (Ne0 − N0 ) dt + βN (N0 )ν dW (t ).

(3.3)

As 1/2 ≤ ν ≤ 1, there is a nonnegative function N0 (t) that satisfies the fourth equation in (3.3) for each Wiener process W(t). In addition, as dE (N0 (t )) = αN (Ne0 − E (N0 (t )))dt, then E (N0 (t )) = Ne0 + (N0 (0) − Ne0 )e−αN t → Ne0 as t → ∞. In particular, E(N0 (t)) is 0 = max{N 0 (0), N 0 } for all t > 0. There also exist nonnegative solutions N(t), P(t) and Z(t) to (3.3). As in section 2, let bounded by Nm e F (t, N, P, Z ) = (F1 , F2 , F3 ) where F1 = (D(N0 (t ) − N) − ah(N)P ), F2 = (abh(N)P − α f (P )Z − δ P − D1 P ), and F3 = (αβ f (P )Z − μZ −  (t )g(P)Z − D2 Z ). Then, we have the following proposition. Proposition 3.1. Solutions of (3.3) exist for t > 0 and remain nonnegative with probability one. Proof. As F is locally Lipschitz with F1 |N=0 ≥ DN0 (t ) ≥ 0, F2 |P=0 = 0 and F3 |Z=0 = 0, it follows from Theorem A.4 of [21] that (3.3) has a unique nonnegative solution defined on some interval [0, tˆ). To show that tˆ = ∞, let S = bβ N + β P + Z. Then,  0 x(0) = S(0). S ≤ bβ DN0 (t ) − mS, where  m = min{D, δ +D1 , μ + D2 }. Consider the scalar linear equation x = bβ DNm − mx with 0 /m + x(0) − bβ DN 0 /m e−mt and E(S(t)) ≤ x(t) for t ≥ 0, we have lim sup 0 /m. Given  > E ( S ( t )) ≤ b β DN Since x(t ) = bβ DNm t→∞ m m 0 /(m). By the Chebyshev–Markov inequality, 0, let M = bβ DNm

lim sup Prob(S(t ) ≤ M) = 1 − lim sup Prob(S(t ) ≥ M) t→∞

t→∞

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

63

Fig. 3. N0 = 6.0 and D = 0.2 in these simulations for model (2.24). (a) provides plankton populations over time when γ = 0.5 and A = 0.5. (b) approximates the Liapunov exponents with γ = 0.01 and A is varied from 0 to 0.99.

≥ 1 − lim sup E (S(t ))/M ≥ 1 −  . t→∞

As  > 0 is arbitrary, this implies that S(t) is bounded with probability unity for all t > 0. Therefore, solutions of (3.3) are defined and nonnegative for all t > 0 with probability unity.  Consider now the stochastic generalization of numerical example (2.24) which has the form:

⎧ dN = (D(N0 (t ) − N) − ah(N)P ) dt ⎪ ⎨ dP = (abh(N)P − α f (P )Z − δ P − D1 P ) dt ⎪ ⎩dZ 0= (αβ f (P0 )Z −0μZ −  (t )g(0P)νZ − D2 Z ) dt dN = αN (Ne − N ) dt + βN (N ) dW (t ),

(3.4)

where h(N) = N/(m1 + N), f (P ) = P/(k + P ), g(P ) = P/(m + P ), and  (t ) = γ (1 + A sin (Tt )). The parameter values for this computational example are selected as a = 0.6, D = 0.2, D1 = 0.03, D2 = 0.01, Ne0 = 6, m1 = 0.1, α = 0.5, δ = 0.01, μ = 0.01, m = 10, k = 10, b = 0.5, β = 0.6, T = 2π /365, γ = 0.01, A = 0.5, αN = 0.5, βN = 1.0, ν = 0.50. The Euler–Maruyama method [33] is applied to compute sample paths of (3.4). This problem’s deterministic solution of zooplankton versus phytoplankton is graphed in Fig. 1(b). One sample path of the solution of (3.4) is graphed in Fig. 4(a) along with the mean zooplankton population size versus mean phytoplankton size for 2000 sample paths up to final time of 2500 in Fig. 4(b). The results for the stochastic case are interesting as one sample path can vary widely in population levels as compared with the deterministic case but the mean plankton levels appear to approach a stable equilibrium, in particular when A = 0 and so  (t ) = γ . To study this, consider the mean values of the levels for (3.4). The mean levels satisfy the deterministic differential equation system:

⎧ dE (N(t )) = (D(E (N0 (t )) − E (N(t ))) − aE (h(N)P )) dt ⎪ ⎨ dE (P (t )) = (abE (h(N)P ) − α E ( f (P )Z ) − δ E (P (t )) − D1 E (P (t ))) dt dE (Z (t )) = (αβ E ( f (P )Z ) − μE (Z (t )) − E ( (t )g(P )Z ) − D2 E (Z (t ))) dt ⎪ ⎩ dE (N0 (t )) = αN (Ne0 − E (N0 (t ))) dt.

(3.5)

Consider now whether a stationary distribution is reached as t → ∞ and so whether the moments approach constant values. It is known that N0 (t) approaches a stationary distribution [32]. We need to consider N, P, and Z. If a stationary distribution is approached, the moments of the solution must approach constant values as t → ∞. Assume that E(N(t)) → N∗ , E(P(t)) → P∗ , E(Z(t)) → Z∗ as t → ∞. As E(N0 (t)) approaches Ne0 for large time t, system (3.4) can be written as:

(t )) d(y (y (t )) + c(t ) =R dt

(3.6)

(t ) = [y1 (t ), y2 (t ), y3 (t )]T , y1 (t ) = N(t ), y2 (t ) = P (t ), and y3 (t ) = Z (t ). Also, c1 (t ) = DN0 (t ), c2 (t ) = c3 (t ) = 0, and where y (t )) → y ∗ . Let limt→∞ E (c(t )) = c∗ . Consider the system E (c1 (t )) → DNe0 = c1∗ . Suppose that E (y

d y1 (y 1 (0) = y ∗ , c(0) = c∗ 1 ) − c(t ) with y =R dt d y2 (y 2 (0) = y ∗ + 0 , c(0) = c∗ 2 ) − c(t ) with y =R dt

(3.7)

2 (t ) − y 1 (t ). Then,  0 is small. Let   (t ) = y  (t ) satisfies where 

d   (y  (y 2 ) − R 1 ) ≈ R ∗ ) with   (0) =  0 = R(y dt

(3.8)

64

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

Fig. 4. Ne0 = 6.0 and γ = 0.01 in these simulations for model (3.4). (a) and (b) provide one sample path and the mean path for plankton population levels with A = 0.5. (c) and (d) are the mean paths for plankton population levels with A = 0.0 and (e) and (f) with A = 0.9.

 (y ∗ ) is the Jacobian matrix evaluated y ∗ . Solving (3.8) for E (  (t )) gives where R

 (y 2 (t )) − E (y 1 (t )) ≈ exp (R ∗ )t )0 . E ( (t )) = E (y

(3.9)

 (y ∗ ) has positive Eq. (3.9) implies that no stable stationary distribution exists if any of the eigenvalues of the Jacobian matrix R real part.  (y ∗ ) has the form For Eq. (3.6), the Jacobian R

⎛

−D − ahN (N∗ )P ∗ ⎝ abhN (N∗ )P∗ 0

ah(N∗ ) abh(N∗ ) − α fP (P ∗ )Z ∗ − δ − D1 αβ fP (P∗ )Z ∗ − γ gP (P∗ )Z ∗

⎞

0 −α f (P ∗ )⎠, j33

(3.10)

 (y ∗ ) were where j33 = αβ f (P ∗ ) − μ − D2 − γ g(P ∗ )Z ∗ . For the numerical example above where A = 0, the eigenvalues of R ∗ calculated as −0.1994, 0.0103 + 0.0692i, 0.0103 − 0.0692i for possible approximate equilibrium mean values of (N , P ∗ , Z ∗ ) = (3.845, 0.7772, 5.325). However, the eigenvalues indicate that a stationary distribution is not reached. Indeed, calculated mean plankton population levels, illustrated in Fig. 4(c) and (d) for 20,000 sample paths up to time 3500, have an oscillatory behavior with oscillations of small magnitude at large time t. Fig. 4(e) and (f) provide mean population levels with γ = 0.01 and A = 0.9. We also simulate the model with Ne0 = 3.0 in Fig. 5. We observe similar behavior as the deterministic model, namely a decreasing of the input nutrient concentration Ne0 reduces the planktonic oscillations. Furthermore, if γ is increased to 0.27, then the mean zooplankton population also approaches zero when Ne0 = 3.0 as presented in Fig. 6. Therefore, increasing toxin production rate can eliminate zooplankton population even in the stochastic environment.

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

65

Fig. 5. Ne0 = 3.0 and γ = 0.01 in these simulations for model (3.4). (a) and (b) provide one sample path and the mean path for plankton population levels with A = 0.5. (c) and (d) are the mean paths for plankton population levels with A = 0.0 and (e) and (f) with A = 0.9.

4. Summary and conclusions Phytoplankton form the base of every food web and they generate about half of the atmosphere’s oxygen. Phytoplankton are consumed by the zooplankton, which transfer energy and organic matter to higher trophic levels. Through their consumption and processing of phytoplankton and other food sources, zooplankton play an important role in aquatic food webs as a resource for consumers of higher trophic levels. There has been a global increase in harmful algal blooms in the last few decades [1] and toxin producing phytoplankton (TPP) are among the contributors in these blooms. The TPP species have shown to have extremely harmful effects on the environment and on other species interacting with them [6,19,34,35]. There are many deterministic models of plankton interactions proposed to study the effects of toxin producing phytoplankton on phytoplankton blooms [3,11–14,20]. However, plankton populations reside in any water body persistently undergo environmental fluctuations. It is unclear whether the biological conclusions made from these deterministic settings remain true in the stochastic systems. Therefore, it is important to investigate the stochastic version of the deterministic models in order to make more precise biological findings. In this study, we propose deterministic and stochastic models of nutrient-phytoplankton-zooplankton interaction to study the effects of toxin producing phytoplankton upon planktonic dynamics. The toxin producing rate is assumed to vary periodically in the deterministic and stochastic systems to model seasonality of the toxin liberation. There are several two trophic levels of either nutrient-phytoplankton or phytoplankton-zooplankton models incorporating TPP in the literature [3,11,12,14,20,36]. The only model of three-trophic-level interactions with TPP consideration is studied in [13], where the toxin liberation is not explicitly modeled.

66

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

Fig. 6. Ne0 = 3.0 and γ = 0.27 in these simulations for model (3.4) with A = 0, A = 0.5, and A = 0.9 in (a) and (b), (c) and (d) and (e) and (f), respectively. These simulations demonstrate that the zooplankton population also goes extinct even in the stochastic environment if the toxin liberation rate by phytoplankton is large.

Based on the model parameters and the nutrient uptake and zooplankton grazing functions, we derive two lumped thresholds ˆ of the R0 and R1 . The threshold R0 depends on the constant input nutrient concentration N0 while R1 depends on the average  toxin liberation rate  (t). It is proven in Theorem 2.7 that the plankton populations in the deterministic setting can coexist if these two thresholds are positive. This analytical result is new in comparing with the study given in [20] for which periodic toxin liberation is also investigated. Both plankton populations are more likely to persist if either the nutrient input concentration N0 ˆ of the toxin liberation rate is large or/and if the phytoplankton outflow rate D1 is small. On the other hand, if the mean value  is large so that the threshold R1 becomes negative, then the zooplankton population will become extinct and consequently the phytoplankton population will be stabilized. Therefore, toxin liberation by phytoplankton can terminate harmful algal blooms by increasing zooplankton mortality. This finding is consistent with the results obtained in [3,11–14,36]. We also develop a stochastic model with Itô differential equations to model environmental variability. Specifically, the input nutrient concentration in the deterministic system is allowed to varying randomly with respect to time. In particular, the input nutrient concentration is modeled with a continuous, mean-reverting process which is biologically reasonable for the nutrient concentration undergoing environmental variability [31]. This stochastic modeling of the environmental randomness is new in the nutrient-plankton investigation. It is shown that the populations cannot reach a stationary distribution even when the toxin liberation is a constant. We simulate the deterministic system with parameter values taken from the literature by varying the input nutrient concentration along with the toxin producing parameters. The stochastic model is also numerically simulated using the same parameter values as the deterministic system. The simulations indicate that increasing the magnitude of the toxin production rate can induce chaotic interactions for the deterministic model (cf. Fig. 1(c) and (d), Fig. 2(b) and (c)). This chaotic behavior is also

S.R.-J. Jang, E.J. Allen / Applied Mathematics and Computation 271 (2015) 52–67

67

demonstrated in a nutrient-phytoplankton study [20] where the toxin production rate is varied periodically. The constant input nutrient concentration N0 , however, can diminish this chaotic interaction if it is smaller (cf. Fig. 1(e) and (f), Fig. 2(e) and (f)). In particular, the planktonic interaction is more stable if N0 is small independent of whether the toxin production rate is constant or varies with large amplitude. Surprisingly, this observation also holds even when the input nutrient concentration randomly varies with respect to time. If the input nutrient concentration N0 is small, then the interaction also depends on the maximum toxin liberation rate. As this toxin liberation rate increases, both the deterministic and the stochastic systems become more stable. Therefore, toxin producing phytoplankton can suppress zooplankton by increasing its mortality and can also eliminate phytoplankton blooms. Acknowledgments We thank the reviewers and the associate editor for providing helpful comments that improved the original manuscript. References [1] G.M. Hallegraeff, A review of harmful algae blooms and the apparent global increase, Phycologia 32 (1993) 79–99. [2] E.J. Buskey, D.A. Stockwell, Effects of a persistent brown tide on zooplankton population in the Laguna Madre of southern Texas, in: Smayda (Ed.), Toxic Phytoplankton Blooms in the Sea, Netherlands, 1993, pp. 659–666. [3] J. Chattopadhayay, R. Sarkar, S. Mandal, Toxin-producing plankton may act as a biological control for planktonic blooms-field study and mathematical modelling, J. Theor. Biol. 215 (2002) 333–344. [4] H. Macintyre, A. Stutes, W. Smith, C. Dorsey, A. Abraham, R. Dockey, Environmental correlates community ompositiion and toxicity during a bloom of Pseudo-nitzschia spp. in the northern Gulf of Mexico, J. Plankton Res. 33 (2011) 273–295. [5] D. Mcgillicuddy, R. Signell, C. Stock, B. Keafer, M. Keller, R. Hetland, D. Anderson, A mechanism for offshore initiation of harmful algal blooms in the coastal Gulf of Maine, J. Plankton Res. 25 (2003) 1131–1138. [6] D.L. Roelke, R.H. Pierce, Effects of inflow harmful algal blooms: some considerations, J. Plankton Res. 33 (2011) 205–209. [7] S. Uye, Impact of copepod grazing on the red tide flagellate Chattonella antiqua., Mar. Biol. 92 (1986) 35–41. [8] S. Ruan, Persistence and coexistence in zooplankton-phytoplankton-nutrient models with instantaneous nutrient recycling, J. Math. Biol. 31 (1993) 633–654. [9] S. Ruan, Oscillations in plankton models with nutrient recycling, J. Theo. Biol. 208 (2001) 15–26. [10] J. Steele, E. Henderson, A simple plankton model, Am. Nat. 117 (1981) 676–691. [11] J. Chattopadhayay, R. Sarkar, A delay differential equation model on harmful algal blooms in the presence of toxic substances, IMA J. Math. Appl. Med. Biol. 19 (2002) 137–161. [12] T. Saha, M. Bandyopadhyay, Dynamical analysis of toxin producing phytoplankton-zooplankton interactions, Nonlin. Anal.: RWA 10 (2009) 314–332. [13] B. Mukhopadhyay, R. Bhattacharyya, Modelling phytoplankton allelopathy in a nutrien-plankton model with spatial heterogeneity, Ecol. Model. 198 (2006) 163–173. [14] S.R.-J. Jang, J. Baglama, L. Wu, Dynamics of phytoplankton-zooplankton systems with toxin producing phytoplankton, Appl. Math. Comput. 227 (2014) 717–740. [15] D.L. DeAngelis, Dynamics of Nutrient Cycling and Food Webs, Chapman and Hall, London, 1992. [16] B.H. Ketchum, The flushing of tidal estuaries, Sew. Ind. Wastes 23 (1951) 198–209. [17] B.H. Ketchum, The relation between circulation and planktonic populations in estuaries, Ecology 35 (1954) 191–200. [18] C.J. Miller, et al., The role of inflow magnitude and frequency on plankton communities from the Guadalupe Estuary, Texas, Est. Coast Shelf Sci. 80 (2008) 67–73. [19] E. Philips, S. Badylak, S. Youn, K. Kelley, The occurrence of potentially toxic dinoflagellates and diatoms in a subtropical lagoon, the Indian River Lagoon, Florida, USA, Harmful Alg. 3 (2004) 39–49. [20] S. Charkraborty, S. Chatterjee, E. Venturino, J. Chattopadhayay, Recurring plankton bloom dynamics modeled via toxin-producing phytoplankton, J. Biol. Phys. 33 (2007) 271–290. [21] H.R. Thieme, Mathematics in Population Biology, Princeton University Press, New Jersey, 2003. [22] L.J.S. Allen, An Introduction to Mathematical Biology, Prentice-Hall, New Jersey, 2006. [23] W. Walter, Ordinary Differential Equations, Graduate Texts in Mathematics, Springer, New York, 1998. [24] H.R. Thieme, Convergence results and a Poincaré–Bendixson trichotomy for asymptotically autonomous differential equations, J. Math. Biol. 30 (1992) 755–763. [25] J.P. LaSalle, An invariance principle in the theory of stability, in: P.R. Mayaguez (Ed.), Proceedings of the International Symposium on Differential Equations and Dynamical Systems, Academic Press, New York, 1965, pp. 277–286. [26] S.R.-J. Jang, J. Baglama, Persistence in variable-yield nutrient-plankton models with nutrient recycling, Math. Comput. Model. 38 (2003) 281–298. [27] G.J. Butler, P. Waltman, Persistence in dynamical systems, J. Diff. Equ. 63 (1986) 255–263. [28] F. Yang, H.I. Freedman, Competing predators for a prey in a chemostat model with periodic nutrient input, J. Math. Biol. 29 (1991) 715–732. [29] A.M. Edwards, J. Brindley, Zooplankton mortality and the dynamical behaviour of plankton population models, Bull. Math. Biol. 61 (1999) 303–339. [30] T.S. Parker, L.O. Chua, Numerical Algorithms for Chaotic Systems, Springer, Berlin, 1989. [31] E.J. Allen, L.J.S. Allen, H. Schurz, A comparison of persistence-time estimation for discrete and continuous stochastic population models that include demographic and environmental variability, Math. Biosci. 196 (2005) 14–38. [32] A. Alfonsi, Strong convergence of some drift implicit Euler scheme, Statistics & Probability Letters 83 (2) (2013) 602–607. [33] P.E. Kloeden, E. Platen, H. Schurz, Numerical Solution of SDE Through Computer Experiments, Springer, Berlin, 1994. [34] B. Hubbart, G. Pitcher, B. Krock, A. Cembella, Toxigenic phytoplankton and concomitant toxicity in the mussel Choromytilus meridionalis off the west coast of South Africa, Harmful Alg. 20 (2012) 30–41. [35] S. Sopanen, O. Setälä, J. Piiparinen, K. Erler, A. Kremp, The toxic dinoflagellate Alexandrium ostenfeldii promotes incapacitation of the calanoid copepods Eurytemora affinis and Acartia bifilosa from the northern Baltic Sea, J. Plankton Res. 33 (2011) 1564–1573. [36] S. Kharea, O. Misraa, J. Dharb, Role of toxin producing phytoplankton on a plankton ecosystem, Nonlin. Anal. Hybrid Sys. 4 (2010) 496–502.