Global Transmission Dynamics of a Schistosomiasis Model and Its Optimal Control

Mouhamadou Diaby¹, Mariama Sène², Abdou Sène³

¹Laboratoire d’Analyse Numérique et d’Informatique (LANI), UFR SAT, Université Gaston Berger de Saint-Louis, Saint-Louis, Sénégal.
²UFR AGRO, Université Gaston Berger de Saint-Louis, Saint-Louis, Sénégal.
³LANI, Pôle STN MAI (Mathématiques Appliquées et Informatique), Université Virtuelle du SENEGAL, Dakar-Fann, Sénégal.

DOI: 10.4236/am.2019.106029                

Abstract

Drug treatment, snail control, cercariae control, improved sanitation and health education are the effective strategies which are used to control the schistosomiasis. In this paper, we consider a deterministic model for schistosomiasis transmission dynamics in order to explore the role of the several control strategies. The global stability of a schistosomiasis infection model that involves mating structure including male schistosomes, female schistosomes, paired schistosomes and snails is studied by constructing appropriate Lyapunov functions. We derive the basic reproduction number R₀ for the deterministic model, and establish that the global dynamics are completely determined by the values of R₀. We show that the disease can be eradicated when R₀ ≤1; otherwise, the system is persistent. In the case where  R₀ >1, we prove the existence, uniqueness and global asymptotic stability of an endemic steady state. Sensitivity analysis and simulations are carried out in order to determine the relative importance of different control strategies for disease transmission and prevalence. Next, optimal control theory is applied to investigate the control strategies for eliminating schistosomiasis using time dependent controls. The characterization of the optimal control is carried out via the Pontryagins Maximum Principle. The simulation results demonstrate that the insecticide is important in the control of schistosomiasis.

Keywords

Schistosomiasis Models, Nonlinear Dynamical Systems, Global Stability, Reproduction Number, Optimal Control, Sensitivity Analysis

1. Introduction

Schistosomiasis (also known as bilharzia, bilharziasis or snail fever) is a vector-borne disease caused by infection of the intestinal or urinary venous system by trematode worms of the genus Schistosoma. More than 220.8 million people required preventive treatment worldwide in 2017, out of which more than 102.3 million people were reported to have been treated [1]. Schistosomiasis is prevalent in tropical and subtropical areas, especially in poor communities without access to safe drinking water and adequate sanitation. Of the 207 million people with schistosomiasis, 85% live in Africa [1]. Of the tropical diseases, only malaria accounts for a greater global burden than schistosomiasis [2]. Therefore, it is vital to prevent and control the schistosomiasis transmission.

Schistosoma requires the use of two hosts to complete its life cycle: the definitive hosts and the intermediate snail hosts. In definitive hosts, schistosoma has two distinct sexes. Mature male and female worms pair and migrate either to the intestines or the bladder where eggs production occurs. One female worm may lay an average of 200 to 2000 eggs per day for up to twenty years. Most eggs leave the bloodstream and body through the intestines. Some of the eggs are not excreted, however, and can lodge in the tissues. It is the presence of these eggs, rather than the worms themselves, that causes the disease. These eggs pass in urine or feces into fresh water into miracidia which infect the intermediate snail hosts. In snail hosts, parasites undergo further asexual reproduction, ultimately yielding large numbers of the second free-living stage, the cercaria. Free-swimming cercariae leave the snail host and move through the aquatic or marine environment, often using a whip-like tail, though a tremendous diversity of tail morphology is seen. Cercariae are infective to the second host and turn it into single schistosoma, and infection may occur passively (e.g., a fish consumes a cercaria) or actively (the cercaria penetrates the fish) and terminates the life cycle of the parasite.

Many effective strategies are used in the real world, such as: based on preventive treatment, snail control, cercariae control, improved sanitation and health education. The WHO strategy for schistosomiasis control focuses on reducing disease through periodic, targeted treatment with praziquantel. This involves regular treatment of all people in at-risk groups [1]. Over the past few decades, different mathematical models [3] [4] [5] [6] have been constructed to describe the transmission dynamics involving two-sex problems. In [3] [4] [5] , a mathematical model is developed for a schistosomiasis infection that involves pair-formation models and studied the existence, uniqueness and the stabilities of exponential solutions. We note that in [4] [5] authors formulate three forms of pair-formation functions (also known as mating functions) that are the harmonic mean function, the geometric mean function and the minimum function. In [7] , Xu et al. have proposed a multi-strain schistosome model with mating structure. Their goal was to study the effect of drug treatment on the maintenance of schistosome genetic diversity. However, in their model they only consider the adult parasite populations. Castillo-Chavez et al. [3] have considered a time delay model but also do not include the snails dynamics. But it is important to take into account the snail dynamics as it is shown in the life cycle of schistosoma. In fact, the parasite offspring is produced directly by infected snails but not by paired parasites as is related in [6]. Recently, Qi et al. [6] have formulated a deterministic mathematical model to study the transmission dynamics of schistosomiasis with a linear mating function incorporating these snail dynamics. This paper gave the expression of a threshold number (and not the basic reproduction number) with a local stability analysis of the disease free equilibrium. The sensitivity analysis of this threshold number is also discussed.

However, no work has been done to investigate the global stability of the equilibria which is more in interest. Here, we take this deterministic schistosomiasis model with mating structure [6] and we propose a complete mathematical analysis. A stability analysis is provided to study the epidemiological consequences of control strategies. We compute the basic reproduction number and we show that when it is less or equal to one then the disease free equilibrium (DFE) is the unique equilibrium of the system and it is globally asymptotically stable, while when the basic reproduction number is greater than one then there is a unique endemic equilibrium which is globally asymptotically stable in the whole space minus the stable manifold of the DFE. A sensitivity analysis of the endemic equilibrium is performed giving a more interest interpretation of the control strategies. Optimal control is a branch of mathematics developed to find optimal ways to control a dynamic system [8] [9]. There are few papers that apply optimal control to schistosomiasis models. Here we propose and analyze one such optimal control problem, where the control function represents the fraction of snails individuals ( $X_s$ and $X_i$ ) that will be submitted to treatment. The objective is to find the optimal treatment strategy through insecticide campaigns that minimizes the number of snails individuals as well as the cost of interventions. This paper is organized as follows. Model formulation is carried out and the basic properties are shown in the next section. In Section 3, we determine the basic reproductive number ${\mathcal{R}}_0$ of the model and also establish local and global stability of the disease-free equilibrium. In the end of this section we show that the disease is uniformly persistent when ${\mathcal{R}}_0>1$. Section 4 is devoted to prove the global asymptotic stability of the endemic equilibrium. In Section 5, a sensitivity analysis of the basic reproductive number and of the endemic equilibrium are explored. The goal is to identify the most sensitive parameter allowing decreasing the disease prevalence. In Section 6 we propose and analyze an optimal control problem. A general conclusion is given in the last section.

2. Mathematical Model

The model that we consider has been presented in [6]. It describes the time evolution of a population divided in three parasites sub-populations and two intermediate snail host sub-populations. The state variables of the model are:

· $X_m\left(t\right)$ the male schistosoma population size.

· $X_f\left(t\right)$ the female schistosoma population size.

· $X_p\left(t\right)$ the pair schistosoma population size.

· $X_s\left(t\right)$ the susceptible (uninfected) snail host population size.

· $X_i\left(t\right)$ the infected snail host population size.

The time evolution of the different populations is governed by the following system of equations:

