Dynamics and Hopf Bifurcation Analysis of a Chemostat Model with Modified Growth Rate and Variable Yield Coefficient

Md. Shariful Islam¹, Touhid Hossain¹, Mir Shariful Islam^2,3,4

¹Department of Mathematics, University of Dhaka, Dhaka, Bangladesh.
²Mechanical Engineering Department, Kyushu University, Fukuoka, Japan.
³International Institute for Carbon-Neutral Energy Research, Kyushu University, Fukuoka, Japan.
⁴Department of Oceanography, University of Dhaka, Dhaka, Bangladesh.

DOI: 10.4236/ojmsi.2022.104023                

Abstract

The objective of this study is to analyze a chemostat model of very simple type with the Haldane expression of growth rate and a variable yield coefficient. The proposed modified model is analyzed qualitatively and quantitatively. Analytic conditions for stability and optimality are determined for washout and no washout equilibrium solutions. One of the main focuses of the study is to determine parameter values for which Hopf Bifurcations occur in a bioreactor. It has been shown that the maximum stable non-washout equilibrium exits at a residence time under suitable parameter values. Hopf bifurcation is observed at three different conditions of the parameters.

Keywords

Chemostat, Residence Time, Hopf Bifurcation, Bioreactor, Growth Rate, Haldane Model, Yield Coefficient

1. Introduction

Researchers have been studying extensively to improve product yield in chemical reactors or bio-reactors for many years. To gain desired scale of production, it is necessary to find effective engineering instruments and mechanism in a bio-reactor. One such engineering instrument is application of oscillatory external force. Since implementation of oscillatory external force is costly, it has very limited use in industry [1] [2]. Many researchers have studied the oscillatory operations in more than one reactors combined in series [3] [4] [5] [6] [7]. In an oscillatory operation, the values of the parameters in a model are determined for which a steady concentration of substrate into one reactor stimulates the successive reactors periodically. This automated periodic forcing of a reactor requires no extra energy or logistics. This is why it is very important to determine the dynamic conditions for which natural oscillations occur in a bio-reactor. So that it can be used to force another bio-reactor periodically without any external forcing [3] [6]. Since we would like to use natural oscillations of a reactor to force another one arranged in series periodically, it is practical to determine sufficient conditions for which natural oscillations occur in the first reactor. Nelson and Sidhu [8] believed that it is better to determine a benchmark for the maximum output from a single reactor before investigating a two-reactor system. Balakrishnan and Yang [9] examined a single bio-reactor model by direct integration on limited range of parameter values. On the other hand, Nelson and Sidhu [8] studied the same model analytically by exploiting theory of bifurcation and path following methods.

2. Model Equations

Our aim is to investigate a microbial system, the process of using systems biology to understand microbes and their environment, where growth of the cell mass (X) depends on substrate species (S). In this study, the specific growth rate is taken from the Haldane kinetic model with variable yield coefficient. We are studying a simple chemostat model to optimize the cell mass concentration as a function of residence time in a reactor and similar kind of study of microbial systems can be found in [1] [3] [6] [9].

2.1. Dimensional Model Description

The microbial system which we consider in this study is represented by the following Equations (1)-(4)

$V\frac{\text{d}S}{\text{d}t}=F\left(S_0-S\right)-VX\frac{\mu \left(S\right)}{Y\left(S\right)}$ (1)

$V\frac{\text{d}X}{\text{d}t}=F\left(X_0-X\right)+VX\mu \left(S\right)$ (2)

with the specific growth rate

$\mu \left(S\right)=\frac{{\mu }_{m}S\left(t\right)}{K_s+S+\frac{S^2}{K_i}}$ (3)

and a Monod expression with variable Yield Coefficients

$Y\left(S\right)=\alpha +\beta S,\left(\alpha ,\beta >0\right)$ (4)

