A MATHEMATICAL MODEL FOR THE SEASONAL TRANSMISSION OF SCHISTOSOMIASIS IN THE LAKE AND MARSHLAND REGIONS OF CHINA

Schistosomiasis, a parasitic disease caused by Schistosoma Japonicum, is still one of the most serious parasitic diseases in China and remains endemic in seven provinces, including Hubei, Anhui, Hunan, Jiangsu, Jiangxi, Sichuan, and Yunnan. The monthly data of human schistosomiasis cases in Hubei, Hunan, and Anhui provinces (lake and marshland regions) released by the Chinese Center for Disease Control and Prevention (China CDC) display a periodic pattern with more cases in late summer and early autumn. Based on this observation, we construct a deterministic model with periodic transmission rates to study the seasonal transmission dynamics of schistosomiasis in these lake and marshland regions in China. We calculate the basic reproduction number R0, discuss the dynamical behavior of solutions to the model, and use the model to fit the monthly data of human schistosomiasis cases in Hubei. We also perform some sensitivity analysis of the basic reproduction number R0 in terms of model parameters. Our results indicate that treatment of at-risk population groups, improving sanitation, hygiene education, and snail control are effective measures in controlling human schistosomiasis in these lakes and marshland regions.

Abstract. Schistosomiasis, a parasitic disease caused by Schistosoma Japonicum, is still one of the most serious parasitic diseases in China and remains endemic in seven provinces, including Hubei, Anhui, Hunan, Jiangsu, Jiangxi, Sichuan, and Yunnan. The monthly data of human schistosomiasis cases in Hubei, Hunan, and Anhui provinces (lake and marshland regions) released by the Chinese Center for Disease Control and Prevention (China CDC) display a periodic pattern with more cases in late summer and early autumn. Based on this observation, we construct a deterministic model with periodic transmission rates to study the seasonal transmission dynamics of schistosomiasis in these lake and marshland regions in China. We calculate the basic reproduction number R 0 , discuss the dynamical behavior of solutions to the model, and use the model to fit the monthly data of human schistosomiasis cases in Hubei. We also perform some sensitivity analysis of the basic reproduction number R 0 in terms of model parameters. Our results indicate that treatment of at-risk population groups, improving sanitation, hygiene education, and snail control are effective measures in controlling human schistosomiasis in these lakes and marshland regions.  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 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 Λ 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 [54]) to describe the transmission rates from cercariae to human and from miracidia to snails, respectively, where positive constants a H , b H and ϕ H represent the human baseline transmission rate, its magnitude of forcing and the initial phase, respectively, positive constants a V , b V and ϕ V in λ V (t) have the similar meanings as constants in λ H (t). We choose bilinear incidence rates (density-dependent or mass action type) λ H (t)S H P and λ V (t)S V M (see [17]).
(5) It is clear that the snail population is seasonally changed in reality, the recruited rate Λ V (t), natural death rate µ V (t) and disease induced death rate α 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 γ H .
(8) We further assume that the human natural death rate µ H , miracidia natural death rate µ M and cercariae natural death rate µ P are all positive constants. Miracidia migration rate λ M from human to snail and cercaria migration rate λ P from snail to human are also supposed to be positive constants.
The model is described by the following system of ordinary differential equations: 3. Basic properties. We denote ω = 12 months. Based on the biological background of model (1), we only consider solutions of model (1) starting at t = 0 with initial values: When I H = 0, M = 0, I V = 0 and P = 0, model (1) has a unique diseasefree periodic solution E 0 = ( S H , 0, 0, S V (t), 0, 0), where S H = Λ H µ H , and S V (t) is the globally asymptotically stable positive ω-periodic solution of equation 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 Let Y (t, s) be the 4 × 4 matrix solution of the following initial value problem where I is the 4 × 4 identity matrix. Further, let C ω be the ordered Banach space of all ω−periodic continuous functions from R to R 4 with maximum norm · and positive cone C + ω := {φ ∈ C ω : φ(t) ≥ 0, ∀t ∈ R}. Suppose φ(s) ∈ C + ω is the initial distribution of infectious individuals, then F (s)φ(s) is the rate of new infection produced by the infectious individuals who were introduced at time s, and Y (t, s)F (s)φ(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 ≥ s.
is the distribution of accumulative new infections at time t produced by all those infected individuals φ(s) introduced at time previous to t. Then, we define a linear operator L : C ω → C ω as follows 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 = ρ(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).
On the positivity and boundedness of solutions of model (1) with nonnegative initial conditions (2), we have the following results.
is nonnegative for all t ≥ 0 and ultimately bounded. In particular, if S 0 Proof. In fact, by the continuous dependence of solutions with respect to initial values, we only need to prove that when S 0 Clearly, m(0) > 0. Assuming that there exists a t 1 > 0 such that m(t 1 ) = 0 and which leads to a contradiction. Similar contradictions can be deduced in the cases of m( 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 ε > 0 there is a T 0 > 0 such that is bounded for t > 0 and by the arbitrariness of ε we also have From the forth and fifth equations of model (1) . By the comparison principle and Lemma 1 in [40], we can obtain 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 . This completes the proof.
Lemma 3.2 implies that Ω is a positively invariant set with respect to model (1).
Proof. By considering the linearization system, we can prove that E 0 is locally asymptotically stable when R 0 < 1, which is equivalent to ρ(Φ F −V (ω)) < 1 by Lemma 3.1. We can choose a small enough positive constant ε such that Then for all t > t 1 , from model (1) we obtain that Considering the following auxiliary system: By Lemma 2.1 in Zhang and Zhao [56], it follows that there exists a positive . We can choose a small enough positive constant η > 0 such that J(t 1 ) ≤ ηφ(t 1 ). Then, from (7) and the comparison principle, we have J(t) ≤ ηe µt φ(t) for all t > t 1 . By ρ(Φ F −V +εN (ω)) < 1, it follows that µ < 0, then lim t→∞ J(t) = 0, which in turn implies that lim t→∞ (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→∞ S H (t) = Λ H µ H and lim t→∞ S V (t) = S V (t). Therefore, E 0 is globally attractive when R 0 < 1. This completes the proof. (1) reduces to an autonomous case. Using the method given by van den Driessche and Watmough [41], we obtain the corresponding basic reproduction numberR As a corollary of Theorem 4.1, we have the following result.
is globally asymptotically stable.
Model (1) can be simplified as a subsystem S on ∂X 0 . This shows that the map P has a global attractor It is clear that on ∂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.
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. By Corollary 5.2, model (1) has at least one positive ω-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.  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 λ H (t) and λ V (t) in (1) by means of the least-square fitting. With the help of the optimization toolbox −F minsearch 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.

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 µ H = 1 74×12 = 1.126 × 10 −3 . The natural death rate of miracidia is 0.9 per day [18], then the natural monthly death rate µ M = 0.9 × 30 = 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 γ 1 = 0.0232 per day [5,33], so the monthly migration rate λ M = 30kγ 1 = 30×300×0.0232 = 209. Existing data show that there are about 13.1% patients received treatments in 2008 [44], we select γ H = 13.1%.
Similarly, the natural death rate of snails is 0.000596 per day [18,36], so the monthly rate µ V = 0.000596×30 = 1.7880×10 −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 α V = 0.012, λ P = 78 and µ P = 0.12, respectively.  [18], [36] (b) The total number of population was 5.699 × 10 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 0 H = 5.699 × 10 7 × 55.70% = 3.174 × 10 7 . The annual average birth rate is 9.19 over one thousand [34], then recruiting number of susceptible humans in January 2008 was Λ H = 3.174×10 7 ×9.19×10 −3 12 = 2.431 × 10 4 . The reported number of infected human cases in January 2008 was 124 [10], which was set as the initial infected human population I 0 H = 124. The snail area was about 7.547 × 10 8 square meters, and the evaluated area of infected snails was about 2.632 × 10 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 Λ 0 V = 7.547 × 10 8 × 0.001 = 7.547 × 10 5 and I 0 V = 2.632 × 10 8 × 0.001 = 2.632 × 10 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 Λ V = 7.547×10 5 ×75% = 5.660×10 5 . We derive the initial miracidia value M 0 and cercariae value P 0 reversely by the parameters λ M and λ P , respectively, so M 0 = 124 × 209 = 25916 and P 0 = 2.632 × 10 5 × 78 = 2.052 × 10 7 .
(d) Parameters a H , ϕ H , a V and ϕ V are obtained by fitting model (1) to data and are given in Table 1.  Table 1. However, the coefficient of determination (R 2 ) is a statistic analysis indicator of correlation between multiple variables, its value is between 0 and 1, the greater the value of R 2 the better of fitting the model (see [29]). Let y i (i = 1, 2 · · · , 84) be the recorded data (that is, the reported monthly human schistosomiasis cases from January 2008 to December 2014 from China CDC [10]), andŷ i (i = 1, 2 · · · , 84) be the estimated data through model (1), then = 0.6463, whereȳ = 1 84 84 1 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 ≥ 84). It can be seen that the number of predicted human schistosomiasis cases will fluctuate periodically and decrease meanwhile after January 2015.  The parameter and initial values are the same as in Figure 5.
The cure rate in 2008 was about γ 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 γ 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.  Table 1. 6.3. Sensitivity analysis. We carry out some sensitivity analysis to investigate the influence of parameters Λ V , γ H , λ M and λ P on R 0 . From Figure 8, it is obvious that R 0 is a increasing function of Λ V , λ M and λ P , respectively, and a decreasing function of γ H . These indicate that R 0 can be less than 1 in Hubei by reducing the recruiting of susceptible snails Λ V , migration rate of cercaria λ P , and migration rate of miracidia λ 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 Λ V increases. (b) in Figure 9. shows that the number of human schistosomiasis cases decrease as the snails death rate µ 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 λ M from humans to physical water, λ 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  Table 1. grass, fish, laundry, wash dishes in the lakes with snails. These will ease the epidemic prevalence.
Increasing cure rate γ 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 λ H (t) and λ 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 = ρ(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 λ H (t) and λ V (t) were estimated by leastsquare 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 • 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.