# Stationary Distribution and Periodic Solution of Stochastic Toxin-Producing Phytoplankton–Zooplankton Systems

Authors: Chunjin Wei; Yingjie Fu
Journal: Complexity
 (2020)

---

## Abstract

In this paper, we investigate the dynamics of autonomous and nonautonomous stochastic toxin-producing phytoplankton–zooplankton system. For the autonomous system, we establish the sufficient conditions for the existence of the globally positive solution as well as the solution of population extinction and persistence in the mean. Furthermore, by constructing some suitable Lyapunov functions, we also prove that there exists a single stationary distribution which is ergodic, what is more important is that Lyapunov function does not depend on existence and stability of equilibrium. For the nonautonomous periodic system, we prove that there exists at least one nontrivial positive periodic solution according to the theory of Khasminskii. Finally, some numerical simulations are introduced to illustrate our theoretical results. The results show that weaker white noise and/or toxicity will strengthen the stability of system, while stronger white noise and/or toxicity will result in the extinction of one or two populations.


---

## Body

## 1. Introduction

As well known, mathematical models describing the plankton dynamics have played an important role in understanding the various mechanisms involved in toxin-producing phytoplankton. There are many scientific works have been carried out to investigate the effects of toxin-producing phytoplankton on plankton ecosystems [1–7]. For example, Upadhyay and Craniopathy [1] proposed three species food chain model with different functional forms to describe the liberation of toxin. The obtained results show that the increase of toxic substances released by toxic-producing phytoplankton has a stabilizing effect. In particular, according to field observations, Chattopadhyay et al. [5] formulated the following toxin-producing phytoplankton–zooplankton model:(1)dpdt=rp1−pk−αfpz,dzdt=βfpz−μz−θgpz,wherept and zt denote the density of toxin-producing plankton (TPP) population and the zooplankton population at time t, respectively, subject to the nonnegative initial condition p0=p0≥0 and z0=z0≥0. r and k represent the intrinsic growth rate and the environmental carrying capacity of TPP population, respectively. α is the rate of predation of zooplankton on TPP population, β is the ratio of biomass consumed by zooplankton for its growth (satisfying the obvious restriction 0<β<α), and μ denotes the mortality rate of zooplankton due to nature death, θ denotes the rate of toxin liberation by TPP population. fp represents the predational response function and gp describes the distribution of toxic substances. All parameters above are positive. As liberation of toxin reduces the growth of zooplankton, causes substantial mortality of zooplankton and in this period toxin-producing phytoplankton population is not easily accessible, hence a more common and intuitively obvious choice is of the saturation functional form to describe the grazing phenomena. For instance, Tapan et al. [8] studies the system (1) when fp=gp=p/γ+p (where γ>0 denotes the half-saturation constant). The obtained result indicates that there is a threshold limit of toxin liberation by the phytoplankton species below which the system does not have any excitable nature and above which the system shows excitability.Clearly, these important and useful works on deterministic phytoplankton–zooplankton model provide a great insight into the dynamics of plankton ecosystems. However, in the real world, the dynamics of plankton ecosystems are inevitably perturbed by various types of environment noises. Development from the deterministic models to the stochastic models can give us new insights into the dynamics of plankton ecosystems [9]. May [10] pointed out that the birth rates, carrying capacity, competition coefficients and other parameters involved in the system can be affected by environmental noise. This may be especially true for plankton ecosystems due to unpredictability of photosynthetically active radiation, nutrient availability, water temperature, water depth, light, eutrophication, acidity, salinity, wind and many other physical factors embedded in aquatic ecosystems. Several scholars have studied the effect of environmental fluctuations on aquatic ecosystems [11–13]. Indeed, stochastic models could be more appropriate way of modeling in comparison with their deterministic counterparts, since they can provide some additional degree of realism. By introducing (stochastic) environmental noise, many investigators have studied stochastic epidemic models [14–23] and stochastic population models [24–42]. They focus on the effect of environmental fluctuations on the dynamic behavior of these models. For instance, Silva [11] investigated a stochastic model of phytoplankton–zooplankton interactions with toxin-producing phytoplankton. Theoretical results show that for certain values of the system parameters, the system posses asymptotic stability around the positive interior equilibrium which depicts the coexistence of all the species. Therefore, it is meaningful to further incorporate the environmental stochasticity into the model (1) with fp=gp=p/γ+p, which could provide us a deeper understanding for real aquatic ecosystems [43].On the other hand, periodic behavior arises naturally in many real world problems, such as in biological, environmental and economic systems [44]. The phenomenon of periodic oscillations has been observed in the growth of populations, such as, a seasonal occurrence of Hematodinium perezi was reported by Newman and Johnson (1975) in the parasitic dinoflagellate in blue crabs on the east coast of the United States from late spring to early winter [45]. Seasonal changes are cyclic, largely predictable, and arguably represent the strongest and most ubiquitous source of external variation influencing human and natural systems [46]. However, to the best of our knowledge, there is little work on the existence of stochastic periodic solution for nonautonomous toxin-producing phytoplankton zooplankton model. Based on the aforementioned, we intend to study toxin-producing phytoplankton–zooplankton model with environmental fluctuations, and then we extend this model into a nonautonomous stochastic model by taking into account seasonal variation in the next section. Then in Section 3, we show the existence and uniqueness of the global positive solution. In Section 4, we obtain the sufficient conditions for the solution of population extinction and persistence in the mean. In Section 5 and Section 6, using Khasminskii’s methods and Lyapunov functions, we derive sufficient conditions for the existence of the single ergodic stationary distribution of the autonomous system (2) and the existence of the nontrivial positive stochastic periodic solution of nonautonomous system (3). In Section 7, some numerical simulations are provided to demonstrate the analytical findings. Finally, some conclusions are given in Section 8.


