1. Introduction

Human schistosomiasis, the third most devastating tropical disease in the world after malaria and intestinal helminthiasis, is a global public health problem [42]. According to World Health Organization (WHO), the number of people need preventive chemotherapy globally in 2013 was 262 million, of which 121.2 million were school-aged children [42,43,44].

The major forms of human schistosomiasis are caused by species of the water-borne flatworm or blood flukes called schistosomes [9,44]. Schistosomiasis in mainland China is caused by Schistosoma Japonicum (S. Japonicum). Though its transmission had been interrupted successively in five of the twelve formerly endemic provinces (see [8,9] and the references therein), schistosomiasis is still one of the most serious parasitic diseases in China and remains endemic in seven provinces, including Hubei, Anhui, Hunan, Jiangsu and Jiangxi in the lake and marshland regions with vast areas of Oncomelania hupensis habita, Sichuan and Yunnan in the mountainous regions with diverse ecologies [8,9]. The monthly data of schistosomiasis cases in Hubei, Hunan and Anhui recorded by Chinese Center for Disease Control and Prevention (China CDC) [10] display a seasonal pattern. The cases in the late summer and early autumn are significantly higher than in the spring and winter (see Figure 1.).

Schistosomiasis is a parasitic disease caused by trematode flatworms of the genus schistosoma [6]. The reproductive cycle of schistosomiasis starts with parasitic eggs released into freshwater through faeces and urine, then some eggs hatch and became miracidia under appropriate conditions, those miracidia swim and penetrate snails as intermediate host. By escaping from the snail, the infective cercariae penetrate the skin of the human host. For more details on the life cycle of schistosome, we refer to [5,7,12,16,17,23,28,30,6]. To focus on the dynamics of S.Japonicum propagating between human and the intermediate host snails, we consider a simplified diagram for the life cycle given in Figure 2.

The earliest mathematical models for schistosomes were developed by Macdonald [32] and Hairston [25]. Since then, a good number of mathematical models involving the transmission dynamics of schistosomes have been proposed (see [8,12,16,17,30] and the references therein). Garira et al. [19] proposed a dynamic model of ordinary differential equations linking the within-host and between-host dynamics of infections with free-living pathogens in the water environment. Wang and Spear [48] explored the impact of infection-induced immunity on the transmission of S. Japonicum in hilly and mountainous environments in China, and underscored the need for improved diagnostic methods for disease control, especially in potentially re-emergent settings. Chen et al. [8] proposed an autonomous mathematical model for controlling schistosomiasis in Hubei Province, China, focusing on the disease spread among people, intermediate hosts snails and cattle. Feng et al. [16] estimated the parameters of a schistosome transmission system, which described the distribution of schistosome parasites in a village in Brazil.

Schistosomiasis often occurs in most tropical and some subtropical regions of the world. Environmental and climatic factors play an important role in the geographical distribution and transmission of schistosomiasis [44]. It was well known that seasonality can cause population fluctuations ranging from annual cycles to multi year oscillations, and even chaotic dynamics [2,22]. From an applied viewpoint, clarifying the mechanisms that link seasonal environmental changes to diseases dynamics may provide help in predicting the long-term health risks, in developing an effective public health program, and in setting objectives and utilizing limited resources more effectively (see [1,31,35] and the references therein). These considerations indicate that seasonal models are needed in order to describe the periodic incidence of schistosomiasis transmission. However, to the best of our knowledge, there are few studies modeling the seasonality influence on the transmission of schistosomiasis in mainland China [54].

More than 82% of infected persons lived in lake and marshland regions (such as Dongting Lake and Poyang Lake) along the Yangtze River, where interruption of transmission has been proven difficult [20,58,59]. The purpose of this paper is to propose a periodic schistosomiasis model to investigate the seasonal transmission dynamics and search for control strategies in these lake and marshland regions in China. We analyze the dynamical behavior, evaluate the basic reproduction number $R_{0},$ use the model to simulate the human cases in Hubei, and forecast the monthly tendency of schistosomiasis after January 2015. Moreover, we perform sensitivity analysis of the basic reproduction number $R_{0}$ in terms of key model parameters. Finally, like Hubei Province, Anhui, Hunan, Jiangsu, and Jiangxi are located in the lake and marshland regions in the central and eastern China, similar control and prevention measures can also be designed and proposed for these provinces.

The paper is organized as follows. In Section 2, we introduce the periodic schistosomiasis model. Some preliminary results are presented in Section 3, such as the positivity and boundedness of solutions and calculation of the basic reproduction number. The extinction and uniform persistence of the disease are discussed in Sections 4 and 5, respectively. Simulations of the schistosomiasis data from Hubei Province are presented in Section 6. Conclussion and discussion are given in Section 7.

---

2. Mathematical modeling

To study the seasonal transmission dynamics of schistosomiasis, we trace the life cycle of schistosome parasites in three different environments: human biological environment, physical water environment, and snail biological environment. The life cycle of schistosomiasis was given in Figure 2 and its transmission diagram among humans, snails, and miracidia and cercariae is illustrated in Figure 3.

We denote the total numbers of humans and snails by $N_{H}(t)$ and $N_{V}(t)$, respectively, and classify each of them into two subclasses: susceptible and infectious, with the numbers of humans denoted by $S_{H}(t)$ and $I_{H}(t)$, and snails sizes denoted by $S_{V}(t)$ and $I_{V}(t)$, respectively. The miracidia and cercariae dynamics are incorporated and their densities are denoted by $M(t)$ and $P(t)$, respectively. The mathematical model is derived based on the following basic assumptions:

(1) There is no vertical transmission of the disease.

(2) Susceptible humans are recruited at a positive constant rate $\Lambda_{H}$.

(3) There are no immigrations of infectious humans.

(4) People living near rivers and lakes are more likely going swimming and fishing in the summer and autumn, they are prone to infection for long contacting with contaminated water. The river, lake, pond water freezes or dry in winter, infected snail seldom or not produce larvae, then infection is not likely to happen. Due to these seasonal phenomena, we use two 12-month periodic functions $\lambda_{H}(t)=a_{H}[1 + b_{H}\sin(\frac{\pi}{6}t+\varphi_H)]$ and $\lambda_{V}(t)=a_{V}[1 + b_{V}\sin(\frac{\pi}{6}t+\varphi_V)]$ (see [54]) to describe the transmission rates from cercariae to human and from miracidia to snails, respectively, where positive constants $a_{H}$, $b_{H}$ and $\varphi_H$ represent the human baseline transmission rate, its magnitude of forcing and the initial phase, respectively, positive constants $a_{V}$, $b_{V}$ and $\varphi_V$ in $\lambda_{V}(t)$ have the similar meanings as constants in $\lambda_{H}(t)$. We choose bilinear incidence rates (density-dependent or mass action type) $\lambda_{H}(t)S_{H}P$ and $\lambda_{V}(t)S_{V}M$ (see [17]).

(5) It is clear that the snail population is seasonally changed in reality, the recruited rate $\Lambda_{V}(t),$ natural death rate $\mu_{V}(t)$ and disease induced death rate $\alpha_{V}(t)$ for the snail population are considered as 12-month periodic continuous functions. For more simulation details, see Section 6.