$\{\begin{array}{l}\frac{\text{d}{X}_m}{\text{d}t}=k_m{X}_i-\left({\mu }_m+{ϵ}_m\right)X_m-\rho X_f,\\ \frac{\text{d}{X}_f}{\text{d}t}=k_f{X}_i-\left({\mu }_f+{ϵ}_f\right)X_f-\rho X_f,\\ \frac{\text{d}{X}_p}{\text{d}t}=\rho X_f-\left({\mu }_p+{ϵ}_p\right)X_p,\\ \frac{\text{d}{X}_s}{\text{d}t}=\Lambda -\left({\mu }_s+{ϵ}_s\right)X_s-\beta X_p{X}_s,\\ \frac{\text{d}{X}_i}{\text{d}t}=\beta X_p{X}_s-\left({\mu }_s+{ϵ}_s+{\alpha }_s\right)X_i.\end{array}$ (1)

The different parameters are:

· $k_m$ and $k_f$ are the recruitment rates of male schistosoma and female schistosoma respectively.

· ${\mu }_m$, ${\mu }_f$, ${\mu }_p$, and ${\mu }_s$ denote the natural death rate for male, female, pair and snail hosts respectively. ${\alpha }_s$ is the disease-induced death rate of snail hosts.

· $\rho$ represents the effective mating rate.

· $\Lambda$ is the recruitment rate of snail hosts.

· $\beta$ is the transmission rate from pairs parasite to susceptible snails.

· ${ϵ}_m$, ${ϵ}_f$, ${ϵ}_p$ and ${ϵ}_s$ are the elimination rates of male shistosoma, female schistosoma, paired schistosoma and snails respectively. These elimination rates represent the control strategies.

As it has been done in [6] , we shall denote

${\mu }_m+{ϵ}_m={\mu }_{mϵ}$, ${\mu }_f+{ϵ}_f={\mu }_{fϵ}$,

${\mu }_p+{ϵ}_p={\mu }_{pϵ}$, ${\mu }_s+{ϵ}_s={\mu }_{sϵ}$.

2.1. Basic Properties

In this section, we give some basic results concerning solutions of system (1) that will be subsequently used in the proofs of the stability results.

Proposition 1 The set $\Gamma =\left\{X_m\ge X_f\ge 0,X_p\ge 0,X_s\ge 0,X_i\ge 0\right\}$ is a positively invariant set for system (1).

Proof. The vector field given by the right-hand side of system (1) points inward on the boundary of $\Gamma$. For example, if $X_s=0$, then, ${\stackrel{\dot{}}{X}}_s=\Lambda >0$. In an analogous manner, the same can be shown for the other system components.

Proposition 2 All solutions of system (1) are forward bounded.

Proof. Let us define $N_X=X_m+X_f+X_p$ and $N_Y=X_s+X_i$. Using system (1), we have $\frac{\text{d}{N}_Y}{\text{d}t}=\Lambda -{\mu }_{sϵ}{N}_Y-{\alpha }_s{X}_i\le \Lambda -{\mu }_{sϵ}{N}_Y$. This implies that the set $\left\{N_Y\le \frac{\Lambda }{{\mu }_{sϵ}}\right\}$ is positively invariant and attracts all the solutions of (1).

We also have:

$\begin{array}{c}\frac{\text{d}{N}_X}{\text{d}t}=\left(k_m+k_f\right)X_i-{\mu }_{mϵ}{X}_m-\left({\mu }_{fϵ}+\rho \right)X_f-{\mu }_{pϵ}{X}_p\\ \le \left(k_m+k_f\right)\frac{\Lambda }{{\mu }_{s\epsilon }}-\min \left\{{\mu }_{mϵ}\mathrm{,}{\mu }_{fϵ}\mathrm{,}{\mu }_{pϵ}\right\}{N}_X-\rho X_f\mathrm{.}\end{array}$

Hence, the set $\left\{N_X\le \frac{\left(k_m+k_f\right)\Lambda }{{\mu }_{sϵ}\gamma }\right\}$, where $\gamma =\min \left\{{\mu }_{mϵ}\mathrm{,}{\mu }_{fϵ}\mathrm{,}{\mu }_{pϵ}\right\}$, is positively invariant set and attracts all the solutions of (1).

Therefore all feasible solutions of system (1) enter the region

$\begin{array}{l}\Omega =\{\left(X_m\mathrm{,}{X}_f\mathrm{,}{X}_p\mathrm{,}{X}_s\mathrm{,}{X}_i\right)\in {ℝ}_{+}^5\mathrm{<mo>\; :\; </mo>}{X}_s+X_i\le \frac{\Lambda }{{\mu }_{sϵ}}\mathrm{,}\\ X_m+X_f+X_p\le \frac{\left(k_m+k_f\right)\Lambda }{{\mu }_{sϵ}\gamma }\}\mathrm{,}\end{array}$

and the set $\Omega$ is a compact positively invariant set for system (1). It is then sufficient to consider solutions in $\Omega$.

3. The Basic Reproduction Number and the Disease-Free Equilibrium

The disease-free equilibrium of system (1) is

${\mathcal{E}}^0=\left(\mathrm{0,0,0,}{X}_s^0\mathrm{,0}\right)=\left(\mathrm{0,0,0,}\frac{\Lambda }{{\mu }_{sϵ}}\mathrm{,0}\right)$. Using the notations of [10] for the model system (1), the matrices F and V for the new infection terms and the remaining transfer terms are, respectively, given by

$F=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& \beta \frac{\Lambda }{{\mu }_{sϵ}}& 0\end{array}\right)$ and $V=\left(\begin{array}{cccc}-k_m& {\mu }_{mϵ}& \rho & 0\\ 0& \rho +{\mu }_{fϵ}& 0& -k_f\\ 0& -\rho & {\mu }_{pϵ}& 0\\ 0& 0& & {\mu }_{sϵ}+{\alpha }_s\end{array}\right)$

The basic reproduction number ${\mathcal{R}}_0$ is equal to the spectral radius of the matrix $F{V}^{-1}$, a simple computation gives:

${\mathcal{R}}_0=\frac{\beta \rho k_f\Lambda }{{\mu }_{sϵ}{\mu }_{pϵ}\left({\mu }_{fϵ}+\rho \right)\left({\mu }_{sϵ}+{\alpha }_s\right)}=\frac{\beta \rho k_f{X}_s^0}{{\mu }_{pϵ}\left({\mu }_{fϵ}+\rho \right)\left({\mu }_{sϵ}+{\alpha }_s\right)}.$

One can remark that there is a mistake in the formula for ${\mathcal{R}}_0$ provided in [6].

The basic reproductive number for system (1) measures the average number of new infections generated by a single infected individual in a completely susceptible population.

As it is well known (see, for instance, [10] ), the local asymptotic stability of the disease-free equilibrium is completely determined by the value of ${\mathcal{R}}_0$ compared to unity, i.e., The disease-free equilibrium ${\mathcal{E}}^0$ of the system (1) is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

Hence ${\mathcal{R}}_0$ determines whether the disease will be prevalent in the given population or will go extinct.

Next, we discuss the global stability of infection-free equilibrium by using suitable Lyapunov function and LaSalle invariance principle for system (1). In recent years, the method of Lyapunov functions has been a popular technique to study global properties of population models. However, it is often difficult to construct suitable Lyapunov functions.

Theorem 3 The disease-free equilibrium ${\mathcal{E}}^0$ of system (1) is globally asymptotically stable (GAS) on the nonnegative orthant ${ℝ}_{+}^5$ whenever ${\mathcal{R}}_0\le 1$.

Proof. We shall use the following notations: $x=\left(X_m,X_f,X_p,X_s,X_i\right)$, and $X_s^0=\frac{\Lambda }{{\mu }_{sϵ}}$. To show the global stability of infection-free equilibrium of system (1), we use the following candidate Lyapunov function:

$V\left(x\right)=\frac{{\mu }_{sϵ}+{\alpha }_s}{k_f}{X}_f+\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{X}_p+{\displaystyle {\int }_{X_s^0}^{X_s}}\frac{X_{\tau }-X_s^0}{X_{\tau }}\text{d}{X}_{\tau }+X_i$ (2)

This function satisfies: $V\left(x\right)\ge 0$ for all $x\in \Omega$, and $V\left(x\right)=0$ if and only if $x=\left(X_m,0,0,X_s^0,0\right)$.

Taking the time derivative of the function V (defined by 2), along the solutions of system (1), we obtain

$\begin{array}{c}\stackrel{\dot{}}{V}=\left(1-\frac{X_s^0}{X_s}\right)\left(\Lambda -{\mu }_{sϵ}{X}_s-\beta X_s{X}_p\right)+\left(\beta X_s{X}_p-\left({\mu }_{sϵ}+{\alpha }_s\right)X_i\right)\\ +\frac{{\mu }_{sϵ}+{\alpha }_s}{k_f}\left(k_f{X}_i-\left({\mu }_{fϵ}+\rho \right)\right)X_f+\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\left(\rho X_f-{\mu }_{pϵ}{X}_p\right)\end{array}$

Using $\Lambda -{\mu }_{sϵ}{X}_s^0=0$, we get

