The basic reproductive ratio of Barbour's two-host schistosomiasis model with seasonal fluctuations.

BACKGROUND: Motivated by the first mathematical model for schistosomiasis proposed by Macdonald and Barbour's classical schistosomiasis model tracking the dynamics of infected human population and infected snail hosts in a community, in our previous study, we incorporated seasonal fluctuations into Barbour's model, but ignored the effect of bovine reservoir host in the transmission of schistosomiasis. Inspired by the findings from our previous work, the model was further improved by integrating two definitive hosts (human and bovine) and seasonal fluctuations, so as to understand the transmission dynamics of schistosomiasis japonica and evaluate the ongoing control measures in Liaonan village, Xingzi County, Jiangxi Province.
METHODS: The basic reproductive ratio R 0 and its computation formulae were derived by using the operator theory in functional analysis and the monodromy matrix theory. The mathematical methods for global dynamics of periodic systems were used in order to show that R 0 serves as a threshold value that determines whether there was disease outbreak or not. The parameter fitting and the ratio calculation were performed with surveillance data obtained from the village of Liaonan using numerical simulation. Sensitivity analysis was carried out in order to understand the impact of R 0 on seasonal fluctuations and snail host control. The modified basic reproductive ratios were compared with known results to illustrate the infection risk.
RESULTS: The Barbour's two-host model with seasonal fluctuations was proposed. The implicit expression of R 0 for the model was given by the spectral radius of next infection operator. The R 0 s for the model ranged between 1.030 and 1.097 from 2003 to 2010 in the village of Liaonan, Xingzi County, China, with 1.097 recorded as the maximum value in 2005 but declined dramatically afterwards. In addition, we proved that the disease goes into extinction when R 0 is less than one and persists when R 0 is greater than one. Comparisons of the different improved models were also made.
CONCLUSIONS: Based on the mechanism and characteristics of schistosomiasis transmission, Barbour's model was improved by considering seasonality. The implicit formula of R 0 for the model and its calculation were given. Theoretical results showed that R 0 gave a sharp threshold that determines whether the disease dies out or not. Simulations concluded that: (i) ignoring seasonality would overestimate the transmission risk of schistosomiasis, and (ii) mollusiciding is an effective control measure to curtail schistosomiasis transmission in Xingzi County when the removal rate of infected snails is small.

Background

Schistosomiasis japonica, a parasitic disease caused by Schistosoma japonicum, has been in existence in the People’s Republic of China (P. R. China) for over 2000 years with considerable public-health and economic significance [1, 2]. A large-scale national schistosomiasis control programme was initiated in the mid-1950s [2, 3], when China’s population was approximately 600 million. An estimated 11.8 million people were infected with S. japonicum [4]. Sustained control efforts have contributed significantly to the dramatic reduction of both transmission intensity and schistosomiasis distribution in China in the past six decades [5-9]. Recent data have shown that progress toward schistosomiasis elimination encountered difficulties and setback due to high re-infection rates, particularly in lake region of the endemic areas [10]. For instance, schistosomiasis re-emerged shortly after the termination of The World Bank Loan Project (WBLP) at the end of 2001 [1, 5, 11]. In order to achieve the goal of schistosomiasis elimination in P. R. China, Chinese government strengthened the national schistosomiasis control programme in 2004. This made schistosomiasis control top priority along with the list of other communicable diseases such as HIV/AIDS, tuberculosis, and hepatitis B in China [4, 12]. Moreover, a revised strategy to effectively control schistosomiasis by using integrated measures in the national control programme has been implemented since 2005 [11].

