# Bifurcation Analysis of a Discrete-Time Chemostat Model

Authors: Khalil Salem Nemer Al-Basyouni; Abdul Qadeer Khan
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

In this article, a discretized two-dimensional chemostat model is investigated. By using the algebraic technique, it is proved that the discrete chemostat model has semitrivial equilibrium solution for all involved parametric values, but it has an interior equilibrium solution under certain parametric conditions. We have studied the local dynamics with topological classifications about the semitrivial and interior equilibrium solution on the basis of the theory of the discrete dynamical system. By bifurcation theory, we studied the existence of bifurcations for the understudied discrete chemostat model and proved that at a semitrivial equilibrium solution, it did not undergo flip bifurcation, but it undergoes flip bifurcation at interior equilibrium solution, and no other bifurcation occurs at this point. It has also explored the existence of periodic points and chaos control by the state feedback method. Finally, with the help of MATLAB software, flip bifurcation diagrams and corresponding maximum Lyapunov exponents were also presented. These numerical simulations not only show the correctness of theoretical results but also reveal the complex dynamics such as the orbits of period−7,−8,−13, and −16.


---

## Body

## 1. Introduction

### 1.1. Motivation and Mathematical Modelling of the Discrete Chemostat Model

In a biological laboratory, chemostat is a basic device which also occupies a central role in the field of mathematics. It is a model of simple lake to study exploitation competition. This model is also used to tract the water wastes process. Oceans are large water bodies which contain major fragments of life form on the Earth including microorganisms. The earliest life form on the planet Earth first appeared in oceans which was a microorganism. Similarly, the first life on the Earth was also migrated from oceans. Chemostat is a device that is used to study physiological interaction of microorganisms under the conditions when nutrients become limited and observed in oligotrophic oceans. These microorganisms can be used to study carbonate chemistry. The effect of accretion on the phytoplankton’s growth rate and general physiology of phytoplankton species should be assessed in the experiment. Chemostat is used for culturing these cyanobacteria to access carbonate chemistry of oceans, effect of the flow rate, temperature,pco2, nutrient concentration, and composition. A few chemostat investigations have been carried out at sea to study the effects of constant low supply of the limiting nutrient iron on the survival and growth of phytoplankton.Chemostat provides most favorable environment for microorganisms to develop into a biomass. That is why, they can interact with raw materials forming a product. It is noted that fermentation takes place in a fermenter. There are two types of fermenter: continuous and batch fermenter. Chemostat is a continuous fermenter. In the continuous fermenter, the chamber of the fermenter is occupied with raw ingredients to be fermented and maintains the temperature and pH for microbes. Fermentation is appropriately settled, and nutrients are added. All the materials are well purified, and pure culture of bacteria is added to the fermenter from a separate vessel. Fermentation continues, and after a certain level of time, the contents of fermenter are taken out. In the fermentation process, the raw materials are continuously added to the fermenter at the same rate that keeps the microorganisms in the development phase, and hence, fermentation products are taken out continuously. An outflow valve allows the development medium to leave at a fix rate, and therefore, the culture’s volume needs to be constant.Hereafter, our aim is to develop a continuous-time mathematical model of the chemostat based on the following main assumptions [1]:(i)
The flow rate of the product should be slow, otherwise the cultivated population of bacteria will disappear.(ii)
The nutrient supply should be satisfactory for the proper development of bacteria because the growth rate of bacteria is directly proportional to nutrients.Since the amount of variables should be measured in keeping way of microorganism’s populace and its nourishment supply, we should be cautious for amassing conditions. Now in the following, we also keep track of assumptions made in the model:(1)
The development compartment is well mixed, and hence, it has no spatial verities in concentration of supplements or microorganisms. These phenomena can be described by following the first-order differential equation:(1)dNdt=KN−FN,where the rate of reproduction isK.(2)
The growth rate of the population depends on nutrient availability, so thatK=KC.Now, in order to write the mathematical equation for change inC, where C is the supplement level in development compartment. Here, again there are numerous influences tending to decrease or increase concentration: inflow of stock supply and depletion by bacteria, as well as outflow of nutrients in the effluent. Now assume that(3)
Supplement consumptions happen persistently as a consequence of reproduction. Hence, we indicate culture development imperatively, for which supplement reduction is basically reaching to be the same as recently.So, based on the previous assumption, the mathematical equation for supplements can be written as follows:(2)dCdt=−αKCN−FC+FC0.From equations (1) and (2), one gets the complete uncorrected, due to nondimensionless, two-dimensional chemostat model represented by the following system of differential equations:(3)dNdt=KCN−FN,dCdt=−αKCN−FC+FC0,whereN,C, and involved chemostat parameters along with corresponding dimensions in model (3) are expressed in the form of Table 1 (see [1]).Hereafter, in order to get the chemostat model with dimensionless form, which is depicted in model (3), one divides FN,FC, and FC0 by V, and hence, model (3) finally takes the following form:(4)dNdt=KCN−FNV,dCdt=−αKCN−FCV+FC0V.In the model similar to chemostat in which substances are transferred from one section to another, an excellent beginning point to write the mathematical equation is the physical rule, where massM is preserved. An identical preservation explanation is that amount of central is preserved. So, model (4) after balancing the mass of species N and C becomes as follows:(5)dNVdt=KCNV−FN,dCVdt=−αKCNV−FC+FC0,whereV denotes the number of microscopic organisms in a compartment and CV represents the mass supplement within the compartment.Hereafter, we will introduce the effect of bacterial growth to model (5), which may depend on supplement availability [1]. It is anticipated that for low nutrient abundance, growth rate typically increases with increasing nutrient concentrations. Eventually, when an excess of nutrient is available, its uptake rate and the resultant reproductive rate of the organisms do not continue to increase indefinitely.Suitable presumptions, in this way, would be ones that join the impact of a soaking reliance due to the following further assumption:(4)
The growth rate increases with supplement availability up to some restricting value.One sort of component that consolidates this impact is that the bacterial development rate and supplement utilizationsKC are expected to be an immersing work of supplement concentration called Michaelis Menten Kinetics, which is given as follows:(6)KC=KmaxCKn+C,where Kmax represents an upper bound for KC. Also, Kmax and Kn depend on the type of bacteria and supplement medium. Now, utilizing equation (6) in (5), one gets the following equation:(7)dNdt=KmaxCKn+CN−FNV,dCdt=−αKmaxCKn+CN−FCV+FC0V.Furthermore, the chemostat model (7) becomes as follows:(8)dN∗N∧dt∗τ=KmaxC∗C∧Kn+C∗C∧N∗N∧−FN∗N∧V,dC∗C∧dt∗τ=−αKmaxC∗C∧Kn+C∗C∧N∗N∧−FC∗C∧V+FC0V,by using the following dimensionless transformations:(9)N=N∗×N∧,C=C∗×C∧,t=t∗×τ,where time independent constants are N∧,C∧,andτ and N∗,C∗,andt∗ has no dimensions. Now, multiplying equation (8) by τ and then dividing the resulting equation by N∧ and C∧, respectively, and finally separating the constant terms, one gets the following equation:(10)dN∗dt∗=τKmaxC∗Kn/C∧+C∗N∗−τFN∗V,dC∗dt∗=−ατKmaxN∧C∗C∗Kn/C∧+C∗N∗−τFC∗V+τFC0VC∧.Now, the chemostat model that is presented in equation (10) becomes(11)dNdt=VFKmaxC1+CN−N,dCdt=−C1+CN−C+C0Kn,by using the following transformations:(12)τ=VF,C∧=Kn,N∧=KnατKmax,and dropping the stars. Now, if one denotes α1=V/FKmax and α2=C0/Kn, then from equation (11), the final form of the two-dimensional chemostat model represents the following system of differential equations:(13)dNdt=α1C1+CN−N,dCdt=−C1+CN−C+α2.Finally, by using the Euler-forward formula, the discrete form of the chemostat model, which is depicted in equation (13), is as follows:(14)Nt+1=1−hNt+α1hCt1+CtNt,Ct+1=α2h+1−hCt−hCt1+CtNt.Table 1
Biological interpretations of chemostat parameters along with the corresponding dimensions.
Chemostat parametersCorresponding dimensionsBiological interpretationsNNo/VBacterial populace thicknessCM/VSupplement concentration in the development compartmentC0M/VSupplement concentration in the storeVVVolume of the development chamberFV/TOutput/intake stream rate


