Risk assessment of waterborne infections in Enugu State, Nigeria: Implications of household water choices, knowledge, and practices

1.   Introduction

Mathematically, there are some complex phenomena in nature that are appealing to be investigated. One example is the interaction between prey and its predator. The first predator-prey interaction model was developed independently by ^[1] and ^[2], and therefore it is known as the Lotka-Volterra predator-prey model. Since then, the model for the interaction between prey and its predator is continuously studied and developed by both mathematicians and biologists ^[3]. One of important features that determine the characteristics of predator-prey interaction is the functional response which explains how many prey are eaten by predators per unit time. For example, the famous Rosenzweig-MacArthur predator-prey model implements the Holling-Type Ⅱ functional response $p(x) = \frac{ax}{b+x}$, where $a$ and $b$ are positive parameters related to the predator's attack rate and the handling time, respectively ^[4,5]. The Rosenzweig-MacArthur model has attracted many scholars to study the effect of various factors on predator-prey interactions, for example the effect of refuge ^[6,7,8,9], harvesting ^[10], quiescence ^[11], disease transmission among predators ^[12] or among prey ^[13] and the mutant population ^[14].

The transition between immature and mature stages from an adaptive dynamics perspective has recently been evaluated in ^[15]. In the predator-prey interaction, the age or size of the prey certainly determines the level of predation. Even in many cases, predators are only able to attack when the prey is still immature. Indeed, small size or immature prey tend to be more easily captured and eaten by predators. Moreover, there are many factors that significantly influence the dynamics of the predator–prey models with stage–structure. For example, ^[16] and ^[17] have studied the harvesting policy in a stage–structure predator–prey model. In ^[18], a ratio-dependent prey–predator model with stage–structure was investigated, and the uniform persistence and impermanence of the model were carried out through the sufficient conditions. Other features which influence the interaction of predator-prey with stage-structure have also been studied in many references, such as maturation delay ^[19,20,21], periodic functional responses ^[22,23], prey refuge ^[24], combination of prey refuge and additional food for predator ^[25], anti-predator ^[26,27,28], as well as diffusive effects ^[29].

Recently, the authors ^[30] have studied a Rosenzweig-MacArthur predator–prey model with stage–structure in prey:

where $x_1 \equiv x_1(t)$ and $x_2 \equiv x_2(t)$ are respectively the densities of immature and mature prey population at time $t$, and $y \equiv y(t)$ is the density of predator population at time $t$. The following assumptions are taken in deriving model (1.1) :

$\rm(H1)$ The immature prey do not produce offspring and their growth rate depends entirely on the reproduction by the mature prey, where the growth is assumed to be logistically with constant intrinsic rate $r > 0$ and constant carrying capacity $K > 0$.

$\rm(H2)$ The immature prey is more vulnerable so predators only consume immature prey. The predation mechanism is assumed to follow the Holling Type Ⅱ functional response where the maximum predation rate is $v_0 > 0$ and the environmental protection for immature prey is denoted by $n_1 > 0$.

$\rm(H3)$ The immature prey grow up and turn into mature prey with a conversion rate $\alpha > 0$.

$\rm(H4)$ The growth of mature prey population depends only on the conversion of immature prey into mature prey. The mature prey do not have a risk to be attacked by predator. The death rate of the mature prey population is $\delta_0 > 0$.

$\rm(H5)$ The predator consumes only the immature prey and the growth rate of predator population is proportional to the predation rate, where the conversion rate of consumed prey into predator births is $c > 0$. The death rate of the predator population is $\delta_1 > 0$.

Notice that the growth of immature prey in model (1.1) is inhibited by the intraspecific competition between the mature and immature prey. In nature, the intraspecific competition also occur between mature prey. Such intraspecific competition will be incoporated in this paper and hence we propose the following stage-structure predator-prey model

System (1.2) describes a predator-prey interaction where the predator only consumes immature prey. One ecological example of such a predator-prey system is the interaction between Moluccan megapode (Eulipoa wallacei) and rats Rattus sp. Rattus sp. is reported to prey upon Eulipoa eggs and chicks, and it cannot attack adult Eulipoa, see ^[31]. Another example of predator-prey system which can be described by model (1.2) is Chinese fire-bellied newt, which is unable to feed on the mature Rana chensinensis, can only prey on its immature, see ^[32]. For the sake of convenience, we describe all variables and parameters in system (1.2) in Table 1.