where F is the flow rate (l·hr⁻¹), K_s is the half-saturation constant (g·l⁻¹), K_i is the inhibitory constant (g·l⁻¹), V is the volume (l), X₀ is the initial cell mass concentration (g·l⁻¹), S₀ is the initial substrate concentration (g·l⁻¹), ${\mu }_m$ is the maximum specific growth rate (/h), α(-) and β(g·l⁻¹) are constant yield coefficients.

2.2. Dimensionless Model

To make the above model dimensionless, we use the transformations $S^{*}=\frac{S}{K_s}$ , $X^{*}=\frac{X}{\alpha K_s}$ , and $t^{*}={\mu }_{m}t$ and we get the following dimensionless model containing the parameters $S_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\tau }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}$ , and ${\gamma }^{\mathrm{<mo>\; *\; </mo>}}$ .

$\frac{\text{d}{S}^{*}}{\text{d}{t}^{*}}=\frac{S_0^{*}-S^{*}}{{\tau }^{*}}-\frac{S^{*}{X}^{*}}{\left(1+S^{*}+{\gamma }^{*}{S}^{*2}\right)\left(1+{\beta }^{*}{S}^{*}\right)}$ (5)

$\frac{\text{d}{X}^{*}}{\text{d}{t}^{*}}=\frac{X_0^{*}-X^{*}}{{\tau }^{*}}+\frac{S^{*}{X}^{*}}{1+S^{*}+{\gamma }^{*}{S}^{*2}}$ (6)

where ${\tau }^{*}=\frac{V{\mu }_m}{F},{\beta }^{*}=\frac{\beta K_s}{\alpha }$ , and ${\gamma }^{*}=\frac{K_s}{K_i}$ .

We consider a chemostat setup consists of a sterile feed means that $X_0^{\mathrm{<mo>\; *\; </mo>}}$ will be zero and assume ${\tau }^{\mathrm{<mo>\; *\; </mo>}}$ (residence time) as our primary bifurcation parameter. Our secondary bifurcation parameters are $S_0^{\mathrm{<mo>\; *\; </mo>}}$ (substrate concentration in the feed), ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ (dimensionless yield coefficient) which depends on specific microbial system and hence not so flexible for tuning, and ${\gamma }^{\mathrm{<mo>\; *\; </mo>}}$ (dimensionless inhibitory constant). We will denote the dimensional and dimensionless variables interchangeably by the same names and notations since they are connected by one to one relation. For example, “the yield coefficient” and “the dimensionless yield coefficient” will be denoted by the same notation $\beta$ .

3. Results

3.1. Steady State Solutions and Their Stability

3.1.1. Steady State Solutions

The system of differential Equations (5) and (6) has one washout steady state solution $\left(S_0\mathrm{,}{X}_0\right)$ , where

$S_0=S_0,X_0=0$ (7)

and two no washout steady state solutions $\left(S_1\mathrm{,}{X}_1\right)$ and $\left(S_2\mathrm{,}{X}_2\right)$ where

$\begin{array}{l}{S}_1=\frac{-1-\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}+\tau }{2\gamma }\\ X_1=\frac{\gamma \left(1+2S_0\gamma +\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}-\tau \right)}{2{\gamma }^2}\\ +\frac{\beta \left(1+\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}-\tau \right)\left(-1+\tau \right)}{2{\gamma }^2}\\ -\frac{\beta \gamma \left(-2+S_0+S_0\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}-S_0\tau \right)}{2{\gamma }^2}\end{array}$ (8)

and

$\begin{array}{l}{S}_2=\frac{-1+\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}+\tau }{2\gamma }\\ X_2=\frac{\gamma \left(1+2S_0\gamma -\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}-\tau \right)}{2{\gamma }^2}\\ -\frac{\beta \left(-1+\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}+\tau \right)\left(-1+\tau \right)}{2{\gamma }^2}\\ +\frac{\beta \gamma \left(2-S_0+S_0\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}+S_0\tau \right)}{2{\gamma }^2}\end{array}$ (9)

