# Impact of Cell Size Effect on Nutrient-Phytoplankton Dynamics

Authors: Tiancai Liao; Hengguo Yu; Chuanjun Dai; Min Zhao
Journal: Complexity
 (2019)

---

## Abstract

In this paper, a nutrient-phytoplankton model, which is described by a system of ordinary differential equations incorporating the effect of cell size, and its corresponding stochastic differential equation version are studied analytically and numerically. A key advantage of considering cell size effect is that it can more accurately reveal the intrinsic law of interaction between nutrient and phytoplankton. The main purpose of this paper is to research how cell size affects the nutrient-phytoplankton dynamics within the deterministic and stochastic environments. Mathematically, we show that the existence and stability of the equilibria in the deterministic model can be determined by cell size: the smaller or larger cell size can lead to the disappearance of the positive equilibrium, but the boundary equilibrium always exists and is globally asymptotically stable; the intermediate cell size is capable to drive the positive equilibrium to appear and be globally asymptotically stable, whereas the boundary equilibrium becomes unstable. In the case of the stochastic model, the stochastic dynamics including the stochastic extinction, persistence in the mean, and the existence of ergodic stationary distribution is found to be largely dependent on cell size and noise intensity. Ecologically, via numerical simulations, it is found that the smaller cell size or larger cell size can result in the extinction of phytoplankton, which is similar to the effect of larger random environmental fluctuations on the phytoplankton. More interestingly, it is discovered that the intermediate cell size is the optimal size for promoting the growth of phytoplankton, but increasing appropriately the cell size can rapidly reduce phytoplankton density and nutrient concentrations at the same time, which provides a possible strategy for biological control of algal blooms.


---

## Body

## 1. Introduction

Phytoplankton blooms, which can negatively affect the aquatic ecosystems, human health, marine fisheries, and local economy, are growing in frequency, magnitude, and duration globally in recent years [1, 2]. For example, in 2007, serious harmful algal blooms broke out in Lake Taihu, resulting in about 2 million people in Wuxi city short of drinking water for more than a week [3]. In 2011, Lake Erie experienced a record-breaking harmful algal bloom, with a peak intensity more than three times the previously observed algal blooms [4]. These have stimulated a number of scholars to study the dynamics of phytoplankton growth by many different ways, in order to explore the possible mechanisms underlying the occurrence or termination of these blooms. Nevertheless, such a mechanism in response to the phenomena of blooms is still under investigation due to the diversity and complexity of influencing factors that can affect the growth of phytoplankton in the real aquatic ecological environments. Hence, the research in seeking for some key factors affecting the growth mechanisms of phytoplankton is currently of great interest.However, in the real aquatic environments, the growth of phytoplankton is generally influenced by many biotic and abiotic factors, such as light [5], cell size [6], climate [7], grazer [8], carbon dioxide [9], nutrient [10], and temperature [11], which make it difficult to determine a clear mechanism of phytoplankton blooms only through experimental studies. Actually, many ecologists, biologists, and biomathematicians increasingly realize that a mathematical model is a powerful tool for exploring biological and physical processes on the dynamic mechanisms of phytoplankton growth in relation to different factors qualitatively and quantitatively [12, 13], as the research results can help us to find out the key factors that may induce the blooms of phytoplankton but are difficult to predict in the experimental analysis, to answer that what the growth mechanism of phytoplankton is, to predict possibly when the phytoplankton blooms will occur, and to determine the optimal strategy for possible control of phytoplankton blooms [14–25]. The application of mathematical models in other research fields, such as investigating other predator-prey dynamics or infectious disease dynamics, can be found in [26–37].In 1949, Riley et al. [38] first used the mathematical model to study the nutrient-plankton dynamics, which leads to the formulation of a growing number of mathematical models to describe the nutrient-phytoplankton dynamics or nutrient-plankton dynamics, and many dynamic mechanisms of phytoplankton growth response to various factors have been revealed [14, 15, 22, 39–47]. For example, Chen et al. [43] showed that the proper control of the ratio for nitrogen and phosphorus can more effectively control and eliminateblue-green algae blooms. Pal et al. [44] indicated that the toxin produced by phytoplankton plays a crucial role in the termination of planktonic blooms when the nutrient concentration is very high. Dai et al. [39] signified that time delay not only induces instability of a positive equilibrium but also promotes the formation of patchiness. Chatterjee et al. [41] observed that the nutrient input rate specially caused by artificial eutrophication has a great influence on controlling the planktonic blooms and maintaining stability around the coexistence equilibrium. These excellent results have made a great contribution to the research processes related to the possible mechanisms for the formation or termination of phytoplankton blooms, which in turn indicate that the modeling studies are an alternative, effective, and feasible method to investigate some important factors that may reveal the nature of the blooms of phytoplankton occurring frequently in various water bodies around the world.Most of the existing mathematical modeling studies on the nutrient-phytoplankton dynamics or nutrient-plankton dynamics usually assume that the nutrient uptake rate, phytoplankton sinking rate, phytoplankton growth rate, and so on are independent of cell size [12, 14, 15, 22, 39–45], which have, in part, been considered to be unrealistic because the factor of cell size is capable to significantly affect the dynamic mechanisms of phytoplankton growth [48]. In fact, cell size is a master functional trait that virtually affects every aspect of phytoplankton biology at the cellular, population, and community levels [49]. Hence, the cell size of phytoplankton not only defines their metabolic activity, growth rates, and numerical abundance but also strongly affects their contributions to biogeochemical cycles via size-dependent sinking and influences community structure and dynamics via size-dependent species interactions [50–53]. A recent remarkable experimental work, in this research direction, was carried out by Marañón et al. [6], where they further determined experimentally that the growth rate, metabolic rate, and nutrient uptake rate of phytoplankton are size-dependent by using cultures of 22 species of marine phytoplankton from five phyla, ranging from 0.1 to 106μm3 in cell volume. Furthermore, by performing an in situ test of Raven’s prediction that there is a reversal of the relationship between cell size and maximum achievable growth rate in unicellular algae at the low end of size classes, Bec et al. [54] found that the biomass specific production and growth rates are similar in both small and large cells but peak at intermediate cell sizes. Nevertheless, despite these facts, there is little theoretical explanation on how cell size affects the growth mechanism of phytoplankton. In recent years, a review of relevant literature shows that several plankton models taking into account the effect of cell size have been reported [55–57]. More specifically, Pu et al. [56] indicated that the coevolutionary dynamics of plankton is closely related to cell size. Zhao et al. [57] found that the cell size can significantly affect the growth and reproduction of phytoplankton. The cell size plays an important role in the dynamics of interacting phytoplankton, which has been clearly demonstrated by these recent studies.In spite of the importance of cell size, a comprehensive understanding towards the dynamic mechanism for phytoplankton growth in relation to cell size is still lacking. Hence, naturally, some questions arise. How cell size affects the nutrient concentration and phytoplankton density? Can the cell size affect the existence and stability of the equilibria? In order to find out the answers to these questions, motivated by these works [12, 55, 56, 58], in the present study, we propose and investigate a nutrient-phytoplankton model, incorporating the effect of cell size, as follows:(1)dNdt=I−μxQxNPa+N−mN,dPdt=μxNPa+N−sxP−dP,subject to the initial conditions N0=N0≥0 and P0=P0≥0, where Nt and Pt are the concentration of nutrient and the density of phytoplankton at time t, respectively. Here, we assume that I is the input rate of nutrient from the environments, m is the removal rate of nutrient, and d is maximum mortality rate of phytoplankton. Let x be the cell size of phytoplankton; we suppose that μx=x/a1x2+a2x+a3 [55] is the maximum specific growth rate of phytoplankton as a function of cell size, and Qx=βx3 [55] signifies the nutrient quota which is proportional to the cube of cell size, where aii=1,2,3 are positive empirical constants and β denotes the phytoplankton nutrient quota coefficient. The nutrient uptake rate of phytoplankton is assumed to depend on μx and Qx. That is, the term μxQxNP/a+N is taken as a Michaelis–Menten function that may provide a more realistic reflection of nutrient uptake dynamics [12], leading to a decrease in nutrient concentration, where a is the half-saturation constant of nutrient. Accordingly, the term μxNP/a+N is the conversion efficiency of nutrient uptake by phytoplankton, which means that phytoplankton absorbs nutrient to reproduce offspring and increase the number of phytoplankton. Suppose that the sinking rate of phytoplankton is proportional to the square of cell size, i.e., sx=αx2 [55], where α is a constant whose value is affected by the density of the water and the algal cell as well as the viscosity of the water. All the parameters mentioned above are assumed to be positive. From a biological viewpoint, we need to ensure that dP/dt>0, so it is assumed that μx−sx−d>0 be always established by default in the whole paper.On the other hand, in the natural world, the aquatic environment in which phytoplankton lives is always uncertain and random because of the environmental noise disturbances [59], such as unpredictable radiation, light availability, and water temperature variation. Hence, the parameters contained in deterministic models are usually assumed to be invariant constants, which have obvious limitations in modeling the natural ecological systems. Actually, most of the natural phenomena do not strictly follow deterministic laws, but rather oscillate randomly about some average behaviors. A fact further pointed out by May [60] is that the birth rate, carrying capacity, and other parameters involved the model should exhibit random fluctuation to a greater or lesser extent because of the environmental noise effects. Consequently, using stochastic differential equation models may reveal the dynamic mechanisms of phytoplankton growth more accurately compared to their deterministic counterparts. For these reasons, some authors recently introduced environmental noise fluctuations into the aquatic plankton models to study noise influences on the interplay and growth of phytoplankton [61–64]. For example, Yu et al. [61] indicated that environmental fluctuations play a key role in the termination of algal blooms. Camara et al. [63] suggested that stochastic environmental constraints have positive and negative effects on the life of Daphnia and algae populations. Obviously, the environmental fluctuations have a vital role in the dynamics of phytoplankton growth. In other words, stochastic environmental disturbance effects should be considered when mathematical models are used to study and model the ecological systems in nature.Consequently, we introduce white noise perturbations into the deterministic model (1). In the existing literature, there are many different ways to incorporate the noise fluctuations into the ecological systems. In this study, by following the method in [65], we assume that stochastic environmental fluctuations mainly affect the growth of phytoplankton μx. In this way, μx changes to a random variable μx¯, and μx¯=μx+δB˙t, where Bt is a standard Brownian motion defined on a complete probability space Ω,ℱ,ℱtt≥0,ℙ, B˙t indicates the white noise, and δ represents the intensity of the white noise. Therefore, by replacing μx in the deterministic model (1) with μx+δB˙t, i.e.,(2)−muxQxNPa+N⟶−μxQxNPa+N−δμxQxNPa+NB˙t,μxNPa+N⟶μxNPa+N+δμxNPa+NB˙t,and then, model (1) becomes(3)dN=I−μxQxNPa+N−mNdt−δμxQxNPa+NdBt,dP=μxNPa+N−sxP−dPdt+δμxNPa+NdBt.In the plankton ecology, the persistence and extinction of phytoplankton, which determine whether the phytoplankton is survival or not in the future, are two important topics. However, in the existing literature, the issues on how cell size affects the stochastic extinction and persistence of phytoplankton and how cell size affects the existence of ergodic stationary distribution of phytoplankton that denotes the weak stability of phytoplankton in a stochastic sense, and so on remain largely unanswered. In this paper, we did some works in this research area by studying the stochastic dynamics of model (3).The main purpose of this paper is to investigate how cell size affects the nutrient-phytoplankton dynamics within the deterministic and stochastic environments by trying to answer the questions proposed in this paper.In order to facilitate the mathematical analysis for model (1) and model (3) in the following paper, we define(4)Xx=μx−sx−dd+sx>0,and we write Xx as X for simplicity and convenience. The rest of this article is organized as follows: in Section 2, we explore the dynamics of model (1). Section 3 is devoted to studying the dynamics of model (3). In Section 4, we carry out the numerical simulations to verify the analytical results. Section 5 focuses on giving a brief discussion and the summary of the main results.


## 2. Dynamics of the Deterministic Model (1)

In this section, we mainly focus on investigating the positivity and boundedness of the solutions, as well as studying the existence and stability of the possible equilibria in model (1).


