\documentclass[reqno]{amsart}
\usepackage{graphicx} 
\usepackage{hyperref} 

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 21, pp. 1--12.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/21\hfil Bifurcation and spatial pattern formation]
{Bifurcation and spatial pattern formation in spreading of disease with
incubation period in a phytoplankton dynamics}

\author[R. S. Baghel, J. Dhar,  R. Jain \hfil EJDE-2012/21\hfilneg]
{Randhir Singh Baghel, Joydip Dhar,  Renu Jain} 

\address{Randhir Singh Baghel \newline
 School of Mathematics and Allied Science, Jiwaji University,
 Gwalior (M.P.)-474011, India}
\email{randhirsng@gmail.com}

\address{Joydip Dhar \newline
Department of Applied Sciences, ABV-Indian Institute of Information
Technology and Management, Gwalior-474010,India}
\email{jdhar@iiitm.ac.in}

\address{Renu Jain \newline
School of Mathematics and Allied Science, Jiwaji University,
Gwalior (M.P.)-474011, India}
\email{renujain3@rediffmail.com}

\thanks{Submitted July 25, 2011. Published February 2, 2012.}
\subjclass[2000]{34C11, 34C23, 34D08, 34D20, 35Q92, 92B05, 92D40}
\keywords{ Phytoplankton dynamics;  reaction-diffusion equation;
 local stability; \hfill\break\indent Hopf-bifurcation; diffusion-driven instability;
 spatial pattern formation}

\begin{abstract}
 In this article, we propose a three dimensional mathematical model of
 phytoplankton dynamics with the help of reaction-diffusion equations
 that studies the bifurcation and pattern formation mechanism.
 We provide an analytical explanation for  understanding phytoplankton 
 dynamics with three population classes:  susceptible, incubated, and infected.
 This model has a Holling  type II response  function  for the  population 
 transformation from susceptible to incubated class in an aquatic ecosystem. 
 Our main goal is to provide a  qualitative  analysis of Hopf bifurcation 
 mechanisms, taking death rate of infected  phytoplankton as bifurcation 
 parameter, and to study further spatial patterns formation  due to spatial 
 diffusion. Here analytical findings are supported by the results of numerical
 experiments. It is observed that the coexistence of all  classes of 
 population depends on the rate of diffusion.  Also we obtained the 
 time evaluation pattern formation of the spatial system.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks

\section{Introduction}

It is well known that the phytoplankton and zooplankton are single cell organisms 
that drift with the currents on the surface of open oceans. 
Further, the phytoplanktons are the staple items for the food web and they are 
the recycler of most of the energy that flows through the ocean ecosystem. 
It has a major role in stabilizing the environment and survival of living 
population as it consumes half of the universal carbon-dioxide and releases oxygen.
 So far, there is a number of studies which show the presence of pathogenic viruses 
in the plankton community \cite{dubey2002model, petrovskii2001wave}.
 A good review of the nature of marine viruses and their ecological as well as 
their biological effects is given in \cite{pozio1980behaviour}. Some researchers 
have shown using an electronic microscope that these viral diseases can affect 
bacteria and phytoplankton in coastal area  and viruses are held responsible 
for the collapse of Emiliania huxleyi bloom in
 Mesocosms \cite{gourley1996instability, ruan1998global}.
 Marine viruses infect not only plankton but cultivated stocks of Crabs,  
Oysters,  Mussels,  Clams shrimp,  Salmon and Catfish,  etc.,  are all susceptible 
to various kinds of viruses. We know that the viruses are nonliving organisms,  
in the sense,  they have no metabolism when out side the host and they can reproduce 
only by infecting the living organisms.  Viral infection of the phytoplankton 
cell is of two types,  namely,  Lysogenic and Lytic. Moreover from literature,  
in lytic viral infection,   when a virus injects its DNA into a cell,  
it hijacks the cell's replication machinery and produces large a number of viruses. 
As a result,  they rupture the host and are released into the environment. 
On the other hand,  in lysogenic viral infection,  the DNA of the viruses 
do not use the machinery of the host themselves,  but their genes are duplicated 
each time as the host cell divides. Many papers have already been developed which 
have used this kind of lysogenic viral infection 
\cite{hilker2006oscillations,malchow2005spatiotemporal, Petrovskii2001, singh2004role}.  
In this kind of infection,  infected phytoplankton do not grow like susceptible 
phytoplankton but their number grows and it is proportional to the number of 
susceptible phytoplankton and infected phytoplankton. Moreover, we do not look on
 the viruses as mere pathogens that destroy others life rather,  they  produce the fuel,  
essential to the running of the marine engine by destroying phytoplankton;
  i.e.,  they produce the essential minerals which are required for the growth of
 phytoplankton. Secondly,  we have introduced incubated class in between a susceptible 
class and  an infected class,  Unlike simple S-I-S models  in  mathematical epidemiology. 
Generally,  we have been using direct shifting from susceptible to infected class,  
but this process is not regular,  rather,   phytoplankton stay for some definite 
period in an incubated class after leaving the susceptible class and joining  
the infected class. The period for which they stay in incubated class is termed as the 
incubation period. The incubation period is defined as the time from the exposure 
to the onset of disease, when they are exposed to infection. 
The incubation period is useful not only for making the rough guesses for finding the 
cause and source of infection,  but also in developing  treatment strategies to extend 
the incubation period,  and for performing  an early projection of  the disease 
prognosis \cite{pozio1980behaviour,Satnuoian2000}.

