Dynamic Behaviors of a Stochastic Eco-Epidemiological Model for Viral Infection in the Toxin-Producing Phytoplankton and Zooplankton System

: It is well known that the evolution of natural populations is almost inevitably disturbed by various environmental factors. Various experiments have shown that the growth of phytoplankton might be affected by nutrient availability, water temperature, and light, while the development of zooplankton is more disturbed by the pH value of the seawater, water temperature, and water move-ment. However, it is not clear how these environmental ﬂuctuations affect the dynamical behavior of the phytoplankton and zooplankton system. In this paper, a stochastic eco-epidemiological model for viral infection in the toxin-producing phytoplankton and zooplankton system is proposed. Firstly, the existence and uniqueness of globally positive solutions for this model is shown. Secondly, the stochastic boundedness of solutions for the model is proved. Moreover, sufﬁcient conditions for the extinction and persistence in the mean for the phytoplankton and zooplankton are obtained by constructing appropriate stochastic Lyapunov functions and using analytical techniques. Numerical simulations are carried out to demonstrate different dynamical behaviors including coexistence, extinction of the whole plankton system, partial persistence and extinction, and their corresponding probability density curves.


Introduction
Plankton, referring to the freely floating and weakly swimming organisms in some bodies of fresh water, are the microscopic organisms of the aquatic environment and play an significant role in the marine ecology. The types of plankton are usually divided into phytoplankton (marine microscopic algae) and zooplankton (small marine animals and immature stages of larger animals mainly feed by the phytoplankton). Phytoplankton, living near the surface of the ocean, have the important effect of balancing the ecosystem by producing oxygen and absorbing carbon dioxide through photosynthesis. The growth of phytoplankton mainly depends on some most important nutrients such as nitrogen and phosphorus, which are naturally available in aquatic environments in various concentrations. Therefore, nutrients, phytoplankton, and zooplankton formed the basis of all aquatic food chains in the marine ecosystem. Due to the eutrophication of water with the presence of too many nutrients, phytoplankton may increase rapidly, and a high concentration of phytoplankton can cause the water to appear blue-green, green, brown, or even red. This phenomenon is popularly known as the algal blooming. Moreover, the nutrients in fresh water are quickly exhausted with algae growing fast. Thus, a low nutrient concentration will inevitably restrict the growth of phytoplankton. Meanwhile, in the process of decay Here, S represents the concentrations of phytoplankton (units of cells/l), I denotes the density of infected cells, and Z (units of animals/l) represents the concentrations of zooplankton. Assume that grazing will follow a mass-action law. In order for infected cells to represent a small fraction of the total biomass, the rate of contagion is made proportional to an interaction term of the form SI S+I , which is roughly I for S large compared to I. The model results are roughly consistent with the behavior observed in field data by employing some methods of nonlinear forecasting. Zooplankton's mortality rate has been increased because of the predation of harmful toxin-producing phytoplankton (both infected and susceptible). Consequently, the linear functional response form, which becomes a naturally evident choice to depict the distribution of toxic substances, was employed in [15]. In addition, some works about the effect of viral disease in recurrent phytoplankton-zooplankton system have been discussed in [16][17][18] and the references therein. Later, Gakkhar and Negi in [19] discussed the effect of toxin substance by the following model where the transmission rate of virus in the phytoplankton is bilinear and θ is the rate of the toxin liberation parameter. The theoretical results illustrated that the system exhibited periodic oscillations in a small region for the rate of infection parameter and rate of toxin liberation parameter. Moreover, the author in [20] further investigated the delay model about the digestion of the zooplankton after feeding on phytoplankton. It was shown that a time delay could result in the appearence of Hopf bifurcation. Moreover, there are some studies about the infection delay [21] and the distribution delay [22], taking the nutrient compartment into consideration [23,24], and involved in virus compartment and environment compartment [25].
Nowadays, it is widely accepted that the evolution of natural populations is rarely deterministic because it is inevitably disturbed by various environmental factors, and the marine ecosystem is no exception. Specifically, phytoplankton rely on light and carbon dioxide for photosynthesis. Therefore, the growth of phytoplankton might be disturbed by some environment factors such as nutrient availability, water temperature, water depth, and light. The growth and development of zooplankton is more affected by stochastic fluctuations of the environment, such as the pH value of the seawater, water temperature, and water movement. On account of the uncertainty of these physical factors (such as nutrient availability, water temperature, light, etc.), the environmental fluctuation is ubiquitous in aquatic ecosystems and prominently affects the dynamics of aquatic ecosystems. In [26], the authors assumed that the growth rate of phytoplankton and the death rate of zooplankton are subjected to the Gaussian white noise, and they studied the following stochastic model: . Sufficient conditions for the existence of globally positive solutions as well as population extinction and persistence in the mean for solutions were established. Furthermore, there exists a single stationary distribution with ergodicity. In [27], the authors established a stochastic phytoplankton-zooplankton model with toxin-producing phytoplankton under regime switching including white and colored noises. They proved the existence of globally positive solutions of the model and derived sufficient conditions for the existence of a unique ergodic stationary distribution. There were some works about introducing the environmental noise into ecosystems (see [26][27][28][29][30][31][32][33][34][35]).
Motivated by these work mentioned above, we consider the following phytoplanktonzooplankton deterministic model with toxin-producing phytoplankton and viral infection where S(t), I(t), and Z(t) stand for the population size at time t of the susceptible phytoplankton, infected phytoplankton, and zooplankton, respectively. K is the carrying capacity of phytoplankton density; r is the intrinsic rate of plankton growth; β is the rate of infection; c and d are the mortality rate of the infected phytoplankton population and zooplankton population, respectively; p 1 , p 2 represent the capture rates of susceptible and infected phytoplankton for zooplankton, respectively; k represents the conversion rate due to the predation of phytoplankton; θ is the rate of toxin liberation by the toxin-producing plankton. All parameters are positive constants. It can be simplified to a dimensionless form with a reduced number of parameters by means of transformation Then, we can obtain the following system (1) In order to describe the effect of stochastic fluctuations for environmental factors on the plankton system, we introduce stochastic fluctuations, which are of the white noise type and the intensities of the white noises are proportional to s(t), i(t), and z(t), which are the same as the modeling idea in [36][37][38]. Therefore, in this paper, we consider the following stochastic eco-epidemiological model for viral infection in a toxin-producing phytoplankton and zooplankton system where B i (t) (1, 2, 3) are mutually independent standard Brownian motions with B i (0) = 0, and σ 2 i > 0 denote the intensities of B i (t), respectively. The rest of the paper is arranged as follows. In Section 2, some useful lemmas are presented. In Section 3, the stochastic boundedness of solutions for model (2) are proved. In Section 4, we focus on the extinction and persistence in the mean of solutions for model (2). In Section 4, numerical simulations are given to validate these theoretical results. Finally, a brief conclusion and discussion are given in Section 5.