### 2.1. Positivity and Boundedness of the Solutions

Now, we first present the positivity of the solutions.Lemma 1.
For any initial valueN0,P0=N0,P0>0, all the solutions of model (1) are positive invariant.Proof.
In model (1), we have Nt>0 and Pt>0 for any t∈0,B, where B is any positive real number. Suppose this is not true, then there exists Bτ∈0,B such that, for any t∈0,Bτ, Nt>0, Pt>0, and either NBτ=0 or PBτ=0. Based on the nutrient and phytoplankton equations in model (1), we can obtain(5)Nt=N0exp∫0tINs−μxQxPsa+Ns−mds,Pt=P0exp∫0tμxNsa+Ns−d−sxds.
SinceNt,Pt are defined and continuous in 0,Bτ, there exists ϱ≥0 such that for all t∈0,Bτ, (6)Nt=N0exp∫0tINs−μxQxPsa+Ns−mds≥N0exp−Bτϱ,Pt=P0exp∫0tμxNsa+Ns−d−sxds≥P0exp−Bτϱ.
Obviously, ifN0,P0=N0,P0>0 and let t⟶Bτ, we have NBτ≥N0exp−Bτϱ>0 and PBτ≥P0exp−Bτϱ>0, which contradict the fact that either NBτ=0 or PBτ=0. Hence, for all t∈0,Bτ, Nt>0 and Pt>0. This completes the proof.
Next, we study the boundedness of the solutions. LetWt=Nt+QxPt, then we have(7)dWdt=dNdt+dPdt=I−μxQxNPa+N−mN+μxQxNPa+N−dQxP−QxsxP=I−mN−d+sxQxP≤I−eW,where e=minm,d+sx. Thus, we have dW/dt+eW≤I. By using the theory of differential inequality [66], for all t≥T≥0, we have(8)0≤Wt≤Ie−Ie−WTexp−et−T.
Hence,(9)limt⟶∞supNt+QxPt≤Ie.
So, we can obtain the theorem as follows.Theorem 1.
All the solutions of model (1) that start in R+2 are uniformly bounded.


### 2.2. Existence and Stability of Equilibria

This section is devoted to investigating the existence and stability of the possible equilibria in model (1). Now, we first consider the existence of possible equilibria in model (1).From the nutrient and phytoplankton equations in model (1), by performing a simple computation, we have that model (1) possesses two equilibria as follows:(a)
The boundary equilibriumE1=I/m,0 always exists(b)
The positive equilibriumE∗=a/X,IX−ma/Qxμx−sx−d, provided X>ma/INext, we deal with the stability of the equilibria in model (1). Through direct calculations, the Jacobian matrix of model (1) at equilibrium EN,P is(10)JE=−μxQxaPa+N2−m−μxQxNa+NμxaPa+N2Nμx−d+sxa+Na+N.Obviously, the variational matrix of model (1) at the boundary equilibrium E1 is(11)JE1=−m−IμxQxI+ma0IX−mad+sxma+I,and the two eigenvalues of JE1 are λ1=−m<0 and λ2=IX−mad+sx/ma+I. Hence, the stability of E1 depends on the sign of λ2. That is, if X<ma/I, E1 is locally asymptotically stable and is unstable if X>ma/I.The variational matrix of model (1) at the positive equilibrium E∗ is given by(12)JE∗=J11J12J21J22=−XXI−maaX+1−m−μxQxX+1XXI−maaX+1Qx0,where J11<0, J12<0, and J21>0. Then, it is not difficult to derive the characteristic equation of the positive equilibrium E∗ which is λ2−trJE∗λ+detJE∗=0, where trJE∗=J11<0 and detJE∗=−J21J12>0. Based on the Routh–Hurwitz criterion, the positive equilibrium E∗ in model (1) is locally asymptotically stable when it exists.Hence, we have the following result.Theorem 2.
For model (1), if X<ma/I, there exists the only boundary equilibrium E1, which is locally asymptotically stable; if X>ma/I, the positive equilibrium E∗ appears and is locally asymptotically stable, but the boundary equilibrium E1 becomes unstable.In what follows, we provide some results concerning the global stability of the equilibria.Theorem 3.
IfX<ma/I holds, the boundary equilibrium of model (1) is globally asymptotically stable.Proof.
LettingE1=I/M,0=N1,P1 and defining fy=μxQxy/a+y with y∈R, we consider the following function:(13)V1N,P=∫N1Nfs−fNfsds+kP,where k is a positive constant that will be determined in the paper later.
The derivative ofV1 with respect to the time t along the solutions of model (1) is(14)dV1dt=fN−fN1fNdNdt+kdPdt=fN−fN1fN−mN−N1−fNP+kQxfNP−sx+dkP=−mfN−fN1N−N1fN+kQx−1fNP+fN1−d+sxkP.
Now, we obtaink=Qx>0, and then, the above equation becomes(15)dV1dt=−mfN−fN1N−N1fN+fN1−d+sxQxP.
Notice thatfy>0 is an increasing function for any y>0, and in view of Lemma 1, if fN1−d+sxQx<0, that is, X<ma/I, one can get dV1/dt≤0 and dV1/dt=0 if and only if N,P=N1,P1. Hence, Lyapunov–LaSalle’s invariance principle implies the global asymptotic stability of E1. This completes the proof.Theorem 4.
IfX>ma/I, the positive equilibrium E∗ in model (1) is globally asymptotically stable.Proof.
We consider such a function:(16)V2N,P=∫N∗Nfs−fN∗fsds+Qx∫P∗Ps−P∗sds,where the function fyy∈R is the same as the definition above.
The derivative ofV2 with respect to the time t along the solutions of model (1) is(17)dV2dt=fN−fN1fNdNdt+QxP−P∗PdPdt=fN−fN∗fN−P∗fN−fN∗−fNP−P∗−mN−N∗+P−P∗fN−fN∗=−P∗fN−fN∗2fN−mN−N∗fN−fN∗fN.
A discussion similar to the proof of Theorem1 shows that dV2/dt≤0 and dV2/dt=0 if and only if N,P=N∗,P∗. Hence, the properties of Lyapunov functional are satisfied, which indicates that the positive equilibrium E∗ is globally asymptotically stable when it exists. Considering the local asymptotical stability of the positive equilibrium, we obtain Theorem 4. This ends the proof.Remark 1.
Theorems2, 3, and 4 indicate that the local stability and global stability of the E∗ (or E1) are equivalent, but the existence and stability of positive equilibrium E∗ and the local or global stability of boundary equilibrium E1 are opposite.


## 2.1. Positivity and Boundedness of the Solutions

Now, we first present the positivity of the solutions.Lemma 1.
For any initial valueN0,P0=N0,P0>0, all the solutions of model (1) are positive invariant.Proof.
In model (1), we have Nt>0 and Pt>0 for any t∈0,B, where B is any positive real number. Suppose this is not true, then there exists Bτ∈0,B such that, for any t∈0,Bτ, Nt>0, Pt>0, and either NBτ=0 or PBτ=0. Based on the nutrient and phytoplankton equations in model (1), we can obtain(5)Nt=N0exp∫0tINs−μxQxPsa+Ns−mds,Pt=P0exp∫0tμxNsa+Ns−d−sxds.
SinceNt,Pt are defined and continuous in 0,Bτ, there exists ϱ≥0 such that for all t∈0,Bτ, (6)Nt=N0exp∫0tINs−μxQxPsa+Ns−mds≥N0exp−Bτϱ,Pt=P0exp∫0tμxNsa+Ns−d−sxds≥P0exp−Bτϱ.
Obviously, ifN0,P0=N0,P0>0 and let t⟶Bτ, we have NBτ≥N0exp−Bτϱ>0 and PBτ≥P0exp−Bτϱ>0, which contradict the fact that either NBτ=0 or PBτ=0. Hence, for all t∈0,Bτ, Nt>0 and Pt>0. This completes the proof.
Next, we study the boundedness of the solutions. LetWt=Nt+QxPt, then we have(7)dWdt=dNdt+dPdt=I−μxQxNPa+N−mN+μxQxNPa+N−dQxP−QxsxP=I−mN−d+sxQxP≤I−eW,where e=minm,d+sx. Thus, we have dW/dt+eW≤I. By using the theory of differential inequality [66], for all t≥T≥0, we have(8)0≤Wt≤Ie−Ie−WTexp−et−T.
Hence,(9)limt⟶∞supNt+QxPt≤Ie.
So, we can obtain the theorem as follows.Theorem 1.
All the solutions of model (1) that start in R+2 are uniformly bounded.


## 2.2. Existence and Stability of Equilibria

This section is devoted to investigating the existence and stability of the possible equilibria in model (1). Now, we first consider the existence of possible equilibria in model (1).From the nutrient and phytoplankton equations in model (1), by performing a simple computation, we have that model (1) possesses two equilibria as follows:(a)
The boundary equilibriumE1=I/m,0 always exists(b)
The positive equilibriumE∗=a/X,IX−ma/Qxμx−sx−d, provided X>ma/INext, we deal with the stability of the equilibria in model (1). Through direct calculations, the Jacobian matrix of model (1) at equilibrium EN,P is(10)JE=−μxQxaPa+N2−m−μxQxNa+NμxaPa+N2Nμx−d+sxa+Na+N.Obviously, the variational matrix of model (1) at the boundary equilibrium E1 is(11)JE1=−m−IμxQxI+ma0IX−mad+sxma+I,and the two eigenvalues of JE1 are λ1=−m<0 and λ2=IX−mad+sx/ma+I. Hence, the stability of E1 depends on the sign of λ2. That is, if X<ma/I, E1 is locally asymptotically stable and is unstable if X>ma/I.The variational matrix of model (1) at the positive equilibrium E∗ is given by(12)JE∗=J11J12J21J22=−XXI−maaX+1−m−μxQxX+1XXI−maaX+1Qx0,where J11<0, J12<0, and J21>0. Then, it is not difficult to derive the characteristic equation of the positive equilibrium E∗ which is λ2−trJE∗λ+detJE∗=0, where trJE∗=J11<0 and detJE∗=−J21J12>0. Based on the Routh–Hurwitz criterion, the positive equilibrium E∗ in model (1) is locally asymptotically stable when it exists.Hence, we have the following result.Theorem 2.
For model (1), if X<ma/I, there exists the only boundary equilibrium E1, which is locally asymptotically stable; if X>ma/I, the positive equilibrium E∗ appears and is locally asymptotically stable, but the boundary equilibrium E1 becomes unstable.In what follows, we provide some results concerning the global stability of the equilibria.Theorem 3.
IfX<ma/I holds, the boundary equilibrium of model (1) is globally asymptotically stable.Proof.
LettingE1=I/M,0=N1,P1 and defining fy=μxQxy/a+y with y∈R, we consider the following function:(13)V1N,P=∫N1Nfs−fNfsds+kP,where k is a positive constant that will be determined in the paper later.
The derivative ofV1 with respect to the time t along the solutions of model (1) is(14)dV1dt=fN−fN1fNdNdt+kdPdt=fN−fN1fN−mN−N1−fNP+kQxfNP−sx+dkP=−mfN−fN1N−N1fN+kQx−1fNP+fN1−d+sxkP.
Now, we obtaink=Qx>0, and then, the above equation becomes(15)dV1dt=−mfN−fN1N−N1fN+fN1−d+sxQxP.
Notice thatfy>0 is an increasing function for any y>0, and in view of Lemma 1, if fN1−d+sxQx<0, that is, X<ma/I, one can get dV1/dt≤0 and dV1/dt=0 if and only if N,P=N1,P1. Hence, Lyapunov–LaSalle’s invariance principle implies the global asymptotic stability of E1. This completes the proof.Theorem 4.
IfX>ma/I, the positive equilibrium E∗ in model (1) is globally asymptotically stable.Proof.
We consider such a function:(16)V2N,P=∫N∗Nfs−fN∗fsds+Qx∫P∗Ps−P∗sds,where the function fyy∈R is the same as the definition above.
The derivative ofV2 with respect to the time t along the solutions of model (1) is(17)dV2dt=fN−fN1fNdNdt+QxP−P∗PdPdt=fN−fN∗fN−P∗fN−fN∗−fNP−P∗−mN−N∗+P−P∗fN−fN∗=−P∗fN−fN∗2fN−mN−N∗fN−fN∗fN.
A discussion similar to the proof of Theorem1 shows that dV2/dt≤0 and dV2/dt=0 if and only if N,P=N∗,P∗. Hence, the properties of Lyapunov functional are satisfied, which indicates that the positive equilibrium E∗ is globally asymptotically stable when it exists. Considering the local asymptotical stability of the positive equilibrium, we obtain Theorem 4. This ends the proof.Remark 1.
Theorems2, 3, and 4 indicate that the local stability and global stability of the E∗ (or E1) are equivalent, but the existence and stability of positive equilibrium E∗ and the local or global stability of boundary equilibrium E1 are opposite.


