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