Dynamic evolution of schistosomiasis distribution under different control strategies: Results from surveillance covering 1991–2014 in Guichi, China

Background Since the founding of the China, the Chinese government, depending on the changing epidemiological situations over time, adopted different strategies to continue the progress towards elimination of schistosomiasis in the country. Although the changing pattern of schistosomiasis distribution in both time and space is well known and has been confirmed by numerous studies, the problem of how these patterns evolve under different control strategies is far from being understood. The purpose of this study is, therefore, to investigate the spatio-temporal change of the distribution of schistosomiasis with special reference to how these patterns evolve under different control strategies. Methodology / Principal findings Parasitological data at the village level were obtained through access to repeated cross-sectional surveys carried out during 1991–2014 in Guichi, a rural district along the Yangtze River in Anhui Province, China. A hierarchical dynamic spatio-temporal model was used to evaluate the evolving pattern of schistosomiasis prevalence, which accounted for mechanism of dynamics of the disease. Descriptive analysis indicates that schistosomiasis prevalence displayed fluctuating high-risk foci during implementation of the chemotherapy-based strategy (1991–2005), while it took on a homogenous pattern of decreasing magnitude in the following period when the integrated strategy was implemented (2006–2014). The dynamic model analysis showed that regularly global propagation of the disease was not present after the effect of proximity to river was taken into account but local pattern transition existed. Maps of predicted prevalence shows that relatively high prevalence (>4%) occasionally occurred before 2006 and prevalence presents a homogenous and decreasing trend over the study area afterwards. Conclusions Proximity to river is still an important determinant for schistosomiasis infection regardless of different types of implemented prevention and control strategies. Between the transition from the chemotherapy-based strategy to the integrated one, we noticed a decreased prevalence. However, schistosomiasis would remain an endemic challenge in these study areas. Further prevention and control countermeasures are warranted.


Methodology / principal findings
Parasitological data at the village level were obtained through access to repeated cross-sectional surveys carried out during 1991-2014 in Guichi, a rural district along the Yangtze River in Anhui Province, China. A hierarchical dynamic spatio-temporal model was used to evaluate the evolving pattern of schistosomiasis prevalence, which accounted for mechanism of dynamics of the disease. Descriptive analysis indicates that schistosomiasis prevalence displayed fluctuating high-risk foci during implementation of the chemotherapy-based strategy (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005), while it took on a homogenous pattern of decreasing magnitude in the following period when the integrated strategy was implemented (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). The dynamic model analysis showed that regularly global propagation of the disease was not present after the effect of proximity to river was taken into account but local pattern transition existed. Maps of predicted prevalence shows that relatively high prevalence (>4%) Introduction Schistosomiasis japonica, caused by the trematode worm Schistosoma japonicum [1], is responsible for human and animal infections in China, the Philippines and parts of Indonesia [2]. It has a documented history of at least 2100 years along the Yangtze River Basin and once had considerable public health and economic significance in the regions [3]. After the founding of the China, the central government made great strides toward fighting schistosomiasis. The initial national schistosomiasis control program began in 1950s [4]; it mainly based on control of the intermediate snail host that spreads the infection to the definitive mammal host, but chemotherapy was also used though the drugs available at the time were not very effective. Strong political will and sufficient long-term financial support [5] resulted in a remarkable decline in both prevalence and intensity of disease [6,7]. Despite the success of mass treatment using praziquantel trough the 10-year World Bank Loan Project (WBLP) that started in 1992, schistosomiasis rebounded in the early 2000s, not only because of the termination of the WBLP [8,9], but also due to common floodings along the Yangtze River basin [10] and other natural, social and economic factors [11,12]. A new national strategy integrating chemotherapy, snail control, health education, improved sanitation and access to safe water, was initiated in 2005 [13,14]. This change of strategy achieved success almost immediately in controlling S. japonicum transmission, as it again resulted in lower numbers of infections, both in domestic animals and humans [15,16].
Uneven schistosomiasis distribution is a well-known phenomenon, even in limited geographical settings [17,18]. Numerous studies have identified this focal pattern [19,20,21,22,23], which is representative for the distribution of schistosomiasis cases in both time and space. However, the question how this pattern evolves is far from understood. The role of, and the causes leading to, epidemiological pattern transition that results in spread of diseases have only recently become a major concern for epidemiological studies [24]. Considering the substantial changes the national strategy of schistosomiasis control in China has undergone, we assume that corresponding changes in the local environmental conditions in the endemic areas have taken place, resulting in pattern transitions of the disease. Understanding and quantifying these transitions would help assessing the effectiveness of control strategies [25,26].
The goal of this study was therefore to investigate the evolution of schistosomiasis distributions patterns under the two latest stages of the national program for schistosomiasis control.