We remark that the intraspecific competition between the mature and immature prey is assumed to have the same strength as that between mature prey. To simplify our analysis, we reduce the number of parameters in system (1.2) by considering the following non-dimensional model

where $\tau = rt$, $u_1 = \frac{x_1}{K}$, $u_2 = \frac{x_2}{K}$, $u_3 = \frac{v_0y}{rK}$, $\omega = \frac{n_1}{K}$, $\alpha_1 = \frac{\alpha}{r}$, $\alpha_2 = \frac{cv_0}{r}$, $\beta_1 = \frac{\delta_0}{r}$ and $\beta_2 = \frac{\delta_1}{r}$.

As far as we are aware, the dynamics of system (1.3) have not been studied. Hence, a Rosenzweig-MacArthur predator-prey model with stage-structure in prey is introduced in this paper and the dynamical properties of the model are investigated. The organization of this paper is as follows. In Section 2, we analyze the boundedness, existence and uniqueness of the solution of the system (1.3). In Section 3 we determine all possible equilibrium points including the conditions for their existence and stability properties. The existence of Hopf-bifurcation and numerical simulations to illustrate our analytical finding are presented in Section 4 and Section 5, respectively. The conclusion of our study is given in the last section.

2.   Basic properties of the model

In this section we present some basic properties of the model (1.3) which include the existence and uniqueness, as well as the positivity and the boundedness of solutions of model (1.3). These are necessary because we deal with a population-related model. Therefore solutions of the model (1.3) must be non-negative. Boundedness of solutions explains that there are natural limits to population growth because resources are limited.

2.1.   Existence and uniqueness of solutions

The existence and uniqueness of solutions of system (1.3) will be investigated in the region $[0, \infty)\times\Omega _{M}$ where

for sufficiently large $M$. The existence of $M$ is assured by the boundedness of the solutions, which will be discussed later. Let ${U} = (u_1, u_2, u_3)^T$, $\bar{U} = (\bar{u}_1, \bar{u}_2, \bar{u}_3)^T$ and consider a mapping

where

For any $U, \bar{U} \in\Omega_M$, we can show that

where $L_1 = (2\alpha_1+M+(1+\alpha_2)M\omega)$, $L_2 = (1+\beta_1+3M)$, $L_3 = (\beta_2+(1+\alpha_2)(\omega+M)M)$, and $L =$max$\{L_1, L_2, L_3\}$. Clearly that $H(U)$ satisfies the Lipschitz condition with respect to $U$, and therefore system (1.3) with any positive initial condition $u_1(0) \ge 0, u_2(0) \ge 0, u_3(0) \ge 0$ has a unique solution ${U(\tau)} = (u_1(\tau), u_2(\tau), u_3(\tau))^T \in\Omega_{M}$. Hence, we have the following theorem related to the existence and uniqueness of the system (1.3).

Theorem 2.1. The stage-structure predator-prey system ($1.3$) subject to any non-negative initial value $\left(u_1(0), u_2(0), u_3(0)\right)$ has a unique solution $\left(u_1(\tau), u_2(\tau), u_3(\tau)\right)^T \in\Omega_{M}$ for all $\tau > 0$.

2.2.   Positive invariance

To show the positive invariance of the system (1.3), we first write the system (1.3) and its initial conditions in matrix form as

where $U(\tau) = (u_1(\tau), u_2(\tau), u_3(\tau))^T, U(0) = (u_1(0), u_2(0), u_3(0))^T \in \Omega_M$ and $H(U)$ defined as in (2.1). From the previous subsection it was shown that $H(U)$ is locally Lipschitz. Furthermore, it is obvious to show that $H_i(U_i) |_{u_i = 0} \ge 0$, for $i = 1, 2, 3$. Then, according to Nagumo theorem, the solution $U(\tau)$ of system (2.4) with initial condition $U(0) = U_0 \in \Omega_M$ remains in $\Omega_M$ for any $\tau \ge 0$. Hence, we have the following theorem.

Theorem 2.2. Every solution of (1.3) with initial conditions $u_1(0) \ge 0, u_2(0) \ge 0, u_3(0) \ge 0$ which exists in $\Omega_M$, remains positive for all $\tau > 0$.