Preliminaries
b} for all a, b ∈ R + . Throughout this paper, we let (Ω, {F t } t≥0 , andP) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F t contains all P-null sets). We assume that model (2) is defined on (Ω, {F t } t≥0 , P).
Firstly, on the existence and uniqueness of globally positive solutions for model (2), we have the following lemma.

Proof.
Since the coefficients of model (2) are locally Lipschitz continuous, for any given initial value (s 0 , i 0 , z 0 ) ∈ R 3 + , there exists a uniquely maximal local positive solution (s(t), i(t), z(t)) on t ∈ [0, τ e ) a.s., where τ e is the explosion time. To show this solution is global, we need to prove τ e = ∞ a.s.. Let k 0 > 0 be sufficiently large for s 0 , i 0 and z 0 lying within the interval [ 1 k 0 , k 0 ]. For each integer k ≥ k 0 , define the stopping time as follows Throughout this paper, we let inf ∅ = ∞ for empty set ∅. Clearly, τ k is increasing as s. If we can show that τ ∞ = ∞ a.s., then τ e = ∞ a.s.. Suppose τ ∞ < ∞ a.s.; then, there exists a pairs of constants T > 0 and ∈ (0, 1) such that P{τ ∞ ≤ T} > . Hence, there is an integer k 1 ≥ k 0 such that Define a C 2 function: The non-negativity of this function can be seen from u − 1 − ln u ≥ 0 for any u > 0. Using the Itô formula, we obtain 2 . By the inequality u ≤ 2(u + 1 − log(u)) − (4 − 2 log 2) for all u > 0 (see [22]), we have 2(u + 1 − log(u)) ≥ 2 + u ≥ u. As a result, it can be obtained that Integrating this inequality from 0 to τ k ∧ T and then taking the expectation on both sides, we obtain by the Gronwall inequality, we conclude that So, we obtain τ ∞ = ∞ a.s. This completes the proof.
Next, we introduce the following lemma, which can be proved by the similar methods given in [33]. Here, we omit the proof.

