Coupled Inverse Fluidized Bed Bioreactor with Advanced Oxidation Processes for Treatment of Vinasse

1.   Introduction

Predator-prey dynamics model is one of the most important topics in mathematical ecology research, which is universal and significant. Since the early model proposed by Lotka ^[1] and Volterra ^[2], the studies of predator-prey interactions have attracted rapidly increasing attention. Generally, a predator-prey model can be described by a system of ordinary differential equations as follows:

where $x(t)$ and $y(t)$ are the population size of the prey and predator at time $t$, respectively. $f_1(x(t))$ is the growth rate of prey population; $f_2(y(t))$ is the growth rate of predator population. $c$ is a positive constant ($c < 1$ due to the biological significance) which represents the coefficient in converting prey into a new predator. $p(x(t), y(t))$ represents the functional response, i.e., how predator consumes the prey species. Vast and rich studies ^[3,4,5,6,7] have been done by taking different kinds of functional responses into consideration to investigate the effect of direct predation. However, in the real world, the presence of predators also has indirect effects on the prey population. The outcomes of recent theoretical studies and experimental shreds of evidence ^[8,9,10,11] have unveiled that the indirect impact of predator on prey species sometimes becomes more crucial in altering the dynamics of ecosystems. The threat of predation alone is enough to reduce the growth, survival and fecundity of prey, which is arbitrated through an alternation of prey's behavioural, morphological or psychological trait ^[12,13]. For instance, the sounds of predator will induce continuous fear on birds which compels them to leave their primary nest and move to a safer zone. This kind of behavior leads to serious damage to the birth rate of the new offspring inevitably. Frequent changes in habitat will not only affect the basic reproduction of the individuals but also cause a serious impact on the survival of adults ^[8]. In 2009, Sheriff et al. ^[14] reported that when pregnant snowshoe hares were exposed to a trained dog (predator), the birth of cubs was significantly reduced. In addition, Zanette et al. ^[15] performed an experiment on song sparrows and the outcome exhibited that the quantity of posterity of the song sparrows was decreased by 40% due to fear induced by predators. Shreds of evidence suggest that the indirect impact of predator should not be ignored when analyzing the dynamics of prey-predator interaction. In this context, Wang et al. ^[16] were the first to establish an ordinary differential equation model including the fear effect. Through their work, we know that the level of fear would not affect the dynamic of the system if the predator consumes the prey according to linear functional response, while for the model with Holling type Ⅱ functional response, the level of fear affected predator-prey interactions in several ways. Subsequently, they proposed a predator-prey model with the cost of fear and adaptive avoidance of predators in ^[17]. Inspired by their research, Sasmal ^[18] investigated a model with fear effect and Allee effect where the prey is taken as infected. Sha et al. ^[19] and Hossain ^[20] focused on eco-epidemiological models with fear effect. Few other works have been found in this direction, such as discrete ^[21], delayed ^[22] and diffusive ^[23,24] models.

Recently, Wang and Zou ^[25] investigated a model with the benefits of the anti-predation response and digestion time which is necessary for biomass transfer from prey to predator after predation. The chance of the prey being vulnerable by predators was thought to be decreased for the presence of anti-predation response, which makes the model more in line with the actual situation. Obviously, such anti-predation responses are beneficial to the survival of the prey. In addition, it is worth noting that in reality, predation comes with risks. The prey will resist when it encounters a predator instead of just escaping. For example, leopards run the risk of being stabbed by the hard spines on the back and tail of the porcupine when preying on the porcupine. Once stabbed, the leopard may die due to wound infection or unable to prey. The death of the predator caused by self-defense of the prey can be regarded as the predation of the predator by the prey, which will reduce the growth of the predator population. The reduction is the cost of predation and it may exist whether the predation is successful or not. Noting that this kind of passive predation does not need time to digest and will not affect the next predation, the linear function response is suitable. What needs to be pointed out is that self-defence is different from conservative and protective behaviour and group defence, as it is the fight for survival after the failure of avoidance strategy. {By fighting against the predator, the chance of being predated will obviously be decreased, which will benefit the prey population. Moreover, the prey's self-defence will definitely consume its energy and this will clearly affect the production rate and death rate of the prey. To avoid making the system too complicated, we only consider the cost of the prey's self-defence for the predator. }Hence, we get the following model with fear factor, digestion delay and self-defence:

where $x(t)$ and $y(t)$ have the same meaning as the system (1.1). $\tau \geq 0$ is the average time needed for biomass transfer after predation from prey to predator. All the parameters are positive constants: $r$ represents the intrinsic growth rate of the prey; $d_i \; (i = 1, 2)$ represent the death rate; $a$ and $h$ donate the intra-specific competition of the prey and predator, respectively; $k$ represents the level of fear; $c_i \; (i = 1, 2)$ represent the decreasing rate of reproduction and predation, respectively; $c$ is the biomass transform efficiency constant and $\theta$ represents the level of self-defence. From the realistic point of view, system (1.2) extend model (3.24) in ^[25] by considering the self-defence of prey and intra-specific of predator. Therefore, for the dynamics of the deterministic model, we will focus on both of these parameters in addition to the fear factor.

