Dynamical Study of an Eco-Epidemiological Delay Model for Plankton System with Toxicity.

In this paper, we analyze the complexity of an eco-epidemiological model for phytoplankton-zooplankton system in presence of toxicity and time delay. Holling type II function response is incorporated to address the predation rate as well as toxic substance distribution in zooplankton. It is also presumed that infected phytoplankton does recover from the viral infection. In the absence of time delay, stability and Hopf-bifurcation conditions are investigated to explore the system dynamics around all the possible equilibrium points. Further, in the presence of time delay, conditions for local stability are derived around the interior equilibria and the properties of the periodic solution are obtained by applying normal form theory and central manifold arguments. Computational simulation is performed to illustrate our theoretical findings. It is explored that system dynamics is very sensitive corresponding to carrying capacity and toxin liberation rate and able to generate chaos. Further, it is observed that time delay in the viral infection process can destabilize the phytoplankton density whereas zooplankton density remains in its old state. Incorporation of time delay also gives the scenario of double Hopf-bifurcation. Some control parameters are discussed to stabilize system dynamics. The effect of time delay on (i) growth rate of susceptible phytoplankton shows the extinction and double Hopf-bifurcation in the zooplankton population, (ii) a sufficiently large value of carrying capacity stabilizes the chaotic dynamics or makes the whole system chaotic with further increment. © Shiraz University 2021.

Introduction

Viral infection in planktonic species affects the bloom dynamics and causes behavioral as well as other changes in the aquatic and marine systems. The capability of regulating the plankton dynamics is still far from understanding. Algal viruses perform a remarkable component in the evolutionary driving force of the aquatic system and responsible for biogeochemical cycles across all the microbial communities. It effectively accounts for the idealization of mathematical biology that handles the new ecological and epidemiological challenges. Viruses have notable prospects as mortality agents for phytoplankton and play a dominant role in extinction and survival behavior among all the planktonic species. Several new developments concerned with the dynamical and behavioral complexity of the prey–predator system have been addressed in the area of ecology and epidemiology (Anderson and May 1986; Das and Chattopadhyay 2015; Thakur and Ojha 2020a). Mathematical models are facilitated to demonstrate the qualitative functioning of the prey–predator system and help to examine the long-term relationship among interacting species of the ecosystem. Many authors assumed only one infected population in their model system, i.e., either prey or predator population is infected, whereas others assumed as both the populations are infected (Gao et al. 2020c; Goyal et al. 2020; Singh et al. 2018; Atangana 2017, 2020, 2018). A number of new observations by using fractional derivative operators have been addressed by Atangana (2017), Atangana (2018), Atangana (2020), Cattani (2018), Gao et al. (2020d), Gao et al. (2020b), İlhan and Kıymaz (2020) and Cattani and Pierro (2013). Eco-epidemiological models are also considered to describe the coronavirus pandemic that describes the real phenomena (Gao et al. 2020a, e). The pioneering work of Kermack and McKendrick (1927) has established a classical SIR model (susceptible, infectious, recovered model) with the idea of plankton disease to study how the population is influenced by infection. Beltrami and Carroll (1994) developed an eco-epidemiological model based on prey–predator in which the prey population seems to be infected by viral contamination and forms an infected group. They found that the system has been destabilized by a minute amount of infection agents otherwise stable tropic configuration noticed. Gakkhar and Negi (2006) studied the role of viral infection and toxic substances on the plankton system and concluded that higher infection rates controlling the plankton blooms. Dhar and Sharma (2010) presented a phytoplankton dynamics along with viral infection and incubation class and found that the absence of incubation class makes the phytoplankton system unstable, whereas the presence of incubation class in the form of delay makes the phytoplankton system stable.

