Effects of zooplankton selectivity on phytoplankton in an ecosystem affected by free-viruses and environmental toxins

1.   Introduction

After the pioneering work of Anderson and May ^[1], research on predator-prey-parasite system has been explored extensively ^[2,3,4]. It is now well documented that predator-prey-virus models exhibit more complex and interesting dynamical behaviors than the conventional ecological or epidemiological models ^[5]. Parasites play a prominent role in trophic interactions and strongly influence the food web characteristics. It can also manipulate the potential of its host and affect the abundance of a predator ^[6]. A number of studies have demonstrated that the predators tend to forage for prey with high parasite burden ^[7]. The different factors regulating the intricate planktonic food web are, namely, selective predation ^[8], nutrient inflow ^[9], grazing by higher predators ^[10], etc. Zooplankton are capable of discriminatory feeding ^[11], and its food selectivity criteria vary to a large extent ^[12]. Zooplankton can selectively ingest its food depending on a myriad of factors, namely size ^[13], digestibility ^[8], toxicity ^[14], availability ^[15], and nutrition value ^[16] of food. A number of studies have documented that Calanoid copepods can discriminate between toxic and nontoxic dinoflagellates ^[17], noxious and innoxious blue-greens ^[18], and live and dead algae of the same species ^[19].

Marine viruses have been recognized to have a decisive lead in altering the physiological and biochemical compositions of their hosts ^[20,21]. During viral infection, the host-cell characteristics such as cell size, the diffused infochemicals and the cell lipid membrane properties are changed, which influence the feeding behavior as well as the growth rate of zooplankton ^[22,23]. Emiliania huxleyi are extremely abundant marine phytoplankton species of the cocolithophore group which are frequently found in most of the marine ecosystems. Under favourable environmental conditions, E. huxleyi outburst into large seasonal blooms ^[24]. Gaint physcodnaviruses (E. huxleyi viruses, EhVs) are known to infect and lyse E. huxleyi cells and to be highly associated with the regulation of populations and the termination of blooms ^[25]. It has been documented that the exposure of E. huxleyi to EhVs causes more production of DMS, DMSP compared to normal bloom conditions ^[26,27]. Some studies have proposed that the infochemicals released from the infected cells act as an activated chemical signal which triggers the selective feeding behavior of the grazer zooplankton ^{[23,28,29,27]}. The role of viral infection in the producer trophic level of marine ecosystems has been investigated by several authors ^{[30,31,32,33]}. Some model based studies explicitly considered the role of free-viruses in eco-epidemiological system ^[6,34], assuming that both the healthy and infected phytoplankton are equally preferable for predation.

Pollutants have a potentially hazardous effect on marine ecosystem. Pollution may arise from different sources, the most well known are influx of domestic sewage, industrial waste, agricultural run-off etc. The environmental toxins affect the plankton communities at different levels $-$ abundance, growth strategies, dominance and succession patterns. During predation these toxins are passed onto higher trophic level in the food chain in a cascading manner. The contact rate between environmental toxins and marine organisms is mainly dependent on their distribution pattern and behavioral attributes in the ecosystems ^[35]. Environmental toxins cause an implicit suppression in the growth of phytoplankton cells among a wide range of phytoplankton population ^[36,37]. In some recent studies, this growth suppression behavior of environmental toxins have been incorporated mathematically ^[38,39]. It has been shown that environmental toxin-induced growth suppression of phytoplankton population can lead to a destabilized system. However, the depletion of environmental toxins from the aquatic system plays a crucial role to stabilize the system. Yu et al. ^[40] have studied a stochastic phytoplankton-zooplankton model by considering toxin-producing phytoplankton and Markov switching in an impulsive polluted environment; they found that environmental fluctuations and exogenous toxicant input have a great influence on the survival fate of plankton.

Recently, Bairagi and Adak ^[41] have studied the predator-prey dynamics, where the prey species are susceptible to infection and predators feed upon both the susceptible and infected prey with some preference. Their findings depict that understanding the consequences of predator's selectivity or preference may play a pivotal role to shape the community structure and ecological properties involving parasites. By dividing the entire phytoplankton population of a pelagic ecosystem into preferred and nonpreferred phytoplankton, Bairagi et al. ^[42] have shown that the selective predation of zooplankton and nutritional value of phytoplankton strongly impact on the complex plankton dynamics, including the bloom phenomenon. In the present investigation, we have captured the idea of selectivity tendency of grazer zooplankton on susceptible phytoplankton cells in the presence of infected ones, for the very first time. The objective of our study is to investigate the role of selectivity behavior on the persistence and dominance of healthy phytoplankton cells. We would also like to investigate the combined effects of free-viruses and environmental toxins on the dynamics of phytoplankton-zooplankton in the system, and how the selectivity behavior of zooplankton regulate the dynamics. The above ecological questions are answered by the model analysis and numerical simulations.

