A Schistosomiasis Model with Mating Structure

Abstract

A schistosomiasis model is proposed including single schistosomes, paired schistosomes, snails, and the latent period of infected snails. A reasonable sex ratio of schistosomes and the minimum mating function are considered. A threshold condition determining the stability of the system is given, and the stability of equilibrium for the model is shown. The impact of the latent period of infected snails on schistosomiasis transmission can be found through numerical simulations. Finally, preferable control strategies are obtained by sensitivity analyses. Killing snails may be the preferred control measure. If we choose chemotherapy, we should use some drugs which are sufficient for reducing egg-associated pathology, since paired schistosomes are mostly harmful to definitive hosts.

1. Introduction

Schistosomiasis is one of the most prevalent parasite diseases. Its transmission needs two hosts: the definitive hosts and the intermediate snail hosts. In definitive hosts, schistosoma has two distinct sexes. Female and male schistosomes pair up they lay eggs that work their way through the intestinal lining or into the ureters and bladder. The eggs themselves are harmless, but in definitive host the body's immune response results in inflammation and eventual damage, as granulomas—fibrous lesions—form around the eggs. These eggs pass in urine or feces into fresh water. In infected water, miracidium hatches from egg and infects the intermediate snail hosts. In snail hosts, parasites asexually reproduce. As a parasite in the body of the snail, the miracidium loses its cilia and undergoes a change into a sac-like primary sporocysts [1]. Many secondary sporocysts are formed in the body of the primary sporocysts. They develop for some time and then rupture the wall of the primary sporocyst and become established in the body of the snail. From the walls of the secondary sporocyst a large number of parthenogenetic eggs are formed, and they develop internally into forms known as cercariae. After about four, weeks cercariae rupture the body-wall of the secondary sporocyst, bore out of the snail's body, and enter the water as very active organisms, characterized by two large suckers and a forked tail. In this stage, they are ready to enter a mammalian host. When such a free-living cercaria comes in contact with a susceptible mammal, it penetrate, the skin of the definitive host and transforms into single schistosoma. Here, the parasite's life cycle is completed.

During the asexual phase in the host mollusc, several groups have observed that the miracidia are already sexually differentiated [2]. Furthermore, male miracidia will only give rise to male later forms, whilst female miracidia will only give rise to female later forms. However, no apparent differences have been registered between physiological alterations in snails infected by a male miracidium and those in snails infected by a female miracidium. There are no distinctive qualitative and quantitative differences in snail host response to infection with a single miracidium and in snail host reaction to infection with numerous miracidia secured by single or multiple exposures [3]. From the above statement, we know that the mating of male schistosoma and female schistosoma in definitive hosts is very important to the transmission of schistosomiasis. Hence, we differentiate female and male parasites only in definitive hosts in this paper.

In the literature, two-sex problems have been studied by many authors [4–6]. All of these models are established for human population diseases, such as sexually transmitted diseases. They constructed pair-formation models and studied the existence, uniqueness and the stabilities of exponential solutions. In [5, 6] authors put forward three forms of pair-formation functions which are also called mating functions: the harmonic mean function, the geometric mean function and the minimum function. All of these three functions can be applied to two-sex models.

As for schistosomiasis, there are few papers considering two-sex problems and mating interactions [4, 7–9]. In [9] Xu et al. have proposed a multi-strain schistosome model with mating structure. Their goal is to infer the impact of drug treatment on the maintenance of schistosome genetic diversity. However, the model of Xu et al. ignores several stages of the parasite's complex life cycle. They consider only the adult parasite populations. Furthermore, in their model they assume that the recruitment rate of single adult parasites at time 𝑡 depends instantaneously on the total number of parasite pairs at time 𝑡. Subsequently, Castillo-Chavez et al. [4] have relaxed this assumption on the parasite's life history by introducing a time delay. The time delay accounts for the average time that must elapse between two adult generations. Their model is in the following: 𝑑𝑚𝑑𝑡=𝑘𝑆_𝑚(𝜏)𝑝(𝑡−𝜏)−𝜇_𝑚+𝜎𝜃𝑚(𝑡)−𝜑𝑚(𝑡),𝑓(𝑡),𝑑𝑓𝑑𝑡=𝑘𝑆_𝑓(𝜏)𝑝(𝑡−𝜏)−𝜇_𝑓+𝜎𝜃𝑓(𝑡)−𝜑𝑚(𝑡),𝑓(𝑡),𝑑𝑝𝑚(𝑡)(𝑡−𝜇𝑑𝑡=𝜑,𝑓)_𝑝+𝜎𝜃𝑝(𝑡). (1) Here 𝑚(𝑡) and 𝑓(𝑡) represent the densities of single male schistosoma and single female schistosoma, respectively, and 𝑝(𝑡) is the density of pairs. The formation of schistosome pairs is described by a mating function 𝜑(𝑚(𝑡),𝑓(𝑡)). 𝑘 is half of the birth rate of a pair per capita. Note that the ratio of male to female offspring is assumed to be 1 : 1 in their model. 𝜇_𝑚, 𝜇_𝑓, and 𝜇_𝑝 are the per capita death rates of single male schistosoma, single female schistosoma, and schistosoma in a mated pair, respectively. 𝜎 is the rate of the chemotherapeutic treatment. 𝜃 denotes the drug resistance of that parasite strain. 𝑆_𝑚(𝜏) and 𝑆_𝑓(𝜏) are functions which keep track of the parasite's survival probabilities in various (nonadult) stages of the life cycle. The dynamical behaviors of their model in [4] are not qualitatively different from those derived from an earlier model in [9], although [9] ignores the impact of time delays associated with the multiple stages in parasite's life cycle. The results of [4] imply that higher treatment rate can allow for coexistence between susceptible and resistant parasite strains. Both in [4, 9], the harmonic mean function is chosen as the mating function 𝜑(𝑚(𝑡),𝑓(𝑡)).

