Hybrid Delayed Feedback Controller Design in a Nutrient-Microorganism System Accompanying Time Delay

1. Introduction

In the sediment of ocean, a great deal of chemical substances and microorganism will appear and they are closely related to each other, which formulate various intricate networks. Thus exploring the law of interaction among them is a very complex problem. In order to reveal the inherent law of chemical substances and microorganism, a lot of works on chemical substances and microorganism have been carried out. For example, Gao et al. [1] explored the spatiotemporal pattern formation and global existence of a nutrient-microorganism system involving nutrient-taxis in the sediment. Baurmann and Feudel [2] investigated the turing patterns of a nutrient-microorganism model. Cao and Wu [3] dealt with the global bifurcation and local bifurcation of both positive constant solutions of a diffusive nutrient-microorganism model. For more concrete contents, one can see [4-8].

In 2004, Baurmann and Feudel [9] proposed the following nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d\mathcal{B}\left(t\right)}{dt}=\alpha {\mathcal{B}}^a\left(t\right)\mathcal{N}\left(t\right)-m\mathcal{B}\left(t\right),}\hfill \\ {\displaystyle \frac{d\mathcal{N}\left(t\right)}{dt}=\varphi (\overline{\mathcal{N}}-\mathcal{N}\left(t\right))-\beta {\mathcal{B}}^a\left(t\right)\mathcal{N}\left(t\right),}\hfill \\ {\displaystyle {\mathcal{B}}^a\left(t\right)=\frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(t\right)+\mathcal{L}}\mathcal{B}\left(t\right),}\hfill \end{array}\right.$$

(1)

where $\mathcal{B}\left(t\right)$ stands for the density of microorganism and $\mathcal{N}\left(t\right)$ stands for the density of one chemical species in the sediment. The bacteria includes two parts (active part and inactive part). ${\mathcal{B}}^a$ stands for the active part bacterium, it grows by feeding on the nutrient with the rate $\alpha$. Noticing that ${\mathcal{B}}^a\left(t\right)\le \mathcal{B}$, we let ${\mathcal{B}}^a\left(t\right)={\displaystyle \frac{\mathcal{B}\left(t\right)}{\mathcal{B}\left(t\right)+L}}\mathcal{B}\left(t\right)$, where $\mathcal{L}$ denotes a constant. m denotes the bacterium mortality. The input of nutrient into the system principally derives from the seawater via the surface of sediment by applying bioirrigation which is describe by $\varphi (\overline{\mathcal{N}}-\mathcal{N}\left(t\right))$, where $\varphi$ stands for the hydraulic conductivity between seawater and sediment, and $\overline{\mathcal{N}}$ represents the concentration of nutrient in the seawater. $\beta$ denotes the capture rate. For more detailed implication of system (1.1), one can see [1]. Applying the following nondimensionalizing variable substitutions:

$$\mathcal{N}={\displaystyle \frac{m}{\alpha }}{u}_2,\overline{\mathcal{N}}={\displaystyle \frac{m}{\alpha }}\beta ,\mathcal{B}={\displaystyle \frac{\varphi }{\beta }}{u}_1,\mathcal{L}={\displaystyle \frac{\varphi }{\beta }}\mathcal{K},\tau ={\displaystyle \frac{1}{\varphi }}t,$$

(2)

we get

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}=\frac{c{u}_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}-c{u}_1\left(t\right),}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}=\beta -u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}},}\hfill \end{array}\right.$$

(3)

where $c={\displaystyle \frac{m}{\varphi }}$. In many cases, the density of microorganism and the density of one chemical species in the sediment depend on the development of the current time and the past history. Here we assume that the density of microorganism depends on the current time and the past history, i.e, there is feedback delay for the variation of the density of microorganism. Relying on this viewpoint, we formulate the following delayed nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}=\frac{c{u}_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+\mathcal{K}}-c{u}_1\left(t\right),}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}=\beta -u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}},}\hfill \end{array}\right.$$

(4)

where $\rho$ is a delay. All the parameters in system (1.4) are positive real numbers.

Generally speaking, delay plays a vital role in affecting the dynamical phenomenon of differential dynamical systems. In numerous situations, delay usually leads to change of stability, onset of bifurcation phenomenon and the emergence of chaos and so on [10-21]. Especially, delay-driven bifurcation phenomenon is a vital dynamical behavior. Chemically, delay-driven bifurcation phenomenon is able to efficaciously describe the balanced relations among the concentrations of various chemical substances. As a result, it is of importance to focus on the delay-driven bifurcation phenomenon for various chemical reaction systems and biological models. Stimulated by this viewpoint above, we will deal with the delay-driven bifurcation phenomenon and bifurcation control aspect for model (1.4). In particular, we will focus on the following key issues: (a) Explore the well-posedness (e.g., boundedness, existence and uniqueness, non-negativeness) of solution of model (1.4). (b) Investigate the onset of bifurcation phenomenon and stability behavior of model (1.4). (c) Design two differential controllers to adjust the stability domain and the time of appearance of bifurcation phenomenon of system (1.4).

The main luminescent spot of this work are summarized as follows: (I) Relying on the earlier studies, a novel stability and bifurcation condition without depending on delay for model (1.4) is formulated. (II) Taking advantage of both different hybrid controllers, stability range and the time of onset of bifurcation phenomenon of model (1.4) are successfully adjusted. (III) The effect of delay on controlling bifurcation behavior and stabilizing the density of microorganism and the density of one chemical species in the sediment of model (1.4) is provided.

This organization of this study is provided as follows. The well-posedness including boundedness, existence and uniqueness and non-negativeness of the solution to model (1.4) is analyzed in “Well-posedness" section. “Bifurcation phenomenon" section studies the Hopf bifurcation issue and stability trait of system (1.4). “Bifurcation control via hybrid delayed feedback controller" section deals with the control issue of Hop bifurcation for system (1.4) by designing a proper hybrid delayed feedback controller including parameter perturbation with delay and state feedback. “Bifurcation control via extended hybrid delayed feedback controller" section focuses on the control issue of Hopf bifurcation for system (1.4) by designing a proper extended hybrid delayed feedback controller including parameter perturbation with delay and state feedback. “Numerical experiments" implements software experiments to check the effectiveness of the gained primary outcomes. A concise conclusion is given to end this article in “Conclusions".

2. Well-Posedness

In this section, we are going to explore the well-posedness of solution of model (1.4)(involve non-negativeness, existence and uniqueness and boundedness) by utilizing fixed point theory, inequality tactics and establishment of an appropriate function.

Theorem 1. 

Denote ${\rm Y}=\{u_1,u_2)\in R^2:max\left\{\right|u_1|,|u_2\left|\right\}\le U\}$, where U denotes a positive constant. For every $(u_{10},u_{20})\in {\rm Y}$, model (1.4) under the initial value $(u_{10},u_{20})$ possesses a unique solution $\mathcal{U}=(u_1,u_2)\in {\rm Y}.$

Proof. 

Let

$$\Gamma \left(\mathcal{U}\right)=({\Gamma }_1\left(\mathcal{U}\right),{\Gamma }_2\left(\mathcal{U}\right)),$$

(5)

where

$$\left\{\begin{array}{c}{\displaystyle {\Gamma }_1\left(\mathcal{U}\right)=\frac{c{u}_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+\mathcal{K}}-c{u}_1\left(t\right),}\hfill \\ {\displaystyle {\Gamma }_2\left(\mathcal{U}\right)=\beta -u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}.}\hfill \end{array}\right.$$