2.   The mathematical model

Pollution of freshwater and marine systems by anthropogenic sources has become a concern over the last several decades. Organic (e.g. triazine herbicides) ^[43] and inorganic compounds (e.g. heavy metals) ^[44] both may have harmful effects to the organisms. For example, photosynthesis of the green algae Selenastrum capricornutum is inhibited due to pollutants originating from industrial and recreational sources ^[45]. In a marine planktonic community comprising mainly diatoms and herbivorous copepods in a low silicate, elevated copper environment ^[46], copper harms diatoms. At low concentrations, the herbicide triazine also affects primary producers directly by inhibiting photosynthesis, while effects on subsequent trophic levels only would be indirect ^[46]. Contamination of water bodies with insecticides such as carbaryl, azadirachtin, or cypermetrin increases the death rate of zooplankton.

In phytoplankton population, viruses are not transmitted from infected phytoplankton to susceptible phytoplankton by contact. Usually, virus particles attach to the host cell and release its genetic material into the cell. The virus then uses the host cells machinery to replicate its own genetic material. Once replication has been completed, the virus particles burst out of the host cell into the extracellular space resulting in the death of the host cell. After escaped from the host cell, virus is ready to infect other susceptible phytoplankton. This kind of transmission converts directly a part of phytoplankton population to an infected class ^[6]. Moreover, in the absence of viral infection, primary consumers like copepods, access over E. huxleyi and other phytoplankton species at random but in the presence of viral infection zooplankton exhibit some sort of grazing preferences ^[23,27].

Let $S$, $I$ and $Z$ be the concentrations of susceptible phytoplankton, infected phytoplankton and zooplankton in an aquatic ecosystem, and let $V$ and $E$ be the concentrations of free-viruses and environmental toxins in the same system. We have the following assumptions:

1. In the absence of free-viruses $V$, environmental toxins $E$ and the grazer zooplankton $Z$, the susceptible phytoplankton $S$ grow logistically as $\frac{dS}{dt} = rS\left(1-\frac{S}{K}\right)$ with $r$ being the intrinsic growth rate and $K$ being the carrying capacity. The environmental toxins $E$ affects the growth rate of $S$, and it reduces to $\frac{r}{1+\gamma\gamma_1SE},$ where $\gamma$ is the contact rate of environmental toxins with phytoplankton and zooplankton, and $\gamma_1$ is the inhibition of the susceptible phytoplankton growth rate due to environmental toxins.

2. The susceptible phytoplankton $S$, becomes infected by direct contact with free viruses, $V$, living in the system following a disc function $\frac{\beta SV}{K_1+V}$ with the transmission parameter $\beta$ and the half saturation constant $K_1$. In the presence of free-viruses, $V$, the susceptible phytoplankton $S$ is consumed by the zooplankton $Z$ with the rate of $\frac{\alpha S}{d+S+\theta I}.$ The parameter $\alpha$ is the predation rate of the zooplankton $Z$; $d$ is the half saturation constant; and $\theta$ measures the selectivity/preference of the zooplankton to infected phytoplankton over the susceptible one or vice-versa. If $\theta > 1$, then zooplankton prefer infected phytoplankton, $\theta = 1$ implies equal preference to susceptible and infected phytoplankton, and $0 < \theta < 1$ means that zooplankton prefer susceptible phytoplankton.

3. It is assumed that the infected phytoplankton, $I$, do not reproduce ^[47,48], and they are removed by cell lysis before having the capability of reproducing ^[49]. The infected phytoplankton do not compete for resource with susceptible phytoplankton, $S$. The infected phytoplankton, $I$, is consumed by the zooplankton, $Z$, with the rate of $\frac{\alpha \theta I}{d+S+\theta I}.$ The infected phytoplankton is assumed to die at a constant rate $\mu$ with an additional mortality rate $\gamma\gamma_2E$ due to the environmental toxins, $E$, in the aquatic system.

4. The zooplankton population increase on consumption of susceptible and infected phytoplankton; experience natural death at a constant rate $\nu$, and additional death due to presence of environmental toxins in the system, $\gamma\gamma_3ZE$.

5. Let $b$ be the virus replication factor through the infected phytoplankton $I$, i.e., lysis of infected phytoplankton, on the average, producing $b$ virus particles ($b\gg1$). The decay rate of virus is assumed to be constant, $\delta$. The virus is removed through the infection of susceptible phytoplankton at the rate $\frac{\beta SV}{K_1+V}.$