In the previous thirty years,  the pattern formation in predator-prey systems has 
been studied by many researchers. The spread of diseases in human populations can
 exhibit large scale patterns, underlining the need for spatially explicit approaches. 
The spatial component of ecological interactions has been identified as an important 
factor in a spatial world and it is a natural phenomenon that a substance goes from 
high density regions to low density regions. Also, spatial patterns are ubiquitous 
in nature. These patterns modify the temporal dynamics and stability properties 
of population densities at a range of spatial scales and their effects must be 
incorporated in temporal ecological models that do not represent space explicitly. 
The spatial component of ecological interactions has been identified as an important
factor in how ecological communities are shaped 
\cite{medvinsky2002spatiotemporal,Murray2003}. Patterns are present in the chemical 
and biological worlds and since Turing \cite{Turing1952a} first introduced his model 
of pattern formation, reaction-diffusion equations have been a primary means of 
predicting them. Similarly structured systems of ordinary differential equations
 govern the spatiotemporal dynamics of ecological population models, yet most of 
the simple models predict spatially homogeneous population distributions 
\cite{liu2009pattern}. It is important to note that the above model takes into 
account the invasion of the prey species by predators but does not include stochastic 
effects or any influences from the environment. Nevertheless, a reaction-diffusion
 equation modeling predator-prey interaction show a wide spectrum of ecologically 
relevant behavior resulting from intrinsic factors alone, and is an intensive area 
of research. An introduction to research in the application of reaction-diffusion 
equations to population dynamics are available in  
\cite{morozov2004bifurcations,lee2004pattern}.

 These numerical simulations reveal that a large variety of different spatiotemporal 
dynamics can be found in this model. The numerical results are consistent with our 
theoretical analysis. It should be noted that, if considered in a somewhat broader 
ecological perspective, our results have an intuitively clear meaning. There has been 
a growing understanding during the past years regarding the dynamics of real ecosystems.
 It is important to reveal different spatial dynamical regimes arising as a result of 
perturbation of the system parameters \cite{sun2009predator, xiao2006chaotic}.

The objective of the paper is to consider the Bifurcation and spatial pattern formation 
in a Phytoplankton dynamics.  In section 2 and 3, we have developed mathematical
 model and analyzed dynamical behavior of system. In section 4, studied the
Hopf bifurcation mechanisms analytical and numerically, in section 5, studied the 
spatial pattern formation and numerical simulation of two dimensional systems and 
finally in section 6, we summarize our results and discuss the relative importance 
of different mechanisms of bifurcation and spatial pattern formation.

\section{Mathematical Model}

Now,  we consider the case of viral infection of phytoplankton population,  
the shifting from susceptible to infected class is not regular. In fact, 
 the susceptible phytoplankton stay for some definite period in the incubated 
class after leaving the susceptible class and joining the infected class. 
Taking the population densities of susceptible and infected phytoplankton  
as $P_s$ and $P_i$ respectively,  at any instant of time $T$. The population 
of susceptible phytoplankton is assumed to be growing logistically with 
intrinsic growth rate $r$ and carrying capacity $K$. Now  taking  $P_{in}$ as 
the population density of population in incubated class.  Here, we will use
 nonlinear Holling Type II functional responses for disease spreading because 
the disease conversion rates become saturated as victim densities increase. 
Let $\alpha$ be the disease contact rate and it is volume-specific encounter 
rate between susceptible and infected phytoplankton, which is equivalent to 
the inverse of the average search time between successful spreading of disease. 
The coefficients $\delta$ and $\beta$ are the total removal of phytoplankton 
from the infected and incubated class because of the death (including recovered) 
from disease and due to natural causes respectively. Again, $\gamma_1$ be the 
fraction of the population recovered from infected phytoplankton population 
and joined in the susceptible phytoplankton population and $\beta_1$ is the 
fraction of the incubated class population which will move to the infected class. 
Therefore, quantitatively $\delta> \gamma_1$ and $\beta > \beta_1$.
 Keeping in view  the above,  the mathematical model for a viral infected 
phytoplankton population with an incubated class is governed by the following 
set of differential equations:
\begin{gather}
\frac{dP_s}{dT} = rP_s\big(1-\frac{P_s}{K}\big)-\frac{\alpha P_sP_i}{P_s+1}
+\gamma_1 P_i,  \label{c2in3}
\\
\frac{dP_{in}}{dT}= \frac{\alpha P_sP_i}{P_s+1} - \beta P_{in}, \label{c2in4}
 \\
\frac{d P_i}{dT} =  \beta_1P_{in} - \delta P_{i}, \label{c2in5} 
\\
P_s(0)>0,  \quad P_i(0)>0,  \quad  P_{in}(0)>0.
\end{gather}
The  Holling type-II the functional response $\frac{\alpha P_sP_i}{P_s+1}$ 
is used \cite{Holling1959} and many other researchers. In this section,  
we have studied the dynamical behavior of the system \eqref{c2in3}-\eqref{c2in5},  with
initial population;  i.e.,  $P_s(0) > 0$,  $P_{in}(0) > 0$,  $P_i(0) > 0$
and the total population at any instant $t$ is $N(t)=P_s(t)+P_{in}(t)+P_i(t)$.