(6)

For each $\mathcal{U},\overline{\mathcal{U}}\in \Phi$, one acquires

$$\begin{array}{ccc}& & \left|\right|\Gamma \left(\mathcal{U}\right)-\Gamma \left(\overline{\mathcal{U}}\right)\left|\right|\hfill \\ & & =\left|{\displaystyle \frac{c{u}_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+\mathcal{K}}}-c{u}_1\left(t\right)-\left[{\displaystyle \frac{c{\overline{u}}_1^2(t-\rho ){\overline{u}}_2\left(t\right)}{{\overline{u}}_1(t-\rho )+\mathcal{K}}}-c{\overline{u}}_1\left(t\right)\right]\right|\hfill \\ & & +\left|\beta -u_2\left(t\right)-{\displaystyle \frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}}-\left[\beta -{\overline{u}}_2\left(t\right)-{\displaystyle \frac{{\overline{u}}_1^2\left(t\right){\overline{u}}_2\left(t\right)}{{\overline{u}}_1\left(t\right)+\mathcal{K}}}\right]\right|\hfill \\ & & \le {\displaystyle \frac{c{\mathcal{U}}^3}{{\mathcal{K}}^2}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|+{\displaystyle \frac{3c{\mathcal{U}}^3}{{\mathcal{K}}^2}}|u_1(t-\rho )-{\overline{u}}_1(t-\rho )|\hfill \\ & & +{\displaystyle \frac{c{\mathcal{U}}^2}{\mathcal{K}}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|+{\displaystyle \frac{c{\mathcal{U}}^3}{\mathcal{K}}}|u_1(t-\rho )-{\overline{u}}_1(t-\rho )|\hfill \\ & & +c|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|\hfill \\ & & +{\displaystyle \frac{{\mathcal{U}}^4}{{\mathcal{K}}^2}}|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+{\displaystyle \frac{{\mathcal{U}}^4}{{\mathcal{K}}^2}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|\hfill \\ & & +{\displaystyle \frac{2{\mathcal{U}}^2}{\mathcal{K}}}|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+{\displaystyle \frac{{\mathcal{U}}^2}{\mathcal{K}}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|\hfill \\ & & \le {\displaystyle \frac{c{\mathcal{U}}^3}{{\mathcal{K}}^2}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|+{\displaystyle \frac{3c{\mathcal{U}}^3}{{\mathcal{K}}^2}}|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|\hfill \\ & & +{\displaystyle \frac{c{\mathcal{U}}^2}{\mathcal{K}}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|+{\displaystyle \frac{c{\mathcal{U}}^3}{\mathcal{K}}}|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|\hfill \\ & & +c|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|\hfill \\ & & +{\displaystyle \frac{{\mathcal{U}}^4}{{\mathcal{K}}^2}}|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+{\displaystyle \frac{{\mathcal{U}}^4}{{\mathcal{K}}^2}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|\hfill \\ & & +{\displaystyle \frac{2{\mathcal{U}}^2}{\mathcal{K}}}|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+{\displaystyle \frac{{\mathcal{U}}^2}{\mathcal{K}}}|u_2\left(t\right)-{\overline{u}}_2\left(t\right)|\hfill \\ & & ={\vartheta }_1|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|+{\vartheta }_2|u_1\left(t\right)-{\overline{u}}_1\left(t\right)|.\hfill \end{array}$$

(7)

where

$$\left\{\begin{array}{c}{\displaystyle {\vartheta }_1=\frac{3c{\mathcal{U}}^3}{{\mathcal{K}}^2}+\frac{c{\mathcal{U}}^3}{\mathcal{K}}+c+\frac{{\mathcal{U}}^4}{{\mathcal{K}}^2}+\frac{2{\mathcal{U}}^2}{\mathcal{K}},}\hfill \\ {\displaystyle {\vartheta }_2=\frac{c{\mathcal{U}}^3}{{\mathcal{K}}^2}+\frac{c{\mathcal{U}}^2}{\mathcal{K}}+1+\frac{{\mathcal{U}}^4}{{\mathcal{K}}^2}+\frac{{\mathcal{U}}^2}{\mathcal{K}}.}\hfill \end{array}\right.$$

(8)

Set

$$\vartheta =max\{{\vartheta }_1,{\vartheta }_2\}.$$

(9)

Then by (2.3), one gets

$$\left|\right|\Gamma \left(\mathcal{U}\right)-\Gamma \left(\overline{\mathcal{U}}\right)\left|\right|\le \vartheta \left|\right|\mathcal{U}-\overline{\mathcal{U}}\left|\right|.$$

(10)

Then $\Gamma \left(\mathcal{U}\right)$ conforms to Lipschitz condition concerning $\mathcal{U}$ (see [22]). By exploiting fixed point theorem, we can lightly acquire that Theorem 2.1 holds. □

Theorem 2. 

Suppose that $\rho =0$, then every solution of system (1.4) beginning with $R_{+}^2$ is non-negative.

Proof. 

Let $U\left(t_0\right)=(u_1\left(t_0\right),u_2\left(t_0\right))$ be the initial condition of system (1.4). If there exists a real constant $t_{0★}>0$ conforming to $t_0<t<t_{0★}$ such that

$$\left\{\begin{array}{c}{\displaystyle u_1\left(t\right)>0,t_0<t<t_{0★},}\hfill \\ {\displaystyle u_1\left(t_{0★}\right)=0,}\hfill \\ {\displaystyle u_1\left(t_{0★}^{+}\right)<0.}\hfill \end{array}\right.$$

(11)

By applying model (1.4), one acquires

$${\displaystyle \frac{d{u}_1\left(t\right)}{dt}}{|}_{u_1\left(t_{0★}\right)=0}=0.$$

(12)

By using Lemma 1 of Das et al. [23], we know that $u_1\left(t_{0★}^{+}\right)=0$, which implies a contradiction (see (2.7)). Then $u_1\left(t\right)\ge 0$ for $t\ge t_0.$ Using the same scheme, one can lightly know that $u_2\left(t\right)\ge 0$ for $t\ge t_0.$ □

Theorem 3. 

If $\rho =0,c<1$, then all solutions to model (1.4) starting with $R_{+}^2$ are uniformly bounded.

Proof. 

Set

$$V\left(t\right)=u_1\left(t\right)+u_2\left(t\right).$$

(13)

Then

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dV\left(t\right)}{dt}}& =& {\displaystyle \frac{d{u}_1\left(t\right)}{dt}}+{\displaystyle \frac{d{u}_2\left(t\right)}{dt}}\hfill \\ & =& {\displaystyle \frac{c{u}_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}}-c{u}_1\left(t\right)+\beta -u_2\left(t\right)-{\displaystyle \frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}}\hfill \\ & =& -c{u}_1\left(t\right)-u_2\left(t\right)+\beta -(1-c){\displaystyle \frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}}\hfill \\ & \le & -c{u}_1\left(t\right)-u_2\left(t\right)+\beta \hfill \\ & \le & -cV\left(t\right)+\beta ,\hfill \end{array}$$

(14)

Making use of Gronwall’s inequality [24], one gains

$$V\left(t\right)\to {\displaystyle \frac{\beta }{c}},ast\to \infty .$$

(15)

The proof of Theorem 2.3 ends.

□

3. Bifurcation Phenomenon

Assume that $E(u_{1★},u_{2★})$ is the equilibrium point of model (1.4), then $u_{1★},u_{2★}$ obey

$$\left\{\begin{array}{c}{\displaystyle \frac{c{u}_{1★}^2{u}_{2★}}{u_{1★}+\mathcal{K}}-c{u}_{1★}=0,}\hfill \\ {\displaystyle \beta -u_{2★}-\frac{u_{1★}^2{u}_{2★}}{u_{1★}+\mathcal{K}}=0.}\hfill \end{array}\right.$$

(16)

The linear equation of model (1.4) at $E(u_{1★},u_{2★})$ is given by

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}=a_1{u}_1\left(t\right)+a_2{u}_2\left(t\right)+a_3{u}_1(t-\rho ),}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}=a_4{u}_1\left(t\right)+a_5{u}_2\left(t\right),}\hfill \end{array}\right.$$

(17)

where

$$\left\{\begin{array}{c}{\displaystyle a_1=-c,}\hfill \\ {\displaystyle a_2=\frac{c{u}_{1★}^2}{u_{1★}+\mathcal{K}},}\hfill \\ {\displaystyle a_3=\frac{2c{u}_{1★}{u}_{2★}}{u_{1★}+\mathcal{K}}-\frac{2c{u}_{1★}^2{u}_{2★}}{{(u_{1★}+\mathcal{K})}^2},}\hfill \\ {\displaystyle a_4=\frac{2u_{1★}{u}_{2★}}{\mathcal{K}}-\frac{u_{1★}^2}{{\mathcal{K}}^2},}\hfill \\ {\displaystyle a_5=\frac{u_{1★}^2}{\mathcal{K}}.}\hfill \end{array}\right.$$

(18)

The characteristic equation of system (3.2) reads as

$$det\left[\begin{array}{cc}\lambda -a_1-a_3{e}^{-\lambda \rho }& -a_2\\ -a_4& \lambda -a_5\end{array}\right]=0,$$

(19)

which results in

$${\lambda }^2+b_1\lambda +b_2+(b_3\lambda +b_4)e^{-\lambda \rho }=0,$$

(20)

where

$$\left\{\begin{array}{c}{\displaystyle b_1=-(a_1+a_5),}\hfill \\ {\displaystyle b_2=a_3{a}_5-a_2{a}_4,}\hfill \\ {\displaystyle b_3=-a_3,}\hfill \\ {\displaystyle b_4=a_3{a}_5.}\hfill \end{array}\right.$$

(21)

When $\rho =0,$ then Equation (3.5) takes the following expression:

$${\lambda }^2+(b_1+b_3)\lambda +b_2+b_4=0.$$

(22)

If

$$\left(M_1\right)b_1+b_3>0,b_2+b_4>0$$

is fulfilled, then the two roots ${\lambda }_1$ and ${\lambda }_2$ of Eq. (3.7) admit negative real parts. So one can know that the equilibrium point $E(u_{1★},u_{2★})$ of model (1.4) under the time delay $\rho =0$ preserves locally asymptotically stable form.

If $\lambda =i\xi$ is the root of Eq. (3.5), then according to Equation (3.5), one gets

$${\left(i\xi \right)}^2+b_1\left(i\xi \right)+b_2+(b_{3}i\xi +b_4)e^{-i\xi \rho }=0,$$

(23)

which generates

$$\left\{\begin{array}{c}{\displaystyle b_{4}cos\xi \rho +b_3\xi sin\xi \rho ={\xi }^2-b_2,}\hfill \\ {\displaystyle b_3\xi cos\xi \rho -b_{4}sin\xi \rho =-b_1\xi .}\hfill \end{array}\right.$$

(24)

In view of (3.9), we gain

$$b_4^2+{\left(b_3\xi \right)}^2={({\xi }^2-b_2)}^2+{\left(b_1\xi \right)}^2,$$

(25)

which results in

$${\xi }^4+(b_1^2-b_3^2-2b_2){\xi }^2-b_4^2=0$$

(26)

Denote

$${\Psi }_1\left(\xi \right)={\xi }^4+(b_1^2-b_3^2-2b_2){\xi }^2-b_4^2.$$

(27)

Notice that ${\Psi }_1\left(0\right)=-b_4^2<0$ and ${lim}_{\xi \to +\infty }{\Psi }_1\left(\xi \right)=+\infty >0$, then we acquire that Eq. (3.11) owns at least one real positive root. Then Equation (3.5) admits at least one pair of purely roots. Without loss of generality, assume that Eq. (3.11) has four real positive roots (say ${\xi }_i,i=1,2,3,4$. By (3.9), we gain

$${\rho }_i^{\left(k\right)}={\displaystyle \frac{1}{{\xi }_i}}\left[arccos\left({\displaystyle \frac{b_4({\xi }_i^2-b_2)-b_1{b}_3{\xi }_i^2}{b_4^2+b_3^2{\xi }_i^2}}\right)+2k\pi \right],$$

(28)

where $i=1,2,3,4;k=0,1,2,\cdots .$ Set ${\rho }_0={min}_{\{i=1,2,3,4;k=0,1,2,\cdots \}}\left\{{\xi }_i^{\left(k\right)}\right\}$ and then if $\rho ={\rho }_0$, (3.5) has a pair of imaginary roots $\pm i{\xi }_0$.

Now we give the assumption as follows:

$$\left(M_2\right){\mathcal{A}}_{1R}{\mathcal{A}}_{2R}+{\mathcal{A}}_{1I}{\mathcal{A}}_{2I}>0,$$

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{A}}_{1R}=b_1+b_{3}cos{\xi }_0{\rho }_0,}\hfill \\ {\displaystyle {\mathcal{A}}_{1I}=2{\xi }_0-b_{3}sin{\xi }_0{\rho }_0,}\hfill \\ {\displaystyle {\mathcal{A}}_{2R}=b_4{\xi }_{0}sin{\xi }_0{\rho }_0-b_3{\xi }_{0}cos{\xi }_0{\rho }_0,}\hfill \\ {\displaystyle {\mathcal{A}}_{2I}=b_4{\xi }_{0}cos{\xi }_0{\rho }_0+b_3{\xi }_{0}sin{\xi }_0{\rho }_0.}\hfill \end{array}\right.$$

(29)

Lemma 1. 

If $\lambda \left(\rho \right)={\iota }_1\left(\rho \right)+i{\iota }_2\left(\rho \right)$ is a root of Eq. (3.5) at $\rho ={\rho }_0$ obeying ${\iota }_1\left({\rho }_0\right)=0,{\iota }_2\left({\rho }_0\right)={\xi }_0$, then $\mathrm{Re}\left({\displaystyle \frac{d\lambda }{d\rho }}\right){|}_{\rho ={\rho }_0,\xi ={\xi }_0}>0.$

Proof. 

Based on Equation (3.5), one gains

$$(2\lambda +b_1){\displaystyle \frac{d\lambda }{d\rho }}+b_3{\displaystyle \frac{d\lambda }{d\rho }}{e}^{-\lambda \rho }-e^{-\lambda \rho }\left({\displaystyle \frac{d\lambda }{d\rho }}\rho +\lambda \right)\left(b_3\lambda +b_4\right)=0,$$

(30)

which generates

$${\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}={\displaystyle \frac{{\mathcal{A}}_1\left(\lambda \right)}{{\mathcal{A}}_2\left(\lambda \right)}}-{\displaystyle \frac{\rho }{\lambda }},$$

(31)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{A}}_1\left(\lambda \right)=2\lambda +b_1+b_3{e}^{-\lambda \rho },}\hfill \\ {\displaystyle {\mathcal{A}}_2\left(\lambda \right)=\lambda (b_3\lambda +b_4)e^{-\lambda \vartheta }.}\hfill \end{array}\right.$$

(32)

Then

$$\mathrm{Re}{\left[{\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}\right]}_{\rho ={\rho }_0,\xi ={\xi }_0}=\mathrm{Re}{\left[{\displaystyle \frac{{\mathcal{A}}_1\left(\lambda \right)}{{\mathcal{A}}_2\left(\lambda \right)}}\right]}_{\rho ={\rho }_0,\xi ={\xi }_0}={\displaystyle \frac{{\mathcal{A}}_{1R}{\mathcal{A}}_{2R}+{\mathcal{A}}_{1I}{\mathcal{A}}_{2I}}{{\mathcal{A}}_{2R}^2+{\mathcal{A}}_{2I}^2}}.$$

(33)

In view of $\left(M_2\right)$, one gains

$$\mathrm{Re}{\left[{\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}\right]}_{\rho ={\rho }_0,\xi ={\xi }_0}>0.$$

(34)

The proof finishes. □

Relying on the analysis above, the following outcome is lightly acquired.

Theorem 4. 

If $\left(M_1\right)$ and $M_2)$ are fulfilled, then the equilibrium point $E(u_{1★},u_{2★})$ of model (1.4) keeps locally asymptotically stable condition if $\rho \in [0,{\rho }_0)$ and model (1.4) produces a cluster of periodic solutions (namely, Hopf bifurcations) near the equilibrium point $E(u_{1★},u_{2★})$ when $\rho ={\rho }_0.$

4. Bifurcation Control via Hybrid Delayed Feedback Controller

In this section, we are to use a hybrid controller (include a state feedback and parameter perturbation) to the first equation of model (1.4). By virtue of this controller, we can dominate the time of emergence of bifurcation phenomenon and stability region of system (1.4). Relying on the ideas of [25-29], we formulate the following fractional controlled nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}=\sigma \left[\frac{c{u}_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+\mathcal{K}}-c{u}_1\left(t\right)\right]+(1-\sigma )[u_1(t-\rho )-u_1\left(t\right)],}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}=\beta -u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}},}\hfill \end{array}\right.$$