6. The environmental toxins $E$ has a constant input rate $A$ and a depletion rate of $\rho$. The additional removal rates are due to being absorbed by both phytoplankton and zooplankton through their contact at a rate $\gamma$ ^[38,40].

The schematic diagram for the interplay among all the considered dynamical variables is depicted in Figure 1. Based on these ecological assumptions, we have the following system of nonlinear differential equations,

The parameters involved in system (2.1) are positive constants. System (2.1) is to be analyzed with the following initial conditions,

We assume $0 < \lambda_2\leq \lambda_1 < 1$, i.e., the per capita contribution of infected phytoplankton to the growth of zooplankton never exceeds that of the susceptible phytoplankton. Biological meanings of parameters in the system (2.1) is given in Table 1.

3.   Mathematical analysis

3.1.   Existence, positive invariance, boundedness and persistence of solutions of system (2.1)

In theoretical ecology, boundedness of a system implies that the system is well behaved. Boundedness of the solutions entails that none of the interacting populations grow exponentially for a long-time interval. The abundance of each population is bounded due to limited resource. Moreover, permanence of a system means the long-term survival of all populations of the system, no matter what the initial populations are. From mathematical point of view, permanence of a system means that strictly positive solutions do not have omega limit points on the boundary of the non-negative cone.

We have the following theorem.

Theorem 3.1. Every solution of system (2.1) with initial conditions (2.2) exists and unique in some interval $[0, \kappa)$, where $S(t) > 0$, $I(t) > 0$, $Z(t) > 0$, $V(t) > 0$, $E(t) > 0$, $\forall t\geq0$. Further, all non-negative solutions of system (2.1) that start in $\mathbb{R}^5_+$ are uniformly bounded, and the feasible region for system (2.1) is given by the following set

Further, if the following inequalities are satisfied:

then system (2.1) is permanent. The quantities $M$, $E_m$, $S_a$, $I_m$, $Z_m$ and $V_m$ are defined in the proof.

For proof of this theorem, see Appendix A.

From the conditions in (3.1), we note that the susceptible phytoplankton may disappear from the aquatic system for very high values of force of infection or its contact rate with environmental toxins; the infected phytoplankton and hence the free-viruses may extinct for very low values of the force of infection or for very high values of its contact rate with environmental toxins; the zooplankton may extinct for very high values of its natural mortality or death rate due to environmental toxins.

3.2.   The ecosystem in the absence of free-viruses and environmental toxins

In the absence of free-viruses and environmental toxins, system (2.1) reduces to,

whose dynamics have been well studied ^[50]. Here, we summarize its dynamics as follows. System (3.2) has the three feasible equilibria.

1. The plankton-free equilibrium point $\widetilde{e}_0 = (0, 0)$ is always a saddle.

2. If $K(\lambda_1\alpha-\nu) < d\nu$, then global stability at the zooplankton-free equilibrium $\widetilde{e}_1 = (K, 0)$.

3. If $K(\lambda_1\alpha-\nu) > d\nu$ and $K < \frac{d(\lambda_1\alpha+\nu)}{\lambda_1\alpha-\nu},$ then global stability at the coexistence equilibrium $\widetilde{e}_* = (\widetilde{S}_*, \widetilde{Z}_*)$, where

4. If $K(\lambda_1\alpha-\nu) > d\nu$ and $K > \frac{d(\lambda_1\alpha+\nu)}{\lambda_1\alpha-\nu},$ then globally stable and unique limit cycle around the coexistence equilibrium $\widetilde{e}_* = (\widetilde{S}_*, \widetilde{Z}_*)$.

3.3.   The ecosystem in the absence of free-viruses

In the absence of free-viruses in the aquatic system, model (2.1) reduces to,

Our model (3.3) differs from the models of ^[38,39] in the sense that here zooplankton are assumed to be affected by environmental toxins. In ^[38], the interaction between phytoplankton and zooplankton is assumed to be modified Beddington type. The effect of toxin liberation by phytoplankton on zooplankton and predation of zooplankton by fish are considered in ^[39]. Yu et al. ^[40] studied a stochastic model by considering that environmental toxins increase the death rate of phytoplankton as well as zooplankton.

System (3.3) has the following three non-negative equilibria:

1. The plankton-free equilibrium $\widehat{e}_0 = (0, 0, A/\rho)$, always feasible.

2. The zooplankton-free equilibrium $\widehat{e}_1 = (K, 0, A/(\gamma K+\rho))$, always feasible.

3. The coexistence equilibrium $\widehat{e}_* = (\widehat{S}_*, \widehat{Z}_*, \widehat{E}_*)$, where