Ethics statement
Approval for oral consent and other aspects of the surveys was granted by the Ethics Committee of Fudan University (ID: IRB#2017-09-0635). Written informed consent was also obtained from adult participants (�18 years). Parental consent was sought for participants aged 5 to 17 years. All the individual information of participants were anonymized by deleting the personal identifiers (such as name, parent name, home address, and telephone number) for the purpose of protecting patients' privacy. All information and consent procedures were conducted in Guichi, Anhui Province, China.
We first conducted an exploratory data analysis of schistosomiasis prevalence in Guichi, a rural district belonging to the city of Chizhou, Anhui Province, using annual village-level prevalence data for the period 1991-2014, and then built a dynamic model to investigate the evolution of schistosomiasis epidemiology during this period.

Study area
Guichi remains one of the most highly endemic areas of schistosomiasis in China and became, for this reason, one of the 20 national sentinel surveillance sites in 2000 [27]. Situated in the middle-lower reaches of the Yangtze River Basin, Guichi spans approximately 2,500 km 2 and is composed of 164 villages. Besides the Yangtze River, another major river, the Qiupu, runs across the region from the Midwest to the North (Fig 1). The humid subtropical, monsoon climate, with an average annual temperature about 16˚C and annual rainfall around 1,600 mm, provide an ideal environment for the survival of O. hupensis, the intermediate snail host [20].

Parasitological data
S. japonicum infection prevalence data for 1991-2014 were provided by the local anti-schistosomiasis station in Guichi. These data were originally collected through village-based field surveys that required residents aged 5 to 65 years to participate in. Indirect Hemagglutination Assay (IHA) test was first used to screen the individuals according to the manufacturer's instructions and then the parasitological test was applied to the seropositive individuals by reading three Kato-Katz [28] thick smears from one stool specimen. All those found to be eggpositive were regarded as having an active S. japonicum infection and therefore treated with praziquantel. For each year, we removed the villages for statistically reliable analysis, in which the number of individuals participating in the serological test are less than 100. With this criteria, the number of sample village ranged from 64 in 1995 to 128 in 2011.

Environmental data
Water body data, including the Yangtze River and Qiupu River, were extracted from the World Wildlife Fund's Conservation Science Data Sets (http://worldwildlife.org). For each sample village, the Euclidian distance (unit: kilometer) to the nearest river was calculated by using ArcGIS software version 10.0 (ESRI Inc., Redlands, CA, USA).

Statistical analysis
With these large spatio-temporal schistosomiasis prevalence data, we carried out an exploratory data analysis to give hints to follow-up spatio-temporal dynamic modeling. A Hovmoller plot, a two-dimensional space-time visualization in which space is collapsed onto one dimension and where the second dimension denotes time [29], was employed.

Dynamic spatio-temporal model
In order to model the dynamic pattern of schistosomiasis, we assume that there is a true unobserved spatio-temporal process hidden behind the yearly observed prevalence in each village, which is incorporated into the framework of a hierarchical model. The basic representation of the hierarchical model is typically composed of data level (whose conditional probability distribution given processes and parameters is independent), process level (which determines change of data level given parameters), and parameter level (which are contained in the previous two levels).
Data level. This describes the relationship between the observations and the latent spatiotemporal process. Specifically, we could present the data model as: where Z t corresponds to the data vector (which includes arbitrary spatial locations) at time t, α is the intercept, X kt are fixed covariates that specified as environmental or socio-economic factors and β kt is the kth covariate's coefficient, Y t is the latent spatio-temporal process, with a linear mapping here, H t that connects the data to the latent process, and finally ε t is a timevarying (but statistically independent in time) continuous mean-zero Gaussian process which is typically measurement error. Although it is in principle possible to parameterize the mapping matrix H t and/or estimate it directly in some cases, we shall typically assume that it is known as a simple incidence matrix (i.e., a matrix of ones and zeros) which depends on the dimension of Y t . In addition, an important assumption is that the data are independent (in time) when conditioned on the latent process Y t and parameters α t , H t , and ε t . Process level. The latent process here is presumed to follow a first-order linear Markovian spatio-temporal process which assumes that the value of the process at a given location at the present time is made up of a weighted combination of the process throughout the spatial domain at previous times, plus an additive, Gaussian, spatially coherent "innovation". This is perhaps best represented in a continuous-spatial context through an integro-difference equation (IDE) as follows: where m(s,x;θ p ) is a transition kernel, depending on parameters θ p that specify "redistribution weights" for the process time over the spatial domain D s , and ω t (�) is a time-varying (but statistically independent in time) continuous mean-zero Gaussian spatial process independent of Y t −1 (�). Elements {m(s,x;θ p )} consist of the "so called" transition matrix M. We assume m(s,x;θ p ) to be a Gaussian-shape kernel as a function of x relative to the location s: where the kernel amplitude is given by θ p,1 , the length-scale (variance) parameter θ p,2 corresponds to kernel scale (aperture) parameter, the mean (shift) parameter θ p,3 and θ p,4 correspond to a shift of kernel relative to location s, and 1 and 2 (� corresponds to s and x) denotes longitude and latitude coordinate, respectively. The kernel parameters θ p in formula (3) are spatially varying which would result in complex dynamic behavior, but they can be spatially invariant as well (more details in S1 Text). In addition, we assume here that θ p dose not vary with time. In general, R D s mðs; x; y p Þdx < 1 is needed for the process to be stable (i.e., nonexplosive) in time.
Parameter level. At this level, we summarize all the parameters to be estimated in above models, which include the intercept α, the fixed-effect coefficient for proximity to river β, and kernel parameters θ p , and variances of the measurement error ε t and the latent process error ω t . All the parameters are estimated using the expectation-maximization (EM) algorithm [30] and inference for the latent process Y t is implemented using kalman filtering and smoothing [31]. Parametric significance (p-value) was calculated using a permutation approach with 99 random simulations of the observed prevalence.
For model validation, we utilize a k-fold cross-validation (k is assumed to be 10 in our study) in terms of prediction performance as follows to determine whether spatially invariant or spatially varying kernel can capture the dynamic behavior of the hidden process Y t : where {Z v (s i ;t j )} are samples of observations by randomly sampling 10% of all the spatial and temporal observation data and fẐ v ðs i ; t j Þg are prediction by the dynamic model based on the rest of the observation. All model fitting was performed using R software, specifically the IDE package (R Development Core Team, 2013). Finally, the latent process, in the form of predicted schistosomiasis prevalence over the study area (displayed as maps), were produced using the optimal model; the spatial resolution for prediction was 1 km. A conceptual framework of the data analysis is shown in Fig 2.

Fig 3 shows
Hovmoller plots for the schistosomiasis prevalence in the study area between 1991 and 2014. Spatially, we see in the left panel that the prevalence increases generally as we move northwards within the county (corresponding to higher values for the latitudes), particularly as we reach the latitude 30.71 degree, but there is no trend with respect to the longitude (right panel) before 2006. Temporally, we find out a potential temporal trend in prevalence with decreasing longitude in the right panel (e.g., the yellow color at the right bottom in the plot slants at a left-top direction from 1991 to 1996 and the same phenomena from 2007 to 2010). In addition, the prevalence seems to be increasing since 2005 but then decrease again after 2010. Fig 4 illustrates the average MSPE of predictions with the spatially invariant and spatially varying kernels at each year, respectively. As shown in this figure, the spatially invariant kernel θ p fit spatio-temporal prevalence data better (less MSPE at almost each year and the whole study period). The EM estimates for all parameters in spatially invariable model are presented in Table 1. There was a statistically negative correlation (p-value = 0.01) between schistosomiasis prevalence and proximity to river. The very small value of shift parameters (θ p,3 = 3.68e (-03) and θ p,4 = 0.38e(-03)) indicates that there was almost no presence of shift in the disease data. Table 2 shows that eigenvalues of the transition matrix M of the spatially invariant model are all less than one, indicating stationary patterns.

Discussion
A dynamic spatio-temporal model was used to investigate the evolution of the schistosomiasis pattern in a high endemic area of schistosomiasis in China. Compared to the descriptive approach, used in most previous studies [23,32,33,34,35,36] which fits means and covariances to spatio-temporal data, we believe that our model represents a more dynamic approach that accounts for prior knowledge (i.e., shifting trend of prevalence shown in Fig 3) to capture the evolution of current and future spatial fields from the characteristics of past ones. A description of this evolution of schistosomiasis prevalence should be helpful for the understanding of Hovmoller plots (Fig 3) for both the latitude and longitude coordinates for the schistosomiasis prevalence data provide empirical recognition of the formation of disease distribution patterns as well as the prior knowledge (i.e. propagation of prevalence) for our dynamic  modeling. In our case, the prevalence during 1991-2005 shows a fluctuating pattern with higher prevalence present at higher latitudes (located further North), which is possibly due to the Yangtze and Qiupu Rivers in the northern part of the study area (Fig 1), and the potential temporal trends in prevalence with decreasing longitude, probably suggesting propagation through time from the eastern longitudes towards the western longitudes. During 2006-2014, the prevalence show a fairly constant pattern with foci at western longitudes which is probably related to the two rivers (Fig 1) which indicates that there is no shift of schistosomiasis prevalence during this period.
The empirical recognition of the disease pattern formation is supported by results of our dynamic modelling. We used a spatially invariant model to fit spatio-temporal prevalence data after carrying out the model validation. We find that there is nearly no shift present in the disease pattern after the effect of proximity to river is taken into account, suggesting that the effect of proximity to river can capture the spatio-temporal variation of schistosomiasis prevalence data and that water contact is possibly still an important determinant for schistosomiasis as a result of agricultural activities and fishing [20]. The stationary pattern means the disease's pattern remain unchanged/stable over time and often exhibit intermittent areas of high density of cases that favor persistence of the disease. This pattern formation and the accompanying transition is more evident in the annual prevalence predictions (Fig 5). During 1991-2005, schistosomiasis prevalence showed heterogeneities, i.e., limited number of geographic foci of high prevalence (except 2005) occurred intermittently over the study area, whereas homogeneities were present during 2006-2014 as large patches of relatively low prevalence constantly occurred.
Although no evidence of regularly global propagation of the disease was seen, local pattern transition existed which can be potentially explained by the change in schistosomiasis control strategies during the study period [9]. In the early 1990s, the local government of Guichi implemented the 10-year WBLP according to the national strategy for schistosomiasis control, which were largely based on chemotherapy [37]. After conclusion of the WBLP in 2001, the local government continued to carry out this strategy [20] but the disease rebounded afterwards due to limited financial support. This rebound can be seen from maps in 2003-2005 (Fig 5). Since 2006, a revised strategy was implemented according to the new national control strategy [38], which aimed at reducing the role of bovines and humans as infection sources. Specific interventions were as follows: (1) all bovines were removed from the study area (e.g., some farm bovines were replaced with tractors, or transferred to non-endemic regions; some bovines were killed); (2) improved sanitation facilities (e.g., sanitary lavatories and piped water were installed); (3) a health education program was initiated, focusing on avoiding snailinfected areas and associated river water; (4) modifying snail microhabitat through agricultural practices (e.g., modifying crop structures in snail-inhabited areas to reduce water use). This integrated strategy substantially changed the local physical environment as well as human behavior. We, therefore, believe this change of control strategy resulted in the local pattern transition of the disease presented here. Although the integrated strategy is more effective than chemotherapy alone in reducing prevalence (as shown in Fig 5), our analysis (Table 1) indicates that it is more likely for the disease to persist under this strategy. Two limitations need discussion. One is that our dynamic model only considered the effect of proximity to river (in formula (1)). We should have used some other factors associated with the disease such as climatic data (e.g., precipitation and temperature) and socio-economic data (e.g., GDP and medical resources). We assumed the climatic factors geographically changed little over the study area (i.e., only at a county level) and the climatic conditions during the study period have not changed substantially [39]. In addition, socio-economic data so far are not available at the village-level. Furthermore, the IDE used in our study models the "residual" after discounting what the covariates explain. If the covariates explain the spatio-temporal variation within the response variable, no trends are left in the "residual" which is the case in our study. Nevertheless, inputting those factors into our dynamic model would be warranted in a future study. The other limitation is that the specificity of serological assays and the sensitivity of stool examination tests are not perfect due to the generally low level of infection, and this could result in an imperfect measure of prevalence data of schistosomiasis. Dynamic modeling with diagnostic errors is also our further study direction.
In summary, we built a hierarchical dynamic spatio-temporal model to investigate the evolution of schistosomiasis pattern in an area known for remaining high transmission. Based on results of pattern analysis, we conclude that local pattern transition was probably due to the change in the two national schistosomiasis control strategies and the integrated strategy is more effective than the chemotherapy-based one, with a caution that the disease would be a long-term endemic. Methodologically, the merit of our study is that it provides a good example for modelling high-dimensional spatio-temporal data in an epidemiological study by using kernel functions to reduce data dimension. As more local and national surveillance systems are established, there would be more large-scale datasets available in public health. Undoubtedly, our study augments the methodology for dealing with such problems.
Supporting information S1 Text. Supporting text containing details of the kernel function. (DOCX)