On the other hand, parameters involved in the system always fluctuate around some average values since the environmental fluctuations, such as temperature, humidity, rainfall and so on. In view of this, many researchers ^{[26,27,28,29]} have introduced stochastic environmental perturbations into deterministic models. Since the death rate $d_i \; (i = 1, 2)$ are most vulnerable by the environmental fluctuations and most influential in determining the dynamics of the system, we will focus on the effects of perturbations on these parameters. In this paper, we perturb $-d_1$ and $-d_2$ by

where $B_1(t)$ and $B_2(t)$ are standard Brownian motion with $B_1(0) = B_2(0) = 0$, $\sigma_1^2$ and $\sigma_2^2$ are the intensity of white noise. Whereupon, we obtain the following stochastic model:

System (1.3) satisfies the initial conditions

where $\phi_1(\varsigma)$ and $\phi_2(\varsigma)$ are both nonnegative continuous functions on $[-\tau, 0]$. In addition, we suppose that $\phi_i(0) > 0 \; (i = 1, 2)$ for biological meaning. Throughout this paper, unless otherwise specified, let $(\Omega, \mathscr{F}, \{\mathscr{F}_t\}_{t\geq0}, P)$ denote a complete probability space with a filtration $\{\mathscr{F}_t\}_{t\geq0}$ satisfying the usual conditions (i.e. it is right continuous and $\mathscr{F}_0$ contains all P-null sets).

The framework of the paper is arranged as follows. In Section 2, we analyze the existence and stability of the equilibrium points and the existence of Hopf bifurcations of system (1.2). In Section 3, we investigate the dynamic behavior of system (1.3), including the existence and uniqueness of the global positive solution and asymptotic property of system (1.3) around the equilibrium points of the corresponding deterministic system (1.2). Several numerical simulations are introduced to support our theoretical results in Section 4. Finally, we give a brief discussion and summarize the main results in Section 5.

2.   Dynamics of the deterministic model

Obviously, system (1.2) has the following three possible non-negative equilibrium points: (ⅰ) The trivial equilibrium $E_0 = (0, 0)$, which always exists; (ⅱ) The predator-free equilibrium $E_1 = (\frac{r-d_1}{a}, 0)$, which exists when $r-d_1 > 0$, matching biological meaning; (ⅲ) The positive (coexistence) equilibrium $E^* = (x^*, y^*)$, where $y^* = \frac{cp}{h(1+c_1k)}\frac{x^*}{1+qx^*}-\frac{\theta}{h}x^*-\frac{d_2}{h}$ and $x^*$ is the positive solution of equation $Ax^5+Bx^4+Cx^3+Dx^2+Ex+F = 0$, where

$A = \frac{ac_2k\theta q^3}{h},$

$B = (2aq+d_1q^2)\frac{c_2kq\theta}{h}-aq^2(q+c_2k\psi_2)+\frac{c_2kq^2\theta^2\psi_1}{h^2}$,

$C = \frac{c_2kq\theta}{h}(a+2d_1q)-(2aq+d_1q^2)(q+c_2k\psi_2)-aq^2(1-\frac{c_2bd_2}{h})-\frac{q^2\theta+2c_2kq\theta}{h}\psi_1 \psi_2 +q^3r$,

$D = \frac{c_2kqd_1\theta}{h}-(q+c_2k\psi_2)(a+2d_1q)+(\frac{q\theta(2c_2kd_2-h)}{h^2}+(q+c_2k\psi_2)\psi_2)\psi_1+\frac{c_2kd_2}{h}+3q^3r-1$,

$E = (1-\frac{c_2kd_2}{h})(\psi_1\psi_2-a-2dq)-(q+c_2k\psi_2)(\frac{d_2\psi_1}{h}+d_1)+3qr$,

$F = -(1-\frac{c_2kd_2}{h})(d_1+\frac{d_2}{h}\psi_1)+r$,

$\psi_1 = \frac{p}{1+c_1k}$, $\psi_2 = \frac{c\psi_1-d_2q-\theta}{h}$.

One can see that $A > 0$, but the sign of $B, C, D, E$ and $F$ is hard to confirm. By Descartes' Rule of Signs ^[30], we have the following existence criterion of the positive root of the equation in view of the sign of the coefficient.

(Ⅰ)The equation has exactly one positive root if

(ⅰ) $B > 0, \quad C > 0, \quad D > 0, \quad E > 0, \quad F < 0;$

(ⅱ) $B > 0, \quad C > 0, \quad D > 0, \quad E < 0, \quad F < 0;$

(ⅲ) $B > 0, \quad C > 0, \quad D < 0, \quad E < 0, \quad F < 0;$

(ⅳ) $B > 0, \quad C < 0, \quad D < 0, \quad E < 0, \quad F < 0;$