and $\widehat{S}_*$ is given by positive root of the following equation:

Regarding local stability of equilibria of the system (3.3), we have the following theorem.

Theorem 3.2. 1. The equilibrium $\widehat{e}_0$ is always unstable.

2. The equilibrium $\widehat{e}_1$ is always stable.

3. The equilibrium $\widehat{e}_*$, if feasible, is locally asymptotically stable if and only if the following conditions hold,

For proof of this theorem, see Appendix B.

3.4.   The ecosystem in the absence of environmental toxins

In the absence of environmental toxins in the aquatic system, model (2.1) reduces to,

In ^[6], authors assumed that both the healthy and infected phytoplankton are equally preferable for predation. They showed that virus replication factor and half saturation constant of viral infection play an important role in persistence of infection in the system. They also observed that for low infection rate, the virus cannot persist in the system and the system becomes infection-free. If the force of infection increases then above a threshold value, all the species coexist and the system becomes endemic. However, for very large values of force of infection, the zooplankton population extinct from the system. The preference behavior of predator over susceptible prey is considered in ^[41]; the explicit role of parasites is not taken into account. In ^[42], the authors showed that the phytoplankton-zooplankton interaction is very complex and the plankton dynamics strongly depends on the selective predation of zooplankton.

System (3.6) exhibits five non-negative equilibria given as,

1. The population-free equilibrium $\overline{e}_0 = (0, 0, 0, 0)$, which is always feasible.

2. The equilibrium with only susceptible phytoplankton $\overline{e}_1 = (K, 0, 0, 0)$, which is always feasible.

3. The disease-free equilibrium $\overline{e}_2 = (\overline{S}_2, 0, \overline{Z}_2, 0)$, where $\overline{S}_2 = \widetilde{S}_*$ and $\overline{Z}_2 = \widetilde{Z}_*$. The equilibrium $\overline{e}_2$ is feasible provided $K(\lambda_1\alpha-\nu) > d\nu$.

4. The zooplankton-free equilibrium $\overline{e}_3 = (\overline{S}_3, \overline{I}_3, 0, \overline{V}_3)$, where $\overline{S}_3 = \frac{K}{r}\left(r-\frac{\beta \overline{V}_3}{K_1+\overline{V}_3}\right),$ $\overline{I}_3 = \frac{K\beta\overline{V}_3}{\mu r(K_1+\overline{V}_3)}\left(1-\frac{\beta \overline{V}_3}{K_1+\overline{V}_3}\right)$ and $\overline{V}_3$ is positive root of the quadratic,

where

As $\overline{a}_2 > 0$, equation (3.7) has exactly one positive root if $\overline{a}_0 < 0$. A necessary condition for $\overline{S}_3$ to be positive is $rK_1+(r-\beta)\overline{V}_3 > 0$. Thus, the equilibrium $\overline{e}_3$ is feasible provided,

In case, the latter is not satisfied, then equation (3.7) has either two or no positive roots.

5. The coexistence equilibrium $\overline{e}_* = (\overline{S}_*, \overline{I}_*, \overline{Z}_*, \overline{V}_*)$. In equilibrium $\overline{e}_*$, we have

and $\overline{V}_*$ is positive root of the equation

Regarding local stability of equilibria of the system (3.6), we have the following theorem.

Theorem 3.3. 1. The equilibrium $\overline{e}_0$ is always unstable.

2. The equilibrium $\overline{e}_1$ is locally asymptotically stable if

3. The equilibrium $\overline{e}_2$, if feasible, is locally asymptotically stable if and only if the following conditions are satisfied,

4. The equilibrium $\overline{e}_3$, if feasible, is locally asymptotically stable provided,

5. The equilibrium $\overline{e}_*$, if feasible, is locally asymptotically stable if and only if the following conditions are satisfied,

For proof of this theorem, see Appendix C.

3.5.   Equilibrium states of the full system (2.1)

System (2.1) exhibits the following five non-negative equilibria:

1. The equilibrium $e_0 = (0, 0, 0, 0, A/\rho)$, which is always feasible.

2. The equilibrium $e_1 = (K, 0, 0, 0, A/(\gamma K+\rho))$, which is always feasible.

3. The equilibrium $e_2 = (S_2, 0, Z_2, 0, E_2)$, where $S_2 = \widehat{S}_*$, $Z_2 = \widehat{Z}_*$ and $E_2 = \widehat{E}_*$.

4. The equilibrium $e_3 = (S_3, I_3, 0, V_3, E_3)$, where

$S_3$ and $V_3$ are positive solutions of the following isoclines:

5. The equilibrium $e_* = (S_*, I_*, Z_*, V_*, E_*)$, where

$S_*$ and $V_*$ are positive solutions of the following isoclines:

Regarding local stability of equilibria of the system (2.1), we have the following theorem.

Theorem 3.4. 1. The equilibrium $e_0$ is always unstable.

2. The equilibrium $e_1$ is locally asymptotically stable provided,

3. The equilibrium $e_2$, if feasible, is locally asymptotically stable if and only if the following conditions hold,

4. The equilibrium $e_3$, if feasible, is locally asymptotically stable if and only if the following conditions hold,

5. The equilibrium $e_*$, if feasible, is locally asymptotically stable if and only if the following conditions hold,

For proof of this theorem, see Appendix D.

3.6.   Existence of Hopf bifurcation

In the context of a biological system, bifurcations describe how the system dynamics may qualitatively drastically change if a parameter varies. A Hopf bifurcation is a critical point where a system's stability switches and a periodic solution arises. We consider here some parameters: selectivity/preference of zooplankton ($\theta$), force of infection ($\beta$), growth suppression of susceptible phytoplankton due to environmental toxins ($\gamma_1$) and carrying capacity of susceptible phytoplankton ($K$) as bifurcation parameters. Analytically we investigate for the possibility of Hopf bifurcation from the coexistence equilibrium $e_*$ by taking the selectivity/preference of zooplankton ($\theta$) as bifurcation parameter, keeping other parameters fixed.

Regarding Hopf bifurcation of the system (2.1) around equilibrium $e_*$ with respect to the parameter $\theta$, we have the following theorem.

Theorem 3.5. When the intensity of selectivity ($\theta$) exceeds a critical value, the system (2.1) enters into Hopf-bifurcation around the coexistence equilibrium $e_*$. The necessary and sufficient conditions for the occurrence of Hopf bifurcation are that there exists $\theta = \theta^*$ such that

1. $\psi(\theta^*)\equiv \{C_3(\theta^*)-C_1(\theta^*)C_2(\theta^*)\}\{C_5(\theta^*)C_2(\theta^*)-C_3(\theta^*)C_4(\theta^*)\}-\{C_5(\theta^*)-C_1(\theta^*)C_4(\theta^*)\}^2 = 0$.

2. $C_1(\theta^*) > 0, \ C_1(\theta^*)C_2(\theta^*)-C_3(\theta^*) > 0, \ C_3(\theta^*)-C_1(\theta^*)\omega^*_0 > 0$, $\omega^*_0 = \frac{C_5(\theta^*)-C_1(\theta^*)C_4(\theta^*)}{C_3(\theta^*)-C_1(\theta^*)C_2(\theta^*)} > 0.$

3. $\frac{d\psi(\theta)}{d\theta}\Bigg\vert_{\theta = \theta^*} = 0.$

For proof of this theorem, see Appendix E.

In Table 2, we listed the equilibria of the systems (2.1), (3.2), (3.3) and (3.6), and the conditions for their stability.

4.   Numerical simulations

Here, we report the simulations to investigate the behavior of system (2.1), performed using the Matlab variable step Runge-Kutta solver ode45. The set of parameter values are chosen within the range prescribed in various previous literature sources ^{[6,34,38,39,41,42]}, and are given in Table 1. Unless it is mentioned, the values of parameters used for numerical simulations are the same as in Table 1. In Table 3, we listed the equilibria of the systems (2.1), (3.2), (3.3) and (3.6), and the corresponding eigenvalues of the system's Jacobian. In Table 3, we choose parameter values in such a way that coexistence equilibrium of each system is stable. For system (2.1), we choose parameter values as in Table 1, while for other systems, the changed parameters are: $\alpha = 0.03$ in system (3.2); $\alpha = 0.057$ in system (3.3); $\alpha = 0.057$ and $\theta = 0.15$ in system (3.6). From the table, we note that the system (3.3) has two coexistence equilibria: $\widehat{e}_{*1}$ and $\widehat{e}_{*2}$, of which the former is stable and the latter is unstable. Also, the equilibria $\widehat{e}_1$ and $\widehat{e}_{*1}$ are simultaneously stable, thus the system (3.3) shows bistability behavior between these two equilibria. Bistability is a phenomenon where the system can converge to two different equilibria in the same parametric region based on the variation of the initial conditions.

4.1.   Dynamics of the subsystems (3.2), (3.3) and (3.6)