If $\tau \ge 1+2\sqrt{\gamma }$ , then the substrate components, $S_1$ and $S_2$ , in (8) and (9) will be real and positive. Conditions for the cell mass concentration component, $X_1$ , in (8) to be real and non-negative are $\gamma >\frac{1}{S_0^2}$ and $\tau \ge 1+2\sqrt{\gamma }$ , and for the cell mass concentration, $X_2$ , in (9) to be real and non-negative are $\gamma >\frac{1}{S_0^2}$ and $1+2\sqrt{\gamma }\le \tau <\frac{1+S_0+S_0^2\gamma }{S_0}$ . These conditions for the cell mass concentration lead us to conclude that both no washout steady state solutions are physically meaningful whenever $\gamma >\frac{1}{S_0^2}$ and $1+2\sqrt{\gamma }\le \tau <\frac{1+S_0+S_0^2\gamma }{S_0}$ . We need to study the effect of variation of the residence time on the reactor’s steady state performance. That means, we need to investigate the relationship between the cell mass concentration X and the residence time $\tau$ in the no washout state. Specially, we want to determine the value of $\tau$ at which X will have optimum values.

From Equation (8), we can show that $\frac{\text{d}X}{\text{d}\tau }=0$ at

${\tau }_{\max }=\frac{-2\beta +4{\beta }^2+2S_0{\beta }^2+\gamma -2S_0\beta \gamma +S_0^2{\beta }^2\gamma }{2\beta \left(-1+S_0\beta \right)}$ (10)

subject to

$S_0>\frac{1}{\beta }$ (11)

and

$\frac{1}{S_0^2}<\gamma <\frac{4{\beta }^2}{1-2S_0\beta +S_0^2{\beta }^2}$ (12)

$X^{″}\left(\tau \right)$ will be negative if critical point ${\tau }_{\max }$ in (10) belongs to the following interval