## 3. Dynamics of the Stochastic Model (3)

In this section, we study the existence and uniqueness of the solution, stochastic extinction, persistence in the mean, and a unique ergodic stationary distribution of model (3). First of all, we define a bounded set Θ as follows:(18)Θ=N,P∈R+2:N+QxP≤Ie,which will be used for the proofs of the following results.


### 3.1. Existence and Uniqueness of the Positive Solution

In model (3), the coefficients do not satisfy the linear growth condition, though they are locally Lipschitz continuous. Hence, the solution of model (3) may explore at a finite time. In this section, following the research [65], we prove that model (3) has a unique global positive solution.Theorem 5.
For any initial valueN0,P0∈Θ, there is a unique positive solution Nt,Pt of model (3) on t≥0, and the solution will remain in R+2 with probability one.Proof.
LetN0,P0∈Θ, based on the nutrient and phytoplankton equations in model (3) and by virtue of the positivity of N and P, we can get(19)dN+QxPdt=I−mN−sx+dQxP≤I−eN+QxP,then limn⟶∞supN+QxP≤I/e, that is,(20)limt⟶∞supNt≤Ie,limt⟶∞supPt≤Ieβx3.
So, for anyc∈0,t, we can obtain(21)Nc,Pc∈0,Iea.s.
Consequently, the setΘ is a positivity invariant to model (3).
Since the coefficients of model (3) are locally Lipschitz continuous, for any given initial value N0,P0, there is a unique local solution Nt,Pt on t∈0,τe, where τe is the explosion time [67]. To show the solution is global, we need to show that τe=∞.
Letn0>0 be sufficiently large for initial values N0 and P0 lying within the interval 1/n0,n0. For each integer n>n0, considering the stopping times,(22)τn=inft∈0,τn,Nt∉1n,norPt∉1n,n,and we set infΦ=∞ (Φ denotes the empty set). Obviously, τn is increasing as n⟶∞. Let τ∞=limn⟶∞τn; hence, τ∞≤τea.s.. Next, we only need to show τ∞=∞. If this is false, there is a pair of constants T>0 and ε∈0,1 such that(23)ℙτn≤T>ε,n≥n1.
Define aC2− function V: R+2⟶R+ by(24)VN,P=−lneNI−lnβx3ePI.
UsingIto^’s formula [68], we have(25)dV=−IN+μxQxPa+N+m+δ2μx2Qx2P22a+N2−μxNa+N+sx+d+δ2μx2N22a+N2dt+δμxPQx−Na+NdB≤μxQxPa+m+δ2μx2Qx2P22a2+sx+d+δ2μx2N22a2dt+δμxPQx−Na+NdB≤μxIae+m+sx+d+δ2μx2I2a2e2dt+δμxPQx−Na+NdB.
Therefore,(26)dVNt,Pt≤φdt+δμxPQx−Na+NdB,where φ=μxI/ae+m+sx+d+δ2μx2I2/a2e2. Then,(27)∫0τn∧TdVNt,Pt≤∫0τn∧Tφdt+∫0τn∧TδμxPQx−Na+NdB,where τn∧T=minτn,T. Taking the expectation of the above inequality, we get(28)EVNτn∧T,Pτn∧T≤VN0,P0+φEτn∧T≤VN0,P0+φT.
LetΩn=τn≤T for n≥n1, by means of (23), then we have ℙΩn≥ε. Note that, for every ϑ∈Ωn, there is at least one of Nτn,ϑ and Pτn,ϑ that equals either 1/n or n. Therefore, VNτn,ϑ,Pτn,ϑ is not less than(29)ϕ=min−lnenI−lnβex3nI,−lnenI−lnβex3In,−lnenI−lnβex3nI,−lnenI−lnβex3nI.
It then follows from (28) that(30)VN0,P0+φT≥E1ΩnϑVNt,Pt≥εϕ,where 1Ωnϑ is the indicator function of Ωn. Letting n⟶∞, we have(31)∞>VN0,P0+φT=∞,which is a contradiction. Hence, we have τ∞=∞. That is, the solution of model (3) will not explore at a finite time with probability one. This completes the proof.


### 3.2. Stochastic Persistence and Extinction of Phytoplankton

From the proof of Theorem5, we know that Θ is the positive invariant set of model (3). In this section, we always assume that the initial value N0,P0∈Θ and first provide the result concerning the persistence of phytoplankton. For the sake of simplicity, we define(32)Λ=δ2μ2xI22ea+I2+sx+d.Then, we have the following result.Theorem 6.
For any given initialN0,P0∈Θ, if one of the following conditions holds:(33)iδ2<2eI+eaImμx,X>2mI+ea2μx−sx−d2eIμxI+ea−I2δ2mμ2x,iiδ2>2eI+eaImμx,X<2mI+ea2μx−sx−d2eIμxI+ea−I2δ2μ2xm,then the solution of model (3) obeys (34)limt⟶∞inf1t∫0tPsds≥eIμx−Λmea+IeQxμxsx+d>0a.s.,that is, the phytoplankton in model (3) is persistent in the mean almost surely.Proof.
Based on the nutrient and phytoplankton equations in model (3), we have(35)dN+θP=I−mN−sx+dQxPdt.
Integrating equation (35) from 0 to t on the both sides, we get(36)Nt+QxPt−N0+QxP0=tI−m∫0tNsds−sx+dQx∫0tPsds.
Note thatNt+QxPt=N0+QxP0 for all t>0, and therefore, we adopt(37)1t∫0tNsds=Im−sx+dQxm1t∫0tPsds.
ByIto^’s formula, we have(38)dlnPt=−δ2μ2x2Na+N2+μxNa+N−sx−ddt+δμxNa+NdB.
LetfN=N/a+N, then it is easy to get dfN/dN>0. Hence, fN is a strictly increasing function for any N∈R+. In light of N≤I/e, we can obtain(39)−δ2μ2x2Na+N2+μxNa+N−sx−d≥−δ2μ2xI22ea+I2+μxeNea+I−sx−d.
Integrating equation (38) from 0 to t on the both sides, we have(40)lnPt−lnP0≥−δ2μ2xI22ea+I2+sx+dt+μxeea+I∫0tNsds+κt,where(41)κt=∫0tδμxNsa+NsdBs.
Dividingt on both sides of equation (40) and together with equation (37), we get(42)lnPtt≥eIμxea+Im−Λ−eQxμxsx+dea+Im1t∫0tPsds+κt+lnP0t.
Notice thatκt is a local continuous martingale with κ0=0 and(43)limt⟶∞κ,κtt≤δμxIae+I2<∞.
By the strong law of large number [69], we can achieve(44)limt⟶∞κt+lnP0t=0a.s.
By means of Lemma 4 in [70], if Λ<eIμx/ea+Im, that is,(45)2mI+ea2μx−sx−d<X2eIμxI+ea−I2δ2μ2xm,which is equivalent to(46)X>2mI+ea2μx−sx−d2eIμxI+ea−I2δ2mμ2xprovidedδ2<2eI+eaImμx,orX<2mI+ea2μx−sx−d2eIμxI+ea−I2δ2μ2xmprovidedδ2>2eI+eaImμx,then it follows equation (42) that(47)limt⟶∞inf1t∫0tPsds≥eIμx−ΛmI+eaeQxμxsx+da.s.
This completes the proof.
Next, we study the stochastic extinction of phytoplankton in model (3).Theorem 7.
For any given initialN0,P0∈Θ, if one of the following conditions holds:(48)iδ2<2I+eaIμx,X<2I+ea2μx−sx−d2IμxI+ea−δ2μ2xI2,iiδ2>2I+eaIμx,X>2I+ea2μx−sx−d2IμxI+ea−δ2I2μ2x,then the solution of model (3) obeys(49)limt⟶∞suplnPtt≤IμxI+ea−Λ<0a.s.,orlimt⟶∞suplnPtt≤12δ2−sx−d<0a.s.,provided δ2>maxI+ea/Iμx,1/2sx+d, which implies that the phytoplankton in model (3) is stochastic extinction with probability one.Proof.
Define the following quadratic function(50)gv=−δ2μ2x2v2+μxv−sx−d,if the symmetric axis 1/δ2μx of the equation gv is larger than I/ae+I, that is, δ2≤I+ea/μxI, we can find that gv takes its maximum value gmax on the interval 0,I/ae+I at v=I/ae+I, where(51)gmax=IμxI+ea−Λ.
Note thatfN=N/a+N is a strictly increasing function for N∈0,I/e when δ2≤I+ea/Iμx; therefore, we have(52)−δ2μ2x2Na+N2+μxNa+N−sx−d≤gmax.
Integrating equation (38) from 0 to t and dividing by t on the both sides, we get(53)lnPtt≤IμxI+ea−Λt+κt+lnP0t.
Taking the limit of equation (53) and considering (44), we have(54)limt⟶∞suplnPtt≤IμxI+ea−Λ<0a.s.,provided Λ>Iμx/I+ea, i.e.,(55)2I+ea2μx−sx−d>X2IμxI+ea−δ2I2μ2x,which equals to(56)δ2<2I+eaIμx,X<2I+ea2μx−sx−d2IμxI+ea−δ2μ2xI2,orδ2>2I+eaIμx,X>2I+ea2μx−sx−d2IμxI+ea−δ2I2μ2x.
On the other hand, if the symmetric axis1/δ2μx of the equation gv is in the interval 0,I/ea+I, i.e., δ2≥I+ea/Iμx, then the equation gv reaches its maximum value 1/2δ2−sx−d at v=2/δ2μx. By performing a similar proof processes as above, we have (57)limt⟶∞suplnPtt≤12δ2−sx−d<0a.s.,provided δ2>maxI+ea/Iμx,1/2sx+d. This completes the proof.Remark 2.
Theorem7 indicates that the phytoplankton in model (3) goes to extinction with probability one under some parameter conditions. In this case, it is natural to ask the question how will the nutrient concentration in model (3) change when the phytoplankton is extinct? In fact, when ℙlimt⟶∞Pt=0=1, then from the model (3), we can obtain the following equation:(58)dNt=I−mNtdt,by direct computation, we have limt⟶∞Nt=I/m. Hence, the nutrient concentration will be a positive constant. That is, in this situation, the noise intensity and cell size will have no impacts on the variation of nutrient, and the dynamics of Nt in model (3) will be similar to the dynamics of E1 in the deterministic model (1).


### 3.3. Existence of Ergodic Stationary Distribution