As we know, in the two models of [4, 9], the ratio of male to female offspring is assumed to be 1 : 1. In fact, some experiments have shown that the number of male schistosomes is bigger than that of female [10–16]. In [10] authors found the phenomena of the natural male bias by experimental observations on the sex ratio of adult schistosoma mansoni. In [12] the sex ratio of female schistosoma to male schistosoma is 1 : 1.81. In [16], authors reported that the male-biased sex ratio is 2.36 males to 1 female in Kenya. In two experiments of [14], the sex ratios of female schistosoma to male schistosoma are (0.89±0.03) : 1 and (0.98±0.03) : 1, respectively. These results tend to male-biased sex ratio. In [15] the sex ratio of female schistosoma to male schistosoma is 1 : 1.38. In artificial experiments of [13], there are 9543 males and 6597 females in 16140 mature schistosoma japonicum; the ratio of female to male is 1 : 1.45. In [11] authors also support male-biased sex ratio. All of these results imply that the sex ratio is not 1 : 1. Furthermore, Mao [13] also denote that the assumption of the 50 : 50 sex ratio in previous schistosomiasis models does not accord with the result of experiments. Hence, in this paper, we consider reasonably that the number of male offspring is bigger than that of female in our model.

On the other hand, the two models in [4, 9] do not include the snail dynamics. In reality, from the life cycle of schistosoma, it is easy to see that the parasite offspring is produced directly by infected snails but not by paired parasites. Nåsell [17] also thought that the number of male and female schistosoma is proportional to the number of infected snails. This implies that the snail dynamics may influence the transmission of schistosomiasis. Hence, it is necessary to add snails to the model. In addition, from the life history of schistosoma, we know that the production of parasite offspring needs about four weeks after the snail hosts are infected by miracidia [4, 9]. Motivated by [4, 9, 13, 17] we established a new model including a more reasonable sex ratio of schistosoma, snail dynamics, the latent period of infected snails, and the mating structure of schistosoma. Here, we take the minimum function as mating function. In this paper, our purpose is to study the dynamics of single schistosoma, pairs, and snails and to put forward preferable control strategies.

Our paper is organized as follows. In Section 2, we establish a mathematical model with single schistosoma, pairs and snails. For convenience of analysis, we first study the model without time delay and define the basic reproductive number. And then the stabilities of the disease-free and endemic equilibria are obtained in Section 3. In Section 4, the stability of the delayed model is investigated. In Section 5, computational simulations and sensitivity analysis are performed, and we give preferable control strategies.

2. Mathematical Model

Considering the mating structure of parasites and snail dynamics, we propose a model with state variables 𝑚(𝑡), 𝑓(𝑡), 𝑝(𝑡), 𝑠(𝑡), and 𝑖(𝑡), where 𝑚(𝑡), 𝑓(𝑡), and 𝑝(𝑡) represent the densities of single male schistosoma, single female schistosoma and pairs, respectively; 𝑠(𝑡) and 𝑖(𝑡) represent the number of susceptible and infectious snail hosts, respectively. The basic assumptions are as follows.

(i) Since the single male and female schistosoma are from the cercaria produced by infected snails, in [17, 18] authors consider that the number of male and female schistosoma is proportional to the number of infected snails. According to [17, 18] we also assume that the recruitment rate of single male and female schistosoma is proportional to the number of infected snail hosts, that is, the recruitment rates of single male and female worms are 𝑘_𝑚𝑖(𝑡) and 𝑘_𝑓𝑖(𝑡), respectively. Following papers [10–16], we assume 𝑘_𝑚≥𝑘_𝑓.

(ii) Natural death rates for single male, single female, and pair worms per capita are 𝜇_𝑚, 𝜇_𝑓, and 𝜇_𝑝, respectively. In reality, pairs of schistosoma may live for a few years while single parasites may only live for a few weeks [9]. Therefore, we assume 𝜇_𝑝<𝜇_𝑚,𝜇_𝑓. As for the investigation of the death rate of single male and female schistosoma, there are many results. In the experiment of [12], the livability of single male parasites is stronger than that of single female parasites, which implies that single female parasites have a higher natural death rate than single male parasites. In [19], the author found that single male schistosoma is larger more muscular. May and Woolhouse [20] and Tchuente et al. [21] also found that the growth of single female schistosoma is restrained when it is not paired up. Standen [22] deduced that single female schistosoma is incapable of going against the blood stream but single male schistosoma can. These results imply that the livability of single male parasites is stronger than that of single female parasites. Cornford and Huot [23] and Cornford and Fitzpatrick [24, 25] had shown that single female schistosoma has less glucose than single male schistosoma. B. G. Atkinson and K. H. Atkinson [26] and Davis et al. [27] found that single female schistosoma has less actin than single male schistosoma. These indicate that single female parasite is not easy to survive. In addition, Mao [13] also reported that the survival rate of single male schistosoma is not lower than that of single female schistosoma. All of the above results imply 𝜇_𝑚≤𝜇_𝑓. Hence, we assume 𝜇_𝑝<𝜇_𝑚≤𝜇_𝑓.