(35)

where $\sigma$ represents feedback gain coefficient. Model (4.1) owns the same equilibrium point $E(u_{1★},u_{2★})$ as that of model (1.4). The linear equation of model (4.1) at $E(u_{1★},u_{2★})$ is given by

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}={\tilde{a}}_1{u}_1\left(t\right)+{\tilde{a}}_2{u}_2\left(t\right)+{\tilde{a}}_3{u}_1(t-\rho ),}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}={\tilde{a}}_4{u}_1\left(t\right)+{\tilde{a}}_5{u}_2\left(t\right),}\hfill \end{array}\right.$$

(36)

where

$$\left\{\begin{array}{c}{\displaystyle {\tilde{a}}_1=-\sigma c-(1-\sigma ),}\hfill \\ {\displaystyle {\tilde{a}}_2=\frac{\sigma c{u}_{1★}^2}{u_{1★}+\mathcal{K}},}\hfill \\ {\displaystyle {\tilde{a}}_3=\frac{2\sigma c{u}_{1★}{u}_{2★}}{u_{1★}+\mathcal{K}}-\frac{2\sigma c{u}_{1★}^2{u}_{2★}}{{(u_{1★}+\mathcal{K})}^2}+1-\sigma ,}\hfill \\ {\displaystyle {\tilde{a}}_4=\frac{2u_{1★}{u}_{2★}}{\mathcal{K}}-\frac{u_{1★}^2}{{\mathcal{K}}^2},}\hfill \\ {\displaystyle {\tilde{a}}_5=\frac{u_{1★}^2}{\mathcal{K}}.}\hfill \end{array}\right.$$

(37)

The characteristic equation of system (4.2) reads as

$$det\left[\begin{array}{cc}\lambda -{\tilde{a}}_1-{\tilde{a}}_3{e}^{-\lambda \rho }& -{\tilde{a}}_2\\ -{\tilde{a}}_4& \lambda -{\tilde{a}}_5\end{array}\right]=0,$$

(38)

which results in

$${\lambda }^2+{\tilde{b}}_1\lambda +{\tilde{b}}_2+({\tilde{b}}_3\lambda +{\tilde{b}}_4)e^{-\lambda \rho }=0,$$

(39)

where

$$\left\{\begin{array}{c}{\displaystyle {\tilde{b}}_1=-({\tilde{a}}_1+{\tilde{a}}_5),}\hfill \\ {\displaystyle {\tilde{b}}_2={\tilde{a}}_3{\tilde{a}}_5-{\tilde{a}}_2{\tilde{a}}_4,}\hfill \\ {\displaystyle {\tilde{b}}_3=-{\tilde{a}}_3,}\hfill \\ {\displaystyle {\tilde{b}}_4={\tilde{a}}_3{\tilde{a}}_5.}\hfill \end{array}\right.$$

(40)

When $\rho =0,$ then Equation (4.5) takes the following expression:

$${\lambda }^2+({\tilde{b}}_1+{\tilde{b}}_3)\lambda +{\tilde{b}}_2{+}_4=0.$$

(41)

If

$$\left(M_3\right){\tilde{b}}_1+{\tilde{b}}_3>0,{\tilde{b}}_2+{\tilde{b}}_4>0$$

is fulfilled, then the two roots ${\lambda }_1$ and ${\lambda }_2$ of Eq. (4.7) admit negative real parts. So one can understand that the equilibrium point $E(u_{1★},u_{2★})$ of model (4.1) under the time delay $\rho =0$ preserves locally asymptotically stable form.

If $\lambda =i\zeta$ is the root of Eq. (4.5), then according to Equation (4.5), one gets

$${\left(i\zeta \right)}^2+{\tilde{b}}_1\left(i\zeta \right)+{\tilde{b}}_2+({\tilde{b}}_{3}i\zeta +{\tilde{b}}_4)e^{-i\zeta \rho }=0,$$

(42)

which generates

$$\left\{\begin{array}{c}{\displaystyle {\tilde{b}}_{4}cos\zeta \rho +{\tilde{b}}_3\zeta sin\zeta \rho ={\zeta }^2-{\tilde{b}}_2,}\hfill \\ {\displaystyle {\tilde{b}}_3\zeta cos\zeta \rho -{\tilde{b}}_{4}sin\zeta \rho =-{\tilde{b}}_1\zeta .}\hfill \end{array}\right.$$

(43)

In view of (4.9), we gain

$${\tilde{b}}_4^2+{\left({\tilde{b}}_3\zeta \right)}^2={({\zeta }^2-{\tilde{b}}_2)}^2+{\left({\tilde{b}}_1\zeta \right)}^2,$$

(44)

which results in

$${\zeta }^4+({\tilde{b}}_1^2-{\tilde{b}}_3^2-2{\tilde{b}}_2){\zeta }^2-{\tilde{b}}_4^2=0$$

(45)

Denote

$${\Psi }_2\left(\zeta \right)={\zeta }^4+({\tilde{b}}_1^2-{\tilde{b}}_3^2-2{\tilde{b}}_2){\zeta }^2-{\tilde{b}}_4^2.$$

(46)

Notice that ${\Psi }_2\left(0\right)=-{\tilde{b}}_4^2<0$ and ${lim}_{\zeta \to +\infty }{\Psi }_2\left(\zeta \right)=+\infty >0$, then we acquire that Eq. (4.11) owns at least one real positive root. Then Equation (4.5) admits at least one pair of purely roots. Without loss of generality, assume that Eq. (4.11) has four real positive roots (say ${\zeta }_j,j=1,2,3,4$. By (4.9), we gain

$${\rho }_j^{\left(l\right)}={\displaystyle \frac{1}{{\zeta }_j}}\left[arccos\left({\displaystyle \frac{{\tilde{b}}_4({\zeta }_j^2-{\tilde{b}}_2)-{\tilde{b}}_1{\tilde{b}}_3{\zeta }_j^2}{{\tilde{b}}_4^2+{\tilde{b}}_3^2{\zeta }_j^2}}\right)+2l\pi \right],$$

(47)

where $j=1,2,3,4;l=0,1,2,\cdots .$ Set ${\rho }_{0★}={min}_{\{j=1,2,3,4;l=0,1,2,\cdots \}}\left\{{\zeta }_j^{\left(l\right)}\right\}$ and then if $\rho ={\rho }_{0★}$, (4.5) has a pair of imaginary roots $\pm i{\zeta }_0$.

Now we give the assumption as follows:

$$\left(M_4\right){\mathcal{B}}_{1R}{\mathcal{B}}_{2R}+{\mathcal{B}}_{1I}{\mathcal{B}}_{2I}>0,$$

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{B}}_{1R}={\tilde{b}}_1+{\tilde{b}}_{3}cos{\zeta }_0{\rho }_{0★},}\hfill \\ {\displaystyle {\mathcal{B}}_{1I}=2{\zeta }_0-{\tilde{b}}_{3}sin{\zeta }_0{\rho }_{0★},}\hfill \\ {\displaystyle {\mathcal{B}}_{2R}={\tilde{b}}_4{\zeta }_{0}sin{\zeta }_0{\rho }_{0★}-{\tilde{b}}_3{\zeta }_{0}cos{\zeta }_0{\rho }_{0★},}\hfill \\ {\displaystyle {\mathcal{B}}_{2I}={\tilde{b}}_4{\zeta }_{0}cos{\zeta }_0{\rho }_{0★}+{\tilde{b}}_3{\zeta }_{0}sin{\zeta }_0{\rho }_{0★}.}\hfill \end{array}\right.$$

(48)

Lemma 2. 

If $\lambda \left(\rho \right)={\tau }_1\left(\rho \right)+i{\tau }_2\left(\rho \right)$ is a root of Eq. (4.5) at $\rho ={\rho }_{0★}$ obeying ${\tau }_1\left({\rho }_{0★}\right)=0,{\tau }_2\left({\rho }_{0★}\right)={\zeta }_0$, then $\mathrm{Re}\left({\displaystyle \frac{d\lambda }{d\rho }}\right){|}_{\rho ={\rho }_{0★},\zeta ={\zeta }_0}>0.$

Proof. 

Based on Equation (4.5), one gains

$$(2\lambda +{\tilde{b}}_1){\displaystyle \frac{d\lambda }{d\rho }}+{\tilde{b}}_3{\displaystyle \frac{d\lambda }{d\rho }}{e}^{-\lambda \rho }-e^{-\lambda \rho }\left({\displaystyle \frac{d\lambda }{d\rho }}\rho +\lambda \right)\left({\tilde{b}}_3\lambda +{\tilde{b}}_4\right)=0,$$

(49)

which generates

$${\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}={\displaystyle \frac{{\mathcal{B}}_1\left(\lambda \right)}{{\mathcal{B}}_2\left(\lambda \right)}}-{\displaystyle \frac{\rho }{\lambda }},$$

(50)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{B}}_1\left(\lambda \right)=2\lambda +{\tilde{b}}_1+{\tilde{b}}_3{e}^{-\lambda \rho },}\hfill \\ {\displaystyle {\mathcal{B}}_2\left(\lambda \right)=\lambda ({\tilde{b}}_3\lambda +{\tilde{b}}_4)e^{-\lambda \vartheta }.}\hfill \end{array}\right.$$