In this section, we shall establish sufficient conditions for the existence and uniqueness of an ergodic stationary distribution for model (3).LetZt be a regular time-homogeneous Markov process in Rd described by the stochastic differential equation:(59)dZt=fZtdt+∑l=1kglZdBltdBt.The diffusion matrix of the processZt is defined as follows:(60)Az=aijz,aijz=∑l=1kglizgljz.Lemma 2 (see [71]).
The Markov processZt has a unique ergodic stationary distribution μ⋅ if there exists a bounded open domain D⊂Rd with regular boundary Γ, having the following properties:B.1: there is a positive number M such that ∑i,j=1dai,jzξiξj≥Mξ2, z∈D and ξ∈RdB.2: there exists a nonnegative C2−function V such that LV is negative for any Rd\DTheorem 8.
IfX<e2I+ea2/e2I+ea2I+δ2I4, then for any given initial N0,P0∈Θ, model (3) has a unique stationary distribution μ⋅ and it has the ergodic property.Proof.
To prove Theorem8, we only need to validate condition B.1 and B.2 in Lemma 2. We first prove the condition B.1. The diffusion matrix of model (3) is given by(61)∑i,j=12aijN,Pξiξj=−δμxQxNPξ1a+NδμxNPξ2a+N−δμxQxNPξ1a+NδμxNPξ2a+N=δ2μ2xQ2xN2P2a+N2ξ12+δ2μ2xN2P2a+N2ξ22≥M0ξ,for anyN,P∈D¯σ⊂Θ,ξ=ξ1,ξ2∈R+2,where M0=minN,P∈D¯σδ2μ2xQ2xN2P2/a+N2,δ2μ2xN2P2/a+N2, and then the condition B.1 in Lemma 2 is satisfied. Next, we prove the condition B.2.
Let(62)M=2λmax2,H,where H=supN,P∈Θ−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN and λ=1/X−I+δ2I4/e2I+ea2>0 provided X<e2I+ea2/e2I+ea2I+δ2I4. Obviously, Mλ/4≥1. Then, we define a C2− function VN,P: Θ⟶R as follows:(63)VN,P=MN+P+N22μ2x+Q2xμ2xP22+lnPμx−sx−1+N+QxP2,and it is easy to see that(64)limn⟶+∞,N,P∈Θ∖DnVN,P=+∞,where Dn=1/n,n×1/n,n. Furthermore, VN,P is continuous function, and so, VN,P must have a minimum point N0,P0 in the interior of Θ. Then, we define nonnegative C2− function VN,P:Θ⟶R in the following form:(65)VN,P=MN+P+N22μ2x+Q2xμ2xP22+lnPμx−sx−1+N+QxP2−VN0,P0:=MV1N,P+V2N,P+V3N,P+V4N,P,where V1N,P=N+P, V2N,P=N2/2μ2x+Q2x/μ2xP2/2, V3=lnP/μx−sx−1, and V4N,P=N+QxP2−VN0,P0.
Making using ofIto^’s formula, we have(66)LV1N,P=I−μxQxNPa+N−mN+μxNPa+N−sx+dP≤I−mN−sx+dP+eμxNPI+ea≤I+eμxNPI+ea,(67)LV2N,P=1μ2xIN−μxQxPN2a+N−mN2+δ2μ2xQ2xN2P22a+N2+Qxμx2μxNP2a+N−sx+dP2+δ2μ2xN2P22a+N2≤INμ2x−mN2μ2x+QxNPμxI+ea−sx+dQ2xP2μ2x+δ2I4e2I+ea2,(68)LV3N,P=1μx−sx−1μxNa+N−sx+d−δ2μ2xN22a+N2≤μxNaμx−sx−d−1X−δ2e2N2μ2x2ea+I2μx−sx−d,(69)LV4N,P=2NI−mN2−sx+dNP+QxIP−Qxsx+dP2+δ2μ2xQxQx+1N2P2a+N2≤−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN.
Hence, we get(70)LVN,P=MLV1N,P+LV2N,P+LV3N,P+V4N,P≤−Mλ+MIμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea−mμ2x+δ2e2μ2x2ea+I2μx−sx−dN2−sx+dQ2xμ2xP2+−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN.
Define a bounded closed set(71)Dε=N,P∈Θ:ε≤N≤1ε,ε≤P≤1ε,where 0<ε<1 is a sufficiently small number. We choose ε sufficiently small such that(72)0<ε<λ4I/μ2x+μx/aμx−sx−d+Qx+eμ2x/μxI+ea,(73)0<ε<Q2xsx+dI+eaμxQx+eμ2x,(74)0<ε<λμxI+ea2Qx+eμ2x,(75)0<ε<2mI+ea2μx−sx−d+δ2e2μ4x2μxI+eaμx−sx−dQx+e2μx,(76)−Mλ2+H¯≤−1,(77)−Mλ−mMμ2x+Mδ2e2μ2x2ea+I2μx−sx−d1ε2+H¯∗,(78)−Mλ−Msx+dQ2xμ2x1ε2+H¯∗,where(79)H¯=supN,P∈Θ−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+MIμ2x+μxaμx−sx−dN,H¯∗=supN,P∈Θ−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+Iμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea.
For convenience, we divideΘ\Dε into four domains:(80)Dε1=N,P∈Θ:0<N<ε,Dε2=N,P∈Θ:0<P<ε,Dε3=N,P∈Θ:N>1ε,Dε4=N,P∈Θ:P>1ε.
Clearly,DεC=∪i4Dεi. Now, we prove that LVN,P≤1 on DεC, which is equivalent to prove it on Dε1, Dε2, Dε3, and Dε4, respectively.Case 1.
WhenN,P∈Dε1, since NP≤εP≤ε1+P2, we have(81)LVN,P≤−Mλ4+MIεμ2x+Mμxεaμx−sx−d+Qx+eμ2xMεμxI+ea−Mλ4−mMμ2x+Mδ2e2μ2x2ea+I2μx−sx−dN2+Qx+eμ2xMεμxI+ea−Msx+dQ2xμ2xP2+−Mλ2+H,which together with conditions (72) and (73), as well as the definition of M in (62), yields(82)LVN,P≤−Mλ4≤−1.Case 2.
WhenN,P∈Dε2, since NP≤εN≤ε1+N2, we get(83)LVN,P≤−Mλ2+Qx+eμ2xMεμxI+ea−MQ2xsx+dμ2xP2+MQx+eμ2xεμxI+ea−Mmμ2x+δ2e2μ2x2I+ea2μx−sx−dN2+−Mλ2+H¯,which leads to LVN,P≤−Mλ/2+H¯≤1 in this domain, in the view of (74)–(76).Case 3.
WhenN,P∈Dε3, we have(84)LVN,P≤−Mλ+M−mμ2x+δ2e2μ2x2ea+I2μx−sx−dN2−sx+dQ2xμ2xP2+−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+Iμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea≤−Mλ−mMμ2x+Mδ2e2μ2x2ea+I2μx−sx−dN2+H¯∗,from (77), and we have LVN,P≤−1 on N,P∈Dε3.Case 4.
WhenN,P∈Dε4, we have(85)LVN,P≤−Mλ+M−mμ2x+δ2e2μ2x2ea+I2μx−sx−dN2−sx+dQ2xμ2xP2+−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+Iμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea≤−Mλ−Msx+dQ2xμ2xP2+H¯∗,in view of (78), and we get LVN,P≤−1. This ends the proof.


## 3.1. Existence and Uniqueness of the Positive Solution

In model (3), the coefficients do not satisfy the linear growth condition, though they are locally Lipschitz continuous. Hence, the solution of model (3) may explore at a finite time. In this section, following the research [65], we prove that model (3) has a unique global positive solution.Theorem 5.
For any initial valueN0,P0∈Θ, there is a unique positive solution Nt,Pt of model (3) on t≥0, and the solution will remain in R+2 with probability one.Proof.
LetN0,P0∈Θ, based on the nutrient and phytoplankton equations in model (3) and by virtue of the positivity of N and P, we can get(19)dN+QxPdt=I−mN−sx+dQxP≤I−eN+QxP,then limn⟶∞supN+QxP≤I/e, that is,(20)limt⟶∞supNt≤Ie,limt⟶∞supPt≤Ieβx3.
So, for anyc∈0,t, we can obtain(21)Nc,Pc∈0,Iea.s.
Consequently, the setΘ is a positivity invariant to model (3).
Since the coefficients of model (3) are locally Lipschitz continuous, for any given initial value N0,P0, there is a unique local solution Nt,Pt on t∈0,τe, where τe is the explosion time [67]. To show the solution is global, we need to show that τe=∞.
Letn0>0 be sufficiently large for initial values N0 and P0 lying within the interval 1/n0,n0. For each integer n>n0, considering the stopping times,(22)τn=inft∈0,τn,Nt∉1n,norPt∉1n,n,and we set infΦ=∞ (Φ denotes the empty set). Obviously, τn is increasing as n⟶∞. Let τ∞=limn⟶∞τn; hence, τ∞≤τea.s.. Next, we only need to show τ∞=∞. If this is false, there is a pair of constants T>0 and ε∈0,1 such that(23)ℙτn≤T>ε,n≥n1.
Define aC2− function V: R+2⟶R+ by(24)VN,P=−lneNI−lnβx3ePI.
UsingIto^’s formula [68], we have(25)dV=−IN+μxQxPa+N+m+δ2μx2Qx2P22a+N2−μxNa+N+sx+d+δ2μx2N22a+N2dt+δμxPQx−Na+NdB≤μxQxPa+m+δ2μx2Qx2P22a2+sx+d+δ2μx2N22a2dt+δμxPQx−Na+NdB≤μxIae+m+sx+d+δ2μx2I2a2e2dt+δμxPQx−Na+NdB.
Therefore,(26)dVNt,Pt≤φdt+δμxPQx−Na+NdB,where φ=μxI/ae+m+sx+d+δ2μx2I2/a2e2. Then,(27)∫0τn∧TdVNt,Pt≤∫0τn∧Tφdt+∫0τn∧TδμxPQx−Na+NdB,where τn∧T=minτn,T. Taking the expectation of the above inequality, we get(28)EVNτn∧T,Pτn∧T≤VN0,P0+φEτn∧T≤VN0,P0+φT.
LetΩn=τn≤T for n≥n1, by means of (23), then we have ℙΩn≥ε. Note that, for every ϑ∈Ωn, there is at least one of Nτn,ϑ and Pτn,ϑ that equals either 1/n or n. Therefore, VNτn,ϑ,Pτn,ϑ is not less than(29)ϕ=min−lnenI−lnβex3nI,−lnenI−lnβex3In,−lnenI−lnβex3nI,−lnenI−lnβex3nI.
It then follows from (28) that(30)VN0,P0+φT≥E1ΩnϑVNt,Pt≥εϕ,where 1Ωnϑ is the indicator function of Ωn. Letting n⟶∞, we have(31)∞>VN0,P0+φT=∞,which is a contradiction. Hence, we have τ∞=∞. That is, the solution of model (3) will not explore at a finite time with probability one. This completes the proof.


## 3.2. Stochastic Persistence and Extinction of Phytoplankton