(iii) The mating function 𝜑(𝑚(𝑡),𝑓(𝑡)) takes the minimum function 𝜑(𝑚(𝑡),𝑓(𝑡))=𝜌min(𝑚(𝑡),𝑓(𝑡)) in which 𝜌 represents the effective mating rate [5, 6].

(iv) The parameter Λ_𝑠 is the recruitment rate of snail hosts. 𝜇_𝑠 is per capita natural death rate of snail hosts. 𝛼_𝑠 is the disease-induced death rate of snail hosts. The transmission rate from pairs to susceptible snails is a constant 𝜉.

(v) 𝜏 is the latent period of infected snails, and 𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏} represents the survival rate of infected snails in the latent period.

(vi) We suppose that infected snails do not recover from schistosomiasis as their life spans are short.

(vii) All parameters are assumed nonnegative in reality.

Then, we have a model with the following form:𝑑𝑚𝑑𝑡=𝑘_𝑚𝑖(𝑡−𝜏)𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}−𝜇_𝑚𝑚(𝑡)−𝜑𝑚(𝑡),𝑓(𝑡),𝑑𝑓𝑑𝑡=𝑘_𝑓𝑖(𝑡−𝜏)𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}−𝜇_𝑓𝑓(𝑡)−𝜑𝑚(𝑡),𝑓(𝑡),𝑑𝑝𝑚(𝑡)(𝑡𝑑𝑡=𝜑,𝑓)−𝜇_𝑝𝑝(𝑡),𝑑𝑠𝑑𝑡=Λ_𝑠−𝜇_𝑠𝑠(𝑡)−𝜉𝑝(𝑡)𝑠(𝑡),𝑑𝑖(𝑡)𝑠(𝑡)−𝜇𝑑𝑡=𝜉𝑝_𝑠+𝛼_𝑠𝑖(𝑡). (2) For convenience of the stability analysis of (2), we first study the case when 𝜏=0, that is, 𝑑𝑚𝑑𝑡=𝑘_𝑚𝑖(𝑡)−𝜇_𝑚𝑚(𝑡)−𝜑𝑚(𝑡),𝑓(𝑡),𝑑𝑓𝑑𝑡=𝑘_𝑓𝑖(𝑡)−𝜇_𝑓𝑓(𝑡)−𝜑𝑚(𝑡),𝑓(𝑡),𝑑𝑝𝑚(𝑡)(𝑡𝑑𝑡=𝜑,𝑓)−𝜇_𝑝𝑝(𝑡),𝑑𝑠𝑑𝑡=Λ_𝑠−𝜇_𝑠𝑠(𝑡)−𝜉𝑝(𝑡)𝑠(𝑡),𝑑𝑖(𝑡)𝑠(𝑡)−𝜇𝑑𝑡=𝜉𝑝_𝑠+𝛼_𝑠𝑖(𝑡). (3)

In the current paper, we take the mating function 𝜑(𝑚(𝑡),𝑓(𝑡))=𝜌min(𝑚(𝑡),𝑓(𝑡)) [5, 6]. The function is linear on each of the following two sets: 𝐾_𝑓=𝑚,𝑓,𝑝,𝑠,𝑖∈𝑅⁵Λ∣𝑓≥𝑚≥0,𝑝≥0,𝑠≥0,𝑖≥0,0≤𝑠+𝑖≤^𝑠𝜇_𝑠,𝐾_𝑚=𝑚,𝑓,𝑝,𝑠,𝑖∈𝑅⁵Λ∣𝑚≥𝑓≥0,𝑝≥0,𝑠≥0,𝑖≥0,0≤𝑠+𝑖≤^𝑠𝜇_𝑠. (4) In general, none of these sets is positively invariant for (3) with 𝜑(𝑚(𝑡),𝑓(𝑡))=𝜌min(𝑚(𝑡),𝑓(𝑡)). But under the conditions 𝑘_𝑚≥𝑘_𝑓,𝜇_𝑝<𝜇_𝑚≤𝜇_𝑓, (5)

the set 𝐾_𝑚 is positively invariant. Thus, the model (3) with an initial surplus of males becomes the so-called female dominance model on the set 𝐾_𝑚 and the system (3) can be rewritten in the following form: 𝑑𝑚𝑑𝑡=𝑘_𝑚𝑖(𝑡)−𝜇_𝑚𝑚(𝑡)−𝜌𝑓(𝑡),𝑑𝑓𝑑𝑡=𝑘_𝑓𝑖(𝑡)−𝜇_𝑓𝑓(𝑡)−𝜌𝑓(𝑡),𝑑𝑝(𝑡)𝑑𝑡=𝜌𝑓−𝜇_𝑝𝑝(𝑡),𝑑𝑠𝑑𝑡=Λ_𝑠−𝜇_𝑠𝑠(𝑡)−𝜉𝑝(𝑡)𝑠(𝑡),𝑑𝑖(𝑡)𝑠(𝑡)−𝜇𝑑𝑡=𝜉𝑝_𝑠+𝛼_𝑠𝑖(𝑡). (6)

