\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2010(2010), No. 98, pp. 1--15.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2010 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2010/98\hfil Allelopathic phytoplankton model] {Stochastically perturbed allelopathic phytoplankton model} \author[S. Abbas, M. Banerjee\hfil EJDE-2010/98\hfilneg] {Syed Abbas, Malay Banerjee} % in alphabetical order \address{Syed Abbas\newline Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur - 208016 India.\newline Dipartimento di Matematica Universita degli Studi di Bologna Piazza di Porta S. Donato, 5 40126 - Bologna, Italy} \email{sabbas.iitk@gmail.com} \address{Malay Banerjee \newline Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur - 208016 India} \email{malayb@iitk.ac.in} \thanks{Submitted May 18, 2010. Published July 21, 2010.} \subjclass[2000]{93E15, 34C23, 37B25, 34K50} \keywords{Stability; bifurcation; Lyapunov functional; stochastic perturbation} \begin{abstract} In this article we have considered a stochastic delay differential equation model for two species competitive phytoplankton system with allelopathic stimulation. We have extended the deterministic model system to a stochastic delay differential equation model system by incorporating multiplicative white noise terms in the growth equations for both species. We have studied the mean square stability of coexisting state using a suitable Lyapunov functional. Numerical simulation results are provided to validate the analytical findings. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction} The freely floating and weakly swimming organisms within aquatic environment are known as plankton. The plant species of plankton community are known as phytoplankton and they are acting as a basic food source within any aquatic environment. They play crucial role towards the regulation of global carbon cycle and which in turn has a significant impact on the climate. Growth and density of phytoplankton within aquatic environment are controlled by variation of available nutrients, intensity of sunlight, temperature gradient, salinity of water, environmental forcing due to seasonal change and many other factors (for details, see \cite{BEEP}). Apart from these factors, the change in phytoplankton population density of one species has the ability to affect the growth of one or more phytoplankton species by producing allelopathic toxins or stimulators, known as allelo-chemicals. Allelochemicals released by a particular species of phytoplankton have positive as well as negative feedbacks on the growth of other phytoplankton species. Allelopathy is the effect of one phytoplankton species on the growth of one or more phytoplankton species by releasing allelo-chemicals into the habitat are known as `allelopathy' \cite{RICE}. These allelo-chemicals are one of the responsible factors for phytoplankton `bloom'. A dramatic increase in the density of one or more phytoplankton species followed by sudden collapse is known as `{\it phytoplankton bloom}'. In last two decades several harmful phytoplankton blooms are reported in various literatures and this issue have received considerable scientific attention from various researchers \cite{DMA, MBTSRP, mb, HB, GMH, rpdbmb, RP, TRR, TSM}. The first mathematical model for allelopathic interaction between two competing phytoplankton species was introduced by Maynard-Smith \cite{MS} and described by Lotka-Volterra type competition model as follows, \begin{equation}\label{1stmodel} \begin{gathered} \frac{du(t)}{dt}=u(t)(k_1-\alpha_1u(t)-\beta_{12}v(t)-\gamma_1u(t)v(t)), \\ \frac{dv(t)}{dt}=v(t)(k_2-\alpha_2v(t)-\beta_{21}u(t)-\gamma_2u(t)v(t)), \end{gathered} \end{equation} subject to the non-negative initial condition $u(0)\equiv u_0\geq0$ and $v(0)\equiv v_0\geq0$. $u(t)$ and $v(t)$ stand for the densities of two phytoplankton species at any instant of time `$t$'. $k_1$, $k_2$ are the cell proliferation rates measured per unit of time, $\alpha_1$, $\alpha_2$ are the intra-specific competition rates for first and second species respectively, $\beta_{12}$, $\beta_{21}$ stand for interspecific competition rates for limited resources \cite{mb,MS}. $\gamma_1$ and $\gamma_2$ are the toxic inhibition rates for the first and second species respectively. All parameters are positive. Basic dynamical results for the above model system was carried out by Chattopadhyay \cite{JC} and obtained the role of toxic inhibition by either species on the stability behavior of two phytoplankton population. Mukhopadhyay {\it et al.} \cite{mct} have considered similar model for two different cases, namely $\gamma_1$, $\gamma_2$ $>0$ as well as $\gamma_1$, $\gamma_2$ $<0$. The later case corresponds to the positive feedback of allelo-chemicals towards the growth of other phytoplankton species. Both the ordinary differential equation model fail to capture the oscillatory growth of phytoplankton population. The oscillatory growth of either phytoplankton species is resulted from the following delay differential equation (DDE) model of two interacting species \begin{gather*} \frac{du(t)}{dt}=u(t)(k_1-\alpha_1u(t)-\beta_{12}v(t) +\gamma_1u(t)v(t)),\\ \frac{dv(t)}{dt}=v(t)(k_2-\alpha_2v(t)-\beta_{21}u(t) +\gamma_2u(t-\tau)v(t)), %\label{ieq0} \end{gather*} where both the species stimulates the growth of the other. $\tau$ is the discrete time delay taking care of the fact that the phytoplankton cells are not capable of releasing toxic substances immediately after cell proliferation, rather delayed by some time time interval required for maturation of the cell before start producing toxic substances. Solution of this delay differential equation model shows oscillatory dynamics, the stable coexisting state becomes unstable through Hopf-bifurcation and small amplitude periodic solution results in around the interior equilibrium point. Recently we have studied the interaction of two phytoplankton population where one phytoplankton species release allelo-chemicals within the aquatic environment, which acts as a stimulator to the growth of other species \cite{AbbasBanerjeeNorbert}. Our study was based upon the following DDE model \begin{gather*} \frac{du(t)}{dt}= u(t)(k_1-\alpha_1 u(t)-\beta_{12}v(t)) \\ \frac{dv(t)}{dt}= v(t)(k_2-\alpha_2 v(t)-\beta_{21}u(t)+\gamma u(t-\tau)v(t)), \end{gather*} subjected to the biologically feasible initial condition $u(t)=\phi(t)>0$ for $t\in[-\tau,0]$, $v(0)=v_0>0$. In \cite{AbbasBanerjeeNorbert} we have established the existence-uniqueness of solutions, local asymptotic stability of various equilibria, persistence of both species and existence of small amplitude periodic solution through Hopf-bifurcation. The purpose of the present paper is to consider the dynamics of a mathematical model of two competing phytoplankton species where one species releases auxin and this auxin stimulate the growth of second phytoplankton species and second one is not a toxic species. After performing some basic dynamical analysis of the ordinary differential equation model system we extend it to a delay differential equation model to understand the effect of time delay (required for maturation of the species before producing toxic substances) on the dynamics of the model system. Finally we construct a two dimensional stochastic delay differential equation model system to understand the effect of environmental fluctuation on the time evolution of both the species. The more on hereditary systems, interested readers may consult \cite{kn,km,ks}. The organization of this paper is as follows: In section 2 we describe the basic model system for two competitive phytoplankton species and study the related dynamical behavior of both delayed and non delayed model systems; in section 3 we will consider the stability of solution in probability for stochastic delay differential equation model system by using the technique of the paper \cite{sha} by Shaikhet and Kolmanovskii and Shaikhet general method of lyapunov functionals construction \cite{bradul,ks1,ks2,pater,sha1,sha2,sha3}. Finally in section 4 we will discuss the basic outcomes of our analytical findings and their ecological interpretations. \section{Basic Model} The basic model considered in this paper is governed by the following two dimensional nonlinear coupled ordinary differential equation model for two competing phytoplankton species \begin{equation}\label{basiceq} \begin{gathered} \frac{du(t)}{dt}= u(t)(k_1-\alpha_1u(t)-\beta_{12}v(t)) \\ \frac{dv(t)}{dt}= v(t)(k_2-\alpha_2v(t)-\beta_{21}u(t)+\gamma u(t)v(t)), \end{gathered} \end{equation} subjected to the biologically feasible initial condition $u(0)\,\equiv\,u_0\,\geq\,0$ and $v(0)\,\equiv\,v_0\,\geq\,0$. All parameters involved with \eqref{basiceq} are positive and have same interpretation as we have mentioned in introduction, $\gamma$ denotes the rate of allelopathic substance released by the first phytoplankton species. The functions present on the right-hand side of system \eqref{basiceq} are continuous functions of the variables in $\mathbb{R}^2_+=[(u, v) : u, v \ge 0]$. Straightforward computation shows that they are locally Lipschitian on $\mathbb{R}^2_+$ and hence the solutions of \eqref{basiceq} with nonnegative initial conditions exist and are unique. It is also easy to verify that these solutions can be extended for all $t > 0$ and stay nonnegative at all future time. By uniqueness principle the solution starting from $(u_0,0)$ will remain on $u$-axis as $t\rightarrow+\infty$ and analogous result hold for the choice of initial condition on $v$-axis. These results combining with the positivity of solutions starting from an interior point of first quadrant imply $\mathbb{R}^2_+$ is the invariant set for the model system \eqref{basiceq}. Here we state the boundedness of solutions for the model system \eqref{basiceq} with positive initial condition, proof of these results can be found in \cite{AbbasBanerjeeNorbert}. \begin{lemma} \label{lem1} If $\alpha_1\alpha_2-\gamma k_1>0$ then any solution $(u(t),v(t))$ of system \eqref{basiceq} with $(u_0,v_0)\in\operatorname{Int}(\mathbb{R}^2_+)$ satisfies $$ \limsup_{t\rightarrow \infty} u(t)\le \frac{k_1}{\alpha_1}\equiv M_1 , \quad \limsup_{t\rightarrow \infty} v(t)\le \frac{\alpha_1 k_2}{\alpha_1 \alpha_2-\gamma k_1}\equiv M_2. $$ \end{lemma} \begin{lemma} \label{lem2} If parameters involved with the model system satisfy the restrictions $\alpha_1\alpha_2>\gamma k_1$, $k_1>\beta_{12} M_2$ and $k_2\alpha_1>k_1\beta_{21}$ then solutions of the system \eqref{basiceq} with $(u_0,v_0)\in\operatorname{Int}(\mathbb{R}^2_+)$ satisfies $$ \liminf_{t\rightarrow \infty} u(t)\ge\frac{k_1-\beta_{12}M_2}{\alpha_1} \equiv m_1 , \quad \liminf_{t\rightarrow \infty} v(t)\ge \frac{k_2\alpha_1-k_1\beta_{21}}{\alpha_1\alpha_2}\equiv m_2. $$ \end{lemma} The model system \eqref{basiceq} have three equilibrium points on the boundary of first quadrant and they exist with out any restriction on the system parameters. $E_0\equiv(0,0)$ is trivial equilibrium point, and two axial equilibrium points are $\displaystyle E_{10}=(\frac{k_1}{\alpha_1},0)$ and $\displaystyle E_{01}=( 0,\frac{k_2}{\alpha_2})$. Coexisting equilibrium point is the point(s) of intersection of two nullclines $\alpha_1u+\beta_{12}v=k_1$ and $\alpha_2v+\beta_{21}u-\gamma uv=k_2$ within the interior of first quadrant. Here we consider the parametric restriction under which interior equilibrium point is unique. If we denote the unique interior equilibrium point by $E_*\equiv(u_*,v_*)$ then \begin{equation} v_*=\frac{k_1-\alpha_1u_*}{\beta_{12}} \end{equation} and $u_*$ is the positive root of the quadratic equation \begin{equation} \alpha_1\gamma x^2+(\beta_{12}\beta_{21}-k_1\gamma-\alpha_1 \alpha_2)x+k_1\alpha_2-\beta_{12}k_2=0. \end{equation} The parametric restriction $k_1\alpha_2<\beta_{12}k_2$ and $u_* 0. \end{equation} Using the algebraic relations $k_1=\alpha_1u_*+\beta_{12}v_*$ and $k_2+\gamma u_*v_*=\alpha_2v_*+\beta_{21}u_*$ it is easy to see, that the stability conditions can be put into the following form, \begin{equation}\label{LAS1} \gamma< \frac{\alpha_1}{v_*}+\frac{\alpha_2}{u_*}\quad \text{and} \quad \alpha_1 \alpha_2 - \beta_{12} \beta_{21} > ( u_* \alpha_1 -v_* \beta_{12}) \gamma . \end{equation} For the model system under consideration, it is difficult to find out the local stability conditions in terms of parameters explicitly. Now we validate the analytical findings with help of a numerical example. For this purpose we choose the parameter values $k_1=2$, $k_2=1$, $\alpha_1=0.07$, $\alpha_2=0.08$, $\beta_{12}=0.05$, $\beta_{21}=0.015$ and $\gamma=0.003$ \cite{AbbasBanerjeeNorbert}. Fig. 1 shows that two nullclines intersect at the interior equilibrium point $E_*(13.85,20.61)$. The two conditions mentioned in \eqref{LAS1} are satisfied for the chosen parameter values. Fig. 2 depicts that the interior equilibrium point is locally asymptotically stable. \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig1} \end{center} \caption{$u$-nullcline and $v$-nullcline intersect at interior equilibrium point $E_*$.} \label{fig1} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig2} \end{center} \caption{Time evolution of $u(t)$ and $v(t)$ obtained from numerical simulation, local asymptotic stability of $E_*(13.85,20.61)$.} \label{fig2} \end{figure} \subsection{Delayed Model System} Now we extend the ordinary differential equation model system to a delay differential equation model by introducing a discrete time delay parameter. The delay parameter is introduced due to the assumption that every phytoplankton cell of the first species is capable to produce toxic substance after a time period measured from their time of birth \cite{AbbasBanerjeeNorbert, mct}. This assumption leads us to the following system of nonlinear coupled delay differential equations, \begin{equation}\label{delaymodel} \begin{gathered} \frac{du(t)}{dt}= u(t)(k_1-\alpha_1 u(t)-\beta_{12}v(t)) \\ \frac{dv(t)}{dt}= v(t)(k_2-\alpha_2 v(t)-\beta_{21}u(t)+\gamma u(t-\tau)v(t)), \end{gathered} \end{equation} subjected to the biologically reasonable initial conditions: \begin{equation}\label{initialcondition} u(t) = \phi(t)>0 \quad\text{for $t \in [-\tau,0]$, } v(0) = v_0>0. \end{equation} Here, $\tau > 0$ is the discrete time delay parameter, and $\phi \in C^0([-\tau,0];\mathbb{R}^{+})$ is a given function. The delayed system \eqref{delaymodel} has the same equilibrium points as the non-delayed model system \cite{gopal}. Linearizing the model system \eqref{delaymodel} around interior equilibrium point $(u_*,v_*)$ we get the following system of linear equations in terms of perturbation variables $x$ and $y$ are as follows \begin{equation}\label{lineareq} \begin{gathered} \frac{dx(t)}{dt}=A x(t)+By(t),\\ \frac{dy(t)}{dt}=Cx(t)+Dy(t)+Ex(t-\tau), \end{gathered} \end{equation} where \begin{gather*} A=k_1-2\alpha_1u_*-\beta_{12}v_*, \quad B=-u_*\beta_{12}, \quad C=-v_*\beta_{21},\\ D=k_2+2\gamma u_*v_*-2\alpha_2 v_*-\beta_{21}u_*,\quad E=\gamma {v_*}^2. \end{gather*} The characteristic equation associated with the linear system \eqref{lineareq} is \begin{equation} \Delta(\lambda, \tau)=\lambda^2-(A+D)\lambda+(AD-BC)-BE e^{-\lambda \tau}=0. \label{eq06} \end{equation} To study the local asymptotic stability properties of $E_*$, we have to consider the nature of the real parts of the roots of characteristic equation \eqref{eq06}. Following theorem states the conditions required for local asymptotic stability of $E_*$ irrespective to the magnitude of delay parameter `$\tau$', detailed proof of this result can be found in \cite{AbbasBanerjeeNorbert}. \begin{theorem} \label{thm1} Interior equilibrium point $E_*$ remains locally asymptotically stable for all $\tau \ge 0$ whenever following conditions are satisfied: \begin{enumerate} \item The real part of the roots of $\Delta(\lambda,0)=0$ are negative.\label{stab-condition1} \item For all real $\mu$ and all $\tau\ge 0,$ $\Delta(i\mu,\tau)\ne 0$. \label{stab-condition2} \end{enumerate} \end{theorem} The conditions of above theorem are satisfied when coefficients of \eqref{chareq} satisfy following restrictions, \begin{equation} \label{LAScond} A+D<0,\quad AD-BC>BE,\quad|AD-BC|>|BE|. \end{equation} The stability properties of $E_*$ under the parametric restriction $(AD-BC)^2<(BE)^2$ depends upon the magnitude of $\tau$. It is established in \cite{AbbasBanerjeeNorbert}, if $A+D<0,$ then in the parametric region $BE0$. In above model system $\sigma_1, \sigma_2$ are two positive real constants known as intensity of environmental fluctuations and $w_1(t)$, $w_2(t))$ are two independent standard Wiener processes \cite{mao}. Consider the transformation $y_1(t)=u(t)-u^*$ and $y_2(t)=v(t)-v^*,$ this will center the system around the positive equilibrium point. Thus we obtain \begin{equation} \begin{gathered} dy_1(t)= (y_1(t)+u^*)[-\alpha_1y_1(t)-\beta_{12}y_2(t)]dt+\sigma_1y_1(t)dw_1(t) \\ \begin{aligned} dy_2(t)&= (y_2(t)+v^*)[(\gamma u^*-\alpha_2)y_2(t)-\beta_{21}y_1(t)+\gamma v^*y_1(t-\tau) \\ &\quad +\gamma y_2(t)y_1(t-\tau)]dt+\sigma_2y_2(t)dw_2(t). \end{aligned}\label{eq2} \end{gathered} \end{equation} To study the local stability of the system \eqref{eq1}, we only need to study the stability of zero solution of the system \eqref{eq2}. Consider the linear part of the system \eqref{eq2}, \begin{equation} \label{eq3} % \label{linearSDDE} \begin{gathered} dy_1(t)= u^*[-\alpha_1y_1(t)-\beta_{12}y_2(t)]dt+\sigma_1y_1(t)dw_1(t) \\ dy_2(t)= v^*[(\gamma u^*-\alpha_2)y_2(t)-\beta_{21}y_1(t)+\gamma v^*y_1(t-\tau)]dt+\sigma_2y_2(t)dw_2(t). \end{gathered} \end{equation} Before proceeding further we establish the existence and uniqueness of solutions for the stochastic model system \eqref{eq1} using \cite[Th. 2.2]{mao}. Two functions $f\equiv (f_1,f_2)^T$ and $g\equiv(g_1,g_2)^T$ are defined as follows \begin{gather*} f_1(t,U) =u(t)(k_1-\alpha_1 u(t)-\beta_{12}v(t)),\\ f_2(t,U) =v(t)(k_2-\alpha_2 v(t)-\beta_{21}u(t)+\gamma u(t-\tau)v(t)), \\ g_1(t,U)=\sigma_1(u(t)-u^*),\quad g_2(t,U)=\sigma_2(v(t)-v^*), \end{gather*} where $U=(u(t),v(t))$. One can easily verify that $f(t,U)$ and $g(t,U)$ are Lipschitz continuous. Using the positivity of $u(t)$, $v(t)$ and as all parameters involved with $f(t,U)$ are positive, we can write $$ f_1(t,U) \le k_1u(t), $$ and \begin{align*} f_2(t,U) &\le v(t)(k_2-\alpha_2v(t)-\beta_{21}u(t)+\gamma M_1 v(t)) \\ &\le v(t)(k_2-(\alpha_2-\gamma M_1) v(t)-\beta_{21}u(t)) \\ &\le k_2v(t), \end{align*} whenever $\alpha_2-\gamma M_1 \ge 0$, (see Lemma 1). Calculating the norm of $f,$ we get \[ |f(t,U)|^2=|f_1(t,U)|^2+|f_2(t,U)|^2 \le k_1^2u^2(t)+k_2^2v^2(t) \le K^2(1+|U|^2), \] where $K=\max\{k_1,k_2\}$. Calculating the norm for $g$, we get \begin{align*} |g(t,U)|^2&= |g_1(t,U)|^2+|g_2(t,U)|^2 \\ &\le \sigma_1^2(u(t)+u^*)^2+\sigma_2^2(v(t)+v^*)^2 \\ &\le 2\sigma_1^2(u^2(t)+{u^*}^2)+2\sigma_2^2(v^2(t)+{v^*}^2)\\ &\le 2\max\{\sigma_1^2,\sigma_2^2\}({v^*}^2+{u^*}^2) \Big(1+\frac{u^2(t)+v^2(t)}{{u^*}^2+{v^*}^2}\Big)\\ &\le 2\max\{\sigma_1^2,\sigma_2^2\}\max\Big\{1,\frac{1}{{u^*}^2 +{v^*}^2}\Big\}({u^*}^2+{v^*}^2)(1+u^2(t)+v^2(t)) \\ &\le M(1+|U|^2). \end{align*} Hence both the conditions of \cite[Th. 2.2]{mao} are satisfied. Hence the solutions of the stochastic delay differential equation \eqref{eq1} exist uniquely. The stability analysis carried our here is based on the paper of Shaikhet \cite{sha}. The interested reader may see the paper of Shaikhet \cite{sha} for detailed analysis. To show the asymptotic mean square stability of the zero solution of system \eqref{eq2}, we need to construct the Lyapunov functional for the system. Consider the neutral form of the system \eqref{eq3}: \begin{equation} \begin{gathered} \frac{dy_1(t)}{dt} = -u^*\alpha_1y_1(t)- u^*\beta_{12}y_2(t)+ \sigma_1y_1(t)\dot{w}_1(t),\\ \begin{aligned} & \frac{d}{dt}\Big[y_2(t)+\gamma {v^*}^2 \int_{t-\tau}^ty_1(s)ds\Big]\\ &= v^*(\gamma u^*-\alpha_2 )y_2(t)+v^*(\gamma v^*-\beta_{21})y_1(t) +\sigma_2y_2(t)\dot{w}_2(t). \end{aligned} \end{gathered} \label{eq4} \end{equation} The auxiliary system corresponding to system \eqref{eq4} is \begin{equation} \begin{gathered} \frac{dx_1(t)}{dt}= -u^*\alpha_1x_1(t) +\sigma_1x_1(t)\dot{w}_1(t), \\ \frac{dx_2(t)}{dt}= v^*(\gamma u^*-\alpha_2 )x_2(t) +\sigma_2x_2(t)\dot{w}_2(t). \end{gathered} \label{eq5} \end{equation} The generating operator $L$ for the system \eqref{eq5} acting on a suitably defined functional $V$ is given by \begin{align*} LV(t,x)&=\frac{\partial V(t,x(t))}{\partial t}+A^{T}(x(t))\frac{\partial V(t,x(t))}{\partial x}\\ &\quad +\frac{1}{2}\operatorname{Trace} \Big(B^{T}(x(t)) \frac{\partial^2 V(t,x(t))}{\partial^2 x}B(x(t))\Big), \end{align*} where $$ A= \begin{pmatrix} -\alpha_1 u^* & 0 \\ 0 & (\gamma u^*-\alpha_2) v^* \end{pmatrix}, \quad B= \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix}, $$ and $$ \frac{\partial V}{\partial x} =\operatorname{col} \Big(\frac{\partial V}{\partial x_1}, \frac{\partial V}{\partial x_2} \Big), \quad \frac{\partial^2 V}{\partial^2 x} =\Big(\frac{\partial^2 V}{\partial x_i \partial x_j} \Big)_{i,j=1,2}. $$ For more details on stochastic calculus, one may consult \cite{fima}. \begin{lemma} \label{lemma1} The zero solution of the system \eqref{eq5} is asymptotically mean square stable if $$ \sigma_1^2 < 2u^* \alpha_1 , \quad \sigma_2^2 < 2v^*(\alpha_2 -\gamma u^*). $$ \end{lemma} \begin{proof} Consider the function $X=x_1^2+x_2^2$. We show that $X$ works as a Lyapunov function for system \eqref{eq5}. Let $L_0$ be the generating operator of the system \eqref{eq5} corresponding to Lyapunov functional $X$. Then we have \begin{align*} L_0X&=2x_1(-u^*\alpha_1x_1(t))+2x_2(-v^*(\alpha_2 -\gamma u^*)x_2(t))+\sigma_1^2x_1^2+\sigma_2^2x_2^2 \\ &= -(2u^*\alpha_1-\sigma_1^2)x_1^2(t)-(2v^*(\alpha_2 -\gamma u^*)-\sigma_2^2)x_2^2(t). \end{align*} Let $$ M=\min\{2u^*\alpha_1-\sigma_1^2, 2v^*(\alpha_2 -\gamma u^*)-\sigma_2^2\}, $$ then we have $$ L_0X \le -M (x_1^2(t)+x_2^2(t))=-M|x(t)|^2, $$ where $x(t)=(x_1(t),x_2(t))$. Thus the zero solution of system \eqref{eq5} is asymptotically mean square stable. \end{proof} To analyze the mean square stability of $E_*$ we consider the functional $V_1$ in order to construct Lyapunov functional as follows, $$ V_1=y_1^2(t)+ \Big(y_2(t)+\gamma {v^*}^2 \int_{t-\tau}^ty_1(s)ds\Big)^2. $$ The Lyapuniv functional for the system \eqref{eq3} is given by $V_1+V_2$. Calculating $LV_1$ and after simplification we get \begin{align*} LV_1 &\le -[2(\alpha_1)u^*-\sigma_1^2-(\beta_{12})u^* -(\beta_{21}-\gamma v^*)v^* +\gamma {v^*}^3(\beta_{21} -\gamma v^*)\tau]y_1^2(t)\\ &\quad - [2(\alpha_2 -\gamma u^*)v^*-\sigma_2^2-(\beta_{12})u^*-(\beta_{21}-\gamma v^*)v^* +\gamma {v^*}^3(\alpha_2 -\gamma u^*)\tau]y_2^2(t)\\ &\quad -\gamma {v^*}^3[(\alpha_2 -\gamma u^*)+(\beta_{21}-\gamma v^*)]\int_{t-\tau}^ty_1^2(s)ds. \end{align*} % 3.7 Next we choose $V_2$ as follows \[ V_2=\gamma {v^*}^3[(\gamma u^*-\alpha_2 )+(\gamma v^*-\beta_{21})]\int_{t-\tau}^t(\theta-t+\tau)y_1^2(\theta)d\theta. \] %3.8 Then calculating $LV$ for $V=V_1+V_2$ we have \begin{equation} \begin{aligned} LV &\le -[2(\alpha_1)u^*-\sigma_1^2-(\beta_{12})u^*-(\beta_{21}-\gamma v^*)v^* +\gamma {v^*}^3(\beta_{21}-\gamma v^*)\tau \\ &\quad +\gamma {v^*}^3((\alpha_2 -\gamma u^*)+(\beta_{21}-\gamma v^*))\tau]y_1^2(t) \\ &\quad -[2(\alpha_2 -\gamma u^*)v^*-\sigma_2^2-(\beta_{12})u^* -(\beta_{21}-\gamma v^*)v^* +\gamma {v^*}^3(\alpha_2 -\gamma u^*)\tau]y_2^2(t). \end{aligned} \label{eq6} \end{equation} Thus we have the functional $V=V_1+V_2$. So we conclude that if both quantities inside the square bracket of \eqref{eq6} are positive, then the zero solution of system \eqref{eq3} is asymptotically mean square stable according to \cite[Th. 1]{sha}. Expressions under square bracket in \eqref{eq6} can be put into the following form: \begin{equation} \begin{gathered} \sigma_1^2 < 2\alpha_1u^*-(\beta_{12})u^*-(\beta_{21}-\gamma v^*)v^*+\gamma {v^*}^3((\alpha_2 -\gamma u^*)+2(\beta_{21}-\gamma v^*))\tau \\ \sigma_2^2 < 2(\alpha_2 -\gamma u^*)v^*-(\beta_{12})u^*+\gamma {v^*}^3(\alpha_2 -\gamma u^*)\tau. \end{gathered}\label{cond1} \end{equation} The conditions \eqref{cond1} are well defined if $$ 2\alpha_1u^*-(\beta_{12})u^*-(\beta_{21}-\gamma v^*)v^* +\gamma {v^*}^3((\alpha_2 -\gamma u^*)+2(\beta_{21}-\gamma v^*))\tau>0 $$ and $$ 2(\alpha_2 -\gamma u^*)v^*-(\beta_{12})u^*+\gamma {v^*}^3(\alpha_2 -\gamma u^*)\tau>0. $$ The function $V=V_1+V_2$ satisfies all the conditions of \cite[Th. 1]{sha}, implying stochastic mean square stability of the model system \eqref{eq4} under the restriction \eqref{cond1}. Now we are in a position to discuss the stability of solution in probability for the model system \eqref{eq2}. \begin{theorem} \label{thm3} If the condition \eqref{cond1} holds, then the zero solution of system \eqref{eq2} is stable in probability. In other words the coexisting equilibrium point $E_4$ is stable in probability for the model system \eqref{eq1}. \end{theorem} \begin{proof} Consider the system \eqref{eq2} in the neutral form \begin{equation} \begin{gathered} \frac{dy_1(t)}{dt}= -(y_1(t)+u^*)[\alpha_1y_1(t)+\beta_{12}y_2(t) ]+\sigma_1y_1(t)\dot{w}_1(t), \\ \begin{aligned} &\frac{d}{dt}\Big[y_2(t)+\gamma {v^*}^2 \int_{t-\tau}^ty_1(s)ds\Big]\\ &= -v^*\Big((\alpha_2-\gamma u^*)y_2(t)+\beta_{21}y_1(t) -\gamma v^*y_1(t) -\gamma y_2(t)y_1(t-\tau)\Big) \\ &\quad -y_2(t)\Big((\alpha_2-\gamma u^*)y_2(t)+\beta_{21}y_1(t)-\gamma v^*y_1(t-\tau) \\ & \quad -\gamma y_2(t)y_1(t-\tau)\Big) +\sigma_2y_2(t)\dot{w}_2(t). \end{aligned} \label{eq10} \end{gathered} \end{equation} For the functional $V_1$ defined in the previous lemma, calculating and simplifying $LV_1$ we obtain \begin{equation} \begin{aligned} LV_1 &= -2(\alpha_1)y_1^3(t)-2(\alpha_2-\gamma u^*)y_2^3(t)-2u^*(\alpha_1)y_1^2(t) - 2v^*(\alpha_2-\gamma u^*)y_2^2(t)\\ &\quad -2(\beta_{12})y_1^2(t)y_2(t) -2[u^*(\beta_{12})+v^*(\beta_{21}-\gamma v^*)+\beta_{21}]y_1(t)y_2(t)\\ &\quad +2v^*\gamma y_2^2(t)y_1(t-\tau) - 2\gamma v^*y_2^2(t)y_1(t-\tau)+2\gamma y_2^2(t)y_1(t-\tau)\\ &\quad - 2\gamma {v^*}^2\Big[v^*(\alpha_2-\gamma u^*)\int_{t-\tau}^ty_2(t)y_1(s)ds +v^*(\beta_{21}-\gamma v^*) \int_{t-\tau}^ty_1(t)y_1(s)ds\\ &\quad -\gamma \int_{t-\tau}^ty_2(t)y_1(t-\tau)y_1(s)ds + (\alpha_2-\gamma u^*)\int_{t-\tau}^ty_2^2(t)y_1(s)ds\\ &\quad +\beta_{21}\int_{t-\tau}^ty_1(t)y_2(t)y_1(s)ds -v^*\gamma \int_{t-\tau}^ty_2(t)y_1(t-\tau)y_1(s)ds\\ &\quad -\gamma \int_{t-\tau}^ty_1(t-\tau)y_2^2(t)y_1(s)ds\Big] +\sigma_1^2 y_1^2(t)+\sigma_2^2 y_2^2(t). \end{aligned} \label{eq11} \end{equation} Assume that there exits a positive quantity $\delta$ such that $\sup_{t\ge \tau}|y_i(s)|<\delta$ for $i=1,2$. Using this estimate we can write \begin{align*} LV_1 & \le \Big[-2u^*(\alpha_1)+2\delta(\alpha_1)+\sigma_1^2+ 2\delta(\beta_{12}) \\ &\quad+ [u^*(\beta_{12})+v^*(\beta_{21}-\gamma v^*)+\beta_{21}]- \gamma {v^*}^2[v^*(\beta_{21}-\gamma v^*)+\beta_{21}\delta]\tau \Big]y_1^2(t) \\ &\quad+ \Big[-2v^*(\alpha_2-\gamma u^*)+2\delta(\alpha_2-\gamma u^*)+\sigma_2^2-4 \delta\gamma v^*-2\gamma \delta \\ &\quad -[u^*(\beta_{12})+v^*(\beta_{21}-\gamma v^*)+\beta_{21}]\\ &\quad -\gamma {v^*}^2[(v^*+\delta)(\alpha_2-\gamma u^*)-2\gamma v^* \delta-\gamma \delta^2]\tau\Big]y_2^2(t) \\ &\quad - \gamma {v^*}^2 \Big[(v^*+\delta)[(\alpha_2-\gamma u^*)+(\beta_{21}-\gamma v^*)-\gamma \delta]\Big] \int_{t-\tau}^ty_1^2(s)ds. \end{align*} For the present problem we choose $V_2$ as \[ V_2=\gamma {v^*}^2 \Big[(v^*+\delta)[(\gamma u^*-\alpha_2)+(\gamma v^*-\beta_{21})+\gamma \delta]\Big]\int_{t-\tau}^t(\theta-t+\tau)y_1^2(\theta)d\theta. \] Then for $V_1+V_2$ we can write, \begin{equation} \begin{aligned} &L(V_1+V_2)\\ &\le \Big[-2u^*(\alpha_1)+2\delta(\alpha_1)+\sigma_1^2 + 2\delta(\beta_{12})+ [u^*(\beta_{12}) +v^*(\beta_{21}-\gamma v^*)+\beta_{21}]\\ &\quad -\gamma {v^*}^2\Big(v^*(\beta_{21}-\gamma v^*)+\beta_{21}\delta +(v^*+\delta)[(\alpha_2-\gamma u^*) +(\beta_{21}-\gamma v^*)\\ &\quad -\gamma \delta]\Big)\tau \Big]y_1^2(t) +\Big[-2v^*(\alpha_2-\gamma u^*)+2\delta(\alpha_2-\gamma u^*)\\ &\quad +\sigma_2^2-4 \delta\gamma v^*-2\gamma \delta - [u^*(\beta_{12})+v^*(\beta_{21}-\gamma v^*)+\beta_{21}]\\ &\quad -\gamma {v^*}^2[(v^*+\delta)(\alpha_2-\gamma u^*)-2\gamma v^* \delta-\gamma \delta^2]\tau\Big]y_2^2(t). \end{aligned} \label{eq14} \end{equation} Therefore, under the condition \eqref{cond1} we have $L(V_1+V_2) \le 0$ for sufficiently small $\delta$. In order to use \cite[Th. 2]{sha} we need the functional to be positive definite. For this purpose we consider $$ W_1=y_1^2(t)+y_2^2(t). $$ Using equation \eqref{eq2} and assuming that $\sup_{t\ge \tau}|y_i(s)| < \delta$, we have \begin{equation} \begin{aligned} LW_1 &\le \Big[-2u^*\alpha_1+2\delta\alpha_1+\sigma_1^2+ 2\delta \beta_{12}+ u^*\beta_{12}+v^*\beta_{21}\Big]y_1^2(t) \\ &\quad+ \Big[-2v^*(\alpha_2-\gamma u^*)+2\delta(\alpha_2-\gamma u^*)+\sigma_2^2- 4\delta\gamma v^* - 2\gamma \delta^2\\ &\quad + u^*(\beta_{12}-\gamma u^*)+v^*\beta_{21}+2\beta_{21}\delta+2{v^*}^2\gamma \Big]y_2^2(t) -2{v^*}^2 \gamma y_1^2(t-\tau). \end{aligned} \label{eq15} \end{equation} Taking the functional $W_2$ of the form $$ W_2=-2{v^*}^2 \gamma \int_{t-\tau}^ty_1^2(s)ds, $$ we obtain \begin{equation} \begin{aligned} &L(W_1+W_2)\\ &\le \Big[-2u^*\alpha_1+2\delta\alpha_1+\sigma_1^2+ 2\delta \beta_{12}+ u^*\beta_{12}+v^*\beta_{21}\Big]y_1^2(t) \\ &\quad + \Big[-2v^*(\alpha_2-\gamma u^*)+2\delta(\alpha_2-\gamma u^*)+\sigma_2^2- 4\delta\gamma v^* \\ &\quad - 2\gamma \delta^2+ u^*(\beta_{12}-\gamma u^*)+v^*\beta_{21}+2\beta_{21}\delta+2{v^*}^2\gamma-2{v^*}^2\gamma \Big]y_2^2(t). \end{aligned} \label{eq16} \end{equation} Finally consider the functional $V=V_1+V_2+\lambda(W_1+W_2)$. It is easy to verify that the functional $V$ satisfies all the conditions of \cite[Th. 2]{sha} for sufficiently small $\lambda$. Hence zero solution of the system \eqref{eq2} is stable in probability under the satisfaction of the condition \eqref{cond1}. \end{proof} Now we present some numerical simulation results of the following stochastic delay differential equation model system for different values of the delay parameter $\tau$ and forcing intensities $\sigma_1$ and $\sigma_2$, \begin{equation} \begin{gathered} du(t)= u(t)[2-.07 u(t)-.05v(t)]dt+\sigma_1(u(t)-u^*)dw_1(t), \\ \begin{aligned} dv(t)&= v(t)[1-.08 v(t)-.015 u(t)+.003 u(t-\tau)v(t)]dt\\ &\quad + \sigma_2(v(t)-v^*)dw_2(t). \end{aligned}\label{example} \end{gathered} \end{equation} We have used the `Euler-Maruyama' (EM) scheme \cite{baker} for numerical simulations. The critical magnitudes of environmental forcing to maintain the stability of co-existing equilibrium point for the above model system is given by (see \eqref{cond1}) \begin{equation} \sigma_1^2<2.21-1.25 \tau, \quad \sigma_2^2<.89-1.01 \tau. \label{40} \end{equation} For feasible values of $\sigma_1$ and $\sigma_2$ we have to take $\tau<0.88$. Relation \eqref{40} defines the bounds for $\sigma_1$ and $\sigma_2$ depending upon the magnitudes of $\tau$. For $\tau=0.2$ we have to take $\sigma_1<1.4$, $\sigma_2<0.83$ and for $\tau=0.5$ the stability in probability for interior equilibrium point $E_*(13.85,20.61)$ demands the restrictions $\sigma_1<1$ and $\sigma_2<0.6$. We can not perform numerical simulation beyond $\tau=0.88$ as in that case the inequalities defined in \eqref{40} are meaningless and the analysis carried out here cannot ensure the stability in probability for the coexistence state. Numerical simulations for the model system \eqref{example} are carried out for two sets of values for $\tau$, $\sigma_1$ and $\sigma_2$ keeping in mind the restriction \eqref{40}. Five sample trajectories are plotted in Fig. 3 and Fig. 4 against time showing stability in probability in both the cases. The parameter values are mentioned in the caption of respective figures. The magnitudes for $\sigma_1$ and $\sigma_2$ with $\tau=0.5$ are chosen greater than their values for $\tau=0.2$. Simulation results depict the fact that the stability in probability is not hampered or altered with the magnitude of forcing intensities once they satisfy the restriction \eqref{40}. \begin{figure}[ht] \begin{center} \includegraphics[width=0.9\textwidth]{fig3} % abbas_sto_1.eps \end{center} %\label{figsto1} \caption{Five simulation results for the model system \eqref{example} for the parametric values $\tau=0.2$, $\sigma_1=0.8$ and $\sigma_2=0.5$.} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.9\textwidth]{fig4} % abbas_sto_2.eps \end{center} \label{figsto1} \caption{Numerical simulation results for the model system \eqref{example} with parametric values $\tau=0.5$, $\sigma_1=1.0$ and $\sigma_2=0.6$.} \end{figure} \section{Conclusion} In this article we have considered the stochastic delay differential equation model for two interacting phytoplankton population where one population produce auxin which enhance the growth of other population. In a recent study \cite{AbbasBanerjeeNorbert} we have studied the dynamics of the same model system with deterministic setup and obtained the parametric restrictions required for the persistence and calculated the threshold magnitude for delay parameter beyond which both the population show significant oscillation about the coexistence steady-state. Here we have extended the model system by introducing multiplicative noise terms to the growth equations of both population where strength of the noise is proportional to distance of $u(t)$ and $v(t)$ from their equilibrium levels. This type of stochastic formulation was firstly proposed in \cite{BerettaKolmanowskiiShaikhet,sha} and later used by several other authors \cite{BandyopadhyayChattopadhyay,Carletti02,Carletti06,Carletti07,pater1} to study the effect of environmental fluctuation on the dynamics of interaction populations. We have obtained the conditions for stability in probability around the coexisting equilibrium point of our stochastic delay differential equation model. The main outcomes of this paper is summarized in Lemma \ref{lemma1} and Th. $3$. It is clear that stability of SDDE system demand some additional restrictions should be satisfied by the magnitude of delay and intensity of environmental driving force. Analytical findings ensure that the magnitude of discrete time delay $\tau$ plays a crucial role to determine the stability in probability of coexisting equilibrium point as well as critical magnitude of environmental forcing intensities. Numerical simulations are carried out to validate the analytical findings and we have observed that all trajectories approach the equilibrium level in the sense of probability. It is interesting to note that the relation \eqref{40} defines the bounds for $\sigma_1$ and $\sigma_2$ in terms of $\tau$ for specific choice of other parameters. This also explains the fact that the stability in probability around coexisting equilibrium point for the model system \eqref{eq1} is possible when the intensities of fluctuation is not arbitrarily large. Finally we would like to remark that the stability of stochastic delay differential equation model constructed by perturbing one or more system parameters and then study of stability behavior remains an open problem for the model considered in this paper and this will be our next project in the near future. \subsection*{Acknowledgments} The authors would like to thanks the anonymous referees for their constructive comments and suggestions which helped us to improve the original manuscript considerably. The first author was supported by grant 2/40(31)/2009-R\&D-II/367 from NBHM. \begin{thebibliography}{99} \bibitem{AbbasBanerjeeNorbert} Abbas, S., Banerjee, M. and Hungerbuhler, N.; Existence, uniqueness and stability analysis of allelopathic stimulatory phytoplankton model, {\it J. Math. Anal. Appl.}, 367, (2010), 249 - 259. \bibitem{DMA} Anderson, D. M.; Toxic algae blooms and red tides: a global perspective, in: T. Okaichi, D.M. Anderson, T. Nemoto (Eds.), Red Tides: Biology, Environmental Science and Toxicology, Elsevier, New York, 1989, pp. 11 - 21. \bibitem{baker} Baker C. T. H. and Buckwar, E.; Numerical analysis of explicit one-step methods for stochastic delay differential equations, {LMS J. Comput. Math.}, 3 (2000), 315 - 335. \bibitem{BandyopadhyayChattopadhyay} Bandyopadhyay, M., Chattopadhyay, J.; Ratio-dependent predatorprey model: effect of environmental fluctuation and stability. Nonlinearity 18 (2005), 913 - 936. \bibitem{MBTSRP} Bandyopadhyay, M., Saha, T., Pal, R.; Deterministic and stochastic analysis of a delayed allelopathic phytoplankton model within fluctuating environment, Nonlinear Analysis : Hybrid systems, 2 (2008), 958 - 970. \bibitem{mb} Bandyopadhyay, M.; Dynamical analysis of a allelopathic phytoplankton model, {\it J. Biol. Sys.}, 14, (2006), 205 - 217. \bibitem{BerettaKolmanowskiiShaikhet} Beretta, E., Kolmanovskii, V.B., Shaikhet, L.; Stability of epidemic model with time delays influenced by stochastic perturbations, Math. Comput. Simul. 45 (1998) 269 - 277. \bibitem{HB} Berglund, H.; Stimulation of growth of two marine algae by organic substances exereted by oxic substances on a two species competitive Enteromorpha linza in unialgal and axenic cultures, { Physiol. Plant}, 22 (1969), 10 - 69. \bibitem{bradul} Bradul, N., Shaikhet, L.; Stability of the positive point of equilibrium of Nicholson's blowflies equation with stochastic perturbations: numerical analysis. Discrete Dynamics in Nature and Society, vol. 2007, Article ID 92959, 25 pages, 2007. \bibitem{Carletti02} Carletti, M.; On the stability properties of a stochastic model for phage-bacteria interaction in open marine environment, Math. Biosci. 175 (2002) 117 - 131. \bibitem{Carletti06} Carletti, M.; Numerical simulation of a Campbell-like stochastic delay model for bacteriophage infection, Math. Med. Biol. 23 (2006) 297 - 310. \bibitem{Carletti07} Carletti, M.; Mean-square stability of a stochastic model for bacteriophage infection with time delays, Math. Biosci. {\bf 210}, (2007), 395 - 414. \bibitem{JC} Chattophadyay, J.; Effect of toxic substances on a two species competitive system, { Ecol. Model}, 84 (1996), 287 - 289. \bibitem{BEEP} Edvarsen, B., Paasche, E.; Bloom dynamics and physiology of Primnseium and Chrysochromulina, in Physiological Ecology of Harmful Algal Bloom, Springer, Berlin, 1998. \bibitem{gopal} Gopalsamy, K.; Stability and oscillation in delay differential equations of population dynamics, {Kulwer Academic, Dordrecht}, 1992. \bibitem{GMH} Hallegraeff, G. M., A review of harmful algae blooms and the apparent global increase, Phycologia 32 (1993) 79 - 99. \bibitem{kn} Kolmanovskii, V. B., Nosov, V. R.; Stability of functional differential equations, { N.-Y., Academic press}, 1986. \bibitem{km} Kolmanovskii, V. B., Myshkis, A. D.; Applied theory of functional differential equations, {Boston, Kluwer Academic Publishers}, 1992. \bibitem{ks1} Kolmanovskii, V., Shaikhet, L.; Some peculiarities of the general method of Lyapunov functionals construction. Applied Mathematics Letters, vol. 15, no. 3, (2002), pp. 355 - 360. \bibitem{ks2} Kolmanovskii, V., Shaikhet, L.; Construction of Lyapunov functionals for stochastic hereditary systems: a survey of some recent results. {\em Mathematical and Computer Modelling,} vol. 36, no. 6, (2002) pp. 691 - 716. \bibitem{ks} Kolmanovskii, V. B., Shaikhet, L. E.; Control of systems with aftereffect, Translations of mathematical monographs, Vol: 157. {American Mathematical Society, Providence, RI}, 1996. \bibitem{fima} Klebaner, Von Fima, C.; Introduction to stochastic calculus with applications, {Imperial College Press}, 2005. \bibitem{MS} Maynard Smith, J.; Models in ecology, {Cambridge University Press, Cambridge}, 1974. \bibitem{mct} Mukhopadhyay, A., Chattophadyay, J., Tapasawi, P. K.; A delay differential equation model of plankton allelopathy, {\it Math. Biosci.}, 149, (1998), 167 - 189. \bibitem{mao} Mao, X.; Exponential stability of stochastic differential equation, {Marcel Dekker}, 1994. \bibitem{pater} Paternoster, B., Shaikhet, L.; About stability of nonlinear stochastic difference equations. {\em Applied Mathematics Letters,} vol. 13, no. 5, (2000), pp. 27 - 32. \bibitem{pater1} Paternoster, B., Shaikhet, L.; Stability of equilibrium points of fractional difference equations with stochastic perturbations. {\em Advances in Difference Equations}, vol. 2008, Article ID 718408, (2008), 21 pages. \bibitem{rpdbmb} Pal, R., Basu D., Banerjee, M.; Modelling of phytoplankton allelopathy with Monod–Haldane-type functional response - A mathematical study, {\it Biosystems}, 95, (2009), 243 - 253. \bibitem{RP} Pratt, R.; Influence of the size of the inoculum on the growth of Chlorella vulgaris in freshly prepared culture medium, {Am. J. Bot.}, 27 (1940), 52 - 67. \bibitem{RICE} Rice, E., Allelopathy; Academic Press, New York, 1984. \bibitem{TRR} Rice, T. R.; Biotic influences affecting population growth of planktonic algae, US Fish Wild Serv. {Fish Bull.}, 54 (1954), 227 - 242. \bibitem{sha} Shaikhet, L.; Stabiliy of predator-prey model with aftereffect by stochastic perturbation, { SACTA}, Vol. 1, No.1 (1998), 3 - 13. \bibitem{sha1} Shaikhet, L.; Stability in probability of nonlinear stochastic hereditary systems, Dynamic Systems and Applications, vol. 4, no. 2, (1995), 199 - 204. \bibitem{sha2} Shaikhet, L.; Modern state and development perspectives of Lyapunov functionals method in the stability theory of stochastic hereditary systems. Theory of Stochastic Processes, vol. 2(18), no. 1-2, (1996), 248 - 259. \bibitem{sha3} Shaikhet, L.; Stability of a positive point of equilibrium of one nonlinear system with aftereffect and stochastic perturbations. Dynamic Systems and Applications, vol.17, (2008), 235 - 253. \bibitem{TSM} Smayda, T.; Novel and nuisance phytoplankton blooms in the sea: Evidance for a global epidemic, in: E. Graneli, B. Sundstrom, L. Edler, D.M. Anderson (Eds.), Toxic Marine Phytoplankton, Elsevier, New York, 1990, 29 - 40. \end{thebibliography} \end{document}