### 1.2. Literature Review, Problem Statement, and Paper Layout

In recent years, many scholars investigated stability and bifurcations of mathematical models. For instance, Ruan and Wolkowicz [2] investigated the Hopf bifurcation analysis of a chemostat model with a distributed delay represented by the integro-differential equations. Yuan and Han [3] investigated dynamical behavior of a chemostat model with two distributed delays. Xu et al. [4] studied the bifurcation analysis of a food chain chemostat model with Michaelis–Menten functional response and double delays. Ma et al. [5] investigated the bifurcation for a fractional-order chemostat model with time delays and incommensurate orders. Rana [6] investigated chaotic dynamics of a discrete predator-prey system. Khan [7] studied bifurcations of a discrete predator-prey model. Mehrjooee et al. [8] studied extended bifurcation analysis of a flexible-link manipulator. Chakraborty et al. [9] investigated stability analysis of a prey-predator model with optimal harvesting and square-root functional response. Liu and Cai [10] studied dynamic behavior of a predator-prey system. Chen [11] studied global stability and permanence of multispecies Lotka–Volterra systems. Fang et al. [12] investigated dynamics of a discrete system with Beddington–DeAngelis function response. By using the parallel fully-implicit ocean model, Thies et al. [13] investigated the bifurcation of 3D ocean flows. Beckmann and Schröter [14] developed mathematical models of circulations in seas and oceans. de Niet et al. [15] studied bifurcation analysis of ocean models. For more articles in this direction, we refer the interesting readers to [16–19] and references cited therein. On the other hand, it is important here to mention that we can also apply the theory of the discrete dynamical system to other recent published works [20–28] to explore the possibilities of bifurcation and chaos.Motivated from aforementioned studies, the purpose of present studies is to investigate the dynamics of a discrete chemostat model (14). More specifically, our key contribution in the article includes(i)
To explore the local dynamical classifications with topological classification at semitrivial and interior equilibrium solutions(ii)
Study of the existence of periodic points for the understudied discrete chemostat model(iii)
Existence of possible bifurcations by bifurcation theory(iv)
Investigation of flip bifurcation at semitrivial equilibrium solutions(v)
Detailed study of flip bifurcation at the interior equilibrium solution, and we proved that no other bifurcation exists at the interior equilibrium solution(vi)
Verification of theoretical results numericallyThe rest of the paper is organized as follows: in Section2, we study equilibrium solutions and linearized form of the discrete chemostat model (14), whereas Section 3 is about the local dynamic behavior at equilibrium solutions. In Section 4, we study the existence of periodic points of the discrete chemostat model (14), whereas brief bifurcation analysis is given in Section 5. Further chaos control by the state feedback method is studied in Section 6. Section 7 contains several numerical simulations, whereas conclusions are given in Section 8.