(51)

Then

$$\mathrm{Re}{\left[{\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}\right]}_{\rho ={\rho }_{0★},\zeta ={\zeta }_0}=\mathrm{Re}{\left[{\displaystyle \frac{{\mathcal{B}}_1\left(\lambda \right)}{{\mathcal{B}}_2\left(\lambda \right)}}\right]}_{\rho ={\rho }_{0★},\zeta ={\zeta }_0}={\displaystyle \frac{{\mathcal{B}}_{1R}{\mathcal{B}}_{2R}+{\mathcal{B}}_{1I}{\mathcal{B}}_{2I}}{{\mathcal{B}}_{2R}^2+{\mathcal{B}}_{2I}^2}}.$$

(52)

In view of $\left(M_4\right)$, one gains

$$\mathrm{Re}{\left[{\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}\right]}_{\rho ={\rho }_{0★},\zeta ={\zeta }_0}>0.$$

(53)

The proof finishes. □

Relying on the analysis above, the following outcome is lightly acquired.

Lemma 3. 

If $\left(M_3\right)$ and $M_4)$ are fulfilled, then the equilibrium point $E(u_{1★},u_{2★})$ of model (4.1) keeps locally asymptotically stable condition if $\rho \in [0,{\rho }_{0★})$ and model (4.1) produces a cluster of periodic solutions (namely, Hopf bifurcations) near the equilibrium point $E(u_{1★},u_{2★})$ when $\rho ={\rho }_{0★}.$