From the proof of Theorem5, we know that Θ is the positive invariant set of model (3). In this section, we always assume that the initial value N0,P0∈Θ and first provide the result concerning the persistence of phytoplankton. For the sake of simplicity, we define(32)Λ=δ2μ2xI22ea+I2+sx+d.Then, we have the following result.Theorem 6.
For any given initialN0,P0∈Θ, if one of the following conditions holds:(33)iδ2<2eI+eaImμx,X>2mI+ea2μx−sx−d2eIμxI+ea−I2δ2mμ2x,iiδ2>2eI+eaImμx,X<2mI+ea2μx−sx−d2eIμxI+ea−I2δ2μ2xm,then the solution of model (3) obeys (34)limt⟶∞inf1t∫0tPsds≥eIμx−Λmea+IeQxμxsx+d>0a.s.,that is, the phytoplankton in model (3) is persistent in the mean almost surely.Proof.
Based on the nutrient and phytoplankton equations in model (3), we have(35)dN+θP=I−mN−sx+dQxPdt.
Integrating equation (35) from 0 to t on the both sides, we get(36)Nt+QxPt−N0+QxP0=tI−m∫0tNsds−sx+dQx∫0tPsds.
Note thatNt+QxPt=N0+QxP0 for all t>0, and therefore, we adopt(37)1t∫0tNsds=Im−sx+dQxm1t∫0tPsds.
ByIto^’s formula, we have(38)dlnPt=−δ2μ2x2Na+N2+μxNa+N−sx−ddt+δμxNa+NdB.
LetfN=N/a+N, then it is easy to get dfN/dN>0. Hence, fN is a strictly increasing function for any N∈R+. In light of N≤I/e, we can obtain(39)−δ2μ2x2Na+N2+μxNa+N−sx−d≥−δ2μ2xI22ea+I2+μxeNea+I−sx−d.
Integrating equation (38) from 0 to t on the both sides, we have(40)lnPt−lnP0≥−δ2μ2xI22ea+I2+sx+dt+μxeea+I∫0tNsds+κt,where(41)κt=∫0tδμxNsa+NsdBs.
Dividingt on both sides of equation (40) and together with equation (37), we get(42)lnPtt≥eIμxea+Im−Λ−eQxμxsx+dea+Im1t∫0tPsds+κt+lnP0t.
Notice thatκt is a local continuous martingale with κ0=0 and(43)limt⟶∞κ,κtt≤δμxIae+I2<∞.
By the strong law of large number [69], we can achieve(44)limt⟶∞κt+lnP0t=0a.s.
By means of Lemma 4 in [70], if Λ<eIμx/ea+Im, that is,(45)2mI+ea2μx−sx−d<X2eIμxI+ea−I2δ2μ2xm,which is equivalent to(46)X>2mI+ea2μx−sx−d2eIμxI+ea−I2δ2mμ2xprovidedδ2<2eI+eaImμx,orX<2mI+ea2μx−sx−d2eIμxI+ea−I2δ2μ2xmprovidedδ2>2eI+eaImμx,then it follows equation (42) that(47)limt⟶∞inf1t∫0tPsds≥eIμx−ΛmI+eaeQxμxsx+da.s.
This completes the proof.
Next, we study the stochastic extinction of phytoplankton in model (3).Theorem 7.
For any given initialN0,P0∈Θ, if one of the following conditions holds:(48)iδ2<2I+eaIμx,X<2I+ea2μx−sx−d2IμxI+ea−δ2μ2xI2,iiδ2>2I+eaIμx,X>2I+ea2μx−sx−d2IμxI+ea−δ2I2μ2x,then the solution of model (3) obeys(49)limt⟶∞suplnPtt≤IμxI+ea−Λ<0a.s.,orlimt⟶∞suplnPtt≤12δ2−sx−d<0a.s.,provided δ2>maxI+ea/Iμx,1/2sx+d, which implies that the phytoplankton in model (3) is stochastic extinction with probability one.Proof.
Define the following quadratic function(50)gv=−δ2μ2x2v2+μxv−sx−d,if the symmetric axis 1/δ2μx of the equation gv is larger than I/ae+I, that is, δ2≤I+ea/μxI, we can find that gv takes its maximum value gmax on the interval 0,I/ae+I at v=I/ae+I, where(51)gmax=IμxI+ea−Λ.
Note thatfN=N/a+N is a strictly increasing function for N∈0,I/e when δ2≤I+ea/Iμx; therefore, we have(52)−δ2μ2x2Na+N2+μxNa+N−sx−d≤gmax.
Integrating equation (38) from 0 to t and dividing by t on the both sides, we get(53)lnPtt≤IμxI+ea−Λt+κt+lnP0t.
Taking the limit of equation (53) and considering (44), we have(54)limt⟶∞suplnPtt≤IμxI+ea−Λ<0a.s.,provided Λ>Iμx/I+ea, i.e.,(55)2I+ea2μx−sx−d>X2IμxI+ea−δ2I2μ2x,which equals to(56)δ2<2I+eaIμx,X<2I+ea2μx−sx−d2IμxI+ea−δ2μ2xI2,orδ2>2I+eaIμx,X>2I+ea2μx−sx−d2IμxI+ea−δ2I2μ2x.
On the other hand, if the symmetric axis1/δ2μx of the equation gv is in the interval 0,I/ea+I, i.e., δ2≥I+ea/Iμx, then the equation gv reaches its maximum value 1/2δ2−sx−d at v=2/δ2μx. By performing a similar proof processes as above, we have (57)limt⟶∞suplnPtt≤12δ2−sx−d<0a.s.,provided δ2>maxI+ea/Iμx,1/2sx+d. This completes the proof.Remark 2.
Theorem7 indicates that the phytoplankton in model (3) goes to extinction with probability one under some parameter conditions. In this case, it is natural to ask the question how will the nutrient concentration in model (3) change when the phytoplankton is extinct? In fact, when ℙlimt⟶∞Pt=0=1, then from the model (3), we can obtain the following equation:(58)dNt=I−mNtdt,by direct computation, we have limt⟶∞Nt=I/m. Hence, the nutrient concentration will be a positive constant. That is, in this situation, the noise intensity and cell size will have no impacts on the variation of nutrient, and the dynamics of Nt in model (3) will be similar to the dynamics of E1 in the deterministic model (1).


## 3.3. Existence of Ergodic Stationary Distribution

In this section, we shall establish sufficient conditions for the existence and uniqueness of an ergodic stationary distribution for model (3).LetZt be a regular time-homogeneous Markov process in Rd described by the stochastic differential equation:(59)dZt=fZtdt+∑l=1kglZdBltdBt.The diffusion matrix of the processZt is defined as follows:(60)Az=aijz,aijz=∑l=1kglizgljz.Lemma 2 (see [71]).
The Markov processZt has a unique ergodic stationary distribution μ⋅ if there exists a bounded open domain D⊂Rd with regular boundary Γ, having the following properties:B.1: there is a positive number M such that ∑i,j=1dai,jzξiξj≥Mξ2, z∈D and ξ∈RdB.2: there exists a nonnegative C2−function V such that LV is negative for any Rd\DTheorem 8.
IfX<e2I+ea2/e2I+ea2I+δ2I4, then for any given initial N0,P0∈Θ, model (3) has a unique stationary distribution μ⋅ and it has the ergodic property.Proof.
To prove Theorem8, we only need to validate condition B.1 and B.2 in Lemma 2. We first prove the condition B.1. The diffusion matrix of model (3) is given by(61)∑i,j=12aijN,Pξiξj=−δμxQxNPξ1a+NδμxNPξ2a+N−δμxQxNPξ1a+NδμxNPξ2a+N=δ2μ2xQ2xN2P2a+N2ξ12+δ2μ2xN2P2a+N2ξ22≥M0ξ,for anyN,P∈D¯σ⊂Θ,ξ=ξ1,ξ2∈R+2,where M0=minN,P∈D¯σδ2μ2xQ2xN2P2/a+N2,δ2μ2xN2P2/a+N2, and then the condition B.1 in Lemma 2 is satisfied. Next, we prove the condition B.2.
Let(62)M=2λmax2,H,where H=supN,P∈Θ−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN and λ=1/X−I+δ2I4/e2I+ea2>0 provided X<e2I+ea2/e2I+ea2I+δ2I4. Obviously, Mλ/4≥1. Then, we define a C2− function VN,P: Θ⟶R as follows:(63)VN,P=MN+P+N22μ2x+Q2xμ2xP22+lnPμx−sx−1+N+QxP2,and it is easy to see that(64)limn⟶+∞,N,P∈Θ∖DnVN,P=+∞,where Dn=1/n,n×1/n,n. Furthermore, VN,P is continuous function, and so, VN,P must have a minimum point N0,P0 in the interior of Θ. Then, we define nonnegative C2− function VN,P:Θ⟶R in the following form:(65)VN,P=MN+P+N22μ2x+Q2xμ2xP22+lnPμx−sx−1+N+QxP2−VN0,P0:=MV1N,P+V2N,P+V3N,P+V4N,P,where V1N,P=N+P, V2N,P=N2/2μ2x+Q2x/μ2xP2/2, V3=lnP/μx−sx−1, and V4N,P=N+QxP2−VN0,P0.
Making using ofIto^’s formula, we have(66)LV1N,P=I−μxQxNPa+N−mN+μxNPa+N−sx+dP≤I−mN−sx+dP+eμxNPI+ea≤I+eμxNPI+ea,(67)LV2N,P=1μ2xIN−μxQxPN2a+N−mN2+δ2μ2xQ2xN2P22a+N2+Qxμx2μxNP2a+N−sx+dP2+δ2μ2xN2P22a+N2≤INμ2x−mN2μ2x+QxNPμxI+ea−sx+dQ2xP2μ2x+δ2I4e2I+ea2,(68)LV3N,P=1μx−sx−1μxNa+N−sx+d−δ2μ2xN22a+N2≤μxNaμx−sx−d−1X−δ2e2N2μ2x2ea+I2μx−sx−d,(69)LV4N,P=2NI−mN2−sx+dNP+QxIP−Qxsx+dP2+δ2μ2xQxQx+1N2P2a+N2≤−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN.
Hence, we get(70)LVN,P=MLV1N,P+LV2N,P+LV3N,P+V4N,P≤−Mλ+MIμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea−mμ2x+δ2e2μ2x2ea+I2μx−sx−dN2−sx+dQ2xμ2xP2+−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN.
Define a bounded closed set(71)Dε=N,P∈Θ:ε≤N≤1ε,ε≤P≤1ε,where 0<ε<1 is a sufficiently small number. We choose ε sufficiently small such that(72)0<ε<λ4I/μ2x+μx/aμx−sx−d+Qx+eμ2x/μxI+ea,(73)0<ε<Q2xsx+dI+eaμxQx+eμ2x,(74)0<ε<λμxI+ea2Qx+eμ2x,(75)0<ε<2mI+ea2μx−sx−d+δ2e2μ4x2μxI+eaμx−sx−dQx+e2μx,(76)−Mλ2+H¯≤−1,(77)−Mλ−mMμ2x+Mδ2e2μ2x2ea+I2μx−sx−d1ε2+H¯∗,(78)−Mλ−Msx+dQ2xμ2x1ε2+H¯∗,where(79)H¯=supN,P∈Θ−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+MIμ2x+μxaμx−sx−dN,H¯∗=supN,P∈Θ−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+Iμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea.
For convenience, we divideΘ\Dε into four domains:(80)Dε1=N,P∈Θ:0<N<ε,Dε2=N,P∈Θ:0<P<ε,Dε3=N,P∈Θ:N>1ε,Dε4=N,P∈Θ:P>1ε.
Clearly,DεC=∪i4Dεi. Now, we prove that LVN,P≤1 on DεC, which is equivalent to prove it on Dε1, Dε2, Dε3, and Dε4, respectively.Case 1.
WhenN,P∈Dε1, since NP≤εP≤ε1+P2, we have(81)LVN,P≤−Mλ4+MIεμ2x+Mμxεaμx−sx−d+Qx+eμ2xMεμxI+ea−Mλ4−mMμ2x+Mδ2e2μ2x2ea+I2μx−sx−dN2+Qx+eμ2xMεμxI+ea−Msx+dQ2xμ2xP2+−Mλ2+H,which together with conditions (72) and (73), as well as the definition of M in (62), yields(82)LVN,P≤−Mλ4≤−1.Case 2.
WhenN,P∈Dε2, since NP≤εN≤ε1+N2, we get(83)LVN,P≤−Mλ2+Qx+eμ2xMεμxI+ea−MQ2xsx+dμ2xP2+MQx+eμ2xεμxI+ea−Mmμ2x+δ2e2μ2x2I+ea2μx−sx−dN2+−Mλ2+H¯,which leads to LVN,P≤−Mλ/2+H¯≤1 in this domain, in the view of (74)–(76).Case 3.
WhenN,P∈Dε3, we have(84)LVN,P≤−Mλ+M−mμ2x+δ2e2μ2x2ea+I2μx−sx−dN2−sx+dQ2xμ2xP2+−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+Iμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea≤−Mλ−mMμ2x+Mδ2e2μ2x2ea+I2μx−sx−dN2+H¯∗,from (77), and we have LVN,P≤−1 on N,P∈Dε3.Case 4.
WhenN,P∈Dε4, we have(85)LVN,P≤−Mλ+M−mμ2x+δ2e2μ2x2ea+I2μx−sx−dN2−sx+dQ2xμ2xP2+−2mN2+2IN−2Qxsx+dP2+2QxI+e−1δ2Iμ2xQx+1P−2sx+Qxm+dPN+Iμ2x+μxaμx−sx−dN+Qx+eμ2xNPμxI+ea≤−Mλ−Msx+dQ2xμ2xP2+H¯∗,in view of (78), and we get LVN,P≤−1. This ends the proof.


## 4. Numerical Results