(ⅴ) $B < 0, \quad C < 0, \quad D < 0, \quad E < 0, \quad F < 0.$

(Ⅱ)The equation has at least one positive root if

(ⅰ) $B > 0, \quad C > 0, \quad D < 0, \quad E > 0, \quad F < 0;$

(ⅱ) $B > 0, \quad C < 0, \quad D < 0, \quad E > 0, \quad F < 0;$

(ⅲ) $B < 0, \quad C > 0, \quad D < 0, \quad E > 0, \quad F < 0.$

(Ⅲ)The equation has an even number of positive roots, including none, for the rest of conditions.

Next, we study the local stability for each of the equilibrium points. The linearized form of system (1.2) at an equilibrium $(\bar{x}, \bar{y})$ is as follows:

where

We can derive the characteristic equation corresponding to system (2.1) as

Theorem 2.1. The trivial equilibrium $E_0$ is locally asymptotically stable if $r < d_1$ and unstable if $r > d_1$. When $r > d_1$, the predator-free equilibrium $E_1$ is locally asymptotically stable if $\frac{cp}{1+c_1k} < \frac{ad_2}{r-d_1}+\frac{\theta q(r-d_1)}{a}+d_2q+\theta$, otherwise is unstable.

Proof. At $E_0 = (0, 0)$, $J_{11} = r-d_1$, $J_{22} = -d_2$, $J_{12} = J_{21} = K_{21} = K_{22} = 0$. The eigenvalues are $\lambda_1 = r-d_1$ and $\lambda_2 = -d_2$. Hence, $E_0$ is locally asymptotically stable if $r < d_1$ and unstable if $r > d_1$.

When $r > d_1$, the predator-free equilibrium $E_1$ exists, at which the characteristic equation becomes

Noting that at $E_1$, $J_{11} = d_1-r < 0$, the stability of $E_1$ is determined by equation

Since $K_{22} > 0$, $F'(\lambda) = 1+K_{22}\tau e^{-\lambda \tau} > 0$. Considering that $\lim\limits_{\lambda \rightarrow -\infty} F(\lambda) = -\infty$, the equation has a negative root if

Calculation shows that Eq (2.3) is equivalent to

The proof is complete.

Remark 2.1. When $\theta = 0$, (2.4) reduces to $d_2(1+c_1k) > (cp-d_2q(1+c_1k))\frac{r-d_1}{a}$, which is consistent with the conclusion obtained in ^[25]. Equation (2.4) indicates that the stability of the predator-free equilibrium $E_1$ is related to the value of $\theta$. Unfortunately, as previously analyzed, the presence of $h$ and $\theta$ makes the conditions for the existence of the positive equilibrium less intuitive. Numerical simulations in Section 4 shows that when $E_1$ loss its stability, the positive equilibrium will exist.

We assume that the positive equilibrium $E^*$ of system (1.2) exists, whose characteristic equation is Eq (2.2). When $\tau = 0$, Eq (2.2) reduces to

The roots of Eq (2.5) are

where

Noting that at $E^*$, $J_{22}+K_{22} = -hy^* < 0$, we will consider the following hypotheses:

(H1): $J_{11}J_{22}+J_{11}K_{22}-J_{12}J_{21}-J_{12}K_{21} < 0$;

(H2): $J_{11}J_{22}+J_{11}K_{22}-J_{12}J_{21}-J_{12}K_{21} > 0$ and $J_{11} < hy^*$;

(H3): $J_{11}J_{22}+J_{11}K_{22}-J_{12}J_{21}-J_{12}K_{21} > 0$ and $J_{11} > hy^*$;

(H4): $J_{11}J_{22}-J_{11}K_{22}-J_{12}J_{21}+J_{12}K_{21} > 0$;

(H5): $J_{11}J_{22}-J_{11}K_{22}-J_{12}J_{21}+J_{12}K_{21} < 0$.

Lemma 2.1. Assume that the positive equilibrium $E^*$ of system (1.2) exists.

(1) If (H1) holds, $\lambda_- < 0$ and $\lambda_+ > 0$, which indicates that the positive equilibrium $E^*$ is unstable.

(2) If (H2) is satisfied, Eq (2.5) only has roots with negative real parts, which indicates that the positive equilibrium $E^*$ is locally asymptotically stable.

(3) If (H3) is satisfied, Eq (2.5) possesses roots with positive real parts, which indicates that the positive equilibrium $E^*$ is unstable.

Now, we check whether the roots of Eq (2.2) will cross the pure imaginary axis as $\tau$ increases. When $\tau\neq0$, plugging $\lambda = \omega i \; (\omega > 0)$ into Eq (2.2) to obtain the critical value where Hopf bifurcation occurs. Separating the real and imaginary parts, we obtain

By trigonometric identity, we obtain

Recalling that $J_{22}+K_{22} < 0$, $J_{11}^2+J_{22}^2-K_{22}^2+2J_{12}J_{21} > 0$. Equation (2.7) has no positive root when

