A mechanistic model of snakebite as a zoonosis: Envenoming incidence is driven by snake ecology, socioeconomics and its impacts on snakes.

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 modelling 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.

Introduction

Snakebite envenoming causes acute life-threatening disease with long lasting consequences for survivors [1]. By most recent estimates, up to 1.8 million people suffer from snakebite envenoming every year, of which 20,000–94, 000 die of the resulting illness [2]. Such a high burden is increasingly recognised as a global health crisis and, in combination with its relative neglect from a research perspective, has led to snakebite’s recent inclusion on the WHO list of class A neglected tropical diseases (NTDs) [3], a first for a non-infectious disease, and to the development of a global action plan to reduce its burden [4].

Mathematical modelling has been useful for identifying processes that affect the incidence of NTDs and managing them to reduce their burden [5,6]. For instance, modelling revealed that only treating confirmed cases of lymphatic filariasis while neglecting vector and alternative host populations facilitated low-level endemic persistence and the evolution of drug resistant parasites, which together hampered long-term mitigation (reviewed in [7]). Models have also been instrumental for testing and implementing interventions for rabies [8]. The spillover of rabies from its zoonotic mammalian host, mostly as a result of a bite from an infected canine, is a relatively simple transmission process that can be represented with epidemiological models.

Rabies, the above example, and snakebite envenoming have some striking similarities. That is, a pathogenic agent—venom—is transmitted to humans after a contact event with a vertebrate host. Despite this similarity, mathematical models of snakebite are scarce (but see [9]), limiting the extent to which mitigation strategies can be assessed prior to field implementation, for example. In contrast, pharmaceutical solutions still dominate mainstream snakebite mitigation and research agendas (e. g., Snakebites: making treatments safe, effective and accessible | Wellcome; [10]). The snakebite roadmap aims to reduce the number of snakebite deaths and disabilities by 50% by 2030 [11], and identifies novel methods and tools to better understand snakebite epidemiology as a priority to help achieve this goal. Improved epidemiological models for snakebite could fill a major current void in understanding snakebite, improving mitigation efforts and maximising the efficacy of post-bite treatment systems (e.g., directing antivenom supplies efficiently) [5,12].

Mechanistic representation of zoonotic spillover predicates that transmission depends on three major principles [13,14]: 1) reservoir and spillover hosts coinciding in time and space (e.g. [15]), 2) the disease prevalence and/or pathogen shedding in reservoir hosts depending on the contact route necessary (e.g. [16]), and 3) the spillover host developing disease when it is infected. Risk mapping studies of snakebites indicate that there are analogous factors relating to these mechanisms for zoonotic spillover. First, models of snake distribution and abundance, and maps of human population density can be used to represent human-snake spatial [17,18] and temporal [19] alignment, which underlies human-snake contact. Second, the equivalent of pathogen prevalence (i.e,. possessing venom) is 100% among venomous snakes, hence the likelihood of transmission (envenoming) is rather a function of how likely it is that venom is delivered during a bite (as opposed to dry bites), which varies considerably between venomous species [20]. And third, host susceptibility relates to the effect of venoms on the human body [21].

Exposure factors and post exposure vulnerabilities (e.g. occupational, cultural and socioeconomic) may mediate these steps to further influence outcomes and ultimately the overall burden in a population [22-24]. For instance, working and living in a agrarian setting could increase spatio-temporal alignment; poverty could limit the use of protective clothing and also exposure to snakebites due to fragile dwellings made of locally found materials, that can increase the probability of envenoming inside houses given a bite; and cultural differences could influence healthcare seeking behaviour, which could determine the toll on the body during an envenoming event. Fig 1, following Plowright et al.,’s [13] framework for zoonotic spillover, represents such a conceptual arrangement as a sieve that starts with snake diversity, distribution, abundance, behaviour, venom toxicity, and its overlap with humans and individual risk/susceptibility factors.

Fig 1

The snakebite sieve, an adaptation of the conceptual framework for zoonoses of Plowright et al. [13], shows how different factors align and result in what we recognise as snakebite burden.

Species and occupational identities are meant to be schematic and do not represent the characteristics of our study system.

Whereas the majority of NTDs comprise transmission of a pathogenic microorganism (or a complex of microorganism species e.g., Leishmania spp.) from one or more reservoir and/or vector hosts to humans, the transmission dynamics of snakebite involve envenoming of humans (or other victim, such as livestock) from one or more venomous snake (‘reservoir’) species. In systems in which only one species dominates the snakebite burden, an appropriately simple single-species model has been shown to successfully predict the geographical variability of snakebite incidence. For instance, Bravo-Vega et al. [9] modelled the frequency of encounter with Bothrops asper as a function of environmental suitability to predict snakebite incidence in Costa Rica. However, such models have not been applied in settings where multiple species bite and influence variability of snakebite envenoming incidence in space and in time. For example, the ‘big four’ species considerably shape the burden of snakebite envenoming in South Asia [25], while Goldstein et al. [26] show that burden dynamics can be complex over time due to the influence of multiple biting species and the environment on the behaviour of both humans and snakes.

In the present study, we sought to capitalise on recent snakebite research developments in incidence mapping [23], snake distributional ecology [18], and zoonotic spillover theory [13] to develop a novel mechanistic epidemiological framework representing the biological components of snakebite. Specifically, we explored the extent to which various types of models typically applied to directly transmitted infectious diseases explain the geography of snakebite and envenoming incidence estimates. We use the island nation of Sri Lanka, a snakebite hotspot, as a model system given that, as of the time of writing, it is the only high snakebite burden region for which high quality/high resolution country-wide data exist to formulate, test and compare models on snakebite and envenoming incidences. In addition to demonstrating successful risk mapping, we show that geographical patterns of incidence arise from dynamic environmental socio-ecological processes that include effects of climate, occupational risk factors, land use and its changes and direct human impacts on snake populations. Recasting snakebite as a zoonosis and formally applying conventional epidemiological models provides a novel way forward in understanding snakebite epidemiology, which in future studies would allow us to better anticipate risks and potentially help achieve ambitious mitigation targets [11] in a rapidly changing world.