Over 40 different species of wild and domestic animals have been identified as definitive hosts of S. japonicum [13]. Bovines are the major reservoirs for S. japonicum in the lake and marshland regions of southern China [14, 15]. A large number of schistosome-infected bovines are distributed in these regions and they excrete large quantities of S. japonicum eggs, the majority of which are deposited near or in the lake [14]. The daily faecal output from a water buffalo (~25 kg) has been estimated to be at least 100 times more than that (~250 g) from an individual human [16, 17]. It was reported that the overall prevalence of S. japonicum was 9.6 and 7.2% in water buffalo and cattle, respectively, in 1995 [18]. High prevalence of S. japonicum infection in bovine reservoir hosts is believed to be the major factor maintaining active transmission in certain areas. A cluster randomized intervention trial which was designed to compare the control (human treatment) and intervention (human and bovine treatment) in some villages concluded that the incidence of human S. japonicum infection is reduced with a decline in the infection rates of water buffaloes [19]. In order to remove bovines as a source of infection, several effective measures, including replacing cattle with farm machinery, isolating marshland, and prohibiting grazing in susceptible areas, have been implemented.

Mathematical modelling is a powerful tool to study the transmission dynamics of schistosomiasis [20]. It was first proposed by Macdonald [21]. Following this pioneering work, many mathematical models have been developed, all of which show great potential in aiding our understanding of the interplay of biology, transmission dynamics and control of schistosomiasis [22-27]. In 1996, Macdonald’s model was improved by Barbour [28]. It tracks dynamics of both infected human and snails in a community. Barbour’s model has played an important role in evaluating possible control strategies [29]. Traditionally, schistosomiasis models assume that all parameters are constant. A real world environment is obviously non-stationary, and should include seasonal variations in snail population. Infection rates vary seasonally due to natural factors (i.e. changes in moisture and temperature) and social factors (i.e. changes in contact rates). Therefore, it is more realistic to assume that the infection rates are periodic rather than constant. In our previous study [30], we constructed the Barbour’s single-host model with seasonal fluctuations (BSHSF model) and calculated the basic reproductive ratio to assess the effect of integrated control measures against schistosomiasis in Liaonan village, Xingzi County, Jiangxi Province. However, the impact of the bovine reservoir host on the transmission of schistosomiasis should not be ignored.

This study aims to further modify the Barbour’s two-host model with seasonal fluctuations (BTHSF model) and give an implicit expression of the basic reproductive ratio of schistosomiasis and computation method. Also, we will further establish some indexes that predict prevalence variation characteristics of schistosomiasis, and evaluate the prevention and control strategies for schistosomiasis.

Methods

Mathematical formulation

The Barbour’s single-host model with seasonal fluctuations (BSHSF model) in Gao et al. [30] was given by the following equations:

All variables of Eq. (1) are described in Table 1. The infection rate a(t) represents the rate at which a single definitive host becomes infected at unit density of infected snails at time t. The infection rate b(t) is the rate at which snails become infected at time t. The two infection rates can be calculated as follows (see [28]):

Table 1

Interpretation of the model (1)

Parameter	Interpretation
a(t)	the rate of incidence for a single definitive host at unit density of infected snails at time t
b(t)	the rate at which an infected definitive host causes snail infections at time t
g	the recovery rate for definitive host infections
Δ	the density of snails
Σ	the density of definitive hosts
μ	per capita removal rate of infected snails

Abbreviations: P, the prevalence of infection in the definitive host population; y, the proportion of infected snails

Here, α1 denotes the rate at which a host has contact with contaminated water per unit time, α2 denotes the density of cercariae, α3 denotes the probability that an encounter with cercaria leads to the establishment of an adult parasite, β 1 is the rate of egg-laying, β 2 is the probability of an egg developing into a miracidium, β 3 is the probability that a miracidium penetrates a snail, β 4 is the probability that miracidial penetration into an uninfected snail develops into infection.

Seasonality was taken into consideration by assuming that the infection rates a and b in model (1) were time-varying and periodic. a(t), b(t) were taken as the absolute value of sine functions with 365 days period, while the average values of functions a(t) and b(t) per year (365 days) were assumed to be equal to the values of parameters a and b in [31], respectively.

