https://orcid.org/0000-0002-5588-7396
Figures
Abstract
The role of climate driving zoonotic diseases’ population dynamics has typically been addressed via retrospective analyses of national aggregated incidence records. A central question in epidemiology has been whether seasonal and interannual cycles are driven by climate variation or generated by socioeconomic factors. Here, we use compartmental models to quantify the role of rainfall and temperature in the dynamics of snakebite, which is one of the primary neglected tropical diseases. We took advantage of space-time datasets of snakebite incidence, rainfall, and temperature for Colombia and combined it with stochastic compartmental models and iterated filtering methods to show the role of rainfall-driven seasonality modulating the encounter frequency with venomous snakes. Then we identified six zones with different rainfall patterns to demonstrate that the relationship between rainfall and snakebite incidence was heterogeneous in space. We show that rainfall only drives snakebite incidence in regions with marked dry seasons, where rainfall becomes the limiting resource, while temperature does not modulate snakebite incidence. In addition, the encounter frequency differs between regions, and it is higher in regions where Bothrops atrox can be found. Our results show how the heterogeneous spatial distribution of snakebite risk seasonality in the country may be related to important traits of venomous snakes’ natural history.
Author summary
Snakebite envenoming is a neglected tropical disease characterized by its high burden on the rural population and high mortality if antivenom is not administered. The ecology of this health problem is not well-understood; however, approaches to address the temporal are growing. So far, we know that rainfall can play an important role in driving snakebite incidence seasonality at a national scale. Moreover, geographical areas with high rainfall are more prone to have high snakebite risk, but the spatial heterogeneity of the temporal association (i.e., if there are different seasonal patterns of rainfall-incidence association in different geographical areas of a country) is just starting to emerge in the literature. By formulating and fitting compartmental models to data, we generated a flexible framework that relies on temporal resolved datasets and a compartmental mathematical model to understand the effect of climatic covariates (such as rainfall and temperature) driving snakebite dynamics in space and time. We applied this framework to Colombia and found that dry seasons cause a decrease in snakebite incidence: Rainfall only drives snakebite dynamics in regions with marked dry seasons. Thus, rainfall is a limiting resource of the system, and its effect is not spatially homogeneous. On the other hand, the temperature had no significant effect driving snakebite incidence. Our modeling approach can also be used to estimate the effect of climate anomalies on snakebite incidence and has the potential to be used as a tool to monitor snakebite incidence.
Citation: Bravo-Vega C, Santos-Vega M, Cordovez JM (2022) Disentangling snakebite dynamics in Colombia: How does rainfall and temperature drive snakebite temporal patterns? PLoS Negl Trop Dis 16(3): e0010270. https://doi.org/10.1371/journal.pntd.0010270
Editor: Stuart Robert Ainsworth, Liverpool School of Tropical Medicine, UNITED KINGDOM
Received: September 1, 2021; Accepted: February 21, 2022; Published: March 31, 2022
Copyright: © 2022 Bravo-Vega et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Data and codes to replicate the results of this study are going to be available in figshare server: DOI: 10.6084/m9.figshare.15029700, (https://figshare.com/s/305c1d006d433debe14d).
Funding: This study was partially supported by: Minciencias, Colombia (https://www.minciencias.gov.co/), Application 727 for doctoral student to CBV, and Universidad de los Andes, Colombia (https://uniandes.edu.co/), Funding program for doctoral students, awarded to CBV. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Snakebite envenoming is the neglected tropical disease (NTD) with the highest mortality rate, but despite its importance, data collection is challenging, and the actual disease burden is underestimated [1,2]. The most reliable estimates are at least 1.8 million cases with 435.000 deaths each year. However, this data can underestimate the real burden because patients prefer seeking traditional medicine instead of allopathic medicine, and mortality data do not quantify the morbidity associated with physical and mental disabilities caused by some snakebites [3–9]. Estimating the real burden is crucial to understand the ecological characteristics that modulate human-snake encounters, where the climate can play an important role in driving the ecology of venomous snakes [10–13]. Although several studies have estimated snakebite risk’s spatial and temporal heterogeneity based on statistical analyses and climatic variables, the nature of these models is empiric: Their main findings are important correlations between snakebite risk and covariates [7,14–17]. Thus, their usage as disease monitoring tools is limited due to their correlational nature, which mostly hides reality’s complexity [1,18,19].
Compartmental modeling has emerged as a tool to understand the processes underlying disease transmission, and it has been used widely for several NTDs [20–25]. This modeling approach subdivides population into a set of compartments, where flows represent movement between the compartments. Usually, these flows are important parameters, e.g., incidence, mortality, and recovery rate. Thus, covariates can affect specific flows in the model, allowing us to understand the relationship between covariates and disease spread and making the model flexible and an excellent tool for generating fundamental knowledge about the basic mechanisms behind disease spread [26]. Contrarily, empirical and correlational approaches use general equations and parameters not based on disease dynamics, contributing with a more general knowledge about the general behavior of disease spread [20,26]. Even so, compartmental models rely strongly on disease’ ecology knowledge, and as in empirical models, the quality of the data. For example, although snakebite envenoming is not an infectious disease, it is caused by the interaction between humans and venomous snakes, and that interaction can be modeled similarly to an infectious disease. A previous study demonstrated that the main assumption of several compartmental models (the law of mass action) could explain snakebite geographical variation in Costa Rica [27]. By combining compartmental models with spatial and temporal surveillance data, it is possible to have reliable estimations for multiple epidemiological parameters and to perform fine political scale extrapolations of the model on other countries with different covariates, making these models a valuable tool to understand, manage and control several NTDs [28,29].
To apply compartmental modeling to snakebites, it is essential to understand the most important venomous snake species ’ biology, which is often poor [11,30]. In the Neotropics, one of the most important venomous snake groups is the genus Bothrops (Wagler, 1824), which is distributed broadly from southern Mexico to northern Argentina and causes the majority of envenomings in this area [31–33]. This genus is viviparous with an average clutch size between 3.5–37 offsprings, mainly during the rainy season [10,12,34–36]. This seasonal dynamic usually increases the venomous snake population during the rainy season, thus increasing snakebite risk during this period. Additionally, rainy seasons can cause flooding, habitat perturbation, and increased prey abundance, making venomous snakes more active, resulting in an additional increase in snakebite risk [37–39]. Previous studies have found that snakebite temporal dynamics tend to be related to seasonal rainfall patterns in tropical countries such as Costa Rica and Sri Lanka, but these studies did not account for climatic heterogeneity [15,16]. Later, in Sri Lanka a study used multivariate Poisson process modeling to evaluate the spatio-temporal association between incidence and climatic and socioeconomic variables, thus finding the spatio-temporal apparition of hotspots in that country [40]. Therefore, the spatial heterogeneity of the temporal association between climate and snakebite incidence has been addressed by statistical modeling rather than more robust modeling approaches such as compartmental modeling, and it has not been studied in neotropical countries.
Colombia is an ideal setting to study this spatial heterogeneity given its location with diverse climatic conditions, where different regions have different rainfall and temperature patterns [41]. In the country, snakebite is a severe public health problem, where recent reports account for around 4500 envenoming cases with at least 40 deaths each year, and two species (Bothrops asper and Bothrops atrox) cause most envenoming’s [42,43]. In this study, we developed and calibrated a compartmental model that can disentangle the role of rainfall and temperature on snakebite temporal patterns. We used compartmental stochastic models combined with iterated filtering statistical inference methods to explore the role of these climatic drivers on snakebite temporal patterns and to understand this association’s geographic distribution.
Materials and methods
We followed four steps to develop and fit a compartmental modeling approach that accounts for the role of temperature and rainfall on snakebite dynamics. i) We developed two compartmental models for national data to determine the effect of temperature and rainfall driving national snakebite incidence. ii) We divided the country into six regions with similar rainfall patterns: This had the aim to aggregate municipality data within these regions to reduce noise caused by the low occurrence of snakebites and heterogeneous climate patterns. iii) We proposed a modeling scheme with different compartmental models to test hypotheses for the association between rainfall and snakebite incidence for the country’s six regions. iv) We fitted our modeling scheme in each region using an iterated filtering approach that maximizes likelihood [44].
Climate and incidence data
Incidence data.
We got publicly available snakebite incidence data between January of 2010 to the end of October of 2016 from the Sistema Nacional de Vigilancia en Salud Pública (SIVIGILA) of Colombia. This dataset is reported in epidemiological weeks, and it details the number of cases for each municipality (the smallest political unit in the country). We only worked with municipalities that reported snakebite for all the study years. First, we looked at the epidemiological calendar of Colombia to convert the timescale of the reported cases from epidemiological weeks to months. We aggregated the weekly reported cases for each month. Then, for weeks that overlapped between two months, we distributed the cases by weighting them based on the number of days of the week that belong to each month. Finally, we removed the temporal trend of the timeseries using a locally estimated scatterplot smoothing (‘stats’ package, R environment v. 4.1.1. [45,46]) to only account for seasonality and inter-annual variation in the data.
Rainfall and temperature data.
We got monthly raster maps for minimum and maximum temperature and rainfall for Colombia between the year 2010 to the end of October 2016 with a resolution of ~ 21.23 km2 from TerraClimate dataset (http://www.climatologylab.org/terraclimate.html) [47]. In addition, we removed climate data from areas where the two most important venomous snakes’ species in the country are absent by using an altitude threshold of 1900 m.a.s.l. for Bothrops asper and of 1500 m.a.s.l. for Bothrops atrox, which are the species that cause most envenomings in the country [33].
Clustering of regions with similar rainfall pattern.
We used a clustering algorithm over rainfall maps to define regions with similar rainfall patterns. Given that we had monthly maps in raster format, each pixel has a rainfall value in a specific month in the area in this pixel; thus, each pixel contains a time series. Before performing the clustering algorithm, we first reduced the resolution of rainfall maps by a factor of 12, where the resolution of the new monthly rainfall maps was 12 times the previous resolution (~ 254 km2). Then, the clustering algorithm groups pixels based on the similarity of the shape of their rainfall time series, where each resulting group of pixels will represent a region with a similar rainfall pattern. Finally, we used a k-shape algorithm to cluster pixels by using a distance matrix constructed with a shape-based distance [48]. These clustering algorithms use the required number of groups as a parameter (i.e., the number of final regions N, so the algorithm will classify each pixel in each one of these N regions), so it is necessary to determine the optimal number of clusters. We evaluated clustering performance between N = 2 and N = 10 regions by using the silhouette (SIL) and Davies-Bouldin star (DB star) indices. These indexes evaluate the capacity of the clustering algorithm to group data by compensating the number of clusters. Thus, we selected the optimal number of clusters by searching the minimum DB star and maximum SIL indexes [49]. The result of this algorithm is a map that defines regions with similar rainfall patterns, and this algorithm was made by using the package dwtclust in R environment v. 4.1.1. [50].
Generating incidence and precipitation data for each region.
We calculated the average rainfall in each region by computing the mean of rainfall maps’ pixels in each region. Then, to compute average incidence in each region, we first rasterized municipality-scale incidence and normalized it by the number of pixels contained in each municipality. This normalization generates raster maps of monthly incidence per pixel. Finally, we summed the incidence per pixel for all pixels contained in each region. Given that we normalized the rasterized incidence by the number of pixels per each municipality, incidence per region will not be inflated after summing incidence per pixel (i.e., National incidence by summing incidences for all municipalities and by summing incidence over all regions is the same). Nevertheless, this incidence distribution has the important assumption of spatial homogeneity within each municipality. Therefore, we did the same process for the total population for each municipality, where we got the total population from the national statistics department DANE. The outputs of this process are: i) Monthly precipitation for each region from 2010 to October of 2016, ii) Monthly reported incidence for each region for the same timespan, iii) Total average population for the same timespan for each region.
Compartmental models
We formulated our compartmental models by assuming that incidence is proportional to the susceptible population (S) and that the parameter representing this proportionality is β: The encounter frequency with venomous snakes. We also implemented a normal random noise to this frequency, the parameter eps. Additionally, we assumed that each time step a proportion γ of the envenomed population (E) recovers from the envenomation (Fig 1A, see S1 Text for a detailed explanation). The first model (Model 1) is our null hypothesis where climate covariates do not modulate this encounter frequency, but the second model does. We first used both models for national data by including as covariates for Model 2 maximum and minimum temperature and rainfall in independent models. In these second models (Model 2R: Rainfall, Model 2MXT: Maximum temperature, Model 2MNT: Minimum temperature), we established the relation between β and climate covariates as a type III functional response (View Fig 1B).
A. General compartmental model. Our general model assumes two populations: Susceptible (S) and Envenomed (E). We defined the flow from S to E as the new envenoming cases per unit of time, which is the incidence. This incidence depends on a gaussian noise (eps), the encounter frequency parameter (β, view [27] for details), and S. By assuming a constant population, we defined a discrete model that only depends on E. B. Compartmental modeling scheme. We defined this scheme to test the effect of climatic covariates on snakebite incidence. First, we tested if the climatic covariate drives snakebite incidence by comparing between a Model where contact rate is constant (Model 1) with a model where this contact rate is modulated by a type III functional response (Model 2). In Model 2, we tested three climatic covariates for national data: Maximum temperature, Minimum temperature and Rainfall. Given that nationally only rainfall modulates snakebite incidence, we did the regional modeling scheme by only using rainfall variables. Then, in regions where Model 1 adjusted better data than Model 2 (Areas where incidence is not modulated by climatic covariates), we tested with Model 3 for rainfall-independent seasonality in incidence timeseries. Finally, in regions where Model 2R (Model 2 with rainfall as covariate) adjusted better data, we tested if the seasonal component of this rainfall (Rainfall-driven seasonality) is the driver behind this association by using Model 4.
Afterwards, we fitted the models to our defined regions to test how rainfall modulates snakebite incidence. We did not use Model 2MXT and Model 2MNT in clustered regions because we did not find any relationship between both variables and incidence in national data (View Fig 2). To test if incidence has rainfall-independent seasonality in regions where rainfall does not drive incidence (Model 1 adjusted better data than Model 2R), we proposed a third model (Model 3), which uses 4 B-splines: Here β is a linear combination of the four B-splines. These splines are functions that make seasonal peaks on different months of the year (View Fig A in S1 Text). Finally, Model 4 tests if the seasonal component of rainfall (Rainfall-driven seasonality) is the mechanism that modulates snakebite incidence in regions where Model 2R fitted data better than Model 1. To extract the seasonal component of rainfall, we used the algorithm based on local regression described in [46]. The four models that compose our modeling scheme can be seen in Fig 1B. (View S1 Text for details).
We show the results of the best model (Model 4) and the models with temperature covariates (Model 2MXT and Model 2MNT) for national data. Climatic data shows the climatic covariates used per each model, and the simulation after parameter adjustment represents the model with the best combination of parameters that capture incidence dynamics. Note that the best model (Model 4) had lower AIC than the other models and note the difference between rainfall-driven seasonality (Dashed light-blue line) and rainfall signal (Dark blue line). In the simulations, Model 4 captures better incidence seasonality than Model 2MXT and Model 2MNT, where the median of the simulations (Dashed red line) is relatively constant.
Fitting compartmental models to snakebite incidence data.
We assumed the process of snakebite envenoming described by compartmental models as a partially-observed Markov process (POMP), where we can declare an observation model to estimate the error of data-reporting [44]. This observation model was a Poisson process, and the mean of this Poisson model is the monthly flow of people between S (Susceptible population) and E (Envenomed population). This flow is the reported incidence, and it is modeled as the expression shown in the arrow between S and E in Fig 1A.
Before estimating the parameters with the maximum likelihood, we detrended rainfall as we did with incidence [46], and we normalized detrended rainfall between 0 and 1. Then, we estimated the parameters that maximize the likelihood between our model and the data by iterated particle filtering, where we defined a parameter space between the biological limits of the values of the parameters of the models. We defined a random walk for each point in the parameter space, so the algorithm computes the likelihood for each combination of parameters. Finally, the algorithm converges when the likelihood reaches a “global” maximum (See S2 Text for details) [44]. To fit models, we used the ‘pomp’ package in R v. 4.1.1. (40).
Best model selection.
After computing the parameter combinations that maximize likelihood per model per region (View next section), we compared models’ performance using the Akaike information criterion (AIC) to weight model capacity to represent data by its number of parameters. Thanks to this criterion, we were able to choose the best model for each region and evade overfitting caused by redundant model complexity, and we used a between-model threshold of 2 AIC units to determine significant differences [51,52]. We first selected the best model between Model 1 and Model 2R for every region, so if Model 2R outperforms Model 1, then rainfall explains the dynamics of snakebite. Then, for regions where Model 2R outperformed Model 1, we used Model 4 to test if rainfall-driven seasonality is the mechanism behind the association in both. Finally, for regions where Model 1 outperformed Model 2R (no association between rainfall and snakebite incidence), we used Model 3 to test any rainfall-independent seasonality in data.
Results
In Colombia, between 201 and 422 envenomings are reported every month, with an average of 3659 cases per year (Min: 3135 cases in 2010, Max: 4089 cases in 2015). We found that the number of reported cases has increased with time (Pearson correlation coefficient: 0.58, p-value < 0.05), a result that the improving reporting system and population growth might cause rather than an increase in the actual incidence. After model selection for national data, Model 2MXT and Model 2MNT had a greater AIC value than the Model 1 (AIC for Model 1: 812.22, Model 2MXT: 819.31, Model 2MNT: 816.38); hence the temperature does not explain incidence variation in the country. On the other hand, the best model for rainfall was Model 4 (AIC: 793.323), where rainfall-driven seasonality is the driver that modulates snakebite incidence dynamics in the country (View Fig 2).
The clustering algorithm determined an optimum number of 6 regions based on the minimum David-Boulin star index value (0.74) and a silhouette index located in the first quantile of the index distribution (Silhouette index for 6 clusters: 0.46, maximum Silhouette index: 0.5 for 12 clusters) [49]. The identified clusters are shown in Fig A in S2 Text, defining the regions: 1. South-west, 2. Andean-Pacific, 3. Orinoco-Amazonian piedmont, 4. Central Amazonas, 5. Eastern Orinoco plains, and 6. Caribbean coast.
We fitted our compartmental modeling scheme for every region after defining these six regions with similar rainfall patterns. Interestingly, a diverse pattern of associations between rainfall and incidence throughout the 6 regions of the country was found. For example, we found that rainfall does not drive snakebit incidence in three regions (region 1, 2, and 4) (View Fig 3). Regions 1 and 4 did not exhibit a seasonal component on its incidence (Model 1), while region 2 has a rainfall-independent seasonal component (Model 3). On the other hand, the incidence in regions 3, 5 and 6 is associated with rainfall-driven seasonality (Model 4), indicating that the seasonal component of rainfall signal is the one that modulates snakebite incidence dynamics (View Fig 4). The AIC values of the adjusted models for each region are shown in Table 1.
Rainfall does not modulate snakebite incidence in three regions of the country, region 1, region 2, and region 4. A rainfall-independent seasonal component (Model 3) modulates snakebite incidence in region 2. This component depends on 4 B-splines (Dashed lines in climatic data plot in region2), not related to rainfall. The best model for regions 1 and 4 was Model 1, where incidence does not have seasonality and is not associated with rainfall. Note the difference between the median of the simulations (dashed red line) for region 2 and regions 1 and 4, wherein region 2 this median has seasonal peaks while in region 1 and 4 is relatively constant. Base map of national boundaries of Colombia was obtained from DIVA-GIS free spatial data (https://www.diva-gis.org/datadown).
Rainfall modulates snakebite incidence in three regions of the country, region 3, region 5, and region 6. For all regions, the best model was Model 4 that includes rainfall-driven seasonality as covariate (Dashed light-blue lines in climatic data). Note how all the medians for the simulations (dashed red lines) have seasonal peaks throughout the years, showing the association between rainfall-driven seasonality and snakebite incidence. Also, it is evident how incidence in region 6 is significantly higher than in the other regions. Base map of national boundaries of Colombia was obtained from DIVA-GIS free spatial data (https://www.diva-gis.org/datadown).
We found that the confidence intervals of the parameters determine the association between rainfall and snakebite incidence (View Table A in S2 Text). The confidence intervals for the parameter K in Model 2R (model with rainfall as a covariate) are undetermined in regions where the best model was Model 1, which cancels the type III functional response (View Table A in S2 Text). Finally, we found that in regions where rainfall-driven seasonality modulates snakebite incidence (region 3, 5 and 6), cases decrease during the dry season, and these regions have the minimum monthly rainfall (region 3: 7.02 mm, region 5: 14.88 mm, region 6: 2.73 mm) during the dry season (View Fig 4).
Discussion
Our study explores the role of rainfall and temperature modulating snakebite incidence dynamics in different Colombian regions. We found that envenoming seasonality is significantly explained by rainfall-driven seasonality at the national level, but not by temperature. However, the association between incidence and rainfall is only present in certain regions with marked seasonal rainfall patterns. Our results suggest that national snakebite incidence has two seasonal peaks modulated by rainfall-driven seasonality, the first occurring between April and June, and the second around October (View Fig 2). This association at the national level can be explained in region 6 (Caribbean coast and Low-Magdalena region) contributing to 54% of the national cases: In this region, incidence also exhibited a significant association with rainfall-driven seasonality (View Fig 4). Therefore, inferring snakebite dynamics based only on national results should be done carefully: In region 2 (Andean and pacific regions), the peak of envenomings occurs during the beginning of the year, and they are modulated by rainfall-independent seasonality, in contrast with the national trend (View Fig 3). Contrary to other studies about snakebite seasonality that have used national incidence data [15,16], our study compares snakebite risk between different regions by analyzing hypotheses represented by compartmental models: National-aggregated analyses can largely neglect the geographical heterogeneity in the association between snakebite and temporal drivers, as it was shown in [40].
Spatial distribution of snakebite risk in Colombia
By looking at the null model (No rainfall, temperature, nor rainfall-independent seasonality as covariates, Model 1), the parameter β represents a “constant” snakebite risk or encounter frequency with venomous snakes as described in [27]. Therefore, the likelihood profile over this parameter (View Table A in S2 Text) and its confidence intervals (View Fig B in S2 Text) can be used to compare risk between regions: The regions with the highest risk (region 1, 3, 4, and 5) have the presence of B. atrox [33,53], while regions with B. asper (regions 2 and 6) are the regions with the lowest risk. This effect can be caused by ecological differences between both species, which makes B. atrox more dangerous or abundant than B.asper, or by economic or sociological differences that can increase the exposure to snakebite of inhabitants of these regions. Sadly, biological information about venomous species in the country is scarce, so it is difficult to determine the cause of these high-risk regions [11,30,33,53,54].
Temporal patterns of snakebite incidence in Colombia
We found that regions with the highest coefficient of variation for rainfall, which determines how seasonal is the temporal pattern of precipitation (regions 3, 5, and 6, View Table B in S2 Text), have their reported incidence modulated by rainfall-driven seasonality, where incidence decreases during dry seasons (Model 2R explained snakebite incidence better than Model 1). Thus, a strong seasonality on rainfall, characterized by marked dry and rainy seasons, determines the association between snakebite incidence and rainfall in Colombia. Therefore, we propose that rainfall acts as a limiting resource in snakebite dynamics in Colombia: The association between snakebite incidence and rainfall in Colombia is mediated via marked dry seasons, which cause a decrease in venomous snakes’ activity, abundance, and finally, snakebite risk [37,55].
Another important time series seasonality measurement is the strength of seasonality, which compares the variance in the noise component of the time series with the variance in the seasonal component [33]. Rainfall seasonal strength was considerably lower in regions where we did not find an association between precipitation and incidence; consequently, the rainfall signal in these places is noisy and does not have a clear temporal periodic trend (View Table B in S2 Text and Fig C in S2 Text). It is important to clarify that region 1 has a part located in the Pacific versant, and another is located in the Amazonian versant, which are two distinct ecological regions with different diversity of venomous snakes, but a similar rainfall pattern [33,53,56]. Thus, these differences in venomous snakes’ diversity and ecology can explain the high noise over incidence time-series in this place. We recommend investigating more deeply the association between rainfall and incidence in this region, where data from Ecuador, which is the neighboring country to this area and has the same species composition in the Pacific and Amazonian versant [53], can help to clarify the dynamics of snakebite incidence in this region.
The mechanisms behind the association between rainfall and snakebite incidence patterns are still unclear. However, we want to propose three not mutually exclusive hypotheses to explain this binding: i) Several neotropical snakes, including some species of the genus Bothrops, have their reproductive cycle related to precipitation pattern: Gravid females give birth to neonates at the beginning of the rainy season, thus increasing the abundance of venomous snakes and the probability of encounter between humans and venomous snakes [10,11]. For example, in Costa Rica, the reproductive cycle of Bothrops asper is known, and the seasonality of snakebite incidence is driven by the population dynamics of this species [11,15]. These dynamics can explain the rainfall-independent seasonality of incidence found in region 2, and the association between incidence and rainfall-driven seasonality found in region 6. Nevertheless, Bothrops asper populations in Colombia are genetically different from populations in Costa Rica, so population dynamics between both populations may vary [57]. In addition, it is known that in Brazilian Amazonas, the reproductive cycle of Bothrops atrox is not seasonal, where births occur during most of the year [12]. This can explain why snakebite incidence is not associated with precipitation or rainfall-independent seasonality in region 4 (central Amazonas), but for regions 3 and 5, where Bothrops atrox is also present, the incidence is modulated by rainfall-driven seasonality. ii) Precipitation can affect the ecology of venomous snakes, either by causing floods which decrease the not flooded area that snakes and humans share, therefore raising the contacts between both populations, or by increasing ecosystem productivity: More prey will be available, so snakes could be more active thus increasing snakebite risk [11]. iii) During rainy seasons, agricultural and cattle productivity increases, causing an increase in the number of farmer workers at risk of encountering a venomous snake [58]. Given that ecological information about venomous snakes is non-existent in Colombia, fieldwork must be done to determine how these three rainfall-related events affect snakebite incidence temporal patterns and determine risky seasons.
Final remarks
Our modeling framework can describe the dynamics of snakebite incidence in Colombia. We believe this modeling framework can be used easily in other countries to monitor snakebite incidence and to improve disease management under changing environments. It is crucial to account for regions with different precipitation patterns to determine spatially heterogeneous snakebite dynamics, which was also addressed recently in [40] by statistical modeling. These results demonstrate how countries affected by snakebite can determine in which regions and in which time they need more antivenom, and how to distribute this scarce resource more accurately. This strategic distribution can help decreasing disease burden because most places affected by snakebite have a deficit in antivenom coverage [4,5,7,19]. In addition, determining the ecological mechanism behind model-estimated snakebite temporal patterns can help to develop prevention strategies to decrease the burden of this severe NTD.
We used a compartmental model to understand the temporal patterns of snakebite incidence by using as explanatory variables rainfall, temperature, and a rainfall-independent seasonal component. Given that these models are robust and their structure is based on the available knowledge behind the causes of the disease, several modifications, hypotheses and extrapolations can be tested: We encourage researchers to use this modeling scheme to understand snakebite epidemiology in other countries. For example, in contrast with statistical modeling, in our compartmental models, the estimated value for the parameters represents characteristics of the ecology of the disease (i.e., encounter frequency β* has units of the number of encounters with venomous snakes that ends in a snakebite per unit of time). Even so, our approach is a first step in the development of a better understanding and prediction of snakebite epidemiology by mathematical modeling: Epidemiological models are more refined models that explicitly use information on snake abundance instead of the mentioned environmental proxies, which can deal with the correlation-causality ambiguity that is still present in our approach.
Epidemiological models have been used in several neglected tropical diseases caused by zoonosis, where the basic biology of the hosts involved in disease transmission is known and modelled [59–62]. Thanks to this synergy between mathematics and biology, more specific control programs have been proposed to control the disease burden over an infected population [29,63–65]. Thus, as it has been proposed before, the role that the biology of venomous snakes plays behind snakebite epidemiology is important, but it is still neglected: We want to encourage the study of the natural history of venomous snakes to fill this enormous vacuum of information that limits the understanding of snakebite, and which can be used to generate more robust epidemiological models by including venomous snakes population dynamics. Finally, an interdisciplinary approach must be undertaken to decrease the high burden of this NTD, and to help contribute to the WHO target to reduce snakebite fatality by 30% by the end of 2050 [66].
Supporting information
S1 Text. Compartmental models and B-splines.
This section contains the development and explanation of the models. We used four models in our modeling framework to determine the association between climatic covariates and snakebite incidence, which are explained in Fig A.
https://doi.org/10.1371/journal.pntd.0010270.s001
(DOCX)
S2 Text. Settings for parameter estimation and model selection, clustered regions with similar rainfall patterns, confidence intervals after parameter adjustment, and seasonal parameters for rainfall over study areas.
We used a Poisson process to model data gathering, and then we used a particle filtering algorithm in two steps to adjust our modeling scheme to the data. This process is explained here. Fig A shows the results of the regions with similar rainfall patterns after clustering algorithm. Table A shows the confidence intervals for all parameters in all models after particle filtering algorithm. Fig B compares the adjusted encounter frequency for MODEL 1 between all regions, showing that regions where Bothrops atrox is distributed share the highest values. Table B shows seasonal parameters for rainfall in our defined regions and Fig C contains the seasonal distribution of rainfall over clusterized regions.
https://doi.org/10.1371/journal.pntd.0010270.s002
(DOCX)
Acknowledgments
We thank to Juan Daniel Umaña, Daniela Angarita, Norma Forero, Gabriela Navas and Carlos Cruz for the teamwork doing preliminary statistical analyses, Alejandro Feged for his support at the beginning of the project, and Mahmood Sasa and Camila Renjifo for their support during the study analysis.
References
- 1. Gutiérrez JM, Warrell DA, Williams DJ, Jensen S, Brown N, Calvete JJ, et al. The Need for Full Integration of Snakebite Envenoming within a Global Strategy to Combat the Neglected Tropical Diseases: The Way Forward. PLoS Negl Trop Dis. 2013;7: e2162. Available: pmid:23785526
- 2. Chippaux JP. Snakebite envenomation turns again into a neglected tropical disease! J Venom Anim Toxins Incl Trop Dis. 2017;23: 1–2. pmid:28804495
- 3.
WHO. WHO | Prevalence of snakebite envenoming. WHO. 2017 [cited 18 Oct 2018]. Available: http://www.who.int/snakebites/epidemiology/en/
- 4. Gutiérrez JM, Calvete JJ, Habib AG, Harrison RA, Williams DJ, Warrell DA. Snakebite envenoming. Nat Rev Dis Prim. 2017;3: 17063. pmid:28905944
- 5. Gutierrez JM. Improving antivenom availability and accessibility: science, technology, and beyond. Toxicon. 2012;60: 676–687. pmid:22781134
- 6. Gutiérrez J. Global Availability of Antivenoms: The Relevance of Public Manufacturing Laboratories. Toxins (Basel). 2018;11: 5. pmid:30586868
- 7. Kasturiratne A, Wickremasinghe AR, de Silva N, Gunawardena NK, Pathmeswaran A, Premaratna R, et al. The Global Burden of Snakebite: A Literature Analysis and Modelling Based on Regional Estimates of Envenoming and Deaths. Winkel K, editor. PLoS Med. 2008;5: e218. pmid:18986210
- 8.
World Health Organization (WHO). Snakebite envenoming. In: World Health Organization (WHO) [Internet]. 2021. Available: https://www.who.int/news-room/fact-sheets/detail/snakebite-envenoming
- 9. Ediriweera DS, Kasturiratne A, Pathmeswaran A, Gunawardena NK, Jayamanne SF, Lalloo DG, et al. Health seeking behavior following snakebites in Sri Lanka: Results of an island wide community based survey. PLoS Negl Trop Dis. 2017;11: e0006073. pmid:29108023
- 10. Solórzano A, Cerdas L. Reproductive Biology and Distribution of the Terciopelo, Bothrops asper Garman (Serpentes: Viperidae), in Costa Rica. Herpetologica. 1989;45: 444–450. Available: http://www.jstor.org/stable/3892835
- 11. Sasa M, Wasko DK, Lamar WW. Natural history of the terciopelo Bothrops asper (Serpentes: Viperidae) in Costa Rica. Toxicon. 2009;54: 904–922. pmid:19563822
- 12. Silva KMP, Silva KB, Sueiro LR, Oliveira MEES, Almeida-Santos SM. Reproductive Biology of Bothrops atrox (Serpentes, Viperidae, Crotalinae) from the Brazilian Amazon. Herpetologica. 2019;75: 198–207.
- 13. Goldstein E, Erinjery JJ, Martin G, Kasturiratne A, Ediriweera DS, de Silva HJ, et al. Integrating human behavior and snake ecology with agent-based models to predict snakebite in high risk landscapes. PLoS Negl Trop Dis. 2021;15: e0009047. pmid:33481802
- 14. Chippaux JP. Estimate of the burden of snakebites in sub-Saharan Africa: A meta-analytic approach. Toxicon. 2011;57: 586–599. pmid:21223975
- 15. Chaves LF, Chuang T-W, Sasa M, Gutiérrez JM. Snakebites are associated with poverty, weather fluctuations, and El Niño. Sci Adv. 2015;1. Available: http://advances.sciencemag.org/content/1/8/e1500249.abstract
- 16. Ediriweera DS, Diggle PJ, Kasturiratne A, Pathmeswaran A, Gunawardena NK, Jayamanne SF, et al. Evaluating temporal patterns of snakebite in Sri Lanka: the potential for higher snakebite burdens with climate change. Int J Epidemiol. 2018 [cited 7 Oct 2018]. pmid:30215727
- 17. Shashar S, Yitshak-Sade M, Sonkin R, Novack V, Jaffe E. The Association Between Heat Waves and Other Meteorological Parameters and Snakebites: Israel National Study. J Emerg Med. 2018;54: 819–826. pmid:29661659
- 18. Williams D, Gutiérrez JM, Harrison R, Warrell DA, White J, Winkel KD, et al. The Global Snake Bite Initiative: an antidote for snake bite. The Lancet. Lancet Publishing Group; 2010. pp. 89–91. pmid:20109867
- 19.
Gutiérrez JM, Fan HW. Improving the control of snakebite envenomation in Latin America and the Caribbean: A discussion on pending issues. Transactions of the Royal Society of Tropical Medicine and Hygiene. Oxford University Press; 2018. pp. 523–526. https://doi.org/10.1093/trstmh/try104 pmid:30219842
- 20.
Brauer F, Castillo-Chávez C. Mathematical models in population biology and epidemiology. Springer; 2012.
- 21. Otero M, Solari HG. Stochastic eco-epidemiological model of dengue disease transmission by Aedes aegypti mosquito. Math Biosci. 2010;223: 32–46. pmid:19861133
- 22. Birget PLG, Koella JC. An epidemiological model of the effects of insecticide-treated bed nets on malaria transmission. PLoS One. 2015;10. pmid:26636568
- 23. Erazo D, Cordovez J. Modeling the effects of palm-house proximity on the theoretical risk of Chagas disease transmission in a rural locality of the Orinoco basin, Colombia. Parasites and Vectors. 2016;9: 1–5. pmid:26728523
- 24. Erazo D, Cordovez J. The role of light in Chagas disease infection risk in Colombia. Parasites and Vectors. 2016;9: 9. pmid:26732186
- 25. Cordovez JM, Rendon LM, Gonzalez C, Guhl F. Using the basic reproduction number to assess the effects of climate change in the risk of Chagas disease transmission in Colombia. Acta Trop. 2014;129: 74–82. pmid:24416781
- 26.
Brauer F, Castillo-Chavez C, Feng Z. Mathematical Models in Epidemiology. New York, NY: Springer New York; 2019. https://doi.org/10.1007/978-1-4939-9828-9
- 27. Bravo-Vega CA, Cordovez JM, Renjifo-Ibáñez C, Santos-Vega M, Sasa M. Estimating snakebite incidence from mathematical models: A test in Costa Rica. Bottazzi ME, editor. PLoS Negl Trop Dis. 2019;13: e0007914. pmid:31790407
- 28. Gaff H, Schaefer E. Optimal control applied to vaccination and treatment strategies for various epidemiological models. Math Biosci Eng. 2009;6: 469–92. Available: http://www.ncbi.nlm.nih.gov/pubmed/19566121 pmid:19566121
- 29.
Feng Z. Applications of epidemiological models to public health policymaking: The role of heterogeneity in model predictions. world scientific. World Scientific Publishing Co.; 2014. https://doi.org/10.1142/8884
- 30. Lynch JD. El contexto de las serpientes en Colombia con un análisis de las amenazas en contra de su conservación. Rev Colomb Cienc. 2012;36: 435–449. Available: http://www.scielo.org.co/pdf/racefn/v36n140/v36n140a09.pdf
- 31. Rojas G, Bogarin G, Gutierrez JM. Snakebite mortality in Costa Rica. Toxicon. 1997;35: 1639–1643. pmid:9428110
- 32.
Sant’Ana Malaque CM, Gutiérrez JM. Snakebite Envenomation in Central and South America. Critical Care Toxicology. Cham: Springer International Publishing; 2015. pp. 1–22. https://doi.org/10.1007/978-3-319-20790-2_146–1
- 33.
Campbell JA, Lamar WW. The venomous reptiles of the Western Hemisphere. United States of America: Comstock Pub. Associates; 2004.
- 34. Almeida-Santos SM, Salomao M da G. Reproduction in neotropical pitvipers, with emphasis on species of the genus Bothrops. Biol vipers. 2002;1: 445–462.
- 35. Marques OA V., Kasperoviczus K, Almeida-Santos SM. Reproductive Ecology of the Threatened Pitviper Bothrops insularis from Queimada Grande Island, Southeast Brazil. J Herpetol. 2013;47: 393–399.
- 36. Hartmann MT, Marques OAV, Almeida-Santos SM. Reproductive biology of the southern Brazilian pitviper Bothrops neuwiedi pubescens (Serpentes, Viperidae). Amphib Reptil. 2004;25: 77–85.
- 37. M Oliveira E, Martins M. When and where to find a pitviper: activity patterns and habitat use of the lancehead, Bothrops atrox, in Central Amazonia, Brazil. Herpetological Natural History. 2001.
- 38. Fraga R de, Magnusson WE, Abrahão CR, Sanaiotti T, Lima AP. Habitat Selection by Bothrops atrox (Serpentes: Viperidae) in Central Amazonia, Brazil. Copeia. 2013;2013: 684–690.
- 39. Reza AS, Islam SN, Ghosh A, Rahman MR, Faiz M, Rahim A. Magnitude of Snake Bite and Drowing During Monsoon Flood in Two Districts of Bangladesh. J Shaheed Suhrawardy Med Coll. 2017;7: 3–5.
- 40. Ediriweera DS, Kasthuriratne A, Pathmeswaran A, Gunawardene NK, Jayamanne SF, Murray K, et al. Evaluating spatiotemporal dynamics of snakebite in Sri Lanka: Monthly incidence mapping from a national representative survey sample. PLoS Negl Trop Dis. 2021;15: e0009447. pmid:34061839
- 41.
Arango C, Dorado J, Guzmán D, Ruiz J. Climatología trimestral de Colombia. Bogotá, D.C: Instituto de Hidrología, Meteorología y estudios ambientales IDEAM; Available: http://www.ideam.gov.co/documents/21021/21789/Climatología+Trimestral+para+Colombia+%28Ruiz%2C+Guzman%2C+Arango+y+Dorado%29.pdf/c2825963-c373-449a-a7cb-8480874478d9
- 42.
María A, Bárcenas R. INFORME DE EVENTO ACCIDENTE OFÍDICO, COLOMBIA, 2017. Dirección de Vigilancia y Análisis del Riesgo en Salud Pública SIVIGILA; 2018. Available: https://www.ins.gov.co/buscador-eventos/Informesdeevento/ACCIDENTEOFÍDICO 2017.pdf
- 43. Nuñez León LJ, Camero-Ramos G, Gutierrez JM. Epidemiology of snakebites in Colombia (2008–2016). Rev Salud Pública. 2020;22: 1–8.
- 44. King AA, Nguyen D, Ionides EL. Statistical inference for partially observed markov processes via the R package pomp. J Stat Softw. 2016;69: 1–43.
- 45.
R Development Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria; 2008. Available: http://www.r-project.org
- 46. Cleveland RB, Cleveland W. S., McRae J.E. STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J Off Stat. 1990;6: 3–73. Available: https://www.wessa.net/download/stl.pdf
- 47. Abatzoglou JT, Dobrowski SZ, Parks SA, Hegewisch KC. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci Data. 2018;5: 1–12. pmid:30482902
- 48.
Paparrizos J, Gravano L. K-shape: Efficient and accurate clustering of time series. Proceedings of the ACM SIGMOD International Conference on Management of Data. New York, New York, USA: Association for Computing Machinery; 2015. pp. 1855–1870. https://doi.org/10.1145/2723372.2737793
- 49. Liu Y, Li Z, Xiong H, Gao X, Wu J. Understanding of Internal Clustering Validation Measures. IEEE Int Conf Data Min. 2010; 911–916.
- 50. Sardá-Espinosa A. Time-series clustering in R Using the dtwclust package. R J. 2019;11.
- 51. Bozdogan H. Model selection and Akaike’s Information Criterion (AIC): The general theory and its analytical extensions. Psychometrika. 1987;52: 345–370.
- 52. Stocks T, Britton T, Höhle M. Model selection and parameter estimation for dynamic epidemic models via iterated filtering: application to rotavirus in Germany. Biostatistics. 2020;21: 400. pmid:30265310
- 53.
Valencia J, Garzón-Tello K, Barragán-Paladínes M. Serpientes venenosas del Ecuador: sistemática, taxonomía, historia natural, conservación. Envenenamiento y aspectos antropológicos. 1st ed. Quito: Fundación Herpetológica Gustavo Orcés; 2016. Available: http://biblioteca.udla.edu.ec/client/en_US/default/search/detailnonmodal?qu=Barragán-Paladines%2C+María+Elena&d=ent%3A%2F%2FSD_ILS%2F0%2FSD_ILS%3A29363~~0&ic=true&te=ILS&ps=300
- 54.
Lynch JD, Angarita-Sierra T, Ruiz-Gómez FJ. Programa nacional para la conservación de las serpientes presentes en Colombia. Bogotá: Ministerio de Ambiente y Desarrollo Sostenible, Colombia, Universidad Nacional de Colombia—Instituto de Ciencias Naturales, Instituto Nacional de Salud.; 2016.
- 55. Wasko DK, Sasa M. Activity Patterns of a Neotropical Ambush Predator: Spatial Ecology of the Fer-de-lance (Bothrops asper, Serpentes: Viperidae) in Costa Rica. Biotropica. 2009;41: 241–249.
- 56. Sánchez-Cuervo AM, Aide TM, Clark ML, Etter A. Land Cover Change in Colombia: Surprising Forest Recovery Trends between 2001 and 2010. PLoS One. 2012;7. pmid:22952816
- 57. Saldarriaga-Córdoba M, Parkinson CL, Daza JM, Wüster W, Sasa M. Phylogeography of the Central American lancehead Bothrops asper (SERPENTES: VIPERIDAE). Chiang T-Y, editor. PLoS One. 2017;12: e0187969. pmid:29176806
- 58. Esquivel A, Llanos-Herrera L, Agudelo D, Prager SD, Fernandes K, Rojas A, et al. Predictability of seasonal precipitation across major crop growing areas in Colombia. Clim Serv. 2018;12: 36–47.
- 59. Koella J, Antia R. Epidemiological models for the spread of anti-malarial resistance. Malar J. 2003;2: 3. pmid:12643812
- 60. Racloz V, Ramsey R, Tong S, Hu W. Surveillance of Dengue Fever Virus: A Review of Epidemiological Models and Early Warning Systems. Anyamba A, editor. PLoS Negl Trop Dis. 2012;6: e1648. pmid:22629476
- 61. Velasco-Hernández JX. An epidemiological model for the dynamics of Chagas’ disease. BioSystems. 1991;26: 127–134. pmid:1841638
- 62. Garchitorena A, Sokolow SH, Roche B, Ngonghala CN, Jocque M, Lund A, et al. Disease ecology, health and the environment: a framework to account for ecological and socio-economic drivers in the control of neglected tropical diseases. Philos Trans R Soc B Biol Sci. 2017;372: 20160128. pmid:28438917
- 63. Sturrock HJW, Picon D, Sabasio A, Oguttu D, Robinson E, Lado M, et al. Integrated Mapping of Neglected Tropical Diseases: Epidemiological Findings and Control Implications for Northern Bahr-el-Ghazal State, Southern Sudan. Lammie PJ, editor. PLoS Negl Trop Dis. 2009;3: e537. pmid:19859537
- 64. Alison Kealey A, Robert Smith R. Neglected Tropical Diseases: Infection, Modeling, and Control. J Health Care Poor Underserved. 2010;21: 53–69. pmid:20173255
- 65.
Mubayi A. Inferring Patterns, Dynamics, and Model-Based Metrics of Epidemiological Risks of Neglected Tropical Diseases. Handbook of Statistics. Elsevier B.V.; 2017. pp. 155–183. https://doi.org/10.1016/bs.host.2017.09.002
- 66.
World health organization. Snakebite envenoming a strategy for prevention and control. Geneva: World health organization; 2019. Available: http://apps.who.int/bookorders.