Methods

We first identified a series of mathematical formulations for human-snake contacts to represent snakebite as a zoonotic disease transmission process. In conventional infectious disease transmission models, disease spread is considered to be frequency-, density-dependent or a mixture of both. When frequency-dependent, the per-capita rate at which susceptible individuals become infected depends on pathogen prevalence in the population. When density-dependent, transmission increases with infected host density. For snakebites, pathogen prevalence for a given reservoir is 100% (all individuals of a venomous species carry venom; Fig 1). Frequency dependence might function for snakebites in the presence of dry bites, such as occurs when venom glands are depleted [20]. However, its causes are still poorly understood to be incorporated in a process-based model. The remaining density-dependent contact formulations (Table 1) summarise different mixing dynamics between humans and snakes, resulting in functional relationships between snake abundance and snakebite incidence that range from linear to asymptotic to bell-shaped [27]. To test the models’ abilities to explain snakebites, we transformed them from continuous to discrete time, representing annual trends without seasonality, and estimated their parameters regressing the estimated snakebite and envenoming incidence rates from Ediriweera et al. [23] against each functional relationship. As Ediriweera et al. [23] estimates are spatially explicit we first fitted a model ignoring the spatial component, and then we fitted the same models with a conditional autoregressive random effect.

Table 1

Disease transmission terms tested representing functional relationships between snakes and humans and resulting in snakebite.

Each term has been previously applied in various settings for infectious disease studies so here we list only the earliest reference for each. DIC = deviance information criterion, pD = potential degrees of freedom, H = humans, S = snakes, β = contact rate. Missing DIC values means that we could not obtain a converged model, so DIC values were not comparable with those of converged models.

