Qualitative analysis of a discrete-time phytoplankton–zooplankton model with Holling type-II response and toxicity

The interaction among phytoplankton and zooplankton is one of the most important processes in ecology. Discrete-time mathematical models are commonly used for describing the dynamical properties of phytoplankton and zooplankton interaction with nonoverlapping generations. In such type of generations a new age group swaps the older group after regular intervals of time. Keeping in observation the dynamical reliability for continuous-time mathematical models, we convert a continuous-time phytoplankton–zooplankton model into its discrete-time counterpart by applying a dynamically consistent nonstandard difference scheme. Moreover, we discuss boundedness conditions for every solution and prove the existence of a unique positive equilibrium point. We discuss the local stability of obtained system about all its equilibrium points and show the existence of Neimark–Sacker bifurcation about unique positive equilibrium under some mathematical conditions. To control the Neimark–Sacker bifurcation, we apply a generalized hybrid control technique. For explanation of our theoretical results and to compare the dynamics of obtained discrete-time model with its continuous counterpart, we provide some motivating numerical examples. Moreover, from numerical study we can see that the obtained system and its continuous-time counterpart are stable for the same values of parameters, and they are unstable for the same parametric values. Hence the dynamical consistency of our obtained system can be seen from numerical study. Finally, we compare the modified hybrid method with old hybrid method at the end of the paper.