Now,  on rescaling the above system \eqref{c2in3}-\eqref{c2in5} using  the change of
variables: 
$ x=\frac{P_s}{K}$,
$y=\frac{P_{in}}{K}$,  $z=\frac{P_i}{K}$,  $t = r T$,   we obtain
\begin{gather}
\frac{dx}{d t} = x (1 - x) - \frac{a x z}{x+1} + c z,  \label{c2in6} 
\\
\frac{dy}{d t} =  \frac{a x z}{x+1} - d y,  \label{c2in7} 
\\
\frac{dz}{d t} =  d_1 y - e z,  \label{c2in8} 
\end{gather} 
where $ a =\frac{\alpha K}{r}$,  $c = \frac{\gamma_1}{r}$,  
 $d =\frac{\beta}{r}$,  $ d_1 = \frac{\beta_1}{r}$,   $e = \frac{\delta}{r}$, 
$x(0) > 0$,  $y(0) > 0$ and $z(0) > 0$.

\section{Boundedness of the system}

Now,  we will study the existence of all possible steady states of the system and 
the boundedness of the solutions.
There are three biologically feasible equilibriums for the 
system \eqref{c2in6}-\eqref{c2in8},  (i) $E_0 = (0,  0,  0)$ is the trivial 
steady state; (ii) $E_1= (1,  0,  0)$ is the disease free steady
state and (iii) $E^*= (x^*,  y^*,  z^*)$ is endemic equilibrium
state,  where
$$ x^* = \frac{e d}{d_1 a - e d}, \quad
 y^* = \frac{e^2 d (a d_1 - 2 e d)}{(e d - d_1 c) (a d_1 - e d)^2}, \quad
z^* = \frac{e d d_1(a d_1 - 2 e d)}{(e d - d_1 c) (a d_1 - e d)^2}.
$$
 Further,   it is clear from the above expression that $E^* \in R_{+}^3$,  if
 $a > \frac{d e}{d_1} > c$. Now,  we will show that all the solutions of the
system \eqref{c2in6}-\eqref{c2in8} are bounded in a region
 $ B \subset R_{+}^3$. We consider the  function
\begin{equation}
 w(\tau) = x(\tau)+y(\tau)+z(\tau) \label{newin1},  
\end{equation}
and substituting the values from \eqref{c2in6}-\eqref{c2in8},  we obtain
\[ 
\frac{d w}{d \tau} = x(1-x) - (d-d_1) y - (e-c)z.
\]
If we choose a positive real number $\eta = min\{d-d_1,  e-c\}$,  then
\[ 
\frac{d w(\tau)}{d \tau} + \eta w(\tau) \le x(1+\eta)- x^2 =f(x).
\]
Moreover,  $f(x)$ has maxima at $x=(1+\eta)/2$ and 
$f(x) \le (1+\eta)^2/4= M$ (say). Hence, 
 $ \dot{w}(\tau) + \eta w(\tau) \le M$. Now,  using comparison theorem,  
as $\tau \to \infty$,  we obtain  
$ sup \ w(\tau) = \frac{M}{\eta}$. Therefore,  
$ 0 \le x(\tau)+y(\tau)+z(\tau) \le \frac{M}{\eta}$.
and let us consider the set $ B = \{ (x, y, z) \in R_{+}^3 : 0
\le x(\tau)+y(\tau)+z(\tau) \le {M}/{\eta} \}$,  hence, 
 The system \eqref{c2in6}-\eqref{c2in8} is
uniformly bounded in the region $B \subset R_{+}^3$.

\subsection{Local stability of the system}

We have already established that the system \eqref{c2in6}-\eqref{c2in8} has three
equilibrium points, namely,
$E_0=(0, 0, 0)$,  $E_1=(1, 0, 0)$ and $E^*=(x^*, y^*, z^*)$ in the previous section. 
The general variational matrix
corresponding to the system is given by
\[
J = \begin{bmatrix}
 & 1- 2 x^* -\frac{a z^*}{(x^*+1)^2} & 0 & -\frac{a x^*}{(x^*+1)} +c \\
 &\frac{a z^*}{(x^*+1)^2}  & -d  & \frac{a x^*}{(x^*+1)} \\
 & 0   & d_1   & -e
\end{bmatrix}
\]
Now,  corresponding to the trivial steady state $E_0=(0, 0, 0)$ the
Jacobian J has the following eigenvalues $\lambda_i = 1,   -d,  -e$.
 Hence,  $E_0$ is repulsive in $x$-direction and attracting in
$y-z$ plane. Clinically,  it means  that when there is no susceptible
population then,  there will be no mass in incubated and in the infected
class. Hence,  $E_0$ is a saddle point. The  corresponding to
the disease free equilibrium point $E_1 = (1, 0, 0)$,  we have
following  $\lambda_1 = -1$ and $\lambda_{2, 3}$ are
the roots of the  quadratic equation
$ \lambda^2 + (d +e) \lambda + (d e - \frac{a d_1}{2}) = 0$
when $d e > d_1 a/2$,  then both the roots  have a negative real parts
and thus $E_1(1, 0, 0)$ is locally stable in this case.

Further,  from the existence of $E^*$ and the stability condition of $E_1$,  it is clear
that the instability of the disease free equilibrium will lead to
the existence of the endemic equilibrium. Now,  we will examine the
local behavior of the flow of the system around the endemic
equilibrium $E^*$. The characteristic equation corresponding to the
equilibrium is
\begin{equation}
P(\lambda) = \lambda^3 + A_1 \lambda^2 + A_2 \lambda + A_3 = 0\label{c2in9},  
\end{equation}
 where
