Generalized polynomial exponential sums and their fourth power mean

1.   Introduction

Schistosomiasis is a neglected tropical disease caused by parasitic helminths of the genus Schistosoma ^[1]. Approximately 250 million people are infected with schistosomiasis ^[2]. Schistosomiasis is commonly found in parts of the Middle East, South America, Southeast Asia and Africa ^[1]. This disease is mostly found in low-income neighborhoods that lack access to clean water and adequate sanitation ^[3]. The three main species infecting humans are Schistosoma haematobium, S. japonicum and S. mansoni ^[4]. S. japonicum causes intestinal schistosomiasis and hepatosplenic schistosomiasis in China, The Philippines and Indonesia, while urogenital schistosomiasis in Africa and parts of the Middle East is caused by S. haematobium infection. In Africa, the Arabian peninsula and Latin America, S. mansoni causes intestinal and liver disease ^[2]. Colley et al. ^[5] stated that S. japonicum is a zoonotic species, while S. mansoni can infect baboons and rodents. On the other hand, S. haematobium is declared to be not a zoonotic species because its definitive host is exclusively humans.

As a parasite, the Schistosoma worm has an elaborate life cycle. To accomplish its life cycle, Schistosoma requires two types of host, namely, intermediate hosts and definitive hosts. Some species of snails can act as intermediate hosts for the parasite. S. japonicum used the Oncomelania snail as an intermediate host while the Bulinus snail could serve as an intermediate host for S. haematobium. Meanwhile, S. mansoni can infect and make Biomphalaria snails intermediate hosts ^[2]. Humans can serve as a definitive host for the parasite ^[6]. In the human body, parasites reproduce sexually. The transmission cycle process occurs when the eggs of adult Schistosoma worm pairs are excreted through the feces and urine by infected humans. After some time, the eggs that survive in the wild hatch and release miracidia which can infect suitable snails. The parasites reproduce asexually in the snail body ^[7]. Inside the snail's body, miracidia develop and then produce cercariae which will be released from the snail's body. Cercariae that successfully infect humans can develop into adult worms ^[2].

Some researchers stated that molluscicide can be used to manage snail populations and schistosomiasis spread ^[8,9,10]. This intervention is very effective in decreasing the prevalence of schistosomiasis. However, it is believed that molluscicide may damage the ecosystem ^[11]. Some researchers recommended the use of snail predators or snail competitors as biological control agents of the snail population ^[12,13,14]. Secor ^[15] stated that several types of fish are predators of host snails. An example is the Cichlid fish in Senegal. In addition, Sokolow et al. ^[13] reported that prawns in Senegal are predators of host snails, while Mkoji et al. ^[16] found that the introduction of Red Swamp Crayfish in schistosomiasis endemic areas in Kenya could reduce snail populations and schistosomiasis cases. Different interventions that can be used to control a snail population are snail habitat modification, such as turning farmland to fish ponds, digging and dredging canals, planting of proper trees such as cotton in high-altitude lands or poplar in marshlands. These interventions can reduce snail habitats, which results in a decrease in snail density and survival ^[17]. Meanwhile, the elimination of snails in some areas in China is achieved after the construction of lake embankments. The construction of the lake embankment prevents flooding which results in a reduced spread of snails around the lake ^[18]. This suggests that snail habitat modification is related to the environmental carrying capacity of the snail population.

To study the dynamics of schistosomiasis spread, Macdonald ^[19] proposed a mathematical model. It is well recognized that the model is the first mathematical model of schistosomiasis. Since then, many researchers have proposed mathematical models which are related to schistosomiasis. Schistosomiasis models that take into account the parasite density in the environment are discussed in ^{[20,21,22,23]}. In these works, the parasite is divided into two distinct classes, i.e., miracidia and cercariae. Miracidia and cercariae are the parasite stages that can infect snails and humans, respectively. To control vector-borne disease spread, we can use any biological control agents, e.g., predator, competitor or parasites of the vector ^[24,25]. A schistosomiasis model considering competitor resistant snails as biological control agents is discussed by Diaby et al. ^[26]. They found that the introduction of competing snails that are resistant to Schistosoma infection could reduce the population of snails that could serve as intermediate hosts for Schistosoma. In 2021, Nur et al. ^[27] extended the model proposed in ^[23] by adding a snail predator population as a biological control agent. Their results showed that the introduction of snail predators in the host snail habitat could reduce the host snail population, cases of schistosomiasis in humans and snails and the number of parasites in the environment. In this paper, we extend the model developed in ^[27] by adding treated human compartments. In addition, to study the impact of snail habitat modification on the spread of schistosomiasis, we assume that the snail population grows logistically. The rest of the paper is organized as follows. Section 2 describes the model formulation and some of the basic properties of the model. The stability analysis of the equilibrium points is discussed in Section 3. In Section 4, several numerical simulation results are presented. The paper ends with some conclusions in Section 5.

2.   Model formulation and basic properties of the model

2.1.   Model formulation

In this subsection, we describe the basis and assumptions used in model formulation. We only consider schistosomiasis caused by S. haematobium. Therefore, in our model, only humans are the definitive hosts, because S. haematobium is not zoonotic. We construct the schistosomiasis model based on the life cycle of the Schistosoma worm ^[6,7,28]. There is a latent period because the parasites need time to migrate through the body, mature and pair up to begin producing eggs. Further, we constructed the model to study the impact of early treatment of exposed humans (humans who have been infected but are still in the prepatent period). Consequently, it is essential to include the exposed human compartment in the model constructed. Thus, we divide the human population into susceptible humans $(S_{h})$, exposed humans $(E_{h})$, infectious humans $(I_{h})$ and treated humans $(T_{h})$. People who have been treated can be reinfected; therefore, in Figure 1 which shows the compartment diagram, there is a flow from the treated human to the susceptible human. Anderson et al. ^[29] stated that the latent period of the snail intermediate host is very dependent on the ambient temperature. High temperatures cause a short latent period. In this article, we assume that the environment has a high temperature so that the latent period in snails is short enough to be negligible. Thus, the snail population is only divided into susceptible snails $(S_{v})$ and infectious snails $(I_{v})$. Similar to ^{[20,21,22,23]}, we only consider two stages of Schistosoma worm development, i.e., miracidia and cercariae. Therefore, the parasite is divided into miracidia $(M)$ and cercariae $(C)$. According to the epidemiology of schistosomiasis ^[7,28], miracidia can infect snails and cercariae have the ability to infect humans. Hence, susceptible humans and susceptible snails may get infected after adequate contact with cercariae and miracidia, respectively. We assume that snails that can serve as intermediate host is the only food for the predator. It should be noted that S. haematobium in the infectious human can produce 100–200 per day per pair. Some of which will be excreted from the infected human body ^[30]. Moreover, an infectious snail can release about 200 cercariae into the environment for the cases of S. haematobium infection ^[7]. Therefore, the decrease of the parasites in the environment on account of direct interaction with susceptible humans and susceptible snails is neglected. In addition, we assume that there is no recovery for infectious snail and no disease-induced death. Figure 1 shows the transition and interaction between compartments. Based on the assumptions and Figure 1, we get the schistosomiasis model as follows:

where

$\phi_{e}, \phi \in [0, 1)$ and the other parameters are positive. The definition of all parameters is given in Table 1.

2.2.   Non-negativity and Boundedness

Theorem 2.1. Solutions of system (2.1) with non-negative initial conditions are non-negative.

Proof. To prove this theorem, we use the method used in ^[31,32]. From system (2.1), we have

where $\zeta(x) = \left\{x = 0\;{\rm{and}}\;S_{h}, E_{h}, I_{h}, T_{h}, S_{v}, I_{v}, C, M, P\in C\left(\mathbb{R}_{+0}, \mathbb{R}_{+0}\right)\right\}$. Based on Lemma 2 presented in ^[33], the invariant region of system (2.1) is $\mathbb{R}_{+0}^{9}$. Hence, solutions of system (2.1) with non-negative initial conditions are non-negative.

Theorem 2.2. Solutions of system (2.1) with non-negative initial values are bounded.

Proof. Here, $N_{h}$ and $N_{v}$ are the total number of humans and the total number of snails, respectively. It is clear that $N_{h} = S_{h}+E_{h}+I_{h}+T_{h}$ and $N_{v} = S_{v}+I_{v}$. From system (2.1), we get

Using the technique of an integrating factor, we obtain the solution of the first equation of (2.2).

It is clear that $0\leq N_{h}(t)\leq \frac{\Omega_{h}}{d_{h}}$ for $t \geq 0$ whenever $0\leq N_{h}(0)\leq \frac{\Omega_{h}}{d_{h}}$. Thus, $N_{h}$ is bounded.

Now, let $\Pi = N_{v}+P$. We will show that $\Pi = N_{v}+P$ is bounded. From the second and third equations of (2.2), we have

where $\xi \geq \tau$, $q = \min \{w_{5}, d_{p}\}$ and $\bar{k} = \frac{\varphi - w_{5}}{\rho}$. According to Gronwall's Inequality Lemma, it is found that $\Pi = N_{v}+P\leq \frac{\varphi \bar{k}}{q}$. Hence, $N_{v}$ and $P$ are bounded. Now, we use the previous results to show that $C$ and $M$ are bounded. Taking into account that $I_{h}\leq N_{h}\leq \frac{\Omega_{h}}{d_{h}}$ and $I_{v}\leq P+N_{v}\leq \frac{\varphi \bar{k}}{q}$, we have

According to Gronwall's Inequality Lemma, we get $C\leq \frac{\sigma \varphi \bar{k}}{q d_{c}}$ and $M \leq \frac{\alpha h_{h}g_{h}\Omega_{h}}{d_{h}d_{m}}$. Thus, the solutions of system (2.1) are bounded.

Therefore, we have the following invariant region of system (2.1).

The solutions of system (2.1) with initial value in $\Xi_{+}$ are always in $\Xi_{+}$.

2.3.   Equilibrium points and basic reproduction number

System (2.1) has five equilibrium points. The first equilibrium point is the snail and predator extinction point ${\bf{Y}}_{0} = \left(S_{h}^{*0}, E_{h}^{*0}, I_{h}^{*0}, T_{h}^{*0}, S_{v}^{*0}, I_{v}^{*0}, C^{*0}, M^{*0}, P^{*0}\right) = \left(\frac{\Omega_{h}}{d_{h}}, 0, 0, 0, 0, 0, 0, 0, 0\right)$, which always exists in $\mathbb{R}^{9}_{+0}$. The second equilibrium point is predator extinction and disease-free point ${\bf{Y}}_{0}^{a} = (S_{h}^{a*}, E_{h}^{a*}, I_{h}^{a*}, T_{h}^{a*}, S_{v}^{a*}, I_{v}^{a*}, C^{a*}, M^{a*}, 0) = \left(\frac{\Omega_{h}}{d_{h}}, 0, 0, 0, \frac{\varphi-w_{5}}{\rho}, 0, 0, 0, 0\right)$, which exists in $\mathbb{R}_{+0}^{9}$ if $\frac{\varphi}{w_{5}} > 1$. The third equilibrium point is predator extinction-endemic point ${\bf{Y}}_{1}^{a} = (S_{h}^{a**}, E_{h}^{a**}, I_{h}^{a**}, T_{h}^{a**}, S_{v}^{a**}, I_{v}^{a**}, C^{a**}, M^{a**}, 0)$, where

It is easy to see that ${\bf{Y}}_{1}^{a}$ exists in $\mathbb{R}_{+0}^{9}$ if $R_{e}^{a} > 1$ and $\frac{\varphi}{w_{5}} > 1$. The fourth equilibrium point is disease-free point ${\bf{Y}}_{0}^{b} = (S_{h}^{b*}, E_{h}^{b*}, I_{h}^{b*}, T_{h}^{b*}, S_{v}^{b*}, I_{v}^{b*}, C^{b*}, M^{b*}, P^{b*}) = \left(\frac{\Omega_{h}}{d_{h}}, 0, 0, 0, \frac{d_{p}}{\tau}, 0, 0, 0, \frac{\tau(\varphi-w_{5})-\rho d_{p}}{\xi \tau}\right)$ which exists in $\mathbb{R}_{+0}^{9}$ if $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} > 1$. The last equilibrium point is endemic point ${\bf{Y}}_{1}^{b} = (S_{h}^{b**}, E_{h}^{b**}, I_{h}^{b**}, T_{h}^{b**}, S_{v}^{b**}, I_{v}^{b**}, C^{b**}, M^{b**}, P^{b*})$, where

It is clear that ${\bf{Y}}_{1}^{b}$ exists in $\mathbb{R}_{+0}^{9}$ if $R_{e}^{b} > 1$ and $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} > 1$.

In epidemiology, there is an important number recognized as the basic reproduction number. In this paper, we employ the next-generation matrix method described in ^[34] to investigate the basic reproduction number. First, we establish $\mathscr{F}$ and $\mathscr{V}$, explained in ^[34]. $\mathscr{F}$ relates to the new infection terms, while $\mathscr{V}$ represents the transition terms. We consider that the infected classes in our model are $E_{h}$, $I_{h}$, $T_{h}$, $I_{v}$, $C$ and $M$. Consequently, we obtain the following $\mathscr{F}$ and $\mathscr{V}$.

We now determine the Jacobian matrix of $\mathscr{F}$ and $\mathscr{V}$ at arbitrary disease-free equilibrium point $\left(S_{h}^{*}, 0, 0, 0, S_{v}^{*}, 0, 0, 0, P^{*}\right)$ and obtain the following results.