## 2. Model Formulation

There are many kinds of approaches to introduce the white noise into the population models. For model (1), we assume that the growth rate r of phytoplankton and the death rate μ of zooplankton are subjected to the Gaussian white noise (we follow the way used in [38]), then we can obtain the following stochastic model:(2)dp=rp1−pk−αpγ+pzdt+pσ1dB1t,dz=βpγ+pz−μz−θpγ+pzdt+zσ2dB2t,whereB1t,B2t are independent Brownian motions defined on the complete probability space. Ω,F,Ftt≥0,P with a filtration Ftt≥0 satisfying the usual conditions (i.e. it is right continuous and increasing while F0 contains all P-null sets). σ12,σ22 represent the intensities of the white noise. Meanwhile, there are evidences suggesting that the toxic substances released by TPP do not remain constant but change over time, which is related to the seasonal changes. Therefore, for better understanding the toxin-producing phytoplankton–zooplankton sustained oscillatory patterns, we further consider the following periodic system with stochastic perturbation by the method of Khasminskii [47].(3)dp=rtp1−pkt−αtpγt+pzdt+pσ1tdB1t,dz=βtpγt+pz−μtz−θtpγt+pzdt+zσ2tdB2t,wherert,kt,αt,βt,μt,θt,γt are all positive θ1-periodic continuous functions. For biological significance, we always assume that βt>θt. Throughout this paper, let R+=0,∞,R+2=x=x1,x2∈R2:xi>0,i=1,2,ux=∂u/∂x,uxx=∂2u/∂x2, a.s. means almost surely.


## 3. Existence and Uniqueness of the Global Positive Solution

As we know, in order for a stochastic differential equation to have a single global solution (i.e. no explosion in a finite time) for any given initial value, the functions involved with stochastic system are generally required to satisfy the linear growth condition and local Lipschitz condition [48, 49]. However, the functions of system (2) do not satisfy the linear growth condition, so the solution of system (2) may explode at a finite time. In this section, we show that there exists a single positive local solution of system (2), then using the Lyapunov analysis method, we prove that this solution is global. Explanation for “explosion time” used in following lemma.Lemma 1.
For
p
0
,
z
0
∈
R
+
2
, there exists a single positive local solution
p
t
,
z
t
of system
(2)
for
t
∈
0
,
τ
e
almost surely, where
τ
e
is the explosion time.Proof.
Using the transformation of variablesut=lnpt,vt=lnzt and applying Itô’s formula, for system (2), we have.(4)
d
u
t
=
r
−
σ
1
2
2
+
r
k
e
u
t
−
α
e
v
t
γ
+
e
u
t
d
t
+
σ
1
d
B
1
t
,
d
v
t
=
β
e
u
t
γ
+
e
u
t
−
μ
−
σ
2
2
2
−
θ
e
u
t
γ
+
e
u
t
d
t
+
σ
2
d
B
2
t
,with the initial valueu0=lnp0,v0=lnz0. The functions involved with drift part of above stochastic differential system satisfy the linear condition and locally Lipschitz condition. Hence there exists a single local solution ut,vt, for t∈0,τe, where τe is any finite positive real number. Clearly, pt,zt=eut,evt is the single positive local solution of stochastic differential system (2) starting from an interior point of the first quadrant.Now we are in a position to show that this single solution is a global solution. To prove this we only need to show thatτe=∞ a.s.Theorem 2.
For any initial value
p
0
,
z
0
∈
R
+
2
, the system
(2)
has a single global positive solution
p
t
,
z
t
for all
t
>
0
, and the solution will remain in
R
+
2
with probability one.Proof.
By Lemma1, we only need to prove that τe=∞ a.s. From the first equation of model (2), we have(5)
d
p
t
≤
p
t
r
−
r
k
p
t
d
t
+
σ
1
p
t
d
B
1
t
.Let(6)Φt=er−σ12/2t+σ1B1t1p0+rk∫0ter−σ12/2s+σ1B1sds,thenΦt is the single solution of the following auxiliary system:(7)dΦt=Φtr−rkΦtdt+σ1ΦtdB1t,dΦ0=P0.Then by the comparison theorem for stochastic equations, we have(8)pt≤Φt,t∈0,τea.s.On the other hand, from the second equation of model (2), we have(9)dzt≤ztβ−μdt+σ1ztdB1t.Consider the following auxiliary system:(10)dMt=Mtβ−μdt+σ2MtdB2t,dM0=z0,ThenMt=z0eβ−μ−σ22/2t+σ2B2t is the single solution of above system. Similarly by the comparison theorem for stochastic differential equations, we have(11)zt≤Mt,t∈0,τea.s.Similarly, we have(12)pt≥Ψt,t∈0,τea.s.,whereΨt=er−σ12/2t−α/γ∫0sMsds+σ1B1t/1/p0+r/k∫0ter−σ12/2s−α/γ∫0sMudu+σ1B1sds is the solution of the following system:(13)dΨt=Ψtr−rkΨt−αMtγdt+σ1ΨtdB1t,dM0=p0,By the same arguments as above, we have(14)dz≥ztβ−βγγ+pt−μ−θdt+σ2dB2t≥ztβ−βγγ+Ψt−μ−θdt+σ2ztdB2t.Similarly, we can obtain(15)zt≥Wt,t∈0,τea.s.,whereWt=z0eβ−μ−θ−σ22/2t−∫0tγβ/γ+Ψtds+σ2B2t is the single solution of the following system:(16)dWt=Wtβ−βγγ+Ψt−μ−θdt+σ2WtdB2t,dW0=z0,By (8), (11), (12), (15), we can get that(17)Ψt≤pt≤Φt,Wt≤zt≤Mt,t≥0a.s.It follows from [34] that Ψt,Φt,Wt and Mt will not be exploded at any finite time, then by the comparison theorem of stochastic differential equations, we can derive that pt,zt will globally exist. This completes the proof.Remark 3.
For any initial valuep0,z0∈R+2, the system (3) has a single global positive solution pt,zt, for all t>0, and the solution will remain in R+2 with probability one.Proof.
The proof is similar to Theorem2, so we omit it.Remark 4.
When we choosefp=pn/γn+pn,gp=pn/γn+pnn=1,2, then the stochastic system (2) has a single global positive solution pt,zt for all t>0, and the solution will remain in R+2 with probability one.Proof.
The proof is similar to Theorem2, so we omit it.


