**3. Occurrence of the Paradox of Enrichment**

According to Jensen and Ginzburg [11], the paradox of enrichment is accepted intuition and must be considered a theory in ecology. In this section, we study the occurrence of the paradox of enrichment on systems (1) and (2) with Holling types I and II. To study this phenomenon, we discuss the stability of the coexistence equilibrium points *E* = (*<sup>x</sup>*, *y*) and *E* = (*<sup>x</sup>*, *y*, *z*) in two and three trophics, respectively.

### *3.1. Occurrence of the Paradox of Enrichment with Holling Type I*

To check this phenomenon with Holling type I, we found the coexistence equilibrium points of systems (1) and (2) and present some theorems that prove the occurrence of the paradox of enrichment in systems (1) and (2) with Holling type I.

The coexistence equilibrium point of system (1) with Holling type I is *E* = (*x*ˆ, *y*ˆ) = *<sup>k</sup>*(*u*+*e*) *e*+*keα* , *ke<sup>α</sup>*−*<sup>u</sup> eα*+*eα*2*<sup>k</sup>*. It exists (positive equilibrium point) under the following condition:

$$
\ker \mathfrak{u} \gg \mathfrak{u} \tag{3}
$$

*E* = (*<sup>x</sup>*, *y*, *z*) represents the coexistence of the equilibrium point of system (2) with Holling type I, which is obtained through the positive solution of the following algebraic system:

$$\begin{aligned} 1 - \frac{\mathbf{x}}{k} - \alpha \mathbf{y} - \beta \mathbf{z} &= 0 \\\\ -\mathbf{u} + \mathbf{e}\_1 \alpha \mathbf{x} - \mathbf{e}\_1 \alpha \mathbf{y} - \mathbf{c}\_1 \mathbf{z} &= 0 \\\\ -\mathbf{w} + \mathbf{e}\_2 \beta \mathbf{x} - \mathbf{e}\_2 \beta \mathbf{z} - \mathbf{c}\_2 \mathbf{y} &= 0 \end{aligned} \tag{4}$$

**Theorem 1.** *The coexistence equilibrium point E*ˆ = (*x*ˆ, *y*ˆ) = *<sup>k</sup>*(*u*+*e*) *e*+*keα* , *ke<sup>α</sup>*−*<sup>u</sup> eα*+*eα*2*<sup>k</sup> of system (1) is globally asymptotically stable within the positive quadrant of the x* − *y plane.*

