Intrinsic and extrinsic drivers of transmission dynamics of hemorrhagic fever with renal syndrome caused by Seoul hantavirus

Seoul hantavirus (SEOV) has recently raised concern by causing geographic range expansion of hemorrhagic fever with renal syndrome (HFRS). SEOV infections in humans are significantly underestimated worldwide and epidemic dynamics of SEOV-related HFRS are poorly understood because of a lack of field data and empirically validated models. Here, we use mathematical models to examine both intrinsic and extrinsic drivers of disease transmission from animal (the Norway rat) to humans in a SEOV-endemic area in China. We found that rat eradication schemes and vaccination campaigns, but below the local elimination threshold, could diminish the amplitude of the HFRS epidemic but did not modify its seasonality. Models demonstrate population dynamics of the rodent host were insensitive to climate variations in urban settings, while relative humidity had a negative effect on the seasonality in transmission. Our study contributes to a better understanding of the epidemiology of SEOV-related HFRS, demonstrates asynchronies between rodent population dynamics and transmission rate, and identifies potential drivers of the SEOV seasonality.

Introduction investigated potential causal relationships among the environmental variables, the population dynamics of the rodent reservoir, and the incidence of HFRS. We subsequently constructed a mathematical model to quantify intrinsic transmission dynamics and extrinsic effects.

Data
Huludao City (40˚56 0 N, 120˚28 0 E) in the Liaoning Province of China (S1 Fig) has a temperate monsoon climate with a hot and rainy summer (mean temperature, 23˚C; mean monthly precipitation, 123 mm) and a cold and rainless winter (mean temperature, -5.7˚C; mean monthly precipitation, 2 mm). The city covers 10,434 square km and has a population of 2.6 million people. HFRS has always posed a severe threat to human health in Huludao City since 1984 [30]. To control HFRS transmission, mass vaccination campaigns against hantavirus were conducted beginning in 2005, and combined with rat extermination programs during 2005 and 2012 [31].
Local demographic data were collected from the Liaoning Statistical Yearbooks. Land cover changes from 1998 to 2015 in Huludao City were derived from the annual European Space Agency (ESA) Climate Change Initiative (CCI) land cover maps, with a 300 m spatial resolution. Meteorological data for Huludao City from 1998 to 2015 were obtained from the Chinese Bureau of Meteorology, and then processed into monthly climate data, including monthly mean maximum temperature, mean temperature, mean minimum temperature, cumulative precipitation, relative humidity, and absolute humidity. Absolute humidity is calculated from temperature and relative humidity [32].
Data of human HFRS cases in Huludao City from 1998 to 2015, were obtained from the Huludao City Center for Disease Control and Prevention, China. All HFRS cases were confirmed according to the standard diagnosis set out by the Ministry of Health of the People's Republic of China, and by detecting antibodies against hantavirus in serum samples. Serum samples were tested for immunoglobulin (Ig) G and IgM antibodies against HTNV and SEOV by indirect immunofluorescent assay (IFA) [33]. Serum-epidemiological surveys on human hantavirus infection were conducted annually between 1998 and 2015. Anonymous (non-personal) information was used.
Rodent population density was investigated indoors and outdoors by the powder-trace method on a monthly basis, obtained from the historical literature [34]. According to the criteria, 400 powdered panels are placed in the special sites for spot checks across the city every month, of which 240 panels are set in restaurant, hotel, station and other public places where rodents are commonly found, and 160 panels are set in general household. Two powdered panels were placed in each room (about 15 m 2 ) in selected spots. Outdoors, each panel was placed along a wall at 5-10 m intervals. All panels were set at night and recovered before sunrise. Rodent population density was calculated as the number of positive powdered panels (with footprints or tail tracks from rodents) divided by the number of effective powdered panels used.

Ethics statement
The study protocol was reviewed by the institutional review board of the Huludao Municipal CDC and ethics approval was not required. We have received consent from home/land owners to collect rodent data on private land and in private homes. The Animal Ethics Committee of the Huludao Municipal CDC also waived approval for this study.