## 1.1. Motivation and Mathematical Modelling of the Discrete Chemostat Model

In a biological laboratory, chemostat is a basic device which also occupies a central role in the field of mathematics. It is a model of simple lake to study exploitation competition. This model is also used to tract the water wastes process. Oceans are large water bodies which contain major fragments of life form on the Earth including microorganisms. The earliest life form on the planet Earth first appeared in oceans which was a microorganism. Similarly, the first life on the Earth was also migrated from oceans. Chemostat is a device that is used to study physiological interaction of microorganisms under the conditions when nutrients become limited and observed in oligotrophic oceans. These microorganisms can be used to study carbonate chemistry. The effect of accretion on the phytoplankton’s growth rate and general physiology of phytoplankton species should be assessed in the experiment. Chemostat is used for culturing these cyanobacteria to access carbonate chemistry of oceans, effect of the flow rate, temperature,pco2, nutrient concentration, and composition. A few chemostat investigations have been carried out at sea to study the effects of constant low supply of the limiting nutrient iron on the survival and growth of phytoplankton.Chemostat provides most favorable environment for microorganisms to develop into a biomass. That is why, they can interact with raw materials forming a product. It is noted that fermentation takes place in a fermenter. There are two types of fermenter: continuous and batch fermenter. Chemostat is a continuous fermenter. In the continuous fermenter, the chamber of the fermenter is occupied with raw ingredients to be fermented and maintains the temperature and pH for microbes. Fermentation is appropriately settled, and nutrients are added. All the materials are well purified, and pure culture of bacteria is added to the fermenter from a separate vessel. Fermentation continues, and after a certain level of time, the contents of fermenter are taken out. In the fermentation process, the raw materials are continuously added to the fermenter at the same rate that keeps the microorganisms in the development phase, and hence, fermentation products are taken out continuously. An outflow valve allows the development medium to leave at a fix rate, and therefore, the culture’s volume needs to be constant.Hereafter, our aim is to develop a continuous-time mathematical model of the chemostat based on the following main assumptions [1]:(i)
The flow rate of the product should be slow, otherwise the cultivated population of bacteria will disappear.(ii)
The nutrient supply should be satisfactory for the proper development of bacteria because the growth rate of bacteria is directly proportional to nutrients.Since the amount of variables should be measured in keeping way of microorganism’s populace and its nourishment supply, we should be cautious for amassing conditions. Now in the following, we also keep track of assumptions made in the model:(1)
The development compartment is well mixed, and hence, it has no spatial verities in concentration of supplements or microorganisms. These phenomena can be described by following the first-order differential equation:(1)dNdt=KN−FN,where the rate of reproduction isK.(2)
The growth rate of the population depends on nutrient availability, so thatK=KC.Now, in order to write the mathematical equation for change inC, where C is the supplement level in development compartment. Here, again there are numerous influences tending to decrease or increase concentration: inflow of stock supply and depletion by bacteria, as well as outflow of nutrients in the effluent. Now assume that(3)
Supplement consumptions happen persistently as a consequence of reproduction. Hence, we indicate culture development imperatively, for which supplement reduction is basically reaching to be the same as recently.So, based on the previous assumption, the mathematical equation for supplements can be written as follows:(2)dCdt=−αKCN−FC+FC0.From equations (1) and (2), one gets the complete uncorrected, due to nondimensionless, two-dimensional chemostat model represented by the following system of differential equations:(3)dNdt=KCN−FN,dCdt=−αKCN−FC+FC0,whereN,C, and involved chemostat parameters along with corresponding dimensions in model (3) are expressed in the form of Table 1 (see [1]).Hereafter, in order to get the chemostat model with dimensionless form, which is depicted in model (3), one divides FN,FC, and FC0 by V, and hence, model (3) finally takes the following form:(4)dNdt=KCN−FNV,dCdt=−αKCN−FCV+FC0V.In the model similar to chemostat in which substances are transferred from one section to another, an excellent beginning point to write the mathematical equation is the physical rule, where massM is preserved. An identical preservation explanation is that amount of central is preserved. So, model (4) after balancing the mass of species N and C becomes as follows:(5)dNVdt=KCNV−FN,dCVdt=−αKCNV−FC+FC0,whereV denotes the number of microscopic organisms in a compartment and CV represents the mass supplement within the compartment.Hereafter, we will introduce the effect of bacterial growth to model (5), which may depend on supplement availability [1]. It is anticipated that for low nutrient abundance, growth rate typically increases with increasing nutrient concentrations. Eventually, when an excess of nutrient is available, its uptake rate and the resultant reproductive rate of the organisms do not continue to increase indefinitely.Suitable presumptions, in this way, would be ones that join the impact of a soaking reliance due to the following further assumption:(4)
The growth rate increases with supplement availability up to some restricting value.One sort of component that consolidates this impact is that the bacterial development rate and supplement utilizationsKC are expected to be an immersing work of supplement concentration called Michaelis Menten Kinetics, which is given as follows:(6)KC=KmaxCKn+C,where Kmax represents an upper bound for KC. Also, Kmax and Kn depend on the type of bacteria and supplement medium. Now, utilizing equation (6) in (5), one gets the following equation:(7)dNdt=KmaxCKn+CN−FNV,dCdt=−αKmaxCKn+CN−FCV+FC0V.Furthermore, the chemostat model (7) becomes as follows:(8)dN∗N∧dt∗τ=KmaxC∗C∧Kn+C∗C∧N∗N∧−FN∗N∧V,dC∗C∧dt∗τ=−αKmaxC∗C∧Kn+C∗C∧N∗N∧−FC∗C∧V+FC0V,by using the following dimensionless transformations:(9)N=N∗×N∧,C=C∗×C∧,t=t∗×τ,where time independent constants are N∧,C∧,andτ and N∗,C∗,andt∗ has no dimensions. Now, multiplying equation (8) by τ and then dividing the resulting equation by N∧ and C∧, respectively, and finally separating the constant terms, one gets the following equation:(10)dN∗dt∗=τKmaxC∗Kn/C∧+C∗N∗−τFN∗V,dC∗dt∗=−ατKmaxN∧C∗C∗Kn/C∧+C∗N∗−τFC∗V+τFC0VC∧.Now, the chemostat model that is presented in equation (10) becomes(11)dNdt=VFKmaxC1+CN−N,dCdt=−C1+CN−C+C0Kn,by using the following transformations:(12)τ=VF,C∧=Kn,N∧=KnατKmax,and dropping the stars. Now, if one denotes α1=V/FKmax and α2=C0/Kn, then from equation (11), the final form of the two-dimensional chemostat model represents the following system of differential equations:(13)dNdt=α1C1+CN−N,dCdt=−C1+CN−C+α2.Finally, by using the Euler-forward formula, the discrete form of the chemostat model, which is depicted in equation (13), is as follows:(14)Nt+1=1−hNt+α1hCt1+CtNt,Ct+1=α2h+1−hCt−hCt1+CtNt.Table 1
Biological interpretations of chemostat parameters along with the corresponding dimensions.
Chemostat parametersCorresponding dimensionsBiological interpretationsNNo/VBacterial populace thicknessCM/VSupplement concentration in the development compartmentC0M/VSupplement concentration in the storeVVVolume of the development chamberFV/TOutput/intake stream rate