**Proof.** Let *<sup>G</sup>*(*<sup>x</sup>*, *y*) = 1*xy* . *G* is a Dulac function. It is continuously differentiable in the positive quadrant of the *x* − *y* plane *A* = {{(*<sup>x</sup>*, *y*)|*x* > 0, *y* > 0} Hsu [43].

$$\begin{aligned} \mathcal{N}\_1(\mathfrak{x}, \mathfrak{y}) &= \mathfrak{x} \left( 1 - \frac{\mathfrak{x}}{k} \right) - \mathfrak{a} \mathfrak{x} \mathfrak{y}, \\\\ \mathcal{N}\_2(\mathfrak{x}, \mathfrak{y}) &= -\mathfrak{u}\mathfrak{y} + \mathfrak{e}\_1 \mathfrak{a} \mathfrak{x} \mathfrak{y} - \mathfrak{e}\_1 \mathfrak{a} \mathfrak{y}^2. \end{aligned}$$

.

> Thus, <sup>Δ</sup>(*GN*1, *GN*2) = *∂*(*GN*1) *∂x* + *∂*(*GN*2) *∂y* = −1 *yk* − *e*1*αx* .

It is observed that Δ(GN1, GN2) is not identically zero and does not change sign in the positive quadrant of the *x* − *y* plane. Per the Bendixson–Dulac criterion, there is no periodic solution inside the positive quadrant of the *x* − *y* plane. *E*2 is globally asymptotically stable inside the positive quadrant of the *x* − *y* plane.

### **Theorem 2.** *The coexistence equilibrium point E* = (*<sup>x</sup>*, *y*, *z*) *of system (2) is globally asymptotically stable.*

**Proof.** The global stability of positive equilibrium point E is proved by using Lyapunov function.

$$V = B\_1(\mathbf{x} - \overline{\mathbf{x}} - \ln\left(\frac{\mathbf{x}}{\overline{\mathbf{x}}}\right)) + B\_2(\mathbf{y} - \overline{\mathbf{y}} - \ln\left(\frac{\mathbf{y}}{\overline{\mathbf{y}}}\right)) + B\_3(z - \overline{z} - \ln\left(\frac{z}{\overline{z}}\right)).\tag{5}$$

. Differentiating *V* with respect to time along the solutions of the system (5)

$$\begin{split} \frac{dV}{dt} &= B\_1(\mathbf{x} - \mathbf{\mathcal{T}}) \Big[ \Big( 1 - \frac{\mathbf{x}}{k} - a\boldsymbol{y} - \beta \mathbf{z} \Big) - \Big( 1 - \frac{\mathbf{x}}{k} - a\mathfrak{J} - \beta \mathfrak{T} \Big) \Big] + B\_2(\mathbf{y} - \mathfrak{J}) \big( -u + \epsilon\_1 a\mathbf{x} - \epsilon\_1 a\mathfrak{y} - \epsilon\_1 \mathbf{z} \Big) \\ &- \Big( -u + \epsilon\_1 a\mathfrak{T} - \epsilon\_1 a\mathfrak{T} - \epsilon\_1 \mathfrak{T} \Big) \big( -w + \epsilon\_2 \beta \mathfrak{x} - \epsilon\_2 \beta \mathfrak{z} - \epsilon\_2 y \big) - \Big( -w + \epsilon\_2 \beta \mathfrak{T} - \epsilon\_2 \beta \mathfrak{T} - \epsilon\_2 \mathfrak{y} \Big) \Big). \end{split} \tag{6}$$

$$\begin{split} \frac{dV}{dt} &= B\_1(\mathbf{x} - \overline{\mathbf{z}}) \left[ \frac{-(\mathbf{x} - \overline{\mathbf{z}})}{k} - a(y - \overline{y}) - \beta(\mathbf{z} - \overline{\mathbf{z}}) \right] + B\_2(y - \overline{y}) [\epsilon\_1 a(\mathbf{x} - \overline{\mathbf{z}}) - \epsilon\_1 a(y - \overline{y}) - \epsilon\_1(\mathbf{z} - \overline{\mathbf{z}})] \\ &+ B\_3(\mathbf{z} - \overline{\mathbf{z}}) [\epsilon\_2 \beta(\mathbf{x} - \overline{\mathbf{z}}) - \epsilon\_2 \beta(\mathbf{z} - \overline{\mathbf{z}}) - \epsilon\_2(y - \overline{y})]. \end{split} \tag{7}$$

*Symmetry* **2018**, *10*, 532

$$\begin{split} \frac{dV}{dt} &= B\_1 \left[ \frac{-(\mathbf{x} - \overline{\mathbf{x}})^2}{k} - a(\mathbf{x} - \overline{\mathbf{x}})(y - \overline{y}) - \beta(\mathbf{x} - \overline{\mathbf{x}})(z - \overline{z})} \right] + B\_2 \begin{bmatrix} e\_1 a(y - \overline{y})(\mathbf{x} - \overline{\mathbf{x}}) - e\_1 a(y - \overline{y})^2 \\ -c\_1 (y - \overline{y})(z - \overline{z}) \end{bmatrix} \\ &+ B\_3 \left[ c\_2 \beta (z - \overline{z}) (\mathbf{x} - \overline{\mathbf{x}}) - c\_2 \beta (z - \overline{z})^2 - c\_2 (z - \overline{z}) (y - \overline{y}) \right]. \end{split} \tag{8}$$