First, we investigate the dynamics of the system (2.1) in the absence of free viruses and environmental toxins i.e., subsystem (3.2), Figure 2. We observe that for system (3.2), the coexistence equilibrium $\widetilde{e}_*$ is stable at $\alpha = 0.03$ but the stability is lost on increasing the values of $\alpha$ or $K$. Moreover, on decreasing the values of $\lambda_1$, the coexistence equilibrium $\widetilde{e}_*$ losses its stability and the zooplankton-free equilibrium $\widetilde{e}_1$ appears in the system. We observe that the system (3.2) is oscillatory at $\alpha = 0.057$, while the introduction of environmental toxins in system (i.e., system (3.3)) makes it stable (Figure 3(a)), but for higher values of input rate of environmental toxins ($A = 20.5$), the system settles to zooplankton-free steady state (Figure 3(b)). Next, we see bistability behavior of the system (3.3) in Figure 4(a). Here, blue region stands for stability domain for the equilibrium $\widehat{e}_1$, and red region stands for stability domain for the equilibria $\widehat{e}_1$ and $\widehat{e}_{*1}$. We observe that the system (3.3) shows bistability between the zooplankton-free equilibrium $\widehat{e}_1$ and the coexistence equilibrium $\widehat{e}_{*1}$ for large input rate of environmental toxins in the aquatic system. Thus, high input rate of environmental toxins leads to extinction of zooplankton or persistence of all the populations depending on the initial conditions. We also plot basin of attractions to confirm bistability behavior of the system. In Figure 4(b), magenta region represents the basin of attraction for the equilibrium $\widehat{e}_1$, and the green region represents the basin of attraction for the equilibrium $\widehat{e}_{*1}$. Next, we observe the dynamics of system (3.6), i.e., the system with free-viruses but without environmental toxins. We observe that the introduction of free-viruses drives the unstable equilibrium of $S$$-$$Z$ system to stable equilibrium of $S$$-$$I$$-$$Z$$-$$V$ system, Figure 5. On increasing the preference of zooplankton on infected phytoplankton, system (3.6) shows limit cycle oscillation around the disease-free equilibrium, $\overline{e}_2$ (see Figure 6). Next, we observe the behavior of system (3.6) by varying the force of infection, $\beta$ (see Figure 7). We find that for low values of force of infection, the system is oscillatory around the disease-free equilibrium, but the oscillation is around the coexistence equilibrium on increasing the values of force of infection; further increase in the values of force of infection results in stable coexistence of all the populations.

4.2.   Dynamics of system (2.1)

4.2.1.   Sensitivity analysis

We employed a global sensitivity analysis to assess the impact of uncertainty and the sensitivity of the outcomes of the numerical simulations to variations in each parameter of the model (2.1) using Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCC) ^[51,52]. To generate the LHS matrices, we assume that all the model parameters are uniformly distributed. We performed PRCC analyses for three different time points in order to investigate which parameters consistently influence free-viruses and environmental toxins in time. Nonlinear and monotone relationships are observed for the free-viruses and environmental toxins with the input parameters of the model (2.1), which is a prerequisite for computing PRCCs. Then, simulations of the model per LHS run were carried out, using the baseline values tabulated in Table 4 and the ranges as 25% from the baseline values (in either direction). The indexes are evaluated at the time points 200,400 and 600 days. For each of the three time points, denoted by 1, 2 and 3, PRCC results are shown in Figs. 8(a) & 8(b) for free-viruses and environmental toxins, respectively. We observe that the PRCC values of free-viruses decrease with time for the parameters, $\alpha$, $d$, $\nu$, $\gamma_1$ and $\gamma_3$ while for the environmental toxins, the PRCC values decrease for the parameters $\beta$ and $\theta$, and increases for the parameter $A$, indicating that these parameters consistently influence free-viruses/environmental toxins in the system; for other parameters there is no trend of either increasing or decreasing in the PRCC's values.

4.2.2.   Existence of Hopf-bifurcation