\begin{gather*}
A_1 =  2 x^* + \frac{a z^*}{(x^*+1)^2} + d +e -1, \\
 A_2 =  2 x^* e+2 x^* d +e d-e-d-\frac{ad_1x^*}{(x^*+1)}+\frac{aez^*}{(x^*+1)^2}
 +\frac{adz^*}{(x^*+1)^2}, \\
 A_3 =  \frac{eadz^*}{(x^*+1)^2}+\frac{ad_1x^*}{(x^*+1)}-\frac{2ad_1{x^*}^2}{(x^*+1)}
 -\frac{acd_1z^*}{(x^*+1)^2}+2 x e d- e d.
\end{gather*}
Hence, using Routh-Hurwitz criteria $E^*$ is locally asymptotically stable,  
if $A_i> 0$,  $i=1, 2, 3$ and $A_1 A_2 > A_3$.

\section{Hopf-bifurcation analysis}

In this section,  we obtain the Hopf-bifurcation criteria of
above system \eqref{c2in6}-\eqref{c2in8},  taking "$e$" 
(i.e.,  total removal of phytoplankton from the infected population ) 
as the bifurcation parameter. Now,  the necessary
and sufficient condition for the existence of the Hopf-bifurcation can be obtained 
if there exists  $e=e_0$,   such that (i) $A_i(e_0)> 0$,  $i=1, 2, 3$, 
 (ii) $A_1(e_0) A_2(e_0) - A_3(e_0) = 0$ and 
(iii) $\frac{d}{d e} (u_i) \ne 0$,  $i=1, 2, 3$,
where $u_i$ is the real part of the eigenvalues of the characteristic
 equation \eqref{c2in9}. After substituting of the values of $A_i$  for $i=1, 2, 3$ 
in equation $A_1 A_2 - A_3=0$ and solving it for $e$,  it reduces to
\begin{equation}
e^7 B_1+e^6 B_2+e^5 B_3+e^4 B_4+e^3 B_5+e^2 B_6+e B_7 +B_8 = 0\label{c2in10}, 
 \end{equation}
 where
\begin{gather*}
B_1 =  (d^6-2d^5d_1a+8d^6d_1a),  \\
\begin{aligned}
B_2 &=  (2d^4{d_1}^2a^2+2d^4{d_1}^2ac-2d^6d_1a+2d^5d_1a+d^7-2d^5d_1c-12d^5{d_1}^2a^2\\
   &\quad -16d^5{d_1}^2ac+8d^2d_1a),
\end{aligned} \\
\begin{aligned}
B_3 &=  (2d^5{d_1}^2ac-2d^4d_1a+9d^3{d_1}^3a^2c+2d^6d_1a-4d^4{d_1}^2ac-2d^4{d_1}^2a^2\\
   & \quad+4d^3{d_1}^3a^3-d^5{d_1}^2a^2+d^4{d_1}^2c^2+4d^2{d_1}^2ac-2d^6d_1c+4d^4{d_1}^3a^3\\
   &\quad +8d^4{d_1}^3ac^2+24d^4{d_1}^3a^2c-12d^6{d_1}^2a^2-16d^6{d_1}^2ac ), 
\end{aligned}\\
\begin{aligned}
B_4 &=  (20d^4{d_1}^3a^2c+2d^6{d_1}^2ac+2d^6{d_1}^2a^2-4d^5{d_1}^2a
     c-8d^5{d_1}^2a^2-22d^3{d_1}^3a^2c\\
 &\quad +2d^4{d_1}^3ac^2+4d^4{d_1}^3a^3-2d{d_1}^4a^3c-3d^2{d_1}^4a^2c^2+6d^5{d_1}^2a^2-2d^3{d_1}^3ac^2\\
 &\quad +d^5{d_1}^2c^2-2d^3{d_1}^3a^3+4d^5{d_1}^2ac-12d^3{d_1}^4a^2c^2-8d^3{d_1}^4a^3c\\
 &\quad +4d^5{d_1}^3a^3+24d^5{d_1}^3a^2c+8d^5{d_1}^3ac^2),
\end{aligned} \\
\begin{aligned}
B_5 &=  (9d^5{d_1}^3a^2c-12d^4{d_1}^3a^2c+10d^4{d_1}^3a^3-12d^2{d_1}^4a^3+4d{d_1}^5a^3c^2-12d^3{d_1}^4a^3c\\
&\quad +7d^2{d_1}^4a^2c^2-8d^4{d_1}^3a^2c-4d^4{d_1}^3a^3+8d{d_1}^4a^3c-d{d_1}^4a^4-3d^3{d_1}^4a^2c^2\\
&\quad +d{d_1}^5a^4c-2d^4{d_1}^3ac^2+d^2{d_1}^4a^4
 +4d^2{d_1}^5a^3c^2-8d^4{d_1}^4a^3c-12d^4{d_1}^4a^2c^2), 
\end{aligned}\\
\begin{aligned}
B_6 &=  (d^4{d_1}^4a^4-4d^3{d_1}^4a^3c-2d^3{d_1}^3a^4-12d^2{d_1}^5a^3c^2
 -2d^5{d_1}^4a^2c^2-6d{d_1}^5a^3c^2\\
&\quad +d^2{d_1}^4a^3c+11d^3{d_1}^4a^2c^2+4d^2{d_1}^5a^4c-d^4{d_1}^4a^2c^2-d^6{d_1}^4c),
\end{aligned} \\
B_7 =  (4d^3{d_1}^5a^3c^2-6d^2{d_1}^5a^3c^2-d{d_1}^6a^4c^2+d^3{d_1}^5a^4c+{d_1}^6a^4c^2), \\
B_8 =  (d{d_1}^6a^4c^2-d^2{d_1}^6a^4c^2).
\end{gather*}

