Evaluating spatiotemporal dynamics of snakebite in Sri Lanka: Monthly incidence mapping from a national representative survey sample.

BACKGROUND: Snakebite incidence shows both spatial and temporal variation. However, no study has evaluated spatiotemporal patterns of snakebites across a country or region in detail. We used a nationally representative population sample to evaluate spatiotemporal patterns of snakebite in Sri Lanka.
METHODOLOGY: We conducted a community-based cross-sectional survey representing all nine provinces of Sri Lanka. We interviewed 165 665 people (0.8% of the national population), and snakebite events reported by the respondents were recorded. Sri Lanka is an agricultural country; its central, southern and western parts receive rain mainly from Southwest monsoon (May to September) and northern and eastern parts receive rain mainly from Northeast monsoon (November to February). We developed spatiotemporal models using multivariate Poisson process modelling to explain monthly snakebite and envenoming incidences in the country. These models were developed at the provincial level to explain local spatiotemporal patterns. PRINCIPAL
FINDINGS: Snakebites and envenomings showed clear spatiotemporal patterns. Snakebite hotspots were found in North-Central, North-West, South-West and Eastern Sri Lanka. They exhibited biannual seasonal patterns except in South-Western inlands, which showed triannual seasonality. Envenoming hotspots were confined to North-Central, East and South-West parts of the country. Hotspots in North-Central regions showed triannual seasonal patterns and South-West regions had annual patterns. Hotspots remained persistent throughout the year in Eastern regions. The overall monthly snakebite and envenoming incidences in Sri Lanka were 39 (95%CI: 38-40) and 19 (95%CI: 13-30) per 100 000, respectively, translating into 110 000 (95%CI: 107 500-112 500) snakebites and 45 000 (95%CI: 32 000-73 000) envenomings in a calendar year.
CONCLUSIONS/SIGNIFICANCE: This study provides information on community-based monthly incidence of snakebites and envenomings over the whole country. Thus, it provides useful insights into healthcare decision-making, such as, prioritizing locations to establish specialized centres for snakebite management and allocating resources based on risk assessments which take into account both location and season.

Introduction

Snakebite envenoming is a neglected tropical disease which mainly affects marginalized people living in rural areas of the tropics. Snakebites can result in mild, moderate or severe envenomation, psychological disturbances, long-term disability and deaths. Annually, 4.5–5.4 million people are bitten by snakes globally, with 1.8–2.7 million envenomings, and 81 000–138 000 deaths [1,2]. Snakebite incidence demonstrates both spatial variation and seasonal patterns in many countries [3,4], but there has been no previous in-depth evaluations of spatiotemporal patterns of snakebite risk at a country level.

Sri Lanka is an agricultural tropical island in the Indian Ocean, extending from 5°55′ to 9°51′N and from 79°42′ to 81°53′E. The geography varies from low flat plains (mean elevation 228 meters) to South-Central highland regions (mean elevation 2400 meters). Annual rainfall varies from less than 1000 mm in some areas to over 4500 mm in others, and the country is divided into three climate zones based on rainfall; wet (South-West and Central parts), dry (Northern and South-Eastern parts) and intermediate (the areas in between). Central, southern and western parts receive rain mainly from Southwest monsoon (May to September) and northern and eastern parts receive from Northeast monsoon (November to February) [5,6].

Sri Lanka has a high snakebite incidence [7]. The country has over 100 snake species. Six are medically important; Naja naja (cobra), Bungarus ceylonicus (Ceylon krait), Bungarus caeruleus (Common krait), Daboia russelii (Russell’s viper), Echis carinatus (Saw scaled viper) and Hypnale hypnale (Hump nosed viper). Naja naja and Daboia russelii are widely found throughout the country. Hypnale hypnale can also be found throughout the country, particularly in the wet zone in rubber, tea, coconut and cocoa plantations, and are responsible for the majority of bites in the country. Bungarus caeruleus are found in the dry zone and Bungarus ceylonicus in wet mountains. Saw scaled vipers are confined to the dry arid zones of the Northern and Eastern parts of the country [8-10].

There is a clear geographical variation in the annual snakebite incidence in Sri Lanka [9]. Snakebite incidence rates have also shown seasonal patterns [10]. Therefore, snakebite is likely to demonstrate spatiotemporal variations in Sri Lanka and, to date, no study has undertaken an in-depth evaluation of spatiotemporal patterns of snakebite in the country. Identifying the high-risk geographical areas with their seasonal patterns will help to implement targeted preventive measures and allocation of resources, including antivenom, to manage the snakebite burden in Sri Lanka [11,12]. The aim of this work was to understand the spatiotemporal dynamics of snakebite in Sri Lanka in order to identify both persistent and time-varying snakebite hotspots.

Methods

Ethics statement

The Ethics Review Committee, Faculty of Medicine, University of Kelaniya provided the ethical approval for the study. Permission to conduct the survey in sampled clusters was obtained from the District and Divisional level public administrators and the respective Grama Niladharis were informed before the data collection. Written informed consent was obtained from the participants. No animals were used in the study.

Snakebite is a bite event from any snake, irrespective of envenoming status.

Envenoming bite is a snakebite with significant envenoming, which is defined as the presence of local tissue necrosis at the site of the bite, neurotoxicity or bleeding manifestations