$\begin{array}{c}\stackrel{\dot{}}{V}=\left(1-\frac{X_s^0}{X_s}\right)\left(-{\mu }_{sϵ}{X}_s+{\mu }_{sϵ}{X}_s^0\right)+\left[\beta X_s^0{X}_p-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\mu }_{pϵ}{X}_p\right]\\ ={\mu }_{sϵ}{X}_s^0\left(1-\frac{X_s^0}{X_s}\right)\left(1-\frac{X_s}{X_s^0}\right)+\frac{\beta \Lambda }{{\mu }_{sϵ}}\left[1-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right){\mu }_{mϵ}{\mu }_{pϵ}}{k_f\rho \Lambda \beta }\right]X_p\\ ={\mu }_{sϵ}{X}_s^0\left(1-\frac{X_s^0}{X_s}\right)\left(1-\frac{X_s}{X_s^0}\right)+\frac{\beta \Lambda }{{\mu }_{sϵ}}\left[1-\frac{1}{{\mathcal{R}}_0}\right]X_p\\ =-\frac{{\mu }_{sϵ}}{X_s}{\left(X_s^0-X_s\right)}^2+\frac{\beta \Lambda }{{\mu }_{sϵ}}\left[1-\frac{1}{{\mathcal{R}}_0}\right]X_p\end{array}$ (3)

Hence, $\stackrel{\dot{}}{V}\le 0$ if ${\mathcal{R}}_0\le 1$, and

$\Omega \cap \left\{\stackrel{\dot{}}{V}=0\right\}=\{\begin{array}{l}\left\{x\in \Omega \mathrm{<mo>\; :\; </mo>}x=\left(X_m\mathrm{,}{X}_f\mathrm{,0,}{X}_s^0\mathrm{,}{X}_i\right)\right\}\text{if}{\mathcal{R}}_0<1\hfill \\ \left\{x\in \Omega \mathrm{<mo>\; :\; </mo>}x=\left(X_m\mathrm{,}{X}_f\mathrm{,}{X}_p\mathrm{,}{X}_s^0\mathrm{,}{X}_i\right)\right\}\text{if}{\mathcal{R}}_0=1\hfill \end{array}$

We will show that the largest invariant set $\mathcal{L}$ contained in $\Omega \cap \left\{\stackrel{\dot{}}{V}=0\right\}$ is reduced to the disease-free equilibrium ${\mathcal{E}}^0$.

Let $x=\left(X_m\mathrm{,}{X}_f\mathrm{,}{X}_p\mathrm{,}{X}_s\mathrm{,}{X}_i\right)\in \mathcal{L}$ and

$x\left(t\right)=\left(X_m\left(t\right)\mathrm{,}{X}_f\left(t\right)\mathrm{,}{X}_p\left(t\right)\mathrm{,}{X}_s\left(t\right)\mathrm{,}{X}_i\left(t\right)\right)$ the solution of (1) issued from this point. By invariance of $\mathcal{L}$, we have $X_s\left(t\right)\equiv X_s^0$ which implies ${\stackrel{\dot{}}{X}}_s\left(t\right)=0=\Lambda -{\mu }_s{X}_s\left(t\right)-\beta X_p\left(t\right)X_s\left(t\right)=\Lambda -{\mu }_s{X}_s^0-\beta X_p\left(t\right)X_s^0$ and hence $X_p\left(t\right)=0$ for all t. But, $X_p\left(t\right)\equiv 0$ implies that ${\stackrel{\dot{}}{X}}_p\left(t\right)=0$ for all t which implies, using system (1), that $X_f\left(t\right)=0$ for all t. In the same way, it can be proved that $X_i\left(t\right)=0$ for all t. Reporting in the first equation of system (1), one obtains that, in $\mathcal{L}$,

${\stackrel{\dot{}}{X}}_m\left(t\right)=-{\mu }_{mϵ}{X}_m\left(t\right)\forall t$

Thus the solution of (1) issued from $x=\left(X_m\mathrm{,}{X}_f\mathrm{,}{X}_p\mathrm{,}{X}_s\mathrm{,}{X}_i\right)\in \mathcal{L}$ is given by $x\left(t\right)=\left(X_m{\text{e}}^{-{\mu }_{mϵ}t}\mathrm{,0,0,}{X}_s^0\mathrm{,0}\right)$ which clearly leaves $\Omega$ and hence $\mathcal{L}$ for $t<0$ if $X_m\ne 0$. Therefore $\mathcal{L}=\left\{{\mathcal{E}}^0\right\}$ and hence ${\mathcal{E}}^0$ is a globally asymptotically stable equilibrium state for system (1) on the compact set $\Omega$ thanks to LaSalle invariance principle [11] , (one can also see [12] , Theorem 3.7.11, page 346). Since the set $\Omega$ is an attractive set, the DFE is actually GAS on the nonnegative orthant ${ℝ}_{+}^5$.

Biologically speaking, Theorem 3 implies that schistosomiasis may be eliminated from the community if ${\mathcal{R}}_0\le 1$. One can remark that ${\mathcal{R}}_0$ does not depend on ${\mu }_{mϵ}={\mu }_m+{ϵ}_m$. Hence it is not helpful to try to control the male schistosoma population and then one can take ${ϵ}_m=0$. Therefore the only way to eliminate schistosomiasis is to increase the killing rates of female schistosoma ( ${ϵ}_f$ ), paired schistosoma ( ${ϵ}_p$ ) and snails ( ${ϵ}_s$ ) in order to have ${\mathcal{R}}_0\le 1$.

In the rest of this section, we show that the disease persists when ${\mathcal{R}}_0>1$. The disease is endemic if the infected fraction of the population persists above a certain positive level. The endemicity of a disease can be well captured and analyzed through the notion of uniform persistence. System (1) is said to be uniformly persistent in $\Omega$ if there exists constant $c>0$, independent of initial conditions in $\stackrel{\circ }{\Omega }$ (the interior of $\Omega$ ), such that all solutions $\left(X_m\left(t\right)\mathrm{,}{X}_f\left(t\right)\mathrm{,}{X}_p\left(t\right)\mathrm{,}{X}_s\left(t\right)\mathrm{,}{X}_i\left(t\right)\right)$ of system (1) satysfy

$\underset{t\to \infty }{\lim \inf }{X}_m\left(t\right)\ge c,\underset{t\to \infty }{\lim \inf }{X}_f\left(t\right)\ge c,\underset{t\to \infty }{\lim \inf }{X}_p\left(t\right)\ge c,$

$\underset{t\to \infty }{\lim \inf }{X}_s\left(t\right)>c,\underset{t\to \infty }{\lim \inf }{X}_i\left(t\right)\ge c,$

provided $\left(X_m\left(0\right)\mathrm{,}{X}_f\left(0\right)\mathrm{,}{X}_p\left(0\right)\mathrm{,}{X}_s\left(0\right)\mathrm{,}{X}_i\left(0\right)\right)\in \stackrel{\circ }{\Omega }$, (see [13] [14] ).

Theorem 4 System (1) is uniformly persistent in $\Omega$ if and only if ${\mathcal{R}}_0>1$.

Proof. When ${\mathcal{R}}_0\le 1$, the infection-free equilibrium ${\mathcal{E}}^0$ is globally asymptotically stable which precludes any sort of persistence and hence ${\mathcal{R}}_0>1$ is a necessary condition for persistence. In order to show that ${\mathcal{R}}_0>1$ is a sufficient condition for uniform persistence, it suffices to verify conditions (1) and (2) of Theorem 4.1 in [15] (one can also see [16] , Theorem 3.5).

We use the notations of [15] with $X=\Omega$ and $Y=\partial \Omega$. Let M be the largest invariant compact set in $Y$. We have already seen that $M=\left\{{\mathcal{E}}^0\right\}$, and so M is isolated. To show that ${\mathcal{W}}^s\left(M\right)$ (the stable set of M) is contained in $Y=\partial \Omega$, we use the following function:

$\mathcal{F}=\frac{{\mu }_{sϵ}+{\alpha }_s}{k_f}{X}_f+\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{X}_p+X_i$

The time derivative of $\mathcal{F}$ along the solutions of system (1) is given by

$\begin{array}{c}\stackrel{\dot{}}{\mathcal{F}}=\beta X_s{X}_p-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\mu }_{pϵ}{X}_p\\ =\left(\beta X_s-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\mu }_{pϵ}\right)X_p\\ =\frac{{\mu }_{pϵ}\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\left(\beta X_s\frac{k_f\rho }{{\mu }_{pϵ}\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}-1\right)X_p\\ =\frac{{\mu }_{pϵ}\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\left({\mathcal{R}}_0\frac{X_s}{X_s^0}-1\right)X_p\end{array}$

Since ${\mathcal{R}}_0>1$, we have $\stackrel{\dot{}}{\mathcal{F}}>0$ for $X_p>0$ and $\frac{X_s^0}{{\mathcal{R}}_0}<X_s\le X_s^0$. Therefore $\stackrel{\dot{}}{\mathcal{F}}>0$ in a neighborhood N of ${\mathcal{E}}^0$ relative to $\Omega \backslash \partial \Omega$. This implies that any solution starting in N must leave N at finite time and hence the stable set of M, ${\mathcal{W}}^s\left(M\right)$ is contained in $\partial \Omega$.

4. Endemic Equilibrium and Its Stability

Endemic equilibrium points are steady-state solutions where the disease persists in the population (all state variables are positive).

In this case system (1) has an endemic equilibrium point given by

${\mathcal{E}}_h=\{\begin{array}{l}{X}_m^{\mathrm{<mo>\; *\; </mo>}}={\mu }_{pϵ}{\mu }_{sϵ}\frac{\rho \left(k_m-k_f\right)+k_m{\mu }_{fϵ}}{\beta \rho k_f{\mu }_{mϵ}}\left({\mathcal{R}}_0-1\right)\mathrm{,}\\ X_f^{\mathrm{<mo>\; *\; </mo>}}=\frac{k_f\Lambda }{\left({\mu }_{sϵ}+{\alpha }_s\right)\left(\rho +{\mu }_{fϵ}\right)}\left(1-\frac{1}{{\mathcal{R}}_0}\right)\mathrm{,}{X}_p^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; =\; </mo>}\frac{\left({\mathcal{R}}_0-1\right){\mu }_{sϵ}}{\beta }\mathrm{,}\\ X_s^{\mathrm{<mo>\; *\; </mo>}}=\frac{\Lambda }{{\mathcal{R}}_0{\mu }_{sϵ}}\mathrm{,}{X}_i^{\mathrm{<mo>\; *\; </mo>}}=\frac{\Lambda }{{\mu }_{sϵ}+{\alpha }_s}\left(1-\frac{1}{{\mathcal{R}}_0}\right)\mathrm{.}\end{array}$

This equilibrium has a biological sense only when ${\mathcal{R}}_0>1$.

Theorem 5 If ${\mathcal{R}}_0>1$, the unique endemic equilibrium ${\mathcal{E}}_h$ is globally asymptotically stable.

Proof. In order to investigate the global stability of the endemic equilibrium, we consider the following function defined on

${\Omega }_1=\left\{x\in \Omega \mathrm{<mo>\; :\; </mo>}{X}_f>\mathrm{0,}{X}_p>\mathrm{0,}{X}_s>0\text{and}{X}_s>0\right\}$ :

$\begin{array}{c}W\left(x\right)=\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)}{k_f}{\displaystyle {\int }_{X_f}^{X_f^{*}}}\frac{u-X_f^{*}}{u}\text{d}u+\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\displaystyle {\int }_{X_p}^{X_p^{*}}}\frac{u-X_p^{*}}{u}\text{d}u\\ +{\displaystyle {\int }_{X_s}^{X_s^{*}}}\frac{u-X_s^{*}}{u}\text{d}u+{\displaystyle {\int }_{X_i}^{X_i^{*}}}\frac{u-X_i^{*}}{u}\text{d}u\end{array}$

