Redefining snakebite envenoming as a zoonosis: disease incidence is driven by snake ecology, socioeconomics and anthropogenic impacts

Snakebite is the only WHO-listed, not infectious neglected tropical disease (NTD), although its eco-epidemiology is similar to that of zoonotic infections: envenoming occurs after a vertebrate host contacts a human. Accordingly, snakebite risk represents the interaction between snake and human factors, but their quantification has been limited by data availability. Models of infectious disease transmission are instrumental for the mitigation of NTDs and zoonoses. Here, we represented snake-human interactions with disease transmission models to approximate geospatial estimates of snakebite incidence in Sri Lanka, a global hotspot. Snakebites and envenomings are described by the product of snake and human abundance, mirroring directly transmitted zoonoses. We found that human-snake contact rates vary according to land cover (surrogate of occupation and socioeconomic status), the impacts of humans and climate on snake abundance, and by snake species. Our findings show that redefining snakebite as zoonosis provides a mechanistic eco-epidemiological basis to understand snakebites, and the possible implications of global environmental and demographic change for the burden of snakebite.

118 We first identified a series of mathematical formulations for human-snake contacts to 119 represent snakebite as a zoonotic disease transmission process. In conventional infectious 120 disease transmission models, disease spread is considered to be frequency-, density-121 dependent or a mixture of both. When frequency-dependent, the per-capita rate at which 122 susceptible individuals become infected depends on pathogen prevalence in the population.
123 When density-dependent, transmission increases with infected host density. For snakebites, 124 pathogen prevalence for a given reservoir is 100% (all individuals of a venomous species 125 carry venom; Figure 1), and as such we discarded all frequency-dependent formulations. The 126 remaining density-dependent contact formulations (Table 1) summarise different mixing 127 dynamics between humans and snakes, resulting in functional relationships between snake 128 abundance and snakebite incidence that range from linear to asymptotic to bell-shaped [27].
129 To test the models' abilities to explain snakebites, we transformed them from continuous to 135 Data 136 Data used for fitting and testing models and estimating parameters were two published 137 datasets (rasters) of the spatial distribution of snakebite and envenoming incidence rates, 138 estimated with model-based geostatistics applied to a country-wide community survey of 139~0.8% of the Sri Lankan population [23], under the assumption that these estimates 140 represent the ground truth. The response variables for regressing the functional relationships 141 were the number of snakebite and envenoming cases, respectively, found by multiplying the 142 mapped incidence rates by human population density (described below). The independent . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; 143 data used to explain incidence rates and which represent the causal snakebite factors (  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; 168 values of adjacent cells by a factor of 5 grid cells along longitude and latitude. Snakebite and 169 envenoming incidence layers were resampled from their original ~1.5 × 3 km to the target 170 resolution using weighted bilinear interpolation. The land cover layer was upscaled from 30 171 m to 5 km by majority vote per pixel (land cover class was assigned to each grid cell based on 172 the most common class among the ~27000 30 m grid cells contained in each 5 × 5 km pixel). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. 193 The tested human-snake contact formulations are given in Table 1 in their original and 194 discretised forms. We held no a priori reason to favour one model over another, so we 195 explored all of their relative abilities to explain snakebite and envenoming incidence. 212 Therefore, the probability that there are any snakebites when during one year (t, t +1) is: is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; 214 To summarise, after taking H s from the right to the left hand side of the equations, snakebite 215 or envenoming incidence is proportional to a function of snakes and/or humans, F (H, S ) = P 216 (t, t +1). More generally the number of snake-bitten people is: 218 To select a more appropriate model in the MCMC sampling process we treated snakebite and 219 envenoming cases as either Poisson: 221 or Negative Binomial, parameterised by P N and dispersion parameter r: 225 in order to use the best statistical distribution for the number of cases.
226 Human-snake contact rates 227 To estimate effects of the different snake species and human-related factors, the contact rate β 228 was decomposed into different aspects. The first aspect included the estimation of one contact 229 rate per snake species, such that S t from equation 2.1, for instance is: is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint 236 related to aggressive behaviour (A s ) or envenoming severity (E s ) that influence human-snake 237 interactions and their outcomes.
238 The second component of the decomposed contact rates are a series of functions of human 239 population density (number of people per grid cell) that attempt to adjust total snake 240 abundance in relation to humans, since high human population density is a well known 241 threatening process for biodiversity [33,34]. The functions of human population density were: 243 In this approach, if β 0, 1, 2 coefficients are drawn from a normal distribution in the MCMC 244 sampling process it is possible to model the negative effect of humans on incidence, whilst 245 retaining a positive relationship between the total number of cases and human population 246 density since β (H ) will always be positive. To estimate human susceptibility in the models we 247 categorised their parameters (β and 2.2) in five land cover classes.
248 To summarise, we estimated parameters and measured the ability to represent snakebite and 249 envenoming of the models with and without the functions of human population density (2.2), 250 with and without its parameters categorised by land cover, and estimating a global β with and 251 without categorisation by land cover.