Tranmission term	Discrete time form F (H, S, θ)	Name and description	Source	DIC, pD*
β HS	1 –exp(-βS)	Simple mass action. H × S is the total number of possible contacts	[30]	17484, 20.02
β H p S q	1 –exp(-β Hp– 1Sq)	Power. Similar to mass action, but p and q take values 0, 1, and increase or decrease the number of contacts of S in relation to H and vice-versa.	[42]	---
β H (S—H/qh)	1—exp(-β (S/H– 1/qh))	Refuge effect on S. Parameter q represents the fraction of snakes exposed to humans depending on the number of humans present	[43]	---
β S (H—S/qs)	1—exp(-β S (1 –S/(qs·H)))	Refuge effect on H. Parameter q represents the fraction of humans exposed to snakes depending on the number of snakes present	[43]	17062, 26.0
βHS·{S1−ε+εSH1−ε+εH	1−exp(−βS·{SH−εH+εHS11−ε+εH)	Separate asymptotic term on H or S, where 0 ≤ ε ≤ 1. When parameter ε = 0, the expression reduces to the simple mass action model, otherwise the effect of H or S becomes asymptotic.	[44]	---
βHS·{1c+S1c+H	1−exp(−βS·{1cH+SH1cH+H^2)	Asymptotic on H or S. Parameter c represents the number of S or H at which 50% of the contacts occur.	[44]	---

Data

Data used for fitting and testing models and estimating parameters were two published datasets (rasters) of the spatial distribution of snakebite and envenoming incidence rates, estimated with model-based geostatistics applied to a country-wide community survey of ~0.8% of the Sri Lankan population [23], under the assumption that these estimates represent the ground truth. The response variables for regressing the functional relationships were the number of snakebite and envenoming cases, respectively, found by multiplying the mapped incidence rates by human population density (described below). The independent data used to explain incidence rates and which represent the causal snakebite factors (Fig 1) were:

Distribution and abundance of reservoir hosts

Raster images of the abundance patterns of the most medically relevant venomous snakes of Sri Lanka (three elapids and four vipers), estimated with point process models as functions of the environment (climate, topography and land cover), and adjusted for species’ relative abundances. Species records to build abundance models span several decades of fieldwork across the island [28].

Distribution and abundance of spillover hosts

Raster image of human population density raster layer from 2010 (closest point in time prior to the snakebite survey period of August 2012 –June 2013; [23]) obtained from the Gridded population of the world (GPW v4) hosted by the Socioeconomic Data and Applications Center (SEDAC, https://sedac.ciesin.columbia.edu).

Spillover host exposure risk factors

Raster image of land cover representing the predominant classes forest, degraded forest, agriculture, urban and tea (see ‘Deriving land cover data’ below for source details). Land cover correlates well with socioeconomic status and predominant occupation [29], which are the primary human-related risk factors [22] and which in the model represent different risk categories via model parameters such as the human-snake contact rate.

Data formatting

Prior to analysis we homogenised and synchronised the resolution of all data to a common grid comprising 5 × 5 km pixels (25 km2) projected to the datum of Sri Lanka (SLD99, EPSG 5235). Synchronising data allowed matching data points to regress the number of snakebite and envenoming cases against human, snake and land cover data according to the model tested. A 5 × 5 km grid was chosen to facilitate computation and retain a reasonable degree of biological detail relevant to the study aims. Human population density and snake abundance estimates were upscaled from their original 1 km resolution by aggregation, summing the values of adjacent cells by a factor of 5 grid cells along longitude and latitude. Snakebite and envenoming incidence layers were resampled from their original ~1.5 × 3 km to the target resolution using weighted bilinear interpolation. The land cover layer was upscaled from 30 m to 5 km by majority vote per pixel (land cover class was assigned to each grid cell based on the most common class among the ~27000 30 m grid cells contained in each 5 × 5 km pixel).

Deriving land cover data

The five categories considered for the analyses were derived using unsupervised isoclustering and visual interpretation of remotely sensed data for the year 2010 (Landsat surface reflectance, original Landsat optical bands and NDVI). The resulting 30 m land cover maps were validated using 600 randomly generated points across Sri Lanka, with which we estimated a classification accuracy of >95%.

Functional relationships for human-snake contacts

The simplest formulation of human-snake contact is the mass action model, whereby the total number of possible different contacts is found by No. Humans × No. Snakes. This formulation is widely used for zoonotic transmission, and the simplest model relevant to snakebites is the susceptible-bitten model [30]. Here, the growth of bitten humans (H) per time unit (dH/dt) is proportional to the number of possible contacts between susceptible humans (H) snakes (S): where β is the human-snake contact rate. This model assumes that time is continuous, however snakebite incidence data has a resolution of one year, for which discrete-time models are more adequate. The continuous (left) and discrete-time (right) models for snakebite based on the SB model are:

The tested human-snake contact formulations are given in Table 1 in their original and discretised forms. We held no a priori reason to favour one model over another, so we explored all of their relative abilities to explain snakebite and envenoming incidence.

Model implementation and selection

Ediriweera et al., [23] report both snakebite and snakebite envenoming incidence nationally for Sri Lanka, which here we refer to as two measures of risk, Snakebite and Envenoming. The former represents all confirmed contacts with snakes with and without envenoming, and the latter are the contacts which resulted in envenoming illness. To model envenoming we used two approaches. First, we modelled snakebite as a contact process with the functional relationships and assumed that envenoming is a subset or secondary event, whose probability of occurring was another function of snake abundance. Second, we treated envenoming cases alone in the same way we treated snakebites (see below).

To estimate model parameters we regressed the number of snakebite or envenoming cases against the functional relationships (Table 1) using Markov Chain Monte Carlo (MCMC) sampling in JAGS [31] via R2Jags package in R [32]. Prior to implementation, the contact processes (Table 1, column 1) were transformed from their continuous-time form into discrete-time, representing the probability that there were any snakebites during the study period (t, t +1). The discrete-time form of Eq 1.4 to represent snakebite incidence is [33]:

Therefore, the probability that there are any snakebites when during one year (t, t +1) is:

To summarise, after taking H from the right to the left hand side of the equations, snakebite or envenoming incidence is proportional to a function of snakes and/or humans, F (H, S) = P (t, t +1). More generally the number of snake-bitten people is:

To select a more appropriate model in the MCMC sampling process we treated snakebite and envenoming cases as either Poisson: or Negative Binomial, parameterised by P and dispersion parameter r: in order to use the best statistical distribution for the number of cases.

Human-snake contact rates

To estimate effects of the different snake species and human-related factors, the contact rate β was decomposed into different aspects. The first aspect included the estimation of one contact rate per snake species, such that S from Eq 2.1, for instance is: the sum of the individual snake species abundances multiplied by their specific human-snake contact rates B, and species aggressiveness or envenoming-severity indices ((A, E); [28]). Each index was included depending on whether the contact process analysed was snakebite (A) or envenoming (A × E). This means that the absolute magnitude of contact rate for snakebite is B × A and for envenoming is B × A × E. Thus B represents other factors not related to aggressive behaviour (A) or envenoming severity (E) that influence human-snake interactions and their outcomes.

The second component of the decomposed contact rates are a series of functions of human population density (number of people per grid cell) that attempt to adjust total snake abundance in relation to humans, since high human population density is a well known threatening process for biodiversity [34,35]. The functions of human population density were:

In this approach, if β coefficients are drawn from a normal distribution in the MCMC sampling process it is possible to model the negative effect of humans on incidence, whilst retaining a positive relationship between the total number of cases and human population density since β (H) will always be positive. To estimate human susceptibility in the models we categorised their parameters (β and 2.2) in five land cover classes.

To summarise, we estimated parameters and measured the ability to represent snakebite and envenoming of the models with and without the functions of human population density (2.2), with and without its parameters categorised by land cover, and estimating a global β with and without categorisation by land cover.

Model of envenoming probability

For the approach of treating envenoming as an event that follows a snakebite, we treated the number of envenoming cases as the subset of snakebites that result in envenoming. Therefore, the number of envenoming cases is: where H is the number of envenoming cases and P is the probability that a snakebite results in envenoming, which we derived from an expert-collated index of species’ envenoming severity (Table 2; [28]). To estimate P we treated the number of envenoming cases as Poisson or Negative binomial, but estimated the probability logistically using two model-formulas: where, in both equations, B is the statistical effect of snake species s, E is the index of envenoming severity (Table 2) and S is snake species s after applying the correction of abundance in relation to humans (β (H)). In model formula 3.2, the term B is a random intercept for land cover class. The criteria to select one formula over the other was minimisation of the deviance information criterion (DIC) and convergence [36].

Table 2

Parameter estimates for the mass action models for snakebites (median and standard deviation).

β parameters are those of the function of human population density to adjust total snake abundance in relation to humans in each land cover class. r are the negative binomial dispersion parameter for each land cover class. B are the estimated contact rates for each snake species after adjusting the point intensities for relative abundances and weighting for aggressiveness behaviour. The star * indicates the parameters whose 95% credible intervals do not contain zero (only relevant to β). Rows with grey background show estimates for the spatial version of the model. H. spp corresponds to Hypnale spp, and T. trig. to T. trigonocephalus. .

θ LC	Agriculture	Degraded	Forest	Tea	Urban
β 0, l	-12.9 (0.19)*	-12.84 (0.2)*	-13.18 (0.21)*	-11.91 (0.26)*	-11.27 (0.61)*
	-12.9 (0.11)*	-12.87 (0.11)*	-13.11 (0.12)*	-11.84 (0.15)*	-11.32 (0.2)*
β 1, l	-0.004 (0.001)*	-0.005 (0.001)*	-0.000 (0.002)	-0.017 (0.002)*	-0.028 (0.002)*
	-0.004 (0.001)*	-0.006 (0.001)*	-0.001 (0.001)	-0.018 (0.001)*	-0.027 (0.002)*
r l	23,70 (2.95)	16.82 (1,13)	9.32 (1,16)	7,18 (0.69)	12,10 (4.52)
	34.47 (2.33)	23.56 (0.91)	13.25 (0.93)	9.39 (0.56)	18.85 (7.02)
θ Spp	B. caeruleus	B. ceylonicus	D. russelli	E. carinatus	H. spp	N. naja	T. trig.
B S	3.48 (0.56)	5.30 (2.85)	1.17 (0.37)	2.99 (1.07)	0.036 (0.02)	8.48 (1.29)	0.17 (0.16)
	4.34 (0.48)	5.58 (2.26)	1.72 (0.39)	2.81 (0.97)	0.047 (0.015)	3.51 (0.69)	0.14 (0.11)

Model selection

We first implemented all functional relationship models and ran short MCMC chains of 10–50 K iterations to discard those with erratic sampling behaviour, poor chain mixing or consistently higher DIC than other functions. Once we obtained a more manageable subset of models we ran the full MCMC chains of up to 750,000 iterations and then checked for convergence with the Gelman diagnostic test [37]. With the latter we made sure that posteriors are unrelated to starting prior values (ratio of between and within chain variances should approach one). Finally, we analysed the spatial pattern and statistical distribution of the residuals by subtracting the median of posterior samples with snakebite and envenoming incidence data.

The criteria to select one functional relationship over another were: 1) ease of parameter estimation and convergence of MCMC chains; 2) adequate reproduction of the spatial pattern of raw number of snakebites and envenomings and their annual incidence rates, by measuring the spatial association of the predictions with a modified T-test for spatially autocorrelated data [38]; 3) minimising the DIC; and 4) adequate representation of the statistical distribution of the snakebite and envenoming data using quantile-quantile plots.

Model with spatial effects

To fit the models with spatial effects we used the mean and standard deviation estimates of the non-spatial version as parameter priors for the main effects β0, β1, and B to ease convergence. Random effects for the above models were incorporated as a log-linear random intercept for the number of predicted cases:

This means that ρ was sampled from a conditional autoregressive normal distribution where each ρ is proportional to the average ρ- of its immediate neighbours in a queen-type neighbourhood. These final models were fitted with Nimble [39] instead of JAGS. Geographical data were manipulated with the raster and rgdal R packages [40,41].

Results

The selected models reproduced well the magnitude and distribution of both snakebite and envenoming incidence rates observed in Ediriweera et al., [23]. National snakebite and envenoming incidence patterns were best predicted by first modelling snakebites with the functional relationship based on simple mass action (Table 1) and then estimating the probability that a snakebite results in envenoming using the land cover random-intercept model (Eq 3.2).

Snakebites model

Of the contact formulations listed in Table 1, the simple mass action (βSH) and refuge effect on humans (βS (H—S/q)) were the best performing models, with the lowest deviance information criterion (DIC; Table 1). The refuge effect model had lower DIC than the simple mass action (Table 1) and higher correlation with the snakebites and incidence data (refuge effect, r = 0.67, d.f. = 128, P = 0; simple mass action, r = 0.61; d.f. = 137; P = 0), suggesting at face value a better fit. We considered that the remaining formulated models were all unsuitable as we failed to obtain reliable parameter estimates due to lack of MCMC convergence.

To select one model over another as the better one, we also took into account the statistical distribution of predicted incidence rates. The statistical distribution of the incidence rates produced by the refuge effect model was very different to the original incidence rate distribution (S1 Fig). Therefore, we chose the simple mass action model. The number of snakebites and incidence patterns produced by this model were qualitatively very similar and had a statistical distribution nearly identical to that of the data (S1 and S2 Figs).

The spatial version of the simple mass action model converged with 1M iterations, and had an even higher correlation with Ediriweera et al. [23]’s estimates, r = 0.87 (d. f. = 81, P = 0) for incidence and r = 0.97 (d. f. = 29, P = 0) for the number of bites (Fig 2).

Fig 2

Snakebite patterns predicted by the spatial model (median of posterior estimates).

Insets in the top right corner of each map show a scatter plot of our model estimates and Ediriweera et al. [23]. Right hand side panels, from top to bottom show the autoregressive random effects, root mean square error for each Sri Lankan district between our model and data used, and the residuals with original survey data analysed by Ediriweera et al. [23].

The decomposition of contact rates that worked best was estimating a contact rate for each snake species and correcting snake abundance in relation to human population density by land cover class and log-transforming human population density:

All the parameters of the spatial and non-spatial simple mass action model converged (Table 2) and did not exceed the very strict threshold of 1.05 of the Gelman test (very similar variances between and within chains; S1 Table). The estimated effects of humans on snake abundance resulted in different responses of incidence rates to human population and snake abundance in each land cover class (Fig 3).

Fig 3

Partial responses of snakebite incidence to human population density and snake abundance per land cover class in the model of snakebites.

The contour lines in each 2D plot indicate the expected incidence level highlighted with the colour scale at that combination of human population and snake abundance.

Model parameters

Parameters of the correction of snake abundance in relation to human population density (β0( and β1(, Eq 4.1) were all significantly negative (Table 2), indicating that humans tend to decrease snake abundance in all land cover classes. Also, given significant differences between land cover classes, the effect of humans on snakes depends on predominant land cover class. The largest effect was estimated for tea and urban cover (snakes decrease faster with human population density), intermediate in forest and degraded forest, and smallest in agriculture (snakes decrease the least with human population density, Table 2).

Regarding individual species’ contact rates, the Indian cobra (Naja naja) had the highest estimated rate, followed by the Ceylon krait (Bungarus ceylonicus) and the common krait (B. caeruleus). The lowest contact rate was estimated for the hump-nosed viper (Hypnale spp). These contact rates showed no relationship with species’ relative abundances, indicating that factors not represented in the decomposed contact rates (i.e., frequency of venom injection, time of activity aligned with human activities) are critical but unobserved determinants of risk patterns. Note, that the estimated rates are not statements of how frequently humans encounter them. Estimated parameters should only be interpreted as indices of species’ importance for snakebites relative to their abundance and that of humans where they occur. Also, these parameters are negatively correlated with those estimated to adjust abundance in relation to humans. Every order of magnitude of increase in the size of the B contact rates must be accompanied by an absolute decrease in log-scale of the value of β parameter of Eq 4.1 (e.g. B × 10 → β0(l)−log(10), or B ÷ 10 → β0(l) + log(10)), allowing the interpretations above for the individual species. The magnitude of these parameters indicates the average number of snakes in logarithmic scale that are lost for every human, resulting in faster declines in urban land cover, and slowest in agriculture (Table 2). The effect of both humans and snakes on snakebite incidence in each land cover class is represented visually in Fig 3.

Envenoming model

The fact that the random-intercept envenoming model (Eq 3.2) was the more adequate, indicated that land cover influenced the probability that a snakebite resulted in envenoming. The correlation between envenoming incidence predicted by the spatial model and the data was r = 0.85 (d. f. = 23, P = 0; Fig 4, right panel) and for the number of envenoming cases was r = 0.93 (d.f. = 9, P = 0; Fig 4, top left panel). Convergence statistics for the spatial model are given in S2 Table, S3 and S4 Figs.

Fig 4

Envenoming patterns predicted by the spatial model of the probability that a bite results in envenoming (median of posterior estimates).

The top right corner of the left side maps show the relationship with the data used to fit the models, and the correlation coefficient adjusted for spatial autocorrelation. Right hand side panels, from top to bottom show the conditional autoregressive random effects, root mean square error for each Sri Lankan district between our model and data used, and the residuals with the original survey data analysed by Ediriweera et al. [23].

The snake species that had the largest significantly positive fitted effect on the probability of envenoming was Bungarus caeruleus, indicating that this species explains most of the spatial heterogeneity of the probability that bites result in envenoming in Sri Lanka (Table 3). D. russellii and N. naja also had significant effects in the outcome of snakebites but were significantly negative, indicating that risk of envenoming after a bite decreased with their predicted abundances (we address this counter-intuitive result in the Discussion). The remaining species’ effects were not significantly different from zero, indicating that their contributions towards envenoming cannot be distinguished from random (Table 3).

Table 3

Parameter estimates of the model for probability of envenoming (median and standard deviation).

Intercept are the estimated effects of each land cover, and B are snake species’ effects. The * symbol indicates that the 95% credible intervals of posterior samples do not contain zero. Rows with grey background show estimates for the spatial version of the model. H. spp corresponds to Hypnale spp, and T. trig. to T. trigonocephalus. .

θ LC	Agriculture	Degraded	Forest	Tea	Urban
Int	0.255 (0.07)*	-0.394 (0.05)*	0.056 (0.07)	-0.72 (0.04)*	0.59 (0.17)*
	-0.05 (0.04)	-0.33 (0.03) *	-0.23 (0.05)*	-0.48 (0.03) *	0.53 (0.11)*
θ Spp	B. caeruleus	B. ceylonicus	D. russelli	E. carinatus	H. spp	N. naja	T.trig.
B S	280.14 (13.99)*	-15.77 (32.12)	-97.34 (28.4)*	5.49 (31.91)	-24.546 (23.74)	-290.34 (29.84)*	-2.88 (31.80)
	197.88 (13.41)*	-3.86 (31.68)	-165.129 (25.61)*	5.22 (31.98)	-34.24 (19.77)	-373.43 (27.66) *	-1.64 (31.59)

Land cover classes also had significant effects on the probability that snakebites resulted in envenoming. Agriculture, had the largest significantly positive effect followed by urban. Degraded forest had the largest significantly negative effect, followed by tea. The only land cover class with a non-significant effect was forest, which suggested that envenoming after a snakebite is more likely to be a function of the biting snake than of the environmental or social context of snake-bitten people in forest environments (Table 3, and see Discussion on the role of land cover).

Discussion

Mathematical models have been critical tools for understanding and controlling the transmission of zoonotic diseases, but despite many ecological and epidemiological similarities few such models exist for snakebite (e. g. [9]). Here, we mapped snakebite and envenoming incidence using a mathematical model that represents human-snake interactions and their outcomes, adapting a mass action model usually applied to the transmission of infectious diseases [30]. We treated venom as pathogenic agent transmitted between venomous snakes and susceptible humans, and tested various functional relationships describing the human-snake contact process and its outcomes. Incidence rates were successfully mapped by estimating contact rates between all medically relevant snakes of Sri Lanka and humans, and by accounting for human and snake factors known to be important determinants of snakebite and envenoming incidence. Human factors (social, economic and cultural) were included categorising parameters by land cover serving as a socio-economic proxy [29]. Snake factors included biological characteristics of the different species, like aggressiveness and severity of the envenoming illness produced by its bite. Furthermore, parameters in the model are sensitive to climate, land cover and topography via their effects on snake population estimates [28]. As such, we have developed a generalisable epidemiological model for snakebite that could be transferrable (given local data) to future environmental conditions.

Predicting spatial and temporal patterns (including future changes) of snakebite risk is possible with purely statistical methods or by focusing only on one or two components of snakebite risk (e.g, [18,45]). An advantage of our approach, however, is that it is process-based and generalisable, explicitly capturing the relevant snakebite processes (e.g., human-snake contact patterns, biological traits of venomous species) to predict epidemiologically meaningful measures of snakebite risk. Such an approach may be more transferrable (i.e., when forecasting burden in other regions) and better suited to forecasting impacts of global change (e.g., including climate change, land-use change and socio-economic development), where a number of complex and potentially interactive mechanisms could push burden in different directions. However, the success of such applications still depends on the availability of relevant data to parameterise the model. The raw material to build our model were snake occurrence records and behavioural trait data [28]. Improving and applying similar frameworks in other regions therefore requires reliable occurrence data in relation to human settlements and improved ecological information on venomous snakes [28,46]., both of which can be regarded as high priorities for future work [47].

To the best of our knowledge, the only previous study of a mathematical model for snakebite comparable to ours, is Bravo-Vega et al. [9]. We considerably extend their approach by decomposing contact rates to incorporate both human and snake factors. The first snake aspect of the decomposed contact rate β represents known (aggressiveness–A−and envenoming–E−indices) and unknown biological aspects (estimated statistically) of multiple species. The unknown estimated snake factors are likely related to the alignment of human-snake activity periods (Table 2 for snakebite model; [26]) and species’ propensities to inject venom during a bite (Table 3 for envenoming model; [20,48]). Also, the estimated parameters summarise snake species’ biologies and how these relate to human social, economic, occupational and cultural aspects relevant for snakebite epidemiology. For instance, rice paddy farmers are more susceptible to Russel’s viper bites (Daboia russelii) because rice is usually harvested barefoot in Sri Lanka [49,50]; common krait (Bungarus caeruleus) bites occur among the poorest of the poor while victims are asleep on the floor [49]; while Hypnale spp. envenoming victims are mostly women who are traditionally in charge of home garden maintenance where this species inhabits leaf litter [51]. These examples show how species’ effects on snakebite are intertwined with human socioeconomic and cultural factors.

The second aspect of the decomposed contact rate was the adjustment of snake abundance as a function of land cover and human population density. Here, land cover may represent predominant occupation and socioeconomic status [29], both of which are known to be important snakebite risk factors [22]. Furthermore, we found that the probability that snakebites result in envenoming (Table 3) is also significantly influenced by land cover, especially in urban and agricultural areas. The latter is supported by empirical evidence, as agricultural workers are at greatest risk of snakebite envenoming [22,52]. Consequently, the estimated effects by land cover class summarise snake responses to humans, and effects of human factors such as agricultural occupation and economic status on snakebite and envenoming incidences. All of these characteristics give our model unparalleled forecasting capabilities to close the gap between public health, ecosystem conservation and sustainability.

In contrast with the conceptual and theoretical strengths of our approach, important questions arise to address in future work, for instance: 1) should models be developed with field data instead of estimates? 2) How does uncertainty and artefacts of incidence estimates affect selection and estimated parameters of our model? 3) Why do highly medically-important species decrease the probability that bites are envenoming according to parameter estimates? For the first question, we infer that a suitable model for field-collected data may be simpler than ours, for which snake abundance estimates in relation to humans should be more robust than currently available [28]. The second point is likely to affect uncertainty of our results, but we interpret the very high similarity between results and estimates as an indication of robustness because most data were independent, apart from human population density, which we discussed above [23]. Finally the negative coefficients estimated for D. russelii and N. naja mean that, in average, the probability that a bite results in envenoming decreases at increasing abundance of these speceis. However the lack of statistical significance of the estimates (credible intervals contain zero, Table 3), reflect a lack of explanatory power of the spatial heterogeneity, not a lack of biological relevance [53]. Therefore, the contradictory estimates from an epidemiological basis result from the higher envenoming incidence in east than west Sri Lanka (Fig 4) because both D. russelli and N. Naja are nearly equally abundant on both sides [28]. An alternative ecological explanation for the estimated negative effects for these two very medically-relevant species could arise if the predicted abundances of N. naja and D. russelli actually serve as proxies for the abundances of non-envenoming species. Abundance surrogacy is relatively common among habitat generalists [54] and non-envenoming snakes are commonly involved in snakebite cases [55], resulting in decreasing risk with the abundance of those non-medically relevant species. In conclusion, given these population-level uncertainties related to the identity of the biting snakes, adequate clinical management in case of a bite should of course be based on expertise and careful clinical evaluation of the envenoming disease.