Lemma 2.
A one-dimensional Brownian motion W(t) t≥0 has the following property for almost every ω ∈ Ω.

Proof. From model
Considering the following auxiliary system It is easy to obtain by Theorem 1.16 in [39] that system (4) has a unique solution and it is easy to verify lim sup t→∞ E[Φ(t)] ≤ 1. Hence, by the comparison theorem of stochastic differential equation [40], we obtain s(t) ≤ Φ(t) for all t ≥ 0 a.s.. and thus That is, s(t) is uniformly bounded in mean. In the following, we will prove that i(t) and z(t) are also uniformly bounded in mean. Define a function Using the Itô formula, we obtain Integrating the last inequality in the above formula from 0 to t, we can obtain Taking the expectation on both sides, we obtain .
that is, W(t) is uniformly bounded in mean a.s.; therefore, i(t) and z(t) are also uniformly bounded in mean a.s.
Moreover, since s(t) is uniformly bounded in mean a.s., by Markov's inequality (see [41]), we assert that for any positive constant η, there is a constant β = 1 η > 0 such that hence, by (3), we have lim sup t→∞ P(s > η) ≤ β, a.s.; thus, s(t) is stochastically ultimately bounded. Similarly, we can show that i(t) and z(t) of model (2) is also stochastically ultimately bounded. This finishes the proof.
Corollary 1. There exist three positive constantss,ī, andz such that for any initial value

