Using Geographical Information Systems to Identify

Populations in Need of Improved Accessibility to

Antivenom Treatment for Snakebite Envenoming in
Costa Rica
Erik Hansson1*, Mahmood Sasa2, Kristoffer Mattisson1, Arodys Robles3, José Marı́a Gutiérrez2
1 Occupational and Environmental Medicine, Lund University, Lund, Sweden, 2 Instituto Clodomiro Picado, Facultad de Microbiogı́a, Universidad de Costa Rica, San José,
Costa Rica, 3 Centro Centroamericano de Población, Universidad de Costa Rica, San José, Costa Rica

Abstract
Introduction: Snakebite accidents are an important health problem in rural areas of tropical countries worldwide, including
Costa Rica, where most bites are caused by the pit-viper Bothrops asper. The treatment of these potentially fatal accidents is
based on the timely administration of specific antivenom. In many regions of the world, insufficient health care systems and
lack of antivenom in remote and poor areas where snakebites are common, means that efficient treatment is unavailable for
many snakebite victims, leading to unnecessary mortality and morbidity. In this study, geographical information systems
(GIS) were used to identify populations in Costa Rica with a need of improved access to antivenom treatment: those living in
areas with a high risk of snakebites and long time to reach antivenom treatment.

Method/Principal Findings: Populations living in areas with high risk of snakebites were identified using two approaches:
one based on the district-level reported incidence, and another based on mapping environmental factors favoring B. asper
presence. Time to reach treatment using ambulance was estimated using cost surface analysis, thereby enabling adjustment
of transportation speed by road availability and quality, topography and land use. By mapping populations in high risk of
snakebites and the estimated time to treatment, populations with need of improved treatment access were identified.

Conclusion/Significance: This study demonstrates the usefulness of GIS for improving treatment of snakebites. By mapping
reported incidence, risk factors, location of existing treatment resources, and the time estimated to reach these for at-risk
populations, rational allocation of treatment resources is facilitated.

Citation: Hansson E, Sasa M, Mattisson K, Robles A, Gutiérrez JM (2013) Using Geographical Information Systems to Identify Populations in Need of Improved
Accessibility to Antivenom Treatment for Snakebite Envenoming in Costa Rica. PLoS Negl Trop Dis 7(1): e2009. doi:10.1371/[Link].0002009
Editor: Jean-Philippe Chippaux, Institut de Recherche pour le Développement, Benin
Received August 22, 2012; Accepted November 30, 2012; Published January 31, 2013
Copyright: ß 2013 Hansson 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.
Funding: The Swedish Society of Medicine supported this study by granting a travel grant for Erik Hansson to go to Costa Rica to perform the study. This grant
consisted of approximately $2000. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The
other authors did not receive any specific funding.
Competing Interests: The authors have declared that no competing interests exist.
* E-mail: [Link]@[Link]

Introduction is the case of Costa Rica, where antivenom is widely available

throughout the public health system, as a result of domestic
Snakebites are a health problem in several parts of the world production and effective acquisition and distribution schemes [5].
[1], especially affecting poor people living in rural areas in tropical The incidence of snakebites is known to vary widely on a sub-
countries [2]. The treatment of snakebite envenoming is based on national level due to environmental and demographic factors [10–
the timely administration of animal-derived antivenoms [3]. 14]. In Central America, an estimated number of 4,000 snakebite
However, a number of factors limit the accessibility of antivenoms cases occur every year [15]. The Central America lancehead
in various areas of the world [3,4]. The relatively high cost of some pitviper Bothrops asper, locally known as ‘terciopelo’ or ‘barba
antivenoms make accessing them difficult in low income countries amarilla’, which is widely distributed in humid, lowland areas (0–
[4]. In addition, long distances to healthcare facilities and incorrect 600 m.a.s.l., but it might be found up to 1,200 m.a.s.l.) [16,17], is
distribution of antivenoms within countries [3,5,6] means that for responsible for 50–80% of the bites and 60–90% of all fatal cases
many of the snakebite victims worldwide, specific treatment with in Central America and northern South America [13]. Epidemi-
safe and effective antivenoms is unavailable, and traditional ological studies indicate that most victims are agricultural workers
healers are instead often consulted [1,7–9]. This is an unfortunate and/or rural residents in general [13]. In the 1990s, the snakebite
situation as antivenom treatment is highly effective at preventing incidence in Costa Rica was approximately 15 per 100,000
morbidity and mortality caused by snakebite envenoming [3]. In inhabitants per year, with a decreasing trend [14]. Snakebite
other regions of the world, such as in various countries in Latin mortality has also been reported to decrease in Costa Rica, an
America, antivenom is readily available in health centers [5]. This observation attributed to, among other things, improvements in

PLOS Neglected Tropical Diseases | [Link] 1 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Author Summary Objectives