Using the standard method, it is easy to see that the disease-free equilibrium 𝐸₀=(0,0,0,Λ_𝑠/𝜇_𝑠,0) always exists. According to the calculational method of the basic reproductive number [28], the basic reproduction number for model (6) is 𝑅₀=³𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓.+𝜌 (7)

The following section shows that the basic reproductive number 𝑅₀ provides a threshold condition for parasite extinction in (6).

3. Stability Analysis of the System (6)

In this section, we will analyze the stability of model (6). The stability of the disease-free equilibrium determines whether the disease will be prevalent in an uninfected population. The following result shows that the parasites will go extinct if 𝑅₀<1.

Theorem 1. The disease-free equilibrium 𝐸₀ of the system (6) is locally asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

Proof. The Jacobian matrix for system (6) is given by 𝐽=−𝜇_𝑚−𝜌00𝑘_𝑚0−𝜇_𝑓+𝜌00𝑘_𝑓0𝜌−𝜇_𝑝(𝑡)0000−𝜉𝑠−𝜇_𝑠(𝑡)0(𝑡)(𝑡)−𝜇−𝜉𝑝00𝜉𝑠𝜉𝑝_𝑠+𝛼_𝑠. (8) Thus, the eigenvalues of 𝐸₀ are −𝜇_𝑠, −𝜇_𝑚, and the roots of the equation are 𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃=0, (9)
where 𝑎₁=𝜇_𝑠+𝛼_𝑠+𝜇_𝑝+𝜇_𝑓𝑎+𝜌>0,₂=𝜇_𝑠+𝛼_𝑠𝜇_𝑝+𝜇_𝑠+𝛼_𝑠𝜇_𝑓+𝜌+𝜇_𝑝𝜇_𝑓𝑎+𝜌>0,₃=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌−𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠. (10)
Note that 𝑎₃>0 if and only if 𝑅₀<1.
By the Routh-Hurwitz criterion, we know 𝐻₁=𝑎₁𝐻>0,₂=𝑎₁𝑎₃1𝑎₂𝐻>0,₃=𝑎₁𝑎₃01𝑎₂00𝑎₁𝑎₃=𝑎₃𝐻₂, (11) which implies that all the eigenvalues of 𝐸₀ have negative real parts if and only if 𝑅₀<1. It follows that the disease-free equilibrium 𝐸₀ is locally asymptotically stable if 𝑅₀<1 and unstable if 𝑅₀>1.

Note that if 𝑅₀>1, (6) has a unique endemic equilibrium 𝐸^∗=(𝑚^∗,𝑓^∗,𝑝^∗,𝑠^∗,𝑖^∗), where 𝑠^∗=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌𝜉𝜌𝑘_𝑓,𝑝^∗=Λ_𝑠𝜌𝑘_𝑓𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓−𝜇+𝜌𝜉^𝑠=𝜇_𝑠𝜉𝑅³₀,𝑓−1_∗=𝜇_𝑝𝜌𝑝^∗,𝑖^∗=𝜇_𝑝𝜇_𝑓+𝜌𝜌𝑘_𝑓𝑝^∗,𝑚^∗=𝑘_𝑚𝜇_𝑓+𝜌−𝜌𝑘_𝑓𝑘_𝑓𝜇_𝑚𝜇_𝑝𝜌𝑝^∗. (12)

Similarly, using the Routh-Hurtwitz criterion we can obtain the stability of the endemic equilibrium 𝐸^∗.

Theorem 2. The endemic equilibrium 𝐸^∗ of the system (6) is locally asymptotically stable if 𝑅₀>1.

4. Stability Analysis of the System (2)

In this section, we analyze the stability of model (2). We recall it with the following form: 𝑑𝑚𝑑𝑡=𝑘_𝑚𝑖(𝑡−𝜏)𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}−𝜇_𝑚𝑚(𝑡)−𝜌𝑓(𝑡),𝑑𝑓𝑑𝑡=𝑘_𝑓𝑖(𝑡−𝜏)𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}−𝜇_𝑓𝑓(𝑡)−𝜌𝑓(𝑡),𝑑𝑝(𝑡)𝑑𝑡=𝜌𝑓−𝜇_𝑝𝑝(𝑡),𝑑𝑠𝑑𝑡=Λ_𝑠−𝜇_𝑠𝑠(𝑡)−𝜉𝑝(𝑡)𝑠(𝑡),𝑑𝑖(𝑡)𝑠(𝑡)−𝜇𝑑𝑡=𝜉𝑝_𝑠+𝛼_𝑠𝑖(𝑡). (13) Define 𝑅₀=𝜉𝜌Λ_𝑠𝑘_𝑓𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝜇_𝑠𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓.+𝜌 (14)