This function satisfies: $W\left(x\right)\ge 0$ for all $x\in {\Omega }_1$, and $W\left(x\right)=0$ if and only if $\left(X_f\mathrm{,}{X}_p\mathrm{,}{X}_s\mathrm{,}{X}_i\right)=\left(X_f^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_p^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_s^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_i^{\mathrm{<mo>\; *\; </mo>}}\right)$. The time derivative of W with respect to the solutions of system (1) is

$\begin{array}{c}\stackrel{\dot{}}{W}=\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)\left(\Lambda -{\mu }_{sϵ}{X}_s-\beta X_s{X}_p\right)\\ +\left(1-\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\right)\left(\beta X_s{X}_p-\left({\mu }_{sϵ}+{\alpha }_s\right)X_i\right)\\ +\frac{\left({\mu }_{s\epsilon }+{\alpha }_s\right)}{k_f}\left(1-\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}\right)\left(k_f{X}_i-\left({\mu }_{fϵ}+\rho \right)X_f\right)\\ +\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\left(1-\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\right)\left(\rho X_f-{\mu }_{pϵ}{X}_p\right)\end{array}$

$\begin{array}{l}=\left(\Lambda -{\mu }_{sϵ}{X}_s\right)\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)+\beta X_s^{\mathrm{<mo>\; *\; </mo>}}{X}_p-\beta X_s{X}_p-\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\beta X_s{X}_p\\ +\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}+\beta X_s{X}_p-\left({\mu }_{sϵ}+{\alpha }_s\right)X_i-\left({\mu }_{sϵ}+{\alpha }_s\right)\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}{X}_i\\ +\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{X}_f^{\mathrm{<mo>\; *\; </mo>}}-\left({\mu }_{sϵ}+{\alpha }_s\right)X_i-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f}{X}_f\\ -\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\mu }_{pϵ}{X}_p+\left({\mu }_{sϵ}+{\alpha }_s\right)X_i\end{array}$

$\begin{array}{l}-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f}{X}_f-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\mu }_{pϵ}{X}_p\\ +\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}{\mu }_{pϵ}{X}_p+\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f}{X}_f\\ -\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}{X}_f\end{array}$

Thus,

$\begin{array}{c}\stackrel{\dot{}}{W}=\left(\Lambda -{\mu }_{sϵ}{X}_s\right)\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)+\beta X_s^{\mathrm{<mo>\; *\; </mo>}}{X}_p-\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\beta X_s{X}_p+\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\\ -\left({\mu }_{sϵ}+{\alpha }_s\right)\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}{X}_i+\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f}{X}_f^{\mathrm{<mo>\; *\; </mo>}}\\ -\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }{\mu }_{pϵ}{X}_p-\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f}\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}{X}_f\\ +\frac{\left({\mu }_{sϵ}+{\alpha }_s\right)\left({\mu }_{fϵ}+\rho \right)}{k_f\rho }\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}{\mu }_{pϵ}{X}_p\mathrm{.}\end{array}$

Using the equilibrium relations ( E Q )

$\mathcal{E}\mathcal{Q}=\{\begin{array}{l}\Lambda -{\mu }_{sϵ}{X}_s^{\mathrm{<mo>\; *\; </mo>}}=\beta X_s^{\mathrm{<mo>\; *\; </mo>}}{X}_p^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\beta \Lambda -{\mu }_{sϵ}{X}_s^{\mathrm{<mo>\; *\; </mo>}}=\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\\ \rho X_f^{\mathrm{<mo>\; *\; </mo>}}={\mu }_{pϵ}{X}_p^{\mathrm{<mo>\; *\; </mo>}},k_f{X}_i^{\mathrm{<mo>\; *\; </mo>}}=\left({\mu }_{sϵ}+\rho \right)X_f^{\mathrm{<mo>\; *\; </mo>}}\mathrm{.}\end{array}$

It follows that

$\begin{array}{c}\stackrel{\dot{}}{W}=\left(\Lambda -{\mu }_{sϵ}{X}_s\right)\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)+\beta X_s^{\mathrm{<mo>\; *\; </mo>}}{X}_p-\left({\mu }_{sϵ}+{\alpha }_s\right)\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{\beta X_s^{\mathrm{<mo>\; *\; </mo>}}{X}_p^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\beta X_s{X}_p\\ +\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}-\left({\mu }_{sϵ}+{\alpha }_s\right)\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}{X}_i+\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\\ -\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}-\left({\mu }_{sϵ}+{\alpha }_s\right)\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_f^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}{X}_f+\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\\ =\left(\Lambda -{\mu }_{sϵ}{X}_s\right)\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)+\beta X_s^{\mathrm{<mo>\; *\; </mo>}}{X}_p+\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\frac{X_f}{X_f^{\mathrm{<mo>\; *\; </mo>}}}-\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\right]\\ -\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\frac{X_s}{X_s^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}+\frac{X_i}{X_i^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}+\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\frac{X_f}{X_f^{\mathrm{<mo>\; *\; </mo>}}}-2\right]\end{array}$

$\begin{array}{l}=\left(\Lambda -{\mu }_{sϵ}{X}_s\right)\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)+\left({\mu }_{sϵ}+{\alpha }_s\right)\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X}_p+\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}-\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}\right]\\ -\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}+\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\frac{X_s}{X_s^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}+\frac{X_i}{X_i^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}+\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_f}\frac{X_p}{X_f^{\mathrm{<mo>\; *\; </mo>}}}-4\right]\\ +\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}-2\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\mathrm{.}\end{array}$

This implies that