The primary aim of the present study is to provide information
Snakebite accidents are a neglected cause of much death to assist decision-making concerning for which primary health care
and suffering worldwide. The situation is especially severe facilities (EBAIS) it is suitable to have antivenom, i.e. those that
in rural areas of low income tropical countries, where long serve a population with a high risk of snakebites and long transport
distances to well-equipped health care facilities mean that times to hospitals or clinics where antivenom is available. The
the time needed to reach antivenom treatment is often
methodology used in this study might be applied in other countries
long. Delays in receiving antivenom treatment of snake-
and regions where snakebite envenoming is a relevant public
bites could lead to severe outcomes, such as death or
permanent disability. In this study we demonstrate how health hazard, in order to identify vulnerable areas that require
Geographical Information Systems (GIS) could be used to interventions aimed at ensuring the access to antivenom and
allocate antivenom rationally and thereby decrease the adequate medical treatment. As part of the primary aim, we will
impact of snakebite in a cost-effective manner. We map demonstrate the benefits of spatial smoothing of the small-area
areas with a high risk of snakebite accidents, based on a snakebite incidence data for increasing the interpretability of this
high reported incidence and environmental conditions data. Aside from the primary aim, we will also describe the
favoring snakebites. We then estimate the time needed to relationship between some district-level environmental and
reach a facility in which antivenom treatment is available demographic factors and snakebite incidence in Costa Rica.
for the population in these high risk areas, thus identifying
areas in need of improved treatment accessibility. Based Materials and Methods
on these maps of the unmet need of antivenom treatment,
allocation of antivenom and other resources needed to Outline of the analysis
treat snakebites can be made more efficiently. Populations in a high risk of snakebite were identified using two
approaches: one based on the reported district-level snakebite
incidence, and another based on identifying populations with
the health care system, including better geographical accessibility environmental and demographic risk factors favoring snakebites.
to health care facilities [12,18]. The use of first aid interventions For the first approach, an incidence of 30 bites per 100,000
and the attendance to traditional healers for treatment of population per year was selected as threshold, on the basis of the
snakebites is very limited in Costa Rica, at least among those overall incidence of snakebite envenomings in Costa Rica (15 cases
seeking hospital treatment [19]. per 100,000 inhabitants per year [14]). In order to reduce the
Identification of areas with high snakebite incidence is an random noise in the small-area incidence, a Bayesian smoothing
important goal for the design of preventive and therapeutic method [26] was employed by fitting a Poisson regression to the
interventions aimed at reducing the impact of snakebites. Such observed incidence in the period 2003–2007. Smoothing can be
interventions could include educational programs, antivenom explained as modelling the underlying risk around which the
distribution and training of health staff in the management of these
observed incidence varies stochastically. As such the smoothed
envenomations [4–6]. Similar issues of rational distribution of
incidence estimates give a more stable picture of the actual risk of
treatment and intervention resources are being considered in the
snakebites than the raw incidence data, making interpretation
struggle against other neglected tropical diseases, with geograph-
easier and giving a better basis for decision-making [26]. Spatial
ical information systems (GIS) increasingly being used as a tool for
smoothing was also used by Leynaud and Reati [21] in their study
improved decision-making [20]. Leynaud and Reati [21] de-
of snakebite epidemiology in Argentina. In this study, we present
scribed the use of the free GIS-software SIGEpi [22] for
an alternative smoothing method that includes also measured
identification of areas with a high risk of snakebites and long
explanatory variables, and allows for estimation of the probability
distances to healthcare in Argentina. It is necessary to further
that an incidence threshold is exceeded [27]. We then compare the
extend the use of this methodology in other countries, in order to
performance of this smoothing method with the one used by
get an integrated view of snakebite envenomings on various
Leynaud and Reati [21]. Measured variables to be included in the
regions of the world.
smoothing were identified in the literature as important factors
The health care system in Costa Rica is considered well
influencing snake presence, i.e. environmental factors [16], and
developed, in terms of geographical accessibility [23] and
snake-human interaction, in this case urban population percentage
insurance coverage [24]. This country has made important
achievements in public health, such as long life expectancy and [13].
low infant mortality rate [24,25]. Antivenom is available in all For the second approach to identifying areas in possible need of
clinics and hospitals and, according to a recent decision by the improved antivenom accessibility, rural areas in humid, lowland
state social insurance administration that runs all public health areas were considered at high risk of snakebites, as this is a habitat
services in the country (Caja Costarricense del Seguro Social, suitable for the medically most relevant snake species B. asper [16].
henceforth CCSS), it can also be distributed to primary health care Among the populations identified as in high risk of snakebites by
teams (Equipos Básicos de Atención Integral en Salud, EBAIS) of either of these two approaches, those with an estimated minimum
which there are approximately 1 per 5,000 inhabitants [24]. This transportation time to hospital or clinic exceeding 2 hours were
opens the possibility to further improve the accessibility to identified as those requiring improved antivenom accessibility,
antivenoms in Costa Rica. However, the distribution of anti- based on the time identified as critical to avoid mortality from
venoms to EBAIS has to be carefully analyzed, in order to ensure snakebites [13].
that mostly EBAIS serving a population with high snakebite risk
and limited access to clinics and hospitals will receive antivenoms. Material
This will prevent unnecessary wastage of this precious drug by The centroid of the census enumeration units were used to
deploying it to regions where snakebites are infrequent, or where estimate the location of the population [23]. A 2 km buffer around
distances to hospitals or clinics where these accidents are treated the census tract centroids were used as a proxy for populated
are short. areas. Census enumeration unit population size in year 2000 and

PLOS Neglected Tropical Diseases | [Link] 2 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