According to ^[34], the basic reproduction number is the dominant eigenvalue of $F_{0}V_{0}^{-1}$.

Obviously, the dominant eigenvalue of $F_{0}V_{0}^{-1}$ is given by

Substituting ${\bf{Y}}_{0}^{a}$ into (2.3) produces the basic reproduction number when the snail survives but the predator goes to extinction.

Substituting ${\bf{Y}}_{0}^{b}$ into (2.3) yields the basic reproduction number when the snail and predator survive.

We notice that $R_{e}^{a} = (R_{0}^{a})^{2}$ and $R_{e}^{b} = (R_{0}^{b})^{2}$. Since $R_{0}^{a}$ and $R_{0}^{b}$ are always positive, it is clear that ${\bf{Y}}_{1}^{a}$ exists in $\mathbb{R}_{+0}^{9}$ if $R_{0}^{a} > 1$ and $\frac{\varphi}{w_{5}} > 1$. Furthermore, ${\bf{Y}}_{1}^{b}$ exists in $\mathbb{R}_{+0}^{9}$ if $R_{0}^{b} > 1$ and $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} > 1$.

3.   Stability analysis

3.1.   Local stability

The general Jacobian matrix of system (2.1) at arbitrary equilibrium point $\mathscr{E} = \left(S_{h}^{*}, E_{h}^{*}, I_{h}^{*}, T_{h}^{*}, S_{v}^{*}, I_{v}^{*}, C^{*}, M^{*}, P^{*}\right)$ is given by

where

First, we present the local stability condition for ${\bf{Y}}_{0}$.

Theorem 3.1. If $\frac{\varphi}{w_{5}} < 1$, then ${\bf{Y}}_{0}$ is locally asymptotically stable. If $\frac{\varphi}{w_{5}} > 1$, then ${\bf{Y}}_{0}$ is unstable.

Proof. The Jacobian matrix of system (2.1) at ${\bf{Y}}_{0}$ is

The eigenvalues of matrix (3.1) are the roots of polynomial (3.2).

It easy to see that all eigenvalues of matrix (3.1) are negative if $\frac{\varphi}{w_{5}} < 1$. Moreover, if $\frac{\varphi}{w_{5}} > 1$, then one eigenvalue of matrix (3.1) is positive. The proof is completed.

Now, we investigate the local stability condition of ${\bf{Y}}_{0}^{a}$. The Jacobian matrix of system (2.1) at ${\bf{Y}}_{0}^{a}$ is given as

Therefore, we get the following characteristic polynomial

where ${\bf{L}}_{1}(\lambda) = \lambda^{5}+l_{1}\lambda^{4}+l_{2}\lambda^{3}+l_{3}\lambda^{2}+l_{4}\lambda+l_{5}$ and

Here $K_{1} = w_{2}, K_{2} = w_{3}, K_{3} = w_{5}+\rho S_{v}^{a*}, K_{4} = d_{c},$ and $K_{5} = d_{m}$. It is clear that matrix (3.3) has three negative eigenvalues, i.e., $\lambda_{1} = -d_{h}$, $\lambda_{2} = -w_{3}$, and $\lambda_{3} = -(\varphi-w_{5})$. Moreover, $\lambda_{4} = \tau S_{v}^{a*}-d_{p} < 0$ if $\frac{\tau S_{v}^{a*}}{d_{p}} = \frac{\tau(\varphi-w_{5})}{\rho d_{p}} < 1$. The other eigenvalues of matrix (3.3) are zeros of ${\bf{L}}_{1}(\lambda)$. It is easy to see that $l_{i} > 0$ for $i = 1, 2, 3, 4$. Further, if $R_{0}^{a} < 1$, then $l_{5} > 0$. Notice that $l_{5} < 0$ if $R_{0}^{a} > 1$. Based on Descartes' rule of signs ^[35], matrix (3.3) has exactly one positive eigenvalue if $R_{0}^{a} > 1$. Thus ${\bf{Y}}_{0}^{a}$ is unstable if $R_{0}^{a} > 1$. We observe that $l_{5} = 0$ if $R_{0}^{a} = 1$. It indicates that one root of ${\bf{L}}_{1}(\lambda)$ is zero if $R_{0}^{a} = 1$. As a result, if $\frac{\tau S_{v}^{a*}}{d_{p}} = \frac{\tau(\varphi-w_{5})}{\rho d_{p}} < 1$ and $R_{0}^{a} = 1$ then one eigenvalue of matrix (3.3) is zero. We now use Routh-Hurwitz test ^[36] to investigate the local stability condition of ${\bf{Y}}_{0}^{a}$. Firstly we establish Routh's array.

Obviously, $l_{1} > 0$ always holds. According to the Routh-Hurwitz condition ^[36,37], all roots of ${\bf{L}}_{1}(\lambda)$ have negative real part if all entries in column 1 of Table 2 have the same sign. We observe that $hz_{1} > 0$ always holds. Clearly, $l_{5}$ is positive if $R_{0}^{a} < 1$. Hence, all eigenvalues of matrix (3.3) have negative real parts if $hz_{3} > 0$, $hz_{4} > 0$, $R_{0}^{a} < 1$, and $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$.

Theorem 3.2. ${\bf{Y}}_{0}^{a}$ is locally asymptotically stable if $hz_{3} > 0$, $hz_{4} > 0$, $R_{0}^{a} < 1$, and $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$. If $R_{0}^{a} > 1$ then ${\bf{Y}}_{0}^{a}$ is unstable.

The method developed in ^[38] is used to study the stability condition of ${\bf{Y}}_{1}^{a}$. Similar approach is used in ^[20,21]. Firstly, we investigate the conditions so that assumption A1 of Theorem 4.1 in ^[38] is fulfilled. Let us choose $\beta_{mv}$ as the bifurcation parameter. We determine the bifurcation point when $R_{0}^{a} = 1$ and obtain $\beta_{mv}^{a*} = \frac{d_{m}d_{h}d_{c}w_{3}w_{2}(w_{5}+\rho N_{v}^{a})}{N_{v}^{a}w_{6}\theta_{ei}w_{1}\beta_{ch}\sigma\Omega_{h}}$. Substituting $\beta_{mv}^{a*}$ into (3.4) produces the following characteristic polynomial of ${\bf{J}}_{{\bf{Y}}_{0}^{a}}\left(\beta_{mv}^{a*}\right)$.

where ${\bf{L}}_{1}^{(1a)}(\lambda) = \lambda^{4}+l_{1}^{(1a)}\lambda^{3}+l_{2}^{(1a)}\lambda^{2}+l_{3}^{(1a)}\lambda^{1}+l_{4}^{(1a)}$ and