A good number of reviews are available on prey–predator dynamics with disease and infection and also their possible ecological and biological impact specified in Biswas et al. (2010), Saifuddin et al. (2016) and Zhao and Jiang (2014). Upadhyay et al. (2008) proposed an eco-epidemiological model based on the Salton sea which contains an infected fish population and tried to explain all the possible ways to the existence of chaos in a detrimental wetland ecosystem. Das et al. (2016) focused on a phytoplankton–zooplankton model system with virally infected species and studied the essential features of plankton dynamics by taking two important parameters, i.e., mortality of phytoplankton and viral infection of zooplankton. Auger et al. (2009) modeled a predator–prey system by using simple Lotka–Volterra equations with disease-affected predator population is considered. Tannoia et al. (2012) discussed a system for transmissible diseases which is disseminating among predators and found that the persistence of oscillation behavior for the system. Bairagi et al. (2007) investigated a comparison-based study of a model with an infected prey–predator population where the predator response function is governed by three different responses function. They observed that when prey is affected with a disease, then species coexistence is not possible whereas some diverse outcomes yield. Later, many investigations have been made by numerous authors to study the eco-epidemiological model in different ecological scenarios where the population is influenced by external toxicity, external disease transmission, prey refuge, Allee effect, etc. (Biswas et al. 2016; Hethcote et al. 2004; Kumar et al. 2019; Venturino 2002). Among them, the toxin-producing phytoplankton–zooplankton system has played a prominent role in marine as well as freshwater ecosystems. Toxicity not only maintains the high diversity of phytoplankton but also reduces the grazing pressure of zooplankton. Some observations show that toxin substance has huge consequences on zooplanktons growth and affects the whole plankton community (Saha and Bandopadhaya 2009; Thakur et al. 2016). Rapid appearance and disappearance of algal bloom are not well understood till now but several scientific inquiries have been done for the virally infected plankton system under the consideration of toxin-producing phytoplankton together with other controlling parameters to understand the significance of toxicity on the survival scenario of plankton. Chattopadhayay et al. (2002a, b) have divulged that the toxicants released by phytoplankton may perform an essential factor for the action of termination of the plankton bloom and act as a bio-controller. Chakraborty and Das (2015) have established two zooplankton and one phytoplankton model system by using Holling type II functional responses that focuses on the mechanism of toxicity with constant harvesting. Upadhyay and Chattopadhyay (2005) have considered a TPP-zooplankton-fish model and observed that TPP in aquatic systems plays a pivotal role by showing stabilizing contribution in population dynamics. Huang et al. (2018) have studied the aquatic toxicity by a simple mathematical model and concluded that the high toxin can lead to the extirpation of the population.

A well-known truth is that time delays exist in every biological process and influence the dynamics of the aquatic as well as the marine ecosystem and its whole community. It can be seen in the form of gestation, maturation, incubation, resource regeneration, etc., into many population models. It reflects the complexity in such a model by showing stability transition phenomenon for equilibria, the occurrence of Hopf-bifurcation, chaotic oscillations and extinction dynamics (Ojha and Thakur 2020; Meng et al. 2011; Wang et al. 2014). One such interesting phenomenon due to time delay is the double Hopf-bifurcation in which the system possesses two critical values of time delay (Jiang et al. 2015). Recently, Thakur and Ojha (2020b) have investigated a three-tier delayed plankton-fish model system and considered the effect of delay on the extinction behavior of fish. They obtained that the increasing value of time delay helps to stabilize the fish dynamics after the extinction through the double Hopf-bifurcation phenomenon. Time delay has an enormous capacity to change system dynamics. It can make the unstable system stable by showing its positive impact, whereas it can also be responsible for making the whole system chaotic (Sharma et al. 2014; Zhang and Rehim 2017). Beretta and Kuang (2002) have substantially studied the biological models by incorporating the effect of time delay. Zhang and Rehim (2017) dealt with a toxic-phytoplankton–zooplankton system in the context of plankton toxicity and found some excitable results of time delay. Sharma et al. (2014) have analyzed the effect of toxin liberation with time delay in the nutrient-plankton system in which planktonic bloom showed oscillation through chaotic behavior. Thakur et al. (2020) have established the plankton-fish model under the consequence of multiple gestation delays and demonstrated that the two equal gestation delays may promote the chaotic phenomenon in the plankton system.

Recently, many studies considered the delay induced infection models to investigate the combined effect of infection and time delay on population models. Gakkhar and Singh (2010) proposed a delayed plankton model, and the results showed that infection in a population may ensure the complexity as well as beneficial for controlling the chaos. It is also very interesting to point out that fluctuation in the toxic parameter pays uniform importance toward the chaos. Furthermore, Agnihotri and Kaur (2019) established a delayed plankton system with viral infection and noticed that system does not get any meaningful behavior without delay but after introducing a time delay in toxin liberation process, the system switches its behavior from stable to unstable periodic oscillations. Meng et al. (2018) considered a delayed eco-epidemiological model under the assumption of predator’s existence in the system where disease disseminates in the prey populations. They found that time delay in the incubation period is harmful to the transmission of infectious diseases. Biswas et al. (2017) discussed optimal harvesting policy and Allee effect for a prey–predator eco-epidemiological model with time delay. They have briefly studied the bistability between all the equilibrium points. Further, the authors concluded that beyond the bifurcation point of time delay, the system produces complex dynamical outcomes together with chaos. Xu and Zhang (2013) explored the global dynamics of disease transmission in a predator population with gestation time delay. It has been shown that the disease-free equilibrium leads to the existence of Hopf-bifurcation with the accumulated value of time delay.