In this section, we study numerically the impact of the cell size effect on the nutrient-phytoplankton dynamics in the deterministic model (1) as well as the stochastic model (3). The parameter values are taken as shown in Table 1. Unless stated otherwise, we will use them for our simulations. Additionally, in the following numerical simulation, we randomly select the initial value satisfying Θ.Table 1
Parameter values in numerical simulations for model (1) and model (3).
Parameters
Unit
Values
I: the nutrient input rate
μ
mol
/
Lay
10 [72]
α: the phytoplankton sinking rate coefficient
1
/
day
/
μ
m
2
0.1 [55]
β: the phytoplankton nutrient quota coefficient
μ
mol nutrient
/
cell
/
μ
m
−
b
1
100 (estimate)
m: the removal rate of nutrient
1
/
day
0.85 (estimate)
d: the death rate of phytoplankton
1
/
day
0.1 [55]
a: the half-saturation constant
μ
mol
/
L
100 (estimate)
a
1: a positive empirical constant
1
/
μ
m
2
0.02 [55]
a
2: a positive empirical constant
1
/
μ
m
0.02 [55]
a
3: a positive empirical constant
Unitless
0.08 [55]


### 4.1. Impact of Cell Size Factor in the Nutrient-Phytoplankton Dynamics of Model (1)

In this section, we study the impact of the cell size effect on the concentration of nutrient and the density of phytoplankton, as well as the stability of the equilibria in model (1).On the basis of previous theoretical analysis, letfx=X−ma/I, and substituting (4) and the required parameters into fx, we get fx<0, i.e.,X<ma/I if(86)H1:=x∈R∣0<x<0.07807216867or2.981219088<x,and fx>0, i.e.,X>ma/I if(87)H2:=x∈R∣0.07807216867<x<2.981219088.Next, we first show the impact of cell size on the nutrient concentration and phytoplankton density.


#### 4.1.1. Impact of the Cell Size Effect on the Nutrient Concentration and Phytoplankton Density

Based on the analytical expressions for the positive equilibrium level of nutrient and phytoplankton, one can know that the nutrient concentration and phytoplankton density are closely related to cell sizex. However, let x∈H1; the increase of cell size x does not affect the nutrient concentration, which is clearly demonstrated in Figures 1(a) and 1(c). Furthermore, in this case, Figures 1(d) and 1(f) show that the density of phytoplankton is negative or zero. That is, the phytoplankton in model (1) cannot survive when the phytoplankton cell is smaller or larger (x∈H1), which may be the main reason why the nutrient concentration does not change with the increase of cell size.Figure 1
(a–c) Bifurcation diagram of nutrient in model (1) with respect to x. (d–f) Bifurcation diagram of phytoplankton in model (1) with respect to x. (a) and (d) for x∈0.01,0.07807216867. (b) and (e) for x∈0.07807216867,2.981219088. (c) and (f) for x∈2.981219088,7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 2
(a–c) Bifurcation diagram of nutrient in model (1) with respect to x∈H2. (b–f) Bifurcation diagram of phytoplankton in model (1) with respect to x∈H2. (a) and (d) for x∈0.07807216867,0.1. (b) and (e) for x∈0.1,0.6. (c) and (f) for x∈0.6,2.981219088.
(a)
(b)
(c)
(d)
(e)
(f)However, as the cell size increases to the domain ofH2, the nutrient concentration is a nonmonotonic concave function with respect to cell size x (see Figure 1(b)) and phytoplankton density is a nonmonotonic convex function of cell size x (see Figure 1(e)). Nevertheless, by further analysis, we find that(1)
Ifx∈0.07807216867,0.1, the nutrient concentration is a monotonic decreasing function with respect to cell size x (see Figure 2(a)), but phytoplankton density is a convex-line increasing function of cell size x (see Figure 2(d)), which signifies that phytoplankton can reproduce rapidly and abundantly and absorb a large amount of nutrients during its growth, resulting in a rapid increase in phytoplankton density and a corresponding rapid decrease in nutrient concentration. In fact, from Figure 1(e) or Figures 2(d) and 2(e), we can know that the peak of phytoplankton density appears at x=0.1, and so, x=0.1 is the optimal size for promoting the rapid reproduction of the phytoplankton population.
(2)
Ifx∈0.1,0.6, both the nutrient concentration and phytoplankton density are concave-line decreasing function of cell size x, which are shown in Figures 2(b) and 2(e), respectively. As the density of phytoplankton reaches its peak at x=0.1, when the cell size of phytoplankton exceeds its optimal size, the phytoplankton density decreases with the increase of cell size, but its density is still very high, resulting in the decrease of nutrient concentration, whereas the reduction rate of the nutrient concentration is decreasing compared with that in Figure 2(a).
(3)
Ifx∈0.6,2.981219088, the nutrient concentration is a concave-line increasing function with respect to x (see Figure 2(c)), but phytoplankton density is a concave-line decreasing function of cell size, as is shown in Figure 2(f). Obviously, in this case, with the increase of cell size, the density of phytoplankton decreases dramatically with the rapid increase of nutrient concentration.


#### 4.1.2. Impact of Cell Size Effect on the Stability of the Equilibria

The analysis indicates that the plane is divided into two regions (i) and (ii) byI−x parameters, as displayed in Figure 3(a). However, when the phytoplankton cell x is smaller or larger (x∈H1), there exists a unique boundary equilibrium E1 in the area (i), which is globally asymptotically stable (see Theorem 3). With the increase of cell size x∈H2, the positive equilibrium E∗ appears in (ii) and it is globally asymptotically stable, but the E1 becomes unstable (see Theorems 2 and 4). Then, we adopt x=0.05∈H1 and x=0.51∈H2; the corresponding results concerning the global stability of E1=11.76470588,0 and E∗=2.414081849,4.754928396 are shown in (i) and (ii) of Figure 3(b), respectively. It is worth noting that the only difference between (i) and (ii) in Figure 3(b) is that the cell size of phytoplankton is different. Actually, from the Figure 3(a), it is clear that the increase of cell size can lead to the stability changes at the boundary equilibrium E1: stable 0<x<0.07807216867 ⟶ unstable 0.07807216867<x<2.981219088 ⟶ stable x>2.981219088, and the existence of positive equilibrium E∗ changes: nonexistence 0<x<0.07807216867 ⟶ existence 0.07807216867<x<2.981219088 ⟶ nonexistence x>2.981219088, which indicates that the cell size plays an important role in determining the existence and stability of the equilibria in model (1).Figure 3
(a) Bifurcation diagram with respect tox and I, where the white area (i) indicates that there exists the only boundary equilibrium E1 that is globally asymptotically stable, and the pink area (ii) implies that the positive equilibrium E∗ is globally asymptotically stable but the boundary equilibrium E1 is unstable. (b) The global stability of the equilibria in model (1) with different cell size x, in which (i) for the boundary equilibrium E1 with x=0.05 and (ii) for the positive equilibrium E∗ with x=0.51.
(a)
(b)


### 4.2. Impact of Cell Size Factor in the Nutrient-Phytoplankton Dynamics of Model (3)

Now, we turn to study the impact of cell size effect on the stochastic dynamics of model (3) by using the Milstein method mentioned in Higham [73]. In this way, model (3) can be rewritten as the following discretization equation:(88)Nk+1=Nk+I−μxQxNkPka+Nk−mNkΔt−δμxQxNkPka+NkΔtξk+δ22μxQxNkPka+Nkξk2−1Δt,Pk+1=Pk+μxNkPka+Nk−sxPk−dPkΔt+δμxNkPka+NkΔtξk+δ22μxNkPka+Nkξk2−1Δt,where ξkk=1,2,…,n are independent Gaussian random variables N0,1.


#### 4.2.1. Impact of Cell Size and Environmental Noise Effects on the Survival of Phytoplankton

First of all, we fixx=0.51 and vary δ=0.1,1.8 to investigate how environmental noise affects the persistence of model (3).When we takeδ=0.1, by simple computations, it is not difficult to get that the condition (i) in Theorem 6 and the condition in Theorem 8 are satisfied. From Theorems 6 and 8, it can be obtained that model (3) is persistent in the mean and exists a unique stationary distribution, as is shown in Figure 4. This result signifies that the small environmental noise disturbance can make model (3) maintain some stability in the random sense and drive the solutions to be perturbed near the positive equilibrium of the deterministic model (1).Figure 4
(a) and (c) are the solutions of nutrient and phytoplankton in the stochastic model (3) and its corresponding deterministic model (1) with initial value N0,P0=8.5,2.4, respectively. (b) and (d) are the probability density function of nutrient and phytoplankton in model (3), respectively.
(a)
(b)
(c)
(d)When we obtain the noise intensityδ=1.8, then the condition (ii) of Theorem 7 is satisfied. Based on Theorem 7, the phytoplankton goes to extinction with probability one, which implies that the larger environmental noise can result in the extinction of phytoplankton. Figure 5(b) clearly confirms this result. Further, in this case, the nutrient concentration is found to increase to N=N1=11.76470588>N∗=2.414081849 (see Figure 5(a)), which agrees well with Remark 2.Figure 5
The solutions of the stochastic model (3) and its corresponding deterministic model (1) with initial value N0,P0=8.5,2.4. (a) For the nutrient. (b) For the phytoplankton.
(a)
(b)Next, we fixδ=0.1 and change x=0.001,0.15,0.45,0.95 to see the effect of cell size on the survival of phytoplankton in model (3). It is clear from Figure 6 that if x=0.15or0.45⊂H2, the phytoplankton is persistent and becomes extinct when x=0.001or0.98⊂H1, which indicates that the smaller or larger cell size can result in the extinction of phytoplankton (see the black and pink curves) and the intermediate cell size is advantage for the persistence of phytoplankton (see the blue and green curves). Moreover, by comparing the green curve with the blue curve, it can be found that the increase of cell size can reduce the random variation of phytoplankton caused by environmental noise disturbance and decrease the density of phytoplankton.Figure 6
The solution of phytoplankton in model (3) with different cell size x, where the noise intensity is fixed to δ=0.1, and the black curve is for x=0.001, the green curve is for x=0.15, the blue curve is for x=0.45, and the pink curve is for x=0.95.


#### 4.2.2. Impact of Cell Size Effect on the Existence of Ergodic Stationary Distribution

In order to explore how cell size affects the existence of stationary distribution of phytoplankton and nutrient in model (3), we fix the noise intensity δ=0.1 and vary the cell size(89)x=0.01,0.02,0.03,0.04,0.05,0.06,which satisfy the parameter condition in Theorem 8. Hence, there exists a stationary distribution of model (3), as is shown in Figure 7.Figure 7
The effect of cell size on the existence of stationary distribution for model (3) with different cell size x=0.01,0.02,0.03,0.04,0.05,0.06. (a) and (c) are the probability density function of phytoplankton and nutrient in model (3), respectively. (b) and (d) are the boxplot of the solution for phytoplankton and nutrient in model (3) with the corresponding cell size, respectively.
(a)
(b)
(c)
(d)In Figure7(a), one can find that, with the increase of cell size, the height of probability density function of phytoplankton decreases, but that of the probability density function of nutrition increases (see Figure 7(c)), which suggests that the increase of cell size is not conducive to stabilize the nutrient-phytoplankton dynamics in a stochastic sense.In Figure7(b), it is shown that the maximum value of the boxplot for phytoplankton increases, while the maximum value of the boxplot for nutrient remains unchanged (see Figure 7(d)), as the cell size increases. This implies that the probability density function of phytoplankton is positively skewed (from left to right), but the probability density function of nutrition is not shifted as the cell size increases.


## 4.1. Impact of Cell Size Factor in the Nutrient-Phytoplankton Dynamics of Model (1)

In this section, we study the impact of the cell size effect on the concentration of nutrient and the density of phytoplankton, as well as the stability of the equilibria in model (1).On the basis of previous theoretical analysis, letfx=X−ma/I, and substituting (4) and the required parameters into fx, we get fx<0, i.e.,X<ma/I if(86)H1:=x∈R∣0<x<0.07807216867or2.981219088<x,and fx>0, i.e.,X>ma/I if(87)H2:=x∈R∣0.07807216867<x<2.981219088.Next, we first show the impact of cell size on the nutrient concentration and phytoplankton density.


### 4.1.1. Impact of the Cell Size Effect on the Nutrient Concentration and Phytoplankton Density

