On the stability, chaos and bifurcation analysis of a discrete-time chemostat model using the piecewise constant argument method

1.   Introduction

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]:

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, $\alpha$ quantifies the intrinsic growth rate of bacterial populations, and $\beta$ represents the rate at which fresh nutrients are supplied to the chemostat. The parameters $\alpha$ and $\beta$ 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:

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:

where $[t]$ represents the integer part of $t$, and $0 < t < \infty$. Furthermore, integrating system (1.3) on an interval $[n, n+1)$ with $n = 0, 1, 2, \cdots$ yields the following system:

Taking $t \to n+1$, we obtain the following discrete-time system:

After replacing $n$ by $t$, we obtain

The model (1.2) can produce negative values for $N(t)$ and $C(t)$, particularly when $h$ is quite big and $\alpha$ 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.

2.   The existence and local stability analysis of FPs

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:

The model (1.6) possesses two FPs $E_1 = (0, \beta)$ and $E_2 = (\frac{\alpha(1+\beta-\alpha \beta)}{1-\alpha}, \frac{1}{-1+\alpha})$. The first FP, $E_1$, is a boundary FP that exists always. The second FP, $E_2$, is the unique positive FP if $\alpha > 1+\frac{1}{\beta}$.

The Jacobian matrix of the model

is the matrix given below:

Thus, the Jacobian matrix $J(N, C)$ of the model (1.6) evaluated at any FP $(N, C)$ is as follows:

The eigenvalues $\xi_{1, 2}$ of the Jacobian matrix $J$ are helpful in determining the stability of FPs. The FP $(x, y)$ is called a sink if $|\xi_{1, 2}| < 1$ and a source if $|\xi_1| > 1$ and $|\xi_2| > 1$. Furthermore, the FP $(x, y)$ is classified as a saddle point (SP) if $|\xi_1| > 1$ and $|\xi_2| < 1$ (or $|\xi_1| < 1$ and $|\xi_2| > 1$). In the case of a non-hyperbolic point (NHBP) $(x, y)$, either $|\xi_1| = 1$ or $|\xi_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 $E_1$.

Theorem 2.1. The boundary FP, $E_1$, is $(1)$ a sink if $\alpha < 1+\frac{1}{\beta}$, $(2)$ never a source, $(3)$ an SP if $\alpha > 1+\frac{1}{\beta}$ and $(4)$ an NHBP if $\alpha = 1+\frac{1}{\beta}$.

Proof. Through computations, it is obtained that

The eigenvalues of $J(E_0)$ are $\xi_1 = e^{-1+\frac{\alpha \beta }{1+\beta }}$ and $\xi_2 = 0$. One can easily check that

□

If $(4)$ is satisfied, then it follows that one of the eigenvalues of the matrix $J(E_1)$ is $1$. Consequently, a TC bifurcation occurs at $E_1$.

Remark 2.1. One key difference we want to highlight here is that the boundary FP, $E_1$, of our discretized model is never a source. But in ^[1], the authors proved that $E_1$ can be a source under certain conditions on the parameters. The Jacobian matrix of the continuous-time model (1.1) evaluated at $E_1$ has eigenvalues $\xi_1 = -1$ and $\xi_2 = \frac{-1+(-1+\alpha)\beta}{1+\beta}$. Clearly, the first eigenvalue is negative. Thus, in the continuous-time model (1.1) the FP, $E_1$, 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 $\alpha > 1+\frac{1}{\beta}$. The positive FP, $E_2$, is $(1)$ a sink if $\frac{1}{-1+\alpha} < \beta < \frac{-1+3\alpha}{1-2\alpha+\alpha^2}$, $(2)$ never a source, $(3)$ an SP if $\beta > \frac{-1+3\alpha}{1-2\alpha+\alpha^2}$and $(4)$ an NHBP if $\beta = \frac{-1+3\alpha}{1-2\alpha+\alpha^2}$.

Proof. The Jacobian matrix computed at $E_1$ is given by

The eigenvalues of $J(E_1)$ are $\xi_1 = 0$ and $\xi_2 = -\frac{1+\beta +\alpha ^2 \beta -2 \alpha (1+\beta)}{\alpha }$. One can easily check that

□

If $(4)$ is true, then it follows that one of the eigenvalues of the matrix $J(E_2)$ is $-1$. Consequently, a PD bifurcation occurs at $E_2$.

Remark 2.2. One can check that the Jacobian matrix of the continuous-time model (1.1) evaluated at $E_2$ has eigenvalues $\xi_1 = -1$ and $\xi_2 = \frac{(1-\alpha)(-1+(-1+\alpha)\beta)}{\alpha}$. Thus, in the continuous-time model (1.1) the positive FP, $E_2$, can never be a source. We also obtain the same result. But in ^[1], the authors proved that $E_2$ 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.

2.1.   Summary