We discuss the Hopf bifurcation near the coexistence equilibrium $e_*$ on varying the intensity of selectivity ($\theta$), inhibition of phytoplankton growth due to environmental toxins ($\gamma_1$), force of infection ($\beta$), and carrying capacity of susceptible phytoplankton ($K$). For the set of parameter values in Table 1, system (2.1) shows stable coexistence of all the populations (see Figure 9) but at $\theta = 1.8$, we observe that the system shows limit cycle oscillations around the coexistence equilibrium (see Figure 10). To get a clear view, we vary $\theta$ from 0 to 6 and draw bifurcation diagram of system (2.1) with respect to $\theta$, Figure 11. We observe that for low values of $\theta$, the system shows stable coexistence of all the populations, but shows limit cycle oscillations around the coexistence equilibrium $e_*$ after a critical value of $\theta$; further increase in the values of $\theta$ leads to instability of the infection-free equilibrium $e_2$. Next, we draw the bifurcation diagram of the system (2.1) by varying the parameter $\gamma_1$ in the interval [0, 0.2], Figure 12. We get that for very low values of $\gamma_1$, the solution trajectories are oscillating around the coexistence equilibrium but the solution trajectories become stable after a threshold value of $\gamma_1$. Now, we vary the force of infection ($\beta$) in the interval [0, 1] and draw bifurcation diagram of the system (2.1), Figure 13. We observe that for very low values of $\beta$, the infection-free equilibrium $e_2$ is unstable, but increase in the values of $\beta$ results in the instability of the coexistence equilibrium, $e_*$. The equilibrium $e_*$ becomes stable for higher values of $\beta$. For very large values of force of infection, the system settles to the zooplankton-free steady state. Next, by varying the values of $K$ in the interval $ ^[5,15] $, we plot bifurcation diagram of the system (2.1), Figure 14. We observe that the stable equilibrium becomes unstable with high-amplitude oscillations with increasing values of $K$. It is apparent from Figure 15 that for large values of input rate of environmental toxins, the zooplankton-free equilibrium, $e_3$, is stable.

Next, we investigate how the stability of the system changes on varying two parameters at a time: $(\theta, A)$ and $(\theta, \beta)$, Figure 16. For this reason, we vary $\theta$ and $A$, and plot stability regions in $(\theta, A)$ plane, Figure 16(a). Here, regions green, red, blue, black and magenta correspond to zooplankton-free, stable coexistence, unstable coexistence, unstable infection-free, and stable infection-free domain, respectively. For low values of $\theta$ ($\theta < 1$) and high values of $A$, we get a zooplankton-free system. By increasing the values of $\theta$, while keeping $A$ at low values, the infection-free equilibrium is unstable while for large values of both $\theta$ and $A$, the infection-free equilibrium is stable. For low values of $\theta$, the coexistence equilibrium is stable while for moderate values of $\theta$, the coexistence equilibrium is unstable. It is apparent from Figure 16(b) that on increasing the values of force of infection, the system becomes zooplankton-free irrespective of the preference of zooplankton for phytoplankton. For low values of $\beta$, on increasing $\theta$, the infection-free equilibrium is stable while on increasing the values of $\theta$ and moderate values of $\beta$, the coexistence equilibrium is stable.

Now, we see the effects of replication factor of free-viruses and input rate of environmental toxins in the aquatic system on the equilibrium abundances of all the densities for the systems (2.1), (3.2), (3.3) and (3.6), Figure 17. The comparisons between equilibrium levels of all the populations for different values of virus replication factor ($b$) and input rate of environmental toxins ($A$) are given in Table 5. We observe that the susceptible phytoplankton is at high equilibrium level if there is no viruses in the system and input rate of toxins is very low; while the susceptible phytoplankton is at low level in the absence of toxins but high value of virus replication factor. The infected phytoplankton is at high equilibrium level if virus replication factor is high, and becomes zero in the absence of free-viruses. The zooplankton population is maximum when the system is free of viruses and toxins but becomes very low in their presence. High values of replication factor of viruses and input rate of toxins maintained the high levels of these two environmental degrading agents. Next, we see the effect of two controllable parameters, decay rate of free-viruses and depletion rate of environmental toxins, on the equilibrium levels of variables in system (2.1), Figure 18. We see that zooplankton are at higher density for high values of these two rate parameters in comparison to the case when these parameters are at low values. Opposite situations occur for other variables of the system. Further, we see the percentage increase/decrease in system's variables by varying the decay rate of free-viruses ($\delta$) and depletion rate of environmental toxins ($\rho$) in the interval [0.01, 0.99], Figure 19. From the figure, it is clear that the concentration of zooplankton changes significantly on increasing the decay rate of free-viruses and depletion rate of environmental toxins; the free-viruses are highly sensitive to its death rate.

5.   Conclusion