$\begin{array}{c}\stackrel{\dot{}}{W}=\left(\Lambda -{\mu }_{sϵ}{X}_s-\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\right)\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)\\ +\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_p}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}-\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}+\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}-1\right]\\ -\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}+\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\frac{X_s}{X_s^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}+\frac{X_i}{X_i^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}+\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_f}\frac{X_p}{X_f^{\mathrm{<mo>\; *\; </mo>}}}-4\right]\end{array}$

And since $\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}=\Lambda -{\mu }_{sϵ}{X}_s^{\mathrm{<mo>\; *\; </mo>}}$, it follows that

$\begin{array}{c}\stackrel{\dot{}}{W}={\mu }_{sϵ}{X}_s^{\mathrm{<mo>\; *\; </mo>}}\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)\left(1-\frac{X_s}{X_s^{\mathrm{<mo>\; *\; </mo>}}}\right)\\ -\left({\mu }_{sϵ}+{\alpha }_s\right)X_i^{\mathrm{<mo>\; *\; </mo>}}\left[\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}+\frac{X_i^{\mathrm{<mo>\; *\; </mo>}}}{X_i}\frac{X_s}{X_s^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_p}{X_p^{\mathrm{<mo>\; *\; </mo>}}}+\frac{X_i}{X_i^{\mathrm{<mo>\; *\; </mo>}}}\frac{X_f^{\mathrm{<mo>\; *\; </mo>}}}{X_f}+\frac{X_p^{\mathrm{<mo>\; *\; </mo>}}}{X_f}\frac{X_p}{X_f^{\mathrm{<mo>\; *\; </mo>}}}-4\right]\end{array}$

From the AM-GM inequality (which says that the algebraic mean is not smaller than the geometric mean), we have

$\frac{X_s^{*}}{X_s}+\frac{X_i^{*}}{X_i}\frac{X_s}{X_s^{*}}\frac{X_p}{X_p^{*}}+\frac{X_i}{X_i^{*}}\frac{X_f^{*}}{X_f}+\frac{X_p^{*}}{X_f}\frac{X_p}{X_f^{*}}-4\ge 0$

$\left(1-\frac{X_s^{\mathrm{<mo>\; *\; </mo>}}}{X_s}\right)\left(1-\frac{X_s}{X_s^{\mathrm{<mo>\; *\; </mo>}}}\right)\le 0.$

Then, $\stackrel{\dot{}}{W}\le 0$ on ${\Omega }_1$ for ${\mathcal{R}}_0>1$. Hence, W is a Lyapunov function on ${\Omega }_1$. Moreover, $\stackrel{\dot{}}{W}=0$ if and only if $X_f=X_f^{*}$, $X_p=X_p^{*}$, $X_s=X_s^{*}$, and $X_i=X_i^{*}$.

To obtain the largest invariant set $\mathcal{L}$ within the region $\left\{x\in {\Omega }_1:\stackrel{\dot{}}{W}=0\right\}$, note that the trajectory of $X_m\left(t\right)$ with an initial condition in $\mathcal{L}$ must be a solution of:

$\frac{\text{d}{X}_m}{\text{d}t}=k_m{X}_i^{*}-{\mu }_{mϵ}{X}_m-\rho X_f^{*}$

Consequently, we have that

$X_m\left(t\right)=\frac{k_m{X}_i^{\mathrm{<mo>\; *\; </mo>}}-\rho X_f^{\mathrm{<mo>\; *\; </mo>}}}{{\mu }_{mϵ}}+{\text{e}}^{-\left({\mu }_{mϵ}+{ϵ}_m\right)t}\left(X_m\left(0\right)-\frac{k_m{X}_i^{\mathrm{<mo>\; *\; </mo>}}-\rho X_f^{\mathrm{<mo>\; *\; </mo>}}}{{\mu }_{mϵ}}\right)\forall t$

Since $X_m\left(t\right)$ must not leave the domain $\mathcal{L}$ for all t, it follows that

$X_m\left(t\right)=\frac{k_m{X}_i^{\mathrm{<mo>\; *\; </mo>}}-\rho X_f^{\mathrm{<mo>\; *\; </mo>}}}{{\mu }_{mϵ}}\forall t$

Hence, the largest invariant set $\mathcal{L}$ contained in $\left\{x\in {\Omega }_1:\stackrel{\dot{}}{W}=0\right\}$ is reduced to $\left\{{\mathcal{E}}_h\right\}$, and therefore by LaSalle’s principle [11] , ${\mathcal{E}}_h$ is globally asymptotically stable over ${\Omega }_1$.

5. Sensitivity Analysis and Numerical Simulations

Sensitivity analysis and simulations are important to determine how best we can reduce the effect of schistosomiasis, by studying the relative importance of different factors responsible for its transmission and prevalence. Generally speaking, initial disease transmission is directly related to the basic reproduction number, and the disease prevalence is directly related to the endemic equilibrium state ${\mathcal{E}}_h$, and more specifically to the magnitude of $X_i^{\mathrm{<mo>\; *\; </mo>}}$, $X_m^{\mathrm{<mo>\; *\; </mo>}}$, $X_f^{\mathrm{<mo>\; *\; </mo>}}$, $X_p^{\mathrm{<mo>\; *\; </mo>}}$. We perform the analysis by deriving the sensitivity indices of the basic reproduction number to the parameters using both local and global methods.

5.1. Local Sensitivity Analysis of ${\mathcal{R}}_0$

We calculate the sensitivity indices of the reproductive number, ${\mathcal{R}}_0$, and the endemic equilibrium point, ${\mathcal{E}}_h$, to the parameters in the model. we can derive an analytical expression for its sensitivity to each parameter using the normalized forward sensitivity index as described by Chitnis et al. [17].

The normalized forward sensitivity index of a variable to a parameter is a ratio of the relative change in the variable to the relative change in the parameter. When a variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives.

Definition 1 The normalised forward sensitivity index of a variable p that depends differentiable on a parameter q is defined as:

${\Upsilon }_q^p=\frac{\partial p}{\partial q}\times \frac{q}{p}$

Sensitivity analysis is commonly used to determine the robustness of model predictions to parameter values (since there are usually errors in data collection and presumed parameter values). Here we use it to discover parameters that have a high impact on ${\mathcal{R}}_0$, and ${\mathcal{E}}_h$, and should be targeted by intervention strategies.

The sensitivity analysis of ${\mathcal{R}}_0$ has already been done in [6]. We just correct here the expressions of the flexibilities of ${\mu }_{sϵ}$, ${\mu }_{fϵ}$, and ${\mu }_{pϵ}$ on the basic reproduction number ${\mathcal{R}}_0$ (the mistakes in [6] are due to the error in the expression of ${\mathcal{R}}_0$ ). The right expressions are:

${\Upsilon }_{{\mu }_{sϵ}}^{{\mathcal{R}}_0}=\frac{\partial {\mathcal{R}}_0}{\partial {\mu }_{sϵ}}\times \frac{{\mu }_{sϵ}}{{\mathcal{R}}_0}=-\frac{2{\mu }_{sϵ}+{\alpha }_s}{{\mu }_{sϵ}+{\alpha }_s}$

${\Upsilon }_{{\mu }_{fϵ}}^{{\mathcal{R}}_0}=\frac{\partial {\mathcal{R}}_0}{\partial {\mu }_{fϵ}}\times \frac{{\mu }_{fϵ}}{{\mathcal{R}}_0}=-\frac{{\mu }_{fϵ}}{{\mu }_{fϵ}+\rho }$

${\Upsilon }_{{\mu }_{pϵ}}^{{\mathcal{R}}_0}=\frac{\partial {\mathcal{R}}_0}{\partial {\mu }_{pϵ}}\times \frac{{\mu }_{pϵ}}{{\mathcal{R}}_0}=-1$

However the conclusions are not affected: the authors remarked that $\left|{\Upsilon }_{{\mu }_{fϵ}}^{{\mathcal{R}}_0}\right|<1=\left|{\Upsilon }_{{\mu }_{pϵ}}^{{\mathcal{R}}_0}\right|<\left|{\Upsilon }_{{\mu }_{sϵ}}^{{\mathcal{R}}_0}\right|$, and so the most sensitive parameter most sensitive parameter is ${\mu }_{sϵ}$ the death rate of snails, followed by ${\mu }_{pϵ}$ the death rate of pair schistosoma. The least sensitive parameter is ${\mu }_{fϵ}$ the death rate of female single schistosoma. Therefore the most efficient way to reduce the value of ${\mathcal{R}}_0$ is to reduce the snail host population.

Sensitivity analysis of ${\mathcal{E}}_h$