2.3.   Boundedness of solutions

Theorem 2.3. All solutions of system (1.3) which start $\in\Omega_M$ are uniformly bounded.

Proof. We begin by assuming that $u_1(\tau)+u_2(\tau) \ge 1$ for all $\tau \ge 0$. Using the first two equations in system (1.3), we have

It is clear that $\frac{d}{d \tau} (u_1 + u_2) \le 0$, and thus $u_1 + u_2$ is monotone decreasing. Let's denote

Based on Barbalat lemma, we can show that if $\kappa > 1$ then

This contradiction leads to $\kappa = 1$. Since $u_1 + u_2$ is monotone decreasing function, we have that

meaning that there exists $\tau_0 > 0$ such that for $\tau > \tau_0$ we get $u_1 + u_2 \le 1 +\epsilon$ or $u_1 \le 1 + \epsilon - u_2$ for any $\epsilon \ge 0$. Then, from the second equation in system (1.3) we get for any $\tau \ge \tau_0$ that

It is obvious that

By taking $\epsilon = 0$, we get

Similarly, we can also show that

Hence $u_1$ and $u_2$ is bounded above. The latter inequality shows that there exists $\tau_1$ such that $u_1 \le \frac{\beta_1}{\alpha_1 + \beta_1} + \epsilon$ for all $\tau \ge \tau_1$. Now define $u = u_1 + \frac{u_3}{\alpha_2}$. Using the first and the third equations in system (1.3), we obtain

By taking $\epsilon = 0$ and denoting $\eta = |\beta_2 - \alpha_1| \frac{\beta_1}{\alpha_1 + \beta_1}$, we get the following inequality

from which we have

Hence, $u_3$ is also bounded above because

Now we suppose that our assumption that $u_1 + u_2 \ge 1$ is violated, i.e. there exists $\tau_2 > 0$ such that we have for the first time $u_1(\tau_2)+u_2(\tau_2) = 1$. Then, we easily have

Consequently, whenever once a solution with $u_1 + u_2$ has entered the interval $(0, 1)$ then it remains bounded there for all $\tau > \tau_2$, i.e.

From the second equation in system (1.3) we get

Using the previous argument, we can prove that

Similarly, we can also show that

Again, by defining $u = u_1 + \frac{u_3}{\alpha_2}$ and using the previous argument, we can easily show that $u$ is bounded above and therefore $u_3$ is also bounded.

At last, if $u_1(0)+u_2(0) < 1$, then by using the previous arguments, we have that $u_1(\tau) + u_2(\tau) < 1$ and $u = u_1 + \frac{u_3}{\alpha_2}$ is bounded for all $\tau > 0$. This completes the proof.

3.   Existence and stability analysis of equilibrium points

It is easy to show that system (1.3) has three non-negative equilibrium points as follows.

● The extinction equilibrium $E_0 = (0, 0, 0)$, which there is no population in the habitat.

● The predator-free equilibrium ${E_1} = (\hat {u}_1, \hat {u}_2, 0)$, which exists if

where $\hat {u}_1 = \beta_1(1-\beta_1) / (\alpha_1+\beta_1)$, $\hat {u}_2 = \alpha_1\hat {u}_1 / \beta_1$.

● The interior equilibrium $E^* = ({u_1^*}, {u_2^*}, {u_3^*})$, i.e. all of species coexist, where

The interior or coexistence equilibrium point exists if

where $\widetilde{\omega} = \frac{\beta_1(\alpha_2-\beta_2)(1-\beta_1)}{\beta_2 (\alpha_1+\beta_1)}$.

From the description above, it can be seen that if the conditions (3.2) are fulfilled then all three equilibrium points exist. The local stability of each equilibrium points of system (1.3) is shown in the following theorem.

Theorem 3.1. For system (1.3), we have the following stability properties of its equilibrium points:

(i) The equilibrium point $E_0$ is locally asymptotically stable if $\beta_1 > 1$.

(ii) The equilibrium point ${E_1}$ is locally asymptotically stable if $\omega > \widetilde{\omega}$.