5. Bifurcation Control via Extended Hybrid Delayed Feedback Controller

In this section, we are to use two hybrid controllers (include a state feedback and parameter perturbation) to the both equations of model (1.4). By virtue of the two controllers, we can dominate the time of emergence of bifurcation phenomenon and stability region of system (1.4). Relying on the ideas of [25-29], we formulate the following fractional controlled nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}={\sigma }_1\left[\frac{c{u}_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+\mathcal{K}}-c{u}_1\left(t\right)\right]+{\sigma }_2[u_1(t-\rho )-u_1\left(t\right)],}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}={\sigma }_3[\beta -u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+\mathcal{K}}]+{\sigma }_4[u_2(t-\rho )-u_2\left(t\right)],}\hfill \end{array}\right.$$

(54)

where ${\sigma }_l(l=1,2,3,4)$ represent feedback gain parameters. Clearly, system (5.1) and system (1.5) admit the mutual equilibrium point $E(u_{1★},u_{2★})$. The linear equation for system (5.1) at $E(u_{1★},u_{2★})$ is given by

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}={\widehat{a}}_1{u}_1\left(t\right)+{\widehat{a}}_2{u}_2\left(t\right)+{\widehat{a}}_3{u}_1(t-\rho ),}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}={\widehat{a}}_4{u}_1\left(t\right)+{\widehat{a}}_5{u}_2\left(t\right)+{\widehat{a}}_6{u}_2(t-\rho ),}\hfill \end{array}\right.$$

(55)

where