Bite refers to both snakebite and envenoming bite to avoid repetition in the description of the methods, as the same statistical methods were applied to model both snakebite and envenoming bites.

Survey month is the month in which the survey was conducted. The survey was conducted for 11 consecutive survey months, August 2012 to June 2013.

Bitten month is the month in which the victim actually experienced a snakebite. The survey reported bites occurring in 23 bitten months, August 2011 to June 2013.

Recall time is the time difference between the bitten month and survey month of a bite. It ranged from 0 to 12 months.

Data sources

Epidemiological data

Nation-wide snakebite data were obtained from the National Snakebite Survey (NSS). This was a community-based survey designed to sample approximately 1% of the Sri Lankan population. Multistage cluster sampling was used to cover all the 9 provinces and 25 districts of the country. According to the previous national census of population and housing survey in 2001, there were estimated 18.8 million population and 3.9 million households in 18 districts with an average household size of 4.2. We estimated that there will be approximately 20 million people in 4.5 million households in Sri Lanka in 2012 (Department of Census and Statistics, Sri Lanka, 2001). Therefore, the NSS was designed to sample 45 000 households (i.e. 1% of households) in the country, consisting of 5 000 households from each of the nine provinces.

The first stage sampling units were the Grama Niladhari (GN) divisions which were considered as “clusters”for data collection. GN divisions are the smallest administrative units in the country, and Sri Lanka is divided into 14 022 GN divisions. We chose 125 clusters each from the nine provinces for the study. A province is comprised of one or more districts and the 125 clusters were proportionally allocated among the districts within the province. The number of clusters per district was selected proportional to district size. Within each district, the required number of clusters was selected by simple random sampling from the list maintained by the Department of Census and Statistics, Sri Lanka. The second stage sampling unit consisted of individual households. In each cluster, 40 households were selected for the survey, where the initial household was randomly selected based on the electoral register, which was taken as the sample frame. Proximity selection was then used to select subsequent households as the "next nearest" until the desired sample size was reached. In cases of non-response (no one in the household when the interviewers visited) the house was visited once again before the data collection in the cluster was completed. If there was still no respondent, the next nearest house in the cluster was selected.

Data were collected using an interviewer-administered two-part questionnaire. Data collection was facilitated by local volunteers. Data collectors were assisted by local field volunteers recruited within each cluster. The respondent of each household was either the head of the household or a responsible adult present in the house. Initially the research assistant screened the households for snakebite within the previous 12 months and obtained socio-demographic data from the households using the first part of the questionnaire. Subsequently, the second part of the questionnaire were administered to the households where snakebites were reported within the previous 12 months in order to obtain details of the bite (i.e. date and time, place, activity of the victim and whether the offending snake identified or not), evidence of significant envenoming, treatment details and outcome. The questionnaire was pre-tested in each province, fine-tuned prior to use and translated into both official languages of Sri Lanka; Sinhala and Tamil. The survey was conducted for 11 consecutive months from August 2012 to June 2013. We did not consider any specific reason to select the time period for the survey [9].

Explanatory variables

Island-wide cluster level data on population density and percentage of agricultural workers were obtained from the Department of Census and Statistics, Sri Lanka. A rasterized elevation map for the country was obtained from DIVA-GIS [13]. The values of these variables at the centroid of each sampled cluster were considered as cluster-level spatial explanatory variables. Monthly average rainfall, temperature and relative humidity in each bitten month were obtained from the Department of Meteorology, Sri Lanka. The values of these variables at the centroid of sampled clusters in each bitten month were considered as temporal explanatory variables. The median area of a cluster is approximately 2.0 km2 (IQR: 0.9–4.2) and we assumed spatial and temporal explanatory variables do not vary within a cluster.

Statistical methods

Spatiotemporal modelling

We constructed separate spatiotemporal prediction models for snakebites and envenoming bites. Explanatory variables included in previously published separate spatial and temporal models for the NSS data were considered as candidate explanatory variables for the spatiotemporal models. Geostatistical binomial logistic [9] and Poisson log-linear [10] models were used to predict spatial variation in snakebite and envenoming incidences and temporal variation in national incidence, respectively.

Separate spatial and temporal models cannot explain locally varying temporal patterns in the country (i.e. persistent and time-varying snakebite hotspots). Here, we developed a spatiotemporal model for each of Sri Lanka’s nine provinces, allowing for cluster level differences within each province (S1 Appendix). This was done as the survey was designed to estimate bite incidences at the provincial level and models developed at district or cluster level would have been underpowered. The NSS was conducted over 11 months and different clusters were surveyed during different months during the study period. Each cluster captured bite events that occurred over a 13-month period which included the survey month and the preceding 12 months. The NSS data contained the location of each sampled individual which was fixed at the cluster centroid, and the month of each recorded bite which varied over the 13 month period.

Our goal is to describe the rate at which bite events occur over time in each of the clusters. Our working assumption was that, in each cluster, bites occur independently of each other. As statistical model for a process of this kind is a multivariate Poisson process, we modelled the data from each cluster as an inhomogenous Poisson process with cluster-level explanatory variables. We used harmonic mathematical functions (sine and cosine terms) to model the annual, biannual, triannual and quadrennial variations in bite incidence. All the province-level models were adjusted for recall bias. This was done by including recall time as an explanatory variable in the fitted models to estimate recall bias and adjusted for this when constructing our estimates of bite incidence. Finally, the model parameters were estimated by maximising the pooled log-likelihood over each province.