Since in general it is not easy to reduce the value of ${\mathcal{R}}_0$ to be less than one and hence to eradicate the disease, one of the control strategy goal could be to reduce the disease prevalence. To this end, we perform a sensitivity analysis of the endemic equilibrium state. Sensitivity analysis of the endemic equilibrium has usually been used to determine the relative importance of different parameters responsible for equilibrium disease prevalence. Equilibrium disease prevalence is related to the magnitude of $\left(X_m^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_f^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_p^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{X}_i^{\mathrm{<mo>\; *\; </mo>}}\right)$, and specifically to the magnitude of $X_i^{\mathrm{<mo>\; *\; </mo>}}$.

The sensitivity indices of $X_i^{\mathrm{<mo>\; *\; </mo>}}$, to the parameters, ${\mu }_{sϵ}$, ${\mu }_{fϵ}$ and ${\mu }_{pϵ}$ are given by

${\Upsilon }_{{\mu }_{sϵ}}^{X_i^{*}}=\frac{\partial X_i^{*}}{\partial {\mu }_{sϵ}}\times \frac{{\mu }_{sϵ}}{X_i^{*}}=\frac{{\mu }_{sϵ}\left(\beta \rho k_f{\Lambda }_s+{\mu }_{pϵ}\left({\mu }_{fϵ}+\rho \right){\left({\alpha }_s+{\mu }_{sϵ}\right)}^2\right)}{\left({\alpha }_s+{\mu }_{sϵ}\right)\left({\mu }_{pϵ}{\mu }_{sϵ}\left({\mu }_{fϵ}+\rho \right)\left({\alpha }_s+{\mu }_{sϵ}\right)-\beta \rho k_f{\Lambda }_s\right)}$

${\Upsilon }_{{\mu }_{fϵ}}^{X_i^{*}}=\frac{\partial X_i^{*}}{\partial {\mu }_{fϵ}}\times \frac{{\mu }_{fϵ}}{X_i^{*}}=\frac{{\mu }_{sϵ}\left(\beta \rho k_f{\Lambda }_s+{\mu }_{pϵ}\left({\mu }_{fϵ}+\rho \right){\left({\alpha }_s+{\mu }_{sϵ}\right)}^2\right)}{\left({\alpha }_s+{\mu }_{sϵ}\right)\left({\mu }_{pϵ}{\mu }_{sϵ}\left({\mu }_{fϵ}+\rho \right)\left({\alpha }_s+{\mu }_{sϵ}\right)-\beta \rho k_f{\Lambda }_s\right)}$

${\Upsilon }_{{\mu }_{pϵ}}^{X_i^{*}}=\frac{\partial X_i^{*}}{\partial {\mu }_{pϵ}}\times \frac{{\mu }_{pϵ}}{X_i^{*}}=\frac{{\mu }_{pϵ}{\mu }_{sϵ}\left({\mu }_{fϵ}+\rho \right)\left({\alpha }_s+{\mu }_{sϵ}\right)}{{\mu }_{pϵ}{\mu }_{sϵ}\left({\mu }_{fϵ}+\rho \right)\left({\alpha }_s+{\mu }_{sϵ}\right)-\beta \rho k_f{\Lambda }_s}$

It follows that

$\frac{{\Upsilon }_{{\mu }_{sϵ}}^{X_i^{*}}}{{\Upsilon }_{{\mu }_{fϵ}}^{X_i^{*}}}=\frac{\beta \rho k_f{\Lambda }_s}{{\mu }_{fϵ}{\mu }_{pϵ}{\left({\alpha }_s+{\mu }_{sϵ}\right)}^2}+\frac{{\mu }_{fϵ}+\rho }{{\mu }_{fϵ}}$

$\frac{{\Upsilon }_{{\mu }_{fϵ}}^{X_i^{*}}}{{\Upsilon }_{{\mu }_{pϵ}}^{X_i^{*}}}=\frac{{\mu }_{fϵ}}{{\mu }_{fϵ}+\rho }$

$\frac{{\Upsilon }_{{\mu }_{sϵ}}^{X_i^{*}}}{{\Upsilon }_{{\mu }_{pϵ}}^{X_i^{*}}}=\frac{\beta \rho k_f{\Lambda }_s}{{\mu }_{pϵ}\left({\mu }_{fϵ}+\rho \right){\left({\alpha }_s+{\mu }_{sϵ}\right)}^2}+1$

This implies that

$\left|{\Upsilon }_{{\mu }_{fϵ}}^{X_i^{*}}\right|<\left|{\Upsilon }_{{\mu }_{pϵ}}^{X_i^{*}}\right|<\left|{\Upsilon }_{{\mu }_{sϵ}}^{X_i^{*}}\right|$

We note that the most sensitive parameter for $X_i^{\mathrm{<mo>\; *\; </mo>}}$ is ${\mu }_{sϵ}$ the death rate of host snails followed by ${\mu }_{pϵ}$ the death rate of pair parasites and ${\mu }_{fϵ}$ the death rate of female parasites.

5.2. Global Sensitivity Analysis of ${\mathcal{R}}_0$

In this subsection we propose the global sensitivity analysis of the model parameter to determine how much the parameters affect the output of the model. Global sensitivity analysis is a collection of more robust procedures, modifying groups of parameters simultaneously, with a specific goal to recognize the impacts of interactions between various parameters. LHS is at present the most productive and refined statistical techniques [18] and Blower presented it of the field of disease modelling in 1994. We use the technique of Latin Hypercube Sampling, which belong to the monte Carlo class of sampling methods [19]. LHS allows for an efficient analysis of parameter variations across simultaneous uncertainty ranges in each parameter. For each parameter, a probability density function is defined and stratified into N equiproportional serial intervals [20]. Here, for each input parameter we have assumed a uniform distribution across the ranges listed in Table 1 due to the absence of data on the distribution function. We then calculated ${\mathcal{R}}_0$ as the model output using $n=1000$ sets of sampled parameters. We used the partial rank correlation coefficient (PRCC) to assess the significance of each parameter with respect to ${\mathcal{R}}_0$. Figure 1 illustrates the results for the range of parameters in Table 1. The sign of the correlation coefficient indicates the direction of the relationship and the value of the correlation indicates the strength of the relationship between input parameters and model output. The global sensitive analysis confirm the local conclusions by showing the influence of the death rate to the model output ${\mathcal{R}}_0$. The death rate ${\mu }_{sϵ}$ have negative PRCC values, all above 0.5 indicating high significance to ${\mathcal{R}}_0$ with indirect proportional relationship, that is, an increase in ${\mu }_{sϵ}$ increases ${\mathcal{R}}_0$. This suggests that this parameter need to be estimated with precision ignored to accurately capture the transmission dynamics of schistosomiasis. The model output is also sensitive to ${\mu }_{pϵ}$ and ${\mu }_{fϵ}$ with PRCC negative indicating a decrease in ${\mathcal{R}}_0$.

Figure 1. Sensitivity indexes of ${\mathcal{R}}_0$ in terms of the model parameters disposed in order of increasing magnitude.

Table 1. Numerical values of the parameters of model system.

6. Optimal Control

In this section, we aim to place the system (1) thereof in an optimal control setting, in order to be able to calculate the optimal intervention strategies. The optimal control represents the most effective way of controlling the disease that can be adopted by authorities in response to its outbreak. We now modify our model (1) with time-dependent treatment effort $u\left(t\right)$ as control for the system. The variable $u\left(t\right)$ represents the amounts of insecticide that is continuously applied during a considered period, as a measure to fight the disease:

$u\left(t\right)\equiv$ level of insecticide campaigns at time t

Our model with snails treatment can be described with the following differential equations:

$\{\begin{array}{l}\frac{\text{d}{X}_m}{\text{d}t}=k_m{X}_i-\left({\mu }_m+{ϵ}_m\right)X_m-\rho X_f,\\ \frac{\text{d}{X}_f}{\text{d}t}=k_f{X}_i-\left({\mu }_f+{ϵ}_f\right)X_f-\rho X_f,\end{array}$ (4)

$\{\begin{array}{l}\frac{\text{d}{X}_p}{\text{d}t}=\rho X_f-\left({\mu }_p+{ϵ}_p\right)X_p,\\ \frac{\text{d}{X}_s}{\text{d}t}=\Lambda -\left({\mu }_s+{ϵ}_s\right)X_s-\beta X_p{X}_s,\\ \frac{\text{d}{X}_i}{\text{d}t}=\beta X_p{X}_s-\left({\mu }_s+{ϵ}_s+{\alpha }_s+u\left(t\right)\right)X_i.\end{array}$

The control variable $u\left(t\right)$ is a bounded, Lebesgue integrable function that is considered in relative terms, varying from 0 to 1. The goal is to maximize the following objective function

$J\left(u\right)=\underset{u}{\min }{\displaystyle {\int }_0^T}\left(c_s{X}_s\left(t\right)+c_i{X}_i\left(t\right)+\frac{1}{2}{c}_{u}u{\left(t\right)}^2\right)\text{d}t$

subject to the system differential Equations (4), where $c_s$, $c_i$ and $c_u$ are the positive balancing constants. We seek to find an optimal control $u^{\mathrm{<mo>\; *\; </mo>}}$ such that

$J\left(u^{\mathrm{<mo>\; *\; </mo>}}\right)=\underset{u}{\min }\left\{J\left(u\right)\right\}$

where the control set is defined as

$\mathcal{U}=\left\{u\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}T\right]\to \left[\mathrm{0,1}\right]\mathrm{,}u\text{is Lebesgue measurable}\right\}$. Here, the running costs of susceptible snails are given by $c_s{X}_s\left(t\right)$, while term $c_i{X}_i\left(t\right)$ determines the costs of infected snails. Notice that $\frac{1}{2}{c}_{u}u{\left(t\right)}^2$ is the cost of eliminating a fraction $N_Y=\left(X_s\left(t\right)+X_i\left(t\right)\right)$ of the snails population. The choice of the cost function as linear in the number of susceptible and infected and quadratic in the control is as generally done [21] [22] [23].

6.1. The Optimality System

This system satisfies standard conditions for the existence of an optimal control and thus by using Pontryagins Maximum Principle as stated in [24] [25] , we derive the necessary conditions for our optimal control and corresponding states. The Hamiltonian H for the control problem is as follows:

$\begin{array}{c}H=c_s{X}_s\left(t\right)+c_i{X}_i\left(t\right)+\frac{1}{2}{c}_{u}u{\left(t\right)}^2+{\lambda }_1\left(t\right)\frac{\text{d}{X}_m}{\text{d}t}+{\lambda }_2\left(t\right)\frac{\text{d}{X}_f}{\text{d}t}\\ +{\lambda }_3\left(t\right)\frac{\text{d}{X}_p}{\text{d}t}+{\lambda }_4\left(t\right)\frac{\text{d}{X}_s}{\text{d}t}+{\lambda }_5\left(t\right)\frac{\text{d}{X}_i}{\text{d}t}\end{array}$

The adjoint variables ${\lambda }_i\left(i=1,2,3,4,5\right)$ are the solution of the following system:

$\{\begin{array}{l}\frac{\text{d}{\lambda }_1}{\text{d}t}={\mu }_m{\lambda }_1\left(t\right),\\ \frac{\text{d}{\lambda }_2}{\text{d}t}=-{\lambda }_2\left(t\right)\left(-{\mu }_f-\rho \right)+\rho {\lambda }_1\left(t\right)-\rho {\lambda }_3\left(t\right),\\ \frac{\text{d}{\lambda }_3}{\text{d}t}={\mu }_p{\lambda }_3\left(t\right)-\beta {\lambda }_4\left(t\right)X_s\left(t\right)-\beta {\lambda }_5\left(t\right)X_s\left(t\right),\\ \frac{\text{d}{\lambda }_4}{\text{d}t}=-c_s-\beta {\lambda }_4\left(t\right)X_p\left(t\right)-\beta {\lambda }_5\left(t\right)X_p\left(t\right),\\ \frac{\text{d}{\lambda }_5}{\text{d}t}=-c_i-k_f{\lambda }_2\left(t\right)-{\lambda }_4\left(t\right)\left(-{\alpha }_s-{\mu }_s-u\left(t\right)\right)-{\lambda }_5\left(t\right)\left(-{\alpha }_s-{\mu }_s-u\left(t\right)\right).\end{array}$ (5)

with the boundary conditions ${\lambda }_i\left(T\right)=0$.

By using Pontryagins Maximum Principle and the existence result for the optimal control from Fleming and Rishel [8] , we obtain

Theorem 6 There exists an optimal strategy $u^{\mathrm{<mo>\; *\; </mo>}}\in \mathcal{U}$ such that

$J\left(u^{\mathrm{<mo>\; *\; </mo>}}\right)=\underset{u\in \mathcal{U}}{\min }\left\{J\left(u\right)\right\}\mathrm{,}$

given by

$u^{*}=\min \left\{1,\max \left\{0,\frac{{\lambda }_4{X}_i+{\lambda }_5{X}_i}{c_u}\right\}\right\}$

where ${\lambda }_i\left(i=1,2,3,4,5\right)$ are the solutions of (5).

Proof. Here the control and the state variables are nonnegative values. The necessary convexity of the objective functional in u is satisfied for this minimising problem. The control variable set $u\in \mathcal{U}$ is also convex and closed by definition. In addition, the integrand of $J\left(u\right)$ with respect to control variables $u^{\mathrm{<mo>\; *\; </mo>}}$ is convex and it is easy to verify the Lipschitz property of the state system with respect to the state variables. Together with a priori boundedness of the state solutions, the existence of an optimal control has been given by in [8] (see corollary 4.1).

System (5) is obtained by differentiating the Hamiltonian function. Furthermore, by equating to zero the derivatives of the Hamiltonian with respect to the control, we obtain

$\frac{\partial H}{\partial t}=c_u\left(t\right)-{\lambda }_4\left(t\right)X_i\left(t\right)-{\lambda }_5{X}_i\left(t\right)=0$

$\Rightarrow u^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)=\frac{{\lambda }_4\left(t\right)X_i\left(t\right)+{\lambda }_5{X}_i\left(t\right)}{c_u}$

Using the property of the control space, we obtain