information whether the unit was urban or rural was available as coverage, and biotic unit [28]) were obtained from a database
attribute information. compiled by Instituto Tecnológico de Costa Rica [29]. The values
Information about snakebite cases 1990–2007 per district of for each of the environmental variables were extracted to the
residence were obtained from hospital discharge reports from the census enumeration unit centroids, and the population-weighted
CCSS. Duplicate cases, with exactly the same age, sex, month of district average value for each of the environmental variables was
bite and district of residence (n = 123) were removed from the calculated. For forest coverage, the proportion of inhabitants living
analysis. Cases reported from health care facilities in which it was within 500 meters from forests larger than 5 ha was calculated. A
highly unlikely that the place of the bite was in the district of few districts (n = 9) were covered by clouds in the satellite image
residence (n = 141), or when the place of residence was unknown from which the forest coverage was derived, and these were
(n = 43) were also removed. These steps left 9,149 snakebites, assigned the mean forest coverage value of the other districts.
divided over 413 districts (the smallest administrative area) The location of lakes, rivers, clinics, hospitals and roads were
(Figure 1). The location of B. asper specimens (n = 241) collected obtained from the same database as the environmental variables
to the serpentarium at Instituto Clodomiro Picado, San José, [29]. A list of primary health care facilities, i.e. Equipos Básicos de
Costa Rica was obtained. Atención Integral en Salud (EBAIS) was obtained from the CCSS
Digital maps of environmental variables (elevation 30630 m [30]. EBAIS were geocoded by matching the name and service
raster, mean annual precipitation, number of dry months, forest area of the facility with the name and district of communities in

Figure 1. Snakebite incidence per 100,000 inhabitants per district 1990–2007.

doi:10.1371/[Link].0002009.g001

PLOS Neglected Tropical Diseases | [Link] 3 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

the country, available as a digital map [29]. In most cases, there The first smoothing method (A) was the same Bayesian model
was a community with the same name as the facility. If there was with random spatial effects and fixed effect explanatory variables
no such community, the facility was located in the main (elevation, number of dry months, precipitation, forest coverage
community of the district. If it was not possible to determine and urban population percentage) as described above. The
which the main community was, the facility was located in one of Bayesian framework makes it possible to estimate the probability
the communities located in the center of the service area. Red that the incidence exceeds a threshold directly from the posterior
Cross ambulance stations were geo-coded using a list and map distribution. For each of the training time periods, the models were
available at the Costa Rican Red Cross Website [31]. run until convergence was reached (after 50,000 iterations) and
then for another 50,000 iterations to obtain samples from the
A Bayesian Poisson regression model of snakebite posterior distribution. From the posterior distribution, the
incidence probability that the incidence exceeded 30 bites per 100,000
inhabitants was calculated.
In order to analyze which factors were important for snakebite
The second smoothing method (B) was the one used by
occurrence, and to smooth the snakebite incidence, a Bayesian
Leynaud and Reati [21]; the tool for automated spatial Bayesian
Poisson regression was used to model the district-level risk of
smoothing of incidence rates (‘‘Suavizador espacial de tasas’’)
snakebites. Factors chosen for inclusion were those that had
available in SIGEpi [22]. The settings of local smoothing and
previously been identified in literature as either influencing snake
neighborhood defined as adjacency were used.
presence [16] or important for snake-human interaction patterns
The abilities of these smoothed estimates, and the unsmoothed
[13]. The number of snakebite cases Y per district i were modeled
incidence, to identify whether the incidence in a district would
as Poisson variates in the form;
exceed 30 bites per 100,000 in the next-coming five-year period
(i.e. the test periods 1995–1999, 1999–2003 and 2003–2007) were
Yi Poisson(mi ) assessed using the ROC Curve function of IBM SPSS 20 [34].
The AUC (area under the curve) of the ROC (receiver operating
characteristics; a plot of sensitivity vs. one minus specificity) is an
X
K   often used tool to assess the discriminatory ability of tests; an AUC
popi
logðmi Þ~az bk Xi zui z log Yall | of 0.7–0.9 indicate reasonable discriminatory ability, and .0.9
k~1
popall very good discriminatory ability [35]. For each of the six
combinations of training time (1 and 5 years) and smoothing
P
K
method (A, B and unsmoothed), the mean AUC (with 95%
where a is the intercept, bk Xi a matrix of five fixed effect
k~1 empirical confidence intervals (C.I.)) were computed by simulating
district-level explanatory variables (urbanity, forest coverage, plausible AUC values from the uncertainty interval for the three
elevation, precipitation, and number of dry months), ui a spatially training-test period pairs and calculating the mean of these.
correlated random effect modelled using a conditional autore-
gressive (CAR) prior structure [32], which assumed dependence Identification of high-risk areas
between districts if they shared a border or corner, and In order to produce smoothed estimates of the underlying
popi snakebite risk in a district, that would be less affected by random
Yall | correspond to the number of snakebites that would
popall noise, and thereby able to more precisely identify the need for
be expected if they were distributed evenly within the population, antivenom accessibility in that district in the upcoming years, the
i.e. the offset. In order to facilitate convergence, continuous above model was applied to the most recently available 5-year-
variable (precipitation, number of dry months and elevation) were period of the snakebite data; i.e. 2003–2007. After reaching
standardized to have mean 0 and standard deviation 1. The convergence after 50,000 iterations, the model was run for another
coefficients for the fixed effects were assigned non-informative 50,000 iterations, during which the probability that the smoothed
normal distribution priors (mean 0 and precision 0.0001), and the incidence exceeded the threshold of 30 bites per 100,000
intercept a non-informative flat prior (range 2‘ to ‘). The inhabitants was calculated (Figure 2).
variance of the spatially correlated random effects was assigned a The incidence threshold exceedance probability estimate that
non-informative gamma prior. The model was fitted in WinBUGS corresponded to 90% sensitivity in detecting districts with an
1.4.3 [33]. incidence above the 30 bites per 100,000 inhabitants threshold in a
When the district-level risk factors of snakebites were analyzed, all future 5-year period was on average 10% for three earlier 5-year
snakebite data (i.e. 1990–2007) were used. Convergence was reached periods (data not shown). Therefore, this was set as the cut-off
after 50 000 iterations, after which another 100 000 iterations were probability for identifying high-risk districts, in need of good
performed to estimate the posterior distribution, from which model antivenom accessibility. Rural population residing in an environ-
parameters with 95% credible intervals (Cr.I.) were obtained. ment suitable for the medically most relevant snake species B. asper
(i.e. below 1200 meters of elevation and in Moist, Wet or Pluvial
Demonstration of spatial smoothing of small-area biotic unit [16,28,29]) (Figure 3) were also identified as living in
snakebite incidence data high-risk areas (Figure 4).
The snakebite incidence dataset was divided into nine temporal
periods, six training periods on which the smoothing was Identification of high-risk areas in need of improved
performed, and three test periods used to assess the ability of the treatment accessibility
estimates produced by the smoothing methods to improve The mountainous terrain of Costa Rica imposes strong
identification of future high-incidence districts. Three of the restrictions on human movement. The Euclidean distance (straight
training periods were five-year periods (1990–1994, 1994–1998 line) approximation of the time needed to move from the place of
and 1998–2002) and three one-year periods (1994, 1998 and snakebite to the healthcare facility (hospital or clinic) might
2002).The three test periods were five years (1995–1999, 1999– therefore not accurately capture the real time spent in this
2003 and 2003–2007). transportation. We aimed at constructing a model of the time

PLOS Neglected Tropical Diseases | [Link] 4 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Figure 2. The probability that the smoothed incidence exceeded 30 bites per 100,000 inhabitants in 2003–2007. Inc = incidence.
doi:10.1371/[Link].0002009.g002

needed to reach antivenom treatment, i.e. hospitals or clinics, distance is not covered by ambulance transportation and it will be
which takes into account the availability and quality of roads, counted twice (see below). Off-road speed in areas covered by
topography, land cover and the mode of transportation used to forest was assumed to be 50% slower. Streams, rivers and lakes,
reach healthcare. We assumed that people would choose that were not crossed by roads, were given a speed one fifth of that
ambulance services after a snakebite accident. Travel time was in open terrain, in order to penalize movement across water. From
therefore calculated as time with ambulance from closest the elevation raster, a slope raster was constructed. The road/off-
ambulance station, to place of residence, and from there to closest road speed raster, and the slope raster were combined to yield a
hospital or clinic. The time estimated should be considered a raster of the time needed to travel one raster cell (tr,a) using
minimum, ideal, time as it is modeled assuming that the Equation 1.
ambulance leaves immediately after the snakebite to meet up the
snakebite victim, without any delays.
The time to treatment was estimated using GRASS 6.4.1 [36]. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
The road data vector layer was converted to a 30630 m raster and tr,! ~tr | 1z tanð!Þzbr z!2 ðEquation1Þ
classified according to the following assumed speeds: primary pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
roads, 60 km/h; secondary and urban roads, 40 km/h; and The first part ( 1ztanð!Þ) of Equation 1 calculates the excess
tertiary, local and other roads, 20 km/h. Off-road speed was set to distance needed to travel the cell due to change in altitude (i.e. the
6 km/h, but the off-road speed will become 3 km/h as this hypotenuse). The second part (br  !2 )of the equation aims at

PLOS Neglected Tropical Diseases | [Link] 5 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Figure 3. Geographical extent of habitat suitable for Bothrops asper and location of specimens collected at ICP. ICP = Instituto
Clodomiro Picado.
doi:10.1371/[Link].0002009.g003

taking into account the reduced speed associated with moving in 2–3 h and .3 h, in order to visualize the need of improved
undulated terrain, as well as to create barriers in extremely steep antivenom accessibility (Figure 6). Finally, the location of EBAIS
slopes. In this part of the equation, we assume that the excess time and roads were added to the map of populations in need of
needed to travel one of the 30 m cells is proportional to the improved antivenom accessibility, to develop small-scale maps of
squared slope ! (in degrees) multiplied by a coefficient br (0.001 areas in need of improved accessibility to antivenom, which can be
for road travel and 0.02 for off-road). The value of this co-efficient used to guide such improvements (Figure 7 and Supporting
was chosen based on simulation of what values produced Information S1).
reasonable estimations.
The time needed to go from an ambulance station to any cell, Results
and the time to go from a hospital or clinic to any cell was
calculated. These two time rasters were then summed to give a The non-spatial descriptive statistics of major parts of our
raster of the total time needed to reach healthcare facilities using dataset correspond very closely to what has been described in
ambulance, i.e. from ambulance station to snakebite victim and detail previously [14], and thus only a brief summary is presented
from there to healthcare (Figure 5). The mean time to reach in this article. A total of 9,333 cases were reported for the period
healthcare (Figure 5) from populated areas in a high risk of 1990–2007, corresponding to an average incidence of 13.8
snakebites (Figure 4) was extracted, and categorized into ,2 h, snakebites per 100,000 inhabitants per year. Seventy-two percent

PLOS Neglected Tropical Diseases | [Link] 6 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Figure 4. Populated areas in high risk of snakebites. Populated areas = areas within 2 km from a census tract centroid in a high risk of
snakebites, i.e. in districts with a high snakebite incidence, and/or in an area suitable for B. asper.
doi:10.1371/[Link].0002009.g004

of the victims were male. Regarding age, 27% of victims were snakebite incidence in the next five-year period (Table 2), but
below 15 years of age, 32% between 15 and 30 years, 21% discriminatory ability was clearly improved by using smoothing.
between 30 and 45 years, and the remaining 22% corresponded to There were only minor differences between the estimates
people older than 45 years. produced by the different smoothing methods A and B in
discriminatory performance. When the training time was only
District-level risk factors of snakebites one the incidence threshold exceedance probability of method A
In multivariate analysis, lowland districts with much precipita- had a tendency to be better than the smoothed incidence
tion and few dry months generally had a higher snakebite produced by method B (p = 0.07). When the training time was
incidence, as did districts with many rural inhabitants, and many five years, both methods had a borderline significantly better
inhabitants residing close to larger forests, although the last discriminatory ability than the unsmoothed incidence observed in
explanatory factor did not reach 95% significance (Table 1). five years (p = 0.02 for method A and p = 0.10 for method B).

Demonstration of spatial smoothing of small-area Identification of high-risk districts

snakebite incidence data There were marked geographical differences in the reported
The unsmoothed incidence observed in a one-year period had incidence, with high rates especially along the southern Pacific
an acceptable ability to identify which districts would have a high coast, and in parts of the Caribbean and northern lowlands

PLOS Neglected Tropical Diseases | [Link] 7 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Figure 5. Clinics, hospitals, ambulance stations and estimated time to hospital or clinic.
doi:10.1371/[Link].0002009.g005

(Figure 1). The districts with a high incidence corresponded to

Table 1. District-level risk factors for snakebite incidence in a large extent to the humid lowland areas identified as suitable
multivariate Poisson regression. habitat for B. asper (Figure 4); 91% of the rural population in
high-incidence districts lived in this type of environment.
However, only 51% of the rural population in humid lowland
95% 95% areas lived in districts with a high snakebite incidence. According
Incidence Cr.I. Cr.I. to the population distribution of the 2000 census, 18% of the total
Covariate Unit Ratio Low High population lived in districts with a high snakebite incidence and an
Elevation +550 m 0.60 0.48 0.75 additional 13% lived in rural, humid lowlands in districts not
Dry months +1.25 months 0.66 0.51 0.86
reporting an incidence above the threshold.
Precipitation +827 mm/year 1.45 1.13 1.80
Identification of high-risk areas in need of improved
Population living within +100%-units 1.53 0.88 2.59
500 m from .5 ha forest treatment accessibility
Urban population +100%-units 0.39 0.29 0.53
Hospitals, clinics, ambulance stations and roads were strongly
aggregated in the great metropolitan area around the cities of
doi:10.1371/[Link].0002009.t001 San José, Heredia, Cartago and Alajuela, where a major part

PLOS Neglected Tropical Diseases | [Link] 8 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Figure 6. Time to hospital or clinic for populated areas at high risk of snakebites. Populated areas = areas within 2 km from a census tract
centroid in a high risk of snakebites.
doi:10.1371/[Link].0002009.g006

of the population lives. There is however also a network of roads to identify vulnerable places where access to hospitals or clinics is
and healthcare facilities in the more peripheral regions of the delayed, thus setting the stage for the design of more effective ways
country (Figure 5). Of the population living in areas with a high to guarantee a more rapid access to antivenom treatment at the
risk of snakebites, 92.5% were estimated to have a minimum local level. Similar maps for other target regions are available as a
transportation time of less than 2 hour to hospitals or clinics, 5% supplement (Supporting Information S1).
were estimated to delay 2–3 hours to hospitals or clinics, and 2.5%
more have transportation times higher than 3 hours. On the south Discussion
Pacific coast around Golfo Dulce, around the Talamanca
highlands in the southeast, and along the northern border, there Key results
are populations in high risk of snakebites and with long The spatial distribution of district-level snakebite incidence in
transportation times to antivenom treatment (Figure 6). Figure 7 Costa Rica largely followed the expected pattern, based on
provides a close-up view of one target region (Golfo Dulce), previous studies and on the distribution of the most important
including all the information presented in the previous maps, as venomous snake in the country, B. asper. Incidence was higher in
well as the location of primary health care facilities (EBAIS). This rural, humid lowlands, notably in the southern part of the country
map demonstrates the type of map that can be prepared and used [14]. Geographical accessibility to antivenom treatment was

PLOS Neglected Tropical Diseases | [Link] 9 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Figure 7. Close-up view of the Golfo Dulce region. Time to reach hospital or clinic for populated areas with a high snakebite risk, health care
facilities and roads. Populated areas = areas within 2 km from a census tract centroid in a high risk of snakebites.
doi:10.1371/[Link].0002009.g007

generally good, however, in some areas there is a need of among snakebite victims receiving treatment at hospitals in 1996
improved treatment accessibility. These are areas where hospitals found that only 2.9% of them had received any type of empirical
and clinics are located relatively far from some of the areas with a treatment before reaching formal healthcare attention [19]. The
likely high snakebite incidence. literacy rate in Costa Rica is 96% [39], and educational campaigns
of various sorts over several decades have raised awareness about
Limitations the importance of seeking formal healthcare after snakebites
The potential bias introduced by using data reported from the among the Costa Rican public [19]. Based on these facts, we
healthcare system to analyze snakebite incidence, due to the use of assume that the degree of utilization of traditional healthcare after
traditional medicine and dysfunctional reporting routines, are well snakebites is very low in Costa Rica and that, on this basis, there
known from other parts of the world [9,37]. However, in contrast will not be much underestimation of the true snakebite incidence
to the situation in many other developing countries in Latin in the statistics available from the Ministry of Health. Nevertheless,
America and elsewhere, traditional medicine is not widely used in our data may still suffer from under- or misreporting because of
Costa Rica [38] and the formal health care system is well reporting errors, such as missing discharge reports.
developed and largely accessible [23]. The percentage of births Assuming equal distribution of snakebite risk within districts is a
attended by skilled personnel was 98.7% in 2007 [39], highlighting strong assumption. Even if districts are the smallest administrative
a highly developed and utilized formal healthcare system. A survey unit, several spatial processes could still lead to large incidence

PLOS Neglected Tropical Diseases | [Link] 10 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

Table [Link] discriminatory abilities of the estimates model. Another reason why it could be incorrect to compare the
produced by the two smoothing methods (A and B) for the estimated times with the time observed in these studies is that they
two training time lengths (1 and 5 years). are at least 15 years old, and there have been improvements in
ambulance and health care facility accessibility, telecommunica-
tions, and possibly road network since then. It should however be
Test statistic 95% Empirical C.I. further emphasized that our model aim at estimating the ideal
time to reach treatment, and that in reality there could be several
Trainingtime Smoothingmet- mean unaccounted-for causes of longer times, such as problems in
(years) hod AUC Low High
communicating with ambulance stations, temporarily impassable
1 None 0.83 0.79 0.86 roads or unavailable ambulance services, etc. Even though there is
5 None 0.93 0.91 0.95 a discrepancy between the estimated and observed times to
1 A 0.95 0.93 0.96 treatment, we consider that the time-to-treatment estimation
model provides important hints about the location of areas where
5 A 0.96 0.95 0.97
the accessibility to antivenom treatment is more difficult. If the
1 B 0.93 0.91 0.95 minimum time to treatment is estimated to be 2 hours in our
5 B 0.96 0.94 0.97 analysis, there is an imperative for improved accessibility as the
actual time to treatment will probably be longer.
A. The smoothed incidence exceedance probability estimate obtained from the
Bayesian Poisson regression model with spatially correlated random effects and It is however important to remember that healthcare accessi-
fixed effect explanatory variables, described in Methods section. B. The bility cannot be reduced to a purely spatial concern. Logistical
smoothed incidence estimate obtained from the built-in Bayesian spatial issues, such as effective communication with ambulance facilities,
smoothing tool ‘‘Suavizador espacial de tasas’’ in SIGEpi [22]. availability of ambulances, and problems with other forms of
doi:10.1371/[Link].0002009.t002
transportation in the communities need to be taken into account as
well. Furthermore, geographical accessibility is just one dimension
variations within districts. One specific example of where problems of health care access [41]; economic, social and cultural
are likely to arise due to within-district variation is in districts with dimensions need to be also considered, something that is easily
high proportion of urban residents. Even though snakebite missed when doing analyses based on maps only. Research
incidence might be high among the rural population of such gathering empirical evidence on the actual, current time needed to
districts, this might go unnoticed due to the large urban population reach treatment, and determinants of this time, would provide
among which there are few snakebites. By mapping census important information for the identification of vulnerable regions
enumeration units with risk factors favoring snakebite occurrence, and for improving access to snakebite treatment in Costa Rica.
the impact of the above two limitations can be reduced as this risk-
factor based approach is not dependent on the quality of the Interpretation
gathered incidence data or arbitrary district divisions. The estimated minimum times to reach antivenom treatment
The information about the location of the health care system were generally short, compared to the actual times to reach health
must be regarded with caution; previous centrally available facilities after snakebites reported in previous studies in Costa Rica
information about primary health care facilities has been found [19,40]. However, our analysis allowed the identification of some
incorrect in an earlier study [23]. Furthermore, the method of areas where accessibility to antivenom treatment needs to be
locating primary health care facilities (EBAIS) by matching facility improved. The specific strategies to be implemented to accomplish
and community names is not infallible as there was not always a this demand a case by case analysis on a local basis, but a feasible
community with the same name as the EBAIS. However, this was alternative might be the distribution of antivenoms to some
mostly a problem in the urban areas, and thus of smaller EBAIS, the strengthening of the training of health staff in
importance for this study. The data about roads was generally old, antivenom use, and the organization of the work in such a way
and as there has likely been some improvement in road that antivenom is available at all times. There is a risk of over-
availability and quality since the data was gathered, our estimates interpreting the messages transmitted through these maps and
of the time needed to reach treatment might be biased towards forget how sensitive it is to data errors and assumptions of, for
overestimation. example, road speeds. Based on these limitations, we advise that
The estimated transportation times to reach hospital or clinic the maps should be interpreted with care, and that the expert
were much lower than those observed in previous studies of the knowledge of actual conditions provided by health care officials at
time to hospital treatment of snakebites in Costa Rica. In a a local level is also taken into account when making decisions
hospital-based study of all snakebites in 1996 [19], the time to about allocation of treatment resources.
reach hospital was recorded for approximately 70% of the
patients. Of these, 61% reached hospital within 3 hours and Generalizability
20% after more than 5 hours. However, the estimates cannot be The snakebite incidence data available in this study, country-
fully compared with the transportation times observed in this study wide, based on the smallest administrative unit and probably
as they were recorded at hospitals, meaning that a major reliable, are not available in many of the countries where this type
proportion of the patients could have received antivenom of study needs to be conducted, owing to the large underestimation
treatment at clinics and subsequently been transferred to a of snakebite incidence and mortality by hospital statistics
hospital, something that would delay the time to reach hospital [37,42,43]. For data available as small area counts, Bayesian
substantially. Saborio et al. [40] found that among children smoothing techniques have a well-known ability to improve
admitted to the hospital in Limon on the Caribbean coast in interpretability [26], as further demonstrated in this study. Thus,
1985–1995, 50% received medical treatment within 3 hours, using Bayesian smoothing, the interpretability of the gathered data
whereas the mean time was 6.8 hours, indicating very long can be increased so that more accurate estimations of area-level
transportation times for some snakebite victims in this area, parts incidence can be made from sparse data. Leynaud and Reati [21]
of which are also estimated to have long times to treatment by our used a spatial Bayesian smoothing technique available in SIGEpi

PLOS Neglected Tropical Diseases | [Link] 11 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

[22]. We compared this technique with a Bayesian smoothing tion could have led to species disappearance [17]. The method of
technique that also allowed for variation in district risk factor determining high-risk areas by environmental determinants
composition, and enabled estimation of the probability that an requires prior knowledge about the habitat of the snake species
incidence threshold value was exceeded. We found that there were and should ideally be complemented by such field studies of actual
benefits of employing such smoothing techniques to improve the snake distribution. Another important application of GIS in the
interpretability of the raw incidence when the data material was struggle for reducing the impact of snakebite envenoming is the
sparse, whereas there was no significant difference when the data ability to analyze geographical accessibility to treatment, an
material was increased (in this case five years of observation time important factor for the outcome of the bite. Leynaud and Reati
instead of one year). [21] used Euclidean distances to hospitals and roads to analyze
Large-scale approaches to identifying areas in need of access to treatment, a common and readily implemented choice.
antivenom could also benefit from using GIS. Available house- We used a more demanding method aiming at estimating the time
hold-based incidence surveys, hospital and mortality records etc., required for a snakebite victim to reach healthcare that takes into
could be mapped and used to construct geostatistical models account topography, land use, road type and location of
which, based on the spatial variation of snakebite burden and its ambulance services. This theoretically allows for a more detailed
relationship with other spatially varying factors, predict snakebite analysis of the time needed to reach treatment, especially in
burden in areas for which there is no data available. By taking this mountainous study areas.
spatial approach, the sparse data available could be better utilised Our study demonstrates that using GIS it is possible to facilitate
than predicting snakebite burden in non-surveyed areas by rational decision-making on localization of treatment resources
extrapolating information to country or Global Burden of Disease against snakebite by overlaying the risk of snakebite accidents,
Region, as was the method used in the most recent review of estimated using reported data and/or presence of risk factors,
global snakebite burden [1]. There has been an attempt to use transport times to existing hospitals or clinics, and the location of
such methods to map snakebite in West Africa [10], and useful
possible additional facilities to which treatment resources could be
methods have been further developed in studies predicting burden
allocated. GIS is a promising tool for devising cost-effective
of other tropical diseases, such as soil-transmitted helminth
interventions aimed at reducing the public health impact of
infection [44] and malaria [45]. However, in order for these
snakebite envenoming.
approaches to be feasible, there is still a large need for more data
on snakebite burden, especially in sub-Saharan Africa, where a
recent systematical review found only a small number of studies Supporting Information
[46]. Large-scale studies such as those recently conducted to Supporting Information S1 Close-up views of other
estimate snakebite mortality and incidence in India [43] and areas at high risk of snakebite and low antivenom
Bangladesh [37], respectively, provide an important source of data accessibility.
for producing snakebite burden maps, especially if the geograph- (DOC)
ical coordinates of the survey clusters are available.
Mapping the availability of treatment is another challenge to
Acknowledgments
implementing our method on a large scale; antivenom availability
in many areas of low income countries is known to be poor [3], but We gratefully acknowledge the help from researchers at Centro
information on antivenom availability on facility or even country Centroamericano de Población at the University of Costa Rica with
level is to our knowledge not easily available, but requires further accessing and analyzing their databases, Fabián Bonilla for gathering data
data collection. on the distribution of B asper and Silvia Vazquez for compiling data about
snakebites from the hospital discharge reports. We are also grateful for the
GIS not only offers the possibility to improve the interpretation
statistical advice provided by Jonas Björk at Lund University and Gilbert
of incidence data through spatial smoothing, but also to identify Brenes at University of Costa Rica, and the advice on travel time
areas which, on the basis of environmental risk factors, could be estimation provided by Eduardo Pérez Molina at PRODUS at University
expected to have a high snakebite incidence. This approach could of Costa Rica.
be especially useful when incidence data are considered unreliable.
We compared the map of areas with an environment considered as Author Contributions
suitable habitat for Bothrops asper, with a map of specimen
collection locations, and found that these had good congruence, Analyzed the data: EH. Contributed reagents/materials/analysis tools: MS
except for the Nicoya Peninsula, where environmental degrada- AR. Wrote the paper: EH MS JMG AR KM.

References
1. Kasturiratne A, Wickremasinghe AR, de Silva N, Gunawardena NK, 8. Gutiérrez JM, Theakston RD, Warrell DA (2006) Confronting the neglected
Pathmeswaran A, et al. (2008) The global burden of snakebite: a literature problem of snake bite envenoming: the need for a global partnership. PLoS Med
analysis and modelling based on regional estimates of envenoming and deaths. 3: e150.
PLoS Med 5: e218. 9. Chippaux JP (1998) Snake-bites: appraisal of the global situation. Bull World
2. Harrison RA, Hargreaves A, Wagstaff SC, Faragher B, Lalloo DG (2009) Snake Health Organ 76: 515–524.
envenoming: a disease of poverty. PLoS Negl Trop Dis 3: e569. 10. Molesworth AM, Harrison R, Theakston RD, Lalloo DG (2003) Geographic
3. WHO (2007) Rabies and envenomings, A neglected public health issue, Report Information System mapping of snakebite incidence in northern Ghana and
of a Consultative Meeting. Available: [Link] Nigeria using environmental indicators: a preliminary study. Trans R Soc Trop
animal_sera/[Link]. Accessed 14 December 2012. Med Hyg 97: 188–192.
4. Gutiérrez JM, Williams D, Fan HW, Warrell DA (2010) Snakebite envenoming from 11. Hansson E, Cuadra S, Oudin A, de Jong K, Stroh E, et al. (2010) Mapping
a global perspective: Towards an integrated approach. Toxicon 56: 1223–1235. snakebite epidemiology in Nicaragua - pitfalls and possible solutions. PLoS Negl
5. Gutiérrez JM, Fan HW, Silvera CL, Angulo Y (2009) Stability, distribution and Trop Dis 4: e896.
use of antivenoms for snakebite envenomation in Latin America: report of a 12. Rojas G, Bogarı́n G, Gutiérrez JM (1997) Snakebite mortality in Costa Rica.
workshop. Toxicon 53: 625–630. Toxicon 35: 1639–1643.
6. Gutierrez JM (2012) Improving antivenom availability and accessibility: Science, 13. Otero-Patiño R (2009) Epidemiological, clinical and therapeutic aspects of
technology, and [Link] 60: 676–687. Bothrops asper bites. Toxicon 54: 998–1011.
7. Habib AG, Gebi UI, Onyemelukwe GC (2001) Snake bite in Nigeria. Afr J Med 14. Sasa M, Vázquez S (2003) Snakebite envenomation in Costa Rica: a revision of
Sci 30: 171–178. incidence in the decade 1990–2000. Toxicon 41: 19–22.

PLOS Neglected Tropical Diseases | [Link] 12 January 2013 | Volume 7 | Issue 1 | e2009
 Mapping Snakebite Treatment Accessibility Needs

15. Gutiérrez JM (2011) Snakebite poisoning in Latin America and the Caribbean: 30. Caja Costarricense de Seguro Social Gerencia Médica Dirección Proyección de
An integral view from a regional perspective. Boletin De Malariologia Y Salud Servicios de Salud, Área Análisis y Proyección de Servicios de Salud (2009)
Ambiental 51: 1–16. Inventario de Áreas de Salud, Sectores, EBAIS, Sede y Puestos de Visita
16. Sasa M, Wasko DK, Lamar WW (2009) Natural history of the terciopelo Periódica en el Ámbito nacional con corte al 31 de Diciembre del 2009.
Bothrops asper (Serpentes: Viperidae) in Costa Rica. Toxicon 54: 904–922. 31. Cruz Roja Costarricense Ubicación de cómites (2011) Available at [Link]
17. Sasa M, Chaves G, Porras LW (2010) The Costa Rican herpetofauna: [Link]/.Accessed 2011-07-15.
conservation status and future perspectives. In Wilson L, Townsend DJH, 32. Besag J, York J, Mollie A (1991) Bayesian Image-Restoration, with 2
Johnson JD, editors. Conservation of Mesoamerican amphibians and reptiles. Applications in Spatial Statistics. Annals of the Institute of Statistical
Utah: Eagle Mountain Press. pp.510–603. Mathematics 43: 1–20.
18. Fernández P, Gutiérrez JM (2008) Mortality due to snakebite envenomation in 33. Imperial College, Medical Research Council, UK (1996–2007) WinBUGS. 1.4.3
Costa Rica (1993–2006). Toxicon 52: 530–533. edition.
19. Arroyo O, Rojas G, Gutiérrez JM (1999) Envenenamiento por mordedura de 34. IBM Corp. (1989–2011) IBM SPSS Statistics. 20.0.0 edition.
serpiente en Costa Rica en 1996: epidemiologı́a y consideraciones clı́nicas. Acta 35. Brooker S, Hay SI, Bundy DA (2002) Tools from ecology: useful for evaluating
Médica Costarricense 41: 23–29. infection risk models? Trends Parasitol 18: 70–74.
20. Brooker S, Hotez PJ, Bundy DA (2010) The global atlas of helminth infection: 36. Open Source Geospatial Foundation GRASS Development Team (2011)
mapping the way forward in neglected tropical disease control. PLoS Negl Trop Geographic Resources Analysis Support System (GRASS) 6.4.1 edition.
Dis 4: e779. 37. Rahman R, Faiz MA, Selim S, Rahman B, Basher A, et al. (2010) Annual
21. Leynaud GC, Reati GJ (2009) Identifying areas of high risk for ophidism in incidence of snake bite in rural Bangladesh. PLoS Negl Trop Dis 4:e860.
Cordoba, Argentina, using SIGEpi software. Pan American Journal of Public 38. Pan American Health Organisation/World Health Organization (2000)
Health 26: 64–69. Traditional Health Systems in Latin America and the Caribbean: Baseline
22. Área de Análisis de Salud y Sistema de Información (AIS), Organizacion Information.
Panamericana de la Salud (OPS/OMS) (2001–2004): SIGEpi - Sistema de 39. WHO (2011) Costa Rica Country Statistics. Available at [Link]
Información Geográfica en Épidemiologı́a y Salud Pública.[Link] edition. ghodata/. Accessed 2011-06-18.
23. Rosero-Bixby L (2004) Spatial access to health care in Costa Rica and its equity: 40. Saborı́o P, González M, Cambronero M (1998) Snake bite accidents in children
a GIS-based study. Soc Sci Med 58: 1271–1284. in Costa Rica: epidemiology and determination of risk factors in the
24. Unger JP, De Paepe P, Buitron R, Soors W (2008) Costa Rica: achievements of a development of abscess and necrosis. Toxicon 36:3 59–366.
heterodox health policy. Am J Public Health 98: 636–643. 41. Cromley EK, McLafferty SL (2002) Analyzing Access to Health Services. In
25. Gonzalez Block M.Á., González McQuire S (2008) Health Systems of Mexico, Cromley EK, McLafferty SL, editors. GIS and Public Health. New York:
Central America and the Caribbean. In Heggenhougen K, editor. International Guilford Press. pp.233–258.
Encyclopedia of Public Health. Elsevier. pp.256–265. 42. Fox S, Rathuwithana AC, Kasturiratne A, Lalloo DG, de Silva HJ (2006)
26. Wakefield JC, Best N, Waller L (2001) Bayesian approaches to disease mapping. Underestimation of snakebite mortality by hospital statistics in the Monaragala
In Elliott WJP, Best N, Briggs D, editors. Spatial Epidemiology, Methods and District of Sri Lanka. Trans R Soc Trop Med Hyg 100: 693–695.
Applications. Oxford Scholarship Online. pp.103–127. 43. Mohapatra B, Warrell DA, Suraweera W, Bhatia P, Dhingra N, et al. (2011)
27. Best N, Richardson S, Thomson A (2005) A comparison of Bayesian spatial Snakebite mortality in India: a nationally representative mortality survey. PLoS
models for disease mapping. Stat Methods Med Res 14: 35–59. Negl Trop Dis 5: e1018.
28. Herrera W, Gómez L.D (1993) Mapa de Unidades Bióticas de Costa Rica. San 44. Pullan RL, Gething PW, Smith JL, Mwandawiro CS, Sturrock HJ, et al. (2011)
José, Costa Rica: US Fish and Wildlife Service, TNC, INCAF, CBCCR, INBio, Spatial modelling of soil-transmitted helminth infections in Kenya: a disease
Fundación Gómez-Dueñas. control planning tool. PLoS Negl Trop Dis 5: e958.
29. Proyecto Atlas Digital Costa Rica (2008) Atlas Digital de Costa Rica. Cartago: 45. Malaria Atlas Project. Available: [Link]. Accessed 2012-01-25.
Tecnológico de Costa Rica. CD-ROM available after contact with the library of 46. Chippaux J-P (2011) Estimate of the burden of snakebites in sub-Saharan Africa:
this institution. A meta-analytic approach. Toxicon 57: 586–599.

PLOS Neglected Tropical Diseases | [Link] 13 January 2013 | Volume 7 | Issue 1 | e2009