## 1.2. Literature Review, Problem Statement, and Paper Layout

In recent years, many scholars investigated stability and bifurcations of mathematical models. For instance, Ruan and Wolkowicz [2] investigated the Hopf bifurcation analysis of a chemostat model with a distributed delay represented by the integro-differential equations. Yuan and Han [3] investigated dynamical behavior of a chemostat model with two distributed delays. Xu et al. [4] studied the bifurcation analysis of a food chain chemostat model with Michaelis–Menten functional response and double delays. Ma et al. [5] investigated the bifurcation for a fractional-order chemostat model with time delays and incommensurate orders. Rana [6] investigated chaotic dynamics of a discrete predator-prey system. Khan [7] studied bifurcations of a discrete predator-prey model. Mehrjooee et al. [8] studied extended bifurcation analysis of a flexible-link manipulator. Chakraborty et al. [9] investigated stability analysis of a prey-predator model with optimal harvesting and square-root functional response. Liu and Cai [10] studied dynamic behavior of a predator-prey system. Chen [11] studied global stability and permanence of multispecies Lotka–Volterra systems. Fang et al. [12] investigated dynamics of a discrete system with Beddington–DeAngelis function response. By using the parallel fully-implicit ocean model, Thies et al. [13] investigated the bifurcation of 3D ocean flows. Beckmann and Schröter [14] developed mathematical models of circulations in seas and oceans. de Niet et al. [15] studied bifurcation analysis of ocean models. For more articles in this direction, we refer the interesting readers to [16–19] and references cited therein. On the other hand, it is important here to mention that we can also apply the theory of the discrete dynamical system to other recent published works [20–28] to explore the possibilities of bifurcation and chaos.Motivated from aforementioned studies, the purpose of present studies is to investigate the dynamics of a discrete chemostat model (14). More specifically, our key contribution in the article includes(i)
To explore the local dynamical classifications with topological classification at semitrivial and interior equilibrium solutions(ii)
Study of the existence of periodic points for the understudied discrete chemostat model(iii)
Existence of possible bifurcations by bifurcation theory(iv)
Investigation of flip bifurcation at semitrivial equilibrium solutions(v)
Detailed study of flip bifurcation at the interior equilibrium solution, and we proved that no other bifurcation exists at the interior equilibrium solution(vi)
Verification of theoretical results numericallyThe rest of the paper is organized as follows: in Section2, we study equilibrium solutions and linearized form of the discrete chemostat model (14), whereas Section 3 is about the local dynamic behavior at equilibrium solutions. In Section 4, we study the existence of periodic points of the discrete chemostat model (14), whereas brief bifurcation analysis is given in Section 5. Further chaos control by the state feedback method is studied in Section 6. Section 7 contains several numerical simulations, whereas conclusions are given in Section 8.


## 2. Equilibrium Solutions and Linearized Form of the Discrete Chemostat Model (14)

The existence of equilibrium solutions and the linearized form of the discrete chemostat model (14) are discussed in this section. We first give the result of the existence of the equilibrium solution of the discrete chemostat model (14) in the region R+2=N,C:N,C≥0 as follows:Lemma 1.
For existence of equilibrium solutions of the discrete chemostat model (14), the following results hold:(i)
∀α1,α2,h, the discrete chemostat model (14), has the semitrivial equilibrium solution F0C0,α2;(ii)
Ifα2>1/α1−1 with α1>1, then the discrete chemostat model (14) has the interior equilibrium solution FNC+α1α2−1/α1−1,1/α1−1.Proof.
If the discrete chemostat model (14) has the equilibrium solution FNCN,C, then(15)N=1−hN+α1hC1+CN,C=α2h+1−hC−hC1+CN.
It should be noted here that equation (15) is satisfied obviously if N=0 and C=α2, and hence, we can conclude that ∀α1,α2,h, discrete chemostat model (14), has the semitrivial equilibrium solution F0C0,α2). Moreover, after some straightforward manipulation, from equation (15), one gets N=α1α2−1/α1−1 and C=1/α1−1. Therefore, we can also conclude that the discrete chemostat model (14) has the interior equilibrium solution FNC+α1α2−1/α1−1,1/α1−1 if α2>1/α1−1 with α1>1.
Now, one has the following map to construct the linearized form of equation (14) at FNCN,C:(16)f1,f2↦Nt+1,Ct+1,where(17)f1≔1−hN+α1hC1+CN,f2≔α2h+1−hC−hC1+CN.
Finally, variation matrixV|FNCN,C at FNCN,C under equation (16) along with equation (17) is as follows:(18)V|FNCN,C≔1−h+α1hC1+Cα1hN1+C2−hC1+C1−h−hN1+C2.