To take account of the uncertainty in estimated incidence maps, we developed PCMs to identify the areas that can be confidently classified as high risk and low-risk, and areas where risk status is highly uncertain. These PCMs quantify the likelihood of exceeding or not exceeding any given threshold level of incidence at a given location and time. Probability values close to one and zero indicate geographical locations with precise classifications, whereas values close to 0.5 indicate locations whose classification is highly uncertain. We selected the monthly mean bite incidences across all the provinces as cut off points to develop PCMs.

All computations used the R programming language version 3.2.3 [14]. Details of the statistical model fitting and prediction of monthly bite incidences are provided in the Technical appendix (S2 Appendix).

Results

Snakebite incidence

Different combinations of spatial variables (elevation from sea level, population density, percentage of agricultural workers in the cluster, average rainfall and temperature) appeared to be important in different provinces. Differences in harmonic terms (i.e. sine and cosine terms) showed different annual, biannual, triannual and quadrennial seasonal snakebite patterns operating in different provinces. Estimates of the recall adjustment factor corresponded to under-reporting of bites in seven of the nine provinces and two provinces where over-reporting was estimated. The two provinces in which the estimated recall factor corresponded to over-reporting had the highest and fourth highest snakebite incidences. The fitted spatiotemporal model is summarized in Table 1 (details of the parameter estimates of the fitted model is provided in S1 Table). The goodness of fit tests (S3 Appendix) and observed versus predicted snakebite incidence for the survey sample (S4 Appendix) did not show evidence of over-dispersion relative to the assumed Poisson process model.

Table 1

Spatiotemporal model for snakebites.

Province	Spatiotemporal model for snakebite incidence
Western	- 1.65 + 0.22*sin(12t) - 0.13*cosine(12t) - 0.11*sin(6t) - 0.20*cosine(6t) + 0.01*sin(4t)+ 0.16*cosine(4t) + 0.37*sin(3t) - 0.21*cosine(3t) - 0.44*rainfall - 0.59*temperature - 0.36*population density + 0.31*agriculture + 0.02*agriculture 9% - 0.01*recall
Central	- 0.41 + 0.29*sin(12t) - 0.12*cosine(12t) - 0.03*sin(6t) - 0.55*cosine(6t) + 0.02*sin(4t)+ 0.08*cosine(4t) + 0.22*sin(3t) - 0.33*cosine(3t) - 1.08*rainfall - 0.39*population density - 0.05*recall
Southern	- 9.18 + 0.06*sin(12t) - 0.26*cosine(12t) - 0.01*sin(6t) - 0.30*cosine(6t) + 0.35*sin(4t)- 0.26*cosine(4t) - 0.12*sin(3t) + 0.21*cosine(3t) - 0.37*rainfall + 0.38*elevation + 0.01*recall
Northern	- 10.30 + 0.19*sin(12t) + 0.07*cosine(12t) - 0.28*sin(6t) - 0.21*cosine(6t) + 0.08*sin(4t) + 0.08*cosine(4t) + 0.03*sin(3t) - 0.36*cosine(3t) + 0.12*agriculture + 0.02*agriculture>9% - 0.06*recall
Eastern	- 1.45 - 0.02*sin(12t) - 0.33*cosine(12t) - 0.12*sin(6t) - 0.13*cosine(6t) + 0.25*sin(4t)- 0.24*cosine(4t) - 0.28*sin(3t) - 0.01*cosine(3t) - 0.91*rainfall - 0.24*population density + 0.64* agriculture>9% - 0.01*recall
North Western	- 3.69 + 0.01*sin(12t) - 0.21*cosine(12t) - 0.02*sin(6t) - 0.14*cosine(6t) - 0.07*sin(4t)- 0.22*cosine(4t) + 0.18*sin(3t) - 0.28*cosine(3t) + 0.67*elevation - 1.37*rainfall - 0.41*agriculture- 0.02*agriculture>9% - 0.05*recall
North Central	+ 0.33 + 0.02*sin(12t) - 0.19*cosine(12t) - 0.15*sin(6t) - 0.21*cosine(6t) + 0.29*sin(4t)+ 0.37*cosine(4t) + 0.29*sin(3t) - 0.01*cosine(3t) - 1.37*rainfall + 0.51*agriculture>9% + 0.01*recall
Uva	- 0.39 - 0.11*sin(12t) - 0.01*cosine(12t) - 0.35*sin(6t) + 0.14*cosine(6t) + 0.02*sin(4t)- 0.34*cosine(4t) - 0.06*sin(3t) - 0.26*cosine(3t) - 1.16*rainfall - 0.01*elevation>195m - 0.35*population density - 0.01*recall
Sabaragamuwa	- 3.37 - 0.20*sin(12t) + 0.04*cosine(12t) - 0.13*sin(6t) - 0.05*cosine(6t) + 0.49*sin(4t)+ 0.02*cosine(4t) + 0.10*sin(3t) - 0.37*cosine(3t) - 0.70*rainfall - 0.23*population density - 0.04*recall

