
In this paper, the dynamics of a discrete-time chemostat model were investigated. The discretization was obtained using the piecewise constant argument method. An analysis was performed to determine the existence and stability of fixed points. In addition, we have shown that the model experiences transcritical and period-doubling bifurcations. Two chaos control techniques, feedback control and hybrid control, were employed to control bifurcation and chaos in the model. Moreover, we provided numerical simulations to substantiate our theoretical results. This study illustrates that the piecewise constant argument method is more dynamically consistent than the forward Euler method.
Citation: Ibraheem M. Alsulami. On the stability, chaos and bifurcation analysis of a discrete-time chemostat model using the piecewise constant argument method[J]. AIMS Mathematics, 2024, 9(12): 33861-33878. doi: 10.3934/math.20241615
[1] | Vinoth Seralan, R. Vadivel, Dimplekumar Chalishajar, Nallappan Gunasekaran . Dynamical complexities and chaos control in a Ricker type predator-prey model with additive Allee effect. AIMS Mathematics, 2023, 8(10): 22896-22923. doi: 10.3934/math.20231165 |
[2] | Limei Liu, Xitong Zhong . Analysis and anti-control of period-doubling bifurcation for the one-dimensional discrete system with three parameters and a square term. AIMS Mathematics, 2025, 10(2): 3227-3250. doi: 10.3934/math.2025150 |
[3] | Abdul Qadeer Khan, Ayesha Yaqoob, Ateq Alsaadi . Neimark-Sacker bifurcation, chaos, and local stability of a discrete Hepatitis C virus model. AIMS Mathematics, 2024, 9(11): 31985-32013. doi: 10.3934/math.20241537 |
[4] | Saud Fahad Aldosary, Rizwan Ahmed . Stability and bifurcation analysis of a discrete Leslie predator-prey system via piecewise constant argument method. AIMS Mathematics, 2024, 9(2): 4684-4706. doi: 10.3934/math.2024226 |
[5] | Yao Shi, Zhenyu Wang . Bifurcation analysis and chaos control of a discrete fractional-order Leslie-Gower model with fear factor. AIMS Mathematics, 2024, 9(11): 30298-30319. doi: 10.3934/math.20241462 |
[6] | Ali Al Khabyah, Rizwan Ahmed, Muhammad Saeed Akram, Shehraz Akhtar . Stability, bifurcation, and chaos control in a discrete predator-prey model with strong Allee effect. AIMS Mathematics, 2023, 8(4): 8060-8081. doi: 10.3934/math.2023408 |
[7] | Tao Xie, Wenqing Zheng . Robustness analysis of Cohen-Grossberg neural network with piecewise constant argument and stochastic disturbances. AIMS Mathematics, 2024, 9(2): 3097-3125. doi: 10.3934/math.2024151 |
[8] | Aeshah A. Raezah, Jahangir Chowdhury, Fahad Al Basir . Global stability of the interior equilibrium and the stability of Hopf bifurcating limit cycle in a model for crop pest control. AIMS Mathematics, 2024, 9(9): 24229-24246. doi: 10.3934/math.20241179 |
[9] | Ahmad Suleman, Rizwan Ahmed, Fehaid Salem Alshammari, Nehad Ali Shah . Dynamic complexity of a slow-fast predator-prey model with herd behavior. AIMS Mathematics, 2023, 8(10): 24446-24472. doi: 10.3934/math.20231247 |
[10] | A. Q. Khan, Ibraheem M. Alsulami, S. K. A. Hamdani . Controlling the chaos and bifurcations of a discrete prey-predator model. AIMS Mathematics, 2024, 9(1): 1783-1818. doi: 10.3934/math.2024087 |
In this paper, the dynamics of a discrete-time chemostat model were investigated. The discretization was obtained using the piecewise constant argument method. An analysis was performed to determine the existence and stability of fixed points. In addition, we have shown that the model experiences transcritical and period-doubling bifurcations. Two chaos control techniques, feedback control and hybrid control, were employed to control bifurcation and chaos in the model. Moreover, we provided numerical simulations to substantiate our theoretical results. This study illustrates that the piecewise constant argument method is more dynamically consistent than the forward Euler method.
Mathematical models are essential tools in the study of complicated events inside biological models. Out of all these models, the chemostat model is particularly notable for its role in explaining the mechanics of cell mass increase in controlled conditions. The chemostat is a continuous bioreactor that is crucial for developing microalgae and different microorganisms in both laboratory and industrial environments. The chemostat model is extensively used in several domains, including microbial ecology and bioprocess engineering. The adaptability and effectiveness of this model arise from its capacity to maintain a consistent environment of nutrient availability while simultaneously eliminating waste products. This allows for the ongoing development of microbial populations under carefully controlled environments.
Ecological theory often uses the chemostat model as a simplified depiction of actual ecosystems, such as lakes or oceans. Researchers obtain insights into the principles regulating population dynamics, resource competition, and ecosystem stability by imitating the dynamics of nutrient inflow and microbial development. This abstraction enables the investigation of ecological topics in a controlled laboratory environment, providing vital understanding of the dynamics of real ecosystems. Moreover, the chemostat model is a great tool for analyzing the effectiveness of antibiotics and developing methods for wastewater treatment. The chemostat is used in antibiotic research to provide a controlled environment for evaluating bacteria susceptibility, resistance mechanisms, and pharmacokinetic factors. Similarly, in the context of wastewater treatment, the chemostat allows researchers to investigate the capacity of microbial communities to break down pollutants and eliminate contaminants from effluent streams.
Consider the following continuous-time chemostat model [1]:
{dNdt=α(C1+C)N−N,dCdt=−(C1+C)N−C+β, | (1.1) |
where N(t) represents the number of bacterial cells per unit volume of the growth medium at time t, reflecting the abundance of bacteria within the chemostat. C(t) denotes the concentration of nutrients at time t available in the growth compartment, influencing bacterial growth rates and population dynamics. Moreover, α quantifies the intrinsic growth rate of bacterial populations, and β represents the rate at which fresh nutrients are supplied to the chemostat. The parameters α and β are positive constants. The detailed derivation of this model is given in [1]. They discretized this model using the forward Euler method and thus obtained the following discrete-time chemostat model:
{Nt+1=(1−h)Nt+αh(Ct1+Ct)Nt,Ct+1=βh+(1−h)Ct−h(Ct1+Ct)Nt, | (1.2) |
where h>0 is the step size. They investigated the existence and stability of fixed points (FPs). Moreover, it is proved that their model experiences only period-doubling (PD) bifurcation at the positive FP.
Dynamic models may be expressed as continuous-time models, which are characterized by differential equations, or as discrete-time models, which are described by difference equations. The study of discrete-time models is growing in popularity due to their effectiveness in non-overlapping generation, ease of numerical solution acquisition, and more complex dynamical behaviors. Recently, discrete-time models have received considerable interest from researchers [2,3,4,5,6,7]. There are several strategies for converting continuous models to discrete ones, including the forward Euler approach [8,9,10,11], the piecewise constant argument method [12,13,14,15], and the nonstandard finite difference scheme [16,17,18].
Akhmet [19] investigated differential equations characterized by a piecewise constant argument of generalized type (EPCAG) and proposed a novel deviation function form, specifying the solution sets and setting criteria for the stability of the zero solution. Alwan et al. [20] established a comparison principle for nonlinear differential equations with piecewise constant arguments (EPCA), employing Lyapunov functions to demonstrate stability characteristics, especially in scenarios where the arguments facilitate stabilization, with implications for logistic growth models. Khan [21] examined a discrete Kolmogorov model with a piecewise constant argument, analyzing local dynamics, chaos, and bifurcations, especially flip bifurcations, and illustrated the efficacy of feedback control in stabilizing chaos. Bozkurt et al. [22] developed a colorectal cancer model using piecewise constant arguments to investigate tumor development and the effects of chemo-immunotherapy, assessing equilibrium stability and demonstrating bifurcation phenomena, including period-doubling and Neimark-Sacker behaviors. Naik et al. [23] used the piecewise constant argument approach on a predator-prey model including the Allee effect, identifying Neimark-Sacker bifurcations and using feedback and hybrid controls to mitigate bifurcations and chaos, therefore emphasizing the Allee effect's significance in population dynamics.
A discrete-time model is dynamically consistent with its continuous counterpart when both demonstrate similar dynamical characteristics, such as boundedness, solution persistence, steady state stability, chaos, and bifurcation. While the forward Euler method is often employed for discretization, it is not dynamically consistent with its continuous counterparts. An important issue is that the discrete model generated by the Euler technique is not completely realistic. This is due to the fact that some initial values and parameters may produce negative values for predator and prey populations. Another issue is that the discretized model has an extra parameter h, the step size. This new parameter also affects the dynamics of the model. Nonetheless, use of the piecewise constant argument method prevents the chance of negative values. Moreover, we have the same parameters in the discretized model as in our original continuous-time model. Thus, we discretize the model (1.1) employing the piecewise constant argument technique. By utilizing piecewise constant arguments to solve nonlinear differential equations and considering the regular time interval for the average growth rate in both populations, we can rewrite system (1.1) as follows:
{1N(t)dNdt=α(C[t]1+C[t])−1,1C(t)dCdt=−(11+C[t])N[t]−1+βC[t], | (1.3) |
where [t] represents the integer part of t, and 0<t<∞. Furthermore, integrating system (1.3) on an interval [n,n+1) with n=0,1,2,⋯ yields the following system:
{Nt=Nnexp((αCn1+Cn−1)(t−n)),Ct=Cnexp((−Nn1+Cn−1+βCn)(t−n)). | (1.4) |
Taking t→n+1, we obtain the following discrete-time system:
{Nn+1=Nnexp(αCn1+Cn−1),Cn+1=Cnexp(−Nn1+Cn−1+βCn). | (1.5) |
After replacing n by t, we obtain
{Nt+1=Ntexp(αCt1+Ct−1),Ct+1=Ctexp(−Nt1+Ct−1+βCt). | (1.6) |
The model (1.2) can produce negative values for N(t) and C(t), particularly when h is quite big and α is considerably small. However, due to the properties of the exponential function, which has a positive range, the model (1.6) excludes the possibility of negative values for N(t) and C(t). This demonstrates that the piecewise constant argument method is more appropriate than the forward Euler method.
The discrete-time chemostat model (1.6) analyzed in this study is novel in its application of the piecewise constant argument method, which overcomes the limitations of the traditional forward Euler method by ensuring biologically realistic, non-negative values for populations and nutrients. The main objective of this study is to investigate the stability and bifurcation in a discrete-time chemostat model via the piecewise constant argument method. For the detailed analysis of stability and bifurcation in discrete-time models, we refer the readers to [24,25,26,27,28,29] and the references therein. The subsequent sections of the paper are arranged as follows: Section 2 is dedicated to investigating the existence and stability of FPs. Section 3 examines the study of bifurcations that include transcritical (TC) and period-doubling (PD) bifurcations. Two chaos control techniques, feedback control and hybrid control, are presented in Section 4. In Section 5, numerical simulation results are presented to support the theoretical analysis and display the new and rich dynamic behavior. Finally, a brief conclusion is presented in Section 6.
Stability analysis in the chemostat model entails evaluating the model's propensity to achieve and sustain a stable equilibrium state over a period of time. The purpose is to ascertain if the microbial populations and nutrient concentrations in the chemostat display stable dynamics or experience oscillations or instability. This study is essential for comprehending the model's ability to withstand environmental changes and disturbances, which in turn helps in optimizing operating conditions and designing resilient bioreactor models. Stability analysis offers valuable information on the long-term dynamics of microbial populations and nutrient levels in continuous culture models. This knowledge helps in making well-informed decisions in the fields of microbial ecology, biotechnology, and bioprocess engineering.
To obtain the FPs (N,P) of the model (1.6), we are required to solve the subsequent algebraic equations:
{N=Nexp(αC1+C−1),C=Cexp(−N1+C−1+βC). | (2.1) |
The model (1.6) possesses two FPs E1=(0,β) and E2=(α(1+β−αβ)1−α,1−1+α). The first FP, E1, is a boundary FP that exists always. The second FP, E2, is the unique positive FP if α>1+1β.
The Jacobian matrix of the model
{Nt+1=φ(Nt,Ct),Ct+1=ϕ(Nt,Ct), |
is the matrix given below:
J(N,C)=[∂φ∂N∂φ∂C∂ϕ∂N∂ϕ∂C]. |
Thus, the Jacobian matrix J(N,C) of the model (1.6) evaluated at any FP (N,C) is as follows:
J(N,C)=[e−1+Cα1+Ce−1+Cα1+CNα(1+C)2−e−1−N1+C+βCC1+Ce−1−N1+C+βC(C+(2+N)C2+C3−(1+C)2β)C(1+C)2]. |
The eigenvalues ξ1,2 of the Jacobian matrix J are helpful in determining the stability of FPs. The FP (x,y) is called a sink if |ξ1,2|<1 and a source if |ξ1|>1 and |ξ2|>1. Furthermore, the FP (x,y) is classified as a saddle point (SP) if |ξ1|>1 and |ξ2|<1 (or |ξ1|<1 and |ξ2|>1). In the case of a non-hyperbolic point (NHBP) (x,y), either |ξ1|=1 or |ξ2|=1. At a NHBP, the model experiences different types of bifurcations depending on the nature of the eigenvalues of the Jacobian matrix. The following result provides the topological classification of E1.
Theorem 2.1. The boundary FP, E1, is (1) a sink if α<1+1β, (2) never a source, (3) an SP if α>1+1β and (4) an NHBP if α=1+1β.
Proof. Through computations, it is obtained that
J(E1)=[e−1+αβ1+β0−β1+β0]. |
The eigenvalues of J(E0) are ξ1=e−1+αβ1+β and ξ2=0. One can easily check that
−1+αβ1+β{<0 if α<1+1β,=0 if α=1+1β,>0 if α>1+1β. | (2.2) |
If (4) is satisfied, then it follows that one of the eigenvalues of the matrix J(E1) is 1. Consequently, a TC bifurcation occurs at E1.
Remark 2.1. One key difference we want to highlight here is that the boundary FP, E1, of our discretized model is never a source. But in [1], the authors proved that E1 can be a source under certain conditions on the parameters. The Jacobian matrix of the continuous-time model (1.1) evaluated at E1 has eigenvalues ξ1=−1 and ξ2=−1+(−1+α)β1+β. Clearly, the first eigenvalue is negative. Thus, in the continuous-time model (1.1) the FP, E1, can never be a source. This comparison strengthens the fact that the piecewise constant argument method is better than the forward Euler method, as it shows more dynamic consistency with its continuous counterpart.
Theorem 2.2. Assume that α>1+1β. The positive FP, E2, is (1) a sink if 1−1+α<β<−1+3α1−2α+α2, (2) never a source, (3) an SP if β>−1+3α1−2α+α2and (4) an NHBP if β=−1+3α1−2α+α2.
Proof. The Jacobian matrix computed at E1 is given by
J(E2)=[1(−1+α)(−1+(−1+α)β)−1α−(−1+α)(−1+(−1+α)β)α]. |
The eigenvalues of J(E1) are ξ1=0 and ξ2=−1+β+α2β−2α(1+β)α. One can easily check that
|−1+β+α2β−2α(1+β)α|{<1 if 1−1+α<β<−1+3α1−2α+α2,=1 if β=−1+3α1−2α+α2,>1 if β>−1+3α1−2α+α2. | (2.3) |
If (4) is true, then it follows that one of the eigenvalues of the matrix J(E2) is −1. Consequently, a PD bifurcation occurs at E2.
Remark 2.2. One can check that the Jacobian matrix of the continuous-time model (1.1) evaluated at E2 has eigenvalues ξ1=−1 and ξ2=(1−α)(−1+(−1+α)β)α. Thus, in the continuous-time model (1.1) the positive FP, E2, can never be a source. We also obtain the same result. But in [1], the authors proved that E2 can be a source under certain conditions on parameters. This comparison strengthens the fact that the piecewise constant argument method is better than the forward Euler method, as it shows more dynamic consistency with its continuous counterpart.
● For α≤1+1β, the model (1.6) has only one FP, E1, which is stable until α<1+1β.
● For α>1+1β, the model possesses two FPs E1 and E2. From these two FPs, E1 is always unstable, but E2 is stable until 1−1+α<β<−1+3α1−2α+α2.
● The boundary FP, E1, always exists and is never destroyed. However, E1 interchanges its stability with the other FP E2 as we vary α about the critical value α=1+1β. Thus, the model experiences TC bifurcation at E1 if α=1+1β.
● The model experiences PD bifurcation at E2 if β=−1+3α1−2α+α2.
Remark 2.3. The stability classifications in our chemostat model reflect the long-term behavior of microbial populations and nutrient levels. A sink indicates a stable equilibrium where populations and nutrients reach a steady state. A saddle point is stable in one direction but unstable in another, resulting in possible alterations in dynamics under changing circumstances. A non-hyperbolic point is a key threshold at which little changes may cause substantial transitions, such as oscillations or instability, signifying the start of complex behaviors such as resource depletion or population surges.
The purpose of this section is to comprehensively examine bifurcation phenomena, specifically focusing on TC bifurcation at E1 and PD at E2. To get a thorough comprehension of bifurcation analysis, we recommend consulting the references [30,31]. Bifurcation plays a critical role in influencing the dynamics of the model, highlighting circumstances where even little alterations in parameters may result in substantial changes in the dynamics of the chemostat model. More precisely, it investigates how changes in factors like nutrition supply, dilution rate, or microbial growth rates result in the appearance of distinct stable conditions, periodic fluctuations, or chaotic patterns in the chemostat. Comprehending the bifurcation behavior allows for the fine-tuning of chemostat operational parameters in order to attain certain objectives, such as increasing biomass yield, limiting nutrient use, or ensuring the stability of microbial communities.
In this subsection, we investigate the TC bifurcation at E1 by considering condition (4) outlined in Theorem 2.1. Adding a sufficiently small perturbation parameter δ into the bifurcation parameter α in the neighborhood of α1=1+1β, the model (1.6) takes the subsequent form:
{Nt+1=Ntexp((α1+δ)Ct1+Ct−1),Ct+1=Ctexp(−Nt1+Ct−1+βCt). | (3.1) |
We shift E1 to (0,0) by taking Ut=Nt and Vt=Ct−β. Consequently, the system (3.1) is expressed as follows:
[Ut+1Vt+1]=[10−β1+β0][UtVt]+[Π1(Ut,Vt,δ)Π2(Ut,Vt,δ)], | (3.2) |
where
Π1(Ut,Vt,δ)=(1β+β2)UtVt+(β1+β)Utδ+(1−2β2β2(1+β)2)UtV2t+(β22(1+β)2)Utδ2+(2(1+β)2)UtVtδ+O((|Ut|+|Vt|+|δ|)4),Π2(Ut,Vt,δ)=(β2(1+β)2)U2t+12βV2t+(β(1+β)2)UtVt+(−β6(1+β)3)U3t−23β2V3t−(β(1+β)3)U2tVt−(1+2β+3β22β(1+β)3)UtV2t+O((|Ut|+|Vt|+|δ|)4). |
Next, the linear part of the system (3.2) is transformed into canonical form by using the subsequent transformation:
[UtVt]=[−1+ββ011][XtYt]. | (3.3) |
As a result, the system (3.2) reduces to the following form:
[Xt+1Yt+1]=[1000][XtYt]+[Γ(Xt,Yt,δ)Υ(Xt,Yt,δ)], | (3.4) |
where
Γ(Xt,Yt,δ)=c1X2t+c2XtYt+c3Xtδ+c4X3t+c5X2tYt+c6X2tδ+c7XtY2t+c8Xtδ2+c9XtYtδ+O((|Xt|+|Yt|+|δ|)4),Υ(Xt,Yt,δ)=d1Y2t+d2Xtδ+d3X3t+d4Y3t+d5X2tYt+d6X2tδ+d7XtY2t+d8Xtδ2+d9XtYtδ+O((|Xt|+|Yt|+|δ|)4),c1=1β+β2, c2=1β+β2, c3=β1+β, c4=1−2β2β2(1+β)2, c5=1−2ββ2(1+β)2,c6=2(1+β)2, c7=1−2β2β2(1+β)2, c8=β22(1+β)2, c9=2(1+β)2, d1=12β,d2=−β1+β, d3=−12β2(1+β)2, d4=−23β2, d5=−2+ββ2(1+β)2, d6=−2(1+β)2,d7=−(2+β)22β2(1+β)2, d8=−β22(1+β)2, d9=−2(1+β)2. |
Next, using the center manifold theory, we obtain the center manifold QC of (3.4) about (0,0) in a small neighborhood of δ=0. Thus, there exists a center manifold QC that can be expressed as follows:
QC={(Xt,Yt,δ)∈R3+|Yt=p1X2t+p2Xtδ+p3δ2+O((|Xt|+|δ|)3)}. |
Through calculations, we obtain that p1=0,p2=d2 and p3=0. Thus, the system (3.4) is restricted to QC in the manner as follows:
F:=Xt+1=Xt+c1X2t+c3Xtδ+c4X3t+(c6+c2d2)X2tδ+c8Xtδ2+O((|Xt|+|δ|)4). | (3.5) |
Through simple computations, we obtain
F(0,0)=0, FXt(0,0)=1, Fδ(0,0)=0,FXtXt(0,0)=2c1=2β+β2≠0,FXtδ(0,0)=c3=β1+β≠0. |
Thus, the model (1.6) experiences TC bifurcation at E1. The next result gives the conditions for the existence and direction of TC bifurcation in the model (1.6) at E1.
Theorem 3.1. If condition (4) in Theorem 2.1 is fulfilled, then (1.6) undergoes TC bifurcation at E1 when α differs in a close neighborhood of α1=1+1β.
Next, we investigate the PD bifurcation at E2 by considering condition (4) given in Theorem 2.2. Adding a sufficiently small perturbation parameter δ into the bifurcation parameter β in the neighborhood of β1=−1+3α1−2α+α2, the model (1.6) takes the subsequent form:
{Nt+1=Ntexp(α1Ct1+Ct−1),Ct+1=Ctexp(−Nt1+Ct−1+β1+δCt). | (3.6) |
We translate the positive FP, E2, to (0,0) by employing the following translation mapping:
Ut=Nt−α(1+β1+δ−α(β1+δ))1−α, Vt=Ct−1−1+α. |
Thus, the system (3.6) can be written as follows:
[Ut+1Vt+1]=[12α−1α−2][UtVt]+[Π1(Ut,Vt,δ)Π2(Ut,Vt,δ)], | (3.7) |
where
Π1(Ut,Vt,δ)=a1V2t+a2UtVt+a3Vtδ+a4V3t+a5UtV2t+a6V2tδ+O((|Ut|+|Vt|+|δ|)4),Π2(Ut,Vt,δ)=b1U2t+b2V2t+b3UtVt+b4Vtδ+b5U3t+b6V3t+b7U2tVt+b8UtV2t+b9V2tδ+b10UtVtδ+O((|Ut|+|Vt|+|δ|)4),a1=(−3+α)(−1+α), a2=−2+1α+α, a3=(−1+α)2, a4=(−1+α)2(13−8α+α2)3α,a5=(−3+α)(−1+α)32α2, a6=(−3+α)(−1+α)32α, b1=−1+α2α2, b2=(−1+α)(4+9α)2α,b3=(−1+α)(1+2α)α2, b4=−(−1+α)2α, b5=−(−1+α)26α3, b6=−(−1+α)2(2+6α+9α2)α2,b7=−(−1+α)2(1+α)α3, b8=−(−1+α)2(2+8α+9α2)2α3,b9=(−1+α)3(1+3α)α2, b10=(−1+α)3α2. |
Next, the linear part of system (3.7) is converted into canonical form by using the following transformation:
[UtVt]=[−α−2α11][XtYt]. | (3.8) |
Thus, the system (3.7) can be written as follows:
[Xt+1Yt+1]=[−1000][XtYt]+[Γ(Xt,Yt,δ)Υ(Xt,Yt,δ)], | (3.9) |
Γ(Xt,Yt,δ)=c1X2t+c2Y2t+c3XtYt+c4Xtδ+c5Ytδ+c6X3t+c7Y3t+c8X2tYt+c9X2tδ+c10XtY2t+c11Y2tδ+c12XtYtδ+O((|Xt|+|Yt|+|δ|)4),Υ(Xt,Yt,δ)=d1X2t+d2Y2t+d3XtYt+d4X3t+d5Y3t+d6X2tYt+d7X2tδ+d8XtY2t+d9Y2tδ+d10XtYtδ+O((|Xt|+|Yt|+|δ|)4),c1=6(−1+α), c2=−5+1α+4α, c3=−10+1α+9α, c4=−(−1+α)2α, c5=−(−1+α)2α,c6=−5(−1+α)2(−1+8α+13α2)6α2, c7=−2(−1+α)2(−2+4α+9α2)3α2,c8=−(−1+α)2(−3+14α+27α2)α2, c9=(−1+α)3(1+9α)2α2,c10=−(−1+α)2(−7+20α+45α2)2α2, c11=(−1+α)3(1+5α)2α2, c12=(−1+α)3(1+7α)α2,d1=4−1α−3α, d2=52−1α−3α2, d3=6−2α−4α, d4=11(−1+α)2(−1+2α+3α2)6α2,d5=2(−1+α)2(−2+α+5α2)3α2, d6=(−1+α)2(−5+7α+14α2)α2, |
d7=−(−1+α)3(−1+5α)2α2, d8=(−1+α)2(−9+8α+24α2)2α2,d9=−(−1+α)3(−1+3α)2α2, d10=−(−1+α)3(−1+4α)α2. |
Next, using the center manifold theory, we determine the center manifold QC of (3.9) about (0,0) in a small neighborhood of δ=0. Thus, there exists a center manifold QC that can be expressed as follows:
QC={(Xt,Yt,δ)∈R3+|Yt=p1X2t+p2Xtδ+p3δ2+O((|Xt|+|δ|)3)}. |
Through computations, we obtain that p1=d1,p2=0 and p3=0. Therefore, (3.9) is restricted to QC as follows:
˜F:=Xt+1=−Xt+c1X2t+c4Xtδ+(c6+c3d1)X3t+(c9+c5d1)X2tδ+O((|Xt|+|δ|)4). | (3.10) |
From [32], it can be seen that the conditions for PD bifurcation to occur are l1≠0 and l2≠0, where
l1=˜Fδ˜FXtXt+2˜FXtδ|(0,0)=2c4=−2(−1+α)2α≠0, | (3.11) |
l2=12(˜FXtXt)2+13˜FXtXtXt|(0,0)=2(c21+c6+c3d1)=−(−1+α)2(1−32α+11α2)3α2. | (3.12) |
From the above discussion, we obtain the next result:
Theorem 3.2. Suppose that condition (4) of Theorem 2.2 is fulfilled. The model (1.6) undergoes PD bifurcation at E2 if l2 given in (3.12) is non-zero and β varies in a small neighborhood of β1=−1+3α1−2α+α2. Moreover, if l2>0 (respectively, l2<0), then a period-2 orbit of the model (1.6) emanates and it is stable (respectively, unstable).
Remark 3.1. Conditions for transcritical and period-doubling bifurcations provide deeper insights into how parameter variations impact system stability. This theoretical advancement is crucial for understanding and predicting the system's long-term behavior under different scenarios.
Optimizing dynamical systems to accomplish specific performance objectives while avoiding chaotic behavior is a much-desired goal. Numerous areas of applied research and engineering widely use chaos control methods. Bifurcations and unstable oscilations have traditionally been considered negatively in the field of mathematical biology because they have a negative impact on biological populations' reproductive abilities. It is feasible to create a controler that can change the bifurcation properties of a dynamic model in order to attain certain dynamic qualities and properly govern chaos in the presence of bifurcations. There are many methods for managing chaos in discrete-time models. This section examines two distinct control methods: state feedback control and hybrid control strategies.
Initially, we apply the state feedback control method described in references [33,34]. The primary objective of the feedback control is to stabilize the positive fixed point E2, which becomes unstable as a result of bifurcation. The suggested approach entails transforming the chaotic system into a piecwise linear system to get an optimal controller that efficiently reduces the upper limit. Subsequently, the problem of optimization is resolved while maintaining compliance with certain limitations. The aforementioned method stabilizes chaotic trajectories at an unstable equilibrium point in the system (1.6). The controlled model via the feedback control method is given by
{Nt+1=Ntexp(αCt1+Ct−1)−Ht,Ct+1=Ctexp(−Nt1+Ct−1+βCt), | (4.1) |
where Ht=κ1(Nt−α(1+β−αβ)1−α)+κ2(Ct−1−1+α) is the feedback controlling force with feedback gains κ1 and κ2. Feedback control aims to adjust the system by adding a corrective term Ht to the state variable Nt+1. Given that Ct+1 is influenced by Nt, introducing a control term in the second equation is not necessary. The values of the feedback gains are chosen to control the stability of E2 by moving eigenvalues of J(E2) into the unit circle. The Jacobian matrix of the model (4.1) at E2 is given by
J(E2)=[1−κ1−κ2+(−1+α)(−1+(−1+α)β)−1α−(−1+α)(−1+(−1+α)β)α]. | (4.2) |
The characteristic equation of J(E2) is
ξ2+(κ1+αβ−2(1+β)+1+βα)ξ+−κ2+κ1(−1+α)(−1+(−1+α)β)α=0. | (4.3) |
Let ξ1 and ξ2 be the roots of (4.3), and we have
ξ1+ξ2=2−κ1+2β−αβ−1+βα, | (4.4) |
ξ1ξ2=−κ2+κ1(−1+α)(−1+(−1+α)β)α. | (4.5) |
Then, the lines of marginal stability are obtained by solving ξ1=±1 and ξ1ξ2=1. These conditions ensure that |ξ1,2|<1. Suppose that ξ1ξ2=1, and then Eq (4.5) implies that
L1:((−1+α)(−1+(−1+α)β)α)κ1−1ακ2−1=0. | (4.6) |
Next, considering ξ1=1 in Eqs (4.4) and (4.5), we get
L2:−(1+(−1+α)2βα)κ1+1ακ2−(−1+α)(−1+(−1+α)β)α=0. | (4.7) |
Next, considering ξ1=−1 in Eqs (4.4) and (4.5), we get
L3:(αβ−2(1+β)+1+βα)κ1−1ακ2−1−3α+β−2αβ+α2βα=0. | (4.8) |
The stable eigenvalues are contained inside the area bounded by the lines L1,L2, and L3. This stability region provides guidance for selecting feedback gain values that will achieve the desired stability. By selecting κ1 and κ2 within this region, it is possible to effectively manage the chaotic dynamics and stabilize the bifurcation.
Subsequently, we use the hybrid control approach [35] to effectively regulate chaos under the influence of both types of bifurcation effects. The hybrid control approach is a technique that integrates state feedback and parameter perturbation to stabilize unstable periodic orbits inside a model's chaotic attractor. The hybrid control method uses a weighted combination of the original model dynamics and the controlled model. The controlled model of (1.6) via the hybrid control method is given by
{Nt+1=ρNtexp(αCt1+Ct−1)+(1−ρ)Nt,Ct+1=ρCtexp(−Nt1+Ct−1+βCt)+(1−ρ)Ct, | (4.9) |
where 0<ρ<1. This interpolation allows for flexible control that adjusts the dynamics without introducing aggressive or unrealistic modifications. The control parameter ρ in the hybrid method serves as a fine-tuning mechanism, effectively balancing stability enhancement with dynamic consistency. The parameter ρ serves as a control parameter, regulating the impact of the original model (1.6) and the modified model (4.9). The negative values of ρ indicate a contrasting effect from the initial model (1.6). If the values of ρ exceed 1, the original model (1.6) may have an exaggerated effect that goes beyond its inherent influence, perhaps resulting in unrealistic or impractical results in the controlled model (4.9). Models (1.6) and (4.9) have identical FPs. The Jacobian matrix of the controlled model evaluated at positive FP, E2, is given by
J(E2)=[1(−1+α)(−1+(−1+α)β)ρ−ρα1+2βρ−αβρ−(1+β)ρα]. | (4.10) |
The eigenvalues of J(E2) are ξ1=1−ρ and ξ2=1−(−1+α)(−1+(−1+α)β)ρα. Thus, we obtain the following result:
Theorem 4.1. Assume that α>1+1β. The FP, E2, of the controlled model (4.9) is a sink if
1−1+α<β<2α−ρ+αρρ−2αρ+α2ρ. |
The purpose of this section is to verify our theoretical results through numerical simulations. Calculations are done using MATHEMATICA, while MATLAB is used for plotting graphs.
Assume that α=2.5 and vary β. For this, the TC bifurcation value is β=1−1+α≈0.66667. The corresponding boundary FP for β≈0.66667 is obtained as E1≈(0,0.66667). By fixing the initial conditions as N0=0.01 and C0=0.7, and varying β∈[0.64,0.70], the bifurcation diagrams, Figures (1a) and (1b), are plotted. These figures illustrate that the model (1.6) has only one FP, E1, for β<0.66667 which is stable. Later, for β>0.66667, the model possesses two FPs E1 and E2, where E1 is unstable and E2 is stable. Thus, at β≈0.66667, the two FPs E1 and E2 collide and exchange their stability. This confirms that the model (1.6) experiences TC bifurcation at E1. This verifies Theorem 2.1. Next, for α=2.5, the PD bifurcation value is β1≈2.88889. The positive FP is obtained as E2≈(5.55556,0.66667). The eigenvalues of the Jacobian matrix J(E2) are ξ1=−1 and ξ2=0. Thus, the model (1.6) undergoes PD bifurcation at E2 when β≈2.88889. Figures (1c) and (1d) depict bifurcation diagrams by fixing N0=5.6 and C0=0.7, and varying β∈[2.84,2.96]. The positive FP, E2, is a sink if β<2.88889. It loses its stability due to PD bifurcation for β≥2.88889. This confirms Theorem 2.2.
Next, our target is to check the effectiveness of the feedback control method. Considering the values α=2.5 and β=2.9, together with the initial conditions N0=5.6 and C0=0.7 for the controlled model (4.1), and the lines of marginal stability are obtained as follows:
L1:κ2=−2.5+5.025κ1, |
L2:κ2=5.025+7.525κ1, |
and
L3:κ2=−0.025+2.525κ1. |
Figure (2a) illustrates the stability region of model (4.1), which is bounded by the lines L1,L2, and L3. In Figures (1c) and (1d), one can see that the original model (1.6) experiences PD bifurcation at E2 for the considered paramatric values. The controlled model (4.1) is analyzed using feedback gains κ1=−0.8 and κ2=−3. Figure 2 shows the graph of xn in Figure (2c), yn in Figure (2d), and the phase portrait in Figure (2b) for the model (4.1). Therefore, it can be concluded that the feedback control technique is effective in managing bifurcation and chaotic behavior.
We shall now evaluate the efficacy of the hybrid control technique. We assume ρ=0.98,α=2.5, and vary β for the controlled model (4.9). If 0.66667<β<2.93424, the positive FP, E2, is stable. One can observe that the stability region has been expanded in the controlled model (4.9). The bifurcation diagrams are depicted in Figures (3a) and (3b) by taking the initial values N0=5.6 and C0=0.7 and varying β∈[2.84,2.99].
Remark 5.1. The numerical simulations clearly demonstrate the shift from stability to instability, highlighting how feedback and hybrid control approaches successfully postpone or eradicate bifurcations, thereby stabilizing the system. This not only validates the theoretical predictions but also underscores the practical importance of these control techniques in preventing chaos, which is crucial for applications such as optimizing bioreactor operations.
This study investigates the dynamics of the discrete-time chemostat model (1.6). The discretization is obtained using the piecewise constant argument method. In [1], the authors investigated the same model by discretizing using the forward Euler method. But the forward Euler method is less dynamically consistent than the piecewise constant argument method. In this study, the existence and stability of all possible FPs are investigated. Furthermore, TC bifurcation at boundary FP, E1, and PD bifurcation at the positive FP, E2, are investigated using center manifold and bifurcation theory. The existence conditions and direction of TC and PD bifurcations are investigated. The thorough bifurcation analysis provides new insights into the system's stability behavior, revealing particular circumstances that cause transcritical and period-doubling bifurcations. These results provide important theoretical insights into the study of discrete-time biological models. Additionally, feedback control and hybrid control methods are employed to control bifurcation and chaos in the model. The application of feedback and hybrid chaos control methods effectively mitigates instability, further underscoring the practical relevance of the model. Some numerical simulations are presented to verify our theoretical results.
Although we obtain similar results as in [1], the following are the main differences observed in both models (1.6) and (1.2):
● The model (1.2) can produce negative values of N and C for specific initial conditions and parametric values. This is not reasonable biologically. But our model (1.6) does not produce negative values of N and C.
● The FPs, E1 and E2, of our discretized model can never be a source. But in [1], the authors proved that E1 and E2 can be a source under certain conditions on parameters. Their result is not consistent with the continuous-time model (1.1).
● We investigated TC bifurcation at E1 which was not reported in [1].
This study illustrates that the piecewise constant argument method is more dynamically consistent than the forward Euler method. Future research may expand this model to include supplementary environmental variables or investigate alternate discretization techniques to improve dynamic accuracy. Furthermore, using these control strategies in experimental or industrial settings may provide practical insights for enhancing microbial growth systems in controlled environments.
The author declares that there are no conflicts of interest in this paper.
[1] |
K. S. N. Al-Basyouni, A. Q. Khan, Bifurcation analysis of a discrete-time chemostat model, Math. Probl. Eng., 2023 (2023), 7518261. https://doi.org/10.1155/2023/7518261 doi: 10.1155/2023/7518261
![]() |
[2] |
D. Zhang, L. Wang, Multistability driven by inhibitory kinetics in a discrete-time size-structured chemostat model, Chaos, 29 (2019), 063112. https://doi.org/10.1063/1.5096661 doi: 10.1063/1.5096661
![]() |
[3] |
M. Zhao, C. Li, J. Wang, Complex dynamic behaviors of a discrete-time predator-prey system, J. Appl. Anal. Comput., 7 (2017), 478–500. https://doi.org/10.11948/2017030 doi: 10.11948/2017030
![]() |
[4] |
A. Khan, I. M. Alsulami, Discrete Leslie's model with bifurcations and control, AIMS Math., 8 (2023), 22483–22506. https://doi.org/10.3934/math.20231146 doi: 10.3934/math.20231146
![]() |
[5] |
S. Akhtar, R. Ahmed, M. Batool, N. A. Shah, J. D. Chung, Stability, bifurcation and chaos control of a discretized Leslie prey-predator model, Chaos Soliton. Fract., 152 (2021), 111345. https://doi.org/10.1016/j.chaos.2021.111345 doi: 10.1016/j.chaos.2021.111345
![]() |
[6] |
P. A. Naik, Z. Eskandari, Z. Avazzadeh, J. Zu, Multiple bifurcations of a discrete-time prey-predator model with mixed functional response, Int. J. Bifurcat. Chaos, 32 (2022), 2250050. https://doi.org/10.1142/s021812742250050x doi: 10.1142/s021812742250050x
![]() |
[7] |
P. A. Naik, Z. Eskandari, A. Madzvamuse, Z. Avazzadeh, J. Zu, Complex dynamics of a discrete-time seasonally forced SIR epidemic model, Math. Method. Appl. Sci., 46 (2023), 7045–7059. https://doi.org/10.1002/mma.8955 doi: 10.1002/mma.8955
![]() |
[8] |
Z. Jing, J. Yang, Bifurcation and chaos in discrete-time predator-prey system, Chaos Soliton. Fract., 27 (2006), 259–277. https://doi.org/10.1016/j.chaos.2005.03.040 doi: 10.1016/j.chaos.2005.03.040
![]() |
[9] |
A. Suleman, R. Ahmed, F. S. Alshammari, N. A. Shah, Dynamic complexity of a slow-fast predator-prey model with herd behavior, AIMS Math., 8 (2023), 24446–24472. https://doi.org/10.3934/math.20231247 doi: 10.3934/math.20231247
![]() |
[10] |
A. Suleman, A. Q. Khan, R. Ahmed, Bifurcation analysis of a discrete Leslie-gower predator-prey model with slow-fast effect on predator, Math. Method. Appl. Sci., 47 (2024), 8561–8580. https://doi.org/10.1002/mma.10032 doi: 10.1002/mma.10032
![]() |
[11] |
R. Ahmed, N. Tahir, N. A. Shah, An analysis of the stability and bifurcation of a discrete-time predator-prey model with the slow-fast effect on the predator, Chaos, 34 (2024), 033127. https://doi.org/10.1063/5.0185809 doi: 10.1063/5.0185809
![]() |
[12] | J. Wiener, Generalized solutions of functional differential equations, World Scientific, 1993. https://doi.org/10.1142/9789814343183 |
[13] |
Q. Zhou, F. Chen, S. Lin, Complex dynamics analysis of a discrete amensalism system with a cover for the first species, Axioms, 11 (2022), 365. https://doi.org/10.3390/axioms11080365 doi: 10.3390/axioms11080365
![]() |
[14] |
R. Ahmed, S. Akhtar, U. Farooq, S. Ali, Stability, bifurcation, and chaos control of predator-prey system with additive Allee effect, Commun. Math. Biol. Neurosci., 2023 (2023), 9. https://doi.org/10.28919/cmbn/7824 doi: 10.28919/cmbn/7824
![]() |
[15] |
P. Amster, G. Robledo, D. Sepulveda, Dynamics of a discrete size-structured chemostat with variable nutrient supply, Discrete Cont. Dyn. B, 28 (2023), 4937–4967. https://doi.org/10.3934/dcdsb.2023048 doi: 10.3934/dcdsb.2023048
![]() |
[16] |
M. S. Shabbir, Q. Din, M. Safeer, M. A. Khan, K. Ahmad, A dynamically consistent nonstandard finite difference scheme for a predator-prey model, Adv. Differ. Equ., 2019 (2019), 381. https://doi.org/10.1186/s13662-019-2319-6 doi: 10.1186/s13662-019-2319-6
![]() |
[17] |
R. Ahmed, A. Ahmad, N. Ali, Stability analysis and Neimark-sacker bifurcation of a nonstandard finite difference scheme for Lotka-Volterra prey-predator model, Commun. Math. Biol. Neurosci., 2022 (2022), 61. https://doi.org/10.28919/cmbn/7534 doi: 10.28919/cmbn/7534
![]() |
[18] |
N. Bairagi, M. Biswas, A predator-prey model with Beddington-DeAngelis functional response: a non-standard finite-difference method, J. Differ. Equ. Appl., 22 (2016), 581–593. https://doi.org/10.1080/10236198.2015.1111345 doi: 10.1080/10236198.2015.1111345
![]() |
[19] |
M. U. Akhmet, Stability of differential equations with piecewise constant arguments of generalized type, Nonlinear Anal. Theor., 68 (2008), 794–803. https://doi.org/10.1016/j.na.2006.11.037 doi: 10.1016/j.na.2006.11.037
![]() |
[20] |
M. S. Alwan, X. Liu, W. C. Xie, Comparison principle and stability of differential equations with piecewise constant arguments, J. Franklin I., 350 (2013), 211–230. https://doi.org/10.1016/j.jfranklin.2012.08.016 doi: 10.1016/j.jfranklin.2012.08.016
![]() |
[21] |
A. Q. Khan, Global dynamics, bifurcation analysis, and chaos in a discrete Kolmogorov model with piecewise-constant argument, Math. Probl. Eng., 2021 (2021), 5259226. https://doi.org/10.1155/2021/5259226 doi: 10.1155/2021/5259226
![]() |
[22] |
F. Bozkurt, A. Yousef, H. Bilgil, D. Baleanu, A mathematical model with piecewise constant arguments of colorectal cancer with chemo-immunotherapy, Chaos Soliton. Fract., 168 (2023), 113207. https://doi.org/10.1016/j.chaos.2023.113207 doi: 10.1016/j.chaos.2023.113207
![]() |
[23] |
P. A. Naik, Y. Javaid, R. Ahmed, Z. Eskandari, A. H. Ganie, Stability and bifurcation analysis of a population dynamic model with Allee effect via piecewise constant argument method, J. Appl. Math. Comput., 70 (2024), 4189–4218. https://doi.org/10.1007/s12190-024-02119-y doi: 10.1007/s12190-024-02119-y
![]() |
[24] |
A. Khan, I. M. Alsulami, Complicate dynamical analysis of a discrete predator-prey model with a prey refuge, AIMS Math., 8 (2023), 15035–15057. https://doi.org/10.3934/math.2023768 doi: 10.3934/math.2023768
![]() |
[25] |
A. A. Khabyah, R. Ahmed, M. S. Akram, S. Akhtar, Stability, bifurcation, and chaos control in a discrete predator-prey model with strong Allee effect, AIMS Math., 8 (2023), 8060–8081. https://doi.org/10.3934/math.2023408 doi: 10.3934/math.2023408
![]() |
[26] |
A. Khan, I. M. Alsulami, S. Hamdani, Controlling the chaos and bifurcations of a discrete prey-predator model, AIMS Math., 9 (2024), 1783–1818. https://doi.org/10.3934/math.2024087 doi: 10.3934/math.2024087
![]() |
[27] |
A. Q. Khan, I. M. Alsulami, U. Sadiq, Stability, chaos, and bifurcation analysis of a discrete chemical system, Complexity, 2022 (2022), 6921934. https://doi.org/10.1155/2022/6921934 doi: 10.1155/2022/6921934
![]() |
[28] |
X. Jiang, X. Chen, M. Chi, J. Chen, On Hopf bifurcation and control for a delay systems, Appl. Math. Comput., 370 (2020), 124906. https://doi.org/10.1016/j.amc.2019.124906 doi: 10.1016/j.amc.2019.124906
![]() |
[29] |
X. Jiang, X. Chen, T. Huang, H. Yan, Bifurcation and control for a predator-prey system with two delays, IEEE T. Circuits II, 68 (2021), 376–380. https://doi.org/10.1109/tcsii.2020.2987392 doi: 10.1109/tcsii.2020.2987392
![]() |
[30] | J. Guckenheimer, P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, New York: Springer, 1983. https://doi.org/10.1007/978-1-4612-1140-2 |
[31] | S. Wiggins, Introduction to applied nonlinear dynamical systems and chaos, New York: Springer, 2003. https://doi.org/10.1007/b97481 |
[32] | Y. A. Kuznetsov, Elements of applied bifurcation theory, New York: Springer, 2004. https://doi.org/10.1007/978-1-4757-3978-7 |
[33] | G. Chen, X. Dong, From chaos to order, World Scientific, 1998. https://doi.org/10.1142/3033 |
[34] |
C. Lei, X. Han, W. Wang, Bifurcation analysis and chaos control of a discrete-time prey-predator model with fear factor, Math. Biosci. Eng., 19 (2022), 6659–6679. https://doi.org/10.3934/mbe.2022313 doi: 10.3934/mbe.2022313
![]() |
[35] |
X. S. Luo, G. Chen, B. H. Wang, J. Q. Fang, Hybrid control of period-doubling bifurcation and chaos in discrete nonlinear dynamical systems, Chaos Soliton. Fract., 18 (2003), 775–783. https://doi.org/10.1016/s0960-0779(03)00028-6 doi: 10.1016/s0960-0779(03)00028-6
![]() |