\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}