Energy and chemical efficient nitrogen removal at a full-scale MBR water reuse facility

Jianfeng Wen; Yanjin Liu; Yunjie Tu and Mark W. LeChevallier; Jianfeng Wen; Yanjin Liu; Yunjie Tu and Mark W. LeChevallier

doi:10.3934/environsci.2015.1.42

AIMS Environmental Science

2015, Volume 2, Issue 1: 42-55. doi: 10.3934/environsci.2015.1.42

Previous Article Next Article

Research article Special Issues

Energy and chemical efficient nitrogen removal at a full-scale MBR water reuse facility

American Water, 1025 Laurel Oak Road, Voorhees, New Jersey 08043, USA

Received: 10 December 2014 Accepted: 12 February 2015 Published: 14 February 2015

With stringent wastewater discharge limits on nitrogen and phosphorus, membrane bioreactor (MBR) technology is gaining popularity for advanced wastewater treatment due to higher effluent quality and smaller footprint. However, higher energy intensity required for MBR plants and increased operational costs for nutrient removal limit wide application of the MBR technology. Conventional nitrogen removal requires intensive energy inputs and chemical addition. There are drivers to search for new technology and process control strategies to treat wastewater with lower energy and chemical demand while still producing high quality effluent. The NPXpress is a patented technology developed by American Water engineers. This technology is an ultra-low dissolved oxygen (DO) operation for wastewater treatment and is able to remove nitrogen with less oxygen requirements and reduced supplemental carbon addition in MBR plants. Jefferson Peaks Water Reuse Facility in New Jersey employs MBR technology to treat municipal wastewater and was selected for the implementation of the NPXpress technology. The technology has been proved to consistently produce a high quality reuse effluent while reducing energy consumption and supplemental carbon addition by 59% and 100%, respectively. Lab-scale kinetic studies suggested that NPXpress promoted microorganisms with higher oxygen affinity. Process modelling was used to simulate treatment performance under NPXpress conditions and develop ammonia-based aeration control strategy. The application of the ammonia-based aeration control at the plant further reduced energy consumption by additional 9% and improved treatment performance with 35% reduction in effluent total nitrogen. The overall energy savings for Jefferson Peaks was $210,000 in four years since the implementation of NPXpress. This study provided an insight in design and operation of MBR plants with NPXpress technology and ultra-low DO operations.

Keywords:

Citation: Jianfeng Wen, Yanjin Liu, Yunjie Tu and Mark W. LeChevallier. Energy and chemical efficient nitrogen removal at a full-scale MBR water reuse facility[J]. AIMS Environmental Science, 2015, 2(1): 42-55. doi: 10.3934/environsci.2015.1.42

Related Papers:

[1]	Xiaoming Su, Jiahui Wang, Adiya Bao . Stability analysis and chaos control in a discrete predator-prey system with Allee effect, fear effect, and refuge. AIMS Mathematics, 2024, 9(5): 13462-13491. doi: 10.3934/math.2024656
[2]	Xiongxiong Du, Xiaoling Han, Ceyu Lei . Dynamics of a nonlinear discrete predator-prey system with fear effect. AIMS Mathematics, 2023, 8(10): 23953-23973. doi: 10.3934/math.20231221
[3]	A. Q. Khan, Ibraheem M. Alsulami . Complicate dynamical analysis of a discrete predator-prey model with a prey refuge. AIMS Mathematics, 2023, 8(7): 15035-15057. doi: 10.3934/math.2023768
[4]	Binfeng Xie, Na Zhang . Influence of fear effect on a Holling type III prey-predator system with the prey refuge. AIMS Mathematics, 2022, 7(2): 1811-1830. doi: 10.3934/math.2022104
[5]	Figen Kangalgil, Seval Isșık . Effect of immigration in a predator-prey system: Stability, bifurcation and chaos. AIMS Mathematics, 2022, 7(8): 14354-14375. doi: 10.3934/math.2022791
[6]	Xuyang Cao, Qinglong Wang, Jie Liu . Hopf bifurcation in a predator-prey model under fuzzy parameters involving prey refuge and fear effects. AIMS Mathematics, 2024, 9(9): 23945-23970. doi: 10.3934/math.20241164
[7]	Ali Yousef, Ashraf Adnan Thirthar, Abdesslem Larmani Alaoui, Prabir Panja, Thabet Abdeljawad . The hunting cooperation of a predator under two prey's competition and fear-effect in the prey-predator fractional-order model. AIMS Mathematics, 2022, 7(4): 5463-5479. doi: 10.3934/math.2022303
[8]	San-Xing Wu, Xin-You Meng . Dynamics of a delayed predator-prey system with fear effect, herd behavior and disease in the susceptible prey. AIMS Mathematics, 2021, 6(4): 3654-3685. doi: 10.3934/math.2021218
[9]	Kottakkaran Sooppy Nisar, G Ranjith Kumar, K Ramesh . The study on the complex nature of a predator-prey model with fractional-order derivatives incorporating refuge and nonlinear prey harvesting. AIMS Mathematics, 2024, 9(5): 13492-13507. doi: 10.3934/math.2024657
[10]	Kwadwo Antwi-Fordjour, Rana D. Parshad, Hannah E. Thompson, Stephanie B. Westaway . Fear-driven extinction and (de)stabilization in a predator-prey model incorporating prey herd behavior and mutual interference. AIMS Mathematics, 2023, 8(2): 3353-3377. doi: 10.3934/math.2023173

Abstract

1. Introduction

The use of mathematical models to study the interaction between different biological elements is an important area in ecosystems. The main purpose of the model is to explore the dynamical behaviors between them, such as predation, symbiosis, parasitism, competition, and so on ^[1,2,3]. Since Lotka ^[4] and Volterra ^[5] separately proposed the Lotka–Volterra model to describe the ecological interaction by differential equation model, many mathematicians and biologists have been deeply interested in studying such models, especially the predator–prey models. In biology and biomathematics, the predation model is one of the important research topics, and it is also an important model of ecosystems ^[6]. There is a large amount of literature on predator–prey models with different types of functional responses. Holling ^[7,8,9] came up with three types of functional response functions to describe the relationship between the quantity of prey captured by predators and the density of prey. Taking into account the living environment and situation of prey and predators in reality, the dwellings of many prey have a certain protective effect, and their ability to resist natural enemies will be enhanced with the increase of the density of prey, or prey will learn to hide themselves under the pursuit of natural enemies, which will reduce the predation ability of natural enemies. Therefore, Andrews ^[10] proposed the Holling–IV functional response function in 1968. Besides, Ajraldi ^[11] proposed a square–root functional response to explain the aggregation behavior that plankton exhibit.

Recently, the discrete population models depicted by differential equations have been extensively studied ^[12,13,14]. Discrete-time systems, compared with continuous–time versions, have obvious advantages. On one side, data collection for biological samples is often on discrete time scales in reality (either in a week, a month, or another certain time period) rather than on a day. On the other side, assume that the population quantity of species does not overlap between successive generations; these could offer more efficient calculation results for numerical simulations ^[15,16]. Therefore, a lot of mathematicians have studied the dynamical behaviors of the discrete-time model corresponding to the continuous–time one. For instance, Agarwal ^[17] formulated bifurcation analysis of a discrete predator–prey model incorporating prey refuge. Ahmed et al. ^[18] explored a class of discrete prey predation models with fast–slow effects. Khan et al. ^[19] studied a class of discrete Rosenzweig–Macarthur predator–prey models and analysed their chaos.

Since the fear effect was first put forward by Wang et al. ^[20], the research of models with the fear effect has been paid more attention by many scholars ^[21,22,23]. The fear effect mainly refers to the mutual relationship between prey and predators, which makes prey populations instinctively fear predators, thereby reducing the birth rate of prey and thus achieving the purpose of reducing prey capture. It is obvious that this is a kind of anti–predation behavior. In ^[20], the expression for the fear effect is $ F(x, y) = 1/(1+ky) $, where the parameter $ k $ represents the level of the fear effect, and the prey population decreases as $ k $ goes up. Compared with the fear effect, the refuge effect is also an anti–predation behavior. It is important to note that the prey cost most of their lives nearby or hiding in shelters for avoiding predation, such as caves, crevices, dense vegetation, shells, or pipes. The notion of prey sanctuary has caught ecologists and mathematicians' attention since Gause et al. ^[24] and Smith ^[25] brought in the parameter refer to refuge. In the ecology, prey refuge could reduce the predation of the prey by the predator, thereby refraining from extinction of the prey population due to the predation factor; see ^{[26,27,28,29,30]} in detail. Many scholars studied a type of continuous predator–prey model with the fear effect or refuge effect ^{[31,32,33,34,35,36]}. For example, Pal et al. ^[32] added the fear effect, the Allee effect, and the influence of external disturbances such as refuges into a continuous predator–prey model. At the same time, the positivity, boundedness, stability of the equilibrium points, and bifurcation phenomena of the continuous predator–prey model were analyzed. Wang et al. ^[34] mainly considered the role of the fear effect and discussed the difference in stability with and without the fear effect; besides, they also studied some bifurcation phenomena. Consider that some populations without overlapping generations are not suitable for continuous models, so we can establish a type of discrete predator–prey model with a fear effect and a prey refuge effect. First of all, a continuous version is put forward as follows:

$\begin{align} \left\{ \begin{array}{l} \frac{dx}{dt} = \frac{rx}{1+ky}(1-\frac{x}{K})-c(x-R)y, \\ \frac{dy}{dt} = y-dy+e(x-R). \end{array} \right. \end{align}$ $

(1.1)

The discretization system (1.1) is obtained

$\begin{align} \left\{ \begin{array}{l} x_{n+1} = x_n+\frac{rx_n}{1+ky_n}(1-\frac{x_n}{K})-c(x_n-R)y_n, \\ y_{n+1} = y_n-dy_n+e(x_n-R)y_n. \end{array} \right. \end{align}$ $

(1.2)

Where $ x_n $ and $ y_n $ are the population densities at $ n $th generation of prey and predator, respectively. $ r, k, K, c, d, e $ and $ R $ severally stand for the growth rate of prey, the level of the fear effect, the carrying capacity of prey, the capture efficiency of the predator of prey, the natural mortality rate of the predator, the rate of conversion of energy, and the quantity of prey using refuge. All of the above parameters are strictly positive.

The dramatic increase in people's demand for resources has led to biological resources being over-exploited. It needs us to maintain sustainable development and maximize economic benefits for the exploitation of ecological resources. Current research refers to the effect of harvesting and mainly focuses on constant harvesting ^[37], proportional harvesting ^[38], and nonlinear harvesting ^[39]. We introduce the type of proportional harvesting of prey and predator. Our system takes the following form:

$\begin{align} \left\{ \begin{array}{l} x_{n+1} = x_n+\frac{rx_n}{1+ky_n}(1-\frac{x_n}{K})-c(x_n-R)y_n-q_1E_1x_n, \\ y_{n+1} = y_n-dy_n+e(x_n-R)y_n-q_2E_2y_n. \end{array} \right. \end{align}$ $

(1.3)

Here $ q_1, q_2 $ represent the catchability coefficients of the prey and predator, and $ E_1, E_2 $ mean the harvesting efforts of the prey and predator, respectively.

The outline of the rest of this paper is in the following manner: Sections 2 and 3 explore the existence and stability of all possible fixed points of system (1.3), respectively. In Sections 4 and 5, we investigate in detail two types of bifurcation analysis (Flip bifurcation and Neimark–Sacker bifurcation) at the interior fixed point of system (1.3). Also, the optimal harvesting strategy is analysed in Section 6. Numerical simulations are applied in Section 7 to illustrate our theoretical results analysed above. Finally, Section 8 contains the conclusion and discussion of our manuscript.

2. The existence of fixed points

In this section, we are about to study the existence of possible fixed points, by solving the nonlinear system given by

$\begin{align} \left\{ \begin{array}{l} x\Rightarrow x+\frac{rx}{1+ky}(1-\frac{x}{K})-c(x-R)y-q_1E_1x, \\ y\Rightarrow y-dy+e(x-R)y-q_2E_2y.\\ \end{array} \right. \end{align}$ $

(2.1)

Theorem 2.1. Fixed points obtained are as follows:

(1) $ A_1(0, 0) $ always exists.

(2) $ A_2(W, 0) $ is feasible if $ r > q_1E_1 $, where $ W = \frac{K(r-q_1E_1)}{r} $.

(3) $ A_3(Q, P) $ exists if $ S < 0 $ is satisfied, where $ Q = \frac{d+q_2E_2+eR}{e} $, $ P = \frac{-V+\sqrt{V^2-4ZS}}{Z} $ and $ Z = ckQ-ckR $, $ V = cQ-cR+q_1E_1kQ $, $ S = q_1E_1Q-rQ(1-\frac{Q}{K}) $.

Proof. Obviously, $ A_1(0, 0) $ always holds true. For $ A_2(W, 0) $, we just need to consider that $ W $ is positive (i.e., $ r > q_1E_1 $). For interior fixed point $ A_3 $, through simple calculation for system (2.1), we obtain Eq (2.2)

$\begin{align} \left\{ \begin{array}{l} Q = Q+\frac{rQ}{1+ky}(1-\frac{Q}{K})-c(Q-R)y-q_1E_1Q, \\ Q = \frac{d+q_2E_2+eR}{e}.\\ \end{array} \right. \end{align}$ $

(2.2)

Appropriate transformation of the first term in Eq (2.2) yields

$\begin{align} (ckQ-ckR)y^2+(cQ-cR+q_1E_1kQ)y+q_1E_1Q-rQ(1-\frac{Q}{K}) = 0. \end{align}$ $

(2.3)

Let's define

$ F(y) = (ckQ-ckR)y^2+(cQ-cR+q_1E_1kQ)y+q_1E_1Q-rQ(1-\frac{Q}{K}) = Zy^2+Vy+S. $

And $ Z = ckQ-ckR $, $ V = cQ-cR+q_1E_1kQ $, and $ S = q_1E_1Q-rQ(1-\frac{Q}{K}) $, notice that

(ⅰ) $ V^2-4ZS > 0 $ implies that there exist two roots for $ F(y) = 0 $;

(ⅱ) Since $ Q = \frac{d+q_2E_2+eR}{e} > R $, $ V = c(Q-R)+q_1E_1kQ > 0 $ and $ Z = ck(Q-R) > 0 $ are valid, it indicates the axis of symmetry $ \overline{y} = -\frac{V}{2Z} < 0 $.

Now we only need to judge the constant term of $ F(y) $ to know whether $ F(y) $ has positive roots. When $ S > 0 $, there are two negative roots, whereas when $ S < 0 $, there is a unique positive root

$ P = \frac{-V+\sqrt{V^2-4ZS}}{Z} > 0. $

At this point, we have completed the proof of Theorem 2.1.

3. Stability analysis

We investigate the local stability of the above fixed points of system (1.3). The Jacobian matrix of (2.1) at $ (x, y) $ is given by

$\begin{align} J(x, y) = \left( \begin{array}{ccc} j_1 & j_2\\ j_3 & j_4\\ \end{array} \right), \end{align}$ $

(3.1)

where

$\begin{align} \begin{array}{l} j_1 = 1+\frac{r}{(1+ky)}-\frac{2rx}{K(1+ky)}-cy-q_1E_1, \\ j_2 = -\frac{krx}{(1+ky)^2}(1-\frac{x}{K})-c(x-R), \\ j_3 = ey, \ \ j_4 = 1-d+e(x-R)-q_2E_2. \end{array} \end{align}$ $

(3.2)

The characteristic equation of this Jacobian matrix (3.1) is

$\begin{align} F(\lambda) = \lambda^2-Tr(J)\lambda+Det(J) = 0, \end{align}$ $

(3.3)

where

$ Tr(J) = j_1+j_4, \ \ Det(J) = j_1j_4-j_2j_3. $

So as to the properties of fixed points of system (1.3), we provide the following lemma.

Lemma 3.1. Assume that $ F(1) > 0 $. Then

(ⅰ) $ |\lambda_1| < 1 $ and $ |\lambda_2| < 1 $ if and only if $ F(-1) > 0 $ and $ Det(J) < 1 $;

(ⅱ) $ |\lambda_1| > 1 $ and $ |\lambda_2| > 1 $ if and only if $ F(-1) > 0 $ and $ Det(J) > 1 $;

(ⅲ) $ (|\lambda_1| > 1 $ and $ |\lambda_2| < 1) $ or $ (|\lambda_1| < 1 $ and $ |\lambda_2| > 1) $ if and only if $ F(-1) < 0 $;

(ⅳ) $ |\lambda_1| = 1 $ and $ |\lambda_2|\neq1 $ if and only if $ F(-1) = 0 $ and $ Tr(J)\neq0, 2 $;

(ⅴ) $ \lambda_1 $ and $ \lambda_2 $ are conjugate complex roots, and $ |\lambda_1| = |\lambda_2| = 1 $ if and only if $ Tr(J)^2-4Det(J) < 0 $ and $ Det(J) = 1 $.

Proposition 3.1. The conclusions of fixed point $ A_1(0, 0) $ are as follows:

(ⅰ) $ A_1(0, 0) $ is sink if $ -2 < r-q_1E_1 < 0 $ and $ 0 < d+eR+q_2E_2 < 2 $;

(ⅱ) $ A_1(0, 0) $ is source (unstable) if $ r-q_1E_1 < -2 $ and $ d+eR+q_2E_2 > 2 $;

(ⅲ) $ A_1(0, 0) $ is saddle if $ r-q_1E_1 < -2 $, $ d+eR+q_2E_2 < 2 $ or $ -2 < r-q_1E_1 < 0 $ and $ d+eR+q_2E_2 > 2 $;

(ⅳ) $ A_1(0, 0) $ is non–hyperbolic if $ r-q_1E_1 = -2 $ and $ d+eR+q_2E_2\neq0\ {\rm, }\ 2 $.

Proof. At fixed point $ A_1(0, 0) $, the Jacobian matrix is

$\begin{align} J_1 = \left( \begin{array}{ccc} 1+r-q_1E_1 & cR\\ 0& 1-d-eR-q_2E_2\\ \end{array} \right). \end{align}$ $

(3.4)

From the above matrix, we can calculate $ Tr(J_1) $ and $ Det(J_1) $ as follows:

$\begin{array}{l} Tr(J_1) = 2+r-q_1E_1-q_2E_2-d-eR, \end{array}$ $

$\begin{array}{l} Det(J_1) = 1-d-eR-q_2E_2+r-rd-erR-rq_2E_2-q_1E_1+dq_1E_1+eRq_1E_1+q_1E_1q_2E_2. \end{array}$ $

At the same time, we can derive the eigenvalues

$ \lambda_{1J_1} = 1+r-q_1E_1, \ \ \lambda_{2J_1} = 1-d-eR-q_2E_2. $

By Lemma 3.1, $ A_1(0, 0) $ is sink if $ \mid \lambda_{1J_1, 2J2} \mid < $1 ($ -2 < r-q_1E_1 < 0 $ and $ 0 < d+eR+q_2E_2 < 2 $). Similarly, $ A_1(0, 0) $ is source if $ r-q_1E_1 < -2 $ and $ d+eR+q_2E_2 > 2 $, saddle if $ r-q_1E_1 < -2 $, $ d+eR+q_2E_2 < 2 $ or $ -2 < r-q_1E_1 < 0 $ and $ d+eR+q_2E_2 > 2 $ and non–hyperbolic if $ r-q_1E_1 = -2 $ and $ d+eR+q_2E_2\neq0{\rm, }2 $.

Proposition 3.2. The conclusions of fixed point $ A_2(W, 0) $ are as follows:

(ⅰ) $ A_2(W, 0) $ is sink if $ -2 < q_1E_1-r < 0 $ and $ -2+d+eR+q_2E_2 < \frac{eK(r-q_1E_1)}{r} < d+eR+q_2E_2 $;

(ⅱ) $ A_2(W, 0) $ is source (unstable) if $ q_1E_1-r < -2 $ and $ \frac{eK(r-q_1E_1)}{r} < -2+d+eR+q_2E_2 $ (or $ \frac{eK(r-q_1E_1)}{r} > d+eR+q_2E_2 $);

(ⅲ) $ A_2(W, 0) $ is saddle if $ q_1E_1-r < -2 $ and $ -2+d+eR+q_2E_2 < \frac{eK(r-q_1E_1)}{r} < d+eR+q_2E_2 $ or $ -2 < q_1E_1-r < 0 $ and $ \frac{eK(r-q_1E_1)}{r} > d+eR+q_2E_2 $ or $ \frac{eK(r-q_1E_1)}{r} < -2+d+eR+q_2E_2 $;

(ⅳ) $ A_2(W, 0) $ is non–hyperbolic if $ q_1E_1-r = -2 $ and $ \frac{eK(r-q_1E_1)}{r}\neq -2+d+eR+q_2E_2 $ and $ \frac{eK(r-q_1E_1)}{r}\neq d+eR+q_2E_2 $.

Proof. Bring $ A_2(\frac{K(r-q_1E_1)}{r}, 0) $ into (3.1), then we obtain

$\begin{align} J_2 = \left( \begin{array}{ccc} 1-r+q_1E_1 & -\frac{kK(r-q_1E_1)}{(1+ky)^2}(1-\frac{r-q_1E_1}{r})-c(\frac{K(r-q_1E_1)}{r}-R)\\ \\ 0& 1-d+\frac{eK(r-q_1E_1)}{r}-eR-q_2E_2 \end{array} \right). \end{align}$ $

(3.5)

At the same time, we can derive the eigenvalues $ \lambda_{1J_2} $ and $ \lambda_{2J_2} $ from the above matrix

$\begin{array}{l} \lambda_{1J_2} = 1-r+q_1E_1, \\ \lambda_{2J_2} = 1+\frac{eK(r-q_1E_1)}{r}-d-eR-q_2E_2. \end{array}$ $

By Lemma 3.1, if $ \lambda_{1J_2} $ = $ \mid1-r+q_1E_1\mid < 1 $ and $ \lambda_{2J_2} $ = $ \mid1+\frac{eK(r-q_1E_1)}{r}-d-eR-q_2E_2\mid < 1 $, i.e., $ -2 < q_1E_1-r < 0 $ and $ -2+d+eR+q_2E_2 < \frac{eK(r-q_1E_1)}{r} < d+eR+q_2E_2 $, then $ A_2(W, 0) $ is sink. Analogously, $ A_2(W, 0) $ is source if $ q_1E_1-r < -2 $ and $ \frac{eK(r-q_1E_1)}{r} < -2+d+eR+q_2E_2 $ (or $ \frac{eK(r-q_1E_1)}{r} > d+eR+q_2E_2 $), saddle if $ q_1E_1-r < -2 $ and $ -2+d+eR+q_2E_2 < \frac{eK(r-q_1E_1)}{r} < d+eR+q_2E_2 $ or $ -2 < q_1E_1-r < 0 $ and $ \frac{eK(r-q_1E_1)}{r} > d+eR+q_2E_2 $ or $ \frac{eK(r-q_1E_1)}{r} < -2+d+eR+q_2E_2 $ and non–hyperbolic if $ q_1E_1-r = -2 $ and $ \frac{eK(r-q_1E_1)}{r}\neq -2+d+eR+q_2E_2 $ and $ \frac{eK(r-q_1E_1)}{r}\neq d+eR+q_2E_2 $.

Proposition 3.3. The conclusions of fixed point $ A_3(Q, P) $ are as follows:

(ⅰ) $ A_3(Q, P) $ is sink if $ 1+Tr(J_3)+Det(J_3) > 0 $ and $ Det(J_3) < 1 $;

(ⅱ) $ A_3(Q, P) $ is source (unstable) if $ 1+Tr(J_3)+Det(J_3) > 0 $ and $ Det(J_3) > 1 $;

(ⅲ) $ A_3(Q, P) $ is saddle if $ 1+Tr(J_3)+Det(J_3) < 0 $;

(ⅳ) $ A_3(Q, P) $ is non–hyperbolic if $ 1+Tr(J_3)+Det(J_3) = 0 $ and $ Tr(J_3)\neq0, 2 $ (or $ Det(J_3) = 1 $ and $ |Tr(J_3)| < 2 $).

Proof. We put $ A_3(\frac{d+q_2E_2+eR}{e}, \frac{-(cx-cR+q_1E_1kx)+\sqrt{(cx-cR+q_1E_1kx)^2-4ck(x-R)[q_1E_1x-rx(1-\frac{x}{R})]}}{2ck(x-R)}) $ into (3.1), and if $ q_1E_1Q-rQ(1-\frac{Q}{K}) < 0 $ and $ 4ac < 0 $ are satisfied, then it yields that

$\begin{align} J_3 = \left( \begin{array}{ccc} 1+\frac{r}{(1+kP)^2}-\frac{2rQ}{K(1+kP)}-cP-q_1E_1 & -\frac{krQ}{(1+kP)^2}(1-\frac{Q}{K})-c(Q-R)) \\ \\ eP& 1-d+e(Q-R)-q_2E_2\\ \end{array} \right). \end{align}$ $

(3.6)

From the above matrix, we can calculate its $ Tr(J_3) $ and $ Det(J_3) $:

$\begin{array}{l} Tr(J_3) = 2+\frac{r}{1+kP}-\frac{2rQ}{K(1+kP)}-cP-q_1E_1-d-q_2E_2+e(Q-R), \\ \\ Det(J_3) = 1+[\frac{r(2Q-K)}{K(1+kP)}+q_1E_1-1][d+q_2E_2-e(Q-R)]-(cP+q_1E_1)+\frac{r(K-2Q)}{K(1+kP)}\\ \qquad\qquad\ +cP(d+q_2E_2)+\frac{kreQP}{(1+kP)^2}(1-\frac{Q}{K}).\\ \end{array}$ $

The characteristic equation can be written as follows:

$\begin{align} F(\lambda) = \lambda^2-Tr(J_3)\lambda+Det(J_3) = 0. \end{align}$ $

(3.7)

By Lemma 3.1, $ A_3(Q, P) $ is sink if $ F(-1) = 1+Tr(J_3)+Det(J_3) > 0 $ and $ Det(J_3) < 0 $. In the same way, $ A_3(Q, P) $ is source if $ 1+Tr(J_3)+Det(J_3) > 0 $ and $ Det(J_3) > 1 $, saddle if $ 1+Tr(J_3)+Det(J_3) < 0 $ and non–hyperbolic if $ 1+Tr(J_3)+Det(J_3) = 0 $ and $ Tr(J_3)\neq0, 2 $ (or $ Det(J_3) = 1 $ and $ |Tr(J_3)| < 2 $).

4. Flip bifurcation

It is well known that flip bifurcation may occur when $ A_3(Q, P) $ is non–hyperbolic. If $ (r, k, K, c, R, q_1, E_1, d, q_2, E_2, e)\in FB $, where $ FB = \bigg\{r, k, K, c, R, q_1, E_1, d, q_2, E_2, e > 0; 4+\frac{kreQP}{(1+kP)^2}(1-\frac{Q}{K})+\frac{2r(K-Q)}{K(1+kP)}+cP(d+q_2E_2)$ = $[2-q_1E_1-\frac{r(2Q-K)}{K(1+kP)}][d+q_2E_2-e(Q-R)]; \frac{r(K-2Q)}{K(1+kP)}\neq(cP+q_1E_1)+[d+q_2E_2-e(Q-R)]$, $ (cP+q_1E_1)+[d+q_2E_2-e(Q-R)]-2\bigg\}. $

In this section we consider the bifurcation case at $ A_3 $. It is not difficult to find that system (1.3) experiences flip bifurcation at $ A_3 $ if $ e $ changes in a small range of $ e = \widehat{e} $. Giving $ \overline{e}^* $ $ (where \ \overline{e}^*\ll 1) $ of the parameter $ e $ in a small range of $ e = \widehat{e} $ to system (1.3), it yields that

$\begin{align} \left\{ \begin{array}{l} x\Rightarrow x+\frac{rx}{1+ky}(1-\frac{x}{K})-c(x-R)y-q_1E_1x, \\ y\Rightarrow y-dy+(\widehat{e}+\overline{e}^*)(x-R)y-q_2E_2y.\\ \end{array} \right. \end{align}$ $

(4.1)

In order to translate $ A_3 $ to the origin, we take the transformation $ u = x-Q $ and $ v = y-P $, then the system (4.1) transforms into the following form:

$\begin{align} \left\{ \begin{array}{l} u\Rightarrow u+\frac{r(u+Q)}{1+k(v+P)}(1-\frac{(u+Q)}{K})-c[(u+Q)-R](v+P)-q_1E_1(u+Q), \\ v\Rightarrow v-d(v+P)+(\widehat{e}+\overline{e}^*)[(u+Q)-R](v+P)-q_2E_2(v+P).\\ \end{array} \right. \end{align}$ $

(4.2)

Taylor expansion of system (4.2) at $ (u, v, \overline{e}^*) = (0, 0, 0) $

$\begin{align} \left\{ \begin{array}{l} u\Rightarrow a_{11}u+a_{12}v+a_{13}u^2+a_{14}uv+o(u, v)^3, \\ v\Rightarrow a_{21}u+a_{22}v+a_{23}u^2+a_{24}uv+a_{25}\overline{e}^*+a_{26}v\overline{e}^*+a_{27}u\overline{e}^*.\\ \end{array} \right. \end{align}$ $

(4.3)

Further, we can obtain

$\begin{array}{l} a_{11} = 1+\frac{r}{1+ky}-\frac{2rx}{K(1+ky)}-cy-q_1E_1, a_{12} = \frac{krx}{(1+ky)^2}(1-\frac{x}{K}), \\ a_{13} = -\frac{r}{K(1+ky)}, a_{14} = -\frac{kr}{2(1+ky)^2}+\frac{rxkK}{[K(1+ky)]^2}-\frac{c}{2}, \\ a_{21} = ey, a_{22} = 1-d+e(x-R)-q_2E_2, a_{23} = 0, a_{24} = \frac{e}{2}, \\ a_{25} = (x-R)y, a_{26} = \frac{(x-R)}{2}, a_{27} = -\frac{Ry}{2}. \end{array}$ $

Subsequently, give an invertible matrix $ T $ as follows:

$ T = \left(

$\begin{array}{ccc} a_{12} &a_{12}\\ -1-a_{12} & -1-a_{12}\\ \end{array}$ \right). $

Then using the following conversion

$\begin{align} \begin{array}{l} \begin{split} \left( \begin{array}{ccc} u\\ v\\ \end{array} \right) = & T \left( \begin{array}{ccc} x\\ y\\ \end{array} \right), \end{split} \end{array} \end{align}$ $

(4.4)

system (4.3) becomes

$ \left(

$\begin{array}{ccc} x\\ y\\ \end{array}$ \right) = \left(

$\begin{array}{ccc} -1 & 0 \\ 0 & \lambda_2 \\ \end{array}$ \right) \left(

$\begin{array}{ccc} u\\ v\\ \end{array}$ \right)+ \left(

$\begin{array}{ccc} f(u, v, \overline{e}^*)\\ g(u, v, \overline{e}^*)\\ \end{array}$ \right), $

where

$\begin{array}{l} f(u, v, \overline{e}^*) = \frac{1}{Det(T)}\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]u^2+[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]uv-a_{12}a_{25}\overline{e}^*-a_{12}a_{26}v\overline{e}^*\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -a_{12}a_{27}u\overline{e}^*+o(u, v, \overline{e}^*)^3\}, \end{array}$ $

$\begin{array}{l} g(u, v, \overline{e}^*) = \frac{1}{Det(T)}\{[(1+a_{11})a_{13}-a_{12}a_{13}]u^2+[(1+a_{11})a_{13}+a_{12}a_{14}]uv+a_{12}a_{25}\overline{e}^*+a_{12}a_{26}v\overline{e}^* \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +a_{12}a_{27}u\overline{e}^*+o(u, v, \overline{e}^*)^3\}. \end{array}$ $

According to the theoretical knowledge related to the center manifold theorem ^[40,41], in a small range $ \overline{e}^* = 0 $, it exists a center manifold $ W^c(0, 0) $ at the fixed point $ (0, 0) $ of system $ (1.3) $ of the form

$\begin{align} W^c(0, 0) = \bigg\{(x, y):y = a_1 \overline{e}^*+a_2x^2+a_3x \overline{e}^*+a_4 \overline{e}^{*2}+o((|x|, |\overline{e}^*|)^3)\bigg\}. \end{align}$ $

(4.5)

With some simple transformations, we obtain

$\begin{align} \left\{ \begin{array}{l} u = a_{12}(x+y), \\ v = -(1+a_{11})x+(\lambda_2-a_{11})y, \end{array} \right. \end{align}$ $

(4.6)

with $ uv = -a_{12}(1+a_{11})x^2+a_{12}(\lambda_2-2a_{11}-1)xy+a_{12}(\lambda_2-a_{11})y^2 $ and $ u^2 = a_{12}^2(x^2+2xy+y^2) $. By combining the above (4.5) and (4.6), we can get the corresponding coefficients of (4.5)

$\begin{array}{l} a_1 = \frac{a_{25}}{1-\lambda_2^2}, \\ a_2 = \frac{1}{(1-\lambda_2^2}\bigg\{[(1+a_{11})a_{13}+a_{12}a_{23}]a_{12}-(1+a_{11})[(\lambda_2-a_{11})a_{14}-a_{12}a_{24}]\bigg\}, \\ a_3 = \frac{2a_2a_{25}+(1+a_{11})a_{26}-a_{12}a_{27}}{(1+\lambda_2)^2}, \\ a_4 = \frac{a_3a_{25}}{(1-\lambda_2)^2}-\frac{a_2a_{25}^2}{(1-\lambda_2)(1+\lambda_2)^2}. \end{array}$ $

Limiting the center manifold $ W^c(0, 0) $ to the map $ G^* $

$\begin{align} \begin{array}{l} G^*(x, \overline{e}^* ) = -x+f(x, y, \overline{e}^*) = -x+h_0 \overline{e}^*+ h_1 x^2+h_2x \overline{e}^*+h_3 \overline{e}^{*2}+h_4 x^2\overline{e}^*+h_5 x \overline{e}^{*2}\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +h_6x^3+h_7\overline{e}^{*3}+O(|x|+|\overline{e}^*|)^3. \end{array} \end{align}$ $

(4.7)

The coefficients of the above system are as follows:

$\begin{array}{l} h_0 = -\frac{a_{25}}{(1+\lambda_2)}, \\ h_1 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^2-[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]a_{12}(1+a_{11})\bigg\}, \\ h_2 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{2[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^2a_1-[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]a_{12}(\lambda_2-2a_{11}-1)a_1+\\ \ \ \ \ \ \ \ \ \ a_{12}a_{26}(1+a_{11})-a_{12}^2a_{27}\bigg\}, \\ h_3 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^2a_1^2+[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]a_{12}(\lambda_2-a_{11})a_1^2\\ \ \ \ \ \ \ \ \ \ -a_{12}a_{26}(\lambda_2-a_{11})a_1-a_{12}^2a_{27}a_1\bigg\}, \\ h_4 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^22(a_3+a_1a_2)+[(\lambda_2-a_{11})a_{14}-a_{12}a_{24}]\\ \ \ \ \ \ \ \ \ \ [a_{12}(\lambda_2-2a_{11}-1)a_3+a_{12}(\lambda_2-a_{11})2a_1a_2]\bigg\}, \\ h_5 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^22a_1a_3+[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]a_{12}(\lambda_2-a_{11})2a_1a_3\bigg\}, \\ h_6 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^22a_2+[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]a_{12}(\lambda_2-2a_{11}-1)2a_2\bigg\}, \\ h_7 = \frac{1}{a_{12}(1+\lambda_2)}\bigg\{[(\lambda_2-a_{11})a_{13}-a_{12}a_{13}]a_{12}^22a_2+[(\lambda_2-a_{11})a_{14}-a_{12}a_{14}]a_{12}(\lambda_2-2a_{11}-1)2a_2\\ \ \ \ \ \ \ \ \ \ -a_{12}a_{26}(\lambda_2-a_{11})a_4-a_{12}^2a_{27}a_4\bigg\}. \end{array}$ $

In order to make system (1.3) generate flip bifurcation, it is sufficient to ensure that both discriminatory quantities $ \Upsilon_1 $ and $ \Upsilon_2 $ are not equal to $ 0 $, i.e.,

$\begin{align} \left\{ \begin{array}{l} \Upsilon_1 = \bigg(2\frac{\partial^2 G^*}{\partial x\partial\overline{e}^*}+\frac{\partial G^*}{\partial\overline{e}^*}\times \frac{\partial^2 G^*}{\partial x^2}\bigg) \bigg|_{(0, 0)} = 2(h_2+h_0h_1)\neq 0, \\ \Upsilon_2 = \bigg(\frac{1}{2} (\frac{\partial^2 G^*}{\partial x^2})^2+\frac{1}{3} (\frac{\partial^3 G^*}{\partial x^3})\bigg) \bigg|_{(0, 0)} = 2(h_6+h_1^2)\neq 0.\\ \end{array} \right. \end{align}$ $

(4.8)

Theorem 4.1. System $ (1.3) $ exists flip bifurcation at $ A_3 $ when $ e $ changes in a small range of $ \widehat{e} $ and $ \Upsilon_1 \neq 0 $ and $ \Upsilon_2 \neq 0 $ are satisfied. Moreover, if $ \Upsilon_2 > 0 $ (or $ \Upsilon_2 < 0 $), then the period–2 point of the bifurcation from $ A_3 $ is stable (or unstable).

By a similar method as above, $ k $ can be selected as the bifurcation parameter, and we can limit the system (1.3) to

$\begin{align} \begin{array}{l} \sigma(x) = -x+\epsilon_1k^*+\epsilon_2x^2+\epsilon_3k^{*2}+\epsilon_4x^2k^*+\epsilon_5xk^{*2}+\epsilon_6x^3+\epsilon_7k^{*2}+o\bigg\{(|x|, |k^*|)^4\bigg\}. \end{array} \end{align}$ $

(4.9)

Similarly, to ensure that flip bifurcation occurs in system (1.3), we set the two discriminants $ \Psi_1 $ and $ \Psi_2 $ to be non-zero

$\begin{align} \left\{ \begin{array}{l} \Psi_1 = \bigg(2\frac{\partial^2 \sigma}{\partial x\partial k^*}+\frac{\partial \sigma}{\partial k^*}\times \frac{\partial^2 \sigma}{\partial x^2}\bigg) \bigg|_{(0, 0)} = 2(\epsilon_2+\epsilon_0h_1)\neq 0, \\ \Psi_1 = \bigg(\frac{1}{2} (\frac{\partial^2 \sigma}{\partial x^2})^2+\frac{1}{3} (\frac{\partial^3 \sigma}{\partial x^3})\bigg) \bigg|_{(0, 0)} = 2(\epsilon_6+\epsilon_1^2)\neq 0.\\ \end{array} \right. \end{align}$ $

(4.10)

The computational processes of $ \Psi_1 $ and $ \Psi_1 $ are shown in Appendix A.

5. Neimark–Sacker bifurcation

Here we focus on the Neimark–Sacker bifurcation of system (1.3) at $ A_3 $. System (1.3) exists Neimark–Sacker bifurcation at $ A_3 $ if and only if all parameters belong to the set $ NSB = \bigg\{r, k, K, c, R, q_1, E_1, d, q_2, E_2$, $e > 0\ | \ e = \frac{\big(\frac{r(2Q-K)}{K(1+kP)}+q_1E_1-1\big)(d+q_1E_1)-(cP+q_1E_1)+\frac{r(K-2Q)}{K(1+kP)}+cP(d+q_2E_2)}{(Q-R)\big(\frac{r(2Q-K)}{K(1+kP)}+q_1E_1-1\big)-\frac{krQP}{(1+kP)^2}(1-\frac{Q}{K})}$, $ |2+\frac{r}{1+kP}-\frac{2rQ}{K(1+kP)}-cP-q_1E_1-d-q_2E_2+e(Q-R)| < 2$; $q_1E_1-r(1-\frac{Q}{k}) < 0 \mbox{ and } r, k, K, c, R, q_1, E_1, d, q_2, E_2 > 0\bigg\} $.

The characteristic equation (3.3) associated with (3.1) at $ A_3 $ is given by

$ F(\lambda) = \lambda^2-Tr(J_3)\lambda+Det(J_3) = 0. $

There exists a pair of conjugate complex roots $ \lambda_1, \lambda_2 $ of the above characteristic equation at $ A_3 $, i.e.,

$ \lambda_{1, 2} = \frac{Tr(J_3)\pm i\sqrt{4Det(J_3)-(Tr(J_3))^2}}{2}. $

Given a parameter $ e^* $ (where $ e^*\ll 1 $) and $ e $ is located in a range of $ e = \overline{e} $ in (3.1), we have

$\begin{align} \begin{array}{l} x\Rightarrow x+\frac{rx}{1+ky}(1-\frac{x}{K})-c(x-R)y-q_1E_1x, \\ y\Rightarrow y-dy+(\widehat{e}+e^*)(x-R)y-q_2E_2y.\\ \end{array} \end{align}$ $

(5.1)

Moving $ A_3 $ to the original point, we take the transformations $ \gamma = x-Q $ and $ \delta = y-P $, system (5.1) transforms into the following form:

$\begin{align} \begin{array}{l} \gamma \Rightarrow \gamma+\frac{r(\gamma+Q)}{1+k(\delta+P)}(1-\frac{(\gamma+Q)}{K})-c[(\gamma+Q)-R](\delta+P)-q_1E_1(\gamma+Q), \\ \delta\Rightarrow \delta-d(\delta+P)+(\widehat{e}+e^*)[(\gamma+Q)-R](\delta+P)-q_2E_2(\delta+P).\\ \end{array} \end{align}$ $

(5.2)

Taylor expansion of the above system (5.2)

$\begin{align} \begin{array}{l} \gamma\Rightarrow c_{11}\gamma+c_{12}\delta+c_{13}\gamma^2+c_{14}\gamma\delta+o(\gamma, \delta)^3, \\ \delta\Rightarrow c_{21}\gamma+c_{22}\delta+c_{23}\gamma^2+c_{24}\gamma\delta+c_{25}e^*+c_{26}\delta e^*+c_{27}\gamma e^*+o(\gamma, \delta)^3.\\ \end{array} \end{align}$ $

(5.3)

Parameters of (5.3) are as follows:

$\begin{array}{l} c_{11} = 1+\frac{r}{1+ky}-\frac{2rx}{K(1+ky)}-cy-q_1E_1, c_{12} = \frac{krx}{(1+ky)^2}(1-\frac{x}{K}), \\ c_{13} = -\frac{r}{K(1+ky)}, c_{14} = -\frac{kr}{2(1+ky)^2}+\frac{rxk}{K(1+ky)^2}-\frac{c}{2}, \\ c_{21} = ey, c_{22} = 1-d+e(x-R)-q_2E_2, c_{23} = 0, c_{24} = \frac{e}{2}, \\ c_{25} = (x-R)y, c_{26} = \frac{(x-R)}{2}, c_{27} = -\frac{Ry}{2}. \end{array}$ $

When $ (\gamma, \delta) = (0, 0) $, the characteristic roots of (5.3) above are as follows:

$ \lambda_{1, 2} = \frac{TrJ(e^*)\pm i\sqrt{4DetJ(e^*)-(TrJ(e^*))^2}}{2}. $

So we have $ |\lambda_{1, 2}(e^*)| = [DetJ(e^*)]^\frac{1}{2} $. When $ e^* = 0 $, then

$\begin{align} \ \theta_1 = \frac{d|\lambda_{1, 2}|}{de}\bigg|_{e = e^*}\neq0. \end{align}$ $

(5.4)

Set $ \alpha_1 = Re(\lambda_{1, 2}) $ and $ \beta_1 = Im(\lambda_{1, 2}) $. We obtain an invertible matrix $ N $

$ N = \left(

$\begin{array}{ccc} c_{12} &0\\ \alpha_1-c_{11}& \beta_1\\ \end{array}$ \right). $

Then use the following conversion:

$\begin{align} \begin{array}{l} \begin{split} \left( \begin{array}{ccc} \gamma\\ \delta\\ \end{array} \right) = & N \left( \begin{array}{ccc} x\\ y\\ \end{array} \right). \end{split} \end{array} \end{align}$ $

(5.5)

System (5.3) becomes

$ \left(

$\begin{array}{ccc} x\\ y\\ \end{array}$ \right) = \left(

$\begin{array}{ccc} \alpha_1 & -\beta_1 \\ \beta & \alpha_1 \\ \end{array}$ \right) \left(

$\begin{array}{ccc} \gamma\\ \delta\\ \end{array}$ \right)+ \left(

$\begin{array}{ccc} F(\gamma, \delta, \overline{e}^*)\\ G(\gamma, \delta, \overline{e}^*)\\ \end{array}$ \right), $

where

$\begin{array}{l} F(\gamma, \delta, e^*) = -\frac{1}{c_{12}}[c_{12}\gamma^2+c_{14}\gamma\delta+O(\gamma, \delta)^3], \\ G(\gamma, \delta, e^*) = -\frac{1}{c_{12}\beta_1}\bigg\{(c_{11}-\alpha_1)c_{13}\gamma^2+(c_{11}-\alpha_1)c_{14}\gamma\delta+c_{12}c_{23}\gamma^2+c_{12}c_{24}\gamma\delta+c_{12}c_{25}e^*\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +c_{12}c_{26}\delta e^*+c_{12}c_{27}\gamma e^*+O(\gamma, \delta)^3\bigg\}. \end{array}$ $

We can know from the transformation (5.5) that $ \gamma = c_{12}x $ and $ \delta = (\alpha_1-c_{11})x-\beta_1y $, therefore

$\begin{array}{l} F(x, y, e^*) = -\frac{1}{c_{12}}\bigg\{c_{12}(c_{12}x)^2+c_{14}c_{12}x[(\alpha_1-c_{11})x-\beta_1y]+O(\gamma, \delta)^3\bigg\}, \\ G(x, y, e^*) = -\frac{1}{c_{12}\beta_1}\bigg\{[(c_{11}-\alpha_1)c_{13}+c_{12}c_{23}](c_{12}x)^2+[(c_{11}-\alpha_1)c_{14}+c_{12}c_{24}](c_{12}x)[(\alpha_1-c_{11})x-\beta_1y]\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +c_{12}c_{25}e^*+c_{12}c_{26}[(\alpha_1-c_{11})x-\beta_1y]e^*+c_{12}c_{27}(c_{12}x)e^*+O(\gamma, \delta)^3\bigg\}. \end{array}$ $

The Neimark–Sacker bifurcation occurs on system (1.3) if the following condition is met:

$\begin{align} l_1 = -Re\bigg\{\frac{(1-2\overline{\lambda})\overline{\lambda}^2}{1-\lambda}\xi_{11}\xi_{20}\bigg\}-\frac{1}{2}|\xi_{11}|^2-|\xi_{02}|^2+Re(\overline{\lambda}\xi_{21})\neq0, \end{align}$ $

(5.6)

where

$\begin{array}{l} \xi_{11} = \frac{1}{4}[F_{xx}+F_{yy}+i(G_{xx}+G_{yy})], \\ \xi_{20} = \frac{1}{8}[F_{xx}+F_{yy}+2G_{xy}+i(G_{xx}-G_{YY}-2F_{xy})], \\ \xi_{02} = \frac{1}{8}[F_{xx}+F_{yy}-2G_{xy}+i(G_{xx}-G_{yy}+2F_{xy})], \\ \xi_{21} = \frac{1}{16}[F_{xxx}+F_{xyy}+G_{yyy}+G_{xxy}+i(G_{xxx}+G_{xyy}-F_{xxy}-F_{yyy})]. \end{array}$ $

By calculation we have

$\begin{array}{l} F_{xx} = 2[c_{12}c_{13}+c_{13}(\alpha_1-c_{11})], \ F_{yy} = 0, \ F_{xy} = -c_{13}\beta_1, \ F_{xxx} = F_{xyy} = F_{xxy} = F_{yyy} = 0, \\ G_{xx} = -\frac{2}{c_{12}\beta_1}\bigg\{[(c_{11}-\alpha_1)c_{13}+c_{12}c_{23}]c_{12}^2-[(c_{11}-\alpha_1)c_{14}+c_{12}c_{24}]c_{12}(\alpha_1-c_{11})\bigg\}, \ G_{yy} = 0, \\ G_{xy} = (c_{11}-\alpha_1)c_{14}+c_{12}c_{24}, \ G_{xxx} = G_{xyy} = G_{xxy} = G_{yyy} = 0. \end{array}$ $

Theorem 5.1. If $ l_1 $ defined is nonzero, then system (1.3) exists Neimark–Sacker bifurcation at $ A_3 $ provided that $ e $ changes in a small range of $ e = e^* $. Moreover, if $ l_1 < 1 $ (or $ l_1 > 1 $), then an attractively invariant closed curve bifurcates at $ A_3 $.

Similarly, consider $ k $ as a bifurcation parameter, it can be obtained that

$\begin{align} \tau_1 = -Re\bigg\{\frac{(1-2\overline{\lambda})\overline{\lambda}^2}{1-\lambda}\varsigma_{11}\varsigma_{20}\bigg\}-\frac{1}{2}|\varsigma_{11}|^2-|\varsigma_{02}|^2+Re(\overline{\lambda}\varsigma_{21})\neq0, \end{align}$ $

($B.3$)

Neimark–Sacker bifurcation occurs in system (1.3). The detailed calculation process of $ \tau_1 $ is shown in Appendix B.

6. Optimal harvesting policy

Our aim is to maximize net returns and maintain ecological balance. Assign the net income function as $ M = \sum\exp(-\delta t)\{(p_1q_1x-h_1)E_1(t)+(p_2q_2x-h_2)E_2(t)\} $, where $ \delta $ expresses the discount rate. Thus

$ \max\sum\limits_{n = 1}^k\exp(-\delta n)\{(p_1q_1x-h_1)E_1(n)+(p_2q_2x-h_2)E_2(n)\}, $

such that

$\begin{align} \left\{ \begin{array}{l} x_{n+1} = x_{n}+\frac{rx_{n}}{1+ky_{n}}(1-\frac{x_{n}}{K})-c(x_{n}-R)y-q_1E_1x_{n}, \\ y_{n+1} = y_{n}-dy_{n}+(\widehat{e}+\overline{e}^*)(x_{n}-R)y_{n}-q_2E_2y_{n}, \\ x_1 = x_0 , y_1 = y_0 \\ 0\leq E_1(t)\leq max E_1 , 0\leq E_2(t)\leq max E_2. \end{array} \right. \end{align}$ $

(6.1)

The Hamiltonian function could be written in the following form:

$\begin{array}{l} H_n = exp(-\delta n)\{(p_1q_1x_{n}-h_1)E_1(n)+(p_2q_2x_{n}-h_2)E_2(n)\}+\lambda_{1(n+1)}[x_{n}+\frac{rx_{n}}{1+ky_{n}}(1-\frac{x_{n}}{K})\\ \ \ \ \ \ \ \ \ \ -c(x_{n}-R)y_{n}-q_1E_1x_{n}]+\lambda_{2(n+1)}[y_{n}-dy_{n}+(\widehat{e}+\overline{e}^*)(x_{n}-R)y_{n}-q_2E_2y_{n}]. \end{array}$ $

According to the Pontryagin maximum principle, we have

$ \frac{\partial H}{\partial H_i} = 0 \ (i = 1, 2);\ \ \frac{d\lambda_{1(n+1)}}{d t } = -\frac{\partial H}{d x}\ \ \mbox{and}\ \ \frac{d\lambda_{2(n+1)}}{d t } = -\frac{\partial H}{d y}. $

Further we can obtain

$\begin{array}{l} \frac{\partial H}{\partial H_1} = 0 \Rightarrow exp(-\delta t)(p_1q_1x_n-h_1)-q_1x_n\lambda_{1(n+1)} = 0\Rightarrow \lambda_1 = exp(-\delta t)(p_1-\frac{h_1}{q_1x_n}), \\ \frac{\partial H}{\partial H_2} = 0 \Rightarrow exp(-\delta t)(p_2q_2x_n-h_2)-q_2x_n\lambda_{2(n+1)} = 0\Rightarrow \lambda_2 = exp(-\delta t)(p_2-\frac{h_2}{q_2x_n}), \\ \lambda_{1(n+1)} = -\frac{\partial H}{d x} = -exp(\delta t)p_1q_1E_1-\lambda_{1(n+1)}[1+\frac{r}{1+ky}-\frac{2rx}{K(1+ky)}-cy-q_1E_1]-\lambda_{2(n+1)}ey, \\ \lambda_{2(n+1)} = -\frac{\partial H}{d y} = -exp(\delta t)p_2q_2E_2-\lambda_{1(n+1)}[-\frac{krx}{(1+ky)^2}(1-\frac{x}{K})-c(x-R)y]\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -\lambda_{2(n+1)}[1-d+e(x-R)-q_2E_2]. \end{array}$ $

By combining the above equations that

$\begin{array}{l} E_1 = -\frac{x}{h_1}\{(p_1-\frac{h_1}{q_1x})[2+\frac{r}{1+ky}-\frac{2rx}{K(1+ky)}-cy]+(p_2-\frac{h_2}{q_2y})ey\}, \\ E_2 = -\frac{y}{h_1}\{(p_1-\frac{h_1}{q_1x})[-\frac{krx}{(1+ky)^2}(1-\frac{x}{K})-c(x-R)y]+(p_2-\frac{h_2}{q_2y})[1-d+e(x-R)]\}. \end{array}$ $

7. Numerical simulations

Example 7.1. To verify the theoretical results numerically. First, we value all parameters in turn $ r = 2.001 $, $ k = 0.57 $, $ K = 4.078 $, $ c = 0.743 $, $ R = 0.829 $, $ d = 0.8 $, $ q_1 = 0.467 $, $ E_1 = 0.201 $, $ q_2 = 0.42 $, $ E_2 = 0.435 $, and $ e\subset(0, 3.5) $ with initial conditions $ (x_0, y_0) = (1, 1.5) $. By simple calculation, we can know that $ \mid \lambda_{1, 2}\mid = 1 $, Neimark–Sacker bifurcation may appear at (2.29572, 1.65075). We can see from Figure 1 that Neimark–Sacker bifurcation occurs when $ e = 0.67 $. It follows that $ \mid \lambda_1\mid = 1 $ and $ \mid \lambda_2\mid \neq 1 $, flip bifurcation may appear at (1.214, 1.8965). We can see from Figure 1 that flip bifurcation occurs when $ e = 2.55 $.

Figure 1. (a) Bifurcation diagram of system (1.3) in the $ (e, x) $–plane for $ r = 2.001 $, $ k = 0.57 $, $ K = 4.078 $, $ c = 0.743 $, $ R = 0.829 $, $ d = 0.8 $, $ q_1 = 0.467 $, $ E_1 = 0.201 $, $ q_2 = 0.42 $, $ E_2 = 0.435 $, and $ e\subset(0, 3.5) $; (b) Bifurcation diagram of system (1.3) in the $ (e, y) $–plane for $ r = 2.001 $, $ k = 0.57 $, $ K = 4.078 $, $ c = 0.743 $, $ R = 0.829 $, $ d = 0.8 $, $ q_1 = 0.467 $, $ E_1 = 0.201 $, $ q_2 = 0.42 $, $ E_2 = 0.435 $, and $ e\subset(0, 3.5) $.

DownLoad: Full-Size Img PowerPoint

From Figure 2(a) and (d), we can know that Neimark–Sacker bifurcation will appear when $ e = 0.67 $, and when $ e $ belongs to (0, 1) and is greater than 0.67, system (1.3) is stable at the internal fixed point $ (2.29572, 1.65075) $, whereas when the bifurcation parameter $ e $ is less than 0.67, the system will lose stability at the internal fixed point $ (2.29572, 1.65075) $. In view of the ecological perspective, when $ e < 0.67 $, a fluctuation will occur, which means that the number of prey and predator populations will continue to fluctuate periodically around a central value, rather than reaching a fixed point. Such fluctuation may reflect seasonal changes in resources or other cyclical environmental factors in the ecosystem. According to Figure 2(b) and (c), with the change of bifurcation parameters $ e $, flip bifurcation will be generated in system (1.3) when its value reaches 2.55. When the bifurcation parameter $ e $ belongs to (1, 3.5) and is less than 2.55, system $ (1.3) $ is stable; conversely, when $ e $ is greater than 2.55, the system will lose stability. Biologically speaking, when $ e > 2.55 $, the number of prey and predators is no longer maintained at a constant level, whereas it presents a state of periodic fluctuation. Such fluctuation could have an impact on the stability of ecosystems and lead to a large fluctuation in the number of prey and predator populations, which can affect the overall ecosystem's sustainability.

Figure 2. (a) Neimark–Sacker bifurcation diagram of system (1.3) in the $ (e, y) $–plane for $ e\subset(0, 0.67) $; (b) Flip bifurcation diagram of system (1.3) in the $ (e, y) $–plane for $ e\subset(0.67, 3.5) $; (c) Flip bifurcation diagram of system (1.3) in the $ (e, x) $–plane for $ e\subset(0.67, 3.5) $; (d) Neimark–Sacker bifurcation diagram of system (1.3) in the $ (e, x) $–plane for $ e\subset(0, 0.67) $.

DownLoad: Full-Size Img PowerPoint

In the process from Figure 3(a) to (f), it is clearly shown that the stability of the system (1.3) varies with the bifurcation parameter $ e $ when it belongs to the range of $ (1, 3.5) $.

Figure 3. Phase portraits of system (1.3) with different $ e $.

DownLoad: Full-Size Img PowerPoint

Example 7.2. We take the parameters $ r = 3.5, K = 8$, $c = 0.69, k = 3.9, d = 3.6, e = 0.77, q_1 = 0.36, q_2 = 0.21$, $p_1 = 0.75, p_2 = 3, h_1 = 1.41, h_2 = 1.8 $, and assign the initial value $ (x_0, y_0, E_{10}, E_{20}) $ = (1.32, 3, 0.05, 0.37). Through Figure 4, with the increase of the refuge coefficient $ R $, harvesting effort $ E_1 $ for prey will increase rapidly, when the refuge coefficient $ R = 0.03 $, the harvesting of the prey will tend to a stable state. Due to the value of $ E_1 $ for the prey population changing too rapidly, as the refuge coefficient $ R $ increases, the prey population will show a downward trend; until $ R = 0.03 $, the number of the prey population will tend to a stable state. For the predator population, the harvesting $ E_2 $ will also increase rapidly with increasing $ R $. Compared with the prey population, the number of predator populations is inherently smaller. Due to the rapid increase in the refuge coefficient $ R $ and the harvesting $ E_2 $, the predator population will be in a state of extinction. In general, whether it is for predators or prey, the harvesting of both of them should be moderate, so as to better promote the harmonious coexistence between humans and nature and ensure the maximum harvesting.

Figure 4. Optimal harvesting with respect to refuge coefficient $ R $.

DownLoad: Full-Size Img PowerPoint

8. Conclusions and discussion

In this article, we mainly study the bifurcation behavior and optimal harvesting strategy of a discrete predator–prey model in possession of fear and refuge effects. Firstly, the existence and stability of three fixed points of system (1.3) are researched, respectively. At the same time, the flip bifurcation and Neimark–Sacker bifurcation are analysed at the internal fixed point in detail. Together with Example 7.1, when the bifurcation parameter $ e = 0.67 $, Neimark–Sacker bifurcation will occur, and when $ e = 2.55 $, flip bifurcation will appear. It is easy to observe that relying on the increasing conversion rate of $ e $, the unstable system turns stable, and then it becomes unstable again. Besides, the predator populations gradually accumulate for the value of the conversion rate at a high level. Finally, we analyze an optimal harvesting problem and theoretically get the value of the optimal harvesting effort, which implies there exists a value of the harvesting effort guaranteeing maximization of the net revenue. From Example 7.2, with the increase of the refuge coefficient $ R $, $ x^{optimal}, y^{optimal}, E_1^{optimal} $ and $ E_2^{optimal} $ change in the range.

Generally speaking, comparing the continuous–time predator–prey models with the corresponding discrete versions, the conclusions on the existence and stability of equilibrium points (or fixed points) can be studied. Also, bifurcation analysis of the models can be discussed. Relative to bifurcation analysis, several bifurcation phenomena, such as Hopf bifurcation, saddle–node bifurcation, and transcritical bifurcation, are mainly investigated in the continuous system. Otherwise, fold bifurcation, flip bifurcation, and Neimark–Sacker bifurcation can be discussed in the discrete system. All in all, whether discrete or continuous system, by studying these characteristics for the system, we can better understand external disturbances, thereby achieving a better role in maintaining ecological balance.

Author contributions

Jie Liu: Conceptualization, Investigation, Methodology, Validation, Writing-original draft, Formal analysis, Software; Qinglong Wang: Conceptualization, Methodology, Formal analysis, Writing-review and editing, Supervision; Xuyang Cao: Validation, Visualization, Data curation; Ting Yu: Validation, Visualisation, Inspection.

Acknowledgments

The authors thank the editor and referees for their careful reading and valuable comments.

The work is supported by the Natural Science Foundation of Hubei Province (No.2023AFB1095) and the National Natural Science Foundation of China (No.12101211) and the Program for Innovative Research Team of the Higher Education Institution of Hubei Province (No.T201812) and the Teaching Research Project of Education Department of Hubei Province (No.2022367) and the Graduate Education Innovation Project of Hubei Minzu University (Nos.MYK2024071, MYK2023042).

Conflict of interest

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

References

[1]	Zhou H, Smith DW (2001) Advanced technologies in water and wastewater treatment. Can J Civil Eng 28: 49-66. doi: 10.1139/cjce-28-S1-49
[2]	Melin T, Jefferson B, Bixio D, et al. (2006) Membrane bioreactor technology for wastewater treatment and reuse. Desalination 187: 271-282. doi: 10.1016/j.desal.2005.04.086
[3]	Verrecht B, Judd SJ, Guglielmi G, et al. (2008) An aeration energy model for an immersed membrane bioreactor. Water Res 42: 4716-4770.
[4]	Schleper C, Nicol GW (2010) Ammonia-oxidising archaea - physiology, ecology and evolution. Adv Microb Physiol 57: 1-36. doi: 10.1016/B978-0-12-381045-8.00001-1
[5]	Giraldo E, Jjemba P, Liu Y, et al. (2011) Ammonia oxidizing archaea, AOA, population and kinetic changes in a full scale simultaneous nitrogen and phosphorus removal MBR. Proc Water Environ Federation 2011: 3156-3168. doi: 10.2175/193864711802721596
[6]	Huang T, Li N, Huang Y (2012) Modelling of nitrogen removal and control strategy in continuous-flow-intermittent-aeration process. Afr J Biotechnol 11: 10626-10631.
[7]	Guglielmi G, Andreottola G (2011) Alternate anoxic/aerobic operation for nitrogen removal in a membrane bioreactor for municipal wastewater treatment. Water Sci Techn 64: 1730-1735. doi: 10.2166/wst.2011.755
[8]	Curko J, Matosic M, Jakopovic HK, et al. (2010) Nitrogen removal in submerged MBR with intermittent aeration. Desalination Water Treat 24: 7-19. doi: 10.5004/dwt.2010.1118
[9]	Parikh C, Trivedi H, Livingston D (2011) A decade of simultaneous nitrification and denitrification experience in over 60 conventional and MBR applications - lessons learned. Proc Water Environ Federation 2011: 3629-3655. doi: 10.2175/193864711802721901
[10]	Holman JB, Wareham DG (2005) COD, ammonia and dissolved oxygen time profiles in the simultaneous nitrification/denitrification process. Biochem Eng J 22: 125-133. doi: 10.1016/j.bej.2004.09.001
[11]	Fu Z, Yang F, An Y, et al. (2009) Simultaneous nitrification and denitrification coupled with phosphorus removal in a modified anoxic/oxic-membrane bioreactor (A/O-MBR). Biochem Eng J 43: 191-196. doi: 10.1016/j.bej.2008.09.021
[12]	Li YZ, He YL, Ohandja DG, et al. (2008) Simutaneous nitrification-denitrification achieved by an improved internal-loop airlift MBR: comparative study. Bioresource Technol 99: 5867-5872. doi: 10.1016/j.biortech.2007.10.001
[13]	Sarioglu M, Insel G, Artan N, et al. (2009) Model evaluation of simultaneous nitrification and denitrification in a membrane bioreactor operated without an anoxic reactor. J Membrane Sci 337: 17-27. doi: 10.1016/j.memsci.2009.03.015
[14]	Wyffels S, Van Hulle SWH, Boeckx P, et al. (2004) Modelling and simulation of oxygen-limited partial nitrification in a membrane-assisted bioreactor (MBR). Biotechnol Bioeng 86: 531-542. doi: 10.1002/bit.20008
[15]	Daebel H, Manser R, Gujer W (2007) Exploring temporal variations of oxygen saturation constants of nitrifying bacteria. Water Res 41: 1094-1102. doi: 10.1016/j.watres.2006.11.011
[16]	Sharp R, Dailey S, Motyl M, et al. (2012) K_DO experiments at NYC DEP's full-scale demonstration facility. Proc Water Environ Federation 2012: 4176-4182. doi: 10.2175/193864712811708482
[17]	Rieger L, Jones RM, Dold PL, et al. (2014) Ammonia-based feedforward and feeback aeration control in activated sludge process. Water Environ Res 86: 63-73. doi: 10.2175/106143013X13596524516987
[18]	Henze M, Gujer W, Mino T, et al. (2000). Activated sludge models: ASM1, ASM2, ASM2d and ASM3. IWA Publishing, London.
[19]	Tchobanoglous G, Burton FL, Stensel HD (2004) Wastewater engineering, treatment and reuse, 4 Eds., Singapore: McGraw-Hill Education, 690.

This article has been cited by:

1.	Xin Yi, Rong Liu, An age-dependent hybrid system for optimal contraception control of vermin, 2025, 10, 2473-6988, 2619, 10.3934/math.2025122
2.	Tiantian Zhao, Xiaoling Han, Complex dynamics of Hastings–Powell model with additive Allee effect and Crowley–Martin functional response, 2025, 198, 09600779, 116431, 10.1016/j.chaos.2025.116431

Reader Comments

Your name:*

Email:*
© 2015 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)