Hence no root of Eq (2.2) will cross the pure imaginary axis for all $\tau\geq0$, that is, the stability of the positive equilibrium $E^*$ will not change as $\tau$ increases.

Lemma 2.2. Assume that the positive equilibrium $E^*$ of system (1.2) exists.

(1) If (H1) and (H5) hold, the positive equilibrium $E^*$ is unstable for all $\tau > 0$.

(2) If (H2) and (H4) are satisfied, the positive equilibrium $E^*$ is locally asymptotically stable for all $\tau > 0$.

(3) If (H3) and (H4) are satisfied, the positive equilibrium $E^*$ is unstable for all $\tau > 0$.

If Eq (2.8) is reversed, Eq (2.7) has a unique positive root

where

Then we obtain a sequence of $\tau_n$ corresponding to $\omega_0$:

Hence we obtain possible value for $\tau$ at which Hopf bifurcation may occur. When $\tau = \tau_0$, $\pm i \omega_0$ is a pair of pure imaginary roots of Eq (2.2). For the sake of confirming whether Hopf bifurcation occurs at the first critical value $\tau = \tau_0$, the transversal condition needs to be verified. Differentiating Eq (2.2) with respect to $\tau$ gives

Therefore, at $\tau = \tau_0$, we have

where $k_1 = \omega_0(J_{11}K_{22}-J_{22}K_{22}-2J_{12}K_{21})$ and $k_2 = -2\omega_0^2 K_{22}-(J_{11}+J_{22})(J_{11}K_{22}-J_{12}K_{21})$. Besides, Eq (2.6) gives

Consequently,

The proof about the transversal condition is complete and a Hopf bifurcation exists when $\tau = \tau_0$.

It is worth noting that Eq (2.7) has a unique positive root and $\Re\left(d(\lambda)/d\tau|_{\tau = \tau_0}\right)^{-1} > 0$ under the premise of $(J_{11}J_{22}+J_{11}K_{22}-J_{12}J_{21}-J_{12}K_{21})(J_{11}J_{22}-J_{12}J_{21}+J_{12}K_{21}-J_{11}K_{22}) < 0$. By the insightful work of Cooke and Grossman ^[31], one can see that (ⅰ) The unstable positive equilibrium of system (1.2) never becomes stable; (ⅱ) If the positive equilibrium is stable for $\tau = 0$, it becomes unstable at $\tau_0$ and remains so as $\tau$ is increased; (ⅲ) Only crossing of the imaginary axis from left to right is possible as $\tau$ increases. Therefore, any roots on the right half plane will remain on the right plane and the stability of the positive equilibrium can only be lost but not regained. Then, we obtain the lemma as follows.

LLemma 2.3. Assume that the positive equilibrium $E^*$ of system (1.2) exists.

(1) If (H1) and (H4) hold, the positive equilibrium $E^*$ is unstable for all $\tau > 0$.

(2) If (H2) and (H5) are satisfied, the positive equilibrium $E^*$ is locally asymptotically stable for $0 < \tau < \tau_0$ and unstable for $\tau > \tau_0$.

(3) If (H3) and (H5) are satisfied, the positive equilibrium $E^*$ is unstable for all $\tau > 0$.

In summary, we obtain the following theorem.

Theorem 2.2. Assume that the positive equilibrium $E^*$ of system (1.2) exists and

(I) If (H1) as well either (H4) or (H5) hold, the positive equilibrium $E^*$ is unstable for all $\tau > 0$.

(II)-(i) If (H2) and (H4) are satisfied, the positive equilibrium $E^*$ is locally asymptotically stable for all $\tau > 0$;

(II)-(ii) If (H2) and (H5) are satisfied, the positive equilibrium $E^*$ is locally asymptotically stable for $0 < \tau < \tau_0$ and unstable for $\tau > \tau_0$. The system (1.2) undergoes a Hopf bifurcation arround $E^*$ at $\tau = \tau_0$.

(III) If (H3) as well either (H4) or (H5) hold, the positive equilibrium $E^*$ is unstable for all $\tau > 0$.

Remark 2.2. When $h = 0$ and $\theta = 0$, $J_{11}J_{22}+J_{11}K_{22}-J_{12}J_{21}-J_{12}K_{21} = -J_{12}K_{21} > 0$, then hypothesis (H1) will never hold, which will simplify Theorem 2.2. Furthermore, Eq (2.8) reduces to $2d_2J_{11}J_{12}K_{21}-J_{12}^2K_{21}^2 > 0$, which happens to be Eq (3.40) in ^[25]. Hence, Theorem 3.4 in ^[25] is a special case of Theorem 2.2 under conditions $h = 0$ and $\theta = 0$. The formula of hypotheses indicates that both the value of $h$ and $\theta$ will affect the stability of the positive equilibrium $E^*$. We shall numerically investigate the impact of them in Section 4.

3.   Dynamics of stochastic model