$\begin{array}{l}1+2\sqrt{\gamma }<{\tau }_{\max }<\text{Root}[\beta -3{\beta }^2-S_0{\beta }^2-6\beta \gamma +16{\beta }^2\gamma +6S_0{\beta }^2\gamma +{\gamma }^2\\ -2S_0\beta {\gamma }^2+S_0^2{\beta }^2{\gamma }^2+\left(-3\beta +6{\beta }^2+3S_0{\beta }^2+6\beta \gamma -6S_0{\beta }^2\gamma \right)\#1\\ +\left(3\beta -3{\beta }^2-3S_0{\beta }^2\right)\#1^2+\left(-\beta +S_0{\beta }^2\right)\#1^3\&,3]\end{array}$

where the upper bound of ${\tau }_{\max }$ depends on $\beta \mathrm{,}\gamma$ , and $S_0$ . For example, if we take $\beta =2,\gamma =0.9,S_0=2$ , we get an upper bound of ${\tau }_{\max }$ at 3.313471. The ${\tau }_{\max }$ does not exceed the upper bound when

$\frac{1}{S_0^2}<\gamma <\frac{4}{S_0^2}$ (13)

or

$\gamma >\frac{4}{S_0^2}\text{and}\beta <2\sqrt{\frac{\gamma }{{\left(-4+S_0^2\gamma \right)}^2}}+\frac{S_0\gamma }{-4+S_0^2\gamma }$ (14)

Therefore, a local maximum value of X has occurred at ${\tau }_{\max }$ .

From Equation (9), we can show that $\frac{\text{d}X}{\text{d}\tau }=0$ at ${\tau }_{\max }$ , subject to the restriction that $S_0>\frac{1}{\beta }$ and $\frac{4{\beta }^2}{1-2S_0\beta +S_0^2{\beta }^2}$ . $X^{″}\left(\tau \right)$ will be negative if the critical point ${\tau }_{\max }$ is satisfied the following conditions.

$\gamma >\frac{4}{S_0^2}\text{and}\beta >2\sqrt{\frac{\gamma }{{\left(-4+S_0^2\gamma \right)}^2}}+\frac{S_0\gamma }{-4+S_0^2\gamma }$ (15)

Therefore, X has a local maximum value at ${\tau }_{\max }$ .

3.1.2. Stability of the Washout Steady State Solution

Jacobian Matrix: At the washout steady state solution (7), the system (5) and (6) has the following Jacobian matrix which has the eigenvalues $-\frac{1}{\tau }$ and $\frac{S_0}{1+S_0+S_0^2\gamma }-\frac{1}{\tau }$ .

$J=\left(\begin{array}{cc}-\frac{1}{\tau }& -\frac{S_0}{\left(1+S_0\beta \right)\left(1+S_0+S_0^2\gamma \right)}\\ 0& \frac{S_0}{1+S_0+S_0^2\gamma }-\frac{1}{\tau }\end{array}\right)$ (16)

Therefore, if $\tau <\frac{1}{S_0}+1+S_0\gamma$ then the washout steady state solution will be stable.

3.1.3. Stability of No Washout Steady State Solution (8)

Jacobian Matrix: At the washout steady state solution (8), the system (5) and (6) has the following Jacobian matrix

$J=\left(\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& 0\end{array}\right)$

where

$J_{11}=\frac{8{\gamma }^2\left(-4S_0{\gamma }^2+S_0\gamma \left(-1+\tau \right)\left(-1+\tau -\omega \right)-2\gamma \left(2-3\tau +\omega \right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}$

$\begin{array}{c}+\frac{8{\gamma }^2\left(\left(1-3\tau +2{\tau }^2\right)\left(1-\tau +\omega \right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{8\beta \gamma \left({\left(-1+\tau \right)}^2\left(-2+5\tau \right)\left(1-\tau +\omega \right)+2{\gamma }^2\left(-4+S_0\left(4-5\tau +2\omega \right)\right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{8\beta \gamma \left(\gamma \left(10+19{\tau }^2+S_0\left(2-5\tau +3{\tau }^2\right)\right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{8\beta {\gamma }^2\left(\left(-1+\tau -\omega \right)+6\omega -\tau \left(29+9\omega \right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\end{array}$

$\begin{array}{l}+\frac{8{\beta }^2\left(4S_0{\gamma }^3+{\left(-1+\tau \right)}^3\left(-1+3\tau \right)\left(1-\tau +\omega \right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{8{\beta }^2\left(-\gamma \left(-1+\tau \right)\left(-6+20\tau -14{\tau }^2-4\omega +8\tau \omega \right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{8{\beta }^2\left(-\gamma \left(-1+\tau \right)\right)\left(S_0\left(1-3\tau +2{\tau }^2\right)\left(1-\tau +\omega \right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{8{\beta }^2\left(+{\gamma }^2\left(2\left(4-6\tau +\omega \right)+S_0\left(-5-8{\tau }^2-3\omega +\tau \left(13+4\omega \right)\right)\right)\right)}{{\tau }^2{\left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}^2{\left(1-\tau +\omega \right)}^2}\end{array}$ (17)

$J_{12}=-\frac{2\gamma }{\tau \left(2\gamma +\beta \left(-1+\tau -\omega \right)\right)}$ (18)

$\begin{array}{c}{J}_{21}=\frac{\left(4\gamma +\left(-1+\tau \right)\left(1-\tau +\omega \right)\right)\left(\beta \gamma \left(2+S_0\left(-1+\tau -\omega \right)\right)\right)}{\gamma {\tau }^2{\left(1-\tau +\omega \right)}^2}\\ +\frac{\left(\gamma {\tau }^2{\left(1-\tau +\omega \right)}^2\right)\left(\beta \left(-1+\tau \right)\left(1-\tau +\omega \right)+\gamma \left(1+2S_0\gamma -\tau +\omega \right)\right)}{\gamma {\tau }^2{\left(1-\tau +\omega \right)}^2}\end{array}$ (19)

and $\omega =\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}$ .

By following the restrictions (11), (12), (13) and (14), it can be shown that $\text{Tr}\left(J\right)=J_{11}$ is negative and $\text{Det}\left(J\right)=-J_{21}{J}_{12}$ is positive at the critical point ${\tau }_{\max }$ . Hence, we can make a conclusion that the steady state solution (8) is stable at ${\tau }_{\max }$ and has practical importance to study.

If $\tau =1+2\sqrt{\gamma }$ or $\tau =\frac{1+S_0+S_0^2\gamma }{S_0}$ , then $Det\left(J\right)=J_{12}{J}_{21}=0$ and Jacobian matrix has a zero eigenvalue. The required conditions to have double zero eigenvalues of the Jacobian matrix are

$\text{Det}\left(J\right)=J_{12}{J}_{21}=0\text{and}\text{Tr}\left(J\right)=J_{11}+J_{22}=0.$

It can be shown that the conditions to have double zero eigenvalues are satisfied when

$S_0=\frac{1}{\beta }+\frac{2}{\sqrt{\gamma }}\text{and}1+2\sqrt{\gamma }=\tau$

Here $S_0$ is positive since $\beta >0$ and $\gamma >0$ . This positivity of $S_0$ in the feed immediately implies that double zero eigenvalues can occur.

3.1.4. Stability of No Washout Steady State Solution (9)

Jacobian Matrix: At the washout steady state solution (9), the system (5) and (6) has the following Jacobian matrix

$J=\left(\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& 0\end{array}\right)$

where

$\begin{array}{c}{J}_{11}=\frac{-8{\gamma }^2\left(4S_0{\gamma }^2+\left(1-3\tau +2{\tau }^2\right)\left(-1+\tau +\omega \right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ +\frac{-8{\gamma }^2\left(\gamma \left(S_0\left(-1+\tau \right)\left(-1+\tau +\omega \right)+2\left(-2+3\tau +\omega \right)\right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ -\frac{8\beta \gamma \left({\left(-1+\tau \right)}^2\left(-2+5\tau \right)\left(-1+\tau +\omega \right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ -\frac{8\beta \gamma \left(\gamma \left(10-29\tau +19{\tau }^2-6\omega +9\tau \omega +S_0\left(2-5\tau +3{\tau }^2\right)\left(-1+\tau +\omega \right)\right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ -\frac{8\beta \gamma \left(+2{\gamma }^2\left(4+S_0\left(-4+5\tau +2\omega \right)\right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\end{array}$

$\begin{array}{l}+\frac{8{\beta }^2\left(-4S_0{\gamma }^3+{\left(-1+\tau \right)}^3\left(-1+3\tau \right)\left(-1+\tau +\omega \right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ -\frac{8{\beta }^2\left(-\gamma \left(-1+\tau \right)\left(6+14{\tau }^2-4\omega +S_0\left(1-3\tau +2{\tau }^2\right)\left(-1+\tau +\omega \right)\right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ +\frac{8{\beta }^2\left(-\gamma \left(-1+\tau \right)4\tau \left(-5+2\omega \right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\\ -\frac{8{\beta }^2\left({\gamma }^2\left(2\left(-4+6\tau +\omega \right)+S_0\left(5+8{\tau }^2-3\omega +\tau \left(-13+4\omega \right)\right)\right)\right)}{{\tau }^2{\left(-1+\tau +\omega \right)}^2{\left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}^2}\end{array}$ (20)

$J_{12}=-\frac{2\gamma }{\tau \left(2\gamma +\beta \left(-1+\tau +\omega \right)\right)}$ (21)

$\begin{array}{c}{J}_{21}=\frac{\left(4\gamma -\left(-1+\tau \right)\left(-1+\tau +\omega \right)\right)\left(\gamma \left(1+2S_0\gamma -\tau -\omega \right)\right)}{\gamma {\tau }^2{\left(-1+\tau +\omega \right)}^2}\\ -\frac{\left(4\gamma -\left(-1+\tau \right)\left(-1+\tau +\omega \right)\right)\beta \left(-1+\tau \right)\left(-1+\tau +\omega \right)}{\gamma {\tau }^2{\left(-1+\tau +\omega \right)}^2}\\ +\frac{\left(4\gamma -\left(-1+\tau \right)\left(-1+\tau +\omega \right)\right)\left(\beta \gamma \left(2+S_0\left(-1+\tau +\omega \right)\right)\right)}{\gamma {\tau }^2{\left(-1+\tau +\omega \right)}^2}\end{array}$ (22)

and $\omega =\sqrt{-4\gamma +{\left(-1+\tau \right)}^2}$ .

By following the condition (15), the $\text{Tr}\left(J\right)$ is never negative and the $\text{Det}\left(J\right)$ is never positive at the critical point. Therefore, the maximum steady state solution (9) is not stable.

If $\tau =1+2\sqrt{\gamma }$ or $\tau =\frac{1+S_0+S_0^2\gamma }{S_0}$ , then $Det\left(J\right)=J_{12}{J}_{21}=0$ and Jacobian matrix has a zero eigenvalue. But the later value of $\tau$ is not relevant. The required conditions to have double zero eigenvalues of the Jacobian matrix are

$\text{Det}\left(J\right)=J_{12}{J}_{21}=0\text{and}\text{Tr}\left(J\right)=J_{11}+J_{22}=0.$

It can be shown that the conditions to have double zero eigenvalues are satisfied when

$S_0=\frac{1}{\beta }+\frac{2}{\sqrt{\gamma }}\text{and}1+2\sqrt{\gamma }=\tau .$

Here $S_0$ is positive since $\beta >0$ and $\gamma >0$ . This positivity of $S_0$ in the feed immediately implies that double zero eigenvalues can occur.

3.1.5. Hopf Bifurcation on the No Washout Steady State Solution (8)

A Hopf bifurcation will take place if $J_{11}=0$ and $J_{12}{J}_{21}<0$ . $J_{12}{J}_{21}$ will be negative when

$\tau >1+2\sqrt{\gamma }\text{and}1+2S_0\gamma +\gamma \sqrt{\frac{-4\gamma +{\left(-1+\tau \right)}^2}{{\gamma }^2}}-\tau >0.$ (23)

These conditions are equivalent to the following three cases. In each case, the values of the residence time correspond to the root of $J_{11}=0$ at which Hopf bifurcations occur .

Case I: $0<S_0\le 2$ , $\gamma >\frac{1}{S_0^2}$ , and $\tau \ge 1+2\sqrt{\gamma }$ . If we take $S_0=1.7,\gamma =2$ , and $\beta =5.25$ , then $J_{11}=0$ at $\tau =4.20144$ (Figure 1) (other values are possible).

A Degenerate Hopf bifurcation occurs at a value of residence time $\tau$ where two Hopf points resulting in a single point by annihilating each other which is known as H2₁ degeneracy and it occurs when $J_{11}=0$ and $\frac{\text{d}{J}_{11}}{\text{d}\tau }=0$ [10].

In this case, it is possible to show that when $J_{11}=0$ , $\frac{\text{d}{J}_{11}}{\text{d}\tau }$ vanishes for suitable values of $\gamma$ and $\beta$ . For example, if we take $\gamma =2$ , $\beta =5.25$ , then $\frac{\text{d}{J}_{11}}{\text{d}\tau }$ vanishes at $\left(S_0\mathrm{,}\tau \right)\approx \left(\mathrm{1.53265,3.90678}\right)$ . Therefore, if the substrate concentration $S_0$ is adequately small ( $S_0<1.53265$ ) or ( $S_0<2.68214\text{g}\cdot {\text{l}}^{-1}$ ) natural oscillations are not possible to occur for $\gamma =2$ , $\beta =5.25$ .

Case II: $S_0>2$ , $\frac{1}{S_0^2}<\gamma \le \frac{1}{4}$ , and $\tau \ge 1+2\sqrt{\gamma }$ . Let $S_0=3$ , $\gamma =\frac{8}{36}$ , $\beta =5.25$ , then $J_{11}=0$ at $\tau =2.281137$ and $\tau =2.467646$ (Figure 2) (other values are possible).

Figure 1. Hopf point unfolding curve. Parameter values: $\gamma =2,\beta =5.25$ .

Figure 2. Hopf point unfolding curve. Parameter values: $\gamma =\frac{8}{36},\beta =5.25$ .

For this situation, it is possible to show that when $J_{11}=0$ , $\frac{\text{d}{J}_{11}}{\text{d}\tau }$ vanishes for appropriate values of $\gamma$ and $\beta$ . For example, if $\gamma =\frac{8}{36}$ , $\beta =5.25$ , then $\frac{\text{d}{J}_{11}}{\text{d}\tau }$ vanishes at $\left(S_0\mathrm{,}\tau \right)\approx \left(\mathrm{2.97798,2.37231}\right)$ . Therefore, if the substrate concentration $S_0$ is sufficiently small ( $S_0<2.97798$ ) or ( $S_0<5.21147\text{g}\cdot {\text{l}}^{-1}$ ) natural oscillations are not possible for $\gamma =\frac{8}{36}$ , $\beta =5.25$ .

Case III: $S_0>2$ , $4\gamma >1$ , and $\tau \ge 1+2\sqrt{\gamma }$ . Let $S_0=3$ , $\gamma =\frac{1}{2}$ , $\beta =5.25$ , then $J_{11}=0$ at $\tau =3.181914$ and $\tau =2.414296$ (Figure 3) (other values are possible).

In this case, it is possible to show that when $J_{11}=0$ , $\frac{\text{d}{J}_{11}}{\text{d}\tau }$ vanishes for appropriate values of $\gamma$ and $\beta$ . For example, if we take $\gamma =\frac{1}{2}$ , $\beta =5.25$ , then

Figure 3. Hopf point unfolding curve. Parameter values: $\gamma =\frac{1}{2},\beta =5.25$ .

$\frac{\text{d}{J}_{11}}{\text{d}\tau }$ vanishes at $\left(S_0,\tau \right)\approx \left(2.43946,2.68646\right)$ . Therefore, if the substrate concentration $S_0$ is small enough ( $S_0<2.43946$ ) or ( $S_0<4.26906\text{g}\cdot {\text{l}}^{-1}$ ) natural oscillations are not possible for $\gamma =\frac{1}{1}$ , $\beta =5.25$ .

In all cases above, we have used parameter values $K_s=1.75{\text{g}}^{-1}$ , $\alpha =0.01$ , $\beta =0.031{\text{g}}^{-1}$ , ${\mu }_m=0.3{\text{h}}^{-1}$ , and ${\beta }^{*}=5.25$ from [6] [9].

The work reported in [3] [6] [9] used $S_0\ge 10\text{g}\cdot {\text{l}}^{-1}$ without confirming conditions under which natural oscillations are possible. The work in [8] reported that natural oscillations are impossible for sufficiently small substrate concentration $S_0<6.84203\text{g}\cdot {\text{l}}^{-1}$ , for $\beta =5.25$ . Both cases above were observed for the specific growth rate equation $\mu \left(S\right)=\frac{{\mu }_{m}S}{K_s+S}$ . In our work, we have used $\mu \left(S\right)=\frac{{\mu }_{m}S}{K_s+S+\frac{S^2}{K_i}}$ . It is possible to gain natural oscillations for even significantly smaller values of $S_0$ .

4. Conclusion

We have studied, analytically, a simple chemostat model in a flow reactor with a variable yield coefficient in which the growth rate is taken to be a Haldane expression. In this study, three steady-state solutions have been discussed, which characterize no washout and washout circumstances in the closed photobioreactor. Under suitable parameter values, a stable steady-state solution attains its maximum value. We also corroborated the parameter ranges for the model with variable yield coefficient and growth rate, which describes natural oscillations in the chemostat. In all three cases of Hopf bifurcation analysis, we have found that natural oscillations can be achieved at significantly lower values of the substrate concentrations. In cases I, II, and III, the value is $S_0>2.68214\text{g}\cdot {\text{l}}^{-1}$ , $S_0>5.21147\text{g}\cdot {\text{l}}^{-1}$ , and $S_0>4.26906\text{g}\cdot {\text{l}}^{-1}$ respectively. A previous study reported that natural oscillation is achieved at a substrate concentration $S_0>6.84203\text{g}\cdot {\text{l}}^{-1}$ and results claim that natural oscillations are possible for sufficiently large values of substrate concentration ( $S_0>10\text{g}\cdot {\text{l}}^{-1}$ ) but the conditions of natural oscillations were not reported.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Silveston, P., Hudgins, R.R. and Renken, A. (1995) Periodic Operation of Catalytic Reactors—Introduction and Overview. Catalysis Today, 25, 91-112. https://doi.org/10.1016/0920-5861(95)00101-K
[2]	Stankiewicz, A. and Kuczynski, M. (1995) An Industrial View on the Dynamic Operation of Chemical Converters. Chemical Engineering and Processing: Process Intensification, 34, 367-377. https://doi.org/10.1016/0255-2701(95)00537-4
[3]	Balakrishnan, A. and Yang, R.Y.K. (1998) Improvement of Chemostat Performance via Nonlinear Oscillations, Part 2. Extension to Other Systems. ACH Models in Chemistry, 135, 1-18.
[4]	Shen, J.Q. and Ray, A.K. (2000) Performance Improvement of Activated Sludge Wastewater Treatment by Nonlinear Natural Oscillations. Chemical Engineering & Technology: Industrial Chemistry-Plant Equipment-Process Engineering-Biotechnology, 23, 1115-1122. https://doi.org/10.1002/1521-4125(200012)23:12<1115::AID-CEAT1115>3.0.CO;2-T
[5]	Ray, A.K. (1995) Performance Improvement of a Chemical Reactor by Non-Linear Natural Oscillations. The Chemical Engineering Journal and the Biochemical Engineering Journal, 59, 169-175. https://doi.org/10.1016/0923-0467(94)02920-2
[6]	Yang, R.Y.K. and Su, J. (1993) Improvement of Chemostat Performance via Nonlinear Oscillations. Bioprocess Engineering, 9, 97-102. https://doi.org/10.1007/BF00369037
[7]	Dae-Hyun, C., Rishiram, R., et al. (2016) Influence of Limiting Factors on Biomass and Lipid Productivities of Axenic Chlorella Vulgaris in Photobioreactor under Chemostat Cultivation. Bioresource Technology, 211, 367-373. https://doi.org/10.1016/j.biortech.2016.03.109
[8]	Nelson, M.I. and Sidhu, H.S. (2005) Analysis of a Chemostat Model with Variable Yield Coefficient. Journal of Mathematical Chemistry, 38, 605-615. https://doi.org/10.1007/s10910-005-6914-2
[9]	Balakrishnan, A. and Yang, R.Y.K. (2002) Self-Forcing of a Chemostat with Self-Sustained Oscillations for Productivity Enhancement. Chemical Engineering Communications, 189, 1569-1585. https://doi.org/10.1080/00986440214998
[10]	Gray, B.F. and Roberts, M.J. (1988) A Method for the Complete Qualitative Analysis of Two Coupled Ordinary Differential Equations Dependent on Three Parameters. Proceedings of the Royal Society A. Mathematical and Physical Sciences, 416, 361-389. https://doi.org/10.1098/rspa.1988.0039