Sin(12t) and cosine(12t), sin(6t) and cosine(6t), sin(4t) and cosine(4t), sin(3t) and cosine(3t) denote the annual, biannual, triannual and quadrennial seasonal patterns respectively; elevation>195m denotes elevation more than 195 meters, agriculture>9% denoted the percentage of agriculture in the cluster is more than 9%.

The spatiotemporal variation in estimated monthly snakebite incidence in the whole country is shown in Fig 1. Snakebite hotspots (i.e. high-risk areas) were mainly scattered in North-Central, North-Western, South-Western and Eastern parts of the country. Most hotspots showed biannual seasonal patterns except in the South-Western inlands where triannual seasonal patterns were noticed. The hotspots were scattered over the country during January. These hotspots then developed along the North-Western and South-Western inlands during February. North-Western hotspots migrated to North-Central regions during April. South-Western hotspots geographically expanded until May. The country as a whole had a relatively low snakebite risk during June and July. In August, hotspots appeared as a band circling the central mountain regions of the country and faded away during the following two months. Thereafter, hotspots appeared in the North-Central inlands in November and then migrated to South-Western inlands during December. Hotspots in the Eastern parts appeared during September and December.

Fig 1

Monthly snakebite incidence in Sri Lanka (bites per 100 000 people).

Average monthly snakebite incidence across all the provinces in Sri Lanka was 39 (95% CI: 38 – 40) bites per 100 000 people per month. This translates into 110 000 (95% CI: 107 500 – 112 500) bites per year in the whole country. The probability of exceeding 39 bites per 100 000 people in each month at any given location is shown in Fig 2. The red islands indicate the most likely places to exceed the average monthly incidence in Sri Lanka. In these areas, the probability of observing more than 39 bites per 100 000 people is at least 0.95. The blue areas have the lowest probability of exceeding the average monthly incidence. In these areas, the probability of observing more than 39 bites per 100 000 people is at most 0.05. Map areas with probability close to 0.5 indicate high uncertainty.

Fig 2

Maps of exceedance probability of more than 39 snakebites per 100 000 people per month.

Envenoming bite incidence

Different provinces showed different combinations of spatial variables (elevation from sea level, population density, percentage of agricultural workers in the cluster, average rainfall and temperature). The harmonic terms (i.e. sine and cosine terms) showed differences in provinces reflecting different annual, biannual, triannual and quadrennial seasonal snakebite patterns. Recall adjustment indicated under-reporting and over-reporting of bites in different provinces. The provinces with over-reporting had either high envenoming incidences or were located adjacent to provinces with high envenoming incidences. The fitted spatiotemporal model is summarized in Table 2 (details of the parameter estimates of the fitted model is provided in S2 Table). The goodness of fit tests (S5 Appendix) and observed versus predicted snakebite incidence for the survey sample (S6 Appendix) did not show evidence of over-dispersion relative to the assumed Poisson process model.

Table 2

Spatiotemporal model for envenoming bites.

Province	Spatiotemporal model for envenoming bite incidence
Western	- 1.36 + 0.27*sin(12t) - 0.51*cosine(12t) - 0.21*sin(6t) - 0.05*cosine(6t) + 0.27*sin(4t)+ 0.34*cosine(4m) - 2.41*temperature - 0.33*population density - 0.19*recall
Central	+ 0.13 + 0.93*sin(12t) - 0.34*cosine(12t) - 0.66*sin(6t) - 0.97*cosine(6t) - 0.30*sin(4t)+ 0.85*cosine(4t) -1.21*rainfall - 0.66*population density + 0.03*recall
Southern	- 3.17 - 0.22*sin(12t) + 0.07*cosine(12t) - 0.28*sin(6t) - 0.38*cosine(6t) + 0.59*sin(4t)- 0.29*cosine(4t) - 2.03*rain + 0.13*elevation>195m + 0.04*agriculture>9% + 0.12*recall
Northern	+ 0.67 + 0.07*sin(12t) + 0.05*cosine(12t) - 0.07*sin(6t) - 0.23*cosine(6t) + 0.02*sin(4t) + 0.16*cosine(4t) - 1.51*rain + 0.29*agriculture - 0.06*recall
Eastern	- 4.13 + 0.21*sin(12t) - 0.12*cosine(12t) - 0.23*sin(6t) + 0.11*cosine(6t) + 0.30*sin(4t)- 0.24*cosine(4t) - 0.51*rainfall - 0.30*population density - 2.35*agriculture + 2.64*agriculture > 9% - 0.01*recall
North Western	-5.09 + 0.26*sin(12t) - 0.66*cosine(12t) + 0.16*sin(6t) - 0.09*cosine(6t) - 0.19*sin(4t)- 0.19*cosine(4t) - 1.02*rainfall + 0.16*population density + 0.01*recall
North Central	-1.92 + 0.13*sin(12t) - 0.13*cosine(12t) - 0.11*sin(6t) - 0.15*cosine(6t) + 0.19*sin(4t)+ 0.44*cosine(4t) - 1.32*rainfall + 0.36*elevation + 0.08*elevation>195m + 0.41*agriculture>9% + 0.05*recall
Uva	-0.14 - 0.12*sin(12t) + 0.16*cosine(12t) - 0.65*sin(6t) + 0.14*cosine(6t) - 0.16*sin(4t)- 0.42*cosine(4t) - 1.49*rain + 0.42*agriculture + 0.01*recall
Sabaragamuwa	-5.67 - 1.36*sin(12t) + 0.06*cosine(12t) - 0.02*sin(6t) + 0.56*cosine(6t) + 0.26*sin(4t)+ 0.08*cosine(4t) - 0.56*rainfall + 0.72*agriculture - 0.05*recall