(iii) The coexistence equilibrium point $E^*$ is locally asymptotically stable if $\varphi_1 > 0, \varphi_3 > 0$ and $\varphi_1\varphi_2 > \varphi_3$ where $\varphi_1, \varphi_2$ and $\varphi_3$ are defined as in the proof.

Proof. The local stability of all equilibrium points can be studied from the linearization of the system (1.3). The Jacobian matrix of the system (1.3) at a point $(u_1, u_2, u_2)$ is given by

By observing the eigenvalues of the Jacobian matrix (3.3) at each equilibrium point, we have the following stability properties.

(ⅰ) The Jacobian matrix of the system (1.3) at $E_0$ has eigenvalues $\lambda_1 = -\beta_2$ and $\lambda_{2, 3} = - \frac{1}{2} B_1 \pm \frac{1}{2}\sqrt{B_1^2-4C_1}$, where $B_1 = \alpha_1+\beta_1 > 0$ and $C_1 = \alpha_1(\beta_1-1)$. If

then $C_1 > 0$. Consequently $\operatorname {Re}(\lambda_{2, 3}) < 0$ and $E_0$ is locally asymptotically stable.

(ⅱ) The Jacobian matrix of the system (1.3) at ${E_1}$ has eigenvalues $\lambda_1 = \frac{\beta_1(\alpha_2-\beta_2)(1-\beta_1)-\beta_2\omega(\alpha_1+\beta_1)}{\omega(\alpha_1+\beta_1)+\beta_1(1-\beta_1)}$ and $\lambda_{2, 3} = - \frac{1}{2} B_2 \pm \frac{1}{2}\sqrt{B_2^2-4C_2}$, where $B_2 = \alpha_1+\beta_1 + \alpha_1 {u}_1 / \beta_1 > 0$ and $C_2 = \alpha_1(1-\beta_1)$. If $E_1$ exists then $C_2 > 0$ and therefore $\operatorname {Re}(\lambda_{2, 3}) < 0$. Thus, ${E_1}$ is locally asymptotically stable if $\lambda_1 < 0$, that is when

Notice that if $E_1$ exists and is locally asymptotically stable, then $E_0$ is unstable and $E^*$ does not exist.

(ⅲ) The characteristic equation of the Jacobian matrix of the system (1.3) at $E^*$ is given by the following cubic equation

where

The stability of $E^*$ can be determined by the Routh-Hurwitz criterion ^[33], i.e. $E^*$ is locally asymptotically stable if $\varphi_i > 0$, $i = 1, 3$ and

From conditions (3.1) and (3.4), we note that if the intrinsic growth rate of immature prey is less than the death rate of mature prey then the extinction point $E_0$ will be locally asymptotically stable. Otherwise, the predator-free point $E_1$ exists whenever the death rate of mature prey is less than the intrinsic growth rate of immature prey. Thus, the existence of predator-free point $E_1$ and the stability of the extinction point $E_0$ are dependent on both intrinsic growth rate of immature prey and the death rate of mature prey. Next, condition (3.2) says that the coexistence equilibrium $E^*$ exists if the death rate of the predator is less than the predation rate.

Based on the existence and stability results of equilibrium points of the system (1.3), we remark that when $\omega > \widetilde{\omega}$ then the coexistence equilibrium point $(E^*)$ does not exist. In this case, parameter $\beta_1$ is a threshold: if it is greater than unity then the asymptotically stable extinction equilibrium point $(E_0)$ is the unique equilibrium point; if it is less than unity then $E_0$ becomes unstable and there appears an asymptotically stable equilibrium point $(E_1)$. Hence the system (1.3) undergoes a forward bifurcation at $\beta_1 = 1$. In fact a forward bifurcation may also be observed when $\beta_1 < 1$, where in this case $\omega$ behaves as the threshold parameter. Since $\beta_1 < 1$, the extinction point $(E_0)$ is unstable. Furthermore, when $\omega > \widetilde{\omega}$, equilibrium point $E_1$ is asymptotically stable and $E^*$ does not appear. On the contrary, if $\omega < \widetilde{\omega}$ then $E_1$ becomes unstable and $E^*$ exists if $\alpha_2 > \beta_2$. If the Routh-Hurwitz criterion for the characteristics equation (3.6) is fulfilled then $E^*$ is asymptotically stable. Consequently, we have the following corollary about the occurrence of bifurcation in the system (1.3).