In this section, we mainly study the dynamics of the stochastic system (1.3). We first prove that there exists a unique global positive solution of model system (1.3) with any given positive initial value. Then, we investigate the asymptotic property of the stochastic system (1.3) around the equilibriums of its corresponding deterministic system (1.2).

3.1.   Existence and uniqueness of the global positive solution

The following result is related to the existence and uniqueness of the global positive solution, which is a prerequisite for researching the long term behavior of system (1.3).

Theorem 3.1. For any given initial value Eq (1.4), system (1.3) has a unique global positive solution $(x(t), y(t))$ for $t\geq -\tau$. Furthermore, the solution will remain in $\mathbb{R}_+^2$ with probability one, namely $(x(t), y(t)) \in \mathbb{R}_+^2$ for all $t\geq -\tau$ almost surely.

Proof. Note that the coefficients of system (1.3) are locally Lipschitz continuous and for any given initial value Eq (1.4), there exists a unique local solution $(x(t), y(t))$ on $t \in [-\tau, \tau_e),$ where $\tau_e$ is the explosion time ^[32]. In order to demonstrate that this solution is global, it's sufficient to prove that $\tau_e = \infty$ a.s.. Let $k_0\geq 0$ be a sufficiently large constant for every component of $(x(0), y(0))$ all lying within the interval $[\frac{1}{k_0}, k_0]\times[\frac{1}{k_0}, k_0]$. For each integer $k\geq k_0$, we define the stopping time as follows:

Throughout this paper, we set $\inf \emptyset = \infty$(as usual $\emptyset$ denotes the empty set). Clearly, $\tau_k$ is increasing as $k\rightarrow \infty$. Set $\tau_\infty = \lim\limits_{k \rightarrow \infty} \tau_k$, hence $\tau_\infty\leq\tau_e$ a.s.. To complete the proof, we only need to show that $\tau_\infty = \infty$ a.s.. If this statement is false, there is a pair of constants $T > 0$ and $\varepsilon\in (0, 1)$ such that $P\{\tau_k \leq T\}\geq \varepsilon$ for any $k\geq k_0$.

We define a $C^2-$function $V$ : $\mathbb{R}_+^2\rightarrow \mathbb{R}_+$ as follows:

Applying $\mathrm{It\hat{o}}$ formula yields

where

where $K$ is a positive number, the remainder of the proof follows that in Li et al. ^[33], here, we omit it. The proof is complete.

3.2.   Asymptotic Property

Due to the interference of stochastic noise, system (1.3) does not exist any equilibriums. However, we are interested in whether system (1.3) has stability and what is the influence of white noise. In this subsection, we discuss the asymptotic behaviors of system (1.3) around the predator-free equilibrium and the positive equilibrium of the corresponding deterministic system (1.2).

Definition 3.1. Let $(x(t), y(t))$ be a solution of model (1.3) with initial value (1.4). $\bar{E} = (\bar{x}, \bar{y})$ is an equilibrium point of the corresponding deterministic system (1.2). If there exists a constant $\bar{K} > 0$ such that

then we say that the solution of system (1.3) will be swinging near the equilibrium $\bar{E} = (\bar{x}, \bar{y})$.

Next, we study the asymptotic property of system (1.3) around the predator-free equilibrium $E_1 = (\frac{r-d_1}{a}, 0)$ when $r-d_1 > 0$ and the positive equilibrium, respectvely.

Theorem 3.2. Let $(x(t), y(t))$ be a solution of model (1.3) with initial value (1.4). If the condition $r-d_1 > \sigma_1^2$ holds, then

where

Proof. Define the function

with

$V_1 = \frac{(x-x_0)^2}{2},$

$V_2 = x-x_0-x_0ln\frac{x}{x_0},$

$V_3 = y+\frac{cp}{1+c_1k}\int_{t-\tau}^t \frac{x(s)y(s)}{1+qx(s)}ds,$

$V_4 = cx+y(t+\tau)+d_2\int_t^{t+\tau}y(s)ds.$

where $\omega_i$(i = 1, 2, 3) are positive constants to be determined later.

Applying $\mathrm{It\hat{o}}$ formula yields:

For

where

Using the inequality $(a+b)^2 \leq 2a^2+2b^2$ for all $a, b\in R$, we have

For

where

For

where

For

where

Therefore, we have

Choosing $\omega_1 = \frac{c}{2x_0}$, $\omega_2 = \frac{c}{2}$, $\omega_3 = \frac{cpx_0}{d_2(1+c_1k)}$, one can see that

Integrating both sides of Eq (3.1) from $0$ to $t$ and taking the expectation yield

where

Divide both sides by $t$ and take the limit superior, and then we obtain

where

This completes the proof.

Remark 3.1. From Theorem 3.2, one can see that the solution of system (1.3) will be swinging near the predator-free equilibrium under certain conditions. In addition, the value of $\sigma_2^2$ has no effect on the swing intensity, but the swing intensity and $\sigma_1^2$ are positively correlated.