Sin(12t) and cosine(12t), sin(6t) and cosine(6t), sin(4t) and cosine(4t) denote the annual, biannual and triannual seasonal patterns respectively; elevation>195m denotes elevation more than 195 meters, agriculture>9% denoted percentage of agriculture in the cluster is more than 9%.

The spatiotemporal variation in monthly envenoming incidence in Sri Lanka is shown in Fig 3. Envenoming hotspots were mainly confined to the North-Central, Eastern and South-Western inlands of the country. Hotspots in North-Central parts showed triannual seasonal patterns that geographically expanded during March to May, July to September and November to December. Hotspots in South-Eastern regions appeared in June and then expanded towards South-Eastern regions until November. Eastern parts of the country showed high envenoming risk throughout the year.

Fig 3

Monthly envenoming bite incidence in Sri Lanka (bites per 100 000 people).

Average monthly envenoming across all the provinces in Sri Lanka was 19 (95% CI: 13 – 30) envenomings per 100 000 people. This translates into 45 000 (95% CI: 32 000 – 73 000) envenomings in a year for the whole country. The probability of exceeding 19 bites per 100 000 people (i.e. the average monthly envenoming) in each month at any given location is shown in Fig 4. The red islands indicate the places most likely to exceed the national average monthly envenoming bite incidence. In these areas, the probability of observing more than 19 envenomings per 100 000 people is at least 0.95. The blue areas have the lowest probability to exceed the national average monthly incidence. In these areas the probability of observing more than 19 envenomings per 100 000 people is at most 0.05. Areas with probability close to 0.5 indicates areas of high uncertainty. Here, the North-Western province shows intermediate risk compared to surrounding provinces and appears as a distinct block. This is a by-product of the separate parameter estimation in each province.

Fig 4

Maps of exceedance probability of more than 19 envenoming bites per 100 000 people per month.

Discussion

Our data showed complex spatiotemporal patterns for both overall snakebites and envenomings in Sri Lanka. Snakebite hotspots were mainly observed in North-Central, North-Western, South-Western and Eastern parts of the country. These hotspots showed biannual seasonal patterns except for South-Western inlands, which had triannual patterns. Envenoming bite hotspots were mainly confined to the North-Central, Eastern and South-Western inlands of the country. Hotspots in North-Central regions showed triannual seasonal patterns. Hotspots in South-Eastern regions had annual seasonality. Hotspots in Eastern regions remained throughout the year. Western and Central provinces did not show clear spatiotemporal variation in the pattern of snakebites or envenomings. This could be due to relatively low bite incidences masking such variation. The majority of recorded bites in the South-Western region were non-envenoming. These patterns are compatible with previously reported local [15-17] and regional data [18-22]. The differences between overall snakebites and envenoming bites are likely due to the differences in the distribution of snake species [23-28], their biological behaviour [29-33], rainfall patterns, operating agricultural seasons [6, 34] and human behaviour [35,36]. The results of our work can be used to find solutions to problems such as snakebite using a bioethical perspective. Such an approach would include discussion among the general public, academia and decision makers to identify vulnerabilities of local communities and to reach consensus on a roadmap to minimize snakebites.

We developed our spatiotemporal models at the provincial level. According to our estimates, national average snakebite and envenoming incidences were 39 and 19 per 100 000 per month, respectively. This translates into 110 000 snakebites and 45 000 envenomings in a year. The previous purely spatial analysis estimated the burden in the country to be 80 500 snakebites and 30 500 envenoming bites per year [9], whereas a purely temporal analysis estimated 119 000 snake bites per year [10]. The underestimation from a purely spatial model is likely to be due to the lack of adjustment for seasonal variations and recall bias. The overestimation from a purely temporal model is likely to be due to lack of adjustment for the spatial variation and subnational level seasonal variations. This highlights the importance of spatiotemporal patterns when estimating disease burdens in countries.

We used harmonic mathematical functions (sine and cosine terms) to represent the seasonally varying weather (rainfall, temperature and humidity), agricultural practices and other seasonally varying phenomena in Sri Lanka [10]. In the present analysis, these harmonic functions alone were able to explain the temporal variations in snakebites and envenoming bites. None of the weather indices improved the model fits. However, other studies have indicated associations between snakebite incidence and weather indices such as rainfall [37-39]. A possible explanation for this apparent discrepancy is that harmonic functions in our models act as proxies for seasonal variation in weather-related and any other unmeasured seasonally varying correlates of snakebite incidence.

Our spatial explanatory variables (elevation from sea level, population density, percentage of agricultural workers in a cluster, long-term averaged rainfall and temperature) showed complex patterns of association with snakebite and envenoming incidences in different parts of the country. In general, incidence dropped when the population density increased, as snake habitats are less likely to coexist in population-dense areas [4]. Percentage of agricultural workers living in communities showed positive association with bite incidence, reflecting the close interaction between human and environmental factors [27]. Long-term averages of rainfall and temperature showed inverse associations with bites. Snakes are ectothermic animals and are less likely to inhabit locations subject to extreme weather conditions [29]. However, direct interpretation of the adjusted effects of variables require care. The regression coefficients in these models represent the effect of changing the values of a single variable fixing the values of all other variables. In the current context, our explanatory variables cannot typically be controlled in this way. Hence, our objective is to explain the joint effect of all the explanatory variables rather than their individual effects.