$K_{1}^{(1a)} = w_{2}, K_{2}^{(1a)} = w_{3}, K_{3}^{(1a)} = w_{5}+\rho S_{v}^{a*}, K_{4}^{(1a)} = d_{c}$, and $K_{5}^{(1a)} = d_{m}$. Clearly, polynomial (3.5) has three negative roots and a simple zero root. Furthermore, the fourth root is also negative, i.e., $\lambda_{4} = \tau S_{v}^{a*}-d_{p} < 0$ if $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$. The other roots of (3.5) are roots of ${\bf{L}}_{1}^{(1a)}(\lambda)$. The Routh's array associated with ${\bf{L}}_{1}^{(1a)}(\lambda)$ is shown in Table 3.

Based on Routh-Hurwitz condition ^[36,37], all roots of ${\bf{L}}_{1}^{(1a)}(\lambda)$ have negative real part if $hr_{2} > 0$, since the other entries in column 1 are always positive. Thus, if $hr_{2} > 0$ and $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$ then ${\bf{J}}_{{\bf{Y}}_{0}^{a}}(\beta_{mv}^{a*})$ has one zero eigenvalue while the other eigenvalues have negative real part, which implies that assumption A1 of Theorem 4.1 ^[38] is met. The right eigenvector of ${\bf{J}}_{{\bf{Y}}_{0}^{a}}(\beta_{mv}^{a*})$ corresponding to zero eigenvalue is

where $v_{8}$ is arbitrarily positive. Obviously, $v_{1} < 0$ and $v_{i} > 0$ for $i = 2, 3, 4, 5, 6, 7$. The left eigenvector of ${\bf{J}}_{{\bf{Y}}_{0}^{a}}(\beta_{mv}^{a*})$ corresponding to zero eigenvalue is $\vec{k} = \left(k_{1}, k_{2}, k_{3}, k_{4}, k_{5}, k_{6}, k_{7}, k_{8}, k_{9}\right)$, where $k_{1} = k_{4} = k_{5} = k_{9} = 0$, $k_{2} = \frac{\theta_{ei}w_{6}\beta_{mv}^{a*}S_{v}^{a*}}{w_{2}w_{3}d_{m}}k_{6}$, $k_{3} = \frac{w_{6}\beta_{mv}^{a*}S_{v}^{a*}}{w_{3}d_{m}}k_{6}$, $k_{7} = \frac{w_{1}\beta_{ch}S_{h}^{a*}\theta_{ei}w_{6}\beta_{mv}^{a*}S_{v}^{a*}}{d_{c}w_{2}w_{3}d_{m}}k_{6}$, $k_{8} = \frac{\beta_{mv}^{a*}S_{v}^{a*}}{d_{m}}k_{6}$, and $k_{6}$ is chosen, such that $\vec{k}.\vec{v} = 1$. It is easy to show that $k_{6} > 0$. Hence, $k_{j} > 0$ for $j = 2, 3, 7, 8$.

Set $x_{1} = S_{h}, x_{2} = E_{h}, x_{3} = I_{h}, x_{4} = T_{h}, x_{5} = S_{v}, x_{6} = I_{v}, x_{7} = C, x_{8} = M, x_{9} = P$ and $f_{i} = \frac{dx_{i}}{dt}$ for $i = 1...9$. Now, we calculate $a$ and $b$, where

The only non-zero terms of $a$ and $b$ obtained are

Since $v_{1} < 0$, we have

Based on Theorem 4.1 in ^[38], forward bifurcation occurs at $R_{0}^{a} = 1$. Hence, the predator extinction-endemic equilibrium point ${\bf{Y}}_{1}^{a}$ is locally asymptotically stable if $R_{0}^{a} > 1$, $hr_{2} > 0$, and $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$.

Theorem 3.3. ${\bf{Y}}_{1}^{a}$ is locally asymptotically stable if $hr_{2} > 0$, $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$, and $R_{0}^{a} > 1$ (near 1).

Now, we investigate the stability condition of ${\bf{Y}}_{0}^{b}$. The Jacobian matrix of system (2.1) at ${\bf{Y}}_{0}^{b}$ is

The eigenvalues of (3.6) are the zeros of ${\bf{B}}(\lambda)$.

where ${\bf{B}}_{1}(\lambda) = \lambda^{7}+b_{1}\lambda^{6}+b_{2}\lambda^{5}+b_{3}\lambda^{4}+b_{4}\lambda^{3}+b_{5}\lambda^{2}+b_{6}\lambda+b_{7}$ and

$K_{1}^{(0b)} = w_{2}, K_{2}^{(0b)} = w_{3}, K_{3}^{(0b)} = (w_{5}+\rho S_{v}^{b*}+\xi P^{b*}), K_{4}^{(0b)} = d_{c},$ and $K_{5}^{(0b)} = d_{m}$. Noticeably, polynomial (3.7) has two negative roots, i.e., $\lambda_{1} = -d_{h}$ and $\lambda_{2} = -w_{4}$. The other roots of polynomial (3.7) are zeros of ${\bf{B}}_{1}(\lambda)$. It is easy to see that $b_{i} > 0$ for $i = 1, 2, 3, 4$. Additionally, $b_{5} > 0$, $b_{6} > 0$ and $b_{7} >$ if $R_{0}^{b} < 1$. It is easily recognized that $b_{7} = 0$ if $R_{0}^{b} = 1$. Hence, if $R_{0}^{b} = 1$ then one eigenvalue of matrix (3.6) is zero. Now, we use Routh's Table ^[36] to analyze the local stability condition of ${\bf{Y}}_{0}^{b}$.

It is clear that $b_{1} > 0$. According to the Routh-Hurwitz criteria ^[36,37], all roots of ${\bf{B}}_{1}(\lambda)$ have negative real parts if all entries in column 1 of Table 4 have the same sign. We recognize that $hs_{1}$ is always positive. It is obvious that $b_{7} > 0$ if $R_{0}^{b} < 1$. Hence, all eigenvalues of (3.6) have negative real part if $hs_{4} > 0$, $hs_{6} > 0$, $hs_{8} > 0$, $hs_{9} > 0$, and $R_{0}^{b} < 1$. Notice that $b_{7} < 0$ if $R_{0}^{b} > 1$. Hence, if $R_{0}^{b} > 1$ then ${\bf{Y}}_{0}^{b}$ is ustable.

Theorem 3.4. ${\bf{Y}}_{0}^{b}$ is locally asymptotically stable if $hs_{4} > 0$, $hs_{6} > 0$, $hs_{8} > 0$, $hs_{9} > 0$, and $R_{0}^{b} < 1$. If $R_{0}^{b} > 1$, then ${\bf{Y}}_{0}^{b}$ is unstable.