## 4. Extinction and Persistence in Mean

In this section, we will investigate the persistence and extinction of the system (2) under certain conditions. We give the definition and lemma which can be used for our main results.Definition 5 [38].
(1)
System (2) is said to be extinct if limt→∞pt=0a.s.;
(2)
System (2) is said to be persistent in mean if limt→∞1/t∫0tpsds>0a.s.LetCΩ×0,∞,R+ be the family functions on Ω which are continuously differentiable with respect to x∈Ω and continuously differentiable with respect to t∈0,∞.Lemma 6 [35].
Suppose thatyt∈CΩ×0,∞,R+, if there are positive constant T,μ0,μ>0 such that for t≥T.(1)
lnyt≤μt−μ0∫0tysds+∑i=1nβiBit, then limsupt→∞1/t∫0tysds≤μ/μ0 a.s.;(2)
lnyt≥μt−μ0∫0tysds+∑i=1nβiBit, then liminft→∞1/t∫0tysds≥μ/μ0 a.s.,
whereβi1≤i≤n is a constant.Theorem 7.
Letpt,zt be the solution of the system (2), then the following statements hold:(i)
Ifσ12<2r,σ22>2β−2μ, then limt→∞zt=0,limt→∞1/t∫0tpsds=kr−σ12/2/r a.s.;(ii)
Ifσ12>2r,σ22>2β−2μ, then limt→∞zt=0,limt→∞pt=0 a.s.Proof.
(i) By the positiveness of the solutionpt,zt of system (2), we have(18)
d
z
=
z
β
p
γ
+
p
−
μ
−
θ
p
γ
+
p
d
t
+
z
σ
2
d
B
2
t
≤
z
β
−
μ
d
t
+
z
σ
2
d
B
2
t
.Consider the following system(19)dΨt=Ψβ−μdt+Ψσ2dB2t,dΨ0=z0.We have(20)dlnΨt=β−μ−σ222dt+σ2dB2t.By the comparison principle of the stochastic differential equation and the theory of diffusion process (see [33] and Lemma A.2 in [34]), we can easily check that(21)S−∞=∫0−∞exp−∫0v2β−μ−σ22/2σ22dsdv=σ222β−μ−σ22/2>−∞,S+∞=∫0+∞exp−∫0v2β−μ−σ22/2σ22dsdv=+∞.Sinceσ22>2β−2μ, we can obtain that limt→∞lnzt=−∞, that is limt→∞zt=0 a.s. For any ϵ>0, there exisit T>0 and a set Ωϵ such that PΩϵ>1−ϵ and zt<ϵ for t>T,ω∈Ωϵ.By the first equation of system (2), we have(22)r−σ122−αϵγ−rkpdt+σ1dB1t≤dlnp≤r−σ122−rkpdt+σ1dB1t.By Lemma6, for ∀ϵ→0, we have liminft→∞1/t∫0tpsds≥kr−σ12/2/r, and limsupt→∞1/t∫0tpsds≤kr−σ12/2/r, that is to say(23)limt→∞1t∫0tpsds=kr−σ12/2r  a.s.  .(ii) Similar to the proof of (i), we obtain that limt→∞zt=0 a.s. for σ22>2β−2μ. In addition, by the first equation of system (2) and r−σ12/2<0, we obtain limt→∞pt=0a.s.. This completes the proof.Remark 8.
From Theorem7, we can see, when σ2 is a constant, if σ12<2r, then the population pt will be persistent in mean, but when we choose σ1 large enough such that σ12>2r, then the population pt will go to extinct. From an ecological point of views, the intensity of white noise has a negative effect on the survival of pt population, which imply that weaker white noise will strengthen the stability of the system, while stronger white noise will lead to population pt extinct.


## 5. Stationary Distribution and Ergodicity