By selecting *B*1 = 1, *B*2 = 1 *e*1 , and *B*3 = 1 *e*2 , so

$$\frac{dV}{dt} = -\frac{1}{k}(\mathbf{x} - \overline{\mathbf{x}})^2 - a(y - \overline{y})^2 - c\_1(y - \overline{y})(z - \overline{z}) - \beta(z - \overline{z})^2 - c\_2(z - \overline{z})(y - \overline{y})\tag{9}$$

We conclude that *dVdt* is a negative definite without any conditions (i.e., no constraints on parameters).

In this section, the main result shows that the dynamic behaviors of systems (1) and (2) are always stable and there is no bifurcation under any conditions through Theorems 1 and 2. Therefore, the paradox of enrichment in systems (1) and (2) with Holling type I does not occur.

The system is stable, or population oscillations with small amplitude are likely to occur if the defense of prey is effective when compared with the predator's attacking [6].

### *3.2. Occurrence of the Paradox of Enrichment with Holling Type II*

We found the coexistence equilibrium points of systems (1) and (2) and present some theorems that prove the occurrence of the paradox of enrichment in these systems with Holling type II.

The coexistence equilibrium point *E*´ = (*x*´, *y*´) is obtained for system (1) with Holling type II through the positive root of the quadratic equation

$$\dot{\mathbf{x}}^2 + \left(\frac{1}{h\_1} - \frac{u}{c\_1} - \frac{1}{k} + \frac{1}{h\_1 a}\right)\dot{\mathbf{x}} - \left(\frac{1}{h\_1 a} + \frac{u}{c\_1 h\_1 a}\right) = 0\tag{10}$$

and

$$\mathcal{Y} = \frac{1}{a} \left( 1 - \frac{\p}{k} \right) (1 + h\_1 a \pounds) \tag{11}$$

The coexistence equilibrium point = *E* = = *x*, = *y*, = *z* of system (2) with Holling type II is obtained through the positive solution of the following algebraic system:

$$1 - \frac{x}{k} - \frac{ay}{1 + h\_1ax} - \frac{\beta z}{1 + h\_2\beta x} = 0$$

$$-u + \frac{c\_1ax}{1 + h\_1ax} - \frac{c\_1a}{1 + h\_1ax}y - c\_1z = 0\tag{12}$$

$$-w + \frac{c\_2\beta x}{1 + h\_2\beta x} - \frac{c\_2\beta}{1 + h\_2\beta x}z - c\_2y = 0$$

To check this phenomenon with Holling type II, we present the following theorems:

**Theorem 3.** *The coexistence equilibrium point E*´ = (*x*´, *y*´) *is asymptotically stable under the following condition:*

$$k < \frac{\circlearrowleft}{\frac{h\_1 a^2 \circlearrowleft y}{\left(1 + h\_1 a \circlearrowleft x\right)^2} - \frac{c\_1 a \circlearrowleft y}{1 + h\_1 a \circlearrowleft x}}\tag{13}$$

**Proof.** The variational matrix of coexistence point *E*´ is as follows:

$$
\dot{V} = \begin{pmatrix}
\dot{\mathfrak{x}} (-\frac{1}{k} + \frac{h\_1 a^2 \mathfrak{y}}{\left(\frac{1 + h\_1 a \mathfrak{x} \mathfrak{x}}{\left(1 + h\_1 a \mathfrak{x} \right)^2}\right)} & -\dot{\mathfrak{x}} (\frac{a}{1 + h\_1 a \mathfrak{x} \mathfrak{x}}) \\
\dot{\mathfrak{y}} (\frac{c\_1 a + h\_1 c\_1 a^2 \mathfrak{y}}{\left(1 + h\_1 a \mathfrak{x} \right)^2}) & -\dot{\mathfrak{y}} (\frac{c\_1 a}{1 + h\_1 a \mathfrak{x}})
\end{pmatrix}
$$

Through the variational matrix, the equilibrium point *E*´ is locally asymptotically stable, providedthe following condition holds:

$$k < \frac{\dot{\mathfrak{x}}}{\frac{h\_1 \alpha^2 \circ \dot{\mathfrak{y}}}{\left(1 + h\_1 \alpha \mathfrak{x}\right)^2} - \frac{c\_1 \alpha \circ \dot{\mathfrak{y}}}{1 + h\_1 \alpha \mathfrak{x}}}$$

**Corollary 1.** *If condition (13) is not satisfied, then the coexistence equilibrium point E*´ = (*x*´, *y*´) *is unstable.*

Through Corollary 1, there is a destabilization of the coexistence equilibrium point *E*´ according to the carrying capacity parameter, so the paradox of enrichment in system (1) with Holling type II would occur. Therefore, the paradox of enrichment occurs in system (1) with Holling type II through some numerical simulations.

**Theorem 4.** *The coexistence equilibrium point* = *E* = = *x*, = *y*, = *z of system (2) is obtained through the positive solution of system (12). It is locally asymptotically stable proven that conditions (15), (16), and (17) hold.*

The variational matrix of *E* is as follows:

=

$$
\begin{bmatrix}
\bar{\boldsymbol{x}} \\
\bar{\boldsymbol{w}} \\
\bar{\boldsymbol{w}}
\end{bmatrix} = \begin{bmatrix}
\bar{\boldsymbol{x}}(-\frac{1}{k} + \frac{h\_{1}a^{2}\bar{\boldsymbol{y}}}{(1+h\_{1}a\bar{\boldsymbol{x}})} + \frac{h\_{2}\beta^{2}\bar{\boldsymbol{z}}}{(1+h\_{2}\beta\bar{\boldsymbol{x}})} & -\bar{\boldsymbol{x}}(\frac{a}{1+h\_{1}a\bar{\boldsymbol{x}}} - \bar{\boldsymbol{x}}(\frac{\beta}{1+h\_{2}\beta\bar{\boldsymbol{x}}}) \\
\bar{\boldsymbol{y}}(\frac{c\_{1}\boldsymbol{x}+h\_{1}c\_{1}a^{2}\bar{\boldsymbol{y}}}{(1+h\_{1}a\bar{\boldsymbol{x}})} & -\bar{\boldsymbol{y}}(\frac{c\_{1}a}{1+h\_{1}a\bar{\boldsymbol{x}}} - \bar{\boldsymbol{c}}(\frac{\bar{\boldsymbol{y}}}{\bar{\boldsymbol{y}}})) \\
\bar{\boldsymbol{z}}(\frac{c\_{2}\beta\bar{\boldsymbol{x}}+c\_{2}a\beta\bar{\boldsymbol{y}}^{2}}{(1+h\_{2}\beta\bar{\boldsymbol{x}})} & -c\_{2}\bar{\boldsymbol{z}} & -\bar{\boldsymbol{z}}(\frac{c\_{2}\beta}{1+h\_{2}\beta\bar{\boldsymbol{x}}}) \\
\bar{\boldsymbol{W}} & \bar{\boldsymbol{W}} & \bar{\boldsymbol{W}}
\end{bmatrix}
$$

where

$$h\_{11} = \overset{=}{x}(-\frac{1}{k} + \frac{h\_1 a^2 \overset{=}{y}}{(1 + h\_1 a \overset{=}{x})^2} + \frac{h\_2 \beta^2 \overset{=}{z}}{(1 + h\_2 \beta \overset{=}{x})^2}), \ h\_{12} = -\overset{=}{x}(\frac{\alpha}{1 + h\_1 a \overset{=}{x}}), h\_{13} = -\overset{=}{x}(\frac{\beta}{1 + h\_2 \beta \overset{=}{x}}),$$

$$h\_{21} = \overset{=}{y}(\frac{e\_1 a + h\_1 e\_1 a^2 \overset{=}{y}}{(1 + h\_1 a \overset{=}{x})^2}), \ h\_{22} = -\overset{=}{y}(\frac{e\_1 a}{1 + h\_1 a \overset{=}{x}}), h\_{23} = -e\_1 \overset{=}{y}$$

$$h\_{31} = \overset{=}{z}(\frac{e\_2 \beta + e\_2 h\_2 \beta^2 \overset{=}{z}}{(1 + h\_2 \beta \overset{=}{x})^2}), \ h\_{32} = -e\_2 \overset{=}{z}, h\_{33} = -\overset{=}{z}(\frac{e\_2 \beta}{1 + h\_2 \beta \overset{=}{x}}),$$

The characteristic equation of the variational matrix *V* is as follows:

$$
\lambda^3 + H\_1 \lambda^2 + H\_2 \lambda + H\_3 = 0 \tag{14}
$$

$$
H\_1 = -(h\_{11} + h\_{22} + h\_{33})
$$

$$
H\_2 = \left( \begin{array}{c} h\_{11}h\_{22} + h\_{23}h\_{32} + h\_{11}h\_{33} + h\_{22}h\_{33} - h\_{12}h\_{21} - h\_{13}h\_{31} \right) \end{array}
$$

$$
H\_3 = \left( h\_{13}h\_{31}h\_{22} + h\_{12}h\_{21}h\_{33} + h\_{11}h\_{23}h\_{32} - h\_{11}h\_{22}h\_{33} - h\_{13}h\_{21}h\_{32} - h\_{11}h\_{22}h\_{33} \right)
$$

=

.

According to Routh–Hurwitz criterion, = *E* = = *x*, = *y*, = *z* is locally asymptotically stable if it holds the following conditions:

$$H\_1 > 0\tag{15}$$

$$H\_3 > 0\tag{16}$$

$$H\_1 H\_2 > H\_3 \tag{17}$$

=

**Theorem 5.** *If one of the conditions (15)–(17) is not satisfied, then the coexistence equilibrium point* = *E* = = *x*, = *y*, = *z is unstable.*

**Proof.** Through the variational matrix of coexistence point *E* = = *x*, = *y*, = *z* , the stability is satisfied according to the Routh–Hurwitz criterion if all the conditions (15), (16), and (17) must be satisfied. However, it is observed from the first condition that

$$H\_1 > 0 \text{ where } H\_1 = -(h\_{11} + h\_{22} + h\_{33})$$

.

> It is observed that *H*1 < 0 when

$$\frac{h\_1 a^2 \bar{y}}{\left(1 + h\_1 a \bar{\mathbf{x}}\right)^2} + \frac{h\_2 \beta^2 \bar{z}}{\left(1 + h\_2 \bar{\beta} \bar{\mathbf{x}}\right)^2} > \frac{\bar{\mathbf{x}}}{k} + \frac{c\_1 a \bar{\mathbf{y}}}{1 + h\_1 a \bar{\mathbf{x}}} + \frac{c\_2 \beta \bar{z}}{1 + h\_2 \bar{\beta} \bar{\mathbf{x}}} \tag{18}$$

=

Thus, the Routh–Hurwitz criterion is not satisfied, so the equilibrium point *E* is unstable.

In this section, we concluded that the linearity and nonlinearity of functional and numerical responses plays important roles in the occurrence of the enrichment paradox. Some studies have shown that functional and numerical responses led to qualitative differences in dynamic behaviors of prey–predator systems [44–47].

Many biological factors control the shape of the functional and numerical responses as foraging theory and densities of prey, as shown in Nowak et al. [7]. Consequently, the shape of the functional and numerical responses affect the dynamic behaviors of prey–predator systems to be steady state, limit cycles, or complex dynamical behaviours.