Furthermore, we considered the impacts of bovine reservoir host on the transmission of schistosomiasis. The BSHSF model (1) was modified and used to create the BTHSF model. The dynamics of the model are governed by the following differential equations:

All variables of Eq. (2) are described in Table 2, where a 1(t) denotes the rate of incidence for one human host at unit density infected by one infectious snail per unit time (one day) at time t, a 2(t) denotes the rate of incidence for one bovine host at unit density infected by one infectious snail per unit time (1 day) at time t, b 1(t) is the rate of incidence for one susceptible snail at unit density infected by one infected human per unit time (1 day) at time t, b 2(t) is the rate of incidence for one susceptible snail at unit density infected by one infected bovine per unit time (1 day) at time t. In this study, we assume that b 1 is proportional to b 1, i.e. b 2 = kb 1, where k is the coefficient ratio. According to Zhao et al. [32]:where A (i = 1,2) denote faeces output per human and bovine, B (i = 1,2) denote the epg of human and bovine, respectively, here epg is the number of eggs per gram of faeces. The number of eggs per bovine per day is more than 105 [33], and the number of eggs per human per day is about ~2800–3000 [34]. Furthermore, many effective measures, such as improved sanitation and health education are implemented to control the spread of schistosomiasis among human beings. In this paper, we choose k =50. We also take the infection rates as:

Table 2

Interpretation of the model (2)

Parameter	Interpretation
a 1(t)	the infection rate from snail to human at time t
a 2(t)	the infection rate from snail to bovine at time t
b 1(t)	the infection rate from human to snail at time t
b 2(t)	the infection rate from bovine to snail at time t
g 1	the recovery rate for human host infections
g 2	the recovery rate for bovine host infections
Δ	the density of snails
Σ1	the density of human hosts
Σ2	the density of bovine hosts
μ	per capita removal rate of infected snails

Abbreviations: P 1, prevalence of infection in the human host population; P 2, the prevalence of infection in the bovine host population; y, the proportion of infected snails

Since infection rates a and b (i = 1,2) denote the average values of a (t) and b (t) per year, respectively, we haveThus we can derive and

The basic reproductive ratio (R0)

For epidemiological models, the basic reproductive ratio, which is often introduced as a threshold parameter, is defined as the expected number of secondary infections produced by a single infective individual in a completely susceptible population during its entire infectious period [35]. R 0 has often been used to predict the trend of the disease transmission and also to assess the effects of control measures. Here, we give the implicit formula R 0 of model (2). It is expected that the infection will persist if R 0 is greater than one. As R 0 increases, the spreading rate of the disease also increases. This implies that more intervention efforts or improved control measures should be implemented. We have also shown that R 0 is a sharp threshold value which determines whether the disease dies out or not in Additional file 1: Appendix B. That is, if R 0 < 1, then the disease will die out, whereas if R 0 > 1, the disease will be endemic.

With reference to Bacaër [36], a biologically meaningful threshold value R 0 of BTHSF model (2) was derived by using operator theory in functional analysis and the monodromy matrix of linear periodic system theory. The numerical computation of the basic reproductive ratio was carried out using the mathematical programming language MATLAB 7.1. Also, sensitivity analysis was carried out in order to assess the impact of seasonal fluctuations on R 0

Study area and estimation of model parameters

Liaonan village, Xingzi County, Jiangxi Province was selected as the study area. The basic reproductive ratio was calculated based on the parameters in model (2). Determination of accurate model parameters was the priority when this threshold was applied. The details of the technique used in estimating the parameters have been described in our previous study [30]. In this paper, part of the parameters in model (2) were determined according to their biological significance. In addition, the remaining parameters were estimated based on the annual report surveillance data of Xingzi County from 2003 to 2010. The meanings of parameters g , Σ (i = 1,2), μ, ω are consistent with those reported in [30] (see Tables 1 and 2).