In this section, we shall consider whether there exists a single stationary distribution of system (2), which means that the zooplankton population can persist and not die out.LetXt be a regular time-homogeneous Markov process in R+d described by the following stochastic differential equation(24)dXt=fXtdt+∑r=1kσrXtdBrt.The diffusion matrix is defined as follows(25)AX=aijx,aijx=∑r=1kσrixσrjx.Next, we shall introduce a lemma which guarantees the existence and uniqueness of a stationary distribution and ergodicity (see Khasminskii [47]).LetU be a given open set in the d-dimensional Euclidean space Rd and C2 denotes the class of functions in Rd which are twice continuously differentiable with respect to x.Lemma 9.
The Markov process
X
t
has a unique stationary distribution
π
·
, if there exists a bounded domain
U
⊂
R
d
with regular boundary
Γ
such that its closure
U
¯
⊂
R
d
, having the following properties:
(B.1) In the open domain
U
and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix
A
t
is bounded away from zero.
(B.2) Ifx∈Rd\U, the mean time τ at which a path issuing from x reaches the set U is finite, and supx∈KExτ<∞ for every compact subset K⊂Rd. Moreover, if f· is a function integrable with respect to the measure π, then(26)
P
l
i
m
T
→
∞
1
T
∫
0
T
f
X
x
t
d
t
=
∫
R
d
f
x
π
d
x
=
1
,
∀
x
∈
R
d
.Remark 10.
To prove condition (B.1), it suffices to verify thatF is uniformly elliptical in U, where Fu=fxux+1/2TrAxuxx,i.e., there exists a positive number M¯such that∑i,j=1daijxξiξj>M¯ξ2,x∈U,ξ∈Rd. To verify condition (B.2), it is sufficient to prove that exists a nonnegative C2−function Vx and a neighborhood U such that for some positive constant M,LVx≤−M for any x∈Rd\U.System (2) can be written into the following form:(27)dptzt=pr−rkp−αzγ+pzβpγ+p−μ−θpγ+pdt+σ1pt0dB1t+0σ2ztdB2t.Suppose the following condition (H) holds:(28)Hλ≜βr1+γ/kr−σ122−μ−σ222−θ>0.Theorem 11.
Assume the condition (H) holds and
σ
i
2
>
0
i
=
1
,
2
, then for any initial value
p
0
,
z
0
∈
R
+
2, system (2) has a single stationary distribution and it has ergodic property.Proof.
By Theorem2, we have obtained that for any initial value p0,z0∈R+2, the system (2) has a single global positive solution pt,zt∈R+2.The following proof is motivated by Liu et al. [38], if we substitute pt=eξt,zt=eηt into system (2) and then by Itô’s formula, we get(29)dξt=r−σ122−rkeξt−αeηtγ+eξtdt+σ1dB1t,dηt=−μ−σ222+βeξtγ+eξt−θeξtγ+eξtdt+σ2dB2t.Define aC2-function V:R2→R+(30)Vξ,η=M−ηt+βr1+γ/klnγ+eξt−ξt+αβ/γ+θr+r/kμβeηt+eξt+meηtq+1q+1=V1+V2,wherem=α/β, q is a constant satisfying 0<q<min1,2μ/σ22, M=2/λmax2,supξ,η∈R2−r/4ke2+qξ−m1+qμ−σ22q/2/4e1+qη+N, N=supξ,η∈R2−r/2ke2+qξ−m1+qμ−σ22q/2/2e1+qη+reξ+meηeξ+σ12q/  2e1+qξ. Clearly, Mλ/4≥1.An application of the operator to−ηt and lnγ+eξt−ξt+αβ/γ+θr+r/k/μβeηt leads to(31)L−η=μ+σ222+θeξγ+eξ−βeξγ+eξ=μ+σ222+θeξγ+eξ−βeξγ+eξ+β1+γ/k−β1+γ/k≤μ+σ222+θ−β1+γ/k+βr1+γ/k·γr−r/keξγ+eξ,and(32)Llnγ+eξ−ξ+αβ/γ+θr+r/  kμβeη=eξγ+eξr−σ122−rkeξ−αeηγ+eξ+−r+σ122+rkeξ+αeηγ+eξ−γσ12eξ2γ+eξ2+αβ/γ+θr+r/  kσ222μβeη+αβ/γ+θr+r/  kμβeη−μ−σ222+βeξγ+eξ−θeξγ+eξ≤eξγ+eξr−rkeξ−αeξ+ηγ+eξ2−γσ12eξ2γ+eξ2+−r+σ122+rkeξ+αeηγ+eξ+αβ/γ+θr+r/  kμβ−μeη+βeξ+ηγ+eξ≤r−rkeξ−γr−r/  keξγ+eξ+−r+σ122+rkeξ+αγeη+αβ/γ+θr+r/  kμβ−μeη+βγeξ+η≤σ122−γr−r/keξγ+eξ+αβ/γ+θr+r/  kγμ⋅eξ+η.Using the above two inequations,LV1 can be estimated as follows(33)LV1≤Mμ+σ222+θ−βr1+γ/kr−σ122+βr1+γ/k·αβ/γ+θr+r/kγμ·eξ+η=−Mλ+Mheξ+η,whereh=β/  r1+γ/  k·  αβ/γ+θr+r/  k/γμ.(34)LV2=eξ+meηqreξ−rke2ξ−mμeη−mθeξ+ηγ+eξ+mβ−  αeξ+ηγ+eξ+qeξ+meηq−1e2ξσ12+m2e2ησ222≤reξ+meηqeξ−rke2+  qξ−m1+qμe1+  qη+qe1+  qξσ12+m1+qe1+  qησ222≤reξ+meηqeξ−rke2+  qξ−m1+qμ−σ22q2e1+  qη+σ12q2e1+  qξ=−r2ke2+  qξ−m1+qμ−σ22q/22e1+  qη+N,whereN=supξ,  η∈R2−γ/2  ke2+  qξ−m1+qμ−σ22q/  2/  2e1+  qη+reξ+meηqeξ+σ12q/  2e1+  qξ. Hence(35)LV=LV1+LV2≤−Mλ+Mheξ+η−r2ke2+  qξ−m1+qμ−σ22q/22e1+  qη+N.To confirm the condition (B.2) of Lemma9, we consider the bounded open subset(36)Uϵ=ξ,  η:ξ≤lnϵ−1,η≤lnϵ−1,ξ,  η∈R2,where0<ϵ<1 is a sufficiently small number. In the set Uc=R2\U, let us choose sufficiently small ϵ such that(37)0<ϵ<λ4h,(38)0<ϵ<m1+qμ−σ22q/24Mh,(39)0<ϵ<r4kMh,(40)−Mλ−r4kϵ+k1≤−1,(41)−Mλ−m1+qμ−σ22q/24ϵ+k2≤−1.For convenience, we divideUc into four domains,(42)Uϵ1=ξ,  η∈R2:−∞<ξ≤ln  ϵ,Uϵ2=ξ,  η∈R2:−∞<η≤ln  ϵ,Uϵ3=ξ,  η∈R2:ξ≥lnϵ−1,Uϵ4=ξ,  η∈R2:η≥lnϵ−1.Clearly,Uc=Uϵ1∪Uϵ2∪Uϵ3∪Uϵ4. Next, we will prove LVξ,  η≤−1 for any ξ,  η∈Uc.Case 1. In domain Uϵ1, owing to −∞<ξ≤lnϵ, then eξ+η≤ϵeη≤ϵ1+e1+  qη, we have(43)LVξ,  η≤−Mλ4+−Mλ4+ϵM  h+−m1+qμ−σ22q/24+ϵM  he1+  qη−r4ke2+  qξ+−Mλ2−r4ke2+  qξ−m1+qμ−σ22q/24e1+  qη+  N,≤−Mλ4+−Mλ4+ϵM  h+−m1+qμ−σ22q/24+ϵM  he1+  qη−r4ke2+  qξ+−Mλ2+supξ,  η∈R2−r4ke2+  qξ−m1+qμ−σ22q/24e1+  qη+  N,sinceMλ/4≥1, combing with the definition of M and Eqs. (37) and (38), we get(44)LVξ,  η≤−Mλ4−r4ke2+  qξ≤−Mλ4≤−1.Case 2. Similarly, for any ξ,  η∈Uϵ2, owing to eξ+η≤ϵeξ≤ϵ1+e2+  qξ, we have(45)LVξ,  η≤−Mλ4+−Mλ4+ϵM  h+−r4k+ϵM  he2+  qξ−m1+qμ−σ22q/24e1+  qη+−Mλ2−r4ke2+  qξ−m1+qμ−σ22q/24e1+  qη+  N,≤−Mλ4+−Mλ4+ϵM  h+−r4k+ϵM  he2+  qξ−m1+qμ−σ22q/24e1+  qη+−Mλ2+supξ,  η∈R2−r4ke2+  qξ−m1+qμ−σ22q/24e1+  qη+  N,by Eqs. (37) and (39), we also obtain(46)LVξ,  η≤−Mλ4−m1+qμ−σ22q/24e1+  qη≤−Mλ4≤−1.Case 3. In domain Uϵ3 we have(47)LV≤−Mλ+Mheξ+η−r2ke2+  qξ−m1+qμ−σ22q/22e1+  qη+N≤−Mλ−r4kϵ+k1,wherek1=supξ,  η∈R2−r/4  ke2+  qξ−m1+qμ−σ22q/  2/  2e1+  qη+Mheξ+η+  N, which gives LV≤−1 in view of Eq. (40).Case 4. In domain Uϵ4 we have(48)LV≤−Mλ+Mheξ+η−r2ke2+  qξ−m1+qμ−σ22q/22e1+  qη+N≤−Mλ−m1+qμ−σ22q/24ϵ+k2,wherek2=supξ,  η∈R2−m1+qμ−σ22q/  2/  4e1+  qη−r/2  ke2+  qξ+Mheξ+η+  N, which gives LV≤−1 according to Eq. (41). Consequently, we can deduce that(49)LVξ,  η≤−1,∀ξ,  η∈Uc.That is, the condition (B.2) holds.On the other hand, there existsM¯=minp,  z∈U¯σ12p2,σ22z2>0 such that(50)∑i,j=12ai,jp,zϵiϵj=σ12p2ξ12+σ22z2ξ22≥M¯ξ2,∀p,  z∈U¯,ξ=ξ1,ξ2∈R2.That is, the condition (B.1) is satisfied. Therefore, according to Lemma9, we know that the system (2) has a single stationary distribution which is ergodic.Remark 12.
From Theorem11, we can see that system (2) exists a single stationary distribution provided that the effects of both environment noise σii=1,2 and the rate of release of toxic substances θ are not too large such that λ>0; The ergodic property reflects the solution of system (2) converges to the single stationary distribution.