Wavelet analysis
Wavelet analysis is widely used in ecology and epidemiology studies to explore the variety in the periodicity of a time series through its decomposition properties [35]. Here, the Morlet wavelet was used to detect the non-stationary characteristics of the incidence of HFRS fluctuations over time, and the bootstrapping method [36] was performed to test the statistical significance of the results. In the significance test, 1000 surrogate datasets were simulated by bootstrapping to test the null hypothesis, where a P value of < 0.05 was considered to be statistically significant.

Convergent cross mapping
Convergent cross mapping (CCM) is used to detect nonlinear causal relationships between time series of HFRS incidence (Y) and environmental factors (X), which is designed to measure causality in nonlinear dynamical systems. CCM can help to identify bidirectional causality (i.e. X and Y are mutually coupled) or unidirectional causality (e.g., X time series variable in a system has a causal influence on Y, but not vice versa) [37,38]. In this study, CCM was used to identify time-delayed effects of a causal interaction between time series of environmental variability, the population dynamics of rodent hosts, and HFRS incidence based on nonlinear state space reconstruction [39]. In a system where x causes y, Sugihara et al. purposed that the state of x(t) can be estimated from the reconstruction of y(t) using the nearest-neighbor forecasting method, also called cross mapping [40]. Pearson's correlation between the estimated x(t) and observed x(t) reflects the causal effect of x on y, called "cross map skill" [39]. In analysis, the embedding dimension (E) were set according to the simplex projection results [40]. Timedelayed effect (tp) were set as 0-6 months for detecting the time lags [41].
Given the seasonality of HFRS incidence, climate variables, and rodent population density, the seasonal component and the response of HFRS risk and rodent population density to the anomalies of other variables were examined to avoid spurious correlations. For each variable, 1500 surrogate time series having the same degree of shared seasonality, but with randomized anomalies, were generated for seasonal surrogate test. For a specific variable, the month of year average (seasonal cycle) was calculated first, and the seasonal anomaly was represented by the difference between the observed value and the seasonal cycle. Then the random shuffling of the time series data of seasonal anomalies was added back to the season average [42]. The above analysis was implemented in the "rEDM" package in R.

The SEOV-related HFRS transmission model
We estimated the epidemiological parameters for HFRS epidemics in Huludao City by fitting the time series from the observed monthly incidence and rodent population density to a discrete-time susceptible-infection model in the Bayesian framework (Table 1) [43]. Norway rat population dynamic model. First, we modelled the Norway rat population dynamic model based on Verhulst-Pearl logistic growth model [44], without consideration of SEOV infection, and identified the effect of climate on population fluctuation. Population fluctuations are regulated by an identical seasonal birth rate (b), natural mortality (m), environmental carrying capacity (K), and potential extrinsic effect (θ r ) in Eq (1). There is no significant sex difference in the Norway rat population in our study area [45]. To reduce the dimensionality of the model, age classes were omitted.
Here, 1 month is used as the time interval in the time-discrete model; b represents an array consisting of 12 values, representing the reproduction rates of the rodent population over 12 months. As Norway rat breeds throughout the year in the regions studied [45], all of the reproduction rates were set at >0. m denotes the natural mortality of the Norway rat population, whose mean life span is generally not more than 1 year [46]. K represents the environmental carrying capacity and is used to model the environmental pressure from intraspecific competition for food availability or other resources. For identifying the extrinsic effect on Norway rat population fluctuation, we constructed θ r to reflect the potential climate effect. Based on the results of the correlation analyses and CCM, minimum temperature and absolute humidity were potential drivers for rodent population dynamics. We constructed equation θ r = δ 1 TMIN t-4 +δ 2 AH t-3 +ε, where TMIN and AH refer to monthly mean minimum temperature and absolute humidity, respectively, and the parameters δ 1 and δ 2 quantify the degree of climate effect and ε is an intercept.
Integrated SEOV-related HFRS transmission model. The mathematical model for host rodent population outlined in Eq (1) can be extended to incorporate rodent-human SEOV transmission. SEOV infections were modeled as occurring among rodents and from rodents to humans in Eqs (2)-(4). The rodent host population (R N ) was divided into two classes: susceptible (R S ) and infected (R I ). A susceptible rodent could be infected when bitten by an infectious rodent or by inhaling air contaminated with SEOV in the excrement from infectious rodents [47].
In the above equations, all newborn rats were regarded as being susceptible to infection because there is no vertical transmission of SEOV [47]. Given that there is evidence that SEOV infection has no evident impact on death and fertility in Norway rat [48], we did not differentiate mortality (m) and birthrate (b) between R S and R I. β R is the transmission rate between rodents and α allows for non-linear contacts from rodents to humans. H N represents the whole human population and H S represents susceptible human population. Observed cases Y I are assumed to be correlated to the unobserved number of cases according to a binomial observation process with constant observation rate ρ, due to non-reporting from mild symptoms overlooked by the health system and asymptomatic infections. The prior information of ρ (0.006, 0.002-0.009) was assessed from our epidemiological survey. β t represents SEOV transmission, which contains seasonal contact rates between rodents and humans β sea and climate-driven transmission potential β cli , β t = β sea β cli . Based on the results of the correlation analyses and CCM, a causal relationship between relative humidity and HFRS incidence was detected. Therefore, we constructed that β cli = δ 3 RH γ t-1, where RH t-1 was relative humidity with 1-month lag, and δ 3 was used to quantify the effect of relative humidity on HFRS transmission dynamics and γ allowed for the non-linearity of the effect. The results of the estimation for γ and δ 3 were the mainly focused where the direction of the climate effect depends on whether γ (δ 3 >0) is positive or negative. All parameter settings are shown in Table 2.
The Markov chain Monte Carlo (MCMC) sampling method was used for parameter estimation with the MATLAB toolbox delayed rejection adaptive metropolis algorithm [49]. The values of the parameters were estimated by sampling from the convergent posterior distributions of the Markov chains. Gelman-Rubin-Brooks MCMC convergence diagnostic was used to test convergence by identifying whether the series came from a stable distribution or not [50]. Prior distributions of the parameters were mainly set according to scientific literature and values of biological relevance, but those parameters that lacked prior information were set to a large variance (100) with no boundaries ( Table 2). The chains were set to perform 1 million iterations with burn-in of 20,000. Deviance information criterion (DIC) was used to measure the goodness of fit for the models with various framework, and a difference value >10 was regarded to be significant. R 2 was used to measure the interpretability of the model for the observed data.