## 3. Local Dynamic Behavior at Equilibrium Solutions of the Chemostat Model (14)

The local dynamic behavior at equilibrium solutionsF0C0,α2 and FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) is deeply studied by the stability theory [29–33] in this section. It is noted that, at F0C0,α2, equation (18) becomes(19)V|F0C0,α2≔1−h+α1α2h1+α20−hα21+α21−h,with(20)λ1≔1−h+α1α2h1+α2,λ2≔1−h.Lemma 2.
The dynamic behavior atF0C0,α2 of the chemostat model (14) can be concluded as follows:(i)
F0C0,α2 of the chemostat model (14) is a sink if(21)α2>1α1−1withα1>1and0<h<2;(ii)
F0C0,α2 of the chemostat model (14) is a source if(22)0<α2<1α1−1withα1>1andh>2;(iii)
F0C0,α2 of the chemostat model (14) is a saddle if(23)α2>1α1−1withα1>1andh>2;(iv)
F0C0,α2 of the chemostat model (14) is nonhyperbolic if(24)h≔2orα2≔1α1−1.Now, atFNC+α1α2−1/α1−1,1/α1−1, equation (18) becomes(25)V|FNC+α1α2−1α1−1,1α1−1≔1hα1−1α2α1−1−1−hα11−h−α1−1hα2α1−1−1α1,whose characteristic equation is as follows:(26)λ2−pα1α2−1α1−1,1α1−1λ+qα1α2−1α1−1,1α1−1=0,where(27)pα1α2−1α1−1,1α1−1=2−h−hα1−1α2α1−1−1α1,qα1α2−1α1−1,1α1−1=1−h1−hα1−1α2α1−1−1α1.Moreover, from equation (27), one has the following equation:(28)λ1,2≔pα1α2−1/α1−1,1/α1−1±Δ2,where(29)Δ=pα1α2−1α1−1,1α1−12−4qα1α2−1α1−1,1α1−1=2−h−hα1−1α2α1−1−1α12−41−h1−hα1−1α2α1−1−1α1=hα11+α2+α12α2−2α11+α22>0.From equations (28) and (29), the straightforward computation yields as follows:(30)λ1≔α1−hα12α2+2hα1α2−hα2−h+hα1α1,λ2≔1−h.Therefore, based on equation (30), the dynamic behavior at FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) can be concluded as follows:Lemma 3.
ForFNC+α1α2−1/α1−1,1/α1−1 of the discrete chemostat model (14), the following expression holds:(i)
FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) is a stable node if(31)0<α2<min2α1+hα1−hhα1−12,1α1−1and0<h<2;(ii)
FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) is an unstable node if(32)α2>2α1+hα1−hhα1−12andh>2;(iii)
FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) is a saddle node if(33)0<α2<min2α1+hα1−hhα1−12,1α1−1andh>2orα2>2α1+hα1−hhα1−12and0<h<2;(iv)
FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) is nonhyperbolic if(34)α2≔2α1+hα1−hhα1−12orh≔2.


## 4. Chemostat Model (14) with Periodic Points

Now, it is investigated thatF0C0,α2 and FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) are periodic points of period-1,2,…,n by the existing theory [34].Theorem 1.
F0C0,α2 and FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) are periodic points of prime period-1.Proof.
From equation (14), one has the following equation:(35)Φ≔f1,f2,where f1 and f2 are depicted in equation (17). The following computation yields the proof of the required results:(36)Φ|F0C0,α2=F0C0,α2,Φ|FNC+α1α2−1α1−1,1α1−1=FNC+α1α2−1α1−1,1α1−1.Theorem 2.
F0C0,α2 and FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) are periodic points of period-2,3,…,n.Proof.
From equation (42), one has the following equation:(37)Φ2=1−hf1+α1hf21+f2f1,α2h+1−hf2−hf21+f2f1⇒Φ2|F0C0,α2=F0C0,α2,Φ3=1−hf12+α1hf221+f22f12,α2h+1−hf22−hf221+f22f12⇒Φ3|F0C0,α2=F0C0,α2,⋮Φn=1−hf1n−1+α1hf2n−11+f2n−1f1n−1,α2h+1−hf2n−1−hf2n−11+f2n−1f1n−1⇒Φn|F0C0,α2=F0C0,α2,which shows that F0C0,α2 of the chemostat model (14) are periodic points of period-2,3,…,n. Similarly, from equation (44), one has the following equation:(38)ΦFNC+α1α2−1/α1−1,1/α1−12=FNC+α1α2−1α1−1,1α1−1,ΦFNC+α1α2−1/α1−1,1/α1−13=FNC+α1α2−1α1−1,1α1−1,⋮ΦFNC+α1α2−1/α1−1,1/α1−1n=FNC+α1α2−1α1−1,1α1−1,which shows that FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) are periodic points of period-2,3,…,n.