For better understanding of the above result, taking an example values of the 
parameters $c=0.01$,  $d=0.11$,  $d_1=0.1$ and $a=7.15$,  we obtain
a positive root $e=0.074$ of the quadratic equation \eqref{c2in10}.
Therefore,  the eigenvalues of the characteristic equation
\eqref{c2in10} at $e=0.074$ are of the form 
$\lambda_{1, 2} = \pm iv$ and $\lambda_3 = - w$,  where $v$ and $w$ are positive real
numbers.

Now,  we will verify the condition (iii) of Hopf-bifurcation. Put $\lambda = u + iv$ in
\eqref{c2in9},  we obtain 
\begin{equation}
(u + iv)^3 + A_1 (u + iv)^2 + A_2 (u + iv) + A_3 = 0. \label{c2in11}
 \end{equation}
On separating the real and imaginary
parts and eliminating $v$ between real and imaginary parts,  we obtain
\begin{equation}
8 u^3 + 8 A_1 u^2 + 2(A_1^2 + A_2) u + A_1 A_2 - A_3 = 0.
\label{c2in12} 
\end{equation}
Now,  we have $u(e_0)=0$ as $A_1(e_0) A_2(e_0) - A_3(e_0) = 0$. Further,  
 $e=e_0$, is the only positive root of $A_1(e_0) A_2(e_0) - A_3(e_0) = 0$,  
and the discriminate of $8 u^2 + 8 A_1 u + 2(A_1^2 + A_2)= 0$ is
 $64A_1^2 -64(A_1^2+A_2) < 0$. Here,
differentiating \eqref{c2in12} with respect to $e$,  we have
$ \left(24 u^2 + 16 A_1 u + 2(A_1^2 +A_2)\right) \frac{d u}{d e}
 + \left(8 u^2 + 4 A_1 u \right) \frac{d A_1}{d e} + 2 u \frac{d A_2}{d e}
 + \frac{d}{d e}(A_1 A_2 - A_3)= 0$.
Now,  since at $e=e_0$,  $u(e_0) = 0$,  we obtain
$$ 
\big[ \frac{d u}{d e} \big]_{e=e_0} 
= \frac{-\frac{d}{d e}(A_1 A_2 - A_3)}{2(A_1^2 +A_2)} \ne 0,
$$
 which ensures that the above system has a Hopf-bifurcation. 
It is shown graphically in figure \ref{f21}.

\begin{figure}
\begin{center}
\includegraphics[width=0.98\textwidth]{fig1}
%\includegraphics[height=5.0in, width=6.0in]{fig1}
\caption{Phase plane representation of three species around the endemic
equilibrium,  taking $c =0.01$,  $d=0.1$1,  $d_1=0.1$,  $a=7.15$: 
(A) $e=0.073$,  (B) $e=0.074$,  (C) $e=0.075$,  (D) $e=0.076$}\label{f21}
\end{center}
\end{figure}

\section{Spatiotemporal model}

Now,  we consider the  phytoplankton dynamics  with movement (i.e.,  diffusion). 
Therefore the population densities;  i.e.,  $P_s$,  $P_i$ and $P_{in}$ 
become space and time dependent. Keeping in view of the above, 
 our mathematical model can stated by the  reaction diffusion system
\begin{gather}
\frac{\partial u}{d t} = u (1 - u) - \frac{a u w}{u+1} + c w+ D_a \nabla^2 u, \label{cin1}
\\
\frac{\partial v}{d t} =  \frac{a u w}{u+1} - d v + D_b \nabla^2 v, \label{cin2}
\\
\frac{\partial w}{d t} =  d_1 v - e w+ D_c \nabla^2 w,\label{cin3}
\end{gather}
where $(x, y)$ is the position in space for two dimensional on a bounded domain  
and $D_a$,  $D_b$ and $D_c$ are diffusion  coefficients for susceptible phytoplankton, 
incubated phytoplankton and infected phytoplankton population,  respectively. 
The no-flux boundary conditions were used.

Now, we will explore the possibility of diffusion-driven instability with respect 
to the steady state solution, i.e., the spatially homogenous solution 
$(u^*,  v^*,  w^*)$ of the reaction diffusion system.

We assume that $E^*$ is stable in the temporal system and unstable in spatiotemporal 
system, which means that the spatially homogeneous equilibrium is unstable with 
respect to spatially homogeneous perturbations.

We obtain the conditions for the diffusion instability to occur in system 
\eqref{cin1}-\eqref{cin3}, one should check how small heterogeneous perturbations 
of the homogeneous steady state develop in the large-time limit. For this purpose, 
we consider the  perturbation
\begin{gather}
u(x,y,t) = u^* + \epsilon \exp((kx+ky)i+\lambda_kt),\label{cin10}
\\
v(x,y,t) = v^* + \eta \exp((kx+ky)i+\lambda_kt),\label{cin11}
\\
w(x,y,t) = w^* + \rho \exp((kx+ky)i+\lambda_kt),\label{cin12}
\end{gather}
where $\epsilon$, $\eta$ and $\rho$ are chosen to be small and $k = (k_x,k_y)$ 
is the wave number.
Substituting \eqref{cin10}-\eqref{cin12} into \eqref{cin1}-\eqref{cin3}, 
linearizing the system around the interior equilibrium $E^*$, we obtain 
the characteristic equation as follows:
\begin{equation}
|J_k-\lambda_k I_2|=0, \label{rs11}
\end{equation}
with
\[
J_k  =
\begin{bmatrix}
 1-{D_a} k^2-2 u^*+\frac{a u^* w^*}{(1+u^*)^2}-\frac{a w^*}{1+u^*} 
 & 0 & c-\frac{a u^*}{1+u^*} \\
 -\frac{a u^* w^*}{(1+u^*)^2}+\frac{a w^*}{1+u^*} & -(d+{D_b} k^2) 
 & \frac{a u^*}{1+u^*} \\
 0 & {d_1} & -(e+{D_c} k^2)
\end{bmatrix},
\]
where  $I_3$ and k are third order identity matrix and  wave number respectively.
The characteristic equation following form:
\begin{equation}
\lambda^3+P_2 \lambda^2+P_1 \lambda+P_0=0, \label{rs03}
\end{equation}
 where