The parameters are stated as follows; g1 = 0.00093 per day, μ = 0.0055 per day, and g2 = 0.00183 per day (the data are given by Wu [31]). The density of snails, Δ(/m 2), and the density of definitive hosts, Σ1, Σ2(/m2), were estimated from annual report surveillance data (see Table 3).

Table 3

The data for Δ, Σ1, Σ2, , and from annual report data in the village of Liaonan, Xingzi County, Jiangxi Province, China

Parameter	2003	2004	2005	2006	2007	2008	2009	2010
Δ	27.997	33.9230	91.4602	133.5238	41.0615	17.4960	11.3121	2.4183
Σ1	0.0300	0.0199	0.0475	0.0167	0.0309	0.0178	0.0309	0.0176
Σ2	0.0013	0.0007	0.0020	0.0006	0.0010	0.0006	0.0013	0.0008
P¯_1	0.0671	0.0533	0.0631	0.0454	0.0454	0.0454	0.04530	0.0457
P¯_2	0.0283	0.1379	0.1756	0.1144	0.1084	0.0676	0.0639	0.0466
y¯	0.0002	0.00005	0.0001	0.0001	0.0002	0.0005	0.0008	0.0015

It is assumed that the coefficients of model (2) are constant. By setting the right-hand side of equations in model (2) to zero, the steady-state (equilibrium) values for the prevalence of infection in humans, bovines and snails were obtained [28] where t ( = b ∑/(μΔ) and t ( = a Δ/g i(i = 1, 2) are the transmission factors for definitive host i. By solving (3) and (4), we obtain the prevalence of infection in humans, bovines and snails from 2003 to 2010.

Results

Calculation of R0

Following the results of Bacaër [36] and Wang & Zhao [37], we show that R 0 is equal to the value of λ, where λ is the positive root of the equation

Hereand ρ(W(365, 0, λ)) is the spectral radius of matrix W(365, 0, λ).

The mathematical details that were used to derive the expression for the basic reproductive ratio can be found in Additional file 1: Appendices A, B (Derivation of R 0 and proof of main results).

Numerical simulation

Based on the annual report data for Xingzi County from 2003 to 2010, the equilibrium values (prevalence of infection in human host, bovine reservoir host and snail) are shown in Table 3. In view of the parameters’ values given above (see Table 3), the infection rates a 1, a 2, b 1 and b 2 of BTH model was derived from Eqs. (3) and (4), and are shown in Table 4.

Table 4

The calculated values of composite parameters a 1, a 2, b 1 and b 2

Parameters	2003	2004	2005	2006	2007	2008	2009	2010
a 1	0.0116	0.0309	0.0051	0.0024	0.0061	0.0050	0.0051	0.0120
a 2	0.0093	0.1726	0.0320	0.0128	0.0305	0.0151	0.0144	0.0241
b 1	0.0082	0.0015	0.0032	0.0260	0.0057	0.0171	0.0085	0.0075
b 2	0.4084	0.0754	0.1606	1.2996	0.2871	0.8565	0.4268	0.3725

In Additional file 1: Appendix A, we illustrate that R 0 is the positive root of the equation ρ(W(365, 0, λ)) = 1. Using the data in Tables 3 and 4, the values of R 0 for BTHSF model were calculated and shown in Table 5. Let R 0′ be the basic reproductive ratio for BSHSF model. It was calculated using established methods (see Table 6).

Table 5

The result of the basic reproductive ratio from annual report data in the village of Liaonan, Xingzi County, Jiangxi Province, China

The basic reproductive ratios	2003	2004	2005	2006	2007	2008	2009	2010
R 0 (BTHSF model)	1.031	1.074	1.097	1.060	1.057	1.038	1.037	1.030
ℜ 0 (BTH model)	1.052	1.142	1.192	1.114	1.107	1.066	1.064	1.050

Table 6

The result of the basic reproductive ratio from annual report data in the village of Liaonan, Xingzi County, Jiangxi Province, China

The basic reproductive ratios	2003	2004	2005	2006	2007	2008	2009	2010
R 0 (BTHSF model)	1.031	1.074	1.097	1.060	1.057	1.038	1.037	1.030
R 0' (BSHSF model)	1.043	1.035	1.040	1.031	1.031	1.031	1.031	1.032

Barbour [28], proposed the following model without seasonality (BTH model):and obtained the basic reproductive ratio

Using the data in Tables 3 and 4, we obtain the values of R 0 in the village of Liaonan from 2003 to 2010.

R 0 for BTHSF model is shown in Table 3, and it clearly shows that the values range between 1.030 and 1.097 from 2003 to 2010 in the village of Liaonan. Moreover, R 0 for BTH model is always greater than that of BTHSF model at the same time. Furthermore, the variation tendency of R 0 for BTHSF model is shown in Fig. 1.

Fig. 1

Changes in the basic reproductive ratio in Liaonan village from 2003 to 2010

The prevalence of two hosts and vector host has strong effect on R 0 for BTHSF model. Except in 2003 and 2010, the values of R 0 for BTHSF model are greater than those for BSHSF model in the other years (see Table 4). The prevalence of infected bovines in 2003 and 2010 are less than the other years (see Table 2). This phenomenon illustrates that if the prevalence ratio of infected bovine hosts and infected human hosts is small, R 0 for BTHSF model may be smaller than for BSHSF model. Moreover, we note that R 0 and the prevalence of infection in bovines for BTHSF model are higher in 2004 and 2005 among the years under study as shown in (Tables 1 and 3).

Sensitivity analysis

In model (2), infection rates a 1(t) and b 1(t) (i = 1,2) were taken as sine functions of 365 days period, with the form and . Thus, πa /2 and πb /2 reflect the amplitude of infection rates a (t) and b (t), respectively. We fix g1 = 0.00093, g2 = 0.00183, Δ = 91.4602, Σ1 = 0.0475, g1 = 0.00093, Σ2 = 0.0020, μ = 0.0055, a 2 = 0.0426, b 2 = 0.1224, and vary a 1 and b 1 in [0,0.04] in (5), with other parameters unchanged as above, by numerical simulations we obtain the curve of R 0 with respect to a 1 and b 1 (see Fig. 2). It shows that the larger a 1 and b 1 are the higher R 0 is. When a 1 and b 1 are near to 0, R 0 is less than 1, which indicates that the control strategy is very effective and disease transmission will be interrupted. This implies that R 0 is very sensitive to the changes in a 1 and b 1 when their values are higher than 0.03. Next, if we fix a 1 = 0.0068, b 1 = 0.0024, and let a 2 and b 2 vary. A graph indicating the relationship between R 0 and a 2 and b 2 was obtained (Fig. 3). This graph shows that R 0 increases with an increase in the amplitude of a 2 and b 2. Finally, as μ varies in [0,0.012] with other parameters unchanged as above, numerical simulations provide the relationship between R 0 and μ (Fig. 4). Figure 4 shows that R 0 is more sensitive to changes in μ when μ is less than 0.002.

Fig. 2

The relationship between a 1, b 1 and R 0. The graph demonstrates the sensitivity of the basic reproductive ratio to the changes of composite parameters a 1 and b 1

Fig. 3

The relationship between a 2, b 2 and R 0. The graph demonstrates the sensitivity of the basic reproductive ratio to the changes of composite parameters a 2 and b 2

Fig. 4

The relationship between μ and R 0. The graph demonstrates the sensitivity of the basic reproductive ratio to the changes of parameter μ

Discussion

This paper proposed Barbour’s two-host model with seasonality (BTHSF model), which is an improvement on the models BTH and BSHSF. Following the idea of Wang & Zhao [37], the implicit formula of R 0 for BTHSF was given. By computation, the annual changing trend of R 0 s in Liaonan village from 2003 to 2010 was given (see Fig. 1). It showed that R 0 peaked in 2005 at 1.09 and declined dramatically until 2010. There are three possible factors responsible for this decline. First, in 2001, the World Bank Loan Project (WBLP) for schistosomiasis control (1992–2001) was terminated. This led to the resurgence of schistosomiasis transmission after 2001 [38, 39]. Secondly, the implementation of integrated control strategies with emphasis on infection source control commenced in 2005. Finally, the high prevalence of S. japonicum in bovines contributes to the high value of R 0.

By numerical simulation, we realize that the basic reproductive number for Barbour’s model with constant infection rates is always higher than that of periodic infection rates. This illustrates that ignoring seasonal fluctuations would lead to overestimating schistosomiasis transmission risk.

Chemotherapy, health education, provision of water, sanitation and hygiene, bovine control and snail host control are the main measures for schistosomiasis control. We know that the implementation of chemotherapy strategy is beneficial to the recovery and well-being of patients and cattle. This would result in increasing the parameters g 1 and g 2 (recovery rates) in model (2). Accordingly, the implementation of control strategies, such as health education, provision of water, sanitation and hygiene, bovine control would lead to decrease in; the host’s contact rate with contaminated water, the probability that an encounter with cercariae results in the development of an adult parasite and the rate at which schistosome eggs are laid. Therefore, the infection rates (a 1(t), a 2(t), b 1(t), b 2(t)) in model (2) would decrease. More so, the implementation of snail control strategy will increase the removal rate of the snail (μ). Due to the implementation of preventive chemotherapy in Xingzi County, we assessed the sensitivity analysis of the average infection rates and removal rate. Numerical results have shown that R 0 is a monotone increasing function of the four infection rates and a monotone decreasing function of the removal rate. This implies that as the infection rates increases or the removal rate decreases, R 0 increases. Our results also indicate that mollusiciding is an effective means of controlling schistosomiasis transmission in Xingzi County when the removal rate of snails is small.

Finally, in line with the theoretical concept of this paper, we introduce several indexes that were explored extensively. Similar to the methods established in [28], R 0 s are estimated from prevalence data and assumed to be at a steady equilibrium state, so they are always more than 1. In reality, R 0 may be less than 1. However, we can study the transmission intensity by using the values of R 0. For example, if we define the relative basic reproductive ratio as:where denotes the basic reproductive ratio of year k. If R  > 1, the transmission status of schistosomiasis in year j is more serious than year i and vice versa.

Next, we define the relative mean basic reproductive ratio as

Here, shows that the basic reproductive ratio in year m is greater than the average ratio from year i to j, and vice versa. Taking Liaonan village as the example, the relative basic reproductive ratio in 2005 and 2006 is defined as:

This illustrates that the epidemic trend in 2005 is relatively strong, compared to 2006. In addition,and

This indicates that the epidemic trends in 2010 are slower than the average level from 2003 to 2010.

Similarly, it is also able to measure the disease transmission situation between two different regions. Let denote the basic reproductive ratios of two different regions in year i, respectively. The index of the basic reproductive ratio is defined as:

γ can be used to compare the effect of control and intervention between two regions at the same time. However, many factors such as heterogeneity of definitive host, drug resistance and optimal control method are yet to be considered in this mathematical model. These issues need to be addressed in the future studies.

Conclusions

We modified and improved Barbour’s model to investigate the role of seasonal fluctuations in the spread of schistosomiasis. The formula for the basic reproductive number R 0 with a single host was generalized with that of two hosts. Theoretical results showed that it is the basic reproductive ratio which is the threshold value that determines whether schistosomiasis should persist or not. The values of R 0 in Liaonan village from 2003 to 2010 were calculated using the mathematical programming language MATLAB 7.1. The results from this study have shown that ignoring seasonality would overestimate the transmission risk of schistosomiasis, and mollusiciding would be the most effective control measure to curtail schistosomiasis transmission in Xingzi County. Also, the indexes established in this paper will lead to the re-evaluation of control strategies between different areas.