Introduction
The study of mathematical models for population dynamics is considered as a key area in abstract ecology from the time when the famous Lotka-Volterra model was presented [1]. The learning of organism movement and spreading has turn out to be a fundamental element for understanding a chain of ecological interrogations associated with the spatiotemporal study of dynamics of populations [2]. Planktons are enormously flexible in abundance, both temporally and spatially. Plankton variability depends on natural along with physical procedure for the spatial structure. Natural processes include, for instance, development, grazing, and behavior, and physical procedures include, for instance, mixing and lateral stirring. Nonlinearity of ecosystems entirely contributes to the spatial organization in plankton allocations [3]. In marine ecology the word plankton refers to the spontaneously moving and faintly swimming organisms. Commonly, plankton is parted into two species, the phytoplankton species and zooplankton species. Phytoplankton species are tiny in their size with a single celled structure [4]. Phytoplankton are beneficial for aquatic life and produce half of the oxygen in the world through the process of photosynthesis. Phytoplankton population is exerting universal-scale effect on atmosphere by transporting CO 2 from water surface to the depth of oceans. Mainly, this process happens due to their death, sinking, and primary production [4]. It is observed that algal species rise abundantly in damped, wet, and marine environments. The stages of speedy growth, slow stagnation, and accelerated decline in the number of cells collectively create an algal bloom. This phenomenon of accelerated variation in the density of phytoplankton population is the central trait in the plankton ecosystem [4]. Despite the fact that the sudden emergence and disappearing of blooms is not clear, the undesirable effect of damaging algal blooms on the health of mankind, aquatic life, and fisheries trade can be easily seen [4].
On the incidence of blooms, phytoplankton and zooplankton interact with each other, and the study of this interaction is the point of focus of many scientific investigations [5]. Phytoplankton produces toxic materials to avert predations by their predators (zooplankton). Furthermore, this is the topic of interest of many researchers from many decades. Mathematical modeling of interactions between plankton species provides us an important optional method in improving the knowledge of any individual related to the biological and physical mechanisms concerning to the ecological study of plankton population [5].
The authors in [6] have considered a plankton-nutrient model related to aquatic environment by consideration of planktonic blooms. In [7] the authors have examined the influence of periodicity and seasonality on planktonic dynamics.
In [8] the authors have presented two mathematical models connected to plankton ecosystem along with a strong representation of viral septic phytoplankton and viruses. The authors in [9] have contemplated the effect of predation on competitory elimination and the coexistence of competitory predators. Moreover, they presented and explored a one-phytoplankton two-zooplankton model along with the consideration of harvesting.
Huppert et al. [10] considered a nutrient-phytoplankton model to examine the dynamical behavior of phytoplankton blooms. In [11] the authors have presented a zooplanktonphytoplankton model with harvesting. Furthermore, they have explained that the extra exploitation may exterminate the population while suitable harvesting guaranties the consolidation of both populations. Moreover, numerous studies have their point of focus on phytoplankton-zooplankton models along with a source of nutrient, the toxic consequence of plankton species, the survival of plankton species, or the harvesting effects [9][10][11][12][13][14][15][16]. It is convenient to introduce the toxin creating lag during the study of the dynamics of phytoplankton-zooplankton models. The authors in [17] have presented a mathematical model including time lag in toxin deliverance by phytoplankton. The work done in [18][19][20][21] motivated us to study the dynamics of a phytoplankton-zooplankton popula-tion model with toxicity. Moreover, the toxic substance is released by phytoplankton and sometimes by other external sources.
We consider the basic phytoplankton-zooplankton model presented by Chattopadhayay et al. [22]. Furthermore, this mathematical model is based on the following conditions.
• We suppose that z(t) and p(t) are the sizes of zooplankton and phytoplankton populations, respectively. • Zooplankton population eats phytoplankton population and then recycles them into their own community. The functional response αp(t)z(t) a+p(t) represents the predation rate of zooplankton population on phytoplankton species. Moreover, this predation increases the growth rate of zooplankton, which is represented by the term βp(t)z(t) a+p(t) . • We assume that zooplankton population becomes infected by eating infected phytoplankton population. Additionally, the infection in phytoplankton may be produced due to external toxic substance (see [22]). • We assume that the infection in phytoplankton may be produced due to external toxic substance (see [22]). • Phytoplankton population has logistic growth [21] in the absence of zooplankton population, where r is their exponential rate of growth, and k is the maximum carrying capacity of environment. Under these conditions we have the following phytoplankton-zooplankton model [22]: (1.1) Kuang [23] have inspected the limit cycle behavior in Gause-type predator-prey systems with Holling type-II response [24]. In addition, he revealed that the study of dynamical properties of predator-prey models using a Holling-type response function is better than the study of dynamics of predator-prey models without using Holling response. Generally, Holling type-II response is modeled and described by using rectangular hyperbola, and its mathematical form is given as where a is any constant. By using Holling type-II response we get the following mathematical form of system (1.1): (1.2) • Next, we assume that the time lag for production and mediation of toxic substance by phytoplankton is zero. • We introduce the catchability coefficients q 1 and q 2 for phytoplankton and zooplankton populations respectively. Generally, functional form for harvesting is expressed by using the hypothesis of catch-per-unit-effort [25]. • Moreover, we introduce E as the parameter for combined effort for harvesting of population [25].
Under these modifications, system (1.2) takes the following mathematical form: where the parameters in system (1.3) are nonnegative and defined as follows: a: constant of partial capturing saturation. α: maximal takeover rate of zooplankton on phytoplankton. β: conversion rate of phytoplankton-zooplankton (β < α). ρ: toxicity rate of phytoplankton per unit biomass. δ: natural rate of death of zooplankton population. Moreover, the term m 1 p 3 (t) appearing in system (1.3) represents the infection produced in phytoplankton population due to an external toxic substance. In addition, d 2 dp 2 (m 1 p 3 ) = 6m 1 p > 0 shows an accelerating growth of toxic substance parallel to phytoplankton population. This is due to fact that approximately each individual in phytoplankton population is increasingly consuming the toxic substances. However, the reduction of grazing by zooplankton due toxicity effect is represented by the term m 2 z 2 (t). Furthermore, the toxicity effect on zooplankton population is less than phytoplankton population, where m 1 and m 2 are the toxicity coefficients with 0 < m 2 < m 1 [25].
Obviously, it is appropriate to explore the dynamics of any biological model by difference equations instead of differential equations when we are dealing with nonoverlapping populations. Furthermore, observation and analysis of chaos in any biological system by using difference equations is better than by using differential equations [26]. Hence it is interesting to study biological models in discrete form. Recently, Ghanbari and Gómez-Aguilar [27] discussed the dynamics of nutrient-phytoplankton-zooplankton system with variable-order fractional derivatives. Moreover, the authors in [28] explored the existence of chaos in a cancer model using fractional derivatives by means of exponential decay and the Mittag-Leffler law. Beigi et al. [29] discussed the use of reinforcement learning for effective vaccination strategies of coronavirus disease 2019 (COVID- 19). The authors in [30] analyzed the role of zooplankton dynamics for Southern Ocean phytoplankton biomass and global biogeochemical cycles. For more detail on the analysis of various dynamical systems, we refer the interested reader to [31][32][33][34][35][36][37]. There are various mathematical techniques for converting the systems of differential equations to their corresponding discrete counterparts. To achieve this goal, the usual way is applying standard difference schemes such as Runge-Kutta methods and Euler approximations. However, numerical inconsistency is experienced with the application of usual finite difference methods. Hence, to avoid this numerical inconsistency, we can apply the nonstandard finite difference method given by Mickens [38].
In general, whenever a nonstandard finite difference scheme is proposed, it is aimed on the preservation of the following properties of the respective continuous-time system: positivity of results, boundedness, stability of equilibrium points, and bifurcations. Moreover, the formation of these type of difference schemes is not straightforward, and there are no usual ways for their construction, which is probably considered as major downside of nonstandard difference schemes. Hence by taking into account the original dynamical properties of model (1.3) a discrete-time model from (1.3) is obtained by using Mickenstype nonstandard scheme such that it remains dynamically consistent [39]. Implementing the Mickens-type nonstandard scheme on model (1.3), we get the following discrete-time mathematical model: where h > 0 is taken as a step size for the nonstandard scheme. Furthermore, (1.4) can be written into the following mathematical form: where β > ρ. Moreover, our model (1.5) loses its biological consistency whenever β < ρ (see [25]), which is impossible biologically. Hence, for the rest of our paper, we assume that a > k and β > ρ.

Boundedness and existence of fixed points for system (1.5)
To obtain steady states of system (1.5), we consider the following two-dimensional system of equations: , 0) remains positive for r > Eq 1 . The existence and uniqueness of (p, z) can be studied as follows. Suppose that p 0 > 0 and z 0 > 0. Then each solution (p n , z n ) of system (1.5) must satisfy p n > 0 and z n > 0 for all n ≥ 0. Then from the first equation of system (1.5) it follows that Consequently, solving (2.2) and then taking the limit, we get In the same way, from second equation of system (1.5) we get Hence, we can obtain the upper bound for zooplankton population: Finally, we have the following theorem about the boundedness of all solutions of (1.5).
Next, we consider the equation system .
(2.5) From (2.5) we get the following pair: From this pair we can write Furthermore, at the upper bound and for each λ Hence F(p) = 0 has at least one positive real root in [0, k]. Furthermore, we can see that In addition, for λ = 0,

Stability analysis of system (1.5) about its fixed points
To discuss the stability of system (1.5) about all its equilibrium points, we compute the variational matrix V (p,z) of system (1.5) about each of its fixed point (p, z). The matrix V (p,z) is given by where Tr = (j 11 + j 22 ) and Dt = j 11 j 22j 12 j 21 .
The next lemma describes the conditions parallel to the Jurry condition for the stability of fixed points; see [40]. If ξ 1 and ξ 2 are characteristic values of (3.1), then the point (p, z) is sink if |ξ 1 | < 1 and |ξ 2 | < 1. Furthermore, it is locally asymptotically stable. The point (p, z) is known as a source (repeller) if |ξ 1 | > 1 and |ξ 2 | > 1, and it provides instability condition for the given system.

Proposition 3.2
Let ξ 1 and ξ 2 be the roots of the characteristic equation of the matrix V (0,0) and suppose that |ξ 2 | < 1 for all parametric values. Let (0, 0) be a population free fixed point of system (1.5). Then (0, 0) is sink or saddle if and only if r < Eq 1 or r > Eq 1 , respectively.
Next, we will explore the local stability of system (1.5) about the zooplankton-free equilibrium ( √ 4k 2 mr+r 2 -4Ek 2 mq 1 -r 2km , 0). Clearly, the first component in the pair , 0) be the variational matrix of the two-dimensional system (1.5) about the zooplankton free equilibrium ( , 0) has the following form: . Moreover, V 1 (x, 0) has the characteristic polynomial Hence we have the following proposition about the local stability of system (1.5) about the zooplankton-free equilibrium (

Proposition 3.3
Let ξ 1 and ξ 2 be the characteristic roots of (3.2), and let r > q 1 E. If

saddle point if and only if one of the following pairs of inequalities
is nonhyperbolic if and only if one of the following conditions is satisfied: Finally, it remains to analyze the local stability of system (1.5) about the only positive fixed point (p, z). Moreover, all parametric conditions for the existence of nonextinction fixed point (p, z) are given in Theorems 2.1 and 2.2. We can calculate the Jacobian matrix V 2 (p, z) of system (1.5) about (p, z) as follows: Let M(ξ ) be the characteristic polynomial for the matrix V 2 (p, z) with Then we have Taking into account the work done in the previous section, by Theorem 2.2 it follows that M(1) = h 2 pz(m 2 (a(βρ)(Ekq 1 + p(2r + km 1 (2a + 3p))) + 2kασ z + kαm 2 z 2 )) Hence the local stability of system (1.5) about (p, z) can be studied with the help of the following proposition.   k)r + Ekq 1 + p(2r + km 1 (2a + 3p))) + akη(βρ)m 2 (a + p)z pz(kασ 2m 2 (aρ((ak)r + Ekq 1 + p(2r + km 1 (2a + 3p))) -2kασ zkαm 2 z 2 )) ,

8)
and hp (ak)r + Ekq 1 + p 2r + km 1 (2a + 3p) (φ + hρm 2 z) = (φ + hρm 2 z)kη(a + p). This section is related to the bifurcation analysis of system (1.5) about (p, z), where all conditions for the existence and positivity of (p, z) are given in Theorems 2.1 and 2.2. Here we will discuss the Neimark-Scaker bifurcation experienced by system (1.5) about (p, z) under some conditions. Bifurcation is the mathematical phenomenon produced in any system due to creation of very small change in stability of the system. Furthermore, it causes some surprising changes in the dynamical standards of any mathematical system. Mathematically, bifurcation arises whenever parameters are varied in a very small neighborhood of an equilibrium point. Moreover, for further study of bifurcation theory and understanding this surprising behavior of a discrete-time mathematical system, we refer to [41][42][43][44][45][46].
Here we use standard theory of bifurcation for the study of Neimark-Sacker bifurcation of system (1.5) at (p, z). Let ξ 1 and ξ 2 be the roots of (3.1). Then both roots are complex with modulus one if (p, z) is a nonhyperbolic fixed point under condition (d) of Proposition 3.4.

Modified hybrid control strategy for controlling bifurcation and chaos
Generally, discrete-time systems are more complex to analyze as compared to a continuous-time one. For survival of life in any environment, it is necessary that the population does not experience any irregular situation. Hence, for controlling accidental uneven and unstable behavior in any mathematical system, chaos control is considered to be an applied tool for evading this complex and chaotic behavior [51][52][53]. In this part of the paper, we study a feedback control method with parameter perturbation to move unstable and irregular trajectories toward the stable trajectories. The most useful and well-known method in the field of chaos is given by Ott et al. [51] to control period-doubling bifurcation, which is known as OGY method. Later on, numerous strategic control methods are developed (see [53]). Here we consider a modified hybrid control method to control the Neimark-Sacker bifurcation and chaos. Furthermore, this mathematical method is well applicable to every discrete-time system experiencing the period-doubling bifurcation and chaos. Originally, a hybrid method was proposed by Liu et al. [52]. Moreover, it was developed to control the period-doubling bifurcation (see [54,55]). Here we reformed the existing hybrid control technique [52] to control the Neimark-Sacker bifurcation and chaos. Furthermore, the newly developed technique has shown better results for almost every discrete dynamical system. Consider the following n-dimensional discrete dynamical system: with Z n ∈ n , n ∈ Z, and the parameter μ ∈ for which system (5.1) experiences the bifurcation. The purpose of proposing the reformed method for controlling the bifurcation is recapturing the extreme range of stable region in (5.1) by lessening the length of unstable region. Hence we present the following generalized hybrid control method by applying state feedback along with parameter perturbation; where ∈ Z, and 0 < L < 1 is a parameter for controlling the bifurcation appearing in (5.2). In addition, g ( ) is the kth value of g(·). By application of (5.2) to system (1.5) we get the following system: The following theorem describes a necessary and sufficient condition for local stability of system (5.3) about (p, z).

Theorem 5.1 The positive constant solution (p, z) of system (5.3) is locally asymptotically stable if and only if
where Tr and Dt are the trace and determinant of (5.4), respectively.
For the understanding of limitation of modified hybrid control technique, we have the following remark.
Remark Like the hybrid method [52], the modified hybrid method (5.2) is feasible and efficient for those discrete-time mathematical models for which the stepsize parameter is taken as a bifurcation parameter.

Numerical simulation
In this section, we numerically study the dynamics of (1.5). This study is a direct verification of our theoretical analysis and analytic results we have proved in the previous sections. Particularly, in this section, we study the existence and direction of Neimark-Sacker bifurcation by using numeric values of the parameters. In addition, in this section, we take the initial conditions in the least neighborhood of the equilibrium point (p, z) for each case study. Additionally, we have Furthermore, for λ = 1.3997, we have (βρ) = 88.34599 and (δ + m 2 f (λ) + q 2 E) = 4.465067496648613. Then (βρ) > (δ + m 2 f (λ) + q 2 E), and we get where f (λ) = -(a + λ)(m 1 r 2 + q 1 E) α = -79.36413532682512 < 0.
Example 7.5 In this example, we compare the generalized hybrid method and hybrid method [52]. From Examples 7.3 and 7.4 we consider two discrete-time models  Table 1 we can observe that |I 1 | > |I 2 | for each variation of the parameter h ∈]0, 1[, where I 1 and I 2 are the controlled intervals corresponding to the controlled systems (7.3) and (7.4), respectively. Hence we can see from Table 1 that the generalized hybrid method (5.2) is much better than the old hybrid method [52].  In addition, our numerical study showed that the generalized hybrid method (5.2) is better than the hybrid method [52]. In addition, it is based on feedback control and it has brought back the stability of system (1.5) for large ranges of parameters. Moreover, from the numerical study we can see that the generalized hybrid method (5.2) is suitable for controlling the Neimark-Sacker bifurcation. At last, in Table 1, we provide a comparison of the generalized hybrid method (5.2) and hybrid method [52]. Moreover, from Fig. 4, Fig. 5, and Table 1 we can see that the generalized hybrid control technique restores the stability of system (1.5) for the maximal range of the control parameter.