We can obtain that the disease-free equilibrium 𝐸₀=(0,0,0,Λ_𝑠/𝜇_𝑠,0) always exists and if 𝑅₀>1 there exists the endemic equilibrium 𝐸^∗=(𝑚^∗,𝑓^∗,𝑝^∗,𝑠^∗̃,𝑖^∗̃), where 𝑠^∗̃=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌𝜉𝜌𝑘_𝑓𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏},𝑝^∗=Λ_𝑠𝜌𝑘_𝑓𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓−𝜇+𝜌𝜉^𝑠=𝜇_𝑠𝜉𝑅₀,𝑓−1_∗=𝜇_𝑝𝜌𝑝^∗,𝑖^∗̃=𝜇_𝑝𝜇_𝑓+𝜌𝜌𝑘_𝑓𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝑝^∗,𝑚^∗=𝑘_𝑚𝜇_𝑓+𝜌−𝜌𝑘_𝑓𝑘_𝑓𝜇_𝑚𝜇_𝑝𝜌𝑝^∗. (15)

Theorem 3. If 𝑅₀<1, the disease-free equilibrium 𝐸₀ of the system (13) is locally asymptotically stable.

Proof. The Jacobian matrix for system (13) is given by ||=−𝜇𝐽−𝜆𝐼^𝑚+𝜎𝜃−𝜌00𝑘_𝑚𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝑒^−𝜆0−𝜇_𝑓𝜎+𝜌+𝜃00𝑘_𝑓𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝑒^−𝜆0𝜌−𝜇_𝑝+𝜎𝜃(𝑡)0000−𝜉𝑠−𝜇_𝑠(𝑡)0(𝑡)(𝑡)−𝜇−𝜉𝑝00𝜉𝑠𝜉𝑝_𝑠+𝛼_𝑠. (16)Thus, the eigenvalues of 𝐸₀ are −𝜇_𝑠, −𝜇_𝑚 and the roots of the equation are 𝜆³+𝑎₁𝜆²+𝑎₂𝜆+𝑎₃=0, (17)
where 𝑎₁=𝜇_𝑠+𝛼_𝑠+𝜇_𝑝+𝜇_𝑓𝑎+𝜌>0,₂=𝜇_𝑠+𝛼_𝑠𝜇_𝑝+𝜇_𝑠+𝛼_𝑠𝜇_𝑓+𝜌+𝜇_𝑝𝜇_𝑓𝑎+𝜌>0,₃=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌−𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝑒^−𝜆. (18) If 𝜆≥0 and 𝜏>0,  0<𝑒^−𝜆≤1, and then 𝑎₃≥𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌−𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}. (19)
Note that 𝑎₃>0 if 𝑅₀<1.
Thus, the left-hand side in (17) is positive for all 𝜆≥0 and 𝜏>0 while the right-hand side is zero. This leads to a contradiction. Then, (17) does not have non-negative real solutions. Following the proof of Theorem 1, we know that (17) has eigenvalues with negative real parts when 𝜏=0. Hence, for the case 𝜏>0, if (17) has roots with non-negative real parts they must be complex roots. Moreover, these complex roots should be obtained from a pair of complex conjugate roots crossing the imaginary axis. Thus, (17) must have a pair of purely imaginary roots.
Suppose 𝜆=𝜔𝑖 (𝜔>0) is a root of (17). Then, we have −𝜔³𝑖−𝑎₁𝜔²+𝑎₂𝜔𝑖+𝑎₃=𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}(cos(𝜔𝜏)−𝑖sin(𝜔𝜏)), (20) where 𝑎₃=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌. (21)
Therefore, −𝜔³+𝑎₂𝜔=−𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}sin(𝜔𝜏)𝑏₃−𝑏₁𝜔²,−𝑎₁𝜔²+𝑎₃=𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}cos(𝜔𝜏). (22) From (22), we can get −𝜔³+𝑎₂𝜔²+−𝑎₁𝜔²+𝑎₃²=𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏2}, (23)
that is, 𝜔⁶+𝑐₁𝜔⁴+𝑐₂𝜔²+𝑐₃=0. (24)
Letting 𝑧=𝜔², we obtain 𝑧³+𝑐₁𝑧²+𝑐₂𝑧+𝑐₃=0, (25)
where 𝑐₁=𝑎²₁−2𝑎₂𝑐>0,₂=𝑎²₂−2𝑎₁𝑎₃𝑐>0,₃=𝑎²₃−𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏2}=𝑎₃𝑎₃+𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}1−𝑅₀. (26)
If 𝑅₀<1, 𝑐₁>0, 𝑐₂>0, and 𝑐₃>0, which implies that (25) does not have positive solutions, then (24) does not have purely imaginary solutions. Hence, If 𝑅₀<1, the disease-free equilibrium 𝐸₀ of the system (13) is locally asymptotically stable.

Remark 4. If 𝑅₀>1, 𝑅₀3<1⟺𝜏>𝜇_𝑠+𝛼_𝑠ln𝑅₀≜𝜏₀. (27) It is easy to see if 𝑅₀<1, we can get 𝑅₀<1. But if 𝑅₀>1 to get 𝑅₀<1 the time delay must satisfy 𝜏>𝜏₀. This implies that to eliminate the disease the time delay must satisfy 𝜏>𝜏₀. In fact, the time delay is decreasing as the global warming. If the time delay is decreased to be smaller than 𝜏₀, 𝑅₀ may be larger than 1, and then the disease may be prevalent.

Now, we turn to the study of the stability of the endemic equilibrium of model (13).

Theorem 5. If 𝑅₀>1, the endemic equilibrium 𝐸^∗ of the system (13) is locally asymptotically stable.