Extinction and Persistence of the Plankton
For the extinction and persistence of the plankton in the mean, we have the following result.
Integrating the last inequality of (5) from 0 to t and then dividing by t on both sides, we get log By the large number theorem [41] , we get lim t→∞ < 0, a.s., that is, s(t) exponentially tends to zero with probability 1 as t → ∞; furthermore, we obtain lim t→∞ 1 t t 0 s(u)du = 0, a.s.. Meanwhile, Integrating this inequality from 0 to t and then dividing by t on both sides, we get By (7), we further obtain lim sup Hence, i(t) exponentially tends to zero with probability 1 as t → ∞ and thus, we obtain Similarly, by the third equation of model (2), we have by (7) and (8), we further get i.e., z(t) exponentially tends to zero with probability 1. Secondly, by applying the Itô formula, we have We further obtain that by the variation-of-constants formula [41] 1 On the one hand, by (9), we have which follows that where k 1 (t) = 1 s 0 exp( and it is easy to obtain sup 0≤t<∞ k 1 (t) < ∞.
And by the inequality in (5), we have Moreover, by (10) and (11), we obtain If the condition a > σ 2 1 2 holds, by taking the superior limit on both sides in the above formula and combining with large number theorem [41] and Lemma 2, we have lim sup On the other hand, by (9), we have By (11), we further obtain As a result of the stochastic ultimate boundedness of solution (s(t), i(t), z(t)) for model (2), we obtain that for an arbitrary small constant 0 < ζ < 1, there exists a constant t ζ > 0 such that P(Ω) ≥ 1 − ζ, whereΩ = {ω : 1 a [(a + 1)i(t, ω) + l 1 z(t, ω)] ≤ ζ for all t ≥ t ζ }. Thus, for any ω ∈Ω, by taking the inferior limit on both sides in (13) and combining with Lemma 2, it is shown that Letting ζ → 0, we get lim inf t→∞ hence, s(t) is persistent in the mean with probability one. Moreover, by the second equation of model (2), we have By (14), we further lim sup It means that lim sup t→0 i(t) = 0 a.s., by combining with the positivity of the solution, we obtain lim t→0 i(t) = 0 a.s.. (15) Similarly, by the last equation of model (2), we get By (14) and (15), we have From the condition l 1 < n, we finally have lim sup t→∞ log z(t) t < 0, and thus, lim t→∞ z(t) = 0. This completes the proof. Remark 1. From Theorem 2, it is suggested that both species in model (2) will die out if σ 2 1 > 2a; that is, when the environment noise intensity for susceptible phytoplankton is larger, the aquatic ecosystem will finally be destroyed, while the environment noise intensity for susceptible phytoplankton is smaller and certain conditions are satisfied, i(t), z(t) in model (2) goes to extinction with probability one, respectively. At the same time, the phytoplankton persists in the mean with probability one; more precisely, phytoplankton tends to its own steady state. It is unfortunate that the persistence of s(t), i(t) and z(t) together could not be obtained in this paper. In [29], the authors investigated a two-dimensional phytoplankton-zooplankton model with regime-switching and toxin-producing phytoplankton but no viral infection in phytoplankton, which is different with model (2). Sufficient conditions for the existence of a unique ergodic stationary distribution were provided in [29]. The question of whether this approach about proving the existence of a unique ergodic stationary distribution in [29] could be used to our model will continue to be discussed in our future work.

Numerical Simulations
In this section, by utilizing Milstein's higher-order method [42], we will make some numerical examples to investigate the behavior of model (2) and illustrate the effects of environmental noises on the survival of plankton, which are performed by using MATLAB and R software.
Example 1. We firstly provide the numerical simulation of the coexistence of susceptible phytoplankton s(t), infected phytoplankton i(t), and zooplankton z(t). Under the parameter values in Table 1 and σ i = 0.01 (i = 1, 2, 3, ) the stochastic paths of s(t), i(t), and z(t) reveal that the susceptible phytoplankton, infected phytoplankton, and zooplankton are persistent in the mean with probability one. Based on 1000 sample paths, we obtain the probability densities for s(t), i(t), and z(t), as shown in Figure 1a-f. Example 2. Choose the parameters a = 0.0015, l 1 = 0.01, l 2 = 0.05, b 1 = 0.0067, n = 0.21, b 2 = 0.00775, σ 1 = 0.012, σ 2 = 0.02, and σ 3 = 0.013, respectively. It is easy to verify that the condition for the first assertion in Theorem 2 is not satisfied; that is, a − σ 2 1 2 = 0.014 > 0 . However, s(t), i(t), and z(t) exponentially tend to zero with probability one. It indicates that the condition of the extinction for the toxin-producing phytoplankton-zooplankton system in Theorem 2 is sufficient and not necessary. As shown in Figure 2, the toxin-producing phytoplankton-zooplankton system would die out. Example 3. Choose n = 0.001, b 2 = 0.95, σ 1 = 0.01, σ 2 = 0.31, and σ 3 = 0.01 but the rest of the parameters are taken from Table 1. By calculating, it can be obtained that a − σ 2 .0019 > 0, and l 1 − n = 0.01 > 0, which indicate that the conditions for the second assertion in Theorem 2 are not satisfied. However, as Figure 3a-f show, s(t) is persistent, while in the meanwhile, i(t) and z(t) would be extinct with probability one, which suggests that the condition of the second assertion in Theorem 2 may be too strong for the persistence only for susceptible phytoplankton but the extinction for infected phytoplankton and zooplankton in the aquatic system.
From the Examples 2 and 3, it can be deduced that the susceptible phytoplankton might be either persistent or extinct if condition a > σ 2 1 2 holds, which also implies that even though the self-growth of phytoplankton is relatively rapid and the effect of environmental noise is small, the phytoplankton and zooplankton system may be extinct or partially extinct. The survival of the phytoplankton and zooplankton system may be subject to other ecological factors, which is a very interesting question deserving our further investigation.     Table 1 and σ i = 0.01 (i = 1, 2, 3). The initial value is (0.08, 0.06, 0.05). (a-f) denote the evolution of a single path and the probability density curve (based on the 1000 simulations) of s(t), i(t), and z(t), respectively.

Example 4.
To further explore the effect of the noises intensities σ 1 , σ 2 , and σ 3 on the dynamical behavior of model (2), hereon, we mainly focus on the effect of environmental variability on the extinction time for z(t) in model (2), which is more meaningful. Now, we define a random variable X to describe the time when z(t) approaches zero, a.s., namely, X inf{t|z(s) < 10 −5 , a.s., for s ≥ t}. We will pay close attention to the dependence of X on three noise intensities σ 1 , σ 2 , and σ 3 . We fix the other parameter values from Table 1, except for b 2 = 0.73 and σ i = 0.01 (i = 1, 2, 3). Under these parameters, s(t) and i(t) are persistent with probability one, while z(t) is extinct with probability one, as shown in Figure 4. By using the Latin hypercube method to generate 1000 random samples of the parameter set (σ 1 , σ 2 , σ 3 ), we compute the partial rank correlation coefficients (PRCCs) for the average value of X over 1000 simulations with respect to each noise intensity, As illustrated in Figure 5, it can be found that these noise intensities have a significantly negative impact for the extinction time of z(t) (under the significance level 0.05); moreover, the absolute value of the PRCC for σ 3 is the largest, which is suggested that σ 3 is more sensitive than the others for the extinction of z(t).            Table 1 except for b 2 = 0.73 and σ i = 0.01 (i = 1, 2, 3). The initial value is (0.08, 0.06, 0.05). (a-f) denote the evolution of a single path and the probability density curve (based on the 1000 simulations) of s(t), i(t), and z(t), respectively.  Table 1, s(t) and i(t) are persistent with probability one, while z(t) is extinct with probability one. The initial values of the solutions are (0.08, 0.06, 0.05).

Discussion and Conclusions
Actually, viral infection commonly exists in aquatic ecosystems. However, there are few models for viral infection in the toxin-producing phytoplankton and zooplankton systems. In this paper, we focus on a stochastic model of viral infection in a toxinproducing phytoplankton and zooplankton system due to considering many stochastic fluctuations in the phytoplankton-zooplankton system (for instance, the randomness of environmental factors).
In this paper, we investigated the stochastic boundness of the solutions for stochastic eco-epidemiological model for viral infection in the toxin-producing phytoplankton and zooplankton system. Sufficient conditions of extinction for the phytoplankton and zooplankton system and persistence in the mean for susceptible phytoplankton but extinction for infected phytoplankton and zooplankton are obtained. It is worth noting that these conditions are sufficient but not necessary, and the persistence for susceptible phytoplankton, infected phytoplankton, and zooplankton together are not obtained in the theoretical analysis, which is left for our further study. In addition, the phytoplankton live in a periodically varying environment; for instance, tidal circulation determines the mobility of aquatic organisms, and light intensity controls the rate of photosynthesis of phytoplankton. Hence, it is natural to incorporate the seasonal fluctuation into model (2) and investigate the effects of stochasticity and seasonality on the dynamic behavior for the phytoplankton and zooplankton system. Moreover, the time delay due to the gestation of zooplankton also could affect the dynamical behavior of the phytoplankton-zooplankton system. Therefore, it is interesting to consider a stochastic model of viral infection in toxin-producing phytoplankton and zooplankton system incorporating time-delays. These models are realistic but more complex, which deserves our further consideration.