(6) There is no immune response in both snail and human populations.

(7) Several effective control strategies, such as drug treatment, improving sanitation and health education, the integrated strategies are considered here. We denote these strategies by the natural recovery and treatment rate $\gamma_{H}.$

(8) We further assume that the human natural death rate $\mu_{H}$, miracidia natural death rate $\mu_{M}$ and cercariae natural death rate $\mu_{P}$ are all positive constants. Miracidia migration rate $\lambda_{M}$ from human to snail and cercaria migration rate $\lambda_{P}$ from snail to human are also supposed to be positive constants.

The model is described by the following system of ordinary differential equations:

$\left\{ \begin{matrix} \begin{align} &{{S}_{{{H}'}}}(t)={{\Lambda }_{H}}-{{\lambda }_{H}}(t){{S}_{H}}P-{{\mu }_{H}}{{S}_{H}}+{{\gamma }_{H}}{{I}_{H}}, \\ &{{I}_{{{H}'}}}(t)={{\lambda }_{H}}(t){{S}_{H}}P-({{\mu }_{H}}+{{\gamma }_{H}}){{I}_{H}}, \\ &{M}'(t)={{\lambda }_{M}}{{I}_{H}}-{{\mu }_{M}}M, \\ &{{S}_{{{V}'}}}(t)={{\Lambda }_{V}}(t)-{{\lambda }_{V}}(t){{S}_{V}}M-{{\mu }_{V}}(t){{S}_{V}}, \\ &{{I}_{{{V}'}}}(t)={{\lambda }_{V}}(t){{S}_{V}}M-{{\alpha }_{V}}(t){{I}_{V}}, \\ &{P}'(t)={{\lambda }_{P}}{{I}_{V}}-{{\mu }_{P}}P. \\ \end{align}&{}&{}&{}&{}&{}&{} \\ \end{matrix} \right.$

(1)

---

3. Basic properties

We denote $\omega=12$ months. Based on the biological background of model (1), we only consider solutions of model (1) starting at $t=0$ with initial values:

$\begin{array}{rl} S^{0}_{H}>0, I^{0}_{H}\geq 0, M^{0}\geq 0, S^{0}_{V}>0, I^{0}_{V}\geq 0, P^{0}\geq 0. \end{array}$

(2)

When $I_{H}=0, M=0, I_{V}=0$ and $P=0,$ model (1) has a unique disease-free periodic solution $E^{0}=(\widehat{S_{H}}, 0, 0, \widehat{S_{V}}(t), 0, 0),$ where $\widehat{S_{H}}=\frac{\Lambda_{H}}{\mu_{H}},$ and $\widehat{S_{V}}(t)$ is the globally asymptotically stable positive $\omega$-periodic solution of equation $S_{V}'(t)=\Lambda_{V}(t)-\mu_{V}(t)S_{V}(t)$ (see Lemma 1 in [40]).

Now, we deduce the basic reproduction number $R_{0}$ for model (1) following the general calculation procedure in Wang and Zhao [50]. Firstly, we can validate that model (1) satisfies the conditions $(A_1)-(A_7)$ given in [50].

Denote

$F(t)=\left( \begin{matrix} 0&0&0&{{\lambda }_{H}}(t)\widehat{{{S}_{H}}} \\ 0&0&0&0 \\ 0&{{\lambda }_{V}}(t)\widehat{{{S}_{V}}}(t)&0&0 \\ 0&0&0&0 \\ \end{matrix} \right),V(t)=\left( \begin{matrix} {{\mu }_{H}}+{{\gamma }_{H}}&0&0&0 \\ -{{\lambda }_{M}}&{{\mu }_{M}}&0&0 \\ 0&0&{{\alpha }_{V}}(t)&0 \\ 0&0&-{{\lambda }_{P}}&{{\mu }_{p}} \\ \end{matrix} \right).$

Let $Y(t, s)$ be the $4\times4$ matrix solution of the following initial value problem

$\left\{ \begin{array}{lll} \displaystyle\frac{dY(t, s)}{dt}&=-V(t)Y(t, s), \forall t\geq s, \\ \displaystyle Y(s, s)&=I, \end{array} \right.$

where $I$ is the $4\times4$ identity matrix. Further, let $C_{\omega}$ be the ordered Banach space of all $\omega-$periodic continuous functions from $\mathbb{R}$ to $\mathbb{R}^4$ with maximum norm $\|\cdot\|$ and positive cone $C^{+}_{\omega}:=\{\phi\in C_{\omega}: \phi(t)\geq 0, \forall t\in R\}.$ Suppose $\phi(s)\in C^{+}_{\omega}$ is the initial distribution of infectious individuals, then $F(s)\phi(s)$ is the rate of new infection produced by the infectious individuals who were introduced at time $s,$ and $Y(t, s)F(s)\phi(s)$ represents the distributions of those infectious individuals who were newly infected at time $s$ and remain in the infected compartment at time $t$ for $t\geq s.$ Naturally,

$\displaystyle \int_{-\infty}^{t}Y(t, s)F(s)\phi(s)ds=\int_{0}^{+\infty}Y(t, t-a)F(t-a)\phi(t-a)da$

is the distribution of accumulative new infections at time $t$ produced by all those infected individuals $\phi(s)$ introduced at time previous to $t$. Then, we define a linear operator $L: C_{\omega} \rightarrow C_{\omega}$ as follows

$\displaystyle (L\phi)(t)=\int_{0}^{+\infty}Y(t, t-a)F(t-a)\phi(t-a)da, \forall t\in R, \phi\in C_{\omega}.$

$L$ is called the next infection operator.

Applying the results obtained in [50], the basic reproduction number $R_{0}$ for model (1) is defined as the spectral radius of operator $L$ ([3,4]); that is, $R_{0}=\rho(L).$

Employing Theorem 2.1 and Theorem 2.2 in Wang and Zhao [50], we can deduce the following results with respect to $R_{0}$ and the locally asymptotical stability of the disease-free periodic solution $E^{0}$ for model (1).

Lemma 3.1. On basic reproduction number $R_{0}$, we have

(ⅰ) $R_{0} < 1$ if and only if $\rho(\Phi_{F-V}(\omega)) < 1$;

(ⅱ) $R_{0}=1$ if and only if $\rho(\Phi_{F-V}(\omega))=1$;

(ⅲ) $R_{0}>1$ if and only if $\rho(\Phi_{F-V}(\omega))>1$.

Then $E^{0}$ is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0}>1$, where $\Phi_{F-V}(t)$ is the monodromy matrix of the linear $\omega-$periodic system $\frac{dz}{zt}=[F(t)-V(t)]z.$

On the positivity and boundedness of solutions of model (1) with nonnegative initial conditions (2), we have the following results.

Lemma 3.2. Let $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))$ be the solution of model (1) with initial conditions (2). Then $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))$ is nonnegative for all $t\geq 0$ and ultimately bounded. In particular, if $S^{0}_{H}>0, I^{0}_{H}>0, M^{0}>0, S^{0}_{V}>0, I^{0}_{V}>0$ and $P^{0}>0$, then $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))$ is also positive for all $t>0.$