With the above motivation, we have considered two interacting components consisting of toxin-producing phytoplankton and zooplankton population in the presence of virally infected phytoplankton. A discrete-time delay is incorporated in the form of incubation period of infected phytoplankton dynamics. Further, we have assumed zooplankton predate both susceptible and infected phytoplankton. This fact has been exhibited through Holling type II function form for the model system. Our main objective of the present article is to study the consequences of some important parameters like toxicity, viral infection, and incubation time delay on the prey–predator dynamics where the prey population is infected by some diseases. Since every species may avoid becoming prey for other living species and try to escape predation by adopting some defensive scheme. Toxin-induced anti-predator defense by phytoplankton is taken that decreases the predation rate as well as the density of zooplankton. It is also considered into account that viral infection is not an instantaneous process, rather than susceptible phytoplankton takes some time lag to convert into infected class. The paper is organized into the following manner: In Sect. 2, the formulation of the plankton system in the presence of viral infection is briefly presented. In Sect. 3, mathematical analysis such as boundedness solution, the existence of all positive equilibria and local stability of the virally infected model are discussed. The model system with time delay is described in Sect. 4 followed by local stability analysis of the delayed system in this section. Using the application of normal form theory and center manifold theorem, the stability and direction of Hopf-bifurcating solution of the delayed system are derived in Sect. 5. Numerical validation of the delayed and non-delayed systems is given in Sect. 6. Conclusions are carried out in the last section of the paper.

Model Formulation

In this section, we consider a plankton interaction model which is made up of two interacting components such as phytoplankton and zooplankton. We assume that the phytoplankton has the capacity to produce toxins as well as some of them infected with the viral. At any time t, phytoplankton is represented by N(t) and zooplankton is represented by Z(t). We would like to impose a brief description of the model system as follows: Implementing the above assumptions gives the following phytoplankton–zooplankton model system as:with associated initial conditions are , and . Eq. (2.1) describes the interaction among susceptible phytoplankton, infected phytoplankton and zooplankton. Susceptible and infected phytoplanktons are taken as prey components. We assume the initial infected phytoplankton population and their growth at any time t are due to viral infection. It is also assumed that infected populations are being recovered at a constant recovery rate . The susceptible phytoplankton grows logistically and per capita growth rate approaches to zero as population size approaches a maximum imposed by limited resources in the environment, whereas the zooplankton population predates the susceptible and infected phytoplankton populations at rates and , respectively. Here, , as infected phytoplankton are more vulnerable to predation than susceptible phytoplankton. Noctiluca scintillans is a marine planktonic dinoflagellate. Although this species is not toxin creating, it has been found to accumulate the poisonous levels of ammonia, which consecutively act as the killing agent in blooms (Okaichi 1976). Experimental works demonstrate that in German Bight disease, Noctiluca scintillans become unable to reproduce (Uhlig and Sahling 1990). The impact of toxins decreases the growth of zooplankton species as per the Holling type II functional response. It is investigated that the toxin liberation parameter plays an important role in the occurrence of bloom. A brief description of parameters and their units of the model system (2.1) are presented in Table 1.

Table 1

Brief description of the notation used for parameters with their units

Parameters	Units	Description
r	h^-1	Intrinsic growth rate of susceptible phytoplankton population
K	l^-1	Carrying capacity of phytoplankton population
C	lh^-1	Rate of infection
η	lh^-1	Maximum predation rate of susceptible phytoplankton population
γ	l^-1	Rate at which infected phytoplankton become susceptible population
ρ	lh^-1	Maximum predation rate of infected phytoplankton population
δ	lh^-1	Total death rate of infected phytoplankton population due to disease
g	lh^-1	Growth rate of zooplankton due to predation of susceptible phytoplankton population
h	lh^-1	Growth rate of zooplankton due to predation of infected phytoplankton
d	lh^-1	Total death rate of zooplankton population
θ	h^-1	The rate of toxin liberation by the toxin producing phytoplankton population
β	l^-1	Half saturation constant for the toxin producing phytoplankton
α	l^-1	Represent the half saturation constant measures the extent to which the environment
		provides protection to susceptible and infected phytoplankton population

We assume that in the absence of viral infection, phytoplankton N(t) grow logistically with carrying capacity K and growth rate r which is expressed as

In the presence of viral infection, the total phytoplankton density can be classified into two categories, S(t) represents susceptible phytoplankton density, and I(t) represents infected phytoplankton density such that

We assume that only susceptible phytoplankton S(t) has the capability to reproduce where the infected phytoplankton does not have the capability to reproduce. However, the infected phytoplankton still contributes to the growth of phytoplankton population toward the same carrying capacity K. Therefore,

Zooplankton predates susceptible phytoplankton with respect to Holling type II functional form, and susceptible phytoplankton becomes infected when they come in contact with infected phytoplankton. Viruses are transmitted by the law of mass action with an infection rate of C. We also assume viruses only infect the phytoplankton population. Moreover, infected phytoplankton does recover with the recovery rate of . Therefore,

The interaction of susceptible phytoplankton with the infected phytoplankton increases the density of infected population whereas recovered ones decrease the density of infected population. Infected phytoplankton will die out their natural death rate . Assuming that zooplankton also predates infected phytoplankton according to Holling type II functional form, the infected phytoplankton dynamics becomes

The consumption of susceptible and infected phytoplankton gives the growth in zooplankton density according to Holling type II functional response at the rate of g and h. In the absence of phytoplankton, zooplankton will die out at the rate of d. Additionally, we assume the extra mortality in zooplankton due to toxin liberation by phytoplankton at the rate of . With this assumption, zooplankton dynamics becomes

Analysis of the Model System (2.1)

Boundedness of Solution

Theorem 1

All solutions of the system (2.1) which initiate in are uniformly bounded for suitably taken positive and .

Proof

To examine the boundedness of the model system (2.1), let us consider a functionDifferentiating with respect to time t, we haveFrom Eq. (2.1), we haveIntroducing positive and rewriting, we haveThis impliesChoosing and ,Using differential inequality theorem in Eq. (3.3), we getwhere .

Moreover,Hence, all the solutions initiating in of Eq. (2.1) are defined in the regionThus, all positive solutions of virally infected model system (2.1) are uniformly bounded for all .

Existence of Equilibria

In this subsection, we have obtained all possible feasible equilibrium point of the model system (2.1).

Existence of . The trivial equilibrium always exit.

Existence of . From Eq. (2.1), it can be seen that . Therefore, positive equilibrium exits.

Existence of . From Eq. (2.1), the susceptible infected equilibrium exits if we get and positive. Thus, we have Solving Eqs. (3.4) and (3.5), we get Hence, the equilibrium point exits if .

Existence of . From Eq. (2.1), the infection-free equilibrium exits if we get and positive. Thus, we have where Hence, the equilibrium point exits if , and .

Existence of non-trivial equilibrium point , where , and satisfy the following algebraic nonlinear equations, let Simplifying Eq. (3.8), we get the following quadratic equation where Roots of above quadratic equation are given by , from which at least one root will be positive provided one of the following condition holds: From Eq. (3.7), we have Clearly, exists if . Now, the value of substituted in Eq. (3.6), we get Clearly, exists if , and .

, ,

, and ,

, .

Stability Analysis

In this subsection, we have discussed the stability condition around all the existing equilibrium point. Hence, we have computed the variational matrix around the positive equilibria, we getwhereWith the help of some straightforward calculation, we construct the following theorem:

The trivial equilibrium point exists. The Jacobian matrix corresponding to the equilibrium point is given by The eigenvalues around are . Therefore, the equilibrium point is stable manifold along direction and unstable manifold along S-direction.

The susceptible phytoplankton equilibrium point 0, 0) exists. The Jacobian matrix corresponding to is given by The eigenvalues around are . Therefore, is locally asymptotically stable provided Thus, the equilibrium point is locally asymptotically stable provided the conditions (3.9) and (3.10) hold. will be saddle point if

The zooplankton-free equilibrium point exists. The Jacobian matrix corresponding to the point is given by The characteristics equation of the Jacobian matrix is given by In the characteristics Eq. (3.11), the quadratic term gives two negative eigenvalue whereas third eigenvalue of the characteristics Eq. (3.11) is , which is either negative or positive, it is depends on or Thus, equilibrium point is locally asymptotically stable if the condition (3.12) holds. Otherwise, is a saddle point if condition (3.13) holds.

The infection-free equilibrium point exists. The Jacobian matrix corresponding to the point is given by One eigenvalue for is and another eigenvalue is roots of sub-matrix The eigenvalue of sub-matrix has negative eigenvalues if tr and det. Therefore, Now, tr if the following condition holds: and det if the following condition holds: Thus, the infection-free equilibrium point is locally asymptotically stable if conditions (3.14) and (3.15) hold.

Variational matrix around the non-trivial equilibrium point is given by where Here, , and . The characteristics equation of the Jacobian matrix is given by where If we assume , , and , then it is easy to verify that , and . Thus, from Routh-Hurwitz criteria, the model system (2.1) is locally asymptotically stable around that equilibrium point .

Theorem 2

The model system is locally asymptotically stable at the interior equilibrium point if the following conditions are satisfied.

,

,

,

,

,

.

Theorem 3

Assume that the positive equilibrium point is locally asymptotically stable in int. Then it is globally asymptotically stable in a subregion of int, providedwhere and .

,

,

,

.

Let us assume a positive definite Lyapunov function around a positive equilibrium point :Differentiating with respect to time t along the solution of the system (2.1), we obtainThe above-mentioned equation can be represented as the sum of quadratics as,Sufficient conditions for to be negative definite areUnder the condition of Theorem 3, we can justify all the above conditions.

Hopf-Bifurcation Analysis

In this subsection, we illustrate the conditions for the Hopf-bifurcation of the system around , when the parameter passes through some critical value. In fact, we have the following theorem.

Theorem 4

The system (2.1) undergoes for Hopf-bifurcation around the positive equilibrium point , whenever the parameter of the system crosses a critical value . If the following conditions are satisfiedthen dynamics of the system (2.1) experiences a Hopf bifurcation around the positive equilibrium points .

, , ,

,

,

Under the conditions of Theorem 2, the system is locally asymptotically stable around equilibrium point . We choose b as a bifurcation parameter and if there exists a critical value says , then , , and . For the occurrence of Hopf-bifurcation at , the characteristics equation must be in the formwhich has three roots , and .

Now, we have to show transversality conditions .

Let , substituting value of in Eq. (3.20) and differentiating with respect to b, we getHere,We havewhere , and .

Hence, the transversality condition holds. This implies that Hopf-bifurcation occurs at .

Model System (2.1) with Time Delay

In this section, we have assumed that the viral dissemination among phytoplankton is not an immediate procedure rather some time lag is essential in the form of an incubation period where the susceptible stay for some definite period of time and after that leaving the susceptible class and joining the infected class. With this purpose, we have introduced discrete-time delay in the infected phytoplankton dynamics. Hence, the corresponding model system takes the following form:with the following initial conditionswhereTo investigate the local stability of system (4.1) about coexistence equilibrium point , we introduce the following perturbations , , . Then, the transformed system of equations is given bywhereThe characteristic equation of system (4.2) at coexistence equilibria is represented aswhere

Lemma 1

If , then Eq. (4.3) has a pair of purely imaginary roots if satisfies

Now, separating real and imaginary parts of Eq. (4.4), we haveEliminating trigonometrical function, this leads to the following sixth degree equationChoosing , this implieswhereNow, letthereforewhere the roots of can be expressed as,From the above-mentioned equations, we have the following lemma (Song and Wei 2004):

Lemma 2

Since , therefore

If then has no positive roots.

If then has exactly two positive roots and (where ). Furthermore, and .

Putting (q=1,2) and solving Eq. (4.5) for , we get the corresponding critical value of time delay as,where and .

If then is the purely imaginary roots of the characteristic Eq. (4.5) and other roots having the nonzero real part.

Let us define where

Theorem 5

We have the following transversality conditions

Proof

Differentiating the characteristic Eq. (4.5) with respect to , we obtainNow, substituting in Eq. (4.13) such that and , it givesIf then and . This is the required transversality condition for occurrence Hopf-bifurcation.

Properties of Periodic Solution of Delayed System (4.1)

In this section, we will elaborate the direction, stability and period of the periodic solutions of the model system (4.1) which is bifurcating from the coexistence equilibrium . We use the normal form and center manifold theory introduced by Hassard et al. (1981) to compute the direction of Hopf-bifurcation which guarantees that the periodic solution exists. Further, we assume without any loss of generality that the critical value of time delay is denoted by at which a pair of purely imaginary roots crosses the imaginary axis and the system undergoes Hopf-bifurcation at the coexistence equilibrium .

Let , , , where and at gives the value of Hopf-bifurcation. Rescale the time by , the system (4.1) can be written into the following continuous real-valued functions as where and are given, respectively,such thatandwhere ,According the Riesz representation theorem, there exists a function of bounded variation for , such thatIn fact, we can pickwhere is the Dirac delta function.

For , defineandNow, the system (5.1) is equivalent towhere for

For , defineand a bilinear inner product is given bywhere . Then, and are adjoint operators. From the above calculation, we know that and are eigenvalues of A(0) and , respectively.

Suppose be the eigenvector of A(0) corresponding to the eigenvalue thenFor , Eq. (5.11) reduces towhere I is identity matrix. This implies,which givesSimilarly, let be the eigenvector of corresponding to the eigenvalue thenwhere I is identity matrix. This implies thatUnder the normalization condition , we havewhich givesNow, following the same manner as given in Hassard et al. (1981), we obtainwhich can also be written asSimplification of Eq. (5.15) giveswhereWe denote and . By computing, we obtainedandwhere and are the constant vectors, which can be determined by using Hassard et al. (1981). Therefore,Solving this system for , we havewhereSimilarly,Solving this system for , we havewhereConsequently, the values of and are determined from Eqs. (5.17) and (5.18). The value of is represented by delay and parameters by Eq. (5.16). Hence, different values are computed as:The above values have been used to determine the bifurcating periodic solution in center manifold at the critical value .

Theorem 6

From Eq. (5.19), the following results have been obtained: Where , and have their usual meanings.

Direction of the Hopf-bifurcation is supercritical (subcritical) if .

Stability of the bifurcating periodic solution is stable (unstable) if .

Period of the bifurcating periodic solution increases (decreases) if .

Numerical Simulation

In this section, the complex dynamical behavior of the system (2.1) and (4.1) with the help of time series, phase portrait and bifurcation diagram is investigated to substantiate our theoretical findings. The main motive of this section is to detect the possible impact of all the sensitive parameters numerically. For this purpose, we have considered the following fixed set of parameter values based on the biological principle:We have observed that the model system (2.1) shows various kinds of dynamics for the different parametric values of carrying capacity K. One can see the increasing value of K from 220 to 300, the system shows stable focus to chaotic attractor through the periodic oscillation of order-1 and order-2. When , the positive equilibrium point converges to an asymptotically stable state, verified by Routh-Hurwitz criterion () and shown by Fig. 1a. As we increase the value of carrying capacity K, one can see that stable equilibrium becomes unstable. Time evolution and phase space representation of limit cycle attractor are given in Fig. 1b for . Similarly, a period-doubling attractor is shown in Fig. 1c for . An excess value of K promotes the occurrence of chaos in the plankton system. Figure 1d shows the time evolution and phase space representation of chaotic attractor when . Under the bifurcation analysis, we have plotted the bifurcation diagram between carrying capacity K and population density S, I and Z. In Fig. 2, successive maxima of S, I and Z are observed in the range [0.0, 300], [0.0, 160] and [5, 30], respectively, as carrying capacity K is taken in the range and the other parameters are the same as given in Eq. (6.1). It is observed that interior equilibrium of the virally infected model system (2.1) is locally asymptotically stable when and unstable via occurrence of Hopf-bifurcation when parameter K passes the threshold value . These results indicate that the small value of carrying capacity is beneficial for a virally infected model system because it makes the system stable. In addition, susceptible phytoplankton grows rapidly, whereas zooplankton grows slowly for large K.

Fig. 1

Time evolution and phase space of model system (2.1) for susceptible phytoplankton, infected phytoplankton and zooplankton at a , b , c , d

Fig. 2

Bifurcation diagram of model system (2.1) for K versus a Max(S), b Max(I), c Max(Z)

In the next scenario, we present the bifurcation diagram for one of the other crucial parameter r. In Fig. 3a–c, successive maxima of S, I and Z are observed in the ranges [40, 220], [20, 90] and [0.0, 25], respectively, as the intrinsic growth rate r is taken in the range and the other parameters are the same as given in Eq. (6.1). A period-doubling cascade is found for the parameter r and the positive equilibrium point is stable when and unstable when , whereas the smaller value of r may lead the extinction in zooplankton. It is also observed that a decreasing value of r switches the stability from plane to plane whereas an increasing value of r leads the periodic solution of order-2.

Fig. 3

a–c Bifurcation diagram of model system (2.1) for r versus a Max(S), b Max(I), c Max(Z)

To understand all the possible aspects of toxin substances release by toxin-producing phytoplankton, we have drawn the time series and phase space diagram for different values of . When , the system has chaotic behavior as shown in Fig. 4a. Periodic behavior is observed at (c.f., Fig. 4b). When is in the interval [6.22, 6.74], a stable focus is noticed in Fig. 4c. For further increment in the toxin rate, a stable focus in the plane is observed as shown in Fig. 4d, and zooplankton goes to extinction. So it is clear that strengthening the toxin liberation shows a stabilizing effect in plankton dynamics. Justification of Fig. 4, we have plotted the bifurcation diagram between toxin liberation rate and successive maxima of S, I and Z in the range [50, 250], [0.0, 140] and [0.0, 25], respectively (c.f., Fig. 5). The amplitude of oscillation decreases as the value of increases.

Fig. 4

Time evolution and phase space of model system (2.1) for susceptible phytoplankton, infected phytoplankton and zooplankton at a , b , c , d

Fig. 5

Bifurcation diagram of model system (2.1) for versus a Max(S), b Max(I), c Max(Z)

In Figs. 6 and 7, we have plotted the bifurcation diagram for infection rate C versus successive maxima of S, I, Z and recovery rate of infected phytoplankton versus successive maxima of S, I, Z, respectively. A simple Hopf-bifurcation scenario is found for both the parameters but in different ranges. It is noticeable that a small infection rate strengthens the phytoplankton density but weakens the zooplankton density whereas a large recovery rate strengthens the phytoplankton density but weakens the zooplankton density. So a similar kind of behavior is observed for the parameters C and but in the opposite way.

Fig. 6

Bifurcation diagram of model system (2.1) for C versus a Max(S), b Max(I), c Max(Z)

Fig. 7

Bifurcation diagram of model system (2.1) for versus a Max(S), b Max(I), c Max(Z)

Now, we study the effect of time delay on model system (2.1), i.e., model system (4.1) corresponding to different values of r while considering the other parameters are the same as given in Eq. (6.1). In Fig. 8a, we have presented the time evaluation of model system (4.1) with which shows stable focus at . In Fig. 8b, we have presented the bifurcation diagram of the model system (4.1) for versus population density S, I, Z at in which the population of S and I bifurcates when the value of time delay increases. The obtained result shows that after crossing the critical value of , i.e., , the population of phytoplankton becomes unstable. These results state that small delay affected the population of S and I by changing their stability but not able to avoid the extinction in Z. If we choose , the model system (4.1) shows stability behavior for all the species at (c.f., Fig 8c). The corresponding bifurcation diagram for varying and fixed is presented in Fig. 8d. It is observed that all the species lose their stability at after crossing the critical value of time delay by showing limit cycle behavior in phytoplankton and extinction behavior in zooplankton. This result indicates that time delay is not able to avoid extinction in zooplankton whereas another parameter supports extinction. Similarly, if we choose , then the model system (4.1) shows limit cycle behavior at (c.f., Fig 8e). In Fig. 8f, we have presented the bifurcation diagram of the model system (4.1) for versus population density S, I, Z at . Double Hopf-bifurcation scenario is noticed while we increase the value of in the range . It is observed from Fig. 8f that the model system (4.1) shows limit cycle behavior initially and after crossing the first critical value of , i.e., , stability occurs and again after crossing the second critical value of , i.e., , the whole system becomes unstable. If we choose , obviously at that state system is chaotic (c.f., Fig. 1d). We have generated the bifurcation diagram for the model system (4.1) with varying time delay parameter with the fixed value of and observed the small interval of changes the chaotic state to stable state and vice versa via limit cycle (c.f., Fig. 9).

Fig. 8

Time evolution of model system (4.1) for susceptible phytoplankton, infected phytoplankton and zooplankton with at a , c , e , bifurcation diagram of model system (4.1) for versus population density S(t), I(t), Z(t) at b , d , f

Fig. 9

Bifurcation diagram of model system (4.1) for versus a Max(S), b Max(I), c Max(Z) at

Discussions and Conclusions

In the last decades, many researchers have focused on the study of plankton dynamics with assumptions in which the role of infective population has been widely recognized. Viruses are one of the most important factors in the aquatic ecosystem that may cause complexity when it interacts with other planktonic species. Another prominent parameter is the mortality of zooplankton due to the toxicity of phytoplankton that has an essential feature in the plankton population. A good number of studies have been characterized that toxin substances released by phytoplankton are one of the major causes in the reduction of grazing power of zooplankton and it also affects the growth of each species of the food chain. The considered model system in our study is motivated by Gakkhar and Negi (2006); Agnihotri and Kaur (2019). Gakkhar and Negi (2006) considered an infectious plankton dynamics in Lotka–Volterra form with toxin liberation process, where the authors conjectured the Holling type I functional form for plankton interaction and Holling type II functional form for studying the consequences of toxin on zooplankton dynamics. Agnihotri and Kaur (2019) extended the work done by Gakkhar and Negi (2006) by taking one crucial parameter, the recovery rate of infected phytoplankton. Further, Holling type I functional form is contemplated for reducing the grazing pressure of zooplankton due to toxicity and assumed time delay in toxin liberation term. We have extended the work done by Agnihotri and Kaur (2019) by using the Holling type II functional form. In this work, we have taken Holling type II functional form for the species interaction as well as toxin liberation term and incorporated all the critical parameters that make the system more realistic. Biologically, Holling type II functional form gives the best description of plankton interaction rather than Holling type I. Additionally, the assumption of two important parameters, i.e., infection rate and recovery rate both in the same model fascinates the plankton system. We have proposed a delay-induced toxin-producing phytoplankton and zooplankton model system with the viral infection where the phytoplankton population is divided into two categories which depend on whether they carry the infection or not. The one who carries the infection forms an infective group and known as infected phytoplankton, and the remaining are known as susceptible phytoplankton. This model consists of Holling type II predation rate of zooplankton over the susceptible phytoplankton as well as infected phytoplankton and also the same type of functional form is considered for the extra reduction in zooplankton due to toxicity. In addition, a feedback time delay is assumed in the infected phytoplankton dynamics due to the incubation period which makes the biological relevance of our model system. Meng et al. (2018) also considered the same kind of time delay in their study with harvesting in the prey population. So the present study is also inspired by the research of Meng et al. (2018) that the affection of incubation time delay makes the system more appropriate to the real biological environment. Analytically, we have investigated the boundedness criteria, existence of all feasible equilibria and local as well as global stability conditions of the delay-free system. It is observed that there exists a smallest critical value of in delayed system beyond that stability switches and periodic oscillations appear. So the local stability of the system with time delay is calculated. Properties of periodic solutions are also examined. Numerical simulations have been performed to explore the impact of model parameters on the system dynamics. Based on the simulation results, one can address that the following parameters have a significant impact on our study: From Fig. 1, we have found that carrying capacity plays a significant role in plankton dynamics by showing different types of attractors including stable focus, periodic solution of order-1, order-2 and chaotic attractor. Further, it is remarkable that a high value of carrying capacity increases the oscillation size of phytoplankton and finally produces chaos which is presented in Fig. 2. From Fig. 3, a decreasing value of r leads the extinction in zooplankton where the increasing value of r leads the period-doubling cascade. We have observed that the small value of growth rate r makes the interior equilibrium stable to unstable. A period halving bifurcation scenario is noticed for the parameter , i.e., toxin liberation rate of phytoplankton (c.f., Fig. 5). These findings indicate that the parameter has a positive impact on the plankton system but when crosses some threshold value, extinction in zooplankton is observed. From Fig. 6, one can see that a minute change in infection rate C not only stabilizes the dynamics but also beneficial for susceptible phytoplankton to increase their density whereas harmful for zooplankton to extinct their density. The small decrement of the infection rate stabilizes the dynamics by showing the weak negative effect in the dynamics of zooplankton. Similarly, we have explored the dynamics with respect to the recovery rate of infected phytoplankton in which a small increment in recovery rate leads the extinction in zooplankton and all the species coexist in periodic mode at the low value of recovery rate (c.f., Fig. 7).

Growth rate of phytoplankton (r).

Carrying capacity of phytoplankton (K).

Toxin liberation rate .

Infection rate of phytoplankton (C).

Recovery rate of infected phytoplankton .

Time delay .

Further, we have studied the impact of time delay for different values of the growth rate r of phytoplankton on the system (4.1). As we know, time delay has an enormous capacity to change the stable equilibria to unstable and vice versa. In our study, we have noticed the stabilization as well as the destabilization effect of time delay when we introduce time delay on infected phytoplankton dynamics. It has been observed that the time delay may destabilize the zooplankton-free equilibrium and cause the convergence of coexistence equilibrium to zooplankton-free equilibrium of plankton system. In Fig. 8b, d, extinction is observed in zooplankton dynamics for and with increasing value of time delay also which shows the negligence impact of time delay in zooplankton extinction behavior whereas phytoplankton density is bifurcated after crossing the threshold value of time delay. In addition, a double Hopf-bifurcation scenario is observed for varying time delay at and this observation is depended on the parametric value of r only (c.f., Fig. 8f). In the absence of time delay, chaos has been observed for (c.f., Fig. 1d) introducing the small value of time delay on this dynamics, one can see the stable focus after the periodic oscillation and finally the system becomes chaotic again for further increment in (c.f., Fig. 9). The system shows the positive impact of time delay that stabilizes the oscillatory and irregular dynamics of the proposed system. If the time delay is varied by some critical value, the interior equilibrium of the system switches from stable to unstable and vice versa. It is remarkable that the occurrence of double Hopf-bifurcation depends on the value of r because as we take the value of r as 3 or smaller than 22, no double Hopf-bifurcation scenario is observed. A similar kind of effect is also observed for the parameter K with time delay. Finally, we conclude that the plankton system with the viral infection has rich dynamics and a number of parameters that are unpredictable and complex in nature. Parameters like carrying capacity and toxin liberation rate of phytoplankton generate chaos, whereas some parameters like infection and recovery rate stabilize the system.