Proof. The eigenvalues of 𝐸^∗ are −𝜇_𝑚 and the roots of the equation are 𝜆⁴+𝑏₁𝜆³+𝑏₂𝜆²+𝑏₃𝜆+𝑏₄=0, (28) where 𝑏₁=𝜇_𝑠+𝛼_𝑠+𝜇_𝑝+𝜇_𝑓+𝜌+𝜇_𝑠+𝜉𝑝^∗𝑏>0,₂=𝜇_𝑠+𝛼_𝑠𝜇_𝑝+𝜇_𝑠+𝛼_𝑠𝜇_𝑓+𝜌+𝜇_𝑠+𝛼_𝑠𝜇_𝑠+𝜉𝑝^∗+𝜇_𝑝𝜇_𝑓+𝜌+𝜇_𝑝𝜇_𝑠+𝜉𝑝^∗+𝜇_𝑓+𝜌𝜇_𝑠+𝜉𝑝^∗,𝑏₃=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌+𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑠+𝜉𝑝^∗+𝜇_𝑝𝜇_𝑓+𝜌𝜇_𝑠+𝜉𝑝^∗−𝜉𝜌𝑘_𝑓𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝑒^−𝜆,𝑏₄=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌𝜇_𝑠+𝜉𝑝^∗−𝜉𝜌𝑘_𝑓𝜇_𝑠𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝑒^−𝜆. (29)
Assume that 𝜆≥0. Since 𝜉𝜌𝑘_𝑓𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}=(𝜇_𝑠+𝛼_𝑠)𝜇_𝑝(𝜇_𝑓+𝜌) and 0<𝑒^−𝜆≤1 for 𝜏>0, we have 𝑏₁>0, 𝑏₂>0, 𝑏₃>0, and 𝑏₄>0 which lead to a contradiction in (28). Then, (28) does not have non-negative real solutions. From Theorem 2, we know that (28) has solutions with negative real parts for 𝜏=0. Hence, for the case 𝜏>0, if (28) has roots with non-negative real parts they must be complex roots. Moreover these complex roots should be obtained from a pair of complex conjugate roots crossing the imaginary axis. Thus, (28) must have a pair of purely imaginary roots.
Suppose 𝜆=𝜔𝑖 (𝜔>0) is a root of (28). Then, we have 𝜔⁴−𝑏₁𝜔³𝑖−𝑏₂𝜔²+𝑏̂₃𝜔𝑖+𝑏̂₄=𝜉𝜌𝑘_𝑓𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝜇_𝑠+𝜔𝑖(cos(𝜔𝜏)−𝑖sin(𝜔𝜏)), (30) where 𝑏̂₃=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌+𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑠+𝜉𝑝^∗+𝜇_𝑝𝜇_𝑓+𝜌𝜇_𝑠+𝜉𝑝^∗,𝑏̂₄=𝜇_𝑠+𝛼_𝑠𝜇_𝑝𝜇_𝑓+𝜌𝜇_𝑠+𝜉𝑝^∗. (31)
Therefore, 𝜔⁴−𝑏₂𝜔²+𝑏̂₄=𝜉𝜌𝑘_𝑓𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝜇_𝑠cos(𝜔𝜏)+𝜔sin(𝜔𝜏),−𝑏₁𝜔³+𝑏̂₃𝜔=𝜉𝜌𝑘_𝑓𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏}𝜔cos(𝜔𝜏)−𝜇_𝑠sin(𝜔𝜏). (32) From (32), we can get 𝜔⁴−𝑏₂𝜔²+𝑏̂₄²+−𝑏₁𝜔³+𝑏̂₃𝜔²=𝜉𝜌𝑘_𝑓𝑠^∗̃𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏2}𝜇²_𝑠+𝜔². (33)
Letting 𝑧=𝜔² again, we obtain 𝑧⁴+𝑐₁𝑧³+𝑐₂𝑧²+𝑐₃𝑧+𝑐₄=0, (34) where 𝑐₁=𝑏²₁−2𝑏₂,𝑐₂=𝑏²₂+2𝑏̂₄−2𝑏₁𝑏̂₃,𝑐₃=𝑏̂²₃−2𝑏₂𝑏̂₄−𝜉𝜌𝑘_𝑓𝑠^∗𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏2},𝑐₄=𝑏̂²₄−𝜉𝜌𝑘_𝑓𝑠^∗𝑒^{−(𝜇𝑠+𝛼𝑠)𝜏2}𝜇²_𝑠. (35)
By similar calculations, we can get 𝑐₁>0, 𝑐₂>0, 𝑐₃>0, and 𝑐₄>0. Thus, (34) does not have positive roots and then (28) cannot have purely imaginary solutions. Hence, If 𝑅₀>1 the endemic equilibrium 𝐸^∗ of the system (13) is locally asymptotically stable.

5. Simulation and Sensitivity Analysis

In this section, we perform some computational simulations to observe the impact of the time delay on schistosomiasis dynamics. Through sensitivity analysis of the death rate of snails, the death rate of single parasites, and the death rate of pairs, we give preferable control strategies.

The parameters are chosen with 𝑘_𝑓=100 per year, 𝑘_𝑚=145 per year [13, 18], Λ_𝑠=150 per year, 𝜉=0.000018 per year [18], 𝜇_𝑠=0.1 per year, 𝛼_𝑠=0.5 per year [29], 𝜇_𝑓=0.2 per year, 𝜇_𝑚=0.1 per year, 𝜇_𝑝=0.02 per year, 𝜌=0.467 [9] per year and 𝜏=18/365~38/365≈0.04~0.1 per year [4]. Then 𝑅₀=12.55, 𝑅₀=149.9 and 𝜏₀=8.43 when we choose 𝜏=30/365=0.082.