## 5. Bifurcation

Here, based on our theoretical findings in Section3, comprehensive bifurcation analysis at equilibrium solutions F0C0,α2 and FNC+α1α2−1/α1−1,1/α1−1 of the chemostat model (14) will be investigated by the bifurcation theory [35, 36]. For this, if equation (24) holds, then from equation (20), one gets λ2|24=−1, but λ1|24=−1−2α1α2/1+α2≠1or−1, which implies that the discrete chemostat model (14) may undergo flip bifurcation if α1,α2,h is located in the following set:(39)F|F0C0,α2≔α1,α2,h,h=2.Theorem 3.
Around the semitrivial equilibrium solutionF0C0,α2 of the discrete chemostat model (14), no flip bifurcation exists if α1,α2,h∈F|F0C0,α2.Proof.
Since the chemostat model (14) is an invariant subject to N=0, so one restricts it on N=0, where model (14) takes the following form:(40)Ct+1=α2h+1−hCt.
From equation (40), we have the map as follows:(41)fh,C≔α2h+1−hC.
Finally, ifh∗=h=2 and C=C∗=α2, then from equation (41), we get the following equations:(42)∂f∂Ch∗=2,C∗=α2≔−1,(43)∂2f∂C2h∗=2,C∗=α2≔0,and(44)∂f∂hh∗=2,C=C∗=α2≔0.
Equations (43) and (44) violate the nondegenerate condition and so if α1,α2,h∈F|F0C0,α2, then around the semitrivial equilibrium solution F0C0,α2, there exists no flip bifurcation.
Now, if equation (40) holds, then from equation (31), one gets λ1|40=−1, but λ2|40=1−h≠1or−1. This implies that, at FNC+α1α2−1/α1−1,1/α1−1 of the discrete chemostat model (14), there may exist flip bifurcation if α1,α2,h is located in the following set:(45)F|FNC+α1α2−1/α1−1,1/α1−1≔α1,α2,h,α2=2α1+hα1−1hα1−12.
In the following theorem, we can conclude that ifα1,α2,h∈F|FNC+α1α2−1/α1−1,1/α1−1, then flip bifurcation must exist for the discrete chemostat model (14).Theorem 4.
Around the interior equilibrium solutionFNC+α1α2−1/α1−1,1/α1−1 of the discrete chemostat model (14), flip bifurcation must exist if α1,α2,h∈F|FNC+α1α2−1/α1−1,1/α1−1.Proof.
Sinceα1,α2,h∈F|FNC+α1α2−1/α1−1,1/α1−1, thus parameter α2 is clearly the bifurcation parameter. Furthermore, the discrete chemostat model (14) takes the following form:(46)Nt+1=1−hNt+α1hCt1+CtNt,Ct+1=α2∗+ϵh+1−hCt−hCt1+CtNt,where α2=α2∗+ϵ and ϵ≪1. Now, the following transformations are utilized in order to transform FNC+α1α2−1/α1−1,1/α1−1 of the discrete chemostat model (46) to F000,0:(47)ut=Nt−N∗,vt=Ct−C∗,with N∗=α1α2−1/α1−1 and C∗=1/α1−1. From equations (47) and (48), we have the following equation:(48)ut+1=ut−hut+N∗+α1hvt+C∗1+vt+C∗ut+N∗,vt+1=α2∗+ϵh+vt−hvt+C∗−hvt+C∗1+vt+C∗ut+N∗.
The Taylor series expansion up to3rd-order gives the following equation:(49)ut+1=κ11^ut+κ12^vt+κ13^utvt+κ14^vt2+κ15^utvt2+κ16^vt3,vt+1=κ21^ut+κ22^vt+κ23^utvt+κ24^vt2+κ25^utvt2+κ26^vt3,where(50)κ11^=1−h+hC∗α11+C∗,κ12^=hN∗α11+C∗2,κ13^=hα11+C∗2,κ14^=−hN∗α11+C∗3,κ15^=−hα11+C∗3,κ16^=hN∗α11+C∗4,κ21^=−hC∗1+C∗,κ22^=1−h+hN∗C∗1+C∗2−hN∗1+C∗,κ23^=−h1+C∗2,κ24^=hN∗1+C∗3,κ25^=h1+C∗3,κ26^=−hN∗1+C∗4.
Now, equation (50) takes the following form:(51)Nt+1Ct+1=−100λ2NtCt+P^ut,vt,ϵQ^ut,vt,ϵ,where(52)P^=hκ13^+2α1/hκ23^α12−hutvt+hκ14^+2α1/hκ24^α12−hvt2+hκ15^+2α1/hα12−hutvt2+hκ16^+2α1/hα12−hvt3,Q^=−hκ13^+α1κ23^α12−hutvt−hκ14^+α1κ24^α12−hvt2−hκ15^+α1κ25^α12−hutvt2−hκ16^+α1κ26^α12−hvt3,ut=−α1Nt−2α1hCt,vt=Nt+Ct,vt2=Nt2+Ct2+2NtCt,utvt=−α1Nt2−α11+2hNtCt−2α1hCt2,utvt2=−α1Nt3−α11+4hNtCt2−2α11+1hNt2Ct−2α1hCt3,vt3=Nt3+3Nt2Ct+3NtCt2+Ct3,by(53)utvt≔−α1−2α1h11NtCt.
Next, for equation (53), one can write the center manifold as follows:(54)McF000,0=Nt,Ct:Ct=s0ϵ+s1Ct2+s2Ntϵ+s3ϵ3+ONt+ϵ3,where the computation yields(55)s0=0,s1=12−hκ13^+α1κ23^−hκ14^−α1κ24^α1+hκ15^+α1κ25^,s2=0=s3.
Lately, in the following equation, we express equation (51) restricted to McF000,0 as follows:(56)fNt=−Nt+h1Nt2+h2Ntϵ+h3Nt2ϵ+h4Ntϵ2+h5Nt3+ONt+ϵ4,where(57)h1=h2−h−κ13^+2α1hκ23^+κ14^+2α1/hκ24^α1−κ15^+2α1h,h2=−h2−hκ13^+2α1hκ23^1+2h,h3=0=h4,h5=12−hs12hκ14^+2α1/hκ24^α1−hκ13^+2α1hκ23^1+2h+hκ16^+2α1/hα1.
For map (56) to undergo flip bifurcation, the following discriminatory quantities are nonzero:(58)Γ1=∂2f∂Nt∂ϵ+12∂f∂ϵ∂2f∂Nt2|F000,0,Γ2=16∂3f∂Nt3+12∂2f∂Nt22|F000,0.
After calculating, one gets the following equations:(59)Γ1=h2≠0,and(60)Γ2=h1+h52≠0.
From equation (60), it should be noted that if Γ2≠0 as α1,α2,h∈F|FNC+α1α2−1/α1−1,1/α1−1, then at FNC+α1α2−1/α1−1,1/α1−1, the discrete chemostat model (14) undergoes flip bifurcation, which further implies that period-2 points bifurcate from FNC+α1α2−1/α1−1,1/α1−1are stable (respectively, unstable) if Γ2>0 (respectively Γ2<0.Remark 1.
In summary, regarding bifurcation analysis for the boundary and interior fixed points of a discrete chemostat model (14), we proved that at the semitrivial equilibrium solution F0C0,α2, the eigenvalues criterion for the existence of flip bifurcation holds, but by computation, the flip bifurcation do not exist if α1,α2,h∈F|F0C0,α2. On the other hand, we also proved that only flip bifurcation occurs at FNC+α1α2−1/α1−1,1/α1−1 when parameters α1,α2,h∈F|FNC+α1α2−1/α1−1,1/α1−1, and no other bifurcation exists at this point.


## 6. Chaos Control

The state feedback control method is used to control chaos motivated from the existing literature [37, 38]. For this, we expressed the controlled system as follows:(61)Nt+1=1−hNt+α1hCt1+CtNt+ut,Ct+1=α2h+1−hCt−hCt1+CtNt,where ut=−k1Nt−N−k2Ct−C, k1,k2 denotes feedback gains and N=α1α2−1/α1−1 and Z=1/α1−1. Moreover, for control system (61), Vcon|FNCN,C becomes(62)Vcon|PNCN,C=l11−k1l12−k2l21l22,where(63)l11=1,l12=hα1−1α1α2−α2−1,l21=−hα1,l22=1−h−hα1−1α1α2−α2−1α1.Now, ifλ1,2 are characteristic roots of Vcon|PNCN,C, then(64)λ1+λ2=l11+l22−k1,(65)λ1λ2=l22l11−k1−l21l12−k2.Lemma 4.
If both eigenvalues ofVcon|PNCN,C satisfy λ1,2<1, then the controlled system (61) is asymptotically stable.Proof.
It should be noted here that marginal stability can be obtained from the conditionsλ1=±1,λ1λ2=1. If λ1λ2=1, then from equation (65), one gets the following equation:(66)L1:1−hα1−hα1−1α1α2−α2−1k1+hk2+hα1−1α1α2−α2−11−h+hα1=0.
Ifλ1=1, then from equations (65) and (66), one gets the following equation:(67)L2:hα1+α1−1α1α2−α2−1k1−hk2−h2α1−1α1α2−α2−1=0.
Finally, ifλ1=−1, then from equations (65) and (66), one gets the following equation:(68)L3:2−hα1−hα1−1α1α2−α2−1k1+hk2−2−h2α1−hα1−1α1α2−α2−1=0.
Therefore, from equations (66)–(68), lines L1,L2, and L3 gives the conditions for the eigenvalues satisfying λ1,2<1. Moreover, the triangular region bounded by L1,L2, and L3 accommodates stable eigenvalues.


## 7. Numerical Simulations

In this section, to verify our theoretical results and reveal new dynamical behaviors of a chemostat model (14), we make the phase portraits, bifurcation diagrams, and maximal Lyapunov exponents for specific parameter values by using MATLAB software. For the completeness of this section, the following three cases are to be considered:Case I: Ifα1=1.23, then 1/α1−1=4.347826086956522, and moreover, if one chooses α2 and h in such a way that satisfied the conclusion i of Lemma 2, that is, α2=4.56>1/α1−1=4.347826086956522 with h=1.9<2. Therefore, in this case, the semitrivial equilibrium solution F0C0,α2=0,4.56 of the chemostat model (14) is a sink whose numerical plots, t vs Nt and Ct, with initial values (0.9, 1.2) are shown in Figures 1(a) and 1(b). But if one chooses α2=4.1<1/α1−1=4.347826086956522 with h=2.1>2, then clearly condition ii of Lemma 2 holds, and so the semitrivial equilibrium solution F0C0,α2=0,4.1 of the chemostat model (14) is a source (Figures 1(c) and 1(d)). Therefore, the simulation agrees with the conclusion of Lemma 2.Case II: Ifα1=3.9,h=1.2<2, then from equation (40), one gets α2=1.117717003567182, and so from equation (32), the fixed point FNC+α1α2−1/α1−1,1/α1−1=FNC+3.0231724137931035,0.3448275862068966 of the chemostat model (14) is a stable node if 0<α2<1.117717003567182,0.3448275862068966, exchange stability if α2=1.117717003567182, and becomes an unstable node if α2>1.117717003567182 with h>2. As a consequence, flip bifurcation occurs if α2=1.117717003567182. For instance, if α2=1.12>1.117717003567182, then from equations (30) and (59), one gets λ1=−1,λ2=−0.2≠1or−1 and Γ1=6.900512820512821≠0. Moreover, from equation (67), one gets Γ2=1645.8589762806776>0, which indicates that stable period-2 point bifurcates from FNC+α1α2−1/α1−1,1/α1−1, and hence, maximum Lyapunov exponents and flip bifurcation diagram are plotted and presented in Figure 2. Moreover, corresponding to Figure 2, 3D flip bifurcation diagrams are plotted in Figure 3. Moreover, Figure 4 gives the complex dynamics of the discrete chemostat model (14) with orbits of period −7,−8,−13, and −16. Finally, for this case, the plots of t vs Nt and Ct are also presented in Figure 5.Case III: Finally, the result of Lemma4 will be verified. For instance, if h=1.8,α1=3.9,andα2=1.22, then from equations (66)–(68), one gets the following equations:(69)L1:−5.468999999999999k1+1.8k2+5.1408=0,(70)L2:9.369k1−1.8k2−4.2282=0,(71)L3:−1.569k1+1.8k2−1.0901999999999998=0.So, the equations (69)–(71) determine the triangular region that gives λ1,2<1 (Figure 6).Figure 1
Phase portraits for the boundary equilibrium solutionF0C0,α2.
(a)(b)(c)(d)Figure 2
(a and b) Flip bifurcation diagrams of the discrete chemostat model (14) with α2∈0.5,1.5. (c) Maximum Lyapunov exponents corresponding to (a and b) with (0.4, 0.3).
(a)(b)(c)Figure 3
3D flip bifurcation diagrams.
(a)(b)(c)(d)(e)Figure 4
Complex dynamics of the discrete chemostat model (14). (a) α2 = 1.13 with (0.4, 0.3), (b) α2 = 1.19 with (0.4, 0.3), (c) α2 = 1.2 with (0.99, 2.6), and (d) α2 = 1.221 with (0.99, 2.6).
(a)(b)(c)(d)Figure 5
Phase portraits for the interior equilibrium solutionFNC+α1α2−1/α1−1,1/α1−1.
(a)(b)Figure 6
Stability region whereλ1,2<1.


## 8. Conclusion

In this work, we have investigated feasible equilibrium solutions, local stability at equilibrium solutions, periodic points, bifurcation, and controlling chaos by the state feedback control method of the discrete chemostat model (14). More precisely, it is investigated that the discrete chemostat model (14) has two equilibrium solutions F0C0,α2 and FNC+α1α2−1/α1−1,1/α1−1. Biologically, the equilibrium solution F0C0,α2 is not an area of interest for the experimentalists: no bacteria are left, and the concentration of the nutrient is the same as the stock solution. On the other hand, the equilibrium solution FNC+α1α2−1/α1−1,1/α1−1 does not always exist biologically, and it looks more inspiring. It depends upon the magnitude of the quantities α1 and α2. If α1<1, then one gets negative values. But biologically negative values would be meaningless, and so population densities and concentrations must always be positive. That is why, the quantities α1 and α2 must be chosen in such a way where α2>1/α1−1 with α1>1. Therefore, mathematically, we can conclude that ∀α1,α2,h, solution F0C0,α2, is designated as the semitrivial equilibrium solution of the discrete chemostat model (14), but if α2>1/α1−1 with α1>1, then it has interior equilibrium solution FNC+α1α2−1/α1−1,1/α1−1. By the stability theory of the discrete dynamical system, it is proved that F0C0,α2 of the chemostat model (14) is a sink if α2>1/α1−1withα1>1/and0<h<2, source if 0<α2<1/α1−1withα1>1andh>2, saddle if α2>1/α1−1withα1>1andh>2, and nonhyperbolic if h=2 or α2=1/α1−1; interior equilibrium solution FNC+α1α2−1/α1−1,1/α1−1 is a stable node if 0<α2<min2α1+hα1−h/hα1−12,1/α1−1 and 0<h<2, unstable node if α2>2α1+hα1−1/hα1−12 and h>2; saddle node if 0<α2<min2α1+hα1−1/hα1−12,1/α1−1 and h>2 or α2>2α1+hα1−1/hα1−12 and h>2 or α2>2α1+hα1−1/hα1−12 and 0<h<2, and nonhyperbolic if α2=2α1+hα1−1/hα1−12 or h=2. It is also explored the existence of periodic points for discrete chemostat model (14). Furthermore, by bifurcation theory, it is proved that at equilibrium solution F0C0,α2, the discrete chemostat model (14) does not undergo flip bifurcation if α1,α2,h∈F|F0C0,α2=α1,α2,h,h=2, but at FNC+α1α2−1/α1−1,1/α1−1, it undergoes flip bifurcation if α1,α2,h∈F|FNC+α1α2−1/α1−1,1/α1−1\\=α1,α2,h,α2=2α1+hα1−1/hα1−12, and no other bifurcation occurs at this point. Furthermore, we have controlled the chaos by the state feedback control method. Finally, with the help of MATLAB software, flip bifurcation diagrams and corresponding maximum Lyapunov exponents are also plotted. These numerical simulations not only show the correctness of theoretical results but also reveal the complex dynamics such as the orbits of the period −7,−8,−13, and −16.



---

*Source: 7518261-2023-02-18.xml*