Proof. In fact, by the continuous dependence of solutions with respect to initial values, we only need to prove that when $S^{0}_{H}>0, I^{0}_{H}>0, M^{0}>0, S^{0}_{V}>0, I^{0}_{V}>0$ and $P^{0}>0$, $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))$ is positive for all $t>0.$ Set

$m(t)=\min\{S_H(t), I_H(t), M(t), S_V(t), I_V(t), P(t) \}, \forall t>0.$

Clearly, $m(0)>0.$ Assuming that there exists a $t_{1}>0$ such that $m(t_1)=0$ and $m(t)>0$ for all $t\in[0, t_{1}).$

If $m(t_1)=S_{H}(t_{1}),$ since $I_{H}(t)>0$ for all $t\in[0, t_{1}),$ from the first equation of model (1), it follows that $S'_{H}(t)\geq-(\lambda_{H}(t)P+\mu_{H})S_{H}$ for all $t\in[0, t_{1}]$. Then

$0=S_{H}(t_{1})\geq S^{0}_{H}\exp\Big(-\int_{0}^{t_{1}}(\lambda_{H}(s)P+\mu_{H})ds\Big)>0,$

which leads to a contradiction.

Similar contradictions can be deduced in the cases of $m(t_1)=I_{H}(t_{1}),$ $m(t_1)=M(t_{1}),$ $m(t_1)=S_{V}(t_{1}),$ $m(t_1)=I_{V}(t_{1})$ and $m(t_1)=P(t_{1}).$ Therefore, $S_{H}(t)>0,$ $I_{H}(t)>0,$ $M(t)>0,$ $S_{V}(t)>0,$ $I_{V}(t)>0,$ and $P(t)>0$ for all $t>0$.

Let $N_{H}(t)=S_{H}(t)+I_{H}(t).$ We have $N'_{H}(t)= \Lambda_{H}-\mu_{H}N_{H}(t)$ which implies that $N_{H}(t)= \frac{\Lambda_{H}}{\mu_{H}}+N^{0}_{H}\exp(-\mu_{H}t),$ where $N^{0}_{H}=S^{0}_{H}+I^{0}_{H}.$ Hence, $N_{H}(t)$ is bounded for all $t\geq 0$ and

$\begin{matrix} \underset{{}}{\mathop{\lim {{\sup }_{t\to \infty }}}}\,{{N}_{H}}(t)=\frac{{{\Lambda }_{H}}}{{{\mu }_{H}}}:={{B}_{H}}, \\ \end{matrix}$

(3)

which implies that $S_{H}(t)$ and $I_{H}(t)$ are also bounded for $t>0$. From the third equation of model (1), we know that for any $\varepsilon>0$ there is a $T_{0}>0$ such that

$M'(t)=\lambda_{M}I_{H}-\mu_{M}M\leq\lambda_{M}(B_{H}+\varepsilon)-\mu_{M}M, \forall t\geq T_{0}.$

Then, we have $M(t)\leq\frac{\lambda_{M}(B_{H}+\varepsilon)}{\mu_{M}}+M^{0}\exp(-\mu_{M}(t-T_{0})).$ Hence, $M(t)$ is bounded for $t>0$ and by the arbitrariness of $\varepsilon$ we also have

$\begin{matrix} \lim {{\sup }_{t\to \infty }}M(t)\le \frac{{{\lambda }_{M}}{{B}_{H}}}{{{\mu }_{M}}}. \\ \end{matrix}$

(4)

Set $N_{V}(t)=S_{V}(t)+I_{V}(t).$ From the forth and fifth equations of model (1) we have $N'_{V}(t)\leq \Lambda_{V}(t)-\mu_{V}(t)N_{V}(t)$. By the comparison principle and Lemma 1 in [40], we can obtain

$\lim {{\sup }_{t\to \infty }}({{N}_{V}}(t)-\widehat{{{S}_{V}}}(t))\le 0,$

(5)

which implies that $S_{V}(t)$ and $I_{V}(t)$ are bounded for $t>0$.

Lastly, from the last equation of model (1), similar to the proof of (4) we obtain $\limsup_{t\rightarrow\infty}P(t)\leq\frac{\lambda_{P}B_{V}}{\mu_{P}}$, where $B_{V}=\sup_{t\in[0, \omega]}\widehat{S_{V}}(t)$. This completes the proof.

Remark 3.3. Denote set $\Omega$ as follows

$\begin{array}{rl} \Omega=&\displaystyle\{(S_{H}, I_{H}, M, S_{V}, I_{V}, P): 0\leq N_{H}\leq B_{H}, 0\leq M\leq \frac{\lambda_{M}B_{H}}{\mu_{M}}, \\ &\displaystyle 0\leq N_{V}\leq B_{V}, 0\leq P\leq\frac{\lambda_{P}B_{V}}{\mu_{P}}\}. \end{array}$

Lemma 3.2 implies that $\Omega$ is a positively invariant set with respect to model (1).

---

4. Extinction of disease

Theorem 4.1. If $R_{0}< 1,$ then the disease-free periodic solution $E^{0}$ of model (1) is globally asymptotically stable.

Proof. By considering the linearization system, we can prove that $E^{0}$ is locally asymptotically stable when $R_{0} < 1$, which is equivalent to $\rho(\Phi_{F-V}(\omega)) < 1$ by Lemma 3.1. We can choose a small enough positive constant $\varepsilon$ such that $\rho(\Phi_{F-V+\varepsilon N}(\omega)) < 1$, where

$N(t)=\left( \begin{array}{cccc} 0&0&0 &\lambda_{H}(t) \\ 0 &0&0& 0 \\ 0&\lambda_{V}(t)&0&0 \\ 0&0&0&0 \\ \end{array} \right).$

(6)

Let $(S_H(t), I_H(t), M(t), S_V(t), I_V(t), P(t))$ be a positive solution of model (1). From (3) and (5), for any given $\varepsilon$, there exists a $t_{1}$ such that $S_{H}(t)\leq \frac{\Lambda_{H}}{\mu_{H}}+\varepsilon$ and $S_{V}(t)\leq \widehat{S_{V}}(t)+\varepsilon$ for all $t>t_{1}$. Then for all $t>t_{1}$, from model (1) we obtain that