The NSS of Sri Lanka reported snakebite events experienced by family members of households in the preceding year from the date of the interview. It is possible that these reported bite counts are affected by recall bias [40]. We used an exponentially decaying function of recall duration to estimate recall bias in the fitted models and adjusted for this when constructing our estimates of snakebite incidence. The results for different provinces indicated a mix of under-reporting or over-reporting of bites. Under-reporting was estimated in provinces with low bite incidences and over-reporting in provinces with high bite incidences. A possible explanation is that in areas with high incidences, where bites are a common feature of day-to-day life, inhabitants are more likely to report a false positive memory of a bite of a previous month and vice versa [41].

Decision-making requires information of known precision. The PCMs identify the areas that can be confidently classified as high risk and low-risk, and areas where risk status is highly uncertain. Red patches on the PCMs indicate the hotspots and blue patches the cold spots, relative to the monthly mean incidences of snakebites and envenoming bites across all the provinces. These PCMs enable decision-makers to identify priority areas for establishment of specialized snakebite management centres and allocate resources, including snakebite antivenom, at operational level. Further, identification of high-risk geographical areas and seasons of snakebite are important to formulate management plans for sustainable use of natural resources. This can be achieved by environmental education processes and reduction of accidental snakebites will benefit the country including the tourism industry.

In this study, we adopted a point process modelling framework instead of the conventional regression modelling approach in order to make best use of the individual level data provided by the NSS. Times of individual events at given locations were modelled using a multivariate Poisson process. The NSS was conducted for 11 months and the survey sampled different numbers of individuals at each cluster. Our methodology allowed us to estimate snakebites and envenoming bites per person incorporating spatial and temporal variables while accounting for recall bias and varying sample sizes at different clusters. The same methods could be applied for spatiotemporal modelling of other diseases for which the times and places of individual events are recorded.

There are limitations in this study. Although NSS had high spatial coverage capturing snakebite events occurring in all the provinces and districts of Sri Lanka, its temporal coverage was limited to a 13-month sampling window at each sampled location. Consequently, we were only able to disentangle seasonal trends from long-term trends by allowing for possible recall bias. Our spatiotemporal models were built from previously published separate spatial and temporal models, and the explanatory variables considered were limited to the variables used in those models. Envenoming bites are a subset of overall bites. The smaller number and fewer data-points available to fit the point process model resulted in higher standard errors and hence wider confidence intervals for envenomings than for overall snakebites. Our spatiotemporal models were constructed at the provincial level and parameter estimates were optimized for each province as the survey was designed to estimate bite incidences at provincial level. This resulted in artificial discontinuities in our fitted incidence maps at province boundaries, which are particularly apparent along the North-Western province. This can be overcome by designing future surveys at sub-province level or at district level to estimate the burden.

The present study aimed at evaluating spatiotemporal patterns of snakebite and envenoming bite in the country at the time of the survey. However, the climate change is known to affect the weather patterns of the country and hence the distribution of snakes [10]. Therefore, we intend to investigate the spatiotemporal patterns of snakebites in the changing climate, in a future study.

In conclusion, this study is the first demonstration of spatiotemporally interacting variations in snakebites and envenomings in a country. The findings can provide insights into healthcare decision-making at local level, taking into account seasonal variations, to prevent and manage snakebites. Sri Lanka could benefit by establishing a notification system for snakebites. This would provide additional documented information and augment future studies. These methods could also be applied for spatiotemporal modelling of other diseases for which the times and places of individual events are recorded.

Figure that illustrates provinces of Sri Lanka.

(DOCX)

Click here for additional data file.

Technical appendix.

(PDF)

Click here for additional data file.

Figure that illustrates reported snakebites versus predicted snakebites in each province.

(DOCX)

Click here for additional data file.

Figure that illustrates reported and predicted snakebites for the survey sample with probability integral transformation.

(DOCX)

Click here for additional data file.

Figure that illustrates reported envenoming bites versus predicted snakebites in each province.

(DOCX)

Click here for additional data file.

Figure that illustrates reported and predicted envenoming bites for the survey sample with probability integral transformation.

(DOCX)

Click here for additional data file.

Table that illustrates parameter estimates and standard errors in each province level for snakebites.

(DOCX)

Click here for additional data file.

Table that illustrated parameter estimates and standard errors in each province level for envenoming bites.

(DOCX)

Click here for additional data file.

Data.

(CSV)

Click here for additional data file.

4 Jan 2021

Dear Dr Ediriweera,,

Thank you very much for submitting your manuscript "Evaluating spatiotemporal dynamics of snakebite in Sri Lanka: monthly incidence mapping from a national representative survey sample" for consideration at PLOS Neglected Tropical Diseases. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Marilia Sá Carvalho

Associate Editor

PLOS Neglected Tropical Diseases

José Gutiérrez

Deputy Editor

PLOS Neglected Tropical Diseases