We utilize center manifold theory ^[38] to investigate the local stability condition of ${\bf{Y}}_{1}^{b}$. It is shown that if $R_{0}^{b} = 1$ then matrix (3.6) has one zero eigenvalue. Let us pick $\beta_{mv}$ as bifurcation parameter. We determine the bifurcation point when $R_{0}^{b} = 1$ and obtain $\beta_{mv}^{b*} = \frac{d_{m}d_{h}d_{c}w_{3}w_{2}(w_{5}+\rho N_{v}^{b}+\xi P^{b*})}{N_{v}^{b}w_{6}\theta_{ei}w_{1}\beta_{ch}\sigma\Omega_{h}}$. Clearly, the characteristic polynomial of ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$ is obtained by substituting $\beta_{mv}^{b*}$ into (3.7). Hence we get (3.8) as a characteristic polynomial of ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$.

where ${\bf{B}}_{1}^{(1b)}(\lambda) = \lambda^{6}+b_{1}^{(1b)}\lambda^{5}+b_{2}^{(1b)}\lambda^{4}+b_{3}^{(1b)}\lambda^{3}+b_{4}^{(1b)}\lambda^{2}+b_{5}^{(1b)}\lambda+b_{6}^{(1b)}$ and

$K_{1}^{(1b)} = w_{2}, K_{2}^{(1b)} = w_{3}, K_{3}^{(1b)} = w_{5}+\rho S_{v}^{b*}+\tau P^{b*}, K_{4}^{(1b)} = d_{c},$ and $K_{5}^{(1b)} = d_{m}$. Thus, ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$ has one zero eigenvalue and two negative eigenvalues, i.e., $-d_{h}$ and $-w_{4}$. The sign of the other roots of (3.8) is studied by investigating the sign of the roots of polynomial ${\bf{B}}_{1}^{(1b)}(\lambda)$. Hence, we utilize Routh-Hurwitz test.

According to Routh-Hurwitz condition ^[36,37], all roots of ${\bf{B}}_{1}^{(1b)}(\lambda)$ have negative real part if $ht_{1}$, $ht_{3}$, $ht_{5}$ and $ht_{6}$ are positive, since the other entries in column 1 of Table 5 are positive, i.e., $b_{1}^{(1b)} > 0$ and $b_{6}^{(1b)} > 0$. Consequently, ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$ has a zero eigenvalue and the other eigenvalues have negative real part if $ht_{1} > 0$, $ht_{3} > 0$, $ht_{5} > 0$ and $ht_{6} > 0$. Therefore, we now investigate the left and the right eigenvector of ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$ corresponding to zero eigenvalue. The right eigenvector of ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$ corresponding to zero eigenvalue is as follows:

where $m_{8}$ is arbitrarily positive. It is not difficult to show that $m_{1} < 0$, $m_{5} < 0$ and $m_{i} > 0$ for $i = 2, 3, 4, 6, 7$. The left eigenvector of ${\bf{J}}_{{\bf{Y}}_{0}^{b}}\left(\beta_{mv}^{b*}\right)$ corresponding to zero eigenvalue is $\vec{h} = \left(h_{1}, h_{2}, h_{3}, h_{4}, h_{5}, h_{6}, h_{7}, h_{8}, h_{9}\right)$, where $h_{1} = h_{4} = h_{5} = h_{9} = 0$, $h_{2} = \frac{\theta_{ei}w_{6}\beta_{mv}^{b*}S_{v}^{b*}}{w_{2}w_{3}d_{m}}h_{6}$, $h_{3} = \frac{w_{6}\beta_{mv}^{b*}S_{v}^{b*}}{w_{3}d_{m}}h_{6}$, $h_{7} = \frac{w_{1}\beta_{ch}S_{h}^{b*}\theta_{ei}w_{6}\beta_{mv}^{b*}S_{v}^{b*}}{d_{c}w_{2}w_{3}d_{m}}h_{6}$, $h_{8} = \frac{\beta_{mv}^{b*}S_{v}^{b*}}{d_{m}}h_{6}$, and $h_{6}$ is chosen, such that $\vec{h}.\vec{m} = 1$. It is easy to show that $h_{6} > 0$. Obviously, $h_{j} > 0$ for $j = 2, 3, 7, 8$.

Suppose $x_{1} = S_{h}, x_{2} = E_{h}, x_{3} = I_{h}, x_{4} = T_{h}, x_{5} = S_{v}, x_{6} = I_{v}, x_{7} = C, x_{8} = M, x_{9} = P$ and $f_{i} = \frac{dx_{i}}{dt}$ for $i = 1...9$. Now, we compute $a$ and $b$, where

The only non-zero terms of $a$ and $b$ are

Since $m_{1} < 0$, we obtain

Based on Theorem 4.1 in ^[38], forward bifurcation occurs at $R_{0}^{b} = 1$. Hence, the endemic equilibrium point ${\bf{Y}}_{1}^{b}$ is locally asymptotically stable if $R_{0}^{b} > 1$, $ht_{1} > 0$, $ht_{3} > 0$, $ht_{5} > 0$ and $ht_{6} > 0$.

Theorem 3.5. ${\bf{Y}}_{1}^{b}$ is locally asymptotically stable if $ht_{1} > 0$, $ht_{3} > 0$, $ht_{5} > 0$, $ht_{6} > 0$ and $R_{0}^{b} > 1$ (near 1).

3.2.   Global stability

To investigate the global behavior of system (2.1), we use the theory of asymptotically autonomous systems ^[39,40], i.e., by studying the limiting systems related to system (2.1). Diaby et al. ^[41] used this theory to study the dynamics of schistosomiasis model. Since $S_{v}+I_{v} = N_{v}$, system (2.1) is equivalent to the following system.

As mentioned in ^[42], global stability can be used as the basis for formulating a limiting system. Nakata ^[43] studied the global dynamics of the total populations, though before constructing a limiting system. Here, we consider studying the global dynamics of snail and snail predator populations before we construct the limiting systems of (3.9), we examine the global behavior of the system formed by the last two equations of (3.9), i.e., system (3.10).

The equilibrium points of system (3.10) are $\mathscr{P}_{0} = (0, 0)$, which always exists, $\mathscr{P}_{1} = (N_{v}^{a}, 0) = \left(\frac{\varphi-w_{5}}{\rho}, 0\right)$ which exists if $\frac{\varphi}{w_{5}} > 1$, and $\mathscr{P}_{1} = (N_{v}^{b}, P^{b*}) = \left(\frac{d_{p}}{\tau}, \frac{\tau(\varphi-w_{5})-\rho d_{p}}{\xi \tau}\right)$, which exists if $\frac{\tau(\varphi-w_{5})}{d_{p}\rho} > 1$.

Obviously, $\mathscr{P}_{0}$ and $\mathscr{P}_{2}$ represent the situation in which both snail and snail predator populations become extinct and survive, respectively. In addition, the condition in which the snail population survives but the snail predator population becomes extinct is denoted by $\mathscr{P}_{1}$. We now present the stability condition of the equilibrium points of the system (3.10).

Theorem 3.6. (i). $\mathscr{P}_{0}$ is global asymptotically stable if $\frac{\varphi}{w_{5}} < 1$.