$\left\{ \begin{array}{lll} \displaystyle I_{H}'(t)\leq\lambda_{H}(t)(\frac{\Lambda_{H}}{\mu_{H}}+\varepsilon)P-(\mu_{H}+\gamma_{H})I_{H}, \\ \displaystyle M'(t)=\lambda_{M}I_{H}-\mu_{M}(t)M, \\ \displaystyle I_{V}'(t)\leq\lambda_{V}(t)(\widehat{S_{V}}(t)+\varepsilon)M-\alpha_{V}(t)I_{V}, \\ \displaystyle P'(t)=\lambda_{P}I_{V}-\mu_{P}(t)P. \end{array} \right.$

(7)

Considering the following auxiliary system:

$\left\{ \begin{array}{lll} \displaystyle \widetilde{I}_{H}'(t)=\lambda_{H}(t)(\frac{\Lambda_{H}}{\mu_{H}}+\varepsilon)\widetilde{P}-(\mu_{H}+\gamma_{H})\widetilde{I}_{H}, \\ \displaystyle \widetilde M'(t)=\lambda_{M}\widetilde{I}_{H}-\mu_{M}(t)\widetilde{M}, \\ \displaystyle \widetilde{I}_{V}'(t)=\lambda_{V}(t)(\widehat{S_{V}}(t)+\varepsilon)\widetilde{M}-\alpha_{V}(t)\widetilde{I}_{V}, \\ \displaystyle \widetilde{P}'(t)=\lambda_{P}\widetilde{I}_{V}-\mu_{P}(t)\widetilde{P}, \end{array} \right.$

that is

$\begin{array}{rl} \displaystyle \frac{dh(t)}{dt}=(F(t)-V(t)+\varepsilon N)h(t), h(t)=(\widetilde{I}_{H}(t), \widetilde{M}(t), \widetilde{I}_{V}(t), \widetilde{P}(t))^{T}. \end{array}$

(8)

By Lemma 2.1 in Zhang and Zhao [56], it follows that there exists a positive $\omega-$periodic function $\phi(t)=(\phi_{1}(t), \phi_{2}(t), \phi_{3}(t), \phi_{4}(t))^{T}$ such that $h(t)=e^{\mu t}\phi(t)$ is a solution of system (8), where $\mu=\frac{1}{\omega}\ln\rho(\Phi_{F-V+\varepsilon N}(\omega)).$

Let $J(t)=({I}_{H}(t), M(t), {I}_{V}(t), P(t))$. We can choose a small enough positive constant $\eta>0$ such that $J(t_{1})\leq \eta \phi(t_{1})$. Then, from (7) and the comparison principle, we have $J(t)\leq \eta e^{\mu t} \phi(t)$ for all $t>t_{1}$. By $\rho(\Phi_{F-V+\varepsilon N}(\omega)) < 1,$ it follows that $\mu < 0,$ then $\lim_{t\rightarrow\infty}J(t)=0,$ which in turn implies that $\lim_{t\rightarrow \infty}(I_{H}(t), M(t), I_{V}(t), P(t))=(0, 0, 0, 0)$. Moreover, by the first and fourth equations in model (1), we obtain that $\lim_{t\rightarrow \infty}S_{H}(t)=\frac{\Lambda_{H}}{\mu_{H}}$ and $\lim_{t\rightarrow \infty}S_{V}(t)=\widehat{S_{V}}(t)$. Therefore, $E^{0}$ is globally attractive when $R_{0} < 1$. This completes the proof.

When $\lambda_{H}(t)\equiv\lambda_{H},$ $\lambda_{V}(t)\equiv\lambda_{V},$ $\Lambda_{V}(t)\equiv \Lambda_{V},$ $\mu_V(t)\equiv\mu_V$ and $\alpha_V(t)\equiv\alpha_V$ are positive constants, model (1) reduces to an autonomous case. Using the method given by van den Driessche and Watmough [41], we obtain the corresponding basic reproduction number

$\begin{array}{rl} \displaystyle \tilde{R}_{0}=\rho(FV^{-1})= \frac{\lambda_{H}\Lambda_{H}\lambda_{M}\lambda_{P}}{\mu_{H}\mu_{M}\mu_{P}(\mu_{H}+\gamma_{H})} \end{array}.$

(9)

As a corollary of Theorem 4.1, we have the following result.

Corollary 4.2. If $\tilde{R}_{0} < 1$, then for the autonomous model (1), the disease-free equilibrium $(\frac{\Lambda_{H}}{\mu_{H}}, 0, 0, \frac{\Lambda_{V}}{\mu_{V}}, 0, 0)$ is globally asymptotically stable.

---

5. Uniform persistence of disease

Theorem 5.1. If $R_{0}>1,$ then model (1) is uniformly persistent; that is, there exists a positive constant $\epsilon$ such that any solution $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))$ of model (1) with initial conditions $S^{0}_{H}>0, I^{0}_{H}> 0, M^{0}> 0, S^{0}_{V}>0, I^{0}_{V}> 0,$ and $P^{0}>0$ satisfies

$\begin{array}{lll} \displaystyle \liminf\limits_{t\rightarrow \infty}(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))\geq(\epsilon, \epsilon, \epsilon, \epsilon, \epsilon, \epsilon). \end{array}$

Proof. By $R_{0}>1,$ which is equivalent to $\rho(\Phi_{F-V}(\omega))>1$ by Lemma 3.1, we can choose a small constant $\theta>0$ such that $\rho(\Phi_{F-V-\theta N}(\omega))>1$, $\frac{\Lambda_{H}}{\mu_{H}}-\theta>0$ and $\widehat{S_{V}}(t)-\theta>0$, where $N(t)$ is given in (6).

For any small enough constant $\varepsilon>0$, we consider the following two perturbed equations

$\begin{array}{rl} \displaystyle U'_{\varepsilon}(t)=\Lambda_{H}-\varepsilon\lambda_{H}(t)U_{\varepsilon}(t)-\mu_{H}U_{\varepsilon}(t) \end{array}$

(10)

and

$\begin{array}{rl} \displaystyle V'_{\varepsilon}(t)=\Lambda_{V}(t)-\varepsilon\lambda_{V}(t)V_{\varepsilon}(t)-\mu_{V}(t)V_{\varepsilon}(t). \end{array}$

(11)

Applying Lemma 2 in [39] and Lemma 1 in [40], equations (10) and (11) admit globally uniformly attractive positive $\omega-$periodic solutions $U^{*}_{\varepsilon}(t)$ and $V^{*}_{\varepsilon}(t)$, respectively. By the continuity of solutions with respect to parameter $\varepsilon,$ for constant $\theta>0$ given above, there exists a constant $\varepsilon_{1}>0$ such that for all $0 < \varepsilon_{1} < \varepsilon$ and $t\in[0, \omega],$

$\begin{array}{lll} \displaystyle U^{*}_{\varepsilon_{1}}(t)>\frac{\Lambda_{H}}{\mu_{H}}-\frac{\theta}{2}, V^{*}_{\varepsilon_{1}}(t)>\widehat{S_{V}}(t)-\frac{\theta}{2}. \end{array}$

(12)

Since model (1) is $\omega-$periodic, we can use the persistence theory of dynamical systems given in [57] to discuss the permanence of model (1). Let

$X=\{(S_{H}, I_{H}, M, S_{V}, I_{V}, P):S_{H}>0, I_{H}\geq0, M\geq0, S_{V}>0, I_{V}\geq0, P\geq0\}$

and

$X_{0}=\{(S_{H}, I_{H}, M, S_{V}, I_{V}, P)\in X:I_{H}>0, M>0, I_{V}>0, P>0\}.$

Then

$\partial X_{0}=X\backslash X_{0}=\{(S_{H}, I_{H}, M, S_{V}, I_{V}, P)\in X: I_{H}MI_{V}P=0\}.$

By Lemma 3.2, $X$ and $X_{0}$ are positively invariant with respect to model (1), and $\partial X_{0}$ is a relatively closed set in $X$.

Define $P: X\rightarrow X$ as the Poincaré map associated with model (1); that is

$P(x^{0})=u(\omega, x^{0}), \forall x^{0}\in X,$