In this paper, we proposed a model that incorporates zooplankton's selectivity on susceptible phytoplankton in the presence of infective ones, and seen the combined effects of marine viruses and environmental toxins on the dynamics of phytoplankton-zooplankton interactions. We studied the dynamics of full as well as submodels. We observed that in the absence of free-viruses and environmental toxins in the ecosystem, the system is stable for low values of capture rate of susceptible phytoplankton by zooplankton and carrying capacity of susceptible phytoplankton. But, increase in either of these two parameters, drives the system to an oscillatory state. If the reproductive gain of zooplankton on consumption of susceptible phytoplankton is very low, then the system will become zooplankton-free. We got that the introduction of environmental toxins in the system makes it stable if the system with susceptible phytoplankton and zooplankton is already oscillatory. Moreover, high input rate of environmental toxins results in the extinction of zooplankton population, which is not reported in the previous studies ^[38,39]. The extinction of zooplankton may be due to the fact that environmental toxins increase their death rate. Previously, the destabilizing effects of input rate of environmental toxins and its contact rate with phytoplankton, and the stabilizing effects of depletion rate of environmental toxins have been reported ^[38,39]. Also, we observed that the presence of free-viruses in the ecosystem can stabilize an unstable equilibrium of the $S$$-$$Z$ system. But, increase in the preference of zooplankton on infected phytoplankton results in the oscillation of the system around disease-free steady state, which was not obtained in the previous study ^[6]. We found that for low values of force of infection, the system is oscillatory around the disease-free equilibrium, but the oscillation is around the coexistence equilibrium on increasing the values of force of infection; further increase in the values of force of infection results in stable coexistence of all the populations.

For system (2.1), we observed that for low values of preference of zooplankton on infected phytoplankton, the coexistence equilibrium is stable, but the coexistence equilibrium becomes unstable after a threshold value of preference of zooplankton; further increase in the values of preference of zooplankton on infected phytoplankton leads to instability of the infection-free equilibrium. Note that selectivity behavior of zooplankton was not considered in the previous study ^[6]. On increasing the values of the parameter related to inhibition of growth rate of susceptible phytoplankton due to environmental toxins, the system get stabilized. It is noted that for very low values of force of infection, the infection-free equilibrium is unstable but increase in the values of force of infection leads to the instability of the coexistence equilibrium; for higher values of force of infection a stable coexistence of all the populations is observed. Moreover, the ecosystem becomes zooplankton-free for very high values of force of infection. We observed that the stable equilibrium becomes unstable with increasing the carrying capacity of susceptible phytoplankton, i.e., increasing nutrient supply, which supports the existing hypothesis of the paradox of enrichment ^[53].

The results of sensitivity analysis suggest that a strategy that reduces the parameters having positive PRCC values with either free-viruses or environmental toxins will adequately reduce the density of phytoplankton in the system. Furthermore, a strategy that increases the parameters having negative PRCC values with free-viruses (environmental toxins) will be effective in curtailing free-viruses (environmental toxins) in the system. Identification of these key parameters is important to the formulation of effective control strategies necessary for combating free-viruses and environmental toxins in the aquatic system. Comparison of the equilibrium values of variables in full and subsystems shows that the replication factor of free-viruses and input rate of environmental toxins in the ecosystems can be used as biological controls to combat the reductions in zooplankton and healthy phytoplankton populations. As different parameters show important effects in reducing the concentration of zooplankton, from modeling and analysis perspective, we can conclude that the concentration of zooplankton can be maintained by increasing the decay rate of free-viruses and depletion rate of environmental toxins. For the control of infected phytoplankton in the system, it is imperative to decrease free-viruses from that system by using antiviral reagents ^[54,55]. Bioremediation technology could be very useful to convert the toxigenic compounds to nontoxic products without further disruption to the local environment, which will enhance the persistence and stability of the populations ^[56]. Although environmental toxins are thought not to be entering marine ecosystems in large quantities yet, experimental evidence reveals the fact that phytoplankton populations are highly vulnerable to damage from environmental toxins ^[57]. Sometimes, for high concentration of environmental toxins, complete growth inhibition in phytoplankton was observed ^[58]. So, there might be a possibility of the population crash. But, it can be controlled through the higher depletion rate of environmental toxins. Regulating the depletion rate of environmental toxins could be an effective control strategy to prevent the crash of the aquatic food web system, which could be achieved with suitable human activities such as safe use of pesticides, chemical toxins, awareness among human etc.

Acknowledgments

The authors thank the associate editor and the anonymous reviewers for valuable comments, which contributed to the improvement in the presentation of the paper. The research work of Saswati Biswas is supported by Council of Scientific and Industrial Research, Government of India, New Delhi in the form of Senior Research Fellowship (Ref. No. 20/12/2015(ⅱ)EU-V). Pankaj Kumar Tiwari is thankful to University Grants Commissions, New Delhi, India for providing financial support in form of D. S. Kothari post-doctoral fellowship (No.F.4-2/2006 (BSR)/MA/17-18/0021). The work of Yun Kang is partially supported by NSF-DMS (1313312 & 1716802); NSF-IOS/DMS (1558127), DARPA (ASC-SIM Ⅱ), and The James S. McDonnell Foundation 21st Century Science Initiative in Studying Complex Systems Scholar Award (UHC Scholar Award 220020472).

Conflict of interest

The authors declare that there is no conflict of interest in this paper.