***********************

Reviewer's Responses to Questions

Key Review Criteria Required for Acceptance?

As you describe the new analyses required for acceptance, please consider the following:

Methods

-Are the objectives of the study clearly articulated with a clear testable hypothesis stated?

-Is the study design appropriate to address the stated objectives?

-Is the population clearly described and appropriate for the hypothesis being tested?

-Is the sample size sufficient to ensure adequate power to address the hypothesis being tested?

-Were correct statistical analysis used to support conclusions?

-Are there concerns about ethical or regulatory requirements being met?

Reviewer #1: The methodology used was consistent with the results achieved.

It's all ok!

Reviewer #2: -Are the objectives of the study clearly articulated with a clear testable hypothesis stated? Need to improve.

-Is the study design appropriate to address the stated objectives? It has some critical points, but it is acceptable considering the place of study.

-Is the population clearly described and appropriate for the hypothesis being tested? I suggest that he explain better the criteria for choosing the interviewed population.

-Is the sample size sufficient to ensure adequate power to address the hypothesis being tested? Yes.

-Were correct statistical analysis used to support conclusions? Yes.

-Are there concerns about ethical or regulatory requirements being met? No.

Reviewer #3: The study shows objectives articulated with a hypothesis

The study design is appropriate

The population is clearly described

The sample size is god

There is a correct statistical analysis

The study was approved by ethical committees

--------------------

Results

-Does the analysis presented match the analysis plan?

-Are the results clearly and completely presented?

-Are the figures (Tables, Images) of sufficient quality for clarity?

Reviewer #1: It's all ok!

Reviewer #2: -Does the analysis presented match the analysis plan? Yes.

-Are the results clearly and completely presented? It can improve a lot.

-Are the figures (Tables, Images) of sufficient quality for clarity? I found the tables a little confusing, but the images are ok.

Reviewer #3: The analysis presented does match the analysis plan; The results are clearly and completely presented and figures (Tables, Images) are relevant, with quality and clarity

--------------------

Conclusions

-Are the conclusions supported by the data presented?

-Are the limitations of analysis clearly described?

-Do the authors discuss how these data can be helpful to advance our understanding of the topic under study?

-Is public health relevance addressed?

Reviewer #1: It's all ok!

Reviewer #2: -Are the conclusions supported by the data presented? Yes.

-Are the limitations of analysis clearly described? Yes.

-Do the authors discuss how these data can be helpful to advance our understanding of the topic under study? Yes.

-Is public health relevance addressed? Yes.

Reviewer #3: From the perspective of my role in the confluence between animals of medical and bioethical interest, I believe that there could be an innovation in the approach of the issue. I'll point out a few points. First, the authors consider it to be auxiliary in the determination of hotspots to target the service. I believe these strategies could be discussed with the community. In the matter of health control and animal accidents the community needs to see as a problem to help in the resolution. Thus, the existence of committees with popular participation, the academy and the management bodies could take advantage of the data in a more applied and real way. In this context I suggest the insertion of paradigms and perspectives of bioethics. Bioethics identifies vulnerabilities, listens to arguments and proposes solutions taken together and that are consensual and fair for all. the second point I highlight tourism, which in my perspective is an important point for the country. therefore, the knowledge of hotspots could be determinant in management plans for sustainable use of natural areas. consequently, the accidents themselves, the relationship of animal residents and the reduction of accidents could be achieved with environmental education processes. I point out that the work as it is brings an important technical contribution with data relevant to decision making, which in itself are sufficient. However, this suggestion is optional and I believe it could help bring a new perspective to the area.

--------------------

Editorial and Data Presentation Modifications?

Use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity. If the only modifications needed are minor and/or editorial, you may wish to recommend “Minor Revision” or “Accept”.

Reviewer #1: Insert the following reference in the correct way below:

França FGR, Braz VS.. Diversity, activity patterns, and habitat use of the snake fauna of Chapada dos Veadeiros National Park in Central Brazil. Biota Neotrop. 2013; 13:74–85. doi:10.1590/S1676-06032013000100008

Reviewer #2: I believe that the larger images would be more readable.

Reviewer #3: From the perspective of my role in the confluence between animals of medical and bioethical interest, I believe that there could be an innovation in the approach of the issue. I'll point out a few points. First, the authors consider it to be auxiliary in the determination of hotspots to target the service. I believe these strategies could be discussed with the community. In the matter of health control and animal accidents the community needs to see as a problem to help in the resolution. Thus, the existence of committees with popular participation, the academy and the management bodies could take advantage of the data in a more applied and real way. In this context I suggest the insertion of paradigms and perspectives of bioethics. Bioethics identifies vulnerabilities, listens to arguments and proposes solutions taken together and that are consensual and fair for all. the second point I highlight tourism, which in my perspective is an important point for the country. therefore, the knowledge of hotspots could be determinant in management plans for sustainable use of natural areas. consequently, the accidents themselves, the relationship of animal residents and the reduction of accidents could be achieved with environmental education processes. I point out that the work as it is brings an important technical contribution with data relevant to decision making, which in itself are sufficient. However, this suggestion is optional and I believe it could help bring a new perspective to the area.

--------------------

Summary and General Comments

Use this section to provide overall comments, discuss strengths/weaknesses of the study, novelty, significance, general execution and scholarship. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. If requesting major revision, please articulate the new experiments that are needed.