where $u(t, x^{0})$ is the unique solution of model (1) with initial values $u(0^{+}, x^{0})=x^{0}$ and $x^{0}=(S^{0}_{H}, I^{0}_{H}, M^{0}, S^{0}_{V}, I^{0}_{V}, P^{0}).$ By Remark 3.3, the Poincaré map $P$ is compact and point dissipative on $X$. Therefore, Theorem 1.1.3 and ($C_{1}$) of Theorem 1.3.1 in [57] hold.

Let $M_{\partial}=\{x^{0}\in\partial X_{0}: P^{n}(x^{0})\in\partial X_{0}, n=1, 2, \cdots\},$ where $P^{n}=P(P^{n-1})$, $n>1$, and $P^{1}=P$. We firstly verify that

$\begin{array}{lll} M_{\partial}=\{(S_{H}, 0, 0, S_{V}, 0, 0): S_{H}>0, S_{V}>0\}. \end{array}$

(13)

If initial conditions $(S^{0}_{H}, I^{0}_{H}, M^{0}, S^{0}_{V}, I^{0}_{V}, P^{0})=(S^{0}_{H}, 0, 0, S^{0}_{V}, 0, 0)$ with $S^{0}_{H}>0$ and $S^{0}_{V}>0,$ then the solution $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t) )\equiv(S_{H}(t), 0, 0, \\ S_{V}(t), 0, 0)$ with $S_{H}(t)>0$ and $S_{V}(t)>0,$ it is clear that $\{(S_{H}, 0, 0, S_{V}, 0, 0): S_{H}>0, S_{V}>0\}\subseteq M_{\partial}.$ On the other hand, if $M_{\partial}\backslash\{(S_{H}, 0, 0, S_{V}, 0, 0): S_{H}>0, S_{V}>0\}\neq \emptyset$, then there exists at least a point $(S^{0}_{H}, I^{0}_{H}, M^{0}, S^{0}_{V}, I^{0}_{V}, P^{0}) \in M_{\partial}$ satisfying $I^{0}_{H}>0$ or $I^{0}_{V}>0$ or $M^{0}>0$ or $P^{0}>0.$

If $I^{0}_{H}>0$, from the second equation of model (1), we have for all $t>0$ that

$I_{H}(t)\geq I_{H}^{0}e^{-(\mu_{H} +\alpha_{H}+\gamma_{H})t}>0.$

Thus, by the third equation of model (1), $M(t)> M^{0}e^{-\mu_{M}t}\geq0$ for all $t>0$. From $S_{H}^{0}>0,$ we can obtain from the forth equation of model (1) that $S_{H}(t)>0$ for all $t>0.$ Therefore, by the fifth equation of model (1), we can easily get $I_{V}(t)>I_{V}^{0}e^{-\int_{0}^{t} \alpha_{V}(\tau)d\tau}\geq0$ for all $t>0$. Furthermore, from the sixth equation of model (1) we can obtain $P(t)>P^{0}e^{-\lambda_{P}t}\geq0$ for all $t>0$. Thus, we finally obtain that $(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))>(0, 0, 0, 0, 0, 0)$ for all $t>0.$ This shows that $(S^{0}_{H}, I^{0}_{H}, M^{0}, S^{0}_{V}, I^{0}_{V}, P^{0})\notin M_\partial,$ which leads to a contradiction.

Similarly, when $I^{0}_{V}>0$ or $M^{0}>0$ or $P^{0}>0,$ we can also prove that $(S_{H}(t), I_{H}(t), \\M(t), S_{V}(t), I_{V}(t), P(t))>(0, 0, 0, 0, 0, 0)$ for all $t>0.$ This shows that $(S^{0}_{H}, I^{0}_{H}, M^{0}, \\ S^{0}_{V}, I^{0}_{V}, P^{0})\notin M_\partial,$ which leads to a contradiction. Therefore, we have $M_\partial\subseteq\{(S_{H}, 0, 0, S_{V}, 0, 0): S_{H}>0, S_{V}>0\}.$ Thus, we finally confirm that claim (13) holds.