$u^{*}\left(t\right)=\{\begin{array}{l}0,\text{if}\frac{l_4\left(t\right)X_i\left(t\right)+l_5{X}_i\left(t\right)}{c_u}\le 0\\ \frac{{\lambda }_4\left(t\right)X_i\left(t\right)+{\lambda }_5{X}_i\left(t\right)}{c_u},\text{if}\frac{l_4\left(t\right)X_i\left(t\right)+l_5{X}_i\left(t\right)}{c_u}\in \left(0,1\right)\\ 1,\text{if}\frac{l_4\left(t\right)X_i\left(t\right)+l_5{X}_i\left(t\right)}{c_u}\ge 1.\end{array}$

Those can be rewritten in compact notation

$u^{*}=\min \left\{1,\max \left\{0,\frac{{\lambda }_4{X}_i+{\lambda }_5{X}_i}{c_u}\right\}\right\}$

6.2. Numerical Examples

The numerical simulations are completed utilizing Matlab and making use of parameter values in [6] to verify the effectiveness of our new model by comparing the disease progression before and after introducing the optimal control variables $u\left(t\right)$. For that, first we solve system (4) with a guess for the controls over the time interval $\left[\mathrm{0,}T\right]$ using a forward fourth-order Runge-Kutta scheme and the transversality conditions ${\lambda }_i\left(T\right)=0,i=1,\cdots ,5$. Then, system (5) is solved by a backward fourth-order Runge-Kutta scheme using the current iteration solution of (4).

The control is updated by using a convex combination of the previous control. The iteration is stopped when the values of the unknowns at the previous iteration are very close to the ones at the present iteration. For more details see, e.g., [25].

We represent the solution curves of the five state variables both in the presence and absence of the control. When viewing the graphs, remember that each of the individuals with control is marked by dashed blue lines. The individuals without control are marked by red lines. It is observed that the application of optimal control reduces a quite larger number of schistosoma (male, female, pair) and snails in the absence of the control. This is occurring as the application of pesticide control reduces the snails population significantly as seen in Figure 2 and Figure 3. Again from the Figures 4-6 it is easy to see that the schistosoma population also much affected due to the use of the insecticide control.

We have considered the schistosomiasis infection in an endemic population (when ${\mathcal{R}}_0$ ). In Figures 4-6, we observe that the fraction of schistosoma (male, female, pair) is lower when control is considered. More precisely, at the end of 15 years, the total number of male, female and pair schistosoma is 10⁵, $0.1\times {10}^5$

Figure 2. The evolution of susceptible snails with and without control. The state is solved forward time with initial conditions $IC\left(0\right)=\left(50000;30000;25000;4500;2500\right)$ while the adjoint system is solved backward in time (in years) with terminal condition $TC=\left(0;0;0;0;0;0\right)$ where $T=5$ and $\left(c_s,c_i,c_u\right)=\left(1,1,0.5\right)$.

Figure 3. The evolution of infected snails with and without control. The state is solved forward time with initial conditions $IC\left(0\right)=\left(50000;30000;25000;4500;2500\right)$ while the adjoint system is solved backward in time in (in years) with terminal condition $TC=\left(0;0;0;0;0;0\right)$ where $T=15$ and $\left(c_s,c_i,c_u\right)=\left(1,1,0.5\right)$.

Figure 4. The evolution of male schistosoma with and without control. The state is solved forward time with initial conditions $IC\left(0\right)=\left(50000;30000;25000;4500;2500\right)$ while the adjoint system is solved backward in time (in years) with terminal condition $TC=\left(0;0;0;0;0;0\right)$ where $T=25$ and $\left(c_s,c_i,c_u\right)=\left(1,1,0.5\right)$.

Figure 5. The evolution of female schistosoma with and without control. The state is solved forward time with initial conditions $IC\left(0\right)=\left(50000;30000;25000;4500;2500\right)$ while the adjoint system is solved backward in time (in years) with terminal condition $TC=\left(0;0;0;0;0;0\right)$ where $T=25$ and $\left(c_s,c_i,c_u\right)=\left(1,1,0.5\right)$.

Figure 6. The evolution of pair schistosoma with and without control. The state is solved forward time with initial conditions IC(0) = (50000; 30000; 25000; 4500; 2500) while the adjoint system is solved backward in time (in years) with terminal condition TC = (0; 0; 0; 0; 0; 0) where T = 25 and (c_s, c_i, c_u) = (1, 1, 0.5).

and $2\times {10}^5$ respectively when control is considered, and $3.5\times {10}^5$, $0.5\times {10}^5$ and $8\times {10}^5$ respectively without control. The schistosoma can survive for very long periods in a dry state, often for more than a few years, that’s why in the stage without control evolution of the number of schistosoma (male, female, pair) take higher levels, once controlled, the number of schistosoma doesn’t develop as shown in Figures 4-6.

In Figure 2, we remark that the number of susceptible snails is more important than in the case without control, this is due to the aim of our approach which focuses on the reduction of the number of susceptible and infected snails population, so in the case without control the number of susceptible snails decrease is less and goes to its stable state, because it also applied a control on schistosoma by the elimination rates ( ${ϵ}_m$, ${ϵ}_f$, ${ϵ}_p$ ).

In Figure 3, it’s the evolution of infected snails which is presented, the major case of our approach, and the graph below shows the effectiveness of the study done in this work. As given above, the numerical simulations suggested 10 years as minimal duration for treatment. We see that if there are control the infected snails population begins to sharply decrease from the very beginning day of treatment and gradually decreases as time goes on.

Acknowledgements

This research was carried out with financial support of CEA-MITIC for postdoc project in Université Gaston Berger de SAINT-LOUIS.

Conflicts of Interest

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

References

[1]	Schistosomiasis. http://www.who.int/mediacentre/factsheets/fs115/en/
[2]	Savioli, L., Stansfield, S., Bundy, D.A., Mitchell, A., Bhatia, R., Engels, D., Montresor, A., Neira, M. and Shein, A.M. (2002) Schistosomiasis and Soil-Transmitted Helminth Infections: Forging Control Efforts. Transactions of the Royal Society of Tropical Medicine and Hygiene, 96, 577-579. https://doi.org/10.1016/S0035-9203(02)90316-0
[3]	Castillo-Chavez, C., Feng, Z. and Xu, D. (2008) A Schistosomiasis Model with Mating Structure and Time Delay. Mathematical Biosciences, 211, 333-341. https://doi.org/10.1016/j.mbs.2007.11.001
[4]	Hadeler, K., Waldstätter, R. and Wörz-Busekros, A. (1988) Models for Pair Formation in Bisexual Populations. Journal of Mathematical Biology, 26, 635-649. https://doi.org/10.1007/BF00276145
[5]	Schmitz, S.-F.H. and Castillo-Chavez, C. (2000) A Note on Pair-Formation Functions. Mathematical and Computer Modelling, 31, 83-91. https://doi.org/10.1016/S0895-7177(00)00025-X
[6]	Qi, L. and Cui, J.-A. (2013) A Schistosomiasis Model with Mating Structure. Abstract and Applied Analysis, 2013, Article ID: 741386. https://doi.org/10.1155/2013/741386
[7]	Xu, D., Curtis, J., Feng, Z. and Minchella, D.J. (2005) On the Role of Schistosome mating Structure in the Maintenance of Drug Resistant Strains. Bulletin of Mathematical Biology, 67, 1207-1226. https://doi.org/10.1016/j.bulm.2005.01.007
[8]	Fleming, W.H. and Rishel, R.W. (2012) Deterministic and Stochastic Optimal Control. Springer Science & Business Media, Berlin, Heidelberg.
[9]	Gamkrelidze, R., Pontrjagin, L.S. and Boltjanskij, V.G. (1964) The Mathematical Theory of Optimal Processes. Macmillan Company, New York.
[10]	Van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[11]	LaSalle, J. (1976) The Stability of Dynamical Systems, Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, PA.
[12]	Bhatia, N.P. and Szegö, G.P. (1967) Dynamical Systems: Stability Theory and Applications. In: Lecture Notes in Mathematics, Springer, Berlin, Heidelberg, New York. https://doi.org/10.1007/BFb0080630
[13]	Thieme, H.R. (1992) Epidemic and Demographic Interaction in the Spread of Potentially Fatal Diseases in Growing Populations. Mathematical Biosciences, 111, 99-130. https://doi.org/10.1016/0025-5564(92)90081-7
[14]	Butler, G., Freedman, H. and Waltman, P. (1986) Uniformly Persistent Systems. Proceedings of the American Mathematical Society, 96, 425-430. https://doi.org/10.2307/2046588
[15]	Hofbauer, J. and So, J.W.-H. (1989) Uniform Persistence and Repellors for Maps. Proceedings of the American Mathematical Society, 107, 1137-1142. https://doi.org/10.2307/2047679
[16]	Lin, X. and So, J.W.-H. (1993) Global Stability of the Endemic Equilibrium and Uniform Persistence in Epidemic Models with Subpopulations. The ANZIAM Journal, 34, 282-295. https://doi.org/10.1017/S0334270000008900
[17]	Chitnis, N., Hyman, J.M. and Cushing, J.M. (2008) Determining Important Parameters in the Spread of Malaria through the Sensitivity Analysis of a Mathematical Model. Bulletin of Mathematical Biology, 70, 1272-1296. https://doi.org/10.1007/s11538-008-9299-0
[18]	Blower, S.M. and Dowlatabadi, H. (1994) Sensitivity and Uncertainty Analysis of Complex Models of Disease Transmission: An HIV Model, as an Example. International Statistical Review/Revue Internationale de Statistique, 66, 229-243. https://doi.org/10.2307/1403510
[19]	Marino, S., Hogue, I.B., Ray, C.J. and Kirschner, D.E. (2008) A Methodology for Performing Global Uncertainty and Sensitivity Analysis in Systems Biology. Journal of Theoretical Biology, 254, 178-196. https://doi.org/10.1016/j.jtbi.2008.04.011
[20]	Wu, J., Dhingra, R., Gambhir, M. and Remais, J.V. (2013) Sensitivity Analysis of Infectious Disease Models: Methods, Advances and Their Application. Journal of the Royal Society Interface, 10, Article ID: 20121018. https://doi.org/10.1098/rsif.2012.1018
[21]	Lee, S., Chowell, G. and Castillo-Chávez, C. (2010) Optimal Control for Pandemic Influenza: The Role of Limited Antiviral Treatment and Isolation. Journal of Theoretical Biology, 265, 136-150. https://doi.org/10.1016/j.jtbi.2010.04.003
[22]	Okosun, K.O., Ouifki, R. and Marcus, N. (2011) Optimal Control Analysis of a Malaria Disease Transmission Model that Includes Treatment and Vaccination with Waning Immunity. Biosystems, 106, 136-145. https://doi.org/10.1016/j.biosystems.2011.07.006
[23]	Tchuenche, J., Khamis, S., Agusto, F. and Mpeshe, S. (2011) Optimal Control and Sensitivity Analysis of an Influenza Model with Treatment and Vaccination. Acta biotheoretica, 59, 1-28. https://doi.org/10.1007/s10441-010-9095-8
[24]	Pontryagin, L.S., Boltyanskii, V., Gamkrelidze, R. and Mishchenko, E.F. (1962) the Mathematical Theory of Optimal Processes. Interscience, New York.
[25]	Lenhart, S. and Workman, J.T. (2007) Optimal Control Applied to Biological Models. CRC Press, Boca Raton, FL.