Dynamical study of infochemical influences on tropic interaction of diffusive plankton system

We propose a model for tropic interaction among the infochemical-producing phytoplankton and non-info chemical-producing phytoplankton and microzooplankton. Volatile information-conveying chemicals (infochemicals) released by phytoplankton play an important role in the food webs of marine ecosystems. Microzooplankton is an ecologically important grazer of phytoplankton for coexistence of a large number of phytoplankton species. Here, we discuss how information transferred by dimethyl sulfide shapes the interaction of phytoplankton. Phytoplankton deterrents may lead to propagation of IPP bloom. The interaction between IPP and microzooplankton follows the Beddington–DeAngelis-type functional response. Analytically, we discuss boundedness, stability and Turing instability of the model system. We perform numerical simulation for temporal (ODE model) as well as a spatial model system. Our numerical investigation shows that microzooplankton grazing refuse of IPP leads to oscillatory dynamics. Increasing diffusion coefficient of microzooplankton shows Turing instability. Time evolution also plays an important role in the stability of system dynamics. The results obtained in this paper are useful to understand the dominance of algal bloom in coastal and estuarine ecosystem.


Introduction
Study of algal bloom is classical and widely recognized for the wealth of marine ecosystem. The research to investigate the mechanisms behind such phenomena is still important to understand the plankton dynamics and ultimately for other marine animals. In food web interaction infochemical, kairomones and pheromones used by an individual can exploit such chemical for searching prey or to avoid predator. Phytoplankton deterrents by producing infochemical may decrease the grazing rate on the producing species in marine environment [1]. Algae species are extremely flexible in their growth from production of toxic substance, biochemical composition and deterrent compounds [2]. Microzooplankton is an important grazer of bacteria and phytoplankton, and their ability to respond on infochemical cues confers a strong selective advantage [3]. The chemicals released by phytoplankton work as a defense mechanism and decrease the grazing pressure, but it also enhances the searching ability of zooplankter [4,5]. Resulting grazing-induced release of infochemical in multi-tropic interaction attracts carnivorous predators that reduce the grazing pressure on primary producer [6]. Recently, the models investigating the effect of toxin-producing phytoplankton (TPP) received much attention [7][8][9][10][11]. But grazing-induced toxic deterrents have received less attention. Inclusion of such toxic deterrents plays an important role for understanding and restructuring the dynamics of plankton bloom. Mathematical models describing the interaction between the toxic phytoplankton and other aquatic species have been discussed by many authors. Bairagi et al. [12] have proposed interacting nutrient-plankton dynamics and suggested that interactions among these species are very complex and situation specific. Model for interacting toxin-producing phytoplankton, non-toxic phytoplankton and zooplankton observed that the populations are independent of the magnitude of the steady-state component [13]. A two zooplankton and one phytoplankton model with toxicity observed that coefficient of toxin plays an important role for existence of Hopf bifurcation around the interior equilibrium [14]. Chatterjee and Pal [15] studied the effect of toxin in nutrient-plankton model and observed that toxic phytoplankton may change the steady-state behavior. De Silva and Jang [16] observed that the mutual interference of zooplankton diminishes harmful blooms.
Many researchers have used the predator-prey interaction models to study the spatiotemporal dynamics in plankton system [10,[17][18][19][20]. Spatiotemporal patterns in plankton dynamics with the sequence of chaotic spiral, strip and spot patterns were studied by Rao [21]. Wang et al. [22] observed that the direction of cross-diffusion influences the spatial distribution. Sekerci and Petrovskii [23] studied the oxygen-plankton model system in pattern formation and observed (for somewhat different initial conditions) irregular spatiotemporal patterns, simple spirals, double spirals or "mushroom-like" structures and "snake-like" structures. Thakur et al. [11] have proposed diffusive three species plankton models with toxic effect for wetland ecosystem. Chaudhuri et al. [24] studied toxic phytoplankton-induced spatiotemporal patterns observed that the populations become inhomogeneous in the presence of toxin-producing phytoplankton. Sekerci and Petrovskii [25] have proposed coupled plankton-oxygen model with effect of global warming. Sekerci [26] has studied theoretically by considering a coupled of the oxygen-plankton model where planktons habitat changes and beachhead as a response to climate change. Recently, some mathematical models exploring the dynamics of interaction between infochemical-producing phytoplankton and microzooplankton have been studied by [1,3,5,[27][28][29][30][31]. Study shows that the inclusion of grazing-induced infochemical term acts as attractor of predatory copepods which can stabilize the system [27]. Inclusion of such term may promote plankton bloom [33]. Hence, DMS acts as principal promoter of multi-trophic interaction [6]. With this motivation, our objective is to understand the effect of grazing-induced DMS and spatiotemporal evolution on system dynamics. The organization of the paper is as follows. Section 2 deals with model formulation, and stability analysis is discussed in Sect. 3. Section 4 describes numerical results. Discussion and conclusion are discussed in Sect. 5 and Sect. 6, respectively.

Model formulation
We present interacting plankton system where microzooplankton M(x; t) grazes on two competing phytoplanktons: P 1 (x;t) infochemical-producing phytoplankton (IPP) and P 2 (x;t) non-infochemical-producing phytoplankton (NIP) at any location x and time t. Both phytoplankton species grow logistically with their intrinsic growth rates r 1 , r 2 and carrying capacities K 1 and K 2 , respectively. The growth of phytoplankton is limited by availability of nutrients and light. Phytoplankton growth is also limited by an interspecific competition term, where 1 is the competitive effect of non-infochemical producers on infochemical producers and 2 is the vice versa. The terms are assumed to represent growth limiting factors such as nutrient depletion and/ or shading by the other species [34]. The relationship between (P 1 , P 2 , M) is defined by Holling type II or BDtype functional response. Microzooplankton predates non-infochemical-producing phytoplankton (NIP) at a rate of w 1 with the maximum conversion rate 1 and predates infochemical-producing phytoplankton (IPP) at a rate of w 2 with the reduction rate 2 , respectively, while k is the half-saturation constant. A key point to note is that although copepods are not explicitly modeled here, the effect of copepod predation is accounted for in the microzooplankton mortality terms in Eq. (2.1). Contrary to studies on toxic phytoplankton [8], DMS has no toxic effects on microzooplankton grazers [31,32], and therefore, IPP has no direct adverse effect on microzooplankton, but their mortality is indirectly increased by a factor proportional to grazing on IPP [1]. The mortality term m represents any mortality incurred by microzooplankton in the absence of infochemicals and the interspecific competition coefficients ( ) of microzooplankton that express the self-limitation of microzooplankton. With this assumption, we propose the model given by: non-vanishing initial conditions (2.1)  (Table 1).

Stability analysis of temporal model system
In this subsection, we have discussed the nonnegative equilibrium points and their stability properties of the model system in the absence of diffusion. The reduced system of ordinary differential equation is as follows: The proof of this lemma can be obtained by simple calculation and hence omitted.
The non-spatial model system (3.1) possesses seven nonnegative real equilibrium points: (v) The planar equilibrium point E 4 = (P 1 , 0,M) exists on the P 1 M-plane. It may be noted that for M to be positive, we must have (3.6)  It may be noted that for M to be positive, we must have This shows that E 5 = (0,P 2 ,M) exists under the condition (3.13) and (3.14). (vii) Existence of interior equilibrium point E 6 = (P * 1 , P * 2 , M * ) has been established below after following [35]. In this case, P * 1 , P * 2 and M * are the positive solutions of following equations: From Eq. (3.15), we get (3.8) (3.12) (3.14) From Eq.(3.20), when P 2 = 0 , then G(P 1 , 0) = 0 has a root P 1b , which is the solution of the following equation Also, we have We noted that holds. We noted that the two isoclines given by Eqs. (3.19) and (3.20) intersect at a unique point (P * 1 , P * 2 ) if in addition to the conditions given in Eqs. (3.22), (3.23), (3.25) and (3.28), the inequality P 1b < P 1a holds. (3.24) (3.28) The local stability of each equilibrium point is now discussed by deriving the variance matrices and using the Routh-Hurwitz criterion. The findings were obtained below: (k +P 2 ) 2 is negative or positive, respectively, if K 2 w 2M < r 2 (K +P 2 ) 2 . (vii) The variational matrix along E 6 (P * 1 , P * 2 , M * ) is given by The characteristic equation for the above matrix W is given by where Theorem 1 Assume that the E 6 (P * 1 , P * 2 , M * ) is positive equilibrium point of the system (3.1). Therefore, the equilibrium point E 6 (P * 1 , P * 2 , M * ) is locally asymptotically stable when The proof of the theorem follows from the Routh-Hurwitz criterion, hence omitted. . (3.32) where Sufficient condition for dV dt to be negative is that the following inequalities hold: , we note that , , (3.34) m 11 > 0,

Stability analysis of spatial model system
In this section, we discuss the stability of interior equilibrium of the diffusive model system. In order to derive the condition of stability, we linearized the model system (2.1) about E 6 (P * 1 , P * 2 , M * ) with small perturbation X(x, t), Y(x, t) and Z(x, t) as P 1 = P * 1 + X (x, t), P 2 = P * 2 + Y(x, t) and M = M * + Z(x, t) . The linearized form of model system is obtained as: Let us assume the solution of system (3. From Eq. (3.42) and using the Routh-Hurwitz criterion, the above theorem follows immediately.

Theorem 4
If the positive equilibrium point E 6 of the model system (3.1) is globally asymptotically stable, then corresponding uniform steady state of model system (2.1) remains globally asymptotically stable.
Proof For stability behavior of the system (2.1), we define a positive definite function V 1 (t) , given by where V (P 1 , P 2 , M) is given in Eq. (3.32

Turing instability
In this subsection, we have derived the required conditions for the existence of Turing instability of the spatial phytoplankton-microzooplankton system (2.1). Due to spatial diffusion, the occurrence of Turing instability changes the stable equilibrium to the unstable one. Mathematically, Turing instability requires at least one of the roots of the characteristic (3.42) which has a nonnegative or nonnegative real part or on the other hand, Re( ) > 0 for some i > 0 . Irrespective of the sign of 1 and 2 , the diffusion-induced instability occurs when 3 = P( i ) < 0. Hence, the condition for diffusive instability is given by a 22 a 33 − a 23 a 32 ) + D 2 (a 11 a 33 − a 13 a 31 ) +D 3 (a 11 a 22 − a 12 a 21 ). P is cubic polynomial in i . The critical value of P( i ) occurs at the i = icr , where (3.46) (3.47) For positive value of critical points i = icr , we require: Now, we have taken the set of parameter values from Lewis et al. [1]. The parameter values given here are fixed unless otherwise stated. These values were taken from parameter ranges found in the literature to give biologically reasonable results.
We have chosen the diffusion coefficients on the basis of the values reported in [18,[36][37][38]. Instabilities of the spatially uniform distribution can obtain if the species move with different velocities regardless of which one is faster [37]. If we take the diffusion coefficient of microzooplankton D 3 which is greater than the diffusion coefficients of infochemical-producing phytoplankton D 1 and non-infochemical-producing phytoplankton D 2 , then Turing instability will prevail in ecological sense (c.f. Fig. 1).
At this set of parameter values as given in Eq. (3.50), the equilibrium point is given by E 6 (P * 1 , P * 2 , M * ) = (25.6430 , 437.9112, 51.8768). The equilibrium point E 6 (P * 1 , P * 2 , M * ) is stable without diffusion. In the presence of D 1 = D 2 = 0.0001, D 3 = 5, 3 = P( i ) changes its sign with i as shown in Fig. 1. At above set of parameter values, the critical value is icr = 13.4515 and corresponding value of P( icr ) = −0.0278 (c.f., Fig. 1). In this case, Turing instability is possible.

Results
In this section, we have performed numerical simulation of the spatial model system (2.1) for one and two dimensions. We have also described the bifurcation diagram and plot of time vs population density for model system (3.1). To investigate the spatiotemporal dynamics of the model system (2.1), we have taken the set of parameter values as given in Eq. (3.50).    and ( ) as bifurcation parameters, respectively. The successive variation of NIP and microzooplankton has been considered in the ranges 0 ≤ P 2 ≤ 600 and 0 ≤ M ≤ 150, respectively, as a function of (w 2 ) which is in the range 0.5 ≤ w 2 ≤ 1.3 . All other parameter values are same as in Eq. (3.50). Under the bifurcation analysis of model system (3.1), we observe that the higher value of w 2 is responsible for the limit cycle oscillations and at the lower value of w 2 , the dynamics shows the stable behavior, and similarly, we observe that the higher value of ( ) is responsible for the stable dynamics and at the lower value of ( ) , the dynamics shows the limit cycle oscillations. From Figs. 2(a)-2(b), it can be seen that initially, the model system shows the stable behavior but when we take the value (w 2 ) > 0.849, the system dynamics becomes oscillatory behavior. Thus, it can say that w 2c = 0.849 is the critical value (w 2 ) for occurrence of bifurcation. Similarly, in Fig.3 we have observed the effect of ( ) on the system dynamics through bifurcation diagram. System dynamics of model (3.1) becomes stable to oscillatory behavior after the critical value = c = 0.00019 . We have observed the effect of (w 2 ) on the dynamical behavior of the model system (3.1) through Fig. 4. In Fig. 4, the plot (time versus population densities) has been obtained. We observe that the increasing value of the (w 2 ) system shows stable to oscillatory behavior of model system (c.f., Fig. 4(a)-4(b)). In Figs. 5 and 6, the spatiotemporal patterns reflect the effect of different parameter values at both spatial and temporal scales. Figure 5 depicts the effect of w 2 on spatiotemporal

Discussion
The microzooplankton grazing selectivity is investigated through a simple mathematical competition model under the consideration of infochemical. The infochemical produced by phytoplankton has directly influenced the feeding as well as selectivity nature of zooplankton. Additionally, grazinginduced DMS infochemicals for multi-trophic interaction are considered for microzooplankton susceptibility. The competition factor between the phytoplanktons is taken that limited their growth. Here, we have taken one more important factor, i.e., intraspecific competition coefficient of microzooplankton , which expresses the self-limitation of microzooplankton in the form M 2 . All the needful comparison with past studies of the study is summarized below: (i) A three-species plankton model with competition is studied by Lewis et al. [1], but the effect of spatial variability on the population dynamics has been avoided. Assuming spatial variability in the form of diffusion fascinated our study. (ii) Thakur et al. [11] investigated a diffusion-based competing plankton system for Sundarban mangrove system and mainly focused on the role of reduction rate of zooplankton, but infochemical term in phytoplankton has been ignored in the model formulation. Therefore, this term has been taken in the third equation of the system (2.1). (iii) The Turing Instability in the case of without diffusion is well obtained for the appropriate set of parameters (c.f., Fig. 1). (iv) Without diffusion, the role of w 2 and is explored.
The increasing value of w 2 makes a stable system to unstable, whereas the increasing value of makes an unstable system to stable one. (v) Spatiotemporal patterns are plotted for different values of w 2 and δ. A stable to oscillatory temporal pattern is observed for different values of w 2 , whereas a stable temporal pattern is observed for different value of δ. A similar kind of temporal pattern is also found by Thakur et al. [10]. (vi) Spatial distribution for different diffusion coefficients of D 3 and time intervals are plotted. An irregular patchy pattern is observed for microzooplankton species.

Conclusions
In this paper, we have investigated a nonlinear diffusive three species food chain models interacting IPP, NIP and microzooplankton with the effect of grazing-induced DMS infochemicals in multi-tropic interaction. We have implicitly incorporated the parameter to understand the sensitivity of predators to infochemicals signals that enhance grazing rates of microzooplankton by copepod. Production of DMS releases the grazing pressure on IPP and allows to reach its densities equivalent to plankton bloom. We first proved the boundedness of the model system which shows that the model system is biologically well behaved. Analytically, we have discussed the stability of non-trivial equilibrium and existence of Turing instability and observed that higher value of diffusion coefficient of microzooplankton leads to occurrence of Turing instability (c.f., Fig. 1). With the help of numerical simulation, we tried to understand the extent to which infochemical influences the tropic interaction and spatiotemporal patterns in heterogeneous marine environment. We have taken weak competition between planktons (i.e., 1 < 2 ) and infochemicals as grazing attractors (i.e., w 1 > w 2 ). Here, we have taken IPP which has lower intrinsic growth rate than the NIP (i.e., r 1 < r 2 ) that corresponds to cost of production of infochemicals. Firstly, we have carried out the bifurcation analysis to understand the microzooplankton grazing selectivity and its intraspecific interference (c.f., Figs. 2 and 3). Figure. 2 shows that increasing value of w 2 leads to oscillatory dynamics and system may become stable with increasing intraspecific competition ( ). Further, we have studied the spatiotemporal pattern to understand the effect of w 2 and . Increasing the value of w 2 and observed stable to oscillatory system dynamics (c.f., Fig. 5) and increasing and observed stable behavior (c.f., Fig. 6). Finally, patch distribution of microzooplankton is studied with the help of snapshots and observed that increasing diffusion coefficients and time evolution show lowering the patch densities (c.f., Figs. 7 and 8). Our result shows that increasing predation of NIP or grazing refuse of IPP (due to grazing-induced infochemicals cue that attracts predatory copepods) destabilizes the system dynamics. Our study may be useful for better understanding of plankton bloom in the ecology of marine food webs.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.