Corollary 3.2. (ⅰ) If $\omega > \widetilde{\omega}$ then the system (1.3) has a forward bifurcation from $E_0$ to $E_1$ around $\beta = 1$.

(ⅱ) If $\beta_1 < 1$, $\alpha_2 > \beta_2$, $\varphi_1 > 0, \varphi_3 > 0$ and $\varphi_1\varphi_2 > \varphi_3$ then the system (1.3) undergoes a forward bifurcation at $\omega = \widetilde{\omega}$ from $E_1$ to $E^*$.

4.   Existence of Hopf-bifurcation

In this section, we study the Hopf-bifurcation around the coexistence equilibrium point $E^* = ({u_1^*}, {u_2^*}, {u_3^*})$ of the system (1.3). We consider $\omega$ and $\alpha_2$ as the bifurcation parameters. $\omega = n_1/K$ and $\alpha_2 = c v_0 / r$ are chosen as the bifurcation parameters because $r$ and $K$ are strongly related to the growth of immature prey, which controls energy input in the predator-prey system. Furthermore $a, v_0$ and $n_1$ are important parameters governing the exchange of energy from prey to predator.

Theorem 4.1. The stage-structure predator-prey system (1.3) undergoes Hopf-bifurcation around the coexistence equilibrium $E^*$ when $\omega$ passes through $\omega^*$ where $\omega^*$ satisfies $\varphi(\omega^*) = \varphi_1(\omega^*) \varphi_2(\omega^*) - \varphi_3(\omega^*) = 0$, $\varphi_i(\omega^*) > 0$ for $i = 1, 2$ and $\varphi_1(\omega^*)\varphi_2'(\omega^*)+\varphi_2(\omega^*)\varphi_1'(\omega^*)-\varphi_3'(\omega^*) \ne 0$.

Proof. At $\omega = \omega^*$, by the condition $\varphi(\omega^*) = 0$, the characteristic equation (3.6) can be written as

If $\varphi_1 = m_1+m_2-m_3 > 0$ and $\varphi_2 = m_1m_2+m_4+m_5-m_3m_6 > 0$, then the roots of the equation (4.1) are $\lambda_1 = -\varphi_1 < 0$ and $\lambda_{2, 3} = i\sqrt{\varphi_2}$. For any $\omega$, the characteristic roots can generally be written as $\lambda_1(\omega) = -\varphi_1(\omega)$, and $\lambda_{2, 3} = \mu_1(\omega)\pm{i}\sigma_1(\omega)$. By substituting $\lambda(\omega) = \mu_1(\omega) + i\sigma_1(\omega)$ into equation (4.1), we have

where

It is obvious that equation (4.2) is satisfied by $\zeta_1(\omega) = 0$ and $\zeta_2(\omega) = 0$. Furthermore, by taking the derivate of equation (4.2) with respect to $\omega$, we get

where

From system (4.3), we can show that

We notice that $\mu_1(\omega^*) = 0$, $\sigma_1(\omega^*) = \sqrt{\varphi_2(\omega^*)}$ and

Therefore the following equation is obtained

Thus, if $\varphi_1(\omega^*)\varphi_2'(\omega^*)+\varphi_2(\omega^*)\varphi_1'(\omega^*)-\varphi_3'(\omega^*) \ne 0$, then transvesality condition is satisfied, and Hopf-bifurcation occurs when $\omega$ passes through $\omega = \omega^*$.

According to Theorem 4.1, there exists a Hopf bifurcation in the stage-structure predator-prey model (1.3) where the Hopf bifurcation is controlled by $\omega$. In fact, using the same argument as in the proof of Theorem 4.1, we can show that the Hopf bifurcation can also be controlled by parameter $\alpha_2$. The possibility of the Hopf bifurcation occurance is stated in the following theorem.

Theorem 4.2. The stage-structure predator-prey system (1.3) undergoes Hopf-bifurcation around the coexistence equilibrium $E^*$ when $\alpha_2$ passes through $\alpha_2^*$ where $\alpha_2^*$ satisfies $\varphi(\alpha_2^*) = \varphi_1(\alpha_2^*) \varphi_2(\alpha_2^*) - \varphi_3(\alpha_2^*) = 0$ provided that $\varphi_i(\alpha_2^*) > 0$ for $i = 1, 2$ and $\varphi_1(\alpha_2^*)\varphi_2'(\alpha_2^*)+\varphi_2(\alpha_2^*)\varphi_1'(\alpha_2^*)-\varphi_3'(\alpha_2^*) \ne 0$.