## 6. The Existence of Periodic Solution of Nonautonomous System

For convenience, we denote(51)fu=maxt∈0,θ1ft,fl=mint∈0,θ1ft,whereft is a continuous θ1-periodic function.In this section, we will recall a basic definition and introduce a lemma which gives a criteria for the existence of a periodic Markov process (see Khasminskii [47]).Definition 13 [47].
A stochastic processξt=ξt,  ω−∞<t<+  ∞ is said to be θ1-periodic if for every finite sequence of numbers t1,t2,···,tn, the joint distribution of random variables ξt1+  h,ξt2+  h,···,ξtn+  hh=kθ1k=1,2,··  · is independent of h.Consider the integral equation(52)Xt=Xt0+∫t0tbs,Xsds+∑r=1k∫t0tσrs,Xsdξrs,wherebs,  x,σis,  xi=1,2,···,  ks∈t0,  T,x∈Rl are continuous functions of s,  x and for some constant B, the following conditions hold.(53)bs,  x−bs,  y+∑r=1kσrs,  x,σrs,  y≤Bx−y,bs,  x+∑r=1kσrs,  x≤B1+x.LetU be a given open set in the d-dimensional Euclidean space Rd. E=Rd×0,  ∞, C2,1 is the family functions on E which are twice continuously differentiable with respect to x∈Rd and continuously differentiable with respect to t∈0,  ∞.Lemma 14 [47].
Suppose that the coefficients of Equation (52) are θ1-periodic in t and satisfy the conditions Equation (53) in every cylinder U×I, and assume further there exists a function Vx,  t∈C2,1, which is θ1-periodic in t and satisfies:(B.1)inf|x|>ℜVx,  t→∞asℜ→∞,(B.2)LVx,  t≤−1 outside some compact set,
then the system (52) exists at least a θ1-periodic Markov process.Theorem 15.
Assume the following conditions hold:(H.1)λ1=1/θ1∫0θ1βlru1+γu/klrs−βurl1+γl/kuσ12s2−μs+θs+σ22s2ds>0,(H.2)λ2=1/θ1∫0θ1μs−qsσ22s2ds>0,
then for any initial valuep0,z0∈R+2, the system (3) has a positive θ1-periodic solution.Proof.
By the same way as in Theorem2, we can obtain that, for any initial value p0,z0∈R+2, the system 3 has a single global positive solution. Next, we only need to verify the conditions (B.1), (B.2) of Lemma 14.Define aC2,1-function V:R2×0,  ∞→R+ by(54)Vξ,η,  t=M1−ηt+βurl1+γl/kulnγu+eξt−ξt+αuβu/γl+θuru+ru/klμlβleηt+ω1t+eω2teξt+mueηtq+1q+1=V1+V1,wherem=αt/βt, 0<q is a constant, M1=2/λ1max2,supξ,  η∈R2−rt/  4kteω2ue2+  qξ−λ2mu1+q/  4eω2ue1+  qη+N1, N1=supξ,  η∈R2−rl/  2kue2+  qξ−λ2mu1+q/  2e1+  qη+rueξ+mueηqeξ+σ12uq/  2e1+  qξ+e1+  qξ  Keω2u. Clearly, M1λ1/4≥1. Here, let(55)ω1′t=−λ1+βlru1+ru/klrt−βurl1+rl/kuσ12t2−μt+θt+σ22t2,ω2′t=−λ2+μt−qσ22t2,whichλ1,λ2 are defined by H.  1,H.  2. It is not difficult to prove that ω1t,ω2t are all θ1-periodic functions, and ω2′t is a bounded function. That is, there exists K>0 such that(56)ω2′t≤K,∀t>0.In order to confirm the condition (B.1) of Lemma14, it is sufficient to verify that(57)infξ,η,  t∈0,+  ∞×R+2\UkVξ,η,  t→∞,k→∞,whereUk=1/k,  k×1/k,  k, which is clearly established since the coefficients of the term e1+  qξt,e1+  qηt of Vξ,η,  t are all positive. Next, we only verify the condition (B.2) of Lemma 14.By Ito’s formula, we have(58)L−  η=μt+σ22t2+θteξγt+eξ−βteξγt+eξ=μt+σ22t2+θteξγt+eξ−βteξγt+eξ+βt1+γt/kt−βt1+γt/kt≤μt+σ22t2+θt−βt1+γt/kt+βtrt1+γt/kt·γtrt−rt/kteξγt+eξ,and(59)Llnγu+eξ−ξ+αuβu/γl+θurl+rl/kuμlβleη=eξγu+eξrt−σ12t2−rtkteξ−αteηγt+eξ+−rt+σ12t2+rtkteξ+αteηγt+eξ−γuσ12teξ2γu+eξ2+αuβu/γl+θuru+ru/klμlβleηt−μt−σ22t2+βteξγt+eξ−θteξγt+eξ+αuβu/γl+θuru+ru/klσ22t2μlβleη≤eξγu+eξrt−rtkteξ−αteξ+ηγu+eξ2−γuσ12eξ2γu+eξ2+−rt+σ12t2+rtkteξ+αteηγt+eξ+αuβu/γl+θuru+ru/klμlβl−μteη+βteξ+ηγt+eξ≤rt−rtkteξ−γtrt−rt/kteξγt+eξ+−rt+σ12t2+rtkteξ+αtγteη+αuβu/γl+θuru+ru/klμlβl−μteη+βtγteξ+η≤σ12t2−γtrt−rt/kteξγt+eξ−αuγl−αtγteη+αuβu2/γl+βuθuru+ru/klμlβlγl·eξ+η≤σ12t2−γtrt−rt/kteξγt+eξ+αuβu2/γl+βuθuru+ru/klμlβlγl·eξ+η.Using the above two inequations,LV1 can be estimated as follows:(60)LV1≤M1μt+σ22t2+θt−βlru1+γu/klrt−βurl1+γl/kuσ12t2+βurl1+γl/ku·αuβu2/γl+βuθuru+ru/klμlβlγl·eξ+η+ω1′t=−M1λ1+M1h1eξ+η,whereh1=βu/  rl1+γl/ku·  αuβu2/γl+βuθuru+ru/kl/μlβlγl, and(61)LV2=eω2teξ+mueηq+1q+1·ω2′t+eξ+mueηqrteξ−rtkte2ξ−muμteη−muθteξ+ηγt+eξ+muβt−αteξ+ηγt+eξ+qeξ+mueηq−1e2ξσ12t+mu2e2ησ22t2≤eω2te1+  qξ+mu1+qe1+  qη·ω2′t+rueξ+mueηqeξ−rlkue2+  qξ−mu1+qμte1+  qη+qe1+  qξσ12t+mu1+qe1+  qησ22t2≤eω2te1+  qξ·ω2′t+mu1+qe1+  qη·ω2′t+rueξ+mueηqeξ−rlkue2+  qξ−mu1+qμte1+  qη+qe1+  qξσ12t+mu1+qe1+  qησ22t2≤eω2ue1+  qξK+rueξ+mueηqeξ−rlkue2+  qξ+σ12uq2e1+  qξ−λ2mu1+qe1+  qη=−eω2url2kue2+  qξ+λ2mu1+q2e1+  qη+N1,whereN1=supξ,  η∈R2−rl/2kue2+  qξ−λ2mu1+q/  2e1+  qη+rueξ+mueηqeξ+σ12uq/  2e1+  qξ+e1+  qξ  Keω2u.Hence(62)LV=LV1+LV2≤−M1λ1+M1h1eξ+η−eω2url2kue2+  qξ+λ2mu1+q2e1+  qη+N1.To confirm the condition(B.2) of Lemma 14, we choose small enough ϵ>0 such that(63)0<ϵ≤14h1minλ1,λ2mu1+qeω2uM1,rleω2ukuM1,(64)−M1λ1+ki′+1≤eω2uminrl4kuϵ,λ2mu1+q4ϵi=1,2.Consider the bounded open subset(65)Uϵ=ξ,  η:ξ≤lnϵ−1,η≤lnϵ−1,ξ,  η∈R2.Denote(66)Uϵ1=ξ,  η∈R2:−∞<ξ≤ln  ϵ,Uϵ2=ξ,η∈R2:−∞<η≤ln  ϵ,Uϵ3=ξ,  η∈R2:ξ≥lnϵ−1,Uϵ4=ξ,  η∈R2:η≥lnϵ−1.Clearly,Uc=Uϵ1∪Uϵ2∪Uϵ3∪Uϵ4. Next, we will prove LV≤−1 for any ξ,  η∈Uc, respectively.Case 1. In domain Uϵ1, owing to −∞<ξ≤lnϵeξ+η≤ϵeη≤ϵ1+e1+  qη, we have(67)LV≤−M1λ14+−M1λ14+ϵM1h1+−λ2mu1+q4eω2u+ϵM1h1e1+  qη−rl4kueω2ue2+  qξ+−M1λ12−rl4kueω2ue2+  qξ−λ2mu1+q4eω2ue1+  qη+N1,≤−M1λ14+−M1λ14+ϵM1h1+−λ2mu1+q4eω2u+ϵM1h1e1+  qη−rl4kueω2ue2+  qξ+−M1λ12+supξ,  η∈R2−rl4kueω2ue2+  qξ−λ2mu1+q4eω2ue1+  qη+N1,sinceM1λ1/4≥1, combined with the definition of M1 and Eq. (63), we get(68)LV≤−M1λ14−rl4kueω2ue2+  qξ≤−M1λ14≤−1.Case 2. Similarly, for any ξ,  η∈Uϵ2, owing to eξ+η≤ϵeξ≤ϵ1+e2+qξ, we have(69)LV≤−M1λ14+−M1λ14+ϵM1h1+−rl4kueω2u+ϵM1h1e2+  qξ−λ2mu1+q4eω2ue1+  qη+−M1λ12−rl4kueω2ue2+qξ−λ2mu1+q4eω2ue1+  qη+N1,≤−M1λ14+−M1λ14+ϵM1h1+rl4kueω2u+ϵM1h1e2+  qξ−λ2mu1+q4eω2ue1+  qη+−M1λ12+supξ,  η∈R2−rl4kueω2ue2+  qξ−λ2mu1+q4eω2ue1+  qη+N1,by Eq. (63), we also obtain(70)LV≤−M1λ14−λ2mu1+qeω2u4e1+  qη≤−M1λ14≤−1.Case 3. In domain Uϵ3 we have(71)LV≤−M1λ1+M1h1eξ+η−eω2url2kue2+  qξ+λ2mu1+q2e1+  qη+N1≤−M1λ1−rleω2u4kuϵ+k1′,wherek1′=supξ,  η∈R2−rleω2u/4kue2+  qξ−λ2mu1+qeω2u/  2e1+  qη+M1h1eξ+η+N1, which gives LV≤−1 in domain Uϵ3 in view of Eq. (64).Case 4. In domain Uϵ4 we have(72)LV≤−M1λ1+M1h1eξ+η−eω2url2kue2+  qξ+λ2mu1+q2e1+  qη+N1≤−M1λ1−λ2mu1+qeω2u4ϵ+k2′,wherek2′=supξ,  η∈R2−λ2mu1+qeω2u/  4e1+  qη−rleω2u/2kue2+  qξ+M1h1eξ+η+N1, which gives LV≤−1 in domain Uϵ4 in view of Eq. (64). Consequently, we can deduce that(73)LV≤−1,∀ξ,  η∈Uc.Thus, the condition (B.2) of Lemma14 is satisfied. By Lemma 14, the system (3) exists a periodic Markov process. This completes the proof.Considering the corresponding deterministic system of system (3)(74)dp=prt−rtktp−αtzγt+pdt,dz=zβtpγt+p−μt−θtpγt+pdt.We can obtain the following result according to Theorem15.Corollary 16.
Assume
r
t
,
k
t
,
α
t
,
β
t
,
γ
t
,
μ
t
,
θ
t
are all positive
θ
1-periodic continuous functions and ∫0θ1βl/ru1+γu/klrs−μs+θsds>0, then for any initial valuep0,z0∈R+2, the system (74) exists a positiveθ1-periodic solution.