(ii). $\mathscr{P}_{1}$ is global asymptotically stable if $\frac{\tau(\varphi-w_{5})}{d_{p}\rho} < 1$.

(iii). $\mathscr{P}_{2}$ is global asymptotically stable if it exists.

Proof. (ⅰ). First, we prove the stability condition of $\mathscr{P}_{0}$. We define the following Lyapunov function.

The derivative of $\mathscr{V}$ with respect to $t$ is given by

It is easy to see that $\frac{d\mathscr{V}}{dt}\leq0$ if $\varphi-w_{5} < 0$. Furthermore, $\frac{d\mathscr{V}}{dt} = 0$ if and only if $N_{v} = 0$ and $P = 0$. Thus, $\mathscr{P}_{0}$ is global asymptotically stable if $\frac{\varphi}{w_{5}} < 1$.

(ⅱ). Now, we prove the stability condition of $\mathscr{P}_{1}$. Consider the following Lyapunov function.

The time derivative of $\mathscr{L}$ is given by

Clearly, $\frac{d\mathscr{L}}{dt}\leq 0$ if $\frac{\tau N_{v}^{a}}{d_{p}} < 1$. Furthermore, $\frac{d\mathscr{L}}{dt} = 0$ if and only if $(N_{v}, P) = (N_{v}^{a}, 0)$. Thus, $\mathscr{P}_{1}$ is global asymptotically stable if $\frac{\tau N_{v}^{a}}{d_{p}} = \frac{\tau (\varphi-w_{5})}{ d_{p}\rho} < 1$.

(ⅲ). We now prove the stability condition of $\mathscr{P}_{2}$. Consider the following Lyapunov function.

The time derivative of $\mathscr{M}$ is given by

From the equilibrium $\mathscr{P}_{2}$, we have $\tau N_{v}^{b} = d_{p}$ and $\varphi-w_{5} = \xi P^{b*}+\frac{\rho d_{p}}{\tau}$. Hence, we obtain

It is clear that $\frac{d\mathscr{M}}{dt} = 0$ if and only if $N_{v} = N_{v}^{b}$. Substituting $N_{v}^{b}$ into system (3.10) shows that $P = P^{b*}$. Thus, $\mathscr{P}_{2}$ is global asymptotically stable if it exists.

According to ^[40], we can investigate the dynamics of Systems (2.1) and (3.9) by studying the dynamics of the limiting systems. Based on the global behavior of the equilibrium points of System (3.10), we divide the process into three parts. In this paper, System (3.9) is called an asymptotically autonomous system with limit Systems (3.11)–(3.13). Theorem 3.6 (ⅰ) guarantees that as $t\to \infty$, $(N_{v}, P)\to (0, 0)$ if $\frac{\varphi}{w_{5}} < 1$. Notice that $S_{v}+I_{v} = N_{v} = 0$ implies that $S_{v} = 0$ and $I_{v} = 0$. Hence, we obtain the following limiting system of (3.9).

We now investigate the stability condition of the equilibrium point of limiting system (3.11), namely $\mathscr{Y}_{0} = \left(S_{h}^{*0}, E_{h}^{*0}, I_{h}^{*0}, T_{h}^{*0}, C^{*0}, M^{*0}\right) = \left(\frac{\Omega_{h}}{d_{h}}, 0, 0, 0, 0, 0\right)$.

Theorem 3.7. $\mathscr{Y}_{0}$ is global asymptotically stable.

Proof. We define the following Lyapunov function.

The derivative of $\mathscr{U}$ with respect to $t$ is given as follows:

From the equilibrium point $\mathscr{Y}_{0}$, we have $\Omega_{h} = d_{h}S_{h}^{*0}$. Thus, we obtain

Since $S_{h}\leq\frac{\Omega_{h}}{d_{h}} = S_{h}^{*0}$, we conclude that $\frac{d\mathscr{U}}{dt}\leq 0$. Further, $\frac{d\mathscr{U}}{dt} = 0$ if and only if $\frac{d\mathscr{U}}{dt}$ is evaluated at $\mathscr{Y}_{0}$. Therefore, $\mathscr{Y}_{0}$ is global asymptotically stable.

Theorem 3.6 (ⅱ) assures that $(N_{v}, P)\to (N_{v}^{a}, 0)$ as $t \to \infty$ if $\frac{\tau (\varphi-w_{5})}{ d_{p}\rho} < 1$. Thus, we have the following limiting system of (3.9).

We now investigate the global stability condition of the equilibrium point of System (3.12), i.e., $\mathscr{Y}_{0}^{a} = (S_{h}^{a*}, E_{h}^{a*}, I_{h}^{a*}, T_{h}^{a*}, I_{v}^{a*}, C^{a*}, M^{a*}) = \left(\frac{\Omega_{h}}{d_{h}}, 0, 0, 0, 0, 0, 0\right)$.

Theorem 3.8. $\mathscr{Y}_{0}^{a}$ is global asymptotically stable if $R_{0}^{a}\leq 1$.

Proof. Consider the following Lyapunov function:

where

The time derivative of $Z$ is

From the equilibrium point $\mathscr{Y}_{0}^{a}$, we have $\Omega_{h} = d_{h}S_{h}^{a*}$. Hence, we obtain

Certainly, $S_{h}\leq S_{h}^{a*}$. Hence, it is easy to see that $\frac{dZ}{dt} \leq 0$ if $R_{0}^{a}\leq1$. Notice that $\frac{dZ}{dt} = 0$ if and only if $S_{h} = S_{h}^{a*}$ which implies that $E_{h} = I_{h} = T_{h} = 0$. Hence, the largest invariant set contained in $\left\{(S_{h}, E_{h}, I_{h}, T_{h}, I_{v}, C, M)|\frac{dZ}{dt} = 0\right\}$ is a singleton set $\left\{\mathscr{Y}_{0}^{a}\right\}$. Therefore, $\mathscr{Y}_{0}^{a}$ is global asymptotically stable if $R_{0}^{a}\leq 1$.

Theorem 3.6 (ⅲ) ensures that $(N_{v}, P)\to (N_{v}^{b}, P^{b*})$ as $t \to \infty$ if $\frac{\tau (\varphi-w_{5})}{ d_{p}\rho} > 1$. Thus, we have the following limiting system of (3.9).

Next, we determine the global stability condition of the equilibrium point of system (3.13), i.e., $\mathscr{Y}_{0}^{b} = (S_{h}^{b*}, E_{h}^{b*}, I_{h}^{b*}, T_{h}^{b*}, I_{v}^{b*}, C^{b*}, M^{b*}) = \left(\frac{\Omega_{h}}{d_{h}}, 0, 0, 0, 0, 0, 0\right)$.

Theorem 3.9. $\mathscr{Y}_{0}^{b}$ is global asymptotically stable if $R_{0}^{b}\leq 1$.

Proof. Consider the Lyapunov function as follows:

where