5.   Numerical simulations

For graphical confirmation of the previous analytical results, several numerical simulations have been carried out. The model (1.3) is solved using the fourth-order Runge-Kutta method with some different initial conditions and some hypothetical values of the parameters. We first consider the following parameter values: $\alpha_1 = 0.5, \alpha_2 = \beta_2 = 0.3$, $\omega = 0.03$ and $\beta_1 = 1.01$. Because $\beta_1 > 1$, Theorem 3.1 says that the extinction of population point ($E_0$) is the only equilibrium point which is locally asymptotically stable. This can be understood from the fact that in this case $\beta_1 = \frac{\delta_0}{r} > 1$, i.e. the death rate of mature prey is larger than the intrinsic growth rate of immature prey. Hence the system will asymptotically be convergent to the extinction equilibrium point. The behavior of this case is depicted in Figure 1.(a). Next, we carry out a simulation using the same previous parameter values, except $\beta_1 = 0.3$. In addition to the unstable extinction of population point ($E_0$), system (1.3) also has the predator-free equilibrium $E_1 = (0.2625, 0.4375, 0.0)$. Since $\alpha_2 = \beta_2$, we have $\widetilde{\omega} = 0$ and therefore $\omega > \widetilde{\omega}$. Based on Theorem 3.1, $E_1$ is asymptotically stable. The case of $\alpha_2 = \beta_2$ is equivalent to the case of $c v_0 = \delta_1$. Hence, from the third equation in the system (1.2), we have that the per capita growth rate of predator is less than the predator death rate. This explains the extinction of predator population. The numerical solutions which describe this situation is plotted in Figure 1.(b). To see a complete view of the dynamics of the system (1.3) with parameter values $\alpha_1 = 0.5, \alpha_2 = \beta_2 = 0.3$, and $\omega = 0.03$, in Figure 2 we plot the equilibrium densities of immature and mature prey as function of $\beta_1$. The equilibrium density of predator in this case is zero and thus it is not shown in the picture. Figure 1 shows that if $\beta_1 < 1$ then $E_0$ is unstable while $E_1$ is stable. On the other hand, when $\beta_1 > 1$, $E_1$ disappears and $E_0$ becomes stable. Thus, the system (1.3) undergoes a forward bifurcation.

We now choose $\alpha_1 = \alpha_2 = 0.5$ and $\beta_1 = \beta_2 = 0.3$. Using Theorem 4.1, we can check that the system undergoes a Hopf bifurcation around $E^*$ where the bifurcation point is at $\omega^* = 0.03989$. To verifiy the existence of Hopf bifurcation numerically, we take two different values of $\omega$, i.e. $\omega = 0.03 < \omega^*$ and $\omega = 0.04 > \omega^*$. For $\omega = 0.03$, we find that the system (1.3) has coexistence equilibrium point $E^* = (0.045, 0.075, 0.072)$, with $\varphi_1 = 1.2617 > 0$ and $\varphi_3 = 0.0338 > 0$, but $\varphi_1 \varphi_2 - \varphi_3 = -0.0323 < 0$. Since the Routh-Hurwitz criterion is not satisfied, the coexistence equilibrium point is unstable. However, for $\omega = 0.04$, the coexistence equilibrium point $E^* = (0.06, 0.1, 0.09)$ is stable because the Routh-Hurwitz criterion is satisfied, i.e. $\varphi_1 = 1.26 > 0$, $\varphi_3 = 0.0324 > 0$ and $\varphi_1 \varphi_2 - \varphi_3 = 0.00036 > 0$. The numerical solutions of the system (1.3) for $\omega = 0.03$ and $\omega = 0.04$ are shown in Figure 3. It can be seen in this figure that our numerical solutions are convergent to a periodic solution when $\omega = 0.03 < \omega^*$ and convergent to $E^*$ when $\omega = 0.04 > \omega^*$. This behavior indicates that the coexistence equilibrium point $(E^*)$ changes its stability, i.e. $E^*$ is unstable and convergent to a limit cycle when $\omega < \omega^*$ and it becomes asymptotically stable if $\omega > \omega^*$. In other words, the system (1.3) experiences a Hopf bifurcation driven by parameter $\omega$. Figure 3 gives an indication that the Hopf bifurcation is supercritical because the limit cycle is stable,