## 7. Numerical Simulation

In this section, we present some numerical simulations to illustrate our theoretical results obtained in previous sections. To this end, based on the method mentioned in Higham [50], we consider the following discretization equations(75)pk+1=pk+pkr−rkpk−αzkγ+pkΔt+σ1pkΔtξk+σ122pkξk2−  1Δt,zk+1=zk+zkβpkγ+pk−θpkγ+pk−  μΔt+σ2zkΔtηk+σ222zkηk2−  1Δt,whereξk and ηkk=1,2...,  n are independent Gaussian random variables which follow standard normal distribution N0,1.We choose the parameters by(76)r=1,k=100,β=0.22,α=1.2,μ=0.1,θ=0.072,γ=25,with the initial valuep0,z0=80,35.Next we use different values ofσ1,σ2 to see the effect of the noise strength on the dynamics of the system (2).(i)
Fixσ2=0.5, and let σ1 vary to see the effect of noise on the dynamics of the system (2). We firstly take σ1=0.5, then σ12<2r,σ22>2β−2μ, by Theorem 7 (i) the population pt will be persistent in mean and population zt will be extinct, see Figure 1; Now let σ1 vary from σ1=0.5 to σ1=1.5, then σ12>2r,σ22>2β−2μ, by Theorem 7 (ii) the population pt and zt will be extinct, see Figure 2. From Figure 1 and Figure 2, we can see that under small noise, system (2) preserves some stability and large noise may lead to the extinction of population.(ii)
Letσ1=0.03,σ2=0.0015, then λ≜β/r1+γ/  kr−σ12/  2−μ−σ12/2−θ=0.00391>0, by Theorem 11, the system (2) exists a single ergodic stationary distribution, see Figure 3. This implies that system (2) still keep some stability and exhibits oscillation around the interior equilibrium of the corresponding deterministic system.Figure 1
The solutions of system (2) and its corresponding deterministic system with σ1=0.5,σ2=0.5 and other parameters are taken as Eq. (76).Figure 2
The solutions of system (2) and its corresponding deterministic system with σ1=1.5,σ2=0.5 and other parameters are taken as Eq. (76).Figure 3
(a) and (b) are the solutions of system (2) and their corresponding deterministic system, respectively; (c) and (d) are the density function diagrams of pt,zt, where σ1=0.03,σ2=0.0015 and other parameters are taken as Eq. (12).
(a)
(b)
(c)
(d)Next, we choose the 2π-periodic functions(77)rt=1+0.3sint,kt=23+0.1sint,βt=1.2+0.05sint,αt=1.6+0.1sint,μt=0.1+0.01sint,γt=20+0.1sint,with the initial valuep0,z0=8,33. Then we use different values of θ,σ1,σ2 to explore how the environmental noise and toxicity affect the dynamics of system (3).(iii)
We firstly adoptθt=0.25+0.005sint and let σ1t=σ2t=0. It is easy to verify that the condition of Corollary 16 holds and deterministic nonautonomous system (74) also exists a periodic solution, shown in the Figure 4(a), and the corresponding solutions of system (74) shown in the Figure 4(b). Next, we increase strengths of environmental forcing to σ1t=0.2+0.5sint,σ2t=0.12+0.1sint, then the conditions (H.1), (H.2) hold. By Theorem 15, we know the system (3) exists a positive 2 π-periodic Markov process, shown in the Figure 5(a), and the corresponding solutions of system (3) shown in the Figure 5(b).(iv)
Now, we selectθt=1.2+0.005sint,σ1t=σ2t=0, then the population zt will be extinct, shown in Figure 6. By comparing Figure 4 with Figure 6, one can observe that the stronger toxicity will result in the extinction of zooplankton population.Figure 4
The solutions and phase diagram of the deterministic periodic system (74) with σ1t=0,σ2t=0 and other parameters are taken as Eq. (77).
(a)
(b)Figure 5
The solutions and phase diagram of the stochastic periodic system (3) with σ1t=0.2+0.5sint,σ2t=0.12+0.1sint and other parameters are taken as Eq. (77).
(a)
(b)Figure 6
The solutions and phase diagram of the deterministic periodic system (72) with θt=1.2+0.005sint,σ1t=0,σ2t=0 and other parameters are taken as Eq. (13).
(a)
(b)