With careful consideration of the strengths and weaknesses of our approach, we encourage applying and testing our framework in data-poor geographical areas for predicting risk in the absence of other data (e.g., national community survey data) or for testing mitigation interventions. Doing so, however, does require some baseline data as inputs. As mentioned above, the first requirement is a collection of geographical occurrence data of venomous snakes to estimate abundance patterns. Methods and concepts for the analysis of this kind of data in relation to the environment are well established and are described elsewhere (e. g. [56]). Here we used point process models (PPMs) for this purpose given important limitations noted for other common distribution modelling methods (e.g., Maxent) [28,57]. Second is to include some key snake biological/behavioural characteristics for the decomposed contact rates. Relevant traits include aggressiveness, overlap of activity periods with humans [26], venom toxicity to humans and propensity to inject venom after a bite [20]. Lastly, socioeconomic and demographic data may also be used if associated risk factors are well known in the study region. Predictions obtained with the suggested approach will not necessarily represent incidence rates or another measure of burden, but are likely to be broadly correlated with them [28].

Achieving the burden reduction goals laid out in the snakebite roadmap (reducing burden by 50% by 2050; [11]) is an exceptionally ambitious target, requiring advances to our basic understanding of snakebite epidemiology and its treatment. As for other zoonotic diseases, global changes are likely already influencing snakebite and envenoming dynamics, and efficient management will need to accommodate for such changes. However, unlike for many zoonoses, few tools currently exist for snakebite that both shed light on its mechanistic underpinnings and provide avenues for burden mapping and prediction under scenarios of global change or for testing interventions. Studying snakebite as a zoonotic disease has considerably improved our understanding of its epidemiology: it is a dynamic process, its burden is the result of the effects of humans on the abundance of snakes and both affect the burden of snakebite envenoming. The ecological footprint of climate and humans, via land use, represents key characteristics of local populations that are related to snakebite levels. All these factors, and the nature of snake models used, make our model a simple, yet effective tool to forecast the impacts of global environmental change on snakebite burden. Such exercises are necessary to develop relevant interventions for the present and into the future to solve the snakebite crisis.