Theorem 3.3. Let $(x(t), y(t))$ be a solution of system (1.3) with initial value (1.4). If the following conditions are satisfied:

then

where

Proof. Recalling that $E^* = (x^*, y^*)$ is the positive equilibrium of system (1.2), we have

Define

An application of $\mathrm{It\hat{o}}$ formula implies

where

Using the inequality $(a+b)^2 \leq 2a^2+2b^2$ for all $a, b\in R$, we have

Set

where

Define

where

where the inequality $a(x-\frac{r}{a})^2\ge0$ is applied.

Define

where

Define

where

Moreover, define

Define

where

Additionally, define

Construct the function

where $\omega_i(i = 1, \cdots, 4)$ are positive constants to be determined as follows.

Therefore, we have

Choosing $\omega_1 = \frac{c^2p^2}{ax^*h(1+c_1k)^2(1+qx^*)^2}$, $\omega_4 = \frac{hq(1+c_1k)}{3cp}$, $\omega_3 = \frac{1}{d_2}(\frac{2\omega_1px^*}{q(1+c_1k)}+\frac{cp}{q(1+c_1k)}+\theta x^*+\omega_4\theta x^*y^*)$, $\omega_2 = \frac{1}{r}(\omega_1(\frac{c_2kx^*y^*}{1+c_2ky^*}+\frac{2py^*}{q(1+c_1k)}+\frac{c_2kry^*}{q(1+c_2ky^*)})+\omega_3cr+\frac{cp}{1+c_1k}\frac{y^*}{1+qx^*}+\theta y^*+\omega_4\theta {y^*}^2)$, then

Integrating both sides from $0$ to $t$ and taking the expectation yield

where

Divide both sides by t and take the limit superior, and then we obtain

where

This completes the proof.

Remark 3.2. From Theorem 3.3, one can see that the solution of system (1.3) will be swinging near the positive equilibrium under certain conditions. In addiction, there is a positive correlation between the swing intensity and the intensity of environmental noise.

4.   Numerical simulations

In this section, we carry out some numerical simulations to illustrate the main analytic results obtained in Sections 2 and 3, respectively. We will focus on the impacts of the response level of fear $k$, the self-defense level $\theta$ and digestion delay $\tau$ on the dynamic of bio-system.

Firstly, we consider the initial value $(x(0), y(0)) = (0.2, 0.1)$ with parameter values, which are fixed unless other specified:

Notice that from Eq (2.4), we can obtain the threshold for $\theta$ and $k$ in comparison with a formula composed of other parameters for the stability of the predator-free equilibrium $E_1$. Simplifying Eq (2.4) gives $\theta^* = \frac{a}{a+q(r-d_1)}(\frac{cp}{1+c_1k}-\frac{ad_2}{r-d_1}-d_2q)$ and $k^* = \frac{1}{c_1}(\frac{acp(r-d_1)}{a^2d_2+\theta q(r-d_1)^2+a(d_2q+\theta)(r-d_1)}-1)$. Then, Theorem 2.1 can be simply expressed as: when $\theta > \theta^*$ and $k > k^*$, the predator-free equilibrium $E_1$ is locally asymptotically stable. Besides, for the sake of convenience, we denote $J_{11}J_{22}+J_{11}K_{22}-J_{12}J_{21}-J_{12}K_{21}$ and $J_{11}J_{22}-J_{12}J_{21}+J_{12}K_{21}-J_{11}K_{22}$ as $\Phi_1$ and $\Phi_2$, respectively. Choosing $\theta = 0.01$, we obtain $k^* = 27.8788$. By Theorem 2.1, when $k = 30 > k^*$, the predator-free equilibrium $E_1$ is locally asymptotically stable (as demonstrated in Figure 1(a)). For $k = 20 < k^*$, the predator-free equilibrium $E_1$ is unstable, but numerical calculation gives that $\Phi_1 = 2.35\times10^{-4}$, $J_{11} = -0.0068 < hy^* = 6.65\times10^{-4}$ and $\Phi_2 = 2.67\times10^{-4}$ (Theorem 2.2-(Ⅱ)-(ⅰ)), and hence, $E^* = (0.9392, 0.0065)$ is locally asymptotically stable for all $\tau > 0$, as demonstrated in Figure 1(b). Furthermore, when $k = 30$ and $\theta$ decreases to $\theta = 0.005 < \theta^* = 0.0084$, the predator-free equilibrium $E_1$ is unstable, while $\Phi_1 = 1.14\times10^{-4}$, $J_{11} = -0.0126 < hy^* = 4.12\times10^{-4}$ and $\Phi_2 = 0.0011$ (Theorem 2.2-(Ⅱ)-(ⅰ)), which indicates that $E^* = (1.3895, 0.0412)$ is locally asymptotically stable for all $\tau > 0$ (as demonstrated in Figure 1(c)), as well. The comparison of the three pictures in Figure 1 can verify Remark 2.1, that is, when $E_1$ is unstable, the positive equilibrium exists.