\[
P_2= -\Big(1-d-e-{D_a} k^2-{D_b} k^2-{D_c} k^2-2 u^*
+\frac{a u^* w^*}{(1+u^*)^2}-\frac{a w^*}{1+u^*}\Big),
\]
\begin{align*}
P_1 &= -\Big(d+e-d e-d D_a k^2+D_b k^2+D_c k^2-d D_c k^2-D_a e k^2-D_b e k^2\\
&\quad -D_a D_b k^4 -D_a D_c k^4-D_b D_c k^4-2 d u^*-2 e u^*-2 D_b k^2 u^*-2 D_c k^2 u^*\\
&\quad +\frac{a {d_1} u^*}{(1+u^*)} +\frac{a d u^* w^*}{(1+u^*)^2}
 +\frac{a e u^* w^*}{(1+u^*)^2}+\frac{a D_b k^2 u^* w^*}{(1+u^*)^2}
 +\frac{a D_c k^2 u^* w^*}{(1+u^*)^2}\\
&\quad -\frac{a d w^*}{1+u^*}-\frac{a e w^*}{1+u^*}
 -\frac{a D_b k^2 w^*}{1+u^*}-\frac{a D_c k^2 w^*}{1+u^*}\Big),
\end{align*}
\begin{align*}
P_0 &= \Big(-d e-d D_c k^2+d D_a e k^2-D_b e k^2+d D_a D_c k^4-D_b D_c k^4+D_a D_b e k^4
\\
& \quad +D_a D_b D_c k^6+2 d e u^*+2 d D_c k^2 u^*+2 D_b e k^2 u^*+2 D_b D_c k^4 u^*
 +\frac{a {d_1} u^*}{1+u^*}\\
& \quad -\frac{a {d_1} D_a k^2 u^*}{1+u^*}-\frac{2 a {d_1} u^{*2}}{1+u^*}+\frac{a c {d_1} u^* w^*}{(1+u^*)^2}-\frac{a d e u^* w^*}{(1+u^*)^2}
 -\frac{a d D_c k^2 u^* w^*}{(1+u^*)^2}\\
 & \quad -\frac{a D_b e k^2 u^* w^*}{(1+u^*)^2}-\frac{a D_b D_c k^4 u^* w^*}{(1+u^*)^2}
 -\frac{a c {d_1} w^*}{1+u^*}+\frac{a d e w^*}{1+u^*}+\frac{a d D_c k^2 w^*}{1+u^*}\\
& \quad +\frac{a D_b e k^2 w^*}{1+u^*}+\frac{a D_b D_c k^4 w^*}{1+u^*}\Big).
\end{align*}
According to the Routh-Hurwitz criterium all the eigenvalues have negative real
 parts if and only if the following conditions hold:
\begin{gather}
P_2>0,  \label{fgh1} \\
P_0>0, \label{fgh2} \\
Q=P_0-P_2P_1<0. \label{fgh3}
\end{gather}
This is best understood in terms of the invariants of the matrix and of its 
inverse matrix
\[
J_k^{-1}= \frac{1}{\det (J_k) }
\begin{pmatrix}
 M_{11} & M_{12} & M_{13} \\
M_{21} & M_{22} & M_{23} \\
 M_{31} & M_{32} & M_{33}
\end{pmatrix},
\]
where  $M_{11}=de-\frac{ad_1u}{(u+1)}$,  $M_{12}=-d_1c-\frac{ad_1u}{u+1}$,  
$M_{13}=d(c-\frac{au}{u+1})$, $M_{21}=\frac{aew}{(u+1)^2}$,  
$M_{22}=-e(1-2u-\frac{aw}{(u+1)^2})$,  
$M_{23}=-(\frac{au}{(u+1)}(1-2u-\frac{aw}{(u+1)^2})
 -\frac{aw}{(u+1)^2}(c-\frac{au}{(u+1)}))$,
$ M_{31}=\frac{ad_1w}{(u+1)^2}$,  $M_{32}=-d_1(1-2u-\frac{aw}{(u+1)^2})$,  
$M_{33}=-d(1-2u-\frac{aw}{(u+1)^2})$. 
Here, matrix $M_{ij}$ is the adjunct of $J_k$.

We obtain the following conditions of the steady-state stability 
(i.e. stability for any value of  k):\\
(i) All diagonal cofactors of matrix $J_{k}$ must be positive.\\
(ii) All diagonal elements of matrix $J_k$ must be negative. 


The two above condition taken together are sufficient to ensure stability of a 
give steady state. It means that instability for some $k>0$ can only be observed 
if at least one of them is violated. Thus we arrive at the
following  necessary condition for the Turing Instability \cite{Qian2003}:
(i) The largest diagonal element of matrix $J_{k}$ must be positive and/or (ii) 
the smallest diagonal cofactor of matrix $J_k$ must be negative.