HFRS epidemics in Huludao City
During 1998 to 2015, 6796 HFRS cases were reported in Huludao City, with an annual incidence from 2.50 (1/100,000, minimum in 2008) to 26.99 (1/100,000, maximum in 2005). The pattern of HFRS incidence showed a significant seasonality with an annual peak in spring from March to May (Fig 1, S2 Fig). The epidemic variability could be divided into three periods: period I (1998-2004), characterized by a high-level of HFRS incidence of over 17.6 (1/ 100,000) and a stable periodicity without limited intervention; period II (2005-2012), vaccination campaign and rat extermination program were implemented simultaneously and the local HFRS incidence decreased dramatically to a relatively low level at 10.0 (1/100,000); period III (2013-2015), the incidence of HFRS gradually rebounded. By using wavelet analysis, we found changes in endemic periodicity after 2006 (Fig 1B), possibly due to the rat extermination program and mass vaccination campaign that were conducted beginning in 2005. Here, to reduce the effect of human interventions on disease dynamics, we focused our analysis on the study period between 1998 and 2004 (period I). In period I, a total of 3914 HFRS cases were reported. The local rodent community mainly comprised Norway rat and the house mouse, with the former accounting for 94.52% of the total. The population density of the rodents ranged from 4.89% to 7.28%, with two relatively low values in January of 1998 and 1999.