By Remark 4, the time delay 𝜏=0.082<8.43=𝜏₀ and 𝑅₀>1 implies that the disease will be prevalent. From Figure 1 we can find, the number of male schistosomes and infected snails increases as the time delay decreases. These phenomena imply that schistosomiasis infection becomes serious as the time delay decreases. As we know, the temperature of snails changes as the environmental temperature changes. Then the environmental temperature can influence the incubation period in infected snails. Currently, the environmental temperature is increasing because of the global warming. Hence, the global warming may reduce the incubation period of infected snails. Consequently, schistosomiasis may become more serious because of the global warming. Hence, the impact of global warming and the incubation period of infected snails on schistosomiasis dynamics can not be ignored.

(a) The trajectory of 𝑚

(b) The trajectory of 𝑖

Figure 1 

The dot line represents 𝜏=0.041, and the real line represents 𝜏=0.082.

As for schistosomiasis control, we can control the disease by killing schistosoma and (or) snails. In reality, there are many drugs for killing schistosoma. When we treat definitive hosts with chemotherapy, different drugs aim at different parasites. For example, male schistosoma is more sensitive to praziquantel (PZQ) [30, 31], but pair schistosoma is more sensitive to artesunate (ART) [32]. Hence, sometimes the death rate of pairs and single schistosoma is increased under some treatments, but the extent of increase is different. Under the treatment of some drugs more pairs are killed, but under the treatment of another drugs, more single schistosomes are killed. However, it is a pity that we lack the data about the rate of killing schistosoma to numerical simulation. Hence, we have done some qualitative analyses from the perspective of sensitivity.

Denote by 𝜀_𝑚, 𝜀_𝑓, 𝜀_𝑝 and 𝜀_𝑠 the killing rates of male schistosoma, female schistosoma, paired schistosoma and snails, respectively. We modify the model by adding these killing rates with the following form: 𝑑𝑚𝑑𝑡=𝑘_𝑚𝑖(𝑡)−𝜇_𝑚+𝜀_𝑚𝑚(𝑡)−𝜌𝑓(𝑡),𝑑𝑓𝑑𝑡=𝑘_𝑓𝑖(𝑡)−𝜇_𝑓+𝜀_𝑓𝑓(𝑡)−𝜌𝑓(𝑡),𝑑𝑝(𝑡)−𝜇𝑑𝑡=𝜌𝑓_𝑝+𝜀_𝑝𝑝(𝑡),𝑑𝑠𝑑𝑡=Λ_𝑠−𝜇_𝑠+𝜀_𝑠𝑠(𝑡)−𝜉𝑝(𝑡)𝑠(𝑡),𝑑𝑖(𝑡)𝑠(𝑡)−𝜇𝑑𝑡=𝜉𝑝_𝑠+𝜀_𝑠+𝛼_𝑠𝑖(𝑡). (36)

For convenience of computation, we denote 𝜇_𝑚+𝜀_𝑚=𝜇_𝑚𝜀,𝜇_𝑓+𝜀_𝑓=𝜇_𝑓𝜀,𝜇_𝑝+𝜀_𝑝=𝜇_𝑝𝜀,𝜇_𝑠+𝜀_𝑠=𝜇_𝑠𝜀. (37)

Similarly, we can obtain that the basic reproduction number for model (36) is as follows 𝑅_0𝜀=³𝜉𝜌Λ_𝑠𝑘_𝑓𝜇_𝑠𝜀𝜇_𝑠𝜀+𝛼_𝑠𝜇_𝑝𝜀𝜇_𝑓𝜀.+𝜌 (38) Note that 𝑅_0𝜀 is independent of the death rate of single male schistosoma. It is easy to see that 𝑅_0𝜀≤𝑅₀, which implies that killing schistosoma and (or) snails can reduce the basic reproduction number and consequently control the disease if 𝑅_0𝜀<1. However, we should choose a rather effective control measure in reality. Thus, we perform some sensitivity analyses by comparing the flexibilities of killing rates. The flexibilities of 𝜇_𝑠𝜀, 𝜇_𝑓𝜀, and 𝜇_𝑝𝜀 on the basic reproduction number 𝑅_0𝜀, respectively, are 𝐸𝑅_0𝜀𝐸𝜇_𝑠𝜀=𝜕𝑅_0𝜀𝜕𝜇_𝑠𝜀×𝜇_𝑠𝜀𝑅_0𝜀1=−32𝜇_𝑠𝜀+𝛼_𝑠𝜇_𝑠𝜀+𝛼_𝑠,𝐸𝑅_0𝜀𝐸𝜇_𝑓𝜀=𝜕𝑅_0𝜀𝜕𝜇_𝑓𝜀×𝜇_𝑓𝜀𝑅_0𝜀1=−3𝜇_𝑓𝜀𝜇_𝑓𝜀,+𝜌𝐸𝑅_0𝜀𝐸𝜇_𝑝𝜀=𝜕𝑅_0𝜀𝜕𝜇_𝑝𝜀×𝜇_𝑝𝜀𝑅_0𝜀1=−3. (39)