Model (1) can be simplified as a subsystem $S_H'(t)=\Lambda_H-\mu_HS_H$ and $S_V'(t)=\Lambda_V(t)-\mu_V(t)S_V(t)$ on $\partial X_{0}.$ This shows that the map $P$ has a global attractor $M_{1}=\{(\frac{\Lambda_{H}}{\mu_{H}}, 0, 0, \widehat{S_{V}}(0), 0, 0)\}$ on $\partial X_{0}.$ It is clear that on $\partial X_{0},$ $\{M_{1}\}$ is isolated, invariant, and does not form a cycle. Therefore, conditions $(a)-(c)$ of $(C'_{2})$ in Theorem 1.3.1 in [57] hold.

Secondly, let $x_{0}=(S^{0}_{H}, I^{0}_{H}, M^{0}, S^{0}_{V}, I^{0}_{V}, P^{0})\in X_{0}.$ By the continuity of solutions with respect to the initial values, for any small enough $\varepsilon_{1}>0,$ there is a $\delta_{1}>0$, if $\|x_{0}-M_{1}\|\leq \delta_{1},$ we have

$\begin{array}{lll} \|u(t, x_{0})-u(t, M_{1})\| < \varepsilon_{1} \mbox{for all} t\in[0, \omega]. \end{array}$

(14)

where $u(t, x_{0})=(S_{H}(t), I_{H}(t), M(t), S_{V}(t), I_{V}(t), P(t))$ is the solution of model (1) with initial values $(S_{H}(0), I_{H}(0), M(0), S_{V}(0), I_{V}(0), P(0))=x_{0}$ and $u(t, M_{1})=(\frac{\Lambda_{H}}{\mu_{H}}, 0, 0, \widehat{S_{V}}(t), 0, 0).$

Now, we claim that

$\begin{array}{*{35}{l}} \lim {{\sup }_{n\to \infty }}\parallel {{P}^{n}}({{x}_{0}})-{{M}_{1}}\parallel \ge {{\delta }_{1}}. \\ \end{array}$

(15)

Suppose (15) is not true, then we have $\limsup_{n\rightarrow\infty}\parallel P^{n}(x_{0})-M_{1}\parallel < \delta_{1}$ for some $x_{0}\in X_{0}.$ For the sake of simplicity, we assume that

$\begin{array}{lll} \parallel P^{n}(x_{0})-M_{1}\parallel < \delta_{1}, \forall n\geq0. \end{array}$

(16)

From (14) we obtain $\|u(t, P^{n}(x_{0}))-u(t, M_{1})\| < \varepsilon_{1}$ for all $n\geq0$ and $t\in[0, \omega]$. Then, for any $t\geq0,$ let $t=n\omega+\widetilde{t}$, where $\widetilde{t}\in [0, \omega)$ and $n=[\frac{t}{\omega}]$ is the greatest integer less than or equal to $\frac{t}{\omega},$ by (16), we have $\|u(t, x_{0})-u(t, M_{1})\|=\|u(\widetilde{t}, P^{n}(x_{0}))-u(\widetilde{t}, M_{1})\| < \varepsilon_{1}.$ It follows that $0\leq I_{H}(t)\leq \varepsilon_{1}$, $0\leq I_{V}(t)\leq \varepsilon_{1}$, $0\leq M(t)\leq \varepsilon_{1}$ and $0\leq P(t)\leq \varepsilon_{1}$ for all $t\geq0.$ Then, by the first and fourth equations of model (1) we obtain that

$\begin{array}{lll} \displaystyle S'_{H}(t)\geq\Lambda_{H}-\varepsilon_{1}\lambda_{H}(t)S_{H}(t)-\mu_{H}S_{H}(t) \end{array}$

and

$\begin{array}{lll} \displaystyle S'_{V}(t)\geq\Lambda_{V}(t)-\varepsilon_{1}\lambda_{V}(t)S_{V}(t)-\mu_{V}(t)S_{V}(t) \end{array}$

for any $t\geq0.$ By the comparison principle, we have for any $t\geq0$ that $S_{H}(t)\geq U_{\varepsilon_{1}}(t)$ and $S_{V}(t)\geq V_{\varepsilon_{1}}(t)$, where $U_{\varepsilon_{1}}(t)$ and $V_{\varepsilon_{1}}(t)$ are the solutions of systems (10) and (11) with parameter $\varepsilon_{1}$ satisfying initial conditions $U_{\varepsilon_{1}}(0)=S^{0}_{H}$ and $V_{\varepsilon_{1}}(t)=S^{0}_{V}$, respectively.

Since systems (10) and (11) with parameter $\varepsilon_{1}$ have globally uniformly attractive positive $\omega-$periodic solutions $U^{*}_{\varepsilon_{1}}(t)$ and $V^{*}_{\varepsilon_{1}}(t)$, respectively, there exists a $t_{2}>0$ such that

$\begin{array}{rl} \displaystyle U_{\varepsilon_{1}}(t)\geq U^{*}_{\varepsilon_{1}}(t)-\frac{\theta}{2}, V_{\varepsilon_{1}}(t)\geq V^{*}_{\varepsilon_{1}}(t)-\frac{\theta}{2} \end{array}$

(17)

for all $t\geq t_{2}.$ From (12) and (17), we obtain $U_{\varepsilon_{1}}(t)>\frac{\Lambda_{H}}{\mu_{H}}-\theta$ and $V_{\varepsilon_{1}}(t)>\widehat{S_{V}}(t)-\theta$ for all $t\geq t_{2}$. Thus, we see that for all $t>t_{2}$

$\left\{ \begin{array}{lll} \displaystyle I_{H}'(t)\geq\lambda_{H}(t)(\frac{\Lambda_{H}}{\mu_{H}}-\theta)P-(\mu_{H}+\gamma_{H})I_{H}, \\ \displaystyle M'(t)=\lambda_{M}I_{H}-\mu_{M}(t)M, \\ \displaystyle I_{V}'(t)\geq\lambda_{V}(t)(\widehat{S_{V}}(t) -\theta)M-\alpha_{V}(t)I_{V}, \\ \displaystyle P'(t)=\lambda_{P}I_{V}-\mu_{P}(t)P. \end{array} \right.$

(18)

Considering the following auxiliary system:

$\left\{ \begin{array}{lll} \displaystyle \widetilde{I}_{H}'(t)=\lambda_{H}(t)(\frac{\Lambda_{H}}{\mu_{H}}-\theta)\widetilde{P}-(\mu_{H}+\gamma_{H})\widetilde{I}_{H}, \\ \displaystyle \widetilde M'(t)=\lambda_{M}\widetilde{I}_{H}-\mu_{M}(t)\widetilde{M}, \\ \displaystyle \widetilde{I}_{V}'(t)=\lambda_{V}(t)(\widehat{S_{V}}(t)-\theta)\widetilde{M}-\alpha_{V}(t)\widetilde{I}_{V}, \\ \displaystyle \widetilde{P}'(t)=\lambda_{P}\widetilde{I}_{V}-\mu_{P}(t)\widetilde{P}, \end{array} \right.$

that is

$\begin{array}{rl} \displaystyle \frac{dh}{dt}=(F(t)-V(t)+\theta N(t))h(t), h(t)=(\widetilde{I}_{H}(t), \widetilde{M}(t), \widetilde{I}_{V}(t), \widetilde{P}(t))^{T}. \end{array}$

(19)

By Lemma 2.1 in Zhang and Zhao [56], there exists a positive $\omega-$periodic function $\psi(t)=(\psi_{1}(t), \psi_{2}(t), \psi_{3}(t), \psi_{4}(t))^{T}$ such that $h(t)=e^{\mu t}\psi(t)$ is a solution of system (19), where $\mu=\frac{1}{\omega}\ln\rho(\Phi_{F-V+\theta N}(\omega)).$

Denote $L(t)=({I}_{H}(t), M(t), {I}_{V}(t), P(t))$. We can choose a small constant $\xi>0$ such that $L(t_{2})\geq\xi \psi(t_{2})$. From (18), the comparison principle implies that $L(t)\geq \xi e^{\mu t} \psi(t)$ for all $t>t_{2}$. By $\rho(\Phi_{F-V+\theta N}(\omega))>1,$ it follows that $\mu>0,$ then $\lim_{t\rightarrow\infty}L(t)=\infty,$ that is $\lim_{t\rightarrow \infty}(I_{H}(t), M(t), I_{V}(t), P(t))=(\infty, \infty, \infty, \infty)$, which is a contradiction with (16). Hence, claim (15) holds. This shows that $W^{s}(M_{1})\cap X_{0}=\emptyset.$ Therefore, condition (d) of $(C'_{2})$ in Theorem 1.3.1 in [57] holds. Consequently, by Theorem 1.3.1 and 3.1.1 in [57], $P$ is uniformly persistent with respect to $(X_{0}, \partial X_{0})$.

Lastly, since model (1) is periodic, we obtain that model (1) is uniformly persistent. From Remark 3.3, model (1) is also permanent. This completed the proof.

As a consequence of Theorem 5.1 and Remark 3.3, from the main results obtained in [38] on the existence of positive periodic solutions for general population dynamical systems, we have the following result.

Corollary 5.2. If $R_{0}>1$, then model (1) admits at least a positive $\omega$-periodic solution.

By Corollary 5.2, model (1) has at least one positive $\omega$-periodic solution when $R_{0}>1.$ Here an important issue is to ascertain its stability, such as local stability and global stability. Unfortunately, it is very difficult to establish the stability of the periodic solutions for model (1) as a six-dimensional system, especially in the use of Floquent theory in periodic linear systems and the Lyapunov method in stability theory. We will discuss these problems in the future research.

Example 5.3. Set $\Lambda_{H}=4\times10^4, \lambda_{H}(t)=0.8\times10^{-11}\times(1+4.5\sin(\frac{\pi t }{6}+4.5)),$$\mu_{H}=5\times10^{-4}, \alpha_{H}=6\times10^{-5}, \gamma_{H}=0.03, \lambda_{M}=20, \mu_{M}=30, \Lambda_{V}(t)=6000,$ $\lambda_{V}(t)= 40\times10^{-8}\times(1+4.5\sin(\frac{\pi t }{6}+4.5)), \mu_{V}(t)=2\times10^{-2}, \alpha_{V}(t)=0.02, \lambda_{P}=60$ and $\mu_{P}=0.2$ in model (1).

The initial values are $\vec{V}_{i}(0)= ((8-0.00005i)\times10^{7}, 100+30i, (1+0.5i)\times10^{4},$$(3+0.5i)\times10^{4}, (1+0.5i)\times10^{4}, (1+0.5i)\times10^{4}), i=1, 2, 3.$ The basic reproduction number $R_{0}=1.8048>1,$ then the model admits at least a positive $12$-periodic solution by Corollary 5.2. Time series of $I_{H}(t)$ are shown in Figure 4, we guess that the periodic solutions of model (1) are stable.

---

6. Application to the control of schistosomiasis in Hubei Province

The monthly reported human schistosomiasis data in Hubei from January 2008 to December 2014 from the China CDC [10] show a seasonal fluctuation, with a peak in late summer to early autumn and a nadir in late winter. We use model (1) to simulate these cases and estimate the values of parameters in $\lambda_{H}(t)$ and $\lambda_{V}(t)$ in (1) by means of the least-square fitting. With the help of the optimization toolbox $-Fminsearch$ in MATLAB, the numerical fitted curve of human schistosomiasis cases is shown in Figure 5. Sensitivity analysis of the main parameters and analysis of control and prevention measures are given in Figure 8 and Figure 9, respectively.

---

6.1. Estimation of model parameters

We explain the parameter values as follows:

(a) The average human lifespan is about 74 years in Hubei in 2008, which is obtained from the National Bureau of Statistics of China [34]. Thus, the monthly average death rate $\mu_{H}= \frac{1}{74\times12} = 1.126\times 10^{-3}.$ The natural death rate of miracidia is 0.9 per day [18], then the natural monthly death rate $\mu_{M}=0.9\times30=27.$ A portion $k=300$ of eggs leave the infective human body with the faeces or urine and enter the fresh water supply where they hatch into miracidia at a rate $\gamma_{1}=0.0232$ per day [5,33], so the monthly migration rate $\lambda_{M}=30k\gamma_{1}=30\times300\times0.0232=209.$ Existing data show that there are about 13.1% patients received treatments in 2008 [44], we select $\gamma_{H}=13.1%$.

Similarly, the natural death rate of snails is $0.000596$ per day [18,36], so the monthly rate $\mu_{V}=0.000596\times30=1.7880\times10^{-2}.$ Based on the daily data in [18,33,36], we obtain the monthly disease induced death rate of snails, migration rate and natural death rate of cercariae as $\alpha_{V}=0.012, \lambda_{P}=78$ and $\mu_{P}=0.12,$ respectively.

(b) The total number of population was $5.699\times10^7$ in Hubei at the end of 2007, and $55.70%$ of them lived in the countryside [34]. Since people who live in the countryside are vulnerable to infected water, the number of the initial susceptible people in January 2008 was $S_{H}^{0}=5.699\times10^7\times55.70%=3.174\times10^7$. The annual average birth rate is 9.19 over one thousand [34], then recruiting number of susceptible humans in January 2008 was $\Lambda_{H}=\frac{3.174\times10^7\times9.19\times10^{-3}}{12}=2.431\times10^4.$ The reported number of infected human cases in January 2008 was 124 [10], which was set as the initial infected human population $I_H^{0}=124$. The snail area was about $7.547\times10^8$ square meters, and the evaluated area of infected snails was about $2.632\times10^8$ square meters at the end of 2007 [8]. The number of average living snails in every square meters in Hubei was between 0.001 and 0.0082 [53], so the initial total number of susceptible and infected snails were estimated as $\Lambda^{0}_{V}=7.547\times10^8\times0.001=7.547\times10^5$ and $I^{0}_{V}=2.632\times10^8\times0.001=2.632\times10^5,$ respectively. After 3-6 months, the survival rate of snails were maintained more than 75% in the natural environment [27], we select the recruiting number of susceptible snails in January 2008 as $\Lambda_{V}=7.547 \times10^5 \times75% =5.660\times10^5.$ We derive the initial miracidia value $M^{0}$ and cercariae value $P^{0}$ reversely by the parameters $\lambda_M$ and $\lambda_P$, respectively, so $M^{0}=124\times209=25916$ and $P^{0}=2.632\times10^5\times78=2.052\times10^7.$

(c) Due to the lack of information about periodic functions $\Lambda_{V}(t)$, $\mu_{V}(t)$ and $\alpha_{V}(t),$ we choose $\Lambda_{V}(t)=\Lambda_{V}$, $\mu_{V}(t)=\mu_{V}$ and $\alpha_{V}(t)=\alpha_{V}+\mu_{V}$ as above constants.

(d) Parameters $a_{H}, \varphi_H,$ $a_{V}$ and $\varphi_V$ are obtained by fitting model (1) to data and are given in Table 1.

---

6.2. Numerical simulations

Based on these known parameter values, the fitted values are $a_{H}=8.0\times10^{-14}, \varphi_H=4.978$, $a_{V}=1.974\times10^{-8}$ and $\varphi_V=4.407,$ see Table 1. However, the coefficient of determination ($\mathbf{R^2}$) is a statistic analysis indicator of correlation between multiple variables, its value is between 0 and 1, the greater the value of $\mathbf{R^2}$ the better of fitting the model (see [29]). Let $y_{i} (i=1, 2\cdots, 84)$ be the recorded data (that is, the reported monthly human schistosomiasis cases from January 2008 to December 2014 from China CDC [10]), and $\hat{y}_{i} ( i=1, 2\cdots, 84)$ be the estimated data through model (1), then

${{\mathbf{R}}^{\mathbf{2}}}=\frac{\sum\nolimits_{1}^{84}{{{({{{\hat{y}}}_{i}}-\bar{y})}^{2}}}}{\sum\nolimits_{1}^{84}{{{({{y}_{i}}-\bar{y})}^{2}}}+\sum\nolimits_{1}^{84}{{{({{{\hat{y}}}_{i}}-\bar{y})}^{2}}}}=0.6463,$

where $\bar{y}=\frac{1}{84}\sum_{1}^{84}{y}_{i}.$ In general, numerical estimation results indicate that our model provides a relatively good match to the reported data. Moreover, we used our model to forecast the disease development trend after December 2014 in Figure 6 $(t\geq84)$. It can be seen that the number of predicted human schistosomiasis cases will fluctuate periodically and decrease meanwhile after January 2015.

The cure rate in 2008 was about $\gamma_{H}=13.1%$ [44], we calculate the basic reproduction number $R_{0}=1.2387>1.$ This shows that the number of infected humans $I_H(t)$ tends to a stable $12$-periodic solution by Theorem 5.1 and Corollary 5.2, see (a) in Figure 7. Data for 2014 show that $20.7%$ of people required for treatments [44]. We set the current cure rate as $\gamma_{H}=20.7%$, then $R_{0}=0.9971 < 1$, and $I_H(t)$ tends to 0 by Theorem 4.1, see (b) in Figure 7. Obviously, the cure rate plays an important role in control of human schistosomiasis. We believe that the human schistosomiasis can be relieved and conclude that the cases will be greatly reduced in the next few years in Hubei if we continue to increase treatments.

---

6.3. Sensitivity analysis

We carry out some sensitivity analysis to investigate the influence of parameters $\Lambda_{V}, \gamma_{H}, \lambda_{M}$ and $\lambda_{P}$ on $R_0$. From Figure 8, it is obvious that $R_0$ is a increasing function of $\Lambda_{V}, \lambda_{M}$ and $\lambda_{P}$, respectively, and a decreasing function of $\gamma_{H}.$ These indicate that $R_0$ can be less than $1$ in Hubei by reducing the recruiting of susceptible snails $\Lambda_{V}$, migration rate of cercaria $\lambda_{P}$, and migration rate of miracidia $\lambda_{M}$ to approximately less than $3000$, $30$ and $90$, respectively, or only increase the cure rate more than 25%.

Traditional strategies in controlling schistosomiasis include chemotherapy, health education, livestock chemotherapy, and snail control in risk areas [8], relying more on treating humans and animals. The sensitivity analysis demonstrate that these are all important measures to control schistosomiasis infection in Hubei, see Figure 8. and Figure 9.

From (a) in Figure 8, we can see that $R_0$ increases as the snails birth rate $\Lambda_V$ increases. (b) in Figure 9. shows that the number of human schistosomiasis cases decrease as the snails death rate $\mu_V$ increases. Thus, by reducing the month new born snails or killing the snails near residential areas as much as possible, the transmission cycle between humans and snails will be broken, and the number of human schistosomiasis can be decreased.

The simulation results (b) and (c) in Figure 8. and (c-e) in Figure 9. also show that virus migration rate $\lambda_{M}$ from humans to physical water, $\lambda_{P}$ from physical water to humans and the baseline transmission rate $a_{H}$ between humans and cercariae have an effect on reduction of the disease. These parameters can be decreased through managing feces and improving sanitation, which aim at killing worm eggs in feces, through hygiene education to give warnings not to swim, dig, water, mow grass, fish, laundry, wash dishes in the lakes with snails. These will ease the epidemic prevalence.

Increasing cure rate $\gamma_{H}$ can reduce human schistosomiasis cases (see (d) in Figure 8. and (f) in Figure 9.). It is recommended that treat groups at risk regularly with praziquantel [26]. Groups targeted for treatment include school-aged children in endemic areas, adults considered to be at risk in endemic areas, and people with occupations involving contact with infested water, such as fishermen, farmers, irrigation workers, women whose domestic tasks bring them in contact with infested water, and entire communities living in highly endemic areas.

---

7. Conclusion and discussion

It was observed that the number of schistosomiasis cases arrives at peak in late summer to early autumn, and reaches nadir in winter and spring in Hubei, Hunan, Anhui and other regions with similar geographic characteristics and environmental factors in China (see Figure 1), which display a seasonal pattern in these epidemic provinces. For the sake of simplicity and convenience, we only list data of three provinces in Figure 1. To investigate the human schistosomiasis transmission dynamics and explore effective control and prevention measures in these lake and marshland regions along the Yangtze River, we developed a nonautonomous model to describe seasonal schistosomiasis incidence rate by incorporating periodic transmission rates $\lambda_{H}(t)$ and $\lambda_{V}(t)$. We deduced the basic reproduction number $R_0$ and analyzed the dynamics of model (1) including global stability of the disease-free periodic solution and uniform persistence of the model. $R_0$ was calculated following the procedure of Wang and Zhao [50], that is $R_0=\rho(L),$ where $L$ is the next infection operator, which was developed from the original definition of Bacaer and Guernaoui [4].

Based on the data from China CDC [10], we used our model (1) to simulate the monthly infected human data from January 2008 to December 2014 in Hubei, the parameters in transmission functions $\lambda_{H}(t)$ and $\lambda_{V}(t)$ were estimated by least-square fitting, see Figure 5. We also predicted the general tendency of the disease by our model after January 2015 in Figure 6. It can be seen that the number of predicted human schistosomiasis cases will decrease and fluctuate periodically after January 2015. The Chinese Government once aimed to reach the criteria of transmission control threshold of less than 1% in the lake and marshland provinces and reach transmission interruption threshold in hilly provinces of Sichuan and Yunnan by the end of 2015 [45]. We believe this periodic model gave a relatively good match with the cases and current situation, see Figure 5 and Figure 6. Prevention and control strategies that we put forward theoretically for Hubei province were demonstrated in Figure 8 and Figure 9. Human schistosomiasis control is based on improving sanitation, hygiene education, snail control, and large-scale treatment of at-risk population groups.

Hubei, Hunan, Jiangxi, Jiangsu and Anhui provinces are located along the Yangtze River in central China, where climate changes clearly all the year round. Rivers and lakes water level rise in rainy spring and summer, then the area of snails increase, farmers and students have more chances to contact with contaminated water for agriculture work or routine life, so epidemics occur naturally in this period. With temperature declining in winter, people have less opportunities to contact with water. Sun et.al [37] estimated that the lowest critical temperature for the infection of snails with miracidia is $3.24\, ^{\circ}\mathrm{C}$, and deduced that the infection rate descends with temperature, so the epidemic outbreak descends in winter. In this way, the infections are subjected to environmental change, fluctuating from season to season [13].

To prevent and control the disease, the most basic work is to increase residents' knowledge of schistosomiasis, including harm of the disease, the transmission through feces of infected people and livestock, how people contract the disease (infection route), the snails as the intermediate host, etc. The best way to prevent infection is to avoid contacting infested water, and once infected, drug treatment of praziquantel is recommended [26].

In schistosomiasis epidemic seasons (April-October), schistosomiasis prevention and control work is very hard. In addition to routine control approaches such as chemotherapy, molluscicide treatment of snail habitats and health education, other major interventions including agriculture mechanization (phasing out the cattle for ploughing and other field work), prohibiting pasture in the grasslands along lake and rivers, building safe grassland for grazing, raising livestock in herds, improving sanitation through supplying safe water, constructing marsh gas pools, building lavatories and latrines, and providing fecal matter containers for fishermen's boats, etc., could decrease the prevalence of schistosomiasis to a very low level (see Figure 9. and [46,49]).

Duo to the increasing migration population and the changes in environments and diet habits, schistosomiasis rebounded in some areas where it had formerly been controlled or eliminated (see [60] and the references therein). Moreover, another threat is that traveling causes new infections of other species of schistosomiasis, for example, an increasing in the cases infected with S. haematobium or S. mansoni is reported in those returning to China after the China-aided projects in Africa and labor services export to Africa [47]. So highly sensitive surveillance and response system for those from overseas is necessary.

The model we set up is used to study the transmission dynamics and control of schistosomiasis in the lake and marshland areas. For mountainous regions, such as Sichuan and Yunnan provinces, the corresponding model needs further research. It is widely acknowledged that the transmission processes of S. japonica is considerably more complex in comparison to other schistosome species because its definitive hosts include more than 40 animal reservoirs, such as cattle, dogs, pigs and rodents [25]. The model should include the role of these hosts. We leave these for future consideration.

---

Acknowledgments

We are grateful to Dr. Daozhou Gao and two anonymous reviewers for their valuable comments and suggestions that greatly improved the presentation of this paper. This research was partially supported by the National Natural Science Foundation of China [11271312,31670656,11501498] and National Science Foundation (DMS-1412454).

---