The time derivative of $\mathscr{Z}$ is

From the equilibrium point $\mathscr{Y}_{0}^{b}$, we have $\Omega_{h} = d_{h}S_{h}^{b*}$. Hence, we obtain

Obviously, $S_{h}\leq S_{h}^{b*}$. Hence, it is easy to see that $\frac{d\mathscr{Z}}{dt} \leq 0$ if $R_{0}^{b}\leq1$. Notice that $\frac{d\mathscr{Z}}{dt} = 0$ if and only if $S_{h} = S_{h}^{b*}$ which implies that $E_{h} = I_{h} = T_{h} = 0$. Hence, the largest invariant set contained in $\left\{(S_{h}, E_{h}, I_{h}, T_{h}, I_{v}, C, M)|\frac{d\mathscr{Z}}{dt} = 0\right\}$ is a singleton set $\left\{\mathscr{Y}_{0}^{b}\right\}$. Hence, $\mathscr{Y}_{0}^{b}$ is global asymptotically stable if $R_{0}^{b}\leq 1$.

Remark $\mathscr{Y}_{0}$ of system (3.11) is qualitatively equivalent to the equilibrium ${\bf{Y}}_{0}$ of system (2.1) when $\frac{\varphi}{w_{5}} < 1$. Hence, we conclude that ${\bf{Y}}_{0}$ is globally asymptotically stable if $\frac{\varphi}{w_{5}} < 1$, since $\mathscr{Y}_{0}$ is globally asymptotically stable if $\frac{\varphi}{w_{5}} < 1$. This result suggests that when the snail and the snail predator populations become extinct, provided by $\frac{\varphi}{w_{5}} < 1$, schistosomiasis will be eradicated.

$\mathscr{Y}_{0}^{a}$ is qualitatively equivalent to the equilibrium ${\bf{Y}}_{0}^{a}$ of system (2.1) if $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$. Thus, ${\bf{Y}}_{0}^{a}$ is globally asymptotically stable if $R_{0}^{a}\leq1$ and $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} < 1$, since $\mathscr{Y}_{0}^{a}$ of is globally asymptotically stable if $R_{0}^{a}\leq1$. This result indicates that when the snail population survives and the population of snail predators becomes extinct, provided by $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} < 1$, schistosomiasis will be successfully eradicated if $R_{0}^{a}\leq 1$.

$\mathscr{Y}_{0}^{b}$ is qualitatively equivalent to the equilibrium point of system (2.1), i.e., ${\bf{Y}}^{b}_{0}$ if $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} > 1$. Thus, ${\bf{Y}}_{0}^{b}$ is globally asymptotically stable if $R_{0}^{b}\leq1$ and $\frac{\tau (\varphi-w_{5})}{d_{p}\rho} > 1$, since $\mathscr{Y}^{b}_{0}$ is globally asymptotically stable if $R_{0}^{b}\leq1$. This result implies that when both snail and snail predator populations survive, provided by $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} > 1$, the disease will die out if $R_{0}^{b}\leq 1$.

4.   Numerical simulations

We conduct numerical simulations for system (2.1) to validate and support the previous theoretical results. The numerical simulations are performed by using parameter values given in Table 6.

4.1.   Snail and snail predator go to extinction

We perform numerical simulation using parameter values given in Table 6 and $\tau = 0.00001$. We set $d_{r} = 0.8$, which indicates the high use of molluscicides. These parameter values yield $\frac{\varphi}{d_{v}+d_{r}} = 0.4996 < 1$. Based on Theorem 3.1, ${\bf{Y}}_{0}$ is locally asymptotically stable. This result is the same as the numerical simulation result shown in Figure 2.

Figure 2(a) shows that snail and snail predator populations lead to zero. It suggests that both populations go to extinction. Furthermore, the infectious human population tends to zero. It indicates that the disease will be eradicated. This result is also supported by the phase portrait displayed in Figure 2(b). Notice that all solution trajectories are towards ${\bf{Y}}_{0}$.

Therefore, the use of molluscicides on a large scale has both positive and negative impacts. This kind of intervention can suppress the spread of schistosomiasis, but it can disrupt the environmental balance, e.g., snail predators may become extinct.

4.2.   Snail survives but snail predator goes to extinction

Now, we present the second numerical simulation results. The parameter values used are presented in Table 6 with $\tau = 0.00001$, which gives $\frac{\varphi}{d_{v}+d_{r}} = 6.9084 \times 10^{2} > 1$ and $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} = 0.2684 < 1$. Thus, the snail population survives, but the snail predator goes to extinction. Moreover, we get a bifurcation point for $\beta_{mv}$, i.e., $\beta_{mv}^{a*} = 5.6887\times10^{-5}$ which corresponds to $R_{0}^{a} = 1$ and $hr_{2} = 1.8276 > 0$. Using $\beta_{mv} = 5.6887\times10^{-6} < \beta_{mv}^{a*}$ implies $R_{0}^{a} = 0.3162 < 1$, $hz_{3} = 1.8280 > 0$, and $hz_{4} = 0.1143 > 0$. Based on Theorem 3.2, ${\bf{Y}}_{0}^{a}$ is asymptotically stable. The numerical simulation result is shown in Figure 3(a). If we set $\beta_{mv} = 5.6887\times10^{-4} > \beta_{mv}^{a*}$, we get $R_{0}^{a} = 3.1626 > 1$. Based on Theorem 3.3, ${\bf{Y}}_{1}^{a}$ is asymptotically stable. The numerical simulation result can be seen in Figure 3(b). Figure 3(a) shows that the disease will be eradicated. On the other hand, Figure 3(b) shows that the disease will be endemic. Furthermore, the snail population does not tend to zero, while the predator population tends to zero. These results tell us that we can eradicate schistosomiasis even though the snail population is extant.

Figure 3(c) displays a phase portrait of system (2.1) when$\frac{\varphi}{d_{v}+d_{r}} = 6.9084 \times 10^{2} > 1$, $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} = 0.2684 < 1$, and $R_{0}^{a} = 0.3162 < 1$. It is clear that all solution trajectories tend to the equilibrium point ${\bf{Y}}_{0}^{a}$. Notice that ${\bf{Y}}_{1}^{a}$ does not exist. On the other hand, when $R_{0}^{a} = 3.1626 > 1$, the equilibrium point ${\bf{Y}}_{1}^{a}$ exists, and all solution trajectories go to ${\bf{Y}}_{1}^{a}$ as illustrated in Figure 3(d).

4.3.   Snail and snail predator survive