Based on the analytical expressions for the positive equilibrium level of nutrient and phytoplankton, one can know that the nutrient concentration and phytoplankton density are closely related to cell sizex. However, let x∈H1; the increase of cell size x does not affect the nutrient concentration, which is clearly demonstrated in Figures 1(a) and 1(c). Furthermore, in this case, Figures 1(d) and 1(f) show that the density of phytoplankton is negative or zero. That is, the phytoplankton in model (1) cannot survive when the phytoplankton cell is smaller or larger (x∈H1), which may be the main reason why the nutrient concentration does not change with the increase of cell size.Figure 1
(a–c) Bifurcation diagram of nutrient in model (1) with respect to x. (d–f) Bifurcation diagram of phytoplankton in model (1) with respect to x. (a) and (d) for x∈0.01,0.07807216867. (b) and (e) for x∈0.07807216867,2.981219088. (c) and (f) for x∈2.981219088,7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 2
(a–c) Bifurcation diagram of nutrient in model (1) with respect to x∈H2. (b–f) Bifurcation diagram of phytoplankton in model (1) with respect to x∈H2. (a) and (d) for x∈0.07807216867,0.1. (b) and (e) for x∈0.1,0.6. (c) and (f) for x∈0.6,2.981219088.
(a)
(b)
(c)
(d)
(e)
(f)However, as the cell size increases to the domain ofH2, the nutrient concentration is a nonmonotonic concave function with respect to cell size x (see Figure 1(b)) and phytoplankton density is a nonmonotonic convex function of cell size x (see Figure 1(e)). Nevertheless, by further analysis, we find that(1)
Ifx∈0.07807216867,0.1, the nutrient concentration is a monotonic decreasing function with respect to cell size x (see Figure 2(a)), but phytoplankton density is a convex-line increasing function of cell size x (see Figure 2(d)), which signifies that phytoplankton can reproduce rapidly and abundantly and absorb a large amount of nutrients during its growth, resulting in a rapid increase in phytoplankton density and a corresponding rapid decrease in nutrient concentration. In fact, from Figure 1(e) or Figures 2(d) and 2(e), we can know that the peak of phytoplankton density appears at x=0.1, and so, x=0.1 is the optimal size for promoting the rapid reproduction of the phytoplankton population.
(2)
Ifx∈0.1,0.6, both the nutrient concentration and phytoplankton density are concave-line decreasing function of cell size x, which are shown in Figures 2(b) and 2(e), respectively. As the density of phytoplankton reaches its peak at x=0.1, when the cell size of phytoplankton exceeds its optimal size, the phytoplankton density decreases with the increase of cell size, but its density is still very high, resulting in the decrease of nutrient concentration, whereas the reduction rate of the nutrient concentration is decreasing compared with that in Figure 2(a).
(3)
Ifx∈0.6,2.981219088, the nutrient concentration is a concave-line increasing function with respect to x (see Figure 2(c)), but phytoplankton density is a concave-line decreasing function of cell size, as is shown in Figure 2(f). Obviously, in this case, with the increase of cell size, the density of phytoplankton decreases dramatically with the rapid increase of nutrient concentration.


### 4.1.2. Impact of Cell Size Effect on the Stability of the Equilibria

The analysis indicates that the plane is divided into two regions (i) and (ii) byI−x parameters, as displayed in Figure 3(a). However, when the phytoplankton cell x is smaller or larger (x∈H1), there exists a unique boundary equilibrium E1 in the area (i), which is globally asymptotically stable (see Theorem 3). With the increase of cell size x∈H2, the positive equilibrium E∗ appears in (ii) and it is globally asymptotically stable, but the E1 becomes unstable (see Theorems 2 and 4). Then, we adopt x=0.05∈H1 and x=0.51∈H2; the corresponding results concerning the global stability of E1=11.76470588,0 and E∗=2.414081849,4.754928396 are shown in (i) and (ii) of Figure 3(b), respectively. It is worth noting that the only difference between (i) and (ii) in Figure 3(b) is that the cell size of phytoplankton is different. Actually, from the Figure 3(a), it is clear that the increase of cell size can lead to the stability changes at the boundary equilibrium E1: stable 0<x<0.07807216867 ⟶ unstable 0.07807216867<x<2.981219088 ⟶ stable x>2.981219088, and the existence of positive equilibrium E∗ changes: nonexistence 0<x<0.07807216867 ⟶ existence 0.07807216867<x<2.981219088 ⟶ nonexistence x>2.981219088, which indicates that the cell size plays an important role in determining the existence and stability of the equilibria in model (1).Figure 3
(a) Bifurcation diagram with respect tox and I, where the white area (i) indicates that there exists the only boundary equilibrium E1 that is globally asymptotically stable, and the pink area (ii) implies that the positive equilibrium E∗ is globally asymptotically stable but the boundary equilibrium E1 is unstable. (b) The global stability of the equilibria in model (1) with different cell size x, in which (i) for the boundary equilibrium E1 with x=0.05 and (ii) for the positive equilibrium E∗ with x=0.51.
(a)
(b)


## 4.1.1. Impact of the Cell Size Effect on the Nutrient Concentration and Phytoplankton Density

Based on the analytical expressions for the positive equilibrium level of nutrient and phytoplankton, one can know that the nutrient concentration and phytoplankton density are closely related to cell sizex. However, let x∈H1; the increase of cell size x does not affect the nutrient concentration, which is clearly demonstrated in Figures 1(a) and 1(c). Furthermore, in this case, Figures 1(d) and 1(f) show that the density of phytoplankton is negative or zero. That is, the phytoplankton in model (1) cannot survive when the phytoplankton cell is smaller or larger (x∈H1), which may be the main reason why the nutrient concentration does not change with the increase of cell size.Figure 1
(a–c) Bifurcation diagram of nutrient in model (1) with respect to x. (d–f) Bifurcation diagram of phytoplankton in model (1) with respect to x. (a) and (d) for x∈0.01,0.07807216867. (b) and (e) for x∈0.07807216867,2.981219088. (c) and (f) for x∈2.981219088,7.
(a)
(b)
(c)
(d)
(e)
(f)Figure 2
(a–c) Bifurcation diagram of nutrient in model (1) with respect to x∈H2. (b–f) Bifurcation diagram of phytoplankton in model (1) with respect to x∈H2. (a) and (d) for x∈0.07807216867,0.1. (b) and (e) for x∈0.1,0.6. (c) and (f) for x∈0.6,2.981219088.
(a)
(b)
(c)
(d)
(e)
(f)However, as the cell size increases to the domain ofH2, the nutrient concentration is a nonmonotonic concave function with respect to cell size x (see Figure 1(b)) and phytoplankton density is a nonmonotonic convex function of cell size x (see Figure 1(e)). Nevertheless, by further analysis, we find that(1)
Ifx∈0.07807216867,0.1, the nutrient concentration is a monotonic decreasing function with respect to cell size x (see Figure 2(a)), but phytoplankton density is a convex-line increasing function of cell size x (see Figure 2(d)), which signifies that phytoplankton can reproduce rapidly and abundantly and absorb a large amount of nutrients during its growth, resulting in a rapid increase in phytoplankton density and a corresponding rapid decrease in nutrient concentration. In fact, from Figure 1(e) or Figures 2(d) and 2(e), we can know that the peak of phytoplankton density appears at x=0.1, and so, x=0.1 is the optimal size for promoting the rapid reproduction of the phytoplankton population.
(2)
Ifx∈0.1,0.6, both the nutrient concentration and phytoplankton density are concave-line decreasing function of cell size x, which are shown in Figures 2(b) and 2(e), respectively. As the density of phytoplankton reaches its peak at x=0.1, when the cell size of phytoplankton exceeds its optimal size, the phytoplankton density decreases with the increase of cell size, but its density is still very high, resulting in the decrease of nutrient concentration, whereas the reduction rate of the nutrient concentration is decreasing compared with that in Figure 2(a).
(3)
Ifx∈0.6,2.981219088, the nutrient concentration is a concave-line increasing function with respect to x (see Figure 2(c)), but phytoplankton density is a concave-line decreasing function of cell size, as is shown in Figure 2(f). Obviously, in this case, with the increase of cell size, the density of phytoplankton decreases dramatically with the rapid increase of nutrient concentration.


## 4.1.2. Impact of Cell Size Effect on the Stability of the Equilibria

The analysis indicates that the plane is divided into two regions (i) and (ii) byI−x parameters, as displayed in Figure 3(a). However, when the phytoplankton cell x is smaller or larger (x∈H1), there exists a unique boundary equilibrium E1 in the area (i), which is globally asymptotically stable (see Theorem 3). With the increase of cell size x∈H2, the positive equilibrium E∗ appears in (ii) and it is globally asymptotically stable, but the E1 becomes unstable (see Theorems 2 and 4). Then, we adopt x=0.05∈H1 and x=0.51∈H2; the corresponding results concerning the global stability of E1=11.76470588,0 and E∗=2.414081849,4.754928396 are shown in (i) and (ii) of Figure 3(b), respectively. It is worth noting that the only difference between (i) and (ii) in Figure 3(b) is that the cell size of phytoplankton is different. Actually, from the Figure 3(a), it is clear that the increase of cell size can lead to the stability changes at the boundary equilibrium E1: stable 0<x<0.07807216867 ⟶ unstable 0.07807216867<x<2.981219088 ⟶ stable x>2.981219088, and the existence of positive equilibrium E∗ changes: nonexistence 0<x<0.07807216867 ⟶ existence 0.07807216867<x<2.981219088 ⟶ nonexistence x>2.981219088, which indicates that the cell size plays an important role in determining the existence and stability of the equilibria in model (1).Figure 3
(a) Bifurcation diagram with respect tox and I, where the white area (i) indicates that there exists the only boundary equilibrium E1 that is globally asymptotically stable, and the pink area (ii) implies that the positive equilibrium E∗ is globally asymptotically stable but the boundary equilibrium E1 is unstable. (b) The global stability of the equilibria in model (1) with different cell size x, in which (i) for the boundary equilibrium E1 with x=0.05 and (ii) for the positive equilibrium E∗ with x=0.51.
(a)
(b)


## 4.2. Impact of Cell Size Factor in the Nutrient-Phytoplankton Dynamics of Model (3)

Now, we turn to study the impact of cell size effect on the stochastic dynamics of model (3) by using the Milstein method mentioned in Higham [73]. In this way, model (3) can be rewritten as the following discretization equation:(88)Nk+1=Nk+I−μxQxNkPka+Nk−mNkΔt−δμxQxNkPka+NkΔtξk+δ22μxQxNkPka+Nkξk2−1Δt,Pk+1=Pk+μxNkPka+Nk−sxPk−dPkΔt+δμxNkPka+NkΔtξk+δ22μxNkPka+Nkξk2−1Δt,where ξkk=1,2,…,n are independent Gaussian random variables N0,1.


### 4.2.1. Impact of Cell Size and Environmental Noise Effects on the Survival of Phytoplankton

First of all, we fixx=0.51 and vary δ=0.1,1.8 to investigate how environmental noise affects the persistence of model (3).When we takeδ=0.1, by simple computations, it is not difficult to get that the condition (i) in Theorem 6 and the condition in Theorem 8 are satisfied. From Theorems 6 and 8, it can be obtained that model (3) is persistent in the mean and exists a unique stationary distribution, as is shown in Figure 4. This result signifies that the small environmental noise disturbance can make model (3) maintain some stability in the random sense and drive the solutions to be perturbed near the positive equilibrium of the deterministic model (1).Figure 4
(a) and (c) are the solutions of nutrient and phytoplankton in the stochastic model (3) and its corresponding deterministic model (1) with initial value N0,P0=8.5,2.4, respectively. (b) and (d) are the probability density function of nutrient and phytoplankton in model (3), respectively.
(a)
(b)
(c)
(d)When we obtain the noise intensityδ=1.8, then the condition (ii) of Theorem 7 is satisfied. Based on Theorem 7, the phytoplankton goes to extinction with probability one, which implies that the larger environmental noise can result in the extinction of phytoplankton. Figure 5(b) clearly confirms this result. Further, in this case, the nutrient concentration is found to increase to N=N1=11.76470588>N∗=2.414081849 (see Figure 5(a)), which agrees well with Remark 2.Figure 5
The solutions of the stochastic model (3) and its corresponding deterministic model (1) with initial value N0,P0=8.5,2.4. (a) For the nutrient. (b) For the phytoplankton.
(a)
(b)Next, we fixδ=0.1 and change x=0.001,0.15,0.45,0.95 to see the effect of cell size on the survival of phytoplankton in model (3). It is clear from Figure 6 that if x=0.15or0.45⊂H2, the phytoplankton is persistent and becomes extinct when x=0.001or0.98⊂H1, which indicates that the smaller or larger cell size can result in the extinction of phytoplankton (see the black and pink curves) and the intermediate cell size is advantage for the persistence of phytoplankton (see the blue and green curves). Moreover, by comparing the green curve with the blue curve, it can be found that the increase of cell size can reduce the random variation of phytoplankton caused by environmental noise disturbance and decrease the density of phytoplankton.Figure 6
The solution of phytoplankton in model (3) with different cell size x, where the noise intensity is fixed to δ=0.1, and the black curve is for x=0.001, the green curve is for x=0.15, the blue curve is for x=0.45, and the pink curve is for x=0.95.