It is easy to see that 𝐸𝑅_0𝜀𝐸𝜇_𝑓𝜀<13=𝐸𝑅_0𝜀𝐸𝜇_𝑝𝜀<𝐸𝑅_0𝜀𝐸𝜇_𝑠𝜀. (40) It implies that the death rate of snails is the most sensitive parameter to the reduction of 𝑅_0𝜀. In addition, the death rate of pair schistosoma is more sensitive than that of female single schistosoma. This is reasonable because the offspring is reproduced by pair schistosoma and the pair schistosoma has stronger effect on the life cycle of schistosome. Hence, to kill the sail hosts may be the first reasonable control measure. Furthermore, to kill pair schistosoma by using some medications may be the other better control measure.

6. Discussion

In this paper, we established a new schistosomiasis model including a more reasonable sex ratio of schistosoma, snail dynamics, the latent period of infected snails and mating structure. By choosing the minimum function as the mating function, we studied the stabilities of model (6) without time delay and model (13) with time delay. If the basic reproduction number of model (6) is less than 1, one can prove the stability of the disease-free equilibrium in system (6). When the basic reproductive number 𝑅₀ is greater than 1, the stability of the endemic equilibrium in system (6) can be obtained. But to get the stability of the disease-free equilibrium in system (13) the time delay must be larger than 𝜏₀. In reality, the environmental temperature because of the global warming is increasing. And then the incubation period of infected snails is shortened as the environmental temperature increases [13]. From Figure 1 we know that schistosomiasis infection become more serious as the incubation period of infected snails decreases. Hence, the impact of the incubation period of infected snails during global warming on schistosomiasis transmission cannot be ignored.

In recent years, the control of schistosomiasis remains one of the highest priorities in parasitology. There are many strategies to control schistosomiasis, such as schistosome vaccine, killing snails and chemotherapy with PZQ. In this paper, through comparing the flexibilities of 𝜇_𝑠𝜀, 𝜇_𝑓𝜀, and 𝜇_𝑝𝜀 on 𝑅_0𝜀, we know that the death rate of snails 𝜇_𝑠𝜀 is the most sensitive to the reduction of 𝑅_0𝜀. Hence killing the snail hosts may be the most reasonable control measure.

Our sensitivity analysis also deduces that 𝑅_0𝜀 is more sensitive to the death rate of pair parasites than to the death rate of female parasite. This implies that killing paired schistosoma is more advantageous than killing single schistosoma to control disease. From the perspective of medicine, eggs are generally accepted as the major cause of pathogenesis of schistosomiasis [14, 30–32]. Nino Incani et al. [14] suggest that the eggs produced by paired schistosoma are the pivotal factor in nosogenesis and transmission, and sexual maturation is the precondition of producing eggs. This means that paired schistosoma is mostly harmful to definitive hosts. Controlling schistosomiasis must attack the ability of the paired schistosoma to produce eggs. Hence, the drug of choice is possibly sufficient for reducing egg-associated pathology.

As we know, PZQ is currently considered the drug of choice for the treatment of schistosomiasis because of its efficacy against all schistosome species, lack of serious side effects, and low cost. However, PZQ has its limitation. PZQ was found to be more active against single male schistosoma [30, 31], not paired schistosomes. In paired schistosoma, male schistosoma are folded together with female schistosoma. The egg production is owed to the paired schistosomes. Reducing male schistosoma in paired schistosoma can help reduce the egg production and consequently the overall parasite density and disease prevalence. In additional, from our analysis, 𝑅_0𝜀 is sensitive to the death rate of pair schistosoma and independent of the death rate of single male schistosoma. Furthermore, recent epidemiological evidence suggests the emergence of PZQ-resistant schistosoma [33]. Hence, PZQ may not be the best drug for schistosomiasis.

Currently, there are some drugs which are possibly sufficient for reducing egg-associated pathology. For example, ART significantly decreased the survival time of both paired male and female schistosomes. It might be responsible for damaging reproductive organs or killing schistosoma [32]. Oxamniquine-praziquantel composite in the treatment of schistosoma mansoni infection is also a controlled trial [31].

In summary, although PZQ remains the drug of choice to treat schistosomiasis, it does not protect from reinfection (especially in children) and is minimally effective against larval stages of the parasite. Resistance can also develop, although its mode of action is poorly understood. Therefore an important priority in developing new control strategies is to search new drug targets, in combination with selection of viable vaccine candidates. Based on our analysis in this paper, we can obtain the following two results. Killing the snail hosts may be the most reasonable control measure. If we choose chemotherapy, we should choose some drugs which are sufficient for reducing egg-associated pathology because paired schistosoma is mostly harmful to definitive hosts.

Acknowledgments

This research is supported by the National Natural Science Foundation of China (11126177 and 11071011), the Natural Science Foundation of Anhui Province (1208085QA15), the Foundation for Young Talents in College of Anhui Province (2012SQRL021), the National Scholarship Foundation of China, Funding Project for Academic Human Resources Development in Institutions of Higher Learning under the Jurisdiction of Beijing Municipality (PHR201107123), the Key Natural Science Foundation of the Anhui Higher Education Institutions of China (KJ2009A49), and the Doctoral Fund of Ministry of Education of China (20113401110001 and 20103401120002).