Quantile-quantile plots of both non-spatial models.

A) Mass-action model and B) refuge effect model. The curved shape of incidence rates for the refuge effect model indicates that the distribution of incidence was very different from the incidence rates used as data.

(TIF)

Click here for additional data file.

Direct comparison of models.

A) Refuge effect and B) Mass action.

(TIF)

Click here for additional data file.

Geweke convergence diagnostics of the spatial random autoregressive effects of the snakebites model.

Estimates are expected to lie within -2 and 2 for sampling convergence.

(TIF)

Click here for additional data file.

Geweke convergence diagnostic of the spatial random autoregressive effects of the envenoming model.

Estimates are expected to lie within -2 and 2 for sampling convergence.

(TIF)

Click here for additional data file.

Upper credible interval of Gelman convergence diagnostic for the spatial mass-action snakebites model.

Upper CI should not exceed 1.05. for sampling convergence.

(XLSX)

Click here for additional data file.

Upper credible interval of Gelman convergence diagnostic for the envenoming spatial model.

Upper CI should not exceed 1.05. for sampling convergence.

(XLSX)

Click here for additional data file.

7 Dec 2021

Dear Dr Martin,

Thank you very much for submitting your manuscript "Redefining snakebite envenoming as a zoonosis: disease incidence is driven by snake ecology, socioeconomics and anthropogenic impacts" 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, including those provided in the attachment.

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. To facilitate this process, we would like to invite you to ensure that all source data is presented in the manuscript and/or electronic supplements, or referenced therein.

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,