## 8. Conclusion

The field studies and laboratory studies point out that the toxic substance plays one of the important role on the growth of the zooplankton population and have a great impact on phytoplankton–zooplankton intersection. It is well established that large number of phytoplankton species produce toxic, such asGymnodinium breve, Gymodinium catenatum, Pyrodinium bahamense, Pfiesteria piscicida, Chrysochromulina polylepis, Prymnesium patelliferum, Noctiluca scintillans and so on. Zooplankton species such as Paracalanus will be greatly affected by harmful phytoplankton species [1–5]. In this paper, we propose autonomous and nonautonomous stochastic toxin-producing phytoplankton–zooplankton model with Holling II functional response and made an attempt to reveal the effects of toxic intensity and environmental fluctuations on the plankton ecosystems. For the autonomous system, we establish sufficient conditions for the existence of the globally positive solution, and obtain sufficient conditions for the solution of population extinction and persistence in the mean. Furthermore, by using Khasminskii’s method and technique of Lyapunov functions, we also prove that there exists a single ergodic stationary distribution, what is important is that Lyapunov function does not depend on existence and stability of equilibrium. For the nonautonomous periodic system, we mainly study the existence of positive periodic solution.The theoretical results and numerical simulations show that the population would extinct as the enhancing of noise and/or toxic intensity, while the reductive speed of the population would slow down as the weakening of noise and/or toxic intensity, and the population would be persistent. That is, toxic intensity and environmental fluctuations have great influence on plankton ecosystems. From Theorems7 and 5.3, we know that the survival of plankton can be significantly affected by the white noise densities σii=1,2 and the release rate of toxic substances θ. That is, if σ1 is sufficiently large, the phytoplankton pt may suffer the danger of extinction. On the contrary, if the intensities of environment noise σii=1,2 and the rate of release of toxic substances θ are not too large such that λ>0, then the system (2) exists a single ergodic stationary distribution. The obtained results also implies that the TPP may provide a possible biological strategy to control the plankton ecosystems.



---

*Source: 4627571-2020-01-20.xml*