$$\left\{\begin{array}{c}{\displaystyle {\widehat{a}}_1=-{\sigma }_{1}c-{\sigma }_2,}\hfill \\ {\displaystyle {\widehat{a}}_2=\frac{{\sigma }_{1}c{u}_{1★}^2}{u_{1★}+\mathcal{K}},}\hfill \\ {\displaystyle {\widehat{a}}_3=\frac{2{\sigma }_{1}c{u}_{1★}{u}_{2★}}{u_{1★}+\mathcal{K}}-\frac{2{\sigma }_{1}c{u}_{1★}^2{u}_{2★}}{{(u_{1★}+\mathcal{K})}^2}+{\sigma }_2,}\hfill \\ {\displaystyle {\widehat{a}}_4=\frac{2{\sigma }_2{u}_{1★}{u}_{2★}}{\mathcal{K}}-\frac{{\sigma }_3{u}_{1★}^2}{{\mathcal{K}}^2},}\hfill \\ {\displaystyle {\widehat{a}}_5=\frac{{\sigma }_3{u}_{1★}^2}{\mathcal{K}}-{\sigma }_4,}\hfill \\ {\displaystyle {\widehat{a}}_6={\sigma }_4.}\hfill \end{array}\right.$$

(56)

The characteristic equation of system (5.2) reads as

$$det\left[\begin{array}{cc}\lambda -{\widehat{a}}_1-{\widehat{a}}_3{e}^{-\lambda \rho }& -{\widehat{a}}_2\\ -{\widehat{a}}_4& \lambda -{\widehat{a}}_5-{\widehat{a}}_6{e}^{-\lambda \rho }\end{array}\right]=0,$$

(57)

which results in

$${\lambda }^2+{\widehat{e}}_1\lambda +{\widehat{e}}_2+({\widehat{e}}_3\lambda +{\widehat{e}}_4)e^{-\lambda \rho }+{\widehat{e}}_5{e}^{-2\lambda \rho }=0,$$

(58)

where

$$\left\{\begin{array}{c}{\displaystyle {\widehat{e}}_1=-({\widehat{a}}_1+{\widehat{a}}_3),}\hfill \\ {\displaystyle {\widehat{e}}_2={\widehat{a}}_1{\widehat{a}}_5-{\widehat{a}}_2{\widehat{a}}_4,}\hfill \\ {\displaystyle {\widehat{e}}_3=-({\widehat{a}}_3+{\widehat{a}}_6),}\hfill \\ {\displaystyle {\widehat{e}}_4={\widehat{a}}_1{\widehat{a}}_6+{\widehat{a}}_3{\widehat{a}}_5,}\hfill \\ {\displaystyle {\widehat{e}}_5={\widehat{a}}_3{\widehat{a}}_6.}\hfill \end{array}\right.$$

(59)

When $\rho =0,$ then Equation (5.5) owns the following expression:

$${\lambda }^2+({\widehat{e}}_1+{\widehat{e}}_3)\lambda +{\widehat{e}}_2+{\widehat{e}}_4+{\widehat{e}}_5=0.$$

(60)

If

$$\left(M_5\right){\widehat{e}}_1+{\widehat{e}}_3>0,{\widehat{e}}_2+{\widehat{e}}_4+{\widehat{e}}_5>0$$

is fulfilled, then the two roots ${\lambda }_1,{\lambda }_2$ of Eq. (5.7) admit negative real parts. So one can conclude that the equilibrium point $E(u_{1★},u_{2★})$ of system (5.1) under the time delay $\rho =0$ keeps locally asymptotically stable level.

By (5.5), we have

$$({\lambda }^2+{\widehat{e}}_1\lambda +{\widehat{e}}_2)e^{\lambda \rho }+({\widehat{e}}_3\lambda +{\widehat{e}}_4)+{\widehat{e}}_5{e}^{-\lambda \rho }=0,$$

(61)

Let $\lambda =i\eta$ be the root of Eq. (5.8). By Eq. (5.8), we get

$$[{\left(i\eta \right)}^2+{\widehat{e}}_{1}i\eta +{\widehat{e}}_2]e^{i\eta \rho }+({\widehat{e}}_{3}i\eta +{\widehat{e}}_4)+{\widehat{e}}_5{e}^{-i\eta \rho }=0,$$

(62)

which results in

$$\left\{\begin{array}{c}{\displaystyle ({\widehat{e}}_2-{\eta }^2+{\widehat{e}}_5)cos\eta \rho -{\widehat{e}}_1\eta sin\eta \rho =-{\widehat{e}}_4,}\hfill \\ {\displaystyle {\widehat{e}}_1\eta cos\eta \rho +({\widehat{e}}_2-{\eta }^2-{\widehat{e}}_5)sin\eta \rho =-{\widehat{e}}_3\eta }\hfill \end{array}\right.$$

(63)

By (5.10), one gains

$$\left\{\begin{array}{c}{\displaystyle sin\eta \rho =\frac{-{\widehat{e}}_3\eta ({\widehat{e}}_2-{\eta }^2+{\widehat{e}}_5)-{\widehat{e}}_1{\widehat{e}}_4\eta }{{({\widehat{e}}_2-{\eta }^2)}^2-{\widehat{e}}_5^2+{\widehat{e}}_1^2{\eta }^2},}\hfill \\ {\displaystyle cos\eta \rho =\frac{-{\widehat{e}}_4({\widehat{e}}_2-{\eta }^2-{\widehat{e}}_5)-{\widehat{e}}_1{\widehat{e}}_3{\eta }^2}{{({\widehat{e}}_2-{\eta }^2)}^2-{\widehat{e}}_5^2+{\widehat{e}}_1^2{\eta }^2}.}\hfill \end{array}\right.$$

(64)

In view of ${sin}^2\eta \rho +{cos}^2\eta \rho =1$, we get

$$\begin{array}{ccc}& & {[{\widehat{e}}_3\eta ({\widehat{e}}_2-{\eta }^2+{\widehat{e}}_5)+{\widehat{e}}_1{\widehat{e}}_4\eta ]}^2\hfill \\ & & +{[{\widehat{e}}_4({\widehat{e}}_2-{\eta }^2-{\widehat{e}}_5)+{\widehat{e}}_1{\widehat{e}}_3{\eta }^2]}^2\hfill \\ & & ={[{({\widehat{e}}_2-{\eta }^2)}^2-{\widehat{e}}_5^2+{\widehat{e}}_1^2{\eta }^2]}^2,\hfill \end{array}$$

(65)

which results in

$${\eta }^8+f_1{\eta }^6+f_2{\eta }^4+f_3{\eta }^2+f_4=0,$$

(66)

where

$$\left\{\begin{array}{c}{\displaystyle f_1=2({\widehat{e}}_1^2-2{\widehat{e}}_2)+{\widehat{e}}_3,}\hfill \\ {\displaystyle f_2=2({\widehat{e}}_2^2-2{\widehat{e}}_5^2)+{({\widehat{e}}_1^2-2{\widehat{e}}_2)}^2+{({\widehat{e}}_1{\widehat{e}}_3-{\widehat{e}}_4)}^2,}\hfill \\ {\displaystyle +2{\widehat{e}}_3({\widehat{e}}_1{\widehat{e}}_4+{\widehat{e}}_2{\widehat{e}}_3+{\widehat{e}}_3{\widehat{e}}_5),}\hfill \\ {\displaystyle f_3=2({\widehat{e}}_2^2-2{\widehat{e}}_5^2)({\widehat{e}}_1^2-2{\widehat{e}}_2)-{({\widehat{e}}_1{\widehat{e}}_4+{\widehat{e}}_2{\widehat{e}}_3+{\widehat{e}}_3{\widehat{e}}_5)}^2}\hfill \\ {\displaystyle -2({\widehat{e}}_2{\widehat{e}}_4-{\widehat{e}}_4{\widehat{e}}_5)({\widehat{e}}_1{\widehat{e}}_3-{\widehat{e}}_4),}\hfill \\ {\displaystyle f_4={({\widehat{e}}_2^2-{\widehat{e}}_5^2)}^2-{({\widehat{e}}_2{\widehat{e}}_4-{\widehat{e}}_4{\widehat{e}}_5)}^2.}\hfill \end{array}\right.$$

(67)

Define

$${\Psi }_3\left(\eta \right)={\eta }^8+f_1{\eta }^6+f_2{\eta }^4+f_3{\eta }^2+f_4.$$

(68)

Suppose that

$$\left(M_6\right)|{\widehat{e}}_2^2-{\widehat{e}}_5^2|<|{\widehat{e}}_2{\widehat{e}}_4-{\widehat{e}}_4{\widehat{e}}_5|$$

holds, notice that ${lim}_{\eta \to +\infty }{\Psi }_3\left(\eta \right)=+\infty >0$, then one understands that Eq. (5.13) admits at least one real positive root. So Equation (5.5) owns at least one pair of purely roots. Without loss of generality, assume that Eq. (5.13) admits eight real positive roots (say ${\eta }_l,l=1,2,3\cdots ,8)$. According to (5.11), we gain

$${\rho }_l^{\left(i\right)}={\displaystyle \frac{1}{{\eta }_l}}\left[arccos\left({\displaystyle \frac{-{\widehat{e}}_4({\widehat{e}}_2-{\eta }_l^2-{\widehat{e}}_5)-{\widehat{e}}_1{\widehat{e}}_3{\eta }_l^2}{{({\widehat{e}}_2-{\eta }_l^2)}^2-{\widehat{e}}_5^2+{\widehat{e}}_1^2{\eta }_l^2}}\right)+2i\pi \right],$$

(69)

where $l=1,2,3\cdots ,8;i=0,1,2,\cdots .$ Let ${\rho }_{0★★}={min}_{\{l=1,2,\cdots ,12;i=0,1,2,\cdots \}}\left\{{\rho }_l^{\left(i\right)}\right\}$ and assume that when $\rho ={\rho }_{0★★}$, (5.5) admits a pair of imaginary roots $\pm i{\eta }_0$.

Next the following condition is given:

$$\left(M_7\right){\mathcal{C}}_{1R}{\mathcal{C}}_{2R}+{\mathcal{C}}_{1I}{\mathcal{C}}_{2I}>0,$$

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{C}}_{1R}={\widehat{e}}_1+{\widehat{e}}_{3}cos{\eta }_0{\rho }_{0★★},}\hfill \\ {\displaystyle {\mathcal{C}}_{1I}=2{\eta }_0-{\widehat{e}}_{3}sin{\eta }_0{\rho }_{0★★},}\hfill \\ {\displaystyle {\mathcal{C}}_{2R}={\widehat{e}}_4{\eta }_{0}sin{\eta }_0{\rho }_{0★★}-{\widehat{e}}_3{\eta }_0^{2}cos{\eta }_0{\rho }_{0★★}+2{\widehat{e}}_5{\eta }_{0}sin2{\eta }_0{\rho }_{0★★},}\hfill \\ {\displaystyle {\mathcal{C}}_{2I}={\widehat{e}}_2{\eta }_0^{2}sin{\eta }_0{\rho }_{0★★}+{\widehat{e}}_4{\eta }_{0}cos{\eta }_0{\rho }_{0★★}+2{\widehat{e}}_5{\eta }_{0}cos2{\eta }_0{\rho }_{0★★}.}\hfill \end{array}\right.$$

(70)

Lemma 4. Lemma 5.1. 

If $\lambda \left(\vartheta \right)={\varphi }_1\left(\rho \right)+i{\varphi }_2\left(\rho \right)$ is a root of Eq. (5.5) near $\rho ={\rho }_{0★★}$ that obeys ${\varphi }_1\left({\rho }_{0★★}\right)=0,{\varphi }_2\left({\rho }_{0★★}\right)={\eta }_0$, then $\mathrm{Re}\left({\displaystyle \frac{d\lambda }{d\vartheta }}\right){|}_{\rho ={\rho }_{0★★},\eta ={\eta }_0}>0.$

Proof. 

By virtue of Equation (5.5), one gains

$$\begin{array}{ccc}& & (2\lambda +{\widehat{e}}_1){\displaystyle \frac{d\lambda }{d\rho }}+{\widehat{e}}_3{\displaystyle \frac{d\lambda }{d\rho }}{e}^{-\lambda \rho }-e^{-\lambda \rho }\left({\displaystyle \frac{d\lambda }{d\rho }}\rho +\lambda \right)\left({\widehat{e}}_3\lambda +{\widehat{e}}_4\right)\hfill \\ & & -2{\widehat{e}}_5{e}^{-2\lambda \rho }\left({\displaystyle \frac{d\lambda }{d\rho }}\rho +\lambda \right)=0,\hfill \end{array}$$

(71)

which results in

$${\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}={\displaystyle \frac{{\mathcal{C}}_1\left(\lambda \right)}{{\mathcal{C}}_2\left(\lambda \right)}}-{\displaystyle \frac{\rho }{\lambda }},$$

(72)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{C}}_1\left(\lambda \right)=2\lambda +{\widehat{e}}_1+{\widehat{e}}_3{e}^{-\lambda \rho },}\hfill \\ {\displaystyle {\mathcal{C}}_2\left(\lambda \right)=e^{-\lambda \rho }\lambda ({\widehat{e}}_3\lambda +{\widehat{e}}_4)+2{\widehat{e}}_5\lambda e^{-2\lambda \rho }.}\hfill \end{array}\right.$$

(73)

Then

$$\mathrm{Re}{\left[{\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}\right]}_{\rho ={\rho }_{0★★},\eta ={\eta }_0}=\mathrm{Re}{\left[{\displaystyle \frac{{\mathcal{C}}_1\left(\lambda \right)}{{\mathcal{C}}_2\left(\lambda \right)}}\right]}_{\rho ={\rho }_{0★★},\eta ={\eta }_0}={\displaystyle \frac{{\mathcal{C}}_{1R}{\mathcal{C}}_{2R}+{\mathcal{C}}_{1I}{\mathcal{C}}_{2I}}{{\mathcal{C}}_{2R}^2+{\mathcal{C}}_{2I}^2}}.$$

(74)

In view of $\left(M_7\right)$, one gains

$$\mathrm{Re}{\left[{\left({\displaystyle \frac{d\lambda }{d\rho }}\right)}^{-1}\right]}_{\rho ={\rho }_{0★★},\eta ={\eta }_0}>0.$$

(75)

The proof completes. □

Depending on the discussion above, the following outcome is lightly acquired.

Theorem 5. 

If $\left(M_5\right)$-$\left(M_7\right)$ are satisfied, then the equilibrium point $E(u_{1★},u_{2★})$ of system (5.1) keeps locally asymptotically stable level if $\rho \in [0,{\rho }_{0★★})$ and model (5.1) produces a cluster of periodic solutions (i.e., Hopf bifurcations) near the equilibrium point $E(u_{1★},u_{2★})$ when $\rho ={\rho }_{0★★}.$

6. Numerical Experiments

Example 1. 

Consider the following delayed nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}=\frac{5u_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+0.05}-5u_1\left(t\right),}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}=8-u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+0.05}.}\hfill \end{array}\right.$$