Ulrich Kuch

Associate Editor

PLOS Neglected Tropical Diseases

Jean-Philippe Chippaux

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 title of the study "Redefining snakebite envenoming as a zoonosis: disease incidence is driven by snake

ecology, socioeconomics and anthropogenic impacts" is not perfectly matching or at least some sort of misleading to what the method/object/design of this study is.

The manuscript itself deals more with the application of epidemiological models of zoonosis on snakebite data in a certain region, while the title somehow suggests a more fundamental attempt/hypothesis/discussion: whether snakebite is a zoonosis. Regrouping snakebite as a zoonosis would of course broaden the general interest and attention to this devastating neglected disease, what would be highly desirable. Nevertheless, the definition of a zoonosis stated by the WHO is: "A zoonosis is any disease or infection that is naturally transmissible from vertebrate animals to humans" and "A zoonosis is an infectious disease that has jumped from a non-human animal to humans. Zoonotic pathogens may be bacterial, viral or parasitic, or may involve unconventional agents and can spread to humans through direct contact or through food, water or the environment." (https://www.who.int/news-room/fact-sheets/detail/zoonoses)

Snakebite or snake venom neither matches well the medical definition of infection/infectious nor exactly acts as a "transmitted"/"jumping" disease, because it has not been a disease/pathogen to the snake. If snakebite would be suggested as a zoonosis(hypothesis) , this would definitely need more discussion here. Additionally, the sheer applicability of epidemiological models used for zoonosis on snakebites is not a proof of snakebite being a zoonosis itself. These zoonosis models roughly describe the probability of a potentially dangerous vertebrate animal- human encounter, highly differentiated towards socioecological factors etc.. These models can be adapted to lot of non-zoonotic infectious diseases, like (non-vertebrate-animal-) vector borne diseases. If we accept that, like the snake venom, the threat of the animal to the human is not a disease but rather a inherent property/ability of the animal, we can adapt such zoonosis models to a lot of incidents as fatalities by hippopotamus or animal-inflicted car accidents (transmission of physical force/energy), and if we delete the „vertebrate“ we can probably even use a "zoonosis" model for further animal-human encounters as arachnidae/insect stings or lepidopterism.

Reviewer #2: The author used appropriate mathematical formulas to represent zoonotic spillovers for the snakebites in Sri Lanka, which is the hotspot of snakebite envenoming. Several formulas were used to best represent the burden attributed to the association between reservoirs (i.e., snake distribution in their habitat, their nature), their spillover hosts (human density/occupation) and locations. Therefore, the study is appropriately designed to explain the mechanistic eco epidemiological aspects and understanding of snakebite risks relating to demographic data.

To support the hypothesis, the author used previous research to calculate the mathematical formulas for the distribution and abundancy of reservoir hosts (snakes) and spillover hosts (humans). With different methods formulated depending on the condition of snakes and humans, the study is well-articulated to the objectives.

Reviewer #3: Accept

Reviewer #4: Snakebite to be considered a zoonotic disease is a very important study and perspective. The only concern is the data being referred is unclear and more than a decade old. Urbanization, deforestation, increase in population and climate change related untimely season changes have caused environmental changes affecting flora and fauna in many places. Urbanization juxtaposing with farming activities and irresponsible garbage dumping seem to have resulted in a higher number of man-animal conflict cases in countries in Southeast Asia, especially where snake distribution thrives on prey base found abundantly close to human habitations.

The basis of the study is whether snakebite should be considered a zoonotic disease. Perhaps more weightage should have been given to bites, circumstances, envenomed cases and medical case history of envenomed patients.

What the authors intend to justify should be succinctly and clearly stated. Many of the statements in the draft are self contradictory and confusing. The analogy regarding low envenoming by D russelii and Naja naja when abundantly found is very different from what is documented regarding both the species in scientific studies and field observations.

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

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: see above/ I do not want to comment

Reviewer #2: The models of magnitude, distribution of snakebites, and envenoming incidence rates were well explained by the author, and this part has matched the analysis plan. The way the results were presented was clear. The author has clearly indicated how the analysis was done and what kind of models were used in the study.

Reviewer #3: Accept

Reviewer #4: What the authors intend to justify should be succinctly and clearly stated. Many of the statements in the draft are self contradictory and confusing. The analogy regarding low envenoming by D russelii and Naja naja when abundantly found is very different from what is documented regarding both the species in scientific studies and field observations.

Given my background, I was not able to appropriately review the formulae for Functional Relationships for human snake contacts.

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

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: Data-based/Evidence-based medicine is absolutely indispensable for an effective disease control, so providing data, as this study does, is absolutely needed and is especially useful when planning interventions/prevention programs or treating a snakebite without informations regarding the snake. Nevertheless, every snakebite patient still should be assessed/examined/interrogated individually and treated individually. Not every single snake might have read the book when, where an whom to bite, so there might emerge danger when treating cases "prejudiced" by statistics. In my personal opinion, this limit should be mentioned with a few words in the discussion.

Reviewer #2: By describing the formulas, the author explained how mathematical modelling can be fit into epidemiological tools for snakebites. The author considered all the possible factors from both the snake and human sides. The snake sides included different species and their envenoming nature, as well as the aggressiveness of the snake. Human factors include social status, occupation, and culture. With these features incorporated, the conclusion has provided a good sketch of generalised epidemiological models for snakebites.

Reviewer #3: Accept

Reviewer #4: The basis of the study is whether snakebite should be considered a zoonotic disease. Perhaps more weightage should have been given to bites, circumstances, envenomed cases and medical case history of envenomed patients to justify why it should be categorized differently. Instead more effort has gone in snake distribution and contradictory observations of the same species in East and West Sri Lanka.

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

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: (No Response)

Reviewer #2: The study is well-written and I would accept it.

Reviewer #3: Minor revisions:

(Line 106 and 107) Inconsistent use of the terms: “snakebite”, “snake bite”.

(Line 108) Consider replacing line with “In future studies would allow us”.

Typos in lines: 302 (“vaulues”) Figure 3 of supplementary material: (“bundances”)

Format issues: Unequal size of text in Table 2.

Table 3. For clarity, maybe explain abbreviations of “H. spp” - Hypnale spp and “T.trig.” with T. trigonocephalus on the Table 2 description.

Reviewer #4: (No Response)

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

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: see above/ I do not want to comment

Reviewer #2: Many formulas were used in the study, but the author was able to explain them thoroughly. The data presentation and figures are clearly tied to the outcomes

Reviewer #3: GENERAL COMMENTS

• It would be interesting to have more information on how some other ecologic factors had an impact on the model. On lines 324-327, for example the effect of human population on snake abundance is described by land cover classes. Yet it’s important to describe that other species exist in varying degree in such land cover classes and not only humans and that these other species may also have an effect by being predators or prey of the vertebrate host of interest (snakes)

• The authors don’t seem to mention or include in the model the effect of seasonal variability, a well-known factor in snake-human interactions. They seem to address this only somewhat tangentially by considering the alignment of human-snake activity periods. If in fact was part of the analysis, this may warrant some clarification.

• Additional information would be useful to better understand some parameters. Why is “Tea” a land cover class, is it the main crop in the region? Also, it is not very clear how they determined what the most “medically relevant species are” (incidence of accidents, severity of cases, both?)

SPECIFIC COMMENTS / MINOR ISSUES

• (Line 117-118) While it’s quite understandable that they would exclude from the model a frequency-dependent analysis, it could be still applicable as some snakes after biting (not only humans but also prey), can become venom depleted for a few days, which, if following a zoonosis model, might be something to consider.

• (Lines 140-141) The authors mention that the reservoir hosts raster was estimated with data from a 10 to 11-month period. This would be important to clarify as snake activity varies greatly throughout the year and the lack of information of even that month might be important.

• (Lines 194-195) Statistical Analysis Software should be informed in the methodology according to convention. (Seems to be R and its packages (JAGS) and on line 279 (NIMBLE). Also, there is no information for the type of geographical information software used.

Reviewer #4: The way the manuscript is drafted, it needs major revisions. The authors themselves have expressed the need for further studies to justify a few observations.

The basis of the study is whether snakebite should be considered a zoonotic disease. Perhaps more weightage should have been given to bites, circumstances, envenomed cases and medical case history of envenomed patients to justify why it should be categorized as a zoonotic disease.

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

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: No

Reviewer #2: No

Reviewer #3: No

Reviewer #4: No

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, we recommend that you deposit your 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. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

Submitted filename: PNTD-D-21-01445_Reviewer comments.docx

Click here for additional data file.

27 Jan 2022

Submitted filename: Response-reviewers.pdf

Click here for additional data file.

7 Apr 2022

Dear Dr Martin,

We are pleased to inform you that your manuscript 'A mechanistic model of snakebite as a zoonosis: envenoming incidence is driven by snake ecology, socioeconomics and its impacts on snakes' 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,

Ulrich Kuch

Associate Editor

PLOS Neglected Tropical Diseases

Jean-Philippe Chippaux

Deputy Editor

PLOS Neglected Tropical Diseases

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

3 May 2022

Dear Dr Martin,

We are delighted to inform you that your manuscript, "A mechanistic model of snakebite as a zoonosis: envenoming incidence is driven by snake ecology, socioeconomics and its impacts on snakes," 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