The numerical simulations presented above indicate that the response level of fear $k$ and the self-defense level $\theta$ have important impacts on the survival of the population, especially on predator. When the coexistence equilibrium $E^*$ exists, taking $\theta = 0.01$ and $k = 20$ for instance, the stability of $E^*$ will still be affected the response level and the self-defense level. By fixing $\theta = 0.01$ and decreasing $k$ to $k = 10$, we obtain $\Phi_1 = 4.71\times10^{-4}$, $J_{11} = -0.0017 < hy^* = 6.87\times10^{-4}$ and $\Phi_2 = -4.55\times10^{-4}$, corresponding to Theorem 2.2-(Ⅱ)-(ⅱ). In this case, the stability of positive equilibrium $E^* = (0.4634, 0.0687)$ is related to the value of $\tau_0$. Further calculations reveal that $\tau_0 = 5.7998$. The numeric solutions are shown in Figure 2, for $\tau < \tau_0$ and $\tau > \tau_0$, respectively in (a) and (b); and the case of Figure 2(b) plotted in the $x-y$ plane is given in Figure 4(a). Analogously, we fix $k = 20$ and change $\theta$ to $\theta = 1\times10^{-4}$. Computations give $\Phi_1 = 3.64\times10^{-4}$, $J_{11} = -0.0043 < hy^* = 7.69\times10^{-4}$ and $\Phi_2 = -9.15\times10^{-6}$, meaning that stability of positive equilibrium $E^* = (0.6942, 0.0769)$ depends on the value of delay $\tau$ (Theorem 2.2-(Ⅱ)-(ⅱ)). We compute to obtain $\tau_0 = 227.4388$ by (2.9). The numeric solutions are illustrated in Figure 3, for $\tau < \tau_0$ and $\tau > \tau_0$, respectively in (a) and (b); and the case of Figure 3(b) plotted in the $x-y$ plane is given in Figure 4(b).

In order to show the role of $\tau$ on results more intuitively, we choose $\theta = 0.01$ and $k = 10$, and consider the set of values $0.1$, $2$, $8$, $15$, and $30$ (as demonstrated in Figure 5). The results indicate that when $\tau > \tau_0 = 5.7998$, the solution of system (1.2) shows more oscillating behaviors and the positive equilibrium is unstable. However, once the parameter is reduced to less than the threshold, the system solution will quickly converge to the equilibrium point, which verifies Theorem 2.2-(Ⅱ)-(ⅱ)) very well.

To explore the impact of self-defence on the dynamics of system (1.2) and to compare with the results obtained in ^[25], we choose the same parameter values as in ^[25]: $h = 0$, $k = 1$, $d_2 = 0.05$, $c = 0.4$, $c_1 = 1$, $c_2 = 1$ and $\tau = 120$, and consider the set of values $0$, $0.01$, $0.03$, $0.036$, and $0.04$ (as demonstrated in Figure 6). The curve of $\theta = 0$ is in line with that in ^[25] and a periodic solution occurs. While as $\theta$ increases, the oscillation amplitude of both the prey and the predator gradually decreases and the equilibrium tends to stable eventually. Hence the level of self-defence acts as a stabilizing factor in system (1.2) (Remark 2.2).

Furthermore, Theorem 2.2 demonstrates that when $\Phi_1 > 0$ and $J_{11} > hy^*$, the positive equilibrium $E^*$ of system (1.2) is unstable for all $\tau > 0$. We set $a = 0.05, q = 0.8, \theta = 0.01$ and $k = 2$ to ensure that the prerequisites are established and the numeric solutions are demonstrated in Figure 7. In such a case, the solutions are periodic and as $\tau$ increases, the amplitude of the oscillation (periodic solution) increases.

For the stochastic model, we shall adopt Milstein's Higher Order Method mentioned in ^[29] to verify our theoretical results. The corresponding discretization equations are

where the time increment $\Delta t > 0$, $\eta$ is the integer part of $\frac{\tau}{\Delta t}-1$, $\xi_k$ and $\varrho_k (k = 1, 2, ..., n)$ are independent Gaussian random variables that follow the standard normal distribution $N(0, 1).$

In Figure 8, we consider the initial value $(x(0), y(0)) = (0.2, 0.1)$ with related parameter values: $\theta = 0.01$, $k = 30$, $\sigma_i = 0.01(i = 1, 2)$ and other parameters are given in Eq (4.1). In this case, $r-d_1 = 0.02 > \sigma_1^2 = 1\times10^{-4}$, which means the condition in Theorem 3.2 is satisfied. Hence, all positive solutions of system (1.3) will fluctuate around the predator-free equilibrium. As clearly shown in Figure 8, the curve of the stochastic model goes around $E_1 = (2, 0)$, which supports Theorem 3.2.

Besides, we take the parameters in Figure 1(b) and $\sigma_i = 0.01(i = 1, 2)$. In this situation, the conditions in Theorem 3.3 are satisfied:

According to Theorem 3.3, all positive solutions of system (1.3) will fluctuate around the positive equilibrium of system (1.2). Figure 9 illustrates that the curve of the stochastic model goes around $E^* = (0.9392, 0.0665)$, which validates the analytical result.

Theorem 3.2 indicates that when the intensities of stochastic perturbation are less than a certain value, the solution of system (1.3) will go around the solution of system (1.2). We are also interested in what happens if the intensities of external disturbance exceed this threshold. So under the condition that $\theta = 0.01$ and $k = 20$ and other parameters are given in Eq (4.1), we choose the different intensities of white noise $\sigma_1 = 0.01, 0.05, 0.2$ (with $\sigma_2 = 0.01$ fixed) and $\sigma_2 = 0.01, 0.05, 0.2$ (with $\sigma_1 = 0.01$ fixed) to explore their influence on the dynamic of system (1.3). The simulation results are observed in Figures 10 and 11. The trajectories of $\sigma_i = 0.01, 0.05 \; (i = 1, 2)$ demonstrate that the small noise will not change the stability of the equilibrium of system (1.2), but the amplitudes of both prey and predator increase with the increase of the intensities of white noise. If the environmental fluctuations keep increasing, the population affected will go to extinction (such as $\sigma_i = 0.2 \; (i = 1, 2)$), which can be observed clearly from Figures 10(a) and Figure 11(b). Notice that Theorem 2.2 claims the co-existence of predator and prey. This discrepancy highlights the impact of stochastic fluctuations on the predator-prey model. What's more, the numerical simulations also verify a basic fact that the survival of predator is based on the survival of the prey while the prey can persist in the absence of the predator.

Ultimately, we focus on the impacts of the response level $k$ and the self-defense level $\theta$ on the prey, respectively. As stated at the beginning of the paper, we firmly believe that the scale of prey must have a significant relationship with anti-predation response and self-defence. Here, under the condition that $\sigma_i = 0.01 \; (i = 1, 2)$ and other factors are given in Eq (4.1), we choose different response level $k = 1, 10, 20$ (with $\theta = 0.01$ fixed) and self-defense level $\theta = 0.001, 0.01, 0.05$ (with $k = 20$ fixed). As is clearly demonstrated in Figure 12, both the response level and the self-defense level exert far-reaching influence on the scale of prey. With the increase of the self-defense ability and anti-predator consciousness, the population of prey will gradually increase and trend to the predator-free equilibrium, which has two implications: (ⅰ) High-level self-defense will prevent the predator from being preyed on by the predator, otherwise the predator will eventually become extinct; (ⅱ) Compared with the negative effects, the benefits of high-level response are more significant. One of the possible reasons is that the increase in escaped prey compensates for the decline in newborns. These are consistent with what we know: predators prefer safer and more efficient predation.

5.   Conclusions

This paper is mainly related to a deterministic predator-prey model with fear factor and self-defence, considering the existence of digestion delay, and its corresponding stochastic model. Rather than simply considering the negative impact of the fear effect, we take the benefits of anti-predation behavior and resistance into consideration. Both avoidance and resistance of prey along with digestion delay all affect the stability of the equilibrium point of the deterministic system and may lead to Hopf bifurcation. In addition, as researchers have not paid much attention to the corresponding stochastic model, we mainly investigate the asymptotic behavior around the solution of the deterministic model. In contrast with existing results, our research shows that the level of fear alone is not enough to completely determine the dynamics of the predator-prey model. The level of self-defence acts as a stabilizing factor, as well. The results obtained can be applied to keep the coexistence of biological populations and maintain ecological balance. Furthermore, numerical results support and expand our analyses, indicating that: (ⅰ) The introduction of disturbance of noise may drift the coexistence equilibrium to the predator-free equilibrium; (ⅱ) Small noise will not change the stability of the equilibrium of system (1.3) and the amplitude of the population is positively correlated with the intensity of withe noise; (ⅲ) Large environmental fluctuations will lead to the extinction of predator and prey; (ⅳ) The scale of prey has a significant relationship with both anti-predation behavior and self-defence, so is predator.

In addition to the conclusions that have been reached, there are some possible extensions of our work. {As mentioned before, we are of great interest in the impact of self-defence on the prey population. A replacement of $\frac{1}{1+c_1k}$ (respectively, $\frac{r}{1+c_2ky(t)}$) by $g(k, \theta)$ (respectively, $f(k, \theta, y(t))$) which is decreasing in both $k$ and $\theta$ may be a meaningful attempt. Furthermore, } as pointed out in Xiao and Chen ^[34], stage structure has a significant influence in the stability of periodic orbits. The probability of self-defense is different when the prey encounters immature or mature predators. Moreover, Lan et al. ^[35] proposed a model with psychological effects and impulsive toxicants in polluted environments. These models are more complex but realistic. We can consider models with factors mentioned above for further investigations.

Acknowledgments

The research was supported by National Natural Sciences Foundation of China (11971405, 22072057).

Conflict of interest

The authors declare there is no conflict of interest.