● For $\alpha\leq 1+\frac{1}{\beta}$, the model (1.6) has only one FP, $E_1$, which is stable until $\alpha < 1+\frac{1}{\beta}$.

● For $\alpha > 1+\frac{1}{\beta}$, the model possesses two FPs $E_1$ and $E_2$. From these two FPs, $E_1$ is always unstable, but $E_2$ is stable until $\frac{1}{-1+\alpha} < \beta < \frac{-1+3\alpha}{1-2\alpha+\alpha^2}$.

● The boundary FP, $E_1$, always exists and is never destroyed. However, $E_1$ interchanges its stability with the other FP $E_2$ as we vary $\alpha$ about the critical value $\alpha = 1+\frac{1}{\beta}$. Thus, the model experiences TC bifurcation at $E_1$ if $\alpha = 1+\frac{1}{\beta}$.

● The model experiences PD bifurcation at $E_2$ if $\beta = \frac{-1+3\alpha}{1-2\alpha+\alpha^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.

3.   Bifurcation analysis

The purpose of this section is to comprehensively examine bifurcation phenomena, specifically focusing on TC bifurcation at $E_1$ and PD at $E_2$. 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.

3.1.   TC bifurcation at $E_1$

In this subsection, we investigate the TC bifurcation at $E_1$ by considering condition $(4)$ outlined in Theorem 2.1. Adding a sufficiently small perturbation parameter $\delta$ into the bifurcation parameter $\alpha$ in the neighborhood of $\alpha_1 = 1+\frac{1}{\beta}$, the model (1.6) takes the subsequent form:

We shift $E_1$ to $(0, 0)$ by taking $U_t = N_t$ and $\ V_t = C_t-\beta$. Consequently, the system (3.1) is expressed as follows:

where

Next, the linear part of the system (3.2) is transformed into canonical form by using the subsequent transformation:

As a result, the system (3.2) reduces to the following form:

where

Next, using the center manifold theory, we obtain the center manifold $Q^C$ of (3.4) about $(0, 0)$ in a small neighborhood of $\delta = 0$. Thus, there exists a center manifold $Q^C$ that can be expressed as follows:

Through calculations, we obtain that $p_1 = 0, p_2 = d_2$ and $p_3 = 0$. Thus, the system (3.4) is restricted to $Q^C$ in the manner as follows:

Through simple computations, we obtain

Thus, the model (1.6) experiences TC bifurcation at $E_1$. The next result gives the conditions for the existence and direction of TC bifurcation in the model (1.6) at $E_1$.

Theorem 3.1. If condition $(4)$ in Theorem 2.1 is fulfilled, then (1.6) undergoes TC bifurcation at $E_1$ when $\alpha$ differs in a close neighborhood of $\alpha_1 = 1+\frac{1}{\beta}$.

3.2.   PD bifurcation at $E_2$

Next, we investigate the PD bifurcation at $E_2$ by considering condition $(4)$ given in Theorem 2.2. Adding a sufficiently small perturbation parameter $\delta$ into the bifurcation parameter $\beta$ in the neighborhood of $\beta_1 = \frac{-1+3\alpha}{1-2\alpha+\alpha^2}$, the model (1.6) takes the subsequent form:

We translate the positive FP, $E_2$, to $(0, 0)$ by employing the following translation mapping:

Thus, the system (3.6) can be written as follows:

where

Next, the linear part of system (3.7) is converted into canonical form by using the following transformation:

Thus, the system (3.7) can be written as follows:

Next, using the center manifold theory, we determine the center manifold $Q^C$ of (3.9) about $(0, 0)$ in a small neighborhood of $\delta = 0$. Thus, there exists a center manifold $Q^C$ that can be expressed as follows:

Through computations, we obtain that $p_1 = d_1, p_2 = 0$ and $p_3 = 0$. Therefore, (3.9) is restricted to $Q^C$ as follows:

From ^[32], it can be seen that the conditions for PD bifurcation to occur are $l_1 \neq 0$ and $l_2 \neq 0$, where

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 $E_2$ if $l_2$ given in (3.12) is non-zero and $\beta$ varies in a small neighborhood of $\beta_1 = \frac{-1+3\alpha}{1-2\alpha+\alpha^2}$. Moreover, if $\ l_2 > 0$ (respectively, $l_2 < 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.

4.   Bifurcation and chaos control

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 $E_2$, 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

where $H_t = \kappa_1 \bigg(N_t-\frac{\alpha(1+\beta-\alpha \beta)}{1-\alpha}\bigg)+\kappa_2 \bigg(C_t-\frac{1}{-1+\alpha}\bigg)$ is the feedback controlling force with feedback gains $\kappa_1$ and $\kappa_2$. Feedback control aims to adjust the system by adding a corrective term $H_t$ to the state variable $N_{t+1}$. Given that $C_{t+1}$ is influenced by $N_t$, introducing a control term in the second equation is not necessary. The values of the feedback gains are chosen to control the stability of $E_2$ by moving eigenvalues of $J(E_2)$ into the unit circle. The Jacobian matrix of the model (4.1) at $E_2$ is given by

The characteristic equation of $J(E_2)$ is

Let $\xi_1$ and $\xi_2$ be the roots of (4.3), and we have

Then, the lines of marginal stability are obtained by solving $\xi_1 = \pm 1$ and $\xi_1 \xi_2 = 1$. These conditions ensure that $|\xi_{1, 2}| < 1$. Suppose that $\xi_1 \xi_2 = 1$, and then Eq (4.5) implies that

Next, considering $\xi_1 = 1$ in Eqs (4.4) and (4.5), we get

Next, considering $\xi_1 = -1$ in Eqs (4.4) and (4.5), we get

The stable eigenvalues are contained inside the area bounded by the lines $L_1, L_2,$ and $L_3$. This stability region provides guidance for selecting feedback gain values that will achieve the desired stability. By selecting $\kappa_1$ and $\kappa_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

where $0 < \rho < 1$. This interpolation allows for flexible control that adjusts the dynamics without introducing aggressive or unrealistic modifications. The control parameter $\rho$ in the hybrid method serves as a fine-tuning mechanism, effectively balancing stability enhancement with dynamic consistency. The parameter $\rho$ serves as a control parameter, regulating the impact of the original model (1.6) and the modified model (4.9). The negative values of $\rho$ indicate a contrasting effect from the initial model (1.6). If the values of $\rho$ 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, $E_2$, is given by

The eigenvalues of $J(E_2)$ are $\xi_1 = 1-\rho$ and $\xi_2 = 1-\frac{(-1+\alpha) (-1+(-1+\alpha) \beta) \rho }{\alpha }$. Thus, we obtain the following result:

Theorem 4.1. Assume that $\alpha > 1+\frac{1}{\beta}$. The FP, $E_2$, of the controlled model (4.9) is a sink if

5.   Numerical simulations

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 $\alpha = 2.5$ and vary $\beta$. For this, the TC bifurcation value is $\beta = \frac{1}{-1+\alpha}\approx 0.66667$. The corresponding boundary FP for $\beta \approx 0.66667$ is obtained as $E_1 \approx (0, 0.66667)$. By fixing the initial conditions as $N_0 = 0.01$ and $C_0 = 0.7$, and varying $\beta \in [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, $E_1$, for $\beta < 0.66667$ which is stable. Later, for $\beta > 0.66667$, the model possesses two FPs $E_1$ and $E_2$, where $E_1$ is unstable and $E_2$ is stable. Thus, at $\beta \approx 0.66667$, the two FPs $E_1$ and $E_2$ collide and exchange their stability. This confirms that the model (1.6) experiences TC bifurcation at $E_1$. This verifies Theorem 2.1. Next, for $\alpha = 2.5$, the PD bifurcation value is $\beta_1\approx 2.88889$. The positive FP is obtained as $E_2\approx (5.55556, 0.66667)$. The eigenvalues of the Jacobian matrix $J(E_2)$ are $\xi_1 = -1$ and $\xi_2 = 0$. Thus, the model (1.6) undergoes PD bifurcation at $E_2$ when $\beta \approx 2.88889$. Figures (1c) and (1d) depict bifurcation diagrams by fixing $N_0 = 5.6$ and $C_0 = 0.7$, and varying $\beta \in [2.84, 2.96]$. The positive FP, $E_2$, is a sink if $\beta < 2.88889$. It loses its stability due to PD bifurcation for $\beta \geq 2.88889$. This confirms Theorem 2.2.

Next, our target is to check the effectiveness of the feedback control method. Considering the values $\alpha = 2.5$ and $\beta = 2.9$, together with the initial conditions $N_0 = 5.6$ and $C_0 = 0.7$ for the controlled model (4.1), and the lines of marginal stability are obtained as follows:

and

Figure (2a) illustrates the stability region of model (4.1), which is bounded by the lines $L_1, L_2,$ and $L_3$. In Figures (1c) and (1d), one can see that the original model (1.6) experiences PD bifurcation at $E_2$ for the considered paramatric values. The controlled model (4.1) is analyzed using feedback gains $\kappa_1 = -0.8$ and $\kappa_2 = -3$. Figure 2 shows the graph of $x_n$ in Figure (2c), $y_n$ 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 $\rho = 0.98, \alpha = 2.5$, and vary $\beta$ for the controlled model (4.9). If $0.66667 < \beta < 2.93424$, the positive FP, $E_2$, 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 $N_0 = 5.6$ and $C_0 = 0.7$ and varying $\beta \in [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.

6.   Conclusions

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, $E_1$, and PD bifurcation at the positive FP, $E_2$, 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, $E_1$ and $E_2$, of our discretized model can never be a source. But in ^[1], the authors proved that $E_1$ and $E_2$ 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 $E_1$ 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.

Conflict of interest

The author declares that there are no conflicts of interest in this paper.