### 4.2.2. Impact of Cell Size Effect on the Existence of Ergodic Stationary Distribution

In order to explore how cell size affects the existence of stationary distribution of phytoplankton and nutrient in model (3), we fix the noise intensity δ=0.1 and vary the cell size(89)x=0.01,0.02,0.03,0.04,0.05,0.06,which satisfy the parameter condition in Theorem 8. Hence, there exists a stationary distribution of model (3), as is shown in Figure 7.Figure 7
The effect of cell size on the existence of stationary distribution for model (3) with different cell size x=0.01,0.02,0.03,0.04,0.05,0.06. (a) and (c) are the probability density function of phytoplankton and nutrient in model (3), respectively. (b) and (d) are the boxplot of the solution for phytoplankton and nutrient in model (3) with the corresponding cell size, respectively.
(a)
(b)
(c)
(d)In Figure7(a), one can find that, with the increase of cell size, the height of probability density function of phytoplankton decreases, but that of the probability density function of nutrition increases (see Figure 7(c)), which suggests that the increase of cell size is not conducive to stabilize the nutrient-phytoplankton dynamics in a stochastic sense.In Figure7(b), it is shown that the maximum value of the boxplot for phytoplankton increases, while the maximum value of the boxplot for nutrient remains unchanged (see Figure 7(d)), as the cell size increases. This implies that the probability density function of phytoplankton is positively skewed (from left to right), but the probability density function of nutrition is not shifted as the cell size increases.


## 4.2.1. Impact of Cell Size and Environmental Noise Effects on the Survival of Phytoplankton

First of all, we fixx=0.51 and vary δ=0.1,1.8 to investigate how environmental noise affects the persistence of model (3).When we takeδ=0.1, by simple computations, it is not difficult to get that the condition (i) in Theorem 6 and the condition in Theorem 8 are satisfied. From Theorems 6 and 8, it can be obtained that model (3) is persistent in the mean and exists a unique stationary distribution, as is shown in Figure 4. This result signifies that the small environmental noise disturbance can make model (3) maintain some stability in the random sense and drive the solutions to be perturbed near the positive equilibrium of the deterministic model (1).Figure 4
(a) and (c) are the solutions of nutrient and phytoplankton in the stochastic model (3) and its corresponding deterministic model (1) with initial value N0,P0=8.5,2.4, respectively. (b) and (d) are the probability density function of nutrient and phytoplankton in model (3), respectively.
(a)
(b)
(c)
(d)When we obtain the noise intensityδ=1.8, then the condition (ii) of Theorem 7 is satisfied. Based on Theorem 7, the phytoplankton goes to extinction with probability one, which implies that the larger environmental noise can result in the extinction of phytoplankton. Figure 5(b) clearly confirms this result. Further, in this case, the nutrient concentration is found to increase to N=N1=11.76470588>N∗=2.414081849 (see Figure 5(a)), which agrees well with Remark 2.Figure 5
The solutions of the stochastic model (3) and its corresponding deterministic model (1) with initial value N0,P0=8.5,2.4. (a) For the nutrient. (b) For the phytoplankton.
(a)
(b)Next, we fixδ=0.1 and change x=0.001,0.15,0.45,0.95 to see the effect of cell size on the survival of phytoplankton in model (3). It is clear from Figure 6 that if x=0.15or0.45⊂H2, the phytoplankton is persistent and becomes extinct when x=0.001or0.98⊂H1, which indicates that the smaller or larger cell size can result in the extinction of phytoplankton (see the black and pink curves) and the intermediate cell size is advantage for the persistence of phytoplankton (see the blue and green curves). Moreover, by comparing the green curve with the blue curve, it can be found that the increase of cell size can reduce the random variation of phytoplankton caused by environmental noise disturbance and decrease the density of phytoplankton.Figure 6
The solution of phytoplankton in model (3) with different cell size x, where the noise intensity is fixed to δ=0.1, and the black curve is for x=0.001, the green curve is for x=0.15, the blue curve is for x=0.45, and the pink curve is for x=0.95.


## 4.2.2. Impact of Cell Size Effect on the Existence of Ergodic Stationary Distribution

In order to explore how cell size affects the existence of stationary distribution of phytoplankton and nutrient in model (3), we fix the noise intensity δ=0.1 and vary the cell size(89)x=0.01,0.02,0.03,0.04,0.05,0.06,which satisfy the parameter condition in Theorem 8. Hence, there exists a stationary distribution of model (3), as is shown in Figure 7.Figure 7
The effect of cell size on the existence of stationary distribution for model (3) with different cell size x=0.01,0.02,0.03,0.04,0.05,0.06. (a) and (c) are the probability density function of phytoplankton and nutrient in model (3), respectively. (b) and (d) are the boxplot of the solution for phytoplankton and nutrient in model (3) with the corresponding cell size, respectively.
(a)
(b)
(c)
(d)In Figure7(a), one can find that, with the increase of cell size, the height of probability density function of phytoplankton decreases, but that of the probability density function of nutrition increases (see Figure 7(c)), which suggests that the increase of cell size is not conducive to stabilize the nutrient-phytoplankton dynamics in a stochastic sense.In Figure7(b), it is shown that the maximum value of the boxplot for phytoplankton increases, while the maximum value of the boxplot for nutrient remains unchanged (see Figure 7(d)), as the cell size increases. This implies that the probability density function of phytoplankton is positively skewed (from left to right), but the probability density function of nutrition is not shifted as the cell size increases.


## 5. Conclusion

In recent years, many field and laboratory evidences indicated that the plankton body size, especially the cell size of phytoplankton, plays a key role in the metabolism, growth, and interaction of phytoplankton [6, 51, 54]. However, little theoretical work has been done in this field. Based on these facts, in this paper, we study analytically and numerically a deterministic nutrient-phytoplankton model taking into account the effects of cell size and its corresponding stochastic version. The consideration of the cell size effect into our proposed model, in this study, makes it different from the traditional nutrient-phytoplankton models [12, 14, 15, 22, 39–45, 61]. The main purpose of this paper is to study the cell size effects on the nutrient-phytoplankton dynamics within the deterministic and stochastic environments.Mathematically, we show that the existence and stability of the equilibria for the deterministic model (1) can be determined by the value of X (i.e., the cell size). In other words, if X<ma/I, there exists the only boundary equilibrium that is locally and globally asymptotically stable (see Theorems 2 and 3); if X<ma/I, the positive equilibrium appears, which is locally and globally asymptotically stable, whereas the boundary equilibrium becomes unstable (see Theorems 2 and 4). For the stochastic model (3), we first prove the existence and uniqueness of the positive solution, and then, the stochastic extinction and persistence in the mean, as well as the existence of ergodic stationary distribution for model (3) are further explored. Ecologically, via numerical simulations, we find that cell size has rich and complex impacts on the nutrient-phytoplankton dynamics in model (1) and model (3) as follow:i If the cell size is smaller or larger (0<x<0.07807216867 or x>2.981219088), the phytoplankton cannot survive (see Figures 1(d) and 1(f)). In these cases, the increase of cell size has no impact on the nutrient concentration (see Figures 1(a) and 1(c)). As the cell size increases to the intermediate size, that is, 0.07807216867<x<2.981219088, the phytoplankton appears and its density as well as nutrient concentration begin to change significantly with cell size.i1 When 0.07807216867<x≤0.1, the increase of cell size can lead to a rapid increase in phytoplankton density and result in a significant reduction in nutrient concentration (see Figures 2(a) and 2(d), respectively). However, the nutrient concentration is still high at this time and the phytoplankton density reached its peak at x=0.1 (see Figures 1(e) and 2(d)), which denotes that the intermediate cell size is the optimum size for the growth of phytoplankton. This result is consistent with the experimental results obtained in [6, 54].i2 When 0.1<x≤0.6, the increase of cell size is capable to decrease the density of phytoplankton and the concentration of nutrient at the same time (see Figures 2(b) and 2(e)). This suggests that this range is most advantageous for reducing the risk of phytoplankton blooms and the possibility of eutrophication of water bodies. Consequently, if the eutrophication of lakes, rivers, and other water bodies, leading to the occurrence of algal blooms, we can invest a number of phytoplankton with cell size x∈0.1,0.6 to possibly control these phenomena.i3 When 0.6<x<2.981219088, the increase of cell size can cause the significant decrease of phytoplankton and the rapid increase of nutrient concentration (see Figures 2(c) and 2(f), respectively). In this case, the density of phytoplankton is very close to zero as the cell size increases. Actually, based on Theorem 2, the phytoplankton will be extinct and the nutrient concentration will eventually increase to a invariant constant as the cell size continuously increases to greater than 2.981219088. Consequently, the increasing cell size of phytoplankton in this region can inhibit the proliferation of phytoplankton (see Figure 2(f)).ii When the cell size is smaller or larger, i.e., 0<x<0.07807216867 or x>2.981219088, there exists the only boundary equilibrium E1 in model (1), which is globally asymptotically stable (see i of Figure 3(b) and Theorem 3), but the positive equilibrium does not exist in these cases (see Theorem 2); when 0.07807216867<x<2.981219088, which lies in the intermediate cell size, the positive equilibrium exists and is globally asymptotically stable, while the boundary becomes unstable (see ii of Figure 3(b) and Theorem 4). Hence, if we can properly control the cell size of phytoplankton to make phytoplankton and nutrient coexist stably, the phytoplankton will not grow rapidly in large quantities; thus, phytoplankton blooms may not occur, which is consistent with the conclusion in i2 and i3.iii With a fixed value of cell size x, the smaller stochastic environmental fluctuations are shown to be unable to affect the persistence of model (3) (see Figure 4 and Theorem 6), but the larger stochastic environmental fluctuations can result in the extinction of phytoplankton with probability one (see Figure 5 and Theorem 7). With a fixed value of noise intensity, the smaller or larger cell size is found to be capable of causing the extinction of phytoplankton, while the intermediate cell size is in favor of the persistence of phytoplankton (see Figure 6). By comparison, the smaller random environmental disturbances and the intermediate cell size have similar positive effects on the persistence of phytoplankton, but the larger cell size, smaller cell size, and larger random environmental disturbances have similar positive effects on the extinction of phytoplankton.iv For the existence of the stationary distribution for model (3), as the cell size increases, the probability density function of phytoplankton has a positive shift and its height is decreasing, but the probability density function of nutrition has hardly moved and its height is increasing (see Figure 7). This indicates that the cell size may be capable to significantly affect the distribution of phytoplankton in water bodies.In view of the above findings, the questions raised in Section1 have been partially well answered. Furthermore, in comparison to the results with respect to cell size obtained in [55–57], the results of this paper are richer and more complex. Though our study cannot directly confirm that cell size is a factor inducing the phytoplankton blooms, we believe that cell size has the potential to influence the formation of this phenomena because it can affect the growth of phytoplankton and the variation of nutrient concentration in the aquatic ecosystems. Hence, the motive of this paper is not to give a feasible theoretical explanation for the size-dependent mechanisms of phytoplankton growth but to provide a new insight into understanding the intrinsic law of the nutrient-phytoplankton dynamics.



---

*Source: 8205696-2019-11-28.xml*