(76)

One can lightly acquire that model (6.1) owns the positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$. By a series of algebraic operation, we can verify that the three assumptions $\left(M_1\right),\left(M_2\right)$ of Theorem 3.1 hold true. Using Matlab software, one gains ${\rho }_0=0.45$. In order to check the onset of Hopf bifurcation and stability issue of Theorem 3.1, both unequal delay values are provided. On the one hand, we select $\rho =0.42$ which demonstrates $\rho =0.42<{\rho }_0=0.45$, namely, ρ is located within the range $[0,{\rho }_0)$. On the other hand, we select $\rho =0.48$ which demonstrates $\rho =0.48>{\rho }_0=0.45$, namely, ρ is located outside of the range $[0,{\rho }_0)$. For $\rho =0.42<{\rho }_0=0.45$, we obtain the simulation Figure 1. Based on Figure 1, we can lightly see that the positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$ is locally stable. In other words, the density of microorganism and the density of one chemical species in the sediment will run towards $6.9928,1.0072$ respectively. For $\rho =0.48>{\rho }_0=0.45$, we obtain the simulation Figure 2. Based on Figure 2, we can lightly see that a family of periodic solutions (i.e, Hopf bifurcations) will take place around positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$. In other words, the density of microorganism and the density of one chemical species in the sediment will retain periodic vibration around the values $6.9928,1.0072$ respectively. In addition, to show visually the bifurcation point, we draw the corresponding bifurcation diagrams which are given in Figure 3 and Figure 4. One can know that the bifurcation value is about $0.45$.