By the Routh-Hurwitz criteria, instability takes place if and only if 
one of the conditions \eqref{fgh1}-\eqref{fgh3} is broken.
 We consider \eqref{fgh2} for instability condition:
\begin{align*}
P_0(k)&= D_aD_bD_ck^6-(D_aD_ba_{33}+D_bD_ca_{11}+D_aD_ca_{22})k^4\\
& \quad+(D_aM_{11}+D_bM_{22}+D_3M_{33})k^2-\det J.
\end{align*}

According to Routh-Hurwitz criterium $P_0(k^2)<0$ is sufficient condition for 
matrix $J_k$ being unstable. Let us assume that $M_{33}<0$. 
If  we choose $D_a=0, D_b=0$ and
\begin{equation}
\begin{aligned}
P_0(k^2)&= D_c M_{33}k^2-det(J_k)\\
 &= -\Big[ d D_c k^2 \Big(1-2 u^*-\frac{a w^*}{(1+u^*)^2}\Big)
+\frac{1}{(1+u^*)^2}(1+u^*) (2 u^*-1)\\
&\quad\times \Big(ed+D_cdk^2 +edu^*+D_cdu^*k^2-ad_1u^*\Big)\\
&\quad +(edaw^*+D_cdaw^*k^2-cad_1)\Big]<0.
\end{aligned}\label{insc}
\end{equation}
Hence, in this system diffusion-driven instability occur.

Now, we  obtain  the  eigenvalues of the characteristic equation \eqref{rs03}
 numerically of the spatial system \eqref{cin1}-\eqref{cin3}. 
Here, we choose  some parametric values of $a=0.4$,  $c=0.01$,  $d=0.11$,  $d_1=0.1$, 
 $e=0.08$. In this set of  values $P_0(k^2)<0$, for all $k>0$, hence 
from \eqref{insc}, we can observe  diffusion driven instability of the system
(see Fig.\ref{rs12}).

\begin{figure}
\begin{center}
\includegraphics[width=0.98\textwidth]{fig2}
%\includegraphics[height = 5.0in,width = 6.0in]{diffusivityfig}
\caption{Plot of max Re$(\lambda(k))$ against $k$. 
The other parametric values are given in text}\label{rs12}
\end{center}
\end{figure}

\subsection{Pattern formation}
Now, we will study numerical system \eqref{cin1}-\eqref{cin3} 
 for the pattern formation of  two dimensional space  with zero-flux boundary
 conditions is used.  We choose the initial spatial distributions of each 
species are random and the numerical results are obtained using finite difference
 method  for figure \ref{f27}-\ref{f28} and we use the parametric values  same 
as above section.

Now, we obtain  that  the spatial distributions of phytoplankton dynamics in 
the time evaluation in figures \ref{f27}-\ref{f28}. By varying coupling parameters, 
 we observed that one parameter value change  then spatial structure change 
over the times of the  spatial system. In figure \ref{f27}-\ref{f28} have  
observed well organized structures for the spatial distribution population 
also observed that time $T$ increase from $10$ to $300$ the density  
of different classes of population become uniform throughout the space. 
Finally,  all  these figure are shown that the qualitative changes spatial 
density distribution of the spatial system for the each species.

\begin{figure}
\begin{center}
\includegraphics[width=0.98\textwidth]{fig3}
\caption{Spatial distribution of susceptible phytoplankton [first column],  
incubated phytoplankton [second column] and infected phytoplankton 
[third column] population density  of the model system 
\eqref{cin1}-\eqref{cin3}. The diffusivity coefficient values are:
$D_a=0.02$,  $D_b=0.2$,  $D_c=5$. Spatial patterns are obtained
at different time levels: for  plot $T =10$ (a, b, c), $T = 40$
 (d, e, f), $T = 80$ (g, h, i)}\label{f27}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=0.98\textwidth]{fig4}
\caption{Spatial distribution of susceptible phytoplankton [first column],  
incubated phytoplankton [second column] and infected phytoplankton  
[third column] population density  of the model system 
\eqref{cin1}-\eqref{cin3}. The diffusivity coefficient  values are:
$D_a=0.02$,  $D_b=0.2$,  $D_c=5$. Spatial patterns are obtained
at different time levels: for  plot $T =100$ (a, b, c), $T = 200$ (d, e, f), 
$T = 300$ (g, h, i)}\label{f28}
\end{center}
\end{figure}

\subsection*{Conclusion}
In this paper, we have studied a phytoplankton dynamics with viral infection. 
We observed that in a  three dimensional phytoplankton system with the Holling-II 
response function, there exit Hopf bifurcation with respect to remove rate 
including nature death   of infected phytoplankton.  In the qualitative analysis,  
we  studied the boundedness,  dynamical behavior and  local stability of the system.  
It is established that the rate of total removal of phytoplankton from the infected 
class;  i.e.,  $e$,  crossed its threshold value,  $e=e_0$,  then susceptible,  
incubated and infected phytoplankton population started  oscillating around the 
endemic equilibrium. The above result has been shown numerically in 
figure \ref{f21} for different values of $e$. In particular, 
 in figure \ref{f21}(A),  we observed that the endemic equilibrium was stable,  
when $e < 0.074$,  but when it crossed the threshold value of $e = 0.074$,  
the above system showed  Hopf-bifurcation,  shown in figure \ref{f21}(B,  C,  D). 
We have also observed  spatially ordered structures of patterns in spatial 
systems and the solutions of the spatial system converges to its  equilibrium 
as time $T$ increase from 10 to 300 in the two-dimensional pattern formation,  
shown in figures \ref{f27}-\ref{f28}. It is numerically established that with 
slight change in a time $T$ parameter of the system \eqref{cin1}-\eqref{cin3}, 
can lead to dramatic changes in the qualitative behavior of the system.