Causal relationships between environmental factors, Norway rat population density and HFRS transmission risk
We first examined the impact of land use and land cover changes on HFRS transmission. Land use changed very slightly in period I (S3 Fig), so the effect of habitat and land use change on rodent hosts and disease transmission were minimal during this time period. We then conducted a cross-correlation analysis of the climate variables, rodent population density, and HFRS incidence. The Norway rat population was found to be positively correlated with mean temperature (Pearson's r = 0.23, P < 0.05), maximum temperature (r = 0.25, P < 0.05), minimum temperature (r = 0.23, P < 0.05), cumulative precipitation (r = 0.28, P < 0.05), and absolute humidity (r = 0.24, P < 0.05), with the maximum correlation coefficients occurring at 3, 3, 4, 4, and 4-month time lags, respectively (S4A Fig). HFRS incidence was strongly correlated with the climate variables. All the tested climate variables showed significant negative correlations with HFRS incidence, and maximum cross-correlations between temperature and HFRS incidence occurred at a lag of 3 months and other climate variables at a lag of 2 months  [51,52], the longest time lag was set as 6 months. Convergent cross mapping (CCM) was used to detect causality between these time series. However, rodent data had weak cross mapping skills for the climate variables. Only a marginally significant, weak causal effect was found between minimum temperature (with a 4-month lag, tp = -4) and absolute humidity (with a 3-month lag, tp = -3) on rodent population density, with cross map skills of 0.18 and 0.12, respectively (S5 Fig). To test the non-stationary (transient) correlation between Norway rat and climate variation, we conducted a further test of wavelet coherence analysis. The results showed scattered and small-sized distribution of significant areas with inconsistent phase differences (S6 Fig). In all, the combined results indicate that there might be a weak causal relationship between Norway rat and climate variables. Only relative humidity with a 1-month lag was identified as a significant causal factor for HFRS transmission (cross map skill, 0.86) by the CCM method (

Model validation
To quantify the effect of climatic factors on rodent population density, we tested the models (see Materials and methods for details) with θ r containing (i) minimum temperature, (ii) absolute humidity, (iii) both variables, and (iv) none of these climate variables. The goodness of fit for all the candidate models is provided in Table 1. The model without the climate influence was regarded as optimal to infer the population dynamics of Norway rat in a local urban setting (DIC = 9.17, R 2 = 0.50, Fig 2A). The reproduction rates of the local Norway rat population, as estimated by our model, peaked in June, July, and October, and reached the lowest level in May and August ( Table 2). The average life span of the Norway rat is estimated to be about 6 months, and the local environmental carrying capacity is estimated to be 8.24. The trace plot of Markov chains and the posterior distribution for the parameters in the optimal rodent model is shown in S8 Fig. All  Based on the causal relationship tested by CCM, we constructed a HFRS dynamic model containing the infected rodent population and seasonal contact rate to quantify the specific effect of relative humidity on risk of HFRS (see Materials and methods for details). We constructed two models: one containing relative humidity and the other without. The model with relative humidity had a slightly better fit when the simulated HFRS cases correlated well with the observations and 75% of the variance was explained (R 2 = 0.75, DIC = 5.57) (Fig 2B), relative to the model without relative humidity (R 2 = 0.72, DIC = 7.89). The optimal models also captured the main seasonal patterns of both the rodent population density and reported cases (Fig 2C and 2D). The parameter γ in β cli = δ 3 RH γ t-1, was estimated to be −1.23 (95% CI, -1.62--0.83) (Fig 3A), indicating that relative humidity had a negative impact on HFRS transmission ( Fig 3B). Additionally, during winter-spring, the estimated seasonal contact rates peaked from February to April (Fig 4), which may result from the behaviors of human and rodent. Most of the patients are farmers (77%) and it's the slack time of winter-spring when people stay at home. Besides, due to the severe weather outside during that time, Norway rats would live and feed in closer proximity to human residence for favorable living conditions and available food, which could be a reason for the high contact rate. The trace plot of Markov chains and the posterior distribution for the parameters used in the optimal HFRS transmission model are shown in S9 Fig.

Discussion
We used a mathematical model to quantify the impact of intrinsic and extrinsic factors on risk of HFRS, and to test the mechanistic understanding of HFRS transmission dynamics in a SEOV endemic area of China. By analyzing time-series data of Norway rat populations and SEOV infections, we found that the population dynamics of Norway rat in urban settings were relatively well-predicted by a simple model which excluded climatic conditions. Unusually, we demonstrated that a potential biotic driver (relative humidity) could enhance predictive ability of the HFRS transmission analysis.
Our results highlight the crucial role of mass immunization campaigns and rat extermination program in HFRS transmission control. It should be noted that SEOV infections began to rise again after 2013 in the post-vaccination era. In HFRS endemic areas where hosts are mainly wild rodents, fluctuations in the population density of wild rodents, such as striped field mouse or bank vole, are influenced by climatic factors and can mediate the effect of climate on the risk of HFRS transmission, mainly because of the close association between climate and host food supply [53]. However, in the SEOV-endemic area, the CCM results showed that the change in population density of the Norway rat was insensitive to climate variability. Our results differed from those of a previous study which found a relationship among the virus-carrying index, climate variables, and HFRS incidence by using structural equation modeling (SEM). This study reported that the Norway rat served as an important mediator of disease transmission when climate variability was found to influence the risk of HFRS [54]. This inconsistency may be due to a difference in sampling frequency as we surveyed on a monthly basis, while Guan et al. (2009) sampled quarterly. Additionally, the SEM analysis used a latent variable to account for measurement error but did not evaluate the contribution of specific variables to the transmission dynamics.
Norway rat is a domestic rodent whose lifestyle differs greatly from those of wild rodents. Norway rat lives in residential areas and feeds on various food items from humans instead of field crops [55]. Food resources and habitats are relatively stable and change only slightly with climate variability. Additionally, Norway rats breed throughout the year in Huludao City, where the winter temperatures rarely drop below −10˚C and residential areas have warmer conditions because of the winter heating policy in China. Winter breeding of wild rodents is limited because of low temperatures [56].
Our findings indicated that an abiotic factor, relative humidity, was a critical indicator of HFRS risk in the SEOV-endemic area of Huludao City, which is consistent with previous studies in which SEOV was the main virus type [57,58]. We assumed that relative humidity may be associated with the survival or infectivity of SEOV in the environment, the activity of the rodent, and human hosts or the transmission process, although the underlying mechanism is not clear. We do not know of any relevant experiments on how environmental conditions affect survival of SEOV outside the host, while PUUV and HTNV have been confirmed to have longer infective periods at low temperatures and high humidities in experimental environments [59]. In our model, interannual variation in SEOV infections was partly explained by changes in relative humidity. A less intense outbreak in 2000 might be associated with a lower relative humidity-driven transmission potential (Fig 2B). Some limitations in our study should be mentioned. First, rodent population density was extremely low in the springs of 1998 and 1999, which reduced the overall explanatory ability of our model. Second, prevalence of SEOV infection in Norway rats was assessed by mathematical model due to lack of available data. Third, to ameliorate the influence of vaccination campaign and other interventions, only the selected time range of data was used to detect the relationship between intrinsic/extrinsic drivers and SEOV transmission. Further efforts should be made to improve model reliability through efforts such as addition of data to the available dataset. More attention should be paid to the association between economic development, such as infrastructure improvement, the spatial distribution of Norway rat, and induced SEOV infections across districts in subsequent studies.
The SEOV-related HFRS is posing a health threat to an increasing number of people. Although the Norway rat is found everywhere, SEOV cases are not. We suspect this may mainly be due to increased surveillance and attention focused on this pathogen, as SEOV previously caused many asymptomatic human infections previously. In conclusion, based on the longitudinal and complete dataset, our study yielded a proper framework for understanding the intrinsic transmission dynamics and extrinsic effects on HFRS risk caused by SEOV, especially at the human-animal-environment interface. The framework can be flexibly adjusted when more information is available. Furthermore, we provided clues to potential environmental drivers on SEOV transmission dynamics, which would be useful for further research related to public health issues. for Norway rat population density and climate variables with time lags (tp) and for HFRS incidence (B). The causal relationship was regarded as significant (P < 0.05) when the solid line (cross map skill of the observed data) exceeded the dashed line (surrogate data). TM, monthly mean temperature; TMAX, monthly mean maximum temperature; TMIN, monthly mean minimum temperature; PREC, monthly cumulative precipitation; RH, monthly mean relative humidity; AH, monthly mean absolutely humidity. (TIF) S6 Fig. Wavelet coherence between Norway rat population density and climate variables from 1998 to 2004. Arrows showed phase differences between Norway rat population density and climate variables. The straight down represents climate variables lead Norway rat population density 1/4 period. Legend is wavelet coherence level. The parts inside white contour shows a significant relationship. The parts outside the transparent cone have affected by edge effect. TM, monthly mean temperature; TMAX, monthly mean maximum temperature; TMIN, monthly mean minimum temperature; PREC, monthly cumulative precipitation; RH, monthly mean relative humidity; AH, monthly mean absolutely humidity; RN, Norway rat population density.