Model of envenoming probability
253 For the approach of treating envenoming as an event that follows a snakebite, we treated the 254 number of envenoming cases as the subset of snakebites that result in envenoming. Therefore, 255 the number of envenoming cases is: is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint 257 where H envenomed is the number of envenoming cases and P env is the probability that a snakebite 258 results in envenoming, which we derived from an expert-collated index of species' 259 envenoming severity ( 264 where, in both equations, B s is the statistical effect of snake species s, E s is the index of 265 envenoming severity (Table 2)  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; 281 of raw number of snakebites and envenomings and their annual incidence rates, by measuring  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint 304 information criterion (DIC; Table 1). The refuge effect model had lower DIC than the simple 305 mass action (Table 1) and higher correlation with the snakebites and incidence data (refuge 306 effect, r = 0.67, d.f. = 128, P = 0; simple mass action, r = 0.61; d.f. = 137; P = 0), suggesting at 307 face value a better fit. We considered that the remaining formulated models were all 308 unsuitable as we failed to obtain reliable parameter estimates due to lack of MCMC 309 convergence.
310 To select one model over another as the better one, we also took into account the statistical 311 distribution of predicted incidence rates. The statistical distribution of the incidence rates 312 produced by the refuge effect model was very different to the original incidence rate 313 distribution (Supplementary materials, Figure S1 and S2). Therefore, we chose the simple 314 mass action model. The number of snakebites and incidence patterns produced by this model 315 were qualitatively very similar and had a statistical distribution nearly identical to that of the 316 data. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint  329 All the parameters of the spatial and non-spatial simple mass action model converged (Table   330 2) and did not exceed the very strict threshold of 1.05 of the Gelman test (very similar 331 variances between and within chains; Supplementary materials, Table S1). The estimated 332 effects of humans on snake abundance resulted in different responses of incidence rates to 333 human population and snake abundance in each land cover class (Figure 3). is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint 334 Table 2. Parameter estimates for the mass action models for snakebites. β i, l parameters are those of the function 335 of human population density to adjust total snake abundance in relation to humans in each land cover class. r l 336 are the negative binomial dispersion parameter for each land cover class. B s are the estimated contact rates for 337 each snake species after adjusting the point intensities for relative abundances and weighting for aggressiveness 338 behaviour. The star * indicates the parameters whose 95% credible intervals do not contain zero (only relevant 339 to β i, l ). Rows with grey background show estimates for the spatial version of the model.  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

Median
The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint 348 Regarding individual species' contact rates, the Indian cobra (Naja naja) had the highest 349 estimated rate, followed by the Ceylon krait (Bungarus ceylonicus) and the common krait (B.  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; 373 The snake species that had the largest significantly positive fitted effect on the probability of 374 envenoming was Bungarus caeruleus, indicating that this species explains most of the spatial 375 heterogeneity of the probability that bites result in envenoming in Sri Lanka (Table 3). D.
376 russellii and N. naja also had significant effects in the outcome of snakebites but were 377 significantly negative, indicating that risk of envenoming after a bite decreased with their 378 predicted abundances (we address this counter-intuitive result in the Discussion). The 379 remaining species' effects were not significantly different from zero, indicating that their 380 contributions towards envenoming cannot be distinguished from random (Table 3).  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint 387 envenoming. Agriculture, had the largest significantly positive effect followed by urban.
388 Degraded forest had the largest significantly negative effect, followed by tea. The only land 389 cover class with a non-significant effect was forest, which suggested that envenoming after a 390 snakebite is more likely to be a function of the biting snake than of the environmental or 391 social context of snake-bitten people in forest environments (Table 3,  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. 457 decrease the probability that bites are envenoming? First question, we infer that a suitable 458 model for field-collected data may be simpler than ours, for which snake abundance estimates 459 in relation to humans should be more robust than currently available [28]. The second point is 460 likely to affect uncertainty of our results, but we interpret the very high similarity between . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; 487 with humans [26], venom toxicity to humans and propensity to inject venom after a bite [20].
488 Lastly, socioeconomic and demographic data may also be used if associated risk factors are 489 well known in the study region. Predictions obtained with the suggested approach will not 490 necessarily represent incidence rates or another measure of burden, but are likely to be 491 broadly correlated with them [28]. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 4, 2021. ; https://doi.org/10.1101/2021.10.01.21264438 doi: medRxiv preprint