\begin{thebibliography}{10}

\bibitem{dubey2002model}
B.~Dubey, B.~Das,  J.~Hussain;
\emph{A model for two competing species with
  self and cross-diffusion}, Indian Journal of Pure and Applied Mathematics
  \textbf{33} (2002), no.~6, 847--860.

\bibitem{gourley1996instability}
S.A.~Gourley;
\emph{Instability in a predator-prey system with delay and spatial
  averaging}, IMA journal of applied mathematics \textbf{56} (1996), no.~2,
  121--132.

\bibitem{hilker2006oscillations}
F. M. Hilker, H.~Malchow, M.~Langlais, S.V. Petrovskii;
\emph{Oscillations
  and waves in a virally infected plankton system:: Part ii: Transition from
  lysogeny to lysis}, Ecological complexity \textbf{3} (2006), no.~3, 200--208.

\bibitem{Holling1959}
C.~Holling;
 \emph{Some characteristics of simple types of predation and
  parasitism}, Can. Entomol. \textbf{91} (1959), 385 -- 398.

\bibitem{liu2009pattern}
P. P. Liu, Z.~Jin,
 \emph{Pattern formation of a predator--prey model},
  Nonlinear Analysis: Hybrid Systems \textbf{3(3)} (2009), no.~3, 177--183.

\bibitem{malchow2005spatiotemporal}
H.~Malchow, F. M. Hilker, R. R. Sarkar, K.~Brauer;
\emph{Spatiotemporal
  patterns in an excitable plankton system with lysogenic viral infection},
  Mathematical and computer modelling \textbf{42} (2005), no.~9-10, 1035--1048.

\bibitem{medvinsky2002spatiotemporal}
A. B. Medvinsky, S. V. Petrovskii, I. A. Tikhonova, H.~Malchow,  B.L. Li,
  \emph{Spatiotemporal complexity of plankton and fish dynamics}, Siam Review
  \textbf{3(44)} (2002), 311--370.

\bibitem{morozov2004bifurcations}
A.~Morozov, S.~Petrovskii,  B. L. Li;
\emph{Bifurcations and chaos in a
  predator-prey system with the allee effect}, Proceedings of the Royal Society
  of London. Series B: Biological Sciences \textbf{271} (2004), no.~1546, 1407.

\bibitem{Murray2003}
J. D. Murray; 
\emph{Mathematical biology ii: Spatial models and biomedical
  applications, 3rd ed., in: Biomathematics}, Springer, New York, 2003.

\bibitem{Petrovskii2001}
S.~Petrovskii, H.~Malchow;
\emph{Wave of chos: new mechanism of pattern
  formation in spatiotemporal population dynamics.}, Theor. Pop. Bio.
  \textbf{59} (2001), 157 -- 174.

\bibitem{petrovskii2001wave}
S. V. Petrovskii, H.~Malchow;
\emph{Wave of chaos: new mechanism of pattern
  formation in spatio-temporal population dynamics}, Theoretical Population
  Biology \textbf{59} (2001), no.~2, 157--174.

\bibitem{pozio1980behaviour}
M. A. Pozio;
\emph{Behaviour of solutions of some abstract functional
  differential equations and application to predator-prey dynamics}, Nonlinear
  Analysis \textbf{4} (1980), no.~5, 917--938.

\bibitem{Qian2003}
H.~Qian, J.~D. Murray; 
\emph{A simple method of parameter space
  determination for diffusion-driven instability with three species}, Appl.
  Math. Lett. \textbf{14} (2003), 405--411.

\bibitem{ruan1998global}
S.~Ruan,  X. Z. He;
 \emph{Global stability in chemostat-type competition
  models with nutrient recycling}, SIAM Journal on Applied Mathematics
  \textbf{31} (1998), 170--192.

\bibitem{lee2004pattern}
H.~K.~Pak, S.~H.~Lee, H.~S. Wi; 
\emph{Pattern dynamics of spatial domains in
  three-trophic population model}, Journal Korean Physical Society
  \textbf{44(3)} (2004), no.~1, 656--659.

\bibitem{Satnuoian2000}
R. A. Satnuoian and M.~Menzinger;
\emph{Non-turing stationary pattern in
  flow-distributed oscillators with general diffusion and flow rates}, Physical
  Review. E \textbf{62(1)} (2000), 113--119.

\bibitem{singh2004role}
B. K. Singh, J.~Chattopadhyay, S.~Sinha;
 \emph{The role of virus infection
  in a simple phytoplankton zooplankton system}, Journal of theoretical biology
  \textbf{231} (2004), no.~2, 153--166.

\bibitem{sun2009predator}
G. Q. Sun, G.~Zhang, Z.~Jin,  L.~Li; 
\emph{Predator cannibalism can give rise
  to regular spatial pattern in a predator--prey system}, Nonlinear Dynamics
  \textbf{58} (2009), no.~1, 75--84.

\bibitem{Turing1952a}
A. M. Turing;
\emph{The chemical basis of morphogenesis}, Philos. Trans. R. Soc.
  London B \textbf{237} (1952), 37--72.

\bibitem{xiao2006chaotic}
J.~Xiao, H.~Li, J.~Yang,  G.~Hu; 
\emph{Chaotic turing pattern formation in
  spatiotemporal systems}, Frontiers of Physics in China \textbf{1} (2006),
  no.~2, 204--208.

\end{thebibliography}


\end{document}