Reviewer #1: It's all ok!

Reviewer #2: Most of my considerations are in the text. I would just like to add, that if possible, make the tables easier to understand.

Reviewer #3: From the perspective of my role in the confluence between animals of medical and bioethical interest, I believe that there could be an innovation in the approach of the issue. I'll point out a few points. First, the authors consider it to be auxiliary in the determination of hotspots to target the service. I believe these strategies could be discussed with the community. In the matter of health control and animal accidents the community needs to see as a problem to help in the resolution. Thus, the existence of committees with popular participation, the academy and the management bodies could take advantage of the data in a more applied and real way. In this context I suggest the insertion of paradigms and perspectives of bioethics. Bioethics identifies vulnerabilities, listens to arguments and proposes solutions taken together and that are consensual and fair for all. the second point I highlight tourism, which in my perspective is an important point for the country. therefore, the knowledge of hotspots could be determinant in management plans for sustainable use of natural areas. consequently, the accidents themselves, the relationship of animal residents and the reduction of accidents could be achieved with environmental education processes. I point out that the work as it is brings an important technical contribution with data relevant to decision making, which in itself are sufficient. However, this suggestion is optional and I believe it could help bring a new perspective to the area.

--------------------

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Paulo Sérgio Bernarde

Reviewer #2: Yes: Gabriela Ferreira Campos Guerra

Reviewer #3: Yes: Marta Luciane Fischer

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see https://journals.plos.org/plosntds/s/submission-guidelines#loc-methods

Submitted filename: PNTD-D-20-01840.docx

Click here for additional data file.

18 Jan 2021

Submitted filename: Response to reviewers.docx

Click here for additional data file.

4 May 2021

Dear Dr Ediriweera:

We are pleased to inform you that your manuscript 'Evaluating spatiotemporal dynamics of snakebite in Sri Lanka: monthly incidence mapping from a national representative survey sample' has been provisionally accepted for publication in PLOS Neglected Tropical Diseases.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Neglected Tropical Diseases.

Best regards,

Marilia Sá Carvalho

Associate Editor

PLOS Neglected Tropical Diseases

José María Gutiérrez

Deputy Editor

PLOS Neglected Tropical Diseases

***********************************************************

Reviewer's Responses to Questions

Key Review Criteria Required for Acceptance?

As you describe the new analyses required for acceptance, please consider the following:

Methods

-Are the objectives of the study clearly articulated with a clear testable hypothesis stated?

-Is the study design appropriate to address the stated objectives?

-Is the population clearly described and appropriate for the hypothesis being tested?

-Is the sample size sufficient to ensure adequate power to address the hypothesis being tested?

-Were correct statistical analysis used to support conclusions?

-Are there concerns about ethical or regulatory requirements being met?

Reviewer #2: (No Response)

Reviewer #3: The authors accepted all suggestions. The manuscript is ready to be published.

**********

Results

-Does the analysis presented match the analysis plan?

-Are the results clearly and completely presented?

-Are the figures (Tables, Images) of sufficient quality for clarity?

Reviewer #2: (No Response)

Reviewer #3: The authors accepted all suggestions. The manuscript is ready to be published.

**********

Conclusions

-Are the conclusions supported by the data presented?

-Are the limitations of analysis clearly described?

-Do the authors discuss how these data can be helpful to advance our understanding of the topic under study?

-Is public health relevance addressed?

Reviewer #2: (No Response)

Reviewer #3: The authors accepted all suggestions. The manuscript is ready to be published.

**********

Editorial and Data Presentation Modifications?

Use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity. If the only modifications needed are minor and/or editorial, you may wish to recommend “Minor Revision” or “Accept”.

Reviewer #2: (No Response)

Reviewer #3: The authors accepted all suggestions. The manuscript is ready to be published.

**********

Summary and General Comments

Use this section to provide overall comments, discuss strengths/weaknesses of the study, novelty, significance, general execution and scholarship. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. If requesting major revision, please articulate the new experiments that are needed.

Reviewer #2: (No Response)

Reviewer #3: The authors accepted all suggestions. The manuscript is ready to be published.

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #2: Yes: Gabriela Ferreira Campos Guerra

Reviewer #3: Yes: ...

Submitted filename: PNTD-D-20-01840_R2.docx

Click here for additional data file.

26 May 2021

Dear Dr Ediriweera,

We are delighted to inform you that your manuscript, "Evaluating spatiotemporal dynamics of snakebite in Sri Lanka: monthly incidence mapping from a national representative survey sample," has been formally accepted for publication in PLOS Neglected Tropical Diseases.

We have now passed your article onto the PLOS Production Department who will complete the rest of the publication process. All authors will receive a confirmation email upon publication.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any scientific or type-setting errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript. Note: Proofs for Front Matter articles (Editorial, Viewpoint, Symposium, Review, etc...) are generated on a different schedule and may not be made available as quickly.

Soon after your final files are uploaded, the early version of your manuscript will be published online unless you opted out of this process. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Neglected Tropical Diseases.

Best regards,

Shaden Kamhawi

co-Editor-in-Chief

PLOS Neglected Tropical Diseases

Paul Brindley

co-Editor-in-Chief

PLOS Neglected Tropical Diseases