As stated in Theorem 4.2, the Hopf bifurcation can also be controlled by parameter $\alpha_2$. To show this behavior we perform simulation using $\alpha_1 = 0.5, \beta_1 = \beta_2 = 0.3$ and $\omega = 0.03$. The critical $\alpha_2^*$ for the occurance of Hopf bifurcation can be determined using Theorem 4.2, where in this case we get $\alpha_2^* = 0.5822$. Using Theorem 3.1, it can be verified that the coexistence equilibrium point $(E^*)$ is unstable for the case of $\alpha = 0.55$, and it is asymptotically stable for $\alpha = 0.59$. As depicted in Figure 4, such stability properties coincide with our numerical simulations. It is clearly seen in Figure 4.(a) that for $\alpha_2 = 0.55 < \alpha_2^*$, $E^*$ is unstable and the solutions converge to a limit cycle (periodic solution). If we take $\alpha_2 = 0.59 > \alpha_2^*$, then the solution is convergent to $E^*$, see Figure 4.(b). From Figure 4, we also see the occurance of (supercritical) Hopf bifurcation controlled by $\alpha_2$.

To see the detail dynamics of the system (1.3) with parameter $\alpha_1 = 0.5, \; \beta_1 = 0.3$ and $\beta_2 = 0.3$, we plot the stability region in $(\omega, \alpha_2)-$plane, see Figure 5. In this picture we see that there are three different regions: the cyan region represents the stability area of $E_1$ ($E_0$ exists but it is unstable, $E^*$ does not exist); the green region represents the stability area of $E^*$ ($E_0$ and $E_1$ exists but they are not stable); and the yellow region corresponds to the stable limit cycle (all the three equilibrium points exists but they are unstable). We identify the occurrence of a Hopf bifurcation as a change from the green region to the yellow region due to the changes of parameters $\omega$ or $\alpha_2$. We also notice the occurrence of a forward bifurcation caused by the changes of $\alpha_2$, see the changes from cyan region to green region. To provide a clearer illustration, we plot a bifurcation diagram for system (1.3) with parameter $\alpha_1 = 0.5, \; \beta_1 = \beta_2 = 0.3$ and $\omega = 0.03$, see Figure 6. From this figure we can see the appearance of transcritical bifurcation which is caused by the changing of the value of $\alpha_2$. When $\alpha_2 < \alpha_2^1 = 0.3343$, the equilibrium point $E_1$ is asymptotically stable. On the contrary, if $\alpha_2^1 < \alpha_2 < \alpha_2^2 = 0.3922$ then equilibrium point $E_1$ becomes ustable and equilibrium point $E^*$ is asymtotically stable. In Figure 6 we also observe that the system (1.3) undergoes 2-Hopf bifurcation. In this case, the equilibrium point $E^*$ is asymptotically stable for $\alpha_2^1 < \alpha_2 < \alpha_2^2 = 0.3922$ or $\alpha_2 > \alpha_2^* = 0.5822$, and it is unstable for $\alpha_2^2 < \alpha_2 < \alpha_2^*$. In the latter case, the system (1.3) is convergent to a periodic solution or limit cycle.

6.   Conclusion

A Rosenzweig-MacArthur predator-prey model with stage-structure in prey has been discussed. It was shown that the proposed model has three equilibrium points, i.e. the extinction of population $(E_0)$, the predator-free point $(E_1)$, and the coexistence point $(E^*)$. All of these equilibrium points are conditionally asymptotically stable. Our analysis also showed that the proposed model exhibits a Hopf-bifurcation which can be driven by $(\omega)$ or $(\alpha_2)$. Furthermore, the proposed model may also undergo a forward bifurcation for suitable parameter values. The analytical findings have been confirmed by our numerical simulations.

Acknowledgments

This research was funded by FMIPA via PNBP-University of Brawijaya according to DIPA-UB No. DIPA-023.17.2.677512/2020, under contract No. 12/UN10.F09/PN/2020.

Conflict of interest

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