Now, we present the third numerical simulation result. The parameter values used are presented in Table 6, and $\tau = 0.0001$, which gives $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} = 2.6841 > 1$, $\beta_{mv}^{b*} = 1.5265\times 10^{-4}$. Substituting $\beta_{mv} = 1.5265\times 10^{-5} < \beta_{mv}^{b*}$ gives $R_{0}^{b} = 0.3162 < 1$, $hs_{4} = 2.2101 > 0$, $hs_{6} = 0.4636 > 0$, $hs_{8} = 0.0263 >$, and $hs_{9} = 0.0003 > 0$. Based on Theorem 3.4, ${\bf{Y}}_{0}^{b}$ is asymptotically stable. This result is confirmed by Figure 4(a). If we use $\beta_{mv} = 1.5265\times 10^{-4} = \beta_{mv}^{b*}$, we get $R_{0}^{b} = 1$, $ht_{1} = 5.1379 > 0$, $ht_{3} = 1.1616 > 0$, $ht_{5} = 0.1338 > 0$, and $ht_{6} = 0.0071 > 0$. Furthermore, $R_{0}^{b} = 3.1622 > 1$ is obtained if we set $\beta_{mv} = 1.5265\times 10^{-3} > \beta_{mv}^{b*}$. Theorem 3.5 states that ${\bf{Y}}^{b}_{1}$ is locally asymptotically stable. This result is similar to the numerical simulation result shown in Figure 4(b). Figure 4(a) shows that the disease will be eradicated. Nevertheless, Figure 4(b) shows that the disease will be endemic. Furthermore, the snail and its predator population do not tend to zero. These results tell us that we can eradicate schistosomiasis even though the snail and snail predator populations are extant.

Figure 4(c) illustrates the phase portrait of the system (1) when $\frac{\tau(\varphi-w_{5})}{\rho d_{p}} = 2.6841 > 1$ and $R_{0}^{b} = 0.3162 < 1$. We can notice that the equilibrium points ${\bf{Y}}_{0}^{a}$, ${\bf{Y}}_{0}^{b}$, and ${\bf{Y}}_{1}^{a}$ exist while ${\bf{Y}}_{1}^{b}$ does not exist. It can be seen that all solution trajectories tend to the equilibrium point ${\bf{Y}}_{0}^{b}$. Otherwise, when $R_{0}^{b} = 3.1622 > 1$, ${\bf{Y}}_{0}^{a}$, ${\bf{Y}}_{0}^{b}$, ${\bf{Y}}_{1}^{a}$ and ${\bf{Y}}_{1}^{b}$ exist. However, it can be witnessed, from Figure 4(d), that all solution trajectories go to ${\bf{Y}}_{1}^{b}$.

4.4.   Effect of control measures

In this section, we investigate the effect of control measures, i.e., $\theta_{et}$, $\rho$, $\tau$ and $\xi$, on schistosomiasis spread dynamics. The parameter values used are listed in Table 6 with $\tau = 0.0001$ and $\beta_{mv} = 1.5265\times 10^{-3}$. After performing numerical simulations, we obtained the following results.

We firstly studied the impact of the average waiting time of a latent human to get treatment on schistosomiasis prevalence. $\theta_{et}$ is varied. The results are shown in Figure 5(a). The figure shows that the prevalence of schistosomiasis decreases as the average waiting time $(1/\theta_{et})$ decreases. In addition, the basic reproduction numbers $(R_{0}^{b})$ obtained are $3.3988$, $3.1624$ and $2.4300$ for $\theta_{et} = 0.0071$, $\theta_{et} = 0.0119$ and $\theta_{et} = 0.0357$, respectively. It is clear that the basic reproduction number decreases as $\theta_{et}$ increases. Agbata et al. ^[44] stated that a screening program could increase the possibility of early detection and treatment. Hence, we can increase $\theta_{et}$ by implementing a screening program followed by early treatment. Therefore, our result shows that a screening program followed by early treatment can reduce the prevalence of schistosomiasis.

Figure 5(b) shows the effect of habitat modification on the spread of schistosomiasis. We varied $\rho$ while the other parameter values are fixed as given in Table 6, $\tau = 0.0001$ and $\beta_{mv} = 1.5265\times 10^{-3}$. The basic reproduction numbers $(R_{0}^{b})$ gained is $3.1622$ for $\rho = 0.001$, $\rho = 0.01$ and $\rho = 0.1$. It can be seen that the number of infectious human at the beginning of intervention decreases as $\rho$ increases. Here, $\rho$ is the competition rate of the snails. According to ^[45], habitat change can affect interspecific and intraspecific competition. It is already known that intraspecific competition is competition of members of the same species for limited resources. Thus, snail habitat modification, which may increase the competition rate in the snail population, can reduce the prevalence of schistosomiasis at the beginning of intervention. Xu et al.^[17] also stated that an intervention that can be used to control the spread of schistosomiasis is snail habitat modification.

Figure 5(c) shows the effects of the conversion rate of snail predators on the dynamics of infectious humans. We varied $\tau$. The other parameter values are fixed as given in Table 6, and $\beta_{mv} = 1.5265\times 10^{-3}$. It is seen that the number of infectious humans decreases as $\tau$ increases. Furthermore, the basic reproduction values $(R_{0}^{b})$ obtained are $3.1622$, $2.5819$ and $2.2360$ for $\tau = 0.0001$, $\tau = 0.00015$ and $\tau = 0.0002$, respectively. The conversion rate is related to the birth rate. Therefore, the birth rate of snail predators has an important role in reducing schistosomiasis cases. This result is similar to the result given in ^[27].

Figure 5(d) shows the dynamics of infectious humans with varying $\xi$. The other parameter values used are given in Table 6 with $\beta_{mv} = 1.5265\times 10^{-3}$ and $\tau = 0.0001$. The basic reproduction number $(R_{0}^{b})$ obtained is $3.1622$ for $\xi = 0.001$, $\xi = 0.005$ and $\xi = 0.01$. The number of infectious humans at the beginning of the outbreak decreases as $\xi$ increases. Here, $\xi$ is predation rate, which is related to the proportion of snails killed per predator per time. Thus, releasing snail predators at snail habitats can reduce the prevalence of schistosomiasis. We recommend using snail predators which can hunt and kill snails effectively, for example, river prawn ^[13], to maximize the effect of intervention.

5.   Conclusions

A schistosomiasis model with treatment, habitat modification and biological control is discussed in this work. Our results show that the basic reproduction number is inversely proportional to $\theta_{et}$. It means that a screening program followed by early treatment can reduce the basic reproduction number and schistosomiasis prevalence. On the other hand, the basic reproduction number is independent of predation rate and competition rate of the snails. However, modifying the snail habitat and releasing snail predators at the snail habitat can reduce the prevalence of schistosomiasis at the beginning of intervention. To maximize this effect, we should use snail predators which can hunt and kill snails effectively.

Acknowledgments

This research is funded by the Directorate of Research and Community Service, Ministry of Research and Technology / National Research and Innovation Agency, Indonesia, via Doctoral Dissertation Research, in accordance with the Research Contract No. 439.17/ UN.10.C10/ TU/ 2021, dated 19 March 2021.

Conflict of interest

The authors declare there is no conflict of interest.