Chaos and Pattern Formation in Spatial Phytoplankton Dynamics

In this paper the dynamics of spatially extended infected phytoplankton with the Holling type II functional response and logistically growing susceptible phytoplankton system is studied. The proposed model is an extension of temporal model available [6], in spatiotemporal domain. The reaction diffusion system exhibits spatiotemporal chaos in phytoplankton dynamics. This is particularly important for the spatially extended systems that are studied in this paper as they display a wide spectrum of ecologically relevant behavior, including chaos. In this system stability of the system is studied with respect to disease contact rate and the growth fraction of infected phytoplankton indirectly rejoin the susceptible phytoplankton population. The results of numerical experiments in one dimension and two dimensions in space as well as time series in temporal models are presented using MATLAB simulation. Moreover, the stability of the corresponding temporal model is studied analytically. Finally, the comparison of the three types of numerical experimentations are discussed in conclusion.


Introduction
Phytoplankton 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 as it consumes half of the universal carbondioxide and releases oxygen. So far, there is a number of studies which show the presence of pathogenic viruses in the plankton community [1,21,23,24]. A good review of the nature of marine viruses and their ecological as well as their biological effects is given in [7]. 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 [2,24]. 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. 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 [1,4,12,22]. So far, no one has used the lytic kind of viral infection in a plankton dynamical model. 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, in this paper we have introduced diffusion of population to capture the spatial distribution (i.e. pattern ) of both susceptible and infected class of population. Nevertheless, reaction-diffusion equations modeling predator-prey interactions show a wide spectrum of ecologically relevant behavior resulting from intrinsic factors alone [17], and is an intensive area of research. For an introduction to research in the application of reaction-diffusion equations to population dynamics, [11] and the references therein. This study is partially motivated by few works, namely, (i) a SIAM review paper [15] that considers the reaction -diffusion system as a model for marine plankton dynamics, (ii) a numerical experimentation of reaction -diffusion model in 1-D and 2-D space coordinates [14], leads to the formation of patterns, (iii) a study on diffusion induced chaos [18] and (iv) a phytoplankton dynamics with susceptible and infective classes [6]. Our mathematical model is an extension of temporal model presented by [6], in spatiotemporal domain. Analysis of the model using similar tools used in [6,14,15]. In this paper numerical experimentation show that, for various initial conditions, the evolution of the system leads to the formation of spiral patterns, followed by irregular patches over the whole domain (spatiotemporal chaos), which are in qualitative agreement with field observations.

Mathematical Model
A mathematical model [6], studied phytoplankton dynamics by considering the population densities of susceptible and infected phytoplankton as s P and i P 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 . Let  be the disease contact rate and  be the removal rate of the diseased phytoplankton population, out of which  fraction of infected phytoplankton rejoin the susceptible phytoplankton population, because, dead infected phytoplankton become nutrients for the growth of susceptible phytoplankton after bacterial decomposition and partially through immunization or natural recovery process in the ecosystem. Finally, their proposed mathematical model was as follows:

Model with Diffusion
In this paper, we will study the phytoplankton dynamics (1)-(2 with movement (i.e., diffusion) and since the phytoplankton population are not uniform in through out the habitat. Therefore the population densities, i.e., s P and i P are become space and time dependent. Keeping in view of the above, our mathematical model can stated by the following reaction diffusion equations: ) are diffusion coefficients and a , b , c , d are similar positive constants as in (1)- (2). Again, for more realistic interaction, we may choose, are the Holling type-II functions. Type II functional responses are the most frequently studied functional responses and well documented in empirical studies [9,13,19]. The shows the 'Holling type II functional response [10]. It is important to note that the above model takes into account the invasion of the phytoplankton dynamics. It is much easier to work with equation that have been scaled to nondimensional form, in above system, we take K and re-scaling the parameters via, . Hence our model reduces to where the parameters  ,  ,  ,  ,  are strictly positive. We assume the system is defined on a bounded domain (habitat), denoted by  , and augmented with appropriate initial conditions and the zero-flux boundary conditions. The particular choice of boundary conditions reflects the assumption that the individual species leave the domain. [8], the numerical solution of (5) -(6) for both dimensions and illustrate the simplicity of the numerical with MATLAB.
In order to provide guidelines on the appropriate choice of parameters for numerical simulation of the full reaction.diffusion system, it is important to consider the local dynamics of the system [15]. It is, naturally, the dynamics in the biologically meaningful region , corresponding to the coexistence of susceptible and infected phytoplankton, given by   . In order for the stationary point to be in the positive quadrant, we must have , The result leads to the stated restrictions on the parameters (assumed throughout). For certain parameter choices [15], the kinetics have a stable limit cycle surrounding the unstable stationary point ) , ( * * v u , i.e., the densities of susceptible and infected phytoplankton cycle periodically in time. An example of kinetics with a limit cycle is illustrated in Figure1.

Numerical Analysis of the System
Almost all of the realistic models in biology are nonlinear [16], without closed form solutions, thus numerical methods have an important part to play in investigating the behavior of their solutions. The rigorous numerical analysis of similar type of system were studied [8,14], using the standard Galerkin finite-element method [3,5] with piecewise linear basis functions. In the coming two subsections, we will study numerically our proposed system (5)-(6) for both dimensions with space coordinate respectively.

Numerical Solutions in One-Dimension
In Fig.2, the numerical solutions of the phytoplankton dynamics (i.e., u, v), are plotted with one space coordinate. The parameter values and initial data are given in the caption. The initial data are small spatial perturbations of the stationary solutions * u and * v of the spatially homogeneous system.
Varying  from 0.5 to 0.05, led to the four basic one-dimensional dynamics, namely stationary, smooth oscillatory, intermittent chaos and chaos covering (almost all) of the domain. In Fig.2(a), the approximations have evolved into the spatially homogeneous stationary states * u and * v . In Fig.2(c), we used the same initial data and parameter values as in [15], resulting in the same intermittent spatial structures at roughly the same positions. Again with same set of parameter as in Fig.2 with , the numerical solutions of phytoplankton dynamics are shown in Fig. 3(a)-(d).

Numerical Solutions in Two-Dimensions
In Figs. 7 4  , the numerical solutions of the phytoplankton dynamics (5)-(6) are plotted with space x and y coordinate. The initial data and parameter values of Figs. 7 4  , were chosen to correspond to similar as in [15].
In Figs. 7 4  , the initial data are small spatial perturbations of the stationary solutions * u and * v of the spatially homogeneous system. In Fig. 4, the evolution of the system led to the formation of spiral patterns, followed by irregular patches covering the whole domain (see Fig. 5), confirming the qualitative behavior reported in [15]. The size of these patches has been related to the characteristic length of observed plankton patterns in the ocean [15] and the references therein.
In Fig. 4, we observe that as t  is reduced, the initial spiral patterns disappear,indicating that the spiral pattern of [15]. The different numerical results may be due to the different domain shape used in [15], namely, a rectangle.
. Spatial patterns are obtained at different values of:

Dynamical Behavior and Stability Analysis
Now, we will study the global behavior of the system (1) and (2) , Further, it is clear from the above expression that Now, we have established the boundedness of the system, considering the function and using (1)-(2), we get If we choose a positive number Now using the comparison theorem, as

M t supW
Again, consider the set Hence, we can state the following lemma: Lemma 1: The system is uniformly bounded in the region . We have already established that the dynamical system (1)-(2) has three equilibria, namely, (i) . Here, the general variation matrix corresponding to the system is given by

Oscillatory Behavior
In this subsection, we will establish the oscillatory behavior of the system using Bendixson-Dulac Theorem. Now, we explore the possibility of a limit cycle around the endemic equilibrium with the help of the Bendixon-Dulec criterion.
Hence, there does not exist any limit cycle around the endemic equilibrium. Now, for the numerical simulation, to make the system (1) and (2)  , rT t = , and re-scaling the parameters via, . Hence our model reduces to .
Taking the same set of parametric values and initial data as in figures 2-3 for the system (12) and (13) Figure 9: Susceptible phytoplankton and infected phytoplankton population density, Time series plots for

Conclusion
In this paper, we have considered a phytoplankton dynamics with spatial movement. We have studied the reaction diffusion model in both one and two dimension space coordinates and also investigated their stability conditions of temporal model. By analyzing the model system (5)-(6) with constant diffusion. We also observed that the point of extinction of both susceptible and infected phytoplankton is unstable while the point of existence of susceptible phytoplankton and extinction of the infected phytoplankton is stable or unstable depending upon the parameter values and also observed that this system stability dependent to gamma as parameter value (i,e. fraction of infected phytoplankton population). Numerical results are also provided, using the standard Galerkin finiteelement method, leading to spiral waves. In one dimension, showed how modest changes in a single parameter of the system with Kinetics, namely  , can lead to dramatic changes in the qualitative dynamics of solutions. The non-dimensional parameter  is the related to half-saturation abundance of susceptible phytoplankton. Thus, the ecological implication is that, how fast the consumption rate of infected phytoplankton saturate with increasing susceptible phytoplankton density can have a profound effect on the fate of the system. In contrast, the experiments in [15], focused on the solution dynamics depend on the current state of the system (initial conditions). The dynamics of the spatially extended system are complicated and will depend on the system parameters, the initial data, and also the specifics of habitat geometry. Nevertheless, there are situations where the local dynamics of solutions gives us important clues to the behavior in the spatially extended situation. The one and two dimensions numerical experiments presented in this paper employ parameter sets that guarantee stable limit cycles in the reaction kinetics. As the diffusion coefficients used are equal, and the direction field of the reaction kinetics does not point out of the limit cycle (see, for example, Fig.1), the limit cycle encloses an 'invariant region' in phase space [20]. Consequently, as the initial data, used to generate Figs. 7 4  , lies entirely in the limit cycle at every point of the domain, the solution is trapped in this region for all time. Numerical experiments for different values T shown in Figs. 6 4  with same value of  and moreover, Fig.4 and Fig.6 shows the system pattern for different initial conditions. Further, it is also observed form Fig.3a-c, Fig.7a-d and Fig. 9a-c in 1-D space, 2-D space and temporal domains respectively that as  increases the system around endemic equilibrium become more stable as compare to other numerical experiments shown in other figures. Further, from Fig.3d, Fig.7e-f and Fig.9d as the level of half-saturation rate (  ) decreases even the stable system with respect to  (i,e. fraction of infected phytoplankton population) become unstable, i.e., system starts oscillating. Hence the the rate of growth of susceptible phytoplankton due to the dead infected phytoplankton (which become nutrients of susceptible phytoplankton after bacterial decomposition) is a major factor for the stability of the system.