Example 2. 

Consider the following controlled delayed nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}=\sigma \left[\frac{5u_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+0.05}-5u_1\left(t\right)\right]+(1-\sigma )[u_1(t-\rho )-u_1\left(t\right)],}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}=8-u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+0.05}.}\hfill \end{array}\right.$$

(77)

One can lightly acquire that model (6.2) owns the positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$. Let $\sigma =0.7$. By a series of algebraic operation, we can verify that the three assumptions $\left(M_3\right),\left(M_4\right)$ of Theorem 4.1 hold true. Using Matlab software, one gains ${\rho }_{0★}=0.84$. In order to check the onset of Hopf bifurcation and stability issue of Theorem 4.1, both unequal delay values are provided. On the one hand, we select $\rho =0.8$ which demonstrates $\rho =0.8<{\rho }_{0★}=0.84$, namely, ρ is located within the range $[0,{\rho }_{0★})$. On the other hand, we select $\rho =0.93$ which demonstrates $\rho =0.93>{\rho }_{0★}=0.84$, namely, ρ is located outside of the range $[0,{\rho }_{0★})$. For $\rho =0.8<{\rho }_{0★}=0.84$, we obtain the simulation Figure 5. Based on Figure 5, we can lightly see that the positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$ is locally stable. In other words, the density of microorganism and the density of one chemical species in the sediment will run towards $6.9928,1.0072$ respectively. For $\rho =0.93>{\rho }_{0★}=0.84$, we obtain the simulation Figure 6. Based on Figure 6, we can lightly see that a family of periodic solutions (i.e, Hopf bifurcations) will take place around positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$. In other words, the density of microorganism and the density of one chemical species in the sediment will retain periodic vibration around the values $6.9928,1.0072$ respectively. In addition, to show visually the bifurcation point, we draw the corresponding bifurcation diagrams which are given in Figure 7 and Figure 8. One can know that the bifurcation value is about $0.84$.

Example 3. 

Consider the following controlled delayed nutrient-microorganism system:

$$\left\{\begin{array}{c}{\displaystyle \frac{d{u}_1\left(t\right)}{dt}={\sigma }_1\left[\frac{5u_1^2(t-\rho )u_2\left(t\right)}{u_1(t-\rho )+0.05}-5u_1\left(t\right)\right]+{\sigma }_2[u_1(t-\rho )-u_1\left(t\right)],}\hfill \\ {\displaystyle \frac{d{u}_2\left(t\right)}{dt}={\sigma }_3\left[8-u_2\left(t\right)-\frac{u_1^2\left(t\right)u_2\left(t\right)}{u_1\left(t\right)+0.05}\right]+{\sigma }_4[u_2(t-\rho )-u_2\left(t\right)].}\hfill \end{array}\right.$$

(78)

One can lightly acquire that model (6.3) owns the positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$. Let ${\sigma }_1=0.5,{\sigma }_2=0.5,{\sigma }_3=0.7,{\sigma }_4=-0.3$. By a series of algebraic operation, we can verify that the three assumptions $\left(M_5\right),\left(M_6\right),\left(M_7\right)$ of Theorem 5.1 hold true. Using Matlab software, one gains ${\rho }_{0★★}=0.62$. In order to check the onset of Hopf bifurcation and stability issue of Theorem 5.1, both unequal delay values are provided. On the one hand, we select $\rho =0.57$ which demonstrates $\rho =0.57<{\rho }_{0★★}=0.62$, namely, ρ is located within the range $[0,{\rho }_{0★★})$. On the other hand, we select $\rho =0.68$ which demonstrates $\rho =0.68>{\rho }_{0★★}=0.62$, namely, ρ is located outside of the range $[0,{\rho }_{0★★})$. For $\rho =0.57<{\rho }_{0★}=0.62$, we obtain the simulation Figure 9. Based on Figure 9, we can lightly see that the positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$ is locally stable. In other words, the density of microorganism and the density of one chemical species in the sediment will run towards $6.9928,1.0072$ respectively. For $\rho =0.68>{\rho }_{0★★}=0.62$, we obtain the simulation Figure 10. Based on Figure 10, we can lightly see that a family of periodic solutions (i.e, Hopf bifurcations) will take place around positive equilibrium point $E(u_{1★},u_{2★})=E(6.9928,1.0072)$. In other words, the density of microorganism and the density of one chemical species in the sediment will retain periodic vibration around the values $6.9928,1.0072$ respectively. In addition, to show visually the bifurcation point, we draw the corresponding bifurcation diagrams which are given in Figure 11 and Figure 12. One can know that the bifurcation value is about $0.62$.

Remark 1. 

Relying on the numerical simulation outcomes in Example 6.1-Example 6.3, we can find that system (6.1) owns the bifurcation value $0.45$, system (6.2) owns the bifurcation value $0.84$, system (6.3) owns the bifurcation value $0.62$, which confirms that we can postpone the time of appearance of bifurcation phenomenon of system (6.1) and expand the stability range of system (6.1) via our designed both different hybrid delayed feedback controllers in varying degrees.

7. Conclusions

During the past decades, mathematical models play a vital role in depicting the evolution of various chemical compositions in chemical reactions. Various chemical reaction models have been formulated and a lot of dynamical behavior of these chemical reaction models have been studied. In the manuscript, based on the forefathers’ research, we build new delayed nutrient-microorganism system. The parameter conditions on the existence and uniqueness, non-negativeness and boundedness of the solution of this model are acquired via fixed point theorem, inequality technique and a reasonable function. The stability and bifurcation phenomenon are discussed by virtue of the stability criterion and bifurcation argument of fractional delayed dynamical system. The critical delay value to guarantee the onset of Hopf bifurcation of this model is acquired. In order to control the the time of emergence of bifurcation phenomenon and stability domain of this model, two different hybrid controllers are provided to achieve the our goal. In these controllers, we can find that delay and control parameters are very important factors in achieve the our goal. The research fruits own enormous theoretical significance in dominating and optimizing the densities of microorganism and nutrient. Also the research fruits can be applied to control the stability, bifurcation behavior and chaos in many differential dynamical models (include integer-order and fractional-order dynamical models) in lots of areas.