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

Khan et al. Advances in Difference Equations (2021) 2021:443 https://doi.org/10.1186/s13662-021-03599-z RESEARCH Open Access Qualitative analysis of a discrete-time phytoplankton–zooplankton model with Holling type-II response and toxicity Muhammad Salman Khan1, Maria Samreen1* , Hassen Aydi2and Manuel De la Sen3 *Correspondence: [email protected]; [email protected] 1Department of Mathematics, Quaid-I-Azam University, 45320, Islamabad, Pakistan Full list of author information is available at the end of the article Abstract 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. Keywords: Phytoplankton–zooplankton model; Boundedness; Local stability analysis; Neimark–Sacker bifurcation; Generalized hybrid control method 1 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 ©The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Khan et al. Advances in Difference Equations (2021) 2021:443 Page 2 of 29 with physical procedure for the spatial structure. Natural processes include, for instance, development,grazing,andbehavior,andphysicalproceduresinclude,forinstance,mixing andlateralstirring.Nonlinearityofecosystemsentirelycontributestothespatialorganization 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 transportingCO2fromwatersurfacetothedepthofoceans.Mainly,thisprocesshappens 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 analgal 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 emergenceand disappearingofbloomsisnot clear,theundesirableeffectofdamagingalgal 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 toavert predationsby 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. Huppertetal. [10]consideredanutrient–phytoplanktonmodeltoexaminethedynamicalbehaviorofphytoplanktonblooms.In[11]theauthorshavepresentedazooplankton– phytoplankton model with harvesting. Furthermore, they have explained that the extra exploitationmayexterminatethepopulationwhilesuitableharvestingguarantiestheconsolidation 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–16].Itis convenienttointroducethetoxin creatinglag duringthe study of the dynamics of phytoplankton–zooplankton models. The authors in [17]havepresentedamathematical model including time lag in toxin deliverance by phytoplankton. The work done in [18–21] motivated us to study the dynamics of a phytoplankton–zooplankton popula- Khan et al. Advances in Difference Equations (2021) 2021:443 Page 3 of 29 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. •Wesupposethatz(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 ris their exponential rate of growth, and kis the maximum carrying capacity of environment. Under these conditions we have the following phytoplankton–zooplankton model [22]: ⎧ ⎨ ⎩ dp dt =rp(t)(1– p(t) k)–αf(p(t))z(t), dz dt =βf(p(t))z(t)–δz(t)–ρg(p(t))z(t). (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–preymodels usingaHolling-typeresponse function isbetterthan thestudyofdynamicsofpredator–preymodelswithoutusingHollingresponse.Generally, Hollingtype-IIresponseismodeledanddescribedbyusingrectangularhyperbola,andits mathematical form is given as ϕ(x)= x a+x, where ais any constant. By using Holling type-II response we get the following mathematical form of system (1.1): ⎧ ⎨ ⎩ dp dt =rp(t)(1– p(t) k)–αp(t) a+p(t)z(t), dz dt =βp(t) a+p(t)z(t)–δz(t)–ρp(t) a+p(t)z(t). (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 q1and q2for 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 Eas the parameter for combined effort for harvesting of population [25]. Khan et al. Advances in Difference Equations (2021) 2021:443 Page 4 of 29 Under these modifications, system (1.2) takes the following mathematical form: ⎧ ⎨ ⎩ dp dt =rp(t)(1– p(t) k)–αp(t) a+p(t)z(t)–m1p3(t)–q1Ep(t), dz dt =βp(t) a+p(t)z(t)–δz(t)–ρp(t) a+p(t)z(t)–m2z2(t)–q2Ez(t), (1.3) 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,thetermm1p3(t)appearinginsystem(1.3)representstheinfectionproduced in phytoplankton population due to an external toxic substance. In addition, d2 dp2(m1p3)= 6m1p>0showsanacceleratinggrowthoftoxicsubstanceparalleltophytoplanktonpopulation.Thisis duetofactthatapproximatelyeach individual in phytoplanktonpopulation isincreasinglyconsumingthetoxicsubstances.However,thereductionofgrazingbyzooplankton due toxicity effect is represented by the term m2z2(t). Furthermore, the toxicity effect on zooplankton population is less than phytoplankton population, where m1and m2are the toxicity coefficients with 0 <m2<m1[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]. Henceitisinterestingtostudybiologicalmodelsindiscreteform.Recently,Ghanbariand 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 (COVID19). The authors in [30] analyzed the role of zooplankton dynamics for Southern Ocean phytoplanktonbiomassandglobalbiogeochemicalcycles.Formoredetailontheanalysis of various dynamical systems, we refer the interested reader to [31–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,numericalinconsistencyisexperiencedwiththeapplicationofusualfinitedifference 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 onthepreservationof thefollowing propertiesof therespectivecontinuous-time system: positivityofresults,boundedness,stabilityofequilibriumpoints,andbifurcations.Moreover, the formation of these type of difference schemes is not straightforward, and there arenousualwaysfortheirconstruction,whichisprobablyconsidered asmajordownside of nonstandard difference schemes. Hence by taking into account the original dynamical properties of model (1.3) a discrete-time model from (1.3)isobtainedbyusingMickenstypenonstandardschemesuchthatitremainsdynamicallyconsistent[39].Implementing Khan et al. Advances in Difference Equations (2021) 2021:443 Page 5 of 29 theMickens-typenonstandardschemeonmodel(1.3),wegetthefollowingdiscrete-time mathematical model: ⎧ ⎨ ⎩ pn+1–pn h=rpn(1– pn+1 k)–αpn+1 a+pnzn–m1p2 npn+1 –q1Epn+1, zn+1–zn h=βpn a+pnzn–δzn+1 –ρpn a+pnzn+1 –m2znzn+1 –q2Ezn+1,(1.4) where h>0 is taken as a step size for the nonstandard scheme. Furthermore, (1.4)canbe written into the following mathematical form: ⎧ ⎪ ⎨ ⎪ ⎩ pn+1 =(1+hr)pn 1+h(r kpn+αzn a+pn+m1p2 n+q1E), zn+1 =(1+hβpn a+pn)zn 1+h(ρpn a+pn+δ+m2zn+q2E),(1.5) 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>kand β>ρ. 2 Boundedness and existence of fixed points for system (1.5) Toobtainsteadystatesofsystem (1.5),weconsiderthefollowingtwo-dimensionalsystem of equations: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ p=(1+hr)p 1+h(r kp+αz a+p+m1p2+q1E), z=(1+hβp a+p)z 1+h(ρp a+p+δ+m2z+q2E).(2.1) Solving (2.1), we can get the following equilibrium points: (0,0) which is an extinction pointforbothpopulations,(√4k2mr+r2–4Ek2mq1–r 2km ,0),whichisanextinctionequilibriumfor zooplankton population, and the unique positive equilibrium (p,z). Additionally, the first componentofthepoint(√4k2mr+r2–4Ek2mq1–r 2km ,0)remainspositiveforr>Eq1.Theexistence and uniqueness of (p,z) can be studied as follows. Suppose that p0>0andz0>0.Then eachsolution(pn,zn)ofsystem(1.5)mustsatisfypn>0andzn>0foralln≥0.Thenfrom the first equation of system (1.5) it follows that pn+1 ≤(1+hr)pn 1+hr kpn. (2.2) Consequently, solving (2.2) and then taking the limit, we get limsup n→∞ pn≤k. (2.3) In the same way, from second equation of system (1.5)weget zn+1 =(1+hβpn a+pn)zn 1+h(ρpn a+pn+δ+m2zn+q2E) ≤(1+hβk a+k)zn 1+h(ρk a+k+m2zn). Khan et al. Advances in Difference Equations (2021) 2021:443 Page 6 of 29 Hence, we can obtain the upper bound for zooplankton population: limsup n→∞ zn≤k(β–ρ) m2(k+a). (2.4) Finally, we have the following theorem about the boundedness of all solutions of (1.5). Theorem 2.1 Assume that 0<p0≤kand0<z0≤k(β–ρ) m2(k+a).Then for all n≥0, every positive solution (pn,zn)of system (1.5)is bounded and contained in the set [0,k]×[0, k(β–ρ) m2(k+a)] whenever β>ρ. Next, we consider the equation system ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ p=(1+hr)p 1+h(r kp+αz a+p+m1p2+q1E), z=(1+hβp a+p)z 1+h(ρp a+p+δ+m2z+q2E).(2.5) From (2.5)wegetthefollowingpair: p=a(β–ρ) (β–ρ)–(δ+m2z+q2E)–a,z=(a+p)(r–rp k–m1p2–q1E) α. From this pair we can write F(p)= a(β–ρ) (β–ρ)–(δ+m2f(p)+q2E)–a–p, (2.6) where f(p)=(a+p)(r–rp k–m1p2–q1E) α, (2.7) with f(0)= a(r–q1E) α>0 and F(0)= a(δ+m2f(0)+q2E) (β–ρ)–(δ+m2f(0)+q2E)>0. Furthermore, at the upper bound and for each λ∈(0,k], if (β–ρ)>(δ+m2f(λ)+q2E), then F(λ)= a(δ+m2f(λ)+q2E) (β–ρ)–(δ+m2f(λ)+q2E)–λ<0, where f(λ)=–(a+λ)(m1r2+q1E) α<0. Khan et al. Advances in Difference Equations (2021) 2021:443 Page 7 of 29 Hence F(p)=0 has at least one positive real root in [0,k]. Furthermore, we can see that F(λ)=–1+ a(β–ρ)(m2f(λ)) ((β–ρ)–(δ+m2f(λ)+q2E))2<0, where f(λ)=–(a+λ)(r k+2m1λ) α+r–rλ k–m1λ2–q1E α<0, whenever r–q1E<r k+m1λλ for every λ∈[0,k]. Hence the equation F(p) =0 has a unique positive solution in [0,k]. Theorem2.2 Assume that 0<p0≤kand0<z0≤k(β–ρ) m2(k+a).Then for r–q1E<r k+m1λλ and r>q1E, there exists a unique positive constant solution (p,z)of system (1.5)in [0,k]×[0, k(β–ρ) m2(k+a)]if and only if for each λ∈(0,k], we have (β–ρ)>δ+m2f(λ)+q2E. In addition,for λ=0, (β–ρ)<δ+m2f(λ)+q2E. 3 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 variationalmatrix V(p,z)ofsystem(1.5)abouteachofitsfixedpoint(p,z).ThematrixV(p,z) is given by V(p,z)=j11 j12 j21 j22. The characteristic polynomial M(ξ)ofV(p,z)is M(ξ)=ξ2–Tr ξ+Dt, (3.1) where Tr =(j11 +j22) Khan et al. Advances in Difference Equations (2021) 2021:443 Page 8 of 29 and Dt =j11j22 –j12j21. ThenextlemmadescribestheconditionsparalleltotheJurryconditionforthestabilityof fixed points; see [40]. Lemma3.1([40]) LetM(ξ)=ξ2–Tr ξ+DtandM(1) >0.Ifξ1,ξ2aretherootsofM(ξ)=0, then: (a)|ξ1|<1and |ξ2|<1if and only if M(–1)>0 and Dt <1; (b)|ξ1|>1and |ξ2|>1if and only if M(–1) >0 and Dt >1; (c)|ξ1|<1and |ξ2|>1or(|ξ1|>1and |ξ2<|1) if and only if M(–1)<0; (d)ξ1and ξ2represent complex conjugates with |ξ1|=1=|ξ2|if and only if Tr2–4Dt <0 and Dt =1. If ξ1and ξ2are characteristic values of (3.1), then the point (p,z)is sink if |ξ1|<1and |ξ2|<1.Furthermore,itislocallyasymptoticallystable.Thepoint(p,z)isknownasasource (repeller)if |ξ1|>1and |ξ2|>1,and it provides instability condition for the given system. The point (p,z)is a saddle point if |ξ1|<1and |ξ2|>1or (|ξ1|>1and |ξ2|<1).Finally, (p,z)is nonhyperbolic if condition (d)is satisfied. Firstly,wewillstudythestabilityofsystem(1.5)aboutpopulationfreeequilibriumpoint (0,0). The variational matrix V(p,z)for system (1.5) evaluated at (0,0) is V(0,0) =1+hr 1+Ehq10 01 1+hδ+Ehq2. Furthermore, V(0,0) isadiagonalmatrix.Hencesystem(1.5)hastwoeigenvaluesrelatedto thepopulationfreeequilibriumpoint(0,0),ξ1=1+hr 1+Ehq1andξ2=1 1+hδ+Ehq2,where,ξ1andξ2 arerootsofthecharacteristicequationofthematrixV(0,0).Itisclearthat|ξ2|=|1 1+hδ+Ehq2|< 1 for all parametric values. Now by considering the condition |ξ2|< 1 we are now able to describe stability conditions for system (1.5)about(0,0). Proposition3.2 Letξ1andξ2betherootsofthecharacteristicequationofthematrixV(0,0) andsupposethat|ξ2|<1forallparametricvalues.Let(0,0)beapopulationfreefixedpoint of system (1.5). Then (0,0) is sink or saddle if and only if r <Eq1or r >Eq1,respectively. Next, we will explore the local stability of system (1.5) about the zooplankton-free equilibrium (√4k2mr+r2–4Ek2mq1–r 2km ,0). Clearly, the first component in the pair (√4k2mr+r2–4Ek2mq1–r 2km ,0) is positive if and only if r>q1E.LetV1(√4k2mr+r2–4Ek2mq1–r 2km ,0) be the variational matrix of the two-dimensional system (1.5) about the zooplankton freeequilibrium(√4k2mr+r2–4Ek2mq1–r 2km ,0).ThenV1(√4k2mr+r2–4Ek2mq1–r 2km ,0)hasthefollowing form: V1(x,0)=k2(1+hr)(1+Ehq1–hm1x2) (k+Ehkq1+hx(r+km1x))2–hk2(1+hr)αx (a+x)(k+Ehkq1+hx(r+km1x))2 0a+x+hβx a+ahδ+x+h(δ+ρ)x+Ehq2(a+x), Khan et al. Advances in Difference Equations (2021) 2021:443 Page 9 of 29 where x=√4k2mr+r2–4Ek2mq1–r 2km .Moreover,V1(x,0) has the characteristic polynomial M(ξ)=ξ2–TrV1(x,0)+DetV1(x,0)(3.2) with TrV1(x,0)=a+x+hβx a+ahδ+x+h(δ+ρ)x+Ehq2(a+x) +k2(1+hr)(1+Ehq1–hm1x2) (k+Ehkq1+hx(r+km1x))2 and DetV1(x,0) =k2(1+hr)(a+x+hβx)(1+Ehq1–hm1x2) (a+ahδ+x+h(δ+ρ)x+Ehq2(a+x))(k+Ehkq1+hx(r+km1x))2. Hencewehavethefollowingpropositionaboutthelocalstabilityofsystem(1.5)aboutthe zooplankton-free equilibrium (√4k2m1r+r2–4Ek2m1q1–r 2km1,0). Proposition 3.3 Let ξ1and ξ2be the characteristic roots of (3.2), and let r >q1E.If (√4k2m1r+r2–4Ek2m1q1–r 2km1,0)=(x,0)is a zooplankton-free constant solution of (1.5), then: (a)(x,0)remains inside the unit disk if and only if 1+Ehq1–hm1x2<(k+Ehkq1+hx(r+km1x))2 k2(1+hr)(3.3) and βx<aδ+(δ+ρ)x+Eq2(a+x). (3.4) (b)(x,0)lies outside the unit disk if and only if 1+Ehq1–hm1x2>(k+Ehkq1+hx(r+km1x))2 k2(1+hr)(3.5) and βx>aδ+(δ+ρ)x+Eq2(a+x). (3.6) (c)(x,0)isasaddlepointifandonlyifoneofthefollowingpairsofinequalities(3.4)–(3.5) or (3.3)–(3.6)is satisfied. (d)(x,0)is nonhyperbolic if and only if one of the following conditions is satisfied: 1+Ehq1–hm1x2=(k+Ehkq1+hx(r+km1x))2 k2(1+hr) or aδ+(δ+ρ–β)x+Eq2(a+x)=0. Khan et al. Advances in Difference Equations (2021) 2021:443 Page 16 of 29 Finally, to write the linear part of (4.3) in the canonical matrix form at ˘ h=0, we consider the similarity transformation P Z=v12 0 –v11 –℘X Y,, (4.6) where =Tr(0) 2, and ℘=4Dt(0)–Tr2(0) 2. Then from (4.6)wehave X Y=1 v12 0 –v11 ℘v12 –1 ℘P Z.. (4.7) By using transformation (4.6)wehavethenextauthoritativeformofsystem(4.3): X Y→–℘ ℘ X Y+˘ F(X,Y) ˘ G(X,Y), (4.8) where ˘ F(X,Y)=v16P3 v12 +v17P2Z v12 +v13P2 v12 +v18PZ2 v12 +v14PZ v12 +Z3v19 v12 +Z2v15 v12 +O|X|+|Y|4 and ˘ G(X,Y)=(–v11)v16 v12℘–v26 ℘P3+(–v11)v17 v12℘–v27 ℘P2Z +(–v11)v13 v12℘–v23 ℘P2+(–v11)v18 v12℘–v28 ℘PZ2 +(–v11)v14 v12℘–v24 ℘PZ +(–v11)v19 v12℘–v29 ℘Z3 +(–v11)v15 v12℘–v25 ℘Z2+O|X|+|Y|4 withP=v12XandZ=(–v11)X–℘Y.Hencebyusingthestandardtheoryofnormalform for analysis of bifurcation we can calculate the first Lyapunov exponent at (X,Y)=(0,0) as follows: =–Re(1–2ξ1)ξ2 2 1–ξ1θ20θ11–1 2|θ11|2–|θ02|2+Re(ξ2θ21)˘ h=0, Khan et al. Advances in Difference Equations (2021) 2021:443 Page 17 of 29 where θ20 =1 8˘ FXX –˘ FYY +2˘ GXY +i(˘ GXX –˘ GYY –2˘ FXY), θ11 =1 4˘ FXX +˘ FYY +i(˘ GXX +˘ GYY), θ02 =1 8˘ FXX –˘ FYY –2˘ GXY +i(˘ GXX –˘ GYY +2˘ FXY), θ21 =1 16(˘ FXXX +˘ FXYY +˘ GXXY +˘ GYYY) +i 16(˘ GXXX +˘ GXYY –˘ FXXY –˘ FYYY). Due to aforementioned analysis, we have the following theorem (see [45–50]). Theorem4.1 Assume that (3.7), (3.8), and (3.9)are satisfied and =0.Then the unique positive fixed point (p,z)of system (1.5)undergoes Neimark–Sacker bifurcation.Additionally,if <0,then for h >˘ h,an attracting invariant closed curve bifurcates from the fixed point (p,z), and if >0,then for h <˘ h,a repelling invariant closed curve bifurcates from the fixed point (p,z). 5 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 unevenandunstablebehaviorinanymathematicalsystem,chaoscontrolisconsideredtobe an applied tool for evading this complex and chaotic behavior [51–53]. In this part of the paper,westudyafeedbackcontrolmethodwithparameterperturbationtomoveunstable 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 wasdevelopedtocontroltheperiod-doublingbifurcation(see[54,55]).Herewereformed theexistinghybrid controltechnique[52]tocontroltheNeimark–Sackerbifurcationand chaos. Furthermore, the newly developed technique has shown better results for almost every discrete dynamical system. Consider the following n-dimensional discrete dynamical system: Zn+1 =g(Zn,μ) (5.1) with Zn∈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 unstableregion.Hencewepresentthefollowinggeneralizedhybridcontrolmethodbyapplying Khan et al. Advances in Difference Equations (2021) 2021:443 Page 18 of 29 state feedback along with parameter perturbation; Zn+k=L3g()(Zn,μ)+1–L3Zn, (5.2) where∈Z,and0<L<1isaparameterforcontrollingthebifurcationappearingin(5.2). In addition, g()is the kth value of g(·). By application of (5.2)tosystem(1.5)wegetthe following system: ⎧ ⎪ ⎨ ⎪ ⎩ pn+1 =L3((1+hr)pn 1+h(r kpn+αzn a+pn+m1p2 n+q1E))+(1–L3)pn, zn+1 =L3((1+hβpn a+pn)zn 1+h(ρpn a+pn+δ+m2zn+q2E))+(1–L3)zn.(5.3) Furthermore, systems (5.3)and(1.5) have the same constant solutions. Additionally, the Jacobian matrix of (5.3)about(p,z) is given as follows: ⎛ ⎝1–hL3p((a–k)r+ekq1+p(2r+km1(2a+3p))) k(1+hr)(a+p)–hL3αp (1+hr)(a+p) hL3z(–β+δ+ρ+eq2+m2z)2 a(β+hβδ–ρ+ehβq2+hβm2z) β+hβδ–ρ+ehβq2+h(β–L3β+L3ρ)m2z β+hβδ–ρ+ehβq2+hβm2z⎞ ⎠. (5.4) The following theorem describes a necessary and sufficient condition for local stability of system (5.3)about(p,z). Theorem5.1 The positive constant solution (p,z)of system (5.3)is locally asymptotically stable if and only if |Tr |<1+Dt <2, 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)isfeasibleand efficient for those discrete-time mathematical models for which the stepsize parameter is taken as a bifurcation parameter. 6 Hybrid control of Neimark–Sacker bifurcation Inthissection,weapplythehybridtechnique[52]tosystem(1.5)tocontroltheNeimark– Sackerbifurcation.Moreover,thismethodisusedascontrolstrategybymanyresearchers for controlling the period-doubling bifurcation, Neimark–Sacker bifurcation, and chaos under the effects of period-doubling bifurcation (see [54,55]). By application of a hybrid method [52]tosystem(1.5)wegetthefollowingsystem: ⎧ ⎪ ⎨ ⎪ ⎩ pn+1 =S1((1+hr)pn 1+h(r kpn+αzn a+pn+m1p2 n+q1E))+(1–S1)pn, zn+1 =S1((1+hβpn a+pn)zn 1+h(ρpn a+pn+δ+m2zn+q2E))+(1–S1)zn,(6.1) Khan et al. Advances in Difference Equations (2021) 2021:443 Page 19 of 29 where0<S1<1isacontrolparameter.Furthermore,systems(6.1)and(1.5)havethesame constant solutions. Additionally, the Jacobian matrix of (6.1)about(p,z)is ⎛ ⎝1–hS1p((a–k)r+ekq1+p(2r+km1(2a+3p))) k(1+hr)(a+p)–hS1αp (1+hr)(a+p) hS1z(–β+δ+ρ+eq2+m2z)2 a(β+hβδ–ρ+ehβq2+hβm2z) β+hβδ–ρ+ehβq2+h(β–S1β+S1ρ)m2z β+hβδ–ρ+ehβq2+hβm2z⎞ ⎠. (6.2) 7 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. Example 7.1 Let a= 2.0099,q1= 0.0189,q2= 1.2994,r= 10.5923,E= 0.9959,k= 1.3997, β= 98.499,α= 2.9999,δ= 0.0384,ρ= 10.5842,m1= 0.6222,m2= 0.4422,p0= 0.105348, z0= 6.8884251, and h∈(0,1)]. In this case the extinction equilibrium and nonextinction equilibrium for zooplankton population are (x,0) = (1.26553,0) and (p,z)= (0.1053484,6.888425), respectively. Then from system (1.5)wehave limsup n→∞ pn≤k=1.3997. Then by using the value k= 1.3997 in the second equation of system (1.5)weget limsup n→∞ zn≤k(β–ρ) m2(k+a)=81.61590194902372. Hence we have (0.1053484,6.888425) ∈[0,k]×[0, k(β–ρ) m2(k+a)]forβ>ρ, which verifies Theorem 2.1. Additionally, we have f(0)= a(r–q1E) α=7.084113606170539>0 and F(0)= a(δ+m2f(0)+q2E) (β–ρ)–(δ+m2f(0)+q2E)=0.10754185654404144 >0. Furthermore, for λ= 1.3997, we have (β–ρ) = 88.34599 and (δ+m2f(λ)+q2E)= 4.465067496648613. Then (β–ρ)>(δ+m2f(λ)+q2E), and we get F(λ)= a(δ+m2f(λ)+q2E) (β–ρ)–(δ+m2f(λ)+q2E)–λ=–1.2921581434559586<0, where f(λ)=–(a+λ)(m1r2+q1E) α=–79.36413532682512<0. Khan et al. Advances in Difference Equations (2021) 2021:443 Page 20 of 29 Moreover, we have F(λ)=–1+ a(β–ρ)(m2f(λ)) ((β–ρ)–(δ+m2f(λ)+q2E))2=–1.0580185510185083<0, where f(λ)=–(a+λ)(r k+2m1λ) α+r–rλ k–m1λ2–q1E α=–10.993342080620321<0, which verifies Theorem 2.2. Example 7.2 Let a= 2.0099,q1= 0.0189,q2= 1.2994,r= 10.5923,E= 0.9959,k= 1.3997, β= 98.499,α= 2.9999,δ= 0.0384,ρ= 10.5842,m1= 0.6222,m2= 0.4422,p0= 0.105348, z0=6.8884251, and h∈(0,1)]. Then system (1.5) takes the form ⎧ ⎪ ⎨ ⎪ ⎩ pn+1 =(1+10.5923h)pn 1+h(10.5923 1.3997 pn+2.9999zn 2.0099+pn+0.6222p2 n+0.0188), zn+1 =(1+h98.499pn 2.0099+pn)zn 1+h(10.5842pn 2.0099+pn+0.0384+0.4422zn+1.2941).(7.1) Additionally, in this case the extinction equilibrium and nonextinction equilibrium for zooplanktonpopulationare(x,0)=(1.26553,0)and(p,z)=(0.1053484,6.888425),respectively. In this case the graphical behavior of both population variables is shown in Fig. 2. In addition, Fig. 2(c) represents the maximum Lyapunov exponent for system (7.1). In Fig. 3,somephaseportraitsaregiven,wherehvaries in ]0,1[. We can easily see thatthere exists the Neimark–Sacker bifurcation when hcertainly passes through h= 0.38022 (see Fig. 3(b)). For the aforementioned values of parameters, the Jacobian matrix V2(p,z)for system (7.1)is V2(0.1053484433,6.88842511) =0.9751639697174742 –0.01143564868241927 36.869495419357726 0.5871690414953413 . Moreover, the characteristic equation M(ξ)=0forV2(0.1053484433,6.88842511) is ξ2–1.8038899298174336ξ+1 =0. (7.2) Solving (7.2), we get ξ1= 0.7811665056064 + 0.6196705420078iand ξ2= 0.7811665056064–0.6196705420078iwith |ξ1|=|ξ1|=1. In addition, we have M(–1)=3.556545701326458 >0 and M(1)=0.43187967890082724 >0. Now from (4.3)wehave ˘ f(P,Z)=0.105348+0.974755P–0.0116237Z–0.256573P2–0.099269PZ Khan et al. Advances in Difference Equations (2021) 2021:443 Page 21 of 29 +0.001282Z2–0.142822P3+0.10288P2Z+0.010039PZ2 –0.000141506Z3+O|P|+|Z|4 and ˘ g(P,Z)=6.88842+0.574993Z+37.95688P–43.12450P2+3.450296PZ –0.0354763Z2+48.99565P3–2.366456P2Z–0.1893441PZ2 +0.00218884Z3+O|P|+|Z|4. Finally,whensystem(4.3)isconvertedintothecanonical form (4.6),weobtainthematrix 1 v12 0 –v11 ℘v12 –1 ℘=–0.01143564869 0 –0.1939974644 –0.6196705420078615 with 1 v12 0 –v11 ℘v12 –1 ℘–1 =–87.4458482512197577 0.0 27.3762776879404903 –1.61376075222132043. Furthermore, from (4.8)wehave ˘ F(X,Y)=22.073323320819P2+8.5402951060818PZ –0.1103358266003Z2 +12.2871550765P3–8.85102538857P2Z–0.8636988319698PZ2 +0.012173994624Z3+O|X|+|Y|4 and ˘ G(X,Y)=61.114564266701P2–8.141785570023PZ +0.0908219975961Z2 –81.22561954155P3+6.52880078485P2Z+0.571452207625PZ2 –0.0072969596695Z3+O|X|+|Y|4. Additionally,plotsfor ˘ F(X,Y)and˘ G(X,Y)withsolutionat(0,0)arepresentedinFigs.1(a) and 1(b), respectively, where P= (–0.1162370758)Xand Z= (–0.1998810610)X– (0.6334408575)Y. Finally, we get θ20 =1 8˘ FXX –˘ FYY +2˘ GXY +i(˘ GXX –˘ GYY –2˘ FXY)=0.005831–0.019066i, θ11 =1 4˘ FXX +˘ FYY +i(˘ GXX +˘ GYY)=–0.01195+0.013959i, θ02 =1 8˘ FXX –˘ FYY –2˘ GXY +i(˘ GXX –˘ GYY +2˘ FXY)=0.02389–0.00182i, θ21 =1 16(˘ FXXX +˘ FXYY +˘ GXXY +˘ GYYY) +i 16(˘ GXXX +˘ GXYY –˘ FXXY –˘ FYYY)= 0.00755+0.00574i, Khan et al. Advances in Difference Equations (2021) 2021:443 Page 22 of 29 Figure 1 Plots for ˘ F(X,Y)and˘ G(X,Y)fora= 2.0099, q1= 0.0189, q2= 1.2994, r= 10.5923, E= 0.9959, k= 1.3997, β= 98.499, α= 2.9999, δ= 0.0384, ρ= 10.5842, m1= 0.6222, m2= 0.4422, and h∈(0,1) and =–Re(1–2ξ1)ξ2 2 1–ξ1θ20θ11–1 2|θ11|2–|θ02|2+Re(ξ2θ21)ˆ h=0 =–0.00138464412<0. Hence the condition for the existence of Neimark–Sacker bifurcation is satisfied (see Theorem 4.1). Example 7.3 This example is related to the study of control of Neimark–Sacker bifurcation by using generalized hybrid technique (5.2). To show the effectiveness of generalized technique, we have used the same values of parameters as in Example 7.1. Consider the following system of difference equations: ⎧ ⎪ ⎨ ⎪ ⎩ pn+1 =L3((1+10.5923h)pn 1+h(10.5923 1.3997 pn+2.9999zn 2.0099+pn+0.6222p2 n+0.0188))+(1–L3)pn, zn+1 =L3((1+h98.499pn 2.0099+pn)zn 1+h(10.5842pn 2.0099+pn+0.0384+0.4422zn+1.2941))+(1–L3)zn,(7.3) where a=2.0099,q1= 0.0189,q2= 1.2994,r=10.5923,E=0.9959,k=1.3997,β= 98.499, α= 2.9999,δ= 0.0384,ρ= 10.5842,m1= 0.6222,m2= 0.4422,h= 0.699909. In addition, 0 < L< 1 is the control parameter. Furthermore, for system (7.3), we have (p,z)= (0.1053484433,6.88842511), which is a unique positive constant solution of the original system (1.5). Additionally, controlled diagrams for zooplankton and phytoplankton populations by using models (7.3)areshowninFigs.4(b) and 4(a), respectively. Finally, we can see that the stability of initial system (1.5) is victoriously regained for large range of control parameter by using the generalized hybrid control method (see Fig. 4). Khan et al. Advances in Difference Equations (2021) 2021:443 Page 23 of 29 Figure 2 Plots of system (1.5)fora= 2.0099, q1= 0.0189, q2= 1.2994, r= 10.5923, E= 0.9959, k= 1.3997, β= 98.499, α= 2.9999, δ= 0.0384, ρ= 10.5842, m1= 0.6222, m2= 0.4422, and h∈(0,1) Example7.4 ThisexampleisrelatedtothestudyofcontrolofNeimark–Sackerbifurcation by using a hybrid technique [52] of control. Consider the system of difference equations ⎧ ⎪ ⎨ ⎪ ⎩ pn+1 =S1((1+10.5923h)pn 1+h(10.5923 1.3997 pn+2.9999zn 2.0099+pn+0.6222p2 n+0.0188))+(1–S1)pn, zn+1 =S1((1+h98.499pn 2.0099+pn)zn 1+h(10.5842pn 2.0099+pn+0.0384+0.4422zn+1.2941))+(1–S1)zn,(7.4) where a=2.0099,q1= 0.0189,q2= 1.2994,r=10.5923,E=0.9959,k=1.3997,β= 98.499, α= 2.9999,δ= 0.0384,ρ= 10.5842,m1= 0.6222,m2= 0.4422,h= 0.699909. In addition, 0<S1< 1 is the control parameter. Furthermore, for system (7.4), we have (p,z)= (0.1053484433,6.88842511), which is a unique positive constant solution of the original system (1.5). Additionally, controlled diagrams for zooplankton and phytoplankton populations for system (7.4) are respectively shown in Figs. 5(b)and5(a). 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 (7.3)and(7.4), respectively. Moreover, in this case, we have taken L,S1∈]0,1[ and a= 2.0099,q1= 0.0189,q2= 1.2994,r= 10.5923,E= 0.9959,k= 1.3997,β= 98.499,α= 2.9999,δ=0.0384,ρ=10.5842,m1=0.6222,m2=0.4422, and h∈(0,1). Form both systems (7.3)and(7.4)weget(p,z)=(0.1053484433,6.88842511) as a nique positivefixedpoint.Additionally,fromTable1wecanobservethat|I1|>|I2|foreachvariationoftheparameterh∈]0,1[,whereI1andI2arethecontrolledintervalscorresponding Khan et al. Advances in Difference Equations (2021) 2021:443 Page 24 of 29 Figure 3 Phase portraits of system (1.5)fora= 2.0099, q1= 0.0189, q2= 1.2994, r= 10.5923, E= 0.9959, k= 1.3997, β= 98.499, α= 2.9999, δ= 0.0384, ρ= 10.5842, m1= 0.6222, m2= 0.4422, and h∈(0,1) to the controlled systems (7.3)and(7.4),respectively.Hencewecanseefrom Table1that the generalized hybrid method (5.2) is much better than the old hybrid method [52]. Example 7.6 In this example, we compare the dynamics of systems (1.3)and(1.5). Forcase(i),wetakea=2.1,q1= 0.09,q2=0.3,r=1.5,c= 0.14,k= 100,β=0.5,α= 0.69,δ= 0.001,ρ=0.1,m2= 0.021, and m1= 0.06. Then we get the fixed point (p,z)= (2.3059,7.38854),whichisauniquepositiveconstantsolutionof(1.3)and(1.5).Moreover, for the initial conditions p0= 2.3059 and z0= 7.38854, Figs. 6(a) and 7(b) are plotted for systems(1.5)and(1.3),respectively. Consequently, we canseethatsystems(1.5)and(1.3) arestableat(p,z) = (2.3059,7.38854) for m1=0.06(seeFigs.6(a) and 7(b)). In addition, Khan et al. Advances in Difference Equations (2021) 2021:443 Page 25 of 29 Figure 4 Controlled diagrams for system (7.1)fora= 2.0099, q1= 0.0189, q2= 1.2994, r= 10.5923, E= 0.9959, k= 1.3997, β= 98.499, α= 2.9999, δ= 0.0384, ρ= 10.5842, m1= 0.6222, m2= 0.4422, h= 0.699909 and L∈(0,1) Figure 5 Controlled diagrams for system (7.4)fora= 2.0099, q1= 0.0189, q2= 1.2994, r= 10.5923, E= 0.9959, k= 1.3997, β= 98.499, α= 2.9999, δ= 0.0384, ρ= 10.5842, m1= 0.6222, m2= 0.4422, h= 0.699909, and S1∈(0,1) Table 1 Comparison of the modified hybrid method (5.2) and hybrid method [52]forL,S1∈]0,1[ and a= 2.0099, q1= 0.0189, q2= 1.2994, r= 10.5923, E= 0.9959, k= 1.3997, β= 98.499, α= 2.9999, δ= 0.0384, ρ= 10.5842, m1= 0.6222, m2= 0.4422, h∈(0,1) h∈]0,1[ Controlled interval I1for (7.3) Controlled interval I2for (7.4) 0.48889569 0 < L< 0.99288325491279 0 < S1< 0.97880134847096 0.58889569 0 < L< 0.983282763587931 0 < S1< 0.95068201684448 0.68889569 0 < L< 0.97635403882826 0 < S1< 0.93072628972275 0.78889569 0 < L< 0.97111704728074 0 < S1< 0.91582972183557 0.88889569 0 < L< 0.96701917628990 0 < S1< 0.90428485868004 0.98889569 0 < L< 0.96372500134675 0 < S1< 0.89507489723917 for case (ii), we take m1= 0.025, a=2.1,q1= 0.09,q2=0.3,r=1.5,c= 0.14,k= 100,β= 0.5,α=0.69,δ=0.001,ρ=0.1,m2=0.021. We get the fixed point (p,z)=(2.8059,8.8854), whichis a unique positive constantsolutionof(1.3)and(1.5).Hencewecanseethatboth systems (1.3)and(1.5)areunstableat(p,z)=(2.8059,8.8854) (see Figs. 6(b)and 7(a)). Finally, Fig. 6(c), (d) shows the existence of Neimark–Sacker bifurcation in system (1.5)for lower values of stepsize h,andFig.7(c), (d) shows that both variables p(t)andz(t)from system (1.3)areunstableat(p,z)=(2.8059,8.8854) (see [25]).