Abstract
Background Post-exposure prophylaxis (PEP) is highly effective at preventing human rabies deaths, however access to PEP is limited in many rabies endemic countries. The 2018 decision by Gavi to add human rabies vaccine to its investment portfolio should expand PEP availability and reduce rabies deaths. We explore how geographic access to PEP impacts the rabies burden in Madagascar and the potential benefits of improved provisioning.
Methodology & Principal Findings We use travel times to the closest clinic providing PEP (N=31) as a proxy for access. We find that travel times strongly predict reported bite incidence across the country. Using resulting estimates in an adapted decision tree framework we extrapolate rabies deaths and reporting and find that geographic access to PEP shapes burden sub-nationally. We estimate 960 human rabies deaths annually (95% Prediction Intervals (PI):790 - 1120), with PEP averting an additional 800 deaths (95% PI: 800 (95% PI: 640 - 970) each year. Under these assumptions, we find that expanding PEP to one clinic per district could reduce deaths by 19%, but even with all major health centers provisioning PEP (1733 additional clinics), we still expect substantial rabies mortality. Our quantitative estimates are most sensitive to assumptions of underlying rabies exposure incidence, but qualitative patterns of the impacts of travel times and expanded PEP access are robust.
Conclusions & Significance PEP is effective at preventing rabies deaths, and in the absence of strong surveillance, targeting underserved populations may be the most equitable way to provision PEP. Our framework could be used to guide PEP expansion and improve targeting of interventions in similar endemic settings where PEP access is geographically restricted. While better PEP access should save many lives, improved outreach and surveillance is needed and if rolled out with Gavi investment could catalyze progress towards achieving zero rabies deaths.
Author Summary Canine rabies causes an estimated 60,000 deaths each year across the world, primarily in low- and middle-income countries where people have limited access to both human vaccines (post-exposure prophylaxis or PEP) and dog rabies vaccines. Given that we have the tools to prevent rabies deaths, a global target has been set to eliminate deaths due to canine rabies by 2030, and recently, Gavi, a multilateral organization that aims to improve access to vaccines in the poorest countries, added human rabies vaccine to it’s portfolio. In this study, we estimated reported bite incidence in relation to travel times to clinics provisioning PEP, and extrapolate human rabies deaths in Madagascar. We find that PEP currently averts around 800 deaths each year, but that the burden remains high (1000 deaths/ year), particularly in remote, hard-to-reach areas. We show that expanding PEP availability to more clinics could significantly reduce rabies deaths in Madagascar, but our results suggest that expansion alone will not eliminate deaths. Combining PEP expansion with outreach, surveillance, and mass dog vaccination programs will be necessary to move Madagascar, and other Low- and Middle-Income countries, forward on the path to rabies elimination.
Introduction
Inequities in access to care are a major driver of disease burden globally [1]. Often, the populations at greatest risk of a given disease are the most underserved [2]. Delivering interventions to these groups is challenging due to financial and infrastructural limitations, and requires careful consideration of how best to allocate limited resources [3].
Canine rabies is estimated to cause approximately 60,000 human deaths annually [4]. Mass vaccination of domestic dogs has been demonstrated to be a highly effective way to control the disease in both animals and humans. While dog vaccination can interrput transmission in the reservoir, human deaths can also be prevented through prompt administration of post-exposure prophylactic vaccines (PEP) following a bite by a rabid animal [5]. However, access to human vaccine is limited in many countries where canine rabies is endemic [6–8], and within countries these deaths are often concentrated in rural, underserved communities [9].
In 2015, a global framework to eliminate deaths due to canine rabies by 2030 (‘Zero by 30’) through a combination of PEP provisioning and dog vaccination was established by the World Health Organization (WHO) and partners [10]. Furthermore, in 2018, Gavi, the Vaccine Alliance, added human rabies vaccines to their proposed investment portfolio [11]. From 2021, Gavi-eligible countries should be able to apply for support to expand access to these vaccines, with potential to greatly reduce deaths due to rabies.
A primary challenge in expanding access effectively is the lack of data on rabies exposures and deaths in humans and incidence in animals in most rabies-endemic countries [12]. Deaths due to rabies are often severely underreported, with many people dying outside of the health system, often in remote and marginalized communities [13]. Instead of directly measuring rabies deaths, the majority of rabies burden studies use bite patient data on reported bites at clinics provisioning PEP and a decision tree framework to extrapolate deaths, assuming that overall reported bite incidence (i.e. both bites due to non-rabid and rabid animals) is proportional to rabies incidence (i.e. the more bites reported in a location, the higher the incidence of rabies exposures), and that reporting to clinics for PEP is uniform across space [8,14,15]. If applied subnationally, these assumptions would likely underestimate rabies deaths in places with poor access to PEP and may overestimate rabies deaths in places with better access to PEP.
In Madagascar, the Institut Pasteur de Madagascar (IPM) provides PEP to 30 Ministry of Health clinics, in addition to it’s own vaccine clinic, where PEP is available at no direct cost to patients [15]. Other than at these 31 anti-rabies medical centers (ARMC), PEP is not available at any other public clinics or through the private sector. In addition, there is limited control of rabies in dog populations and the disease is endemic throughout the country [16,17]. Due to the spatially restricted nature of PEP provisioning and lack of direct costs for PEP, geographic access is likely to be a major driver of disease burden within the country. Previously, we estimated the burden of rabies in Madagascar nationally using data from a single district to extrapolate to the country, and did not account for spatial variation in access [15]. Here, we provide revised estimates of human rabies deaths by incorporating the impact of access to PEP at the sub-national level on preventing human rabies deaths and explore the potential impact of expanding provisioning of human rabies vaccines on further reducing these deaths. This framework may usefully apply to other countries where PEP availability is currently geographically restricted in considering how to most effectively and equitably provision these life-saving vaccines.
Methods
Estimating geographic access to ARMC
To estimate mean and population weighted travel times to the nearest clinic at the scale of administrative units, we used two rasters: 1) the friction surface from the Malaria Atlas Project [18] at an ∼1 km^2 scale (Fig S1.1A) and 2) the population estimates from the the 2015 UN adjusted population projections from World Pop ([19], originally at an ∼100m^2 resolution, Fig S1.1B), which we aggregated to the friction surface.
From GPS locations of the 31 ARMCs we estimated the travel time to the nearest clinic at an approximately 1 x 1 km scale as described in [18]. We then extracted the mean and population-weighted mean travel times for the district and commune (the administrative unit below the district), and euclidean distance, i.e. the minimum distance from the administrative unit centroid to any ARMC. We used shapefiles from the UN Office for the Coordination of Humanitarian Affairs for the administrative boundaries (as of October 31, 2018). We compared travel times and distance estimates to driving times collected by IPM during field missions and patient reported travel times from a subset of Moramanga ARMC patients (see Fig S1.2 for raw data).
Estimating bite incidence
We used two datasets on bite patients reporting to ARMC:
A national database of individual bite patient forms submitted to IPM from ARMC across the country between 2014 - 2017. These forms were submitted with frequencies ranging from monthly to annually, included the patient reporting date and were resolved to the district level (patient residence).
33 months of data (between October 2016 and June 2019) on patients reporting to the Moramanga District ARMC resolved to commune level.
For the national data, some clinics did not submit any data, or had substantial periods (months to a whole year), with no submitted data. To correct for this, we exclude periods of 15 consecutive days with zero submitted records (see Supplementary Appendix, section S2). For each clinic we divided the total number of bites reported in a given year by the estimated proportion of forms which were not submitted (under-submission). Due to yearly variation in submissions, we took the average of annual bite incidence estimates aggregated to district level. We validated this approach by comparing estimated vial demand given the total reported bites corrected for under-submission to vials provisioned to clinics for 2014-2017 (see Supplementary Appendix, section S2). At both the commune and district level, we assigned clinic catchments by determining which were closest in terms of travel times for the majority of the population within the administrative unit. For national data, we excluded any districts in a catchment of a clinic which submitted less than 10 forms and any years for which we estimated less than 25% of forms were submitted.
Modeling reported bite incidence
We modeled the number of reported bites as a function of travel time (T) using a Poisson regression:
where yi is the average number of bites reported to a clinic annually and μi the expected number of bite patients presenting at the ARMC as a function of travel time (Ti) and human population size (Pi) (an offset which scales the incidence to the expected number of bites) for a given source location (district or commune). We fit this model to both the national data (district level) and the Moramanga data (commune level). To more directly compare estimates between datasets, we also modeled the national data with a latent commune-level travel time covariate (Tj):
As travel times are correlated with population size (Fig S3.1), we also compared how well bites were predicted by population size alone, and in combination with travel times. For the models with population size, we removed the offset and used either population size alone 
or population size and travel times 
as predictors.
For the models fit to the national data, we also modeled variation between clinics with a catchment random effect: B0,k ∼ norm(μ, σ0)), where μ is the mean and σ0 is standard deviation and B0,k is the catchment level intercept.
We tested whether the catchment random effect captured overdispersion in the data (i.e. variance > mean – the expectation given a Poisson distribution) rather than any catchment specific effects by extending these models with an overdispersion parameter, σe [20]:
where 
is the sum of the all parameters for a given model. We fit all models in a Bayesian regression framework via MCMC using the R package ‘rjags’ [21]. We used model estimates to generate fitted and out-of-fit predictions, and examined the sensitivity of estimates to adjustments for under-submission of forms (Supplementary Appendix, section S3).
Modeling human rabies deaths
We estimate rabies deaths as a function of the number of bites predicted by our model and estimates of endemic rabies exposure incidence using an adapted decision tree framework (Fig 1). To model uncertainty in parameter estimates we used triangular distributions, as with previous studies [8,22], for two key parameters: Ei, the annual exposure incidence per administrative unit, and prabid, the proportion of reported bites that are rabies exposures.
For a given administrative unit i, human deaths due to rabies (Di) are calculated from model predicted reported bites (Bi). To get Ri, the number of reported bites that were rabies exposures, we multiply Bi by prabid, the proportion of reported bites that are rabies exposures. prabid is constrained such that reported rabies exposures cannot exceed total rabies exposures (Ei) and reporting of patients to clinics cannot exceed a maximum (pmax). Ri is subtracted from Ei to get the number of unreported bites (Ui) and then multiplied by the probability of death given a rabies exposure (pdeath) to get deaths due to rabies (Di). Similarly, deaths averted by PEP, Ai, are estimated by multiplying Ri by pdeath, i.e. those who would have died given exposure, but instead received PEP. Both Ei and prabid are drawn from a triangular distribution. Parameter values and sources are in Table 1.
For Ei, we center the distribution at the lower end of our estimated exposure incidence from the Moramanga District (42 exposures/100,000 persons), with a range applied assuming 1% rabies incidence in dogs (estimated across a range of human-to-dog ratios between 5 - 25) and that on average a rabid dog exposes 0.39 persons [4] (see Fig S4.1). As there is little data on dog population size and human exposure incidence in Madagascar[16,23], the range we used encompasses both observed human-to-dog ratios across Africa [14,24] and recent subnational estimates from Madagascar [25], and generates similar exposure incidences as observed previously across Africa [26,27].
Parameters used in the decision tree to estimate human rabies deaths at the administrative level.
For prabid, we use a range between 0.2 - 0.6 estimated from a study of bite patients in the Moramanga District [15]. So that rabid reported bites cannot exceed the total expected number of rabies exposures or a maximum reporting (even with minimal travel times, people may not report for PEP for other reasons), we constrain prabid :
where ρmax is the maximum reporting, estimated from the Moramanga ARMC data for the commune of Moramanga Ville, the closest commune to the ARMC (average of 3.12 minutes travel time to the clinic), and where we find that approximately 2% of rabies exposures go unreported [15].
We assume that all rabies exposed patients who report to an ARMC receive and complete PEP, and PEP is completely effective at preventing rabies. Fig 1 describes the decision tree, the key inputs, and outputs (Ai, deaths averted by PEP, and Di, deaths due to rabies). Table 1 list all inputs and their sources.
Estimating the impact of expanding PEP provisioning
We developed a framework to rank clinics by how much their PEP provision improves access for underserved communities, estimating incremental reductions in burden and increases in vaccine demand. Specifically, we aggregated our model-predicted estimates of annual bites to the clinic level. As multiple clinics may serve a single district or commune, we allocated bites to clinics according to the proportion of the population in each administrative unit which were closest. For each clinic, we simulate throughput by randomly assigning patient presentation dates, and then assume perfect compliance (i.e. patients report for all doses) to generate subsequent vaccination dates. We use these dates to estimate vial usage given routine vial sharing practices in Madagascar [15], but assuming adoption of the WHO-recommended abridged intradermal regimen (2 x 0.1 ml injections on days 0, 3, and 7 [29]). For both burden and vial estimates, we take the mean of 1000 simulations as each clinic is added.
We simulate expansion first to each district (N = 114) and then to each commune in the country for all communes with a clinic. We select the CSB II (Centre de Sante Niveau II, major health centers with immunization capacity) in the highest density grid cell of the administrative unit as candidates for expansion. For the 85 communes without a CSB II, we chose the CSB I (secondary health posts) in the highest density grid cell. 94 communes lacked any CSB I or II. Finally, we explore a scenario where all additional CSB II (totalling 1733) provision PEP.
We tested three metrics for ranking additional clinics: 1) The proportion of people living >3 hours from an existing ARMC for which travel times were reduced; 2) This proportion (1), weighted by the magnitude of the change in travel times and 3) The mean reduction in travel times for people living >3 hours from an existing ARMC. We simulated expansion of ARMC to each district using these three metrics, and chose the metric which decreased burden the most compared to simulations (N = 10) where clinics were added randomly to districts for the full expansion of PEP. For the full simulation of expanded access, once clinics reduced travel times for less than 0.01% of the population (< 2400 living greater than x hrs away, starting with x = 3 hrs), we reduced the travel time threshold by 25%.
Sensitivity analysis
To test the effect of our model assumptions on estimates of rabies burden and vial demand, we did a univariate sensitivity analysis of both parameters from the models of bite incidence and the decision tree (see Table S6.1 & S6.2 for parameter ranges used). We also examined how systematic variation in rabies incidence with human population size affected burden estimates. Specifically, if human-to-dog ratios are positively correlated with human populations (i.e. dog ownership/populations are higher in more populated, urban areas), we might expect higher rabies exposure incidence as population size increases. Alternatively, if human-to-dog ratios inversely correlate with population size (i.e. dog ownership is more common in less populated, rural areas), we might expect exposure incidence to scale negatively with population size. We use scaling factors to scale incidence either positively or negatively with observed population sizes at the district and commune levels, while constraining them to the range of exposure incidence used in the main analyses (15.6 - 76 exposures per 100,000 persons, Fig S4.2) and simulated baseline burden, as well as expanded PEP access.
Data and analyses
All analyses were done in R version 4.0.2 (2020-06-22) [30] and using the packages listed in the supplementary references (Supplementary appendix, section S7). All processed data, code, and outputs are archived on Zenodo (http://doi.org/10.5281/zenodo.4064312 and https://doi.org/10.5281/zenodo.4064304), and maintained at https://github.com/mrajeev08/MadaAccess. The raw bite patient data at the national level are maintained in two secure REDCap (project-redcap.org) databases, one for IPM and another for all peripheral ARMC. These databases were last queried on September 19, 2019 for these analyses. The IPM GIS unit provided the data on geolocated clinics across the country. Anonymized raw bite patient data and full data on geolocated clinics are available from IPM following insitutional data transfer protocols. Anonymized raw data collected from the Moramanga District were retrieved from the Wise Monkey Portal (wisemonkeyfoundation.org) on the same date and are shared in the archived repository.
Ethics statement
Data collection from the Moramanga District was approved by the Princeton University IRB (7801) and the Madagascar Ministry of Public Health Ethics Committee (105-MSANP/CE). Oral informed consent was obtained from all interviewed participants. Data collected from bite patients at the national level are maintained jointly by the Ministry of Health and IPM as a routine part of PEP provisioning.
Results
Estimating access to ARMC
Estimates of travel times to the closest ARMC varied greatly (Fig 2A), with approximately 36% of the population estimated to live over 3 hrs from a clinic (Fig 2B). We found that travel time estimates from the friction surface underestimated both driving times across the country and patient-reported travel times to the Moramanga ARMC (Fig 2C), but were correlated with ground truthed driving and patient-reported times (Fig 3C, Fig S1.4). Patient-reported travel times were also significantly more variable than estimates from the friction surface (Fig S1.2). The friction surface assumes that the fastest available mode of transport is used across each grid cell (i.e. the presence of a road indicates that all travel through that grid cell is by vehicle). However, Moramanga ARMC patients reported using multiple modes of transport, with some individuals walking days to a clinic (Fig S1.3). Travel times weighted by population at the grid cell level were a better predictor than unweighted travel times or distance (Table S1.1). Therefore we use population-weighted travel time as a proxy for access at the commune/district level in subsequent analyses.
(A) Estimated at an ∼ 1 km2 scale. (B) Distribution of the population across travel times. (C) Correlation between ground-truthed travel times (mean of patient reported travel times to the Moramanga ARMC at the commune level and reported driving times between GPS points) and friction surface travel time estimates. The vertical lines show the 95% quantiles for reported travel times and the point size shows the number of observations for each commune. The best fit lines (red and grey) and 95% confidence interval (shading) from a linear model where observed travel times are predicted by estimated travel times for each data source are also shown. The dashed black line is the 1:1 line.
(A) at the district level for the national data and (B) commune level for the Moramanga data: circles with a black outline represent the total number of patients reporting to each ARMC for which we have data. Color corresponds to the clinic catchment. Circles with a white outline are the total number of bites reported for that administrative unit (plotted as the centroid). Lines show which ARMC those patients reported to, with the line width proportional to number of patients from that district reporting to the ARMC; flows of less than 5 patients were excluded. Out-of-catchment reporting is indicated where points and line colors are mismatched. For panel (A) districts in catchments excluded due to lack of forms submitted by the clinic are colored in grey. For (B) the inset of Madagascar shows the location of the enlarged area plotted, which shows the district of Moramanga (outlined in black), all communes included in it’s catchment (red polygons), and other communes where bites were reported to colored by their catchment (C) The estimated average annual bite incidence from the national and Moramanga data plotted at the district scale (both datasets) and at the (D) commune scale (Moramanga dataset). Colors correspond to the clinic catchment, shape indicates the dataset, and the size of the point indicates the number of observations (i.e. the number of years for which data was available for the national data; note for Moramanga 33 months of data were used). The point lines indicate the range of estimated bite incidence for each district.
Estimating bite incidence
Most patients from each district reported to their closest ARMC by the weighted travel time metric (Fig 3). Accordingly, we assigned catchments based on which clinic was the closest for the majority of the population. While there are discrepancies between commune and district catchment assignments (Fig S2.1A), over 75% of the population in a given district or commune were closest to a single clinic (Fig S2.1C). We excluded any clinics which submitted less than 10 forms (excluded 11 catchments, Fig 3A grey polygons) and corrected for periods where clinics did not submit any forms (see Supplementary Appendix Section S2). After additionally excluding any year with less than 25% of forms submitted, our final dataset consisted of estimates of average bite incidence for 83 of 114 districts (Fig 3C), and 58 communes within the catchment of the Moramanga District (Fig 3D).
Bite incidence estimates generally increased with decreasing weighted travel times at both scales, although there was considerable variation between catchments for the magnitude of this relationship (Fig 3C and D). For the national data, there were two outliers, Toamasina II (the sub-urban district surrounding the city of Toamasina) and Soanierana Ivongo, with higher bite incidence when compared to other districts with similar travel times. While the estimates from the Moramanga data showed higher reported incidence at low travel times at the commune level compared to the district estimates, when aggregated to the district, bite incidence estimates fell within the ranges observed from the national dataset.
Modeling reported bite incidence
Travel times were a strong and consistent predictor of reported bite incidence in both datasets and across scales (Fig 4). Population size alone was the poorest fit to the data as estimated by DIC (Table S3.1), and models with population size as an additional covariate did not generate realistic predictions to the observed data or when used to predict out of fit (Figs S3.2 and S3.3).
(A) The estimated relationship between travel time in hours (x-axis) and mean annual reported bite incidence (y-axis). The lines are the mean estimates and the envelopes are the 95% prediction intervals generated by drawing 1000 independent samples from the parameter posterior distributions for three candidate models: model with travel times at the 1) commune- and 2) district-level fitted to the national data with an overdispersion parameter (σe) and 3) travel times at the commune level fitted to the Moramanga data with a fixed intercept and unadjusted for overdispersion. The points show the data: National data (circles) at the district level used to fit the District and Commune models, and Moramanga data (triangles) at the commune level used to fit the Moramanga model. (B) The posterior distribution of parameters from the respective models for the model intercept, travel time effect, and for overdispersion (national data only).
For the national data, including a catchment random effect improved predictions (Fig S3.2 & Fig S3.3). However, after accounting for overdispersion, catchment effects were not clearly identifiable (Table S3.1) and the models resulted in similar predictions (Fig S3.6 & S3.7), indicating that catchment effects could not be differentiated from random variation in the data. Similarly, while the commune model fit to the Moramanga data generated stronger travel time effects (Fig 4B), after accounting for data overdispersion, the posterior estimates of the parameters overlapped for the commune and district models fit to the national data (Fig S3.4), and the model estimates were in general less robust to overdispersion than for the national data, particularly at low travel times (Fig S3.5).
As the predictions from the model fit to the Moramanga data without accounting for overdispersion fall within the prediction intervals for the models fit to the national data (Fig 4A), for subsequent predictions, we used the parameter estimates from models fit to the national data, which encompass the range of travel time effects observed in our datasets. Moreover, our out-of-fit predictions to the data across scales suggest that the commune model is able to capture travel time impacts at the commune level (Fig S3.3), therefore we use both the district and commune model to disaggregate burden to the finest scale possible. Finally, we examined the sensitivity of models to how we corrected for underreporting of data, and found that parameter estimates of travel time impacts were similar across models and performed similarly in prediction (Fig S3.8 and Fig S3.9).
Estimating human rabies deaths
Overall, we estimate close to 1000 rabies deaths (95% PI: 800 - 1100) annually in Madagascar. Our estimates vary only slightly depending on the scale of the model (Table 2), but disaggregating deaths to the commune level shows considerable variation in predicted burden within districts (Fig 5A). Under the current system of 31 ARMCs in Madagascar, we estimate that use of PEP prevents approximately 800 (95% PI: 600 - 1000) deaths due to rabies each year. In general, the incidence of rabies deaths increases with travel times to clinics, and there is significant sub-national variation when deaths are modeled at the district and commune scale, with the least accessible communities having most deaths (Fig 5B & C).
Model predictions of average annual reported bite incidence, total deaths, and deaths averted at the national level for the two models (commune level and district level models with travel time predictor and an overdispersion parameter); 95% prediction interval in parentheses.
(A) for each district (y-axis) in Madagascar. Diamonds show the predicted incidence for the district model and squares show predicted incidence for the commune model fit to the National data for all communes in a district. Points are colored and districts ordered by travel times. The vertical lines show the average national incidence of human rabies deaths for the commune (grey) and district (black) models. Incidence mapped to the (B) commune- and (C) district-level from the respective models; grey X’s show locations of current ARMC.
Estimating the impact of PEP provisioning
For a subset of additional clinics (n = 83, up to one per district), we compared three methods of ranking for expanding PEP provisioning. We found that targeted expansion of PEP to clinics based on the proportion of the population they reduced travel times for resulted in fewest deaths (Fig S5.1). Therefore, we used this approach for simulating expansion of PEP to a larger set of clinics (N = 1733). Here we report results from the commune model, as estimates were consistent across models (Fig 6 and Supplementary appendix, section S5).
(A) Decrease in deaths due to rabies, (B) increase in total # of vials as additional ARMC are added at the national level, and (C) increase in vials needed per death averted based on the two models of reported bites. Lines are the mean of 1000 simulations with envelopes representing 95% prediction intervals. The points show when all additional csb (N = 1733) clinics have been added).
We estimated that strategic PEP expansion to these additional 83 clinics (1 per district) reduced rabies deaths by 19% (95% PI: 14 - 23%) (Fig 6A). With one clinic per commune (where available, N = 1696), we see a further reduction of 38% (95% PI: 30 - 46%). However, reductions in burden saturate as more clinics are added (Fig S5.2). Even when all CSB II provisioning PEP, our model still predicts 600 (95% PI: 400 - 800) deaths per annum, and average reporting of rabies exposures remains approximately 66% (95% PI: 32 - 77%) (Fig S5.5); as more clinics are added, reported bite incidence saturates (Fig S5.4), and patients shift which clinic they report to (S5.7 & S5.8).
Vial demand also outpaces reductions in burden (Fig 6B), with more vials needed per death averted (Fig 6C). Our model predicts an increase from 33500 vials (95% PI: 22900 - 49400) per annum under current provisioning but with the abridged intradermal regimen (i.e. visits on days 0, 3, 7), to ∼56900 vials (95% PI: 40200 - 77800) with 250 clinics providing PEP, and ∼86400 vials (95% PI: 61600 - 118000) if all CSB IIs provision PEP. In these scenarios, clinic catchment populations and throughput decrease, with clinics seeing fewer patients per day (S5.6).
Sensitivity Analyses
To quantify which parameters contribute to the uncertainty in our estimates, we performed univariate sensitivity analyses, varying parameters in both the bite incidence model (Table S6.1) and the decision tree (Table S6.2). These analyses show that assumptions of the underlying rabies exposure incidence (Ei) contribute most uncertainty to burden estimates (Fig S6.1) and impacts of PEP access (Fig S6.2). Uncertainty in bite model parameters contribute less to uncertainty in estimates of burden. For the estimates of vial demand, uncertainty around the model intercept (i.e. the baseline reported bite incidence) has most impact, rather than the travel time effect or the overdispersion parameter (Fig S6.3). Finally, scaling of incidence with population size (Fig S4.2) modulates the impact of travel times on deaths, with positive scaling of rabies incidence with population size (i.e. more rabies in more populated places) dampening and negative scaling exacerbating the relationship between access and deaths (Fig S6.4A). However, the impact of adding clinics remains broadly the same (Fig S6.4B). Overall, while uncertainty in the underlying rabies exposure incidence results in considerable variation in our burden estimates, the projected impact of travel times and of access to PEP are qualitatively similar across parameter assumptions.
Discussion
Main findings
We find that the burden of rabies in Madagascar is likely concentrated in areas with poor PEP access. We estimate that current PEP provisioning (at 31 clinics) averts 45%% of deaths that would otherwise occur, and that expanding PEP access should reduce mortality, with provisioning in one clinic per district (N = 83), or per commune (N = 1733), expected to reduce mortality by 33% and 16%, respectively. However, improved PEP provisioning is unlikely to eliminate rabies deaths; with over 600 deaths expected even with PEP at all CSB II (N = 1733). This is partly because travel times remain high (> 2 hrs as estimated by the friction surface for over 10% of the population, Fig S5.4) even after expanding PEP to all major health centers. With reduced travel times, over 20% of exposures will still not seek PEP (Fig S5.5), resulting in ∼1.65 rabies deaths per 100,000 people. PEP is expected to remain cost-effective as provisioning expands, to a maximum of 450 USD per death averted (assuming 5 USD per vial), similar to other estimates [4]. While our quantitative predictions depend on assumptions of underlying rabies exposure incidence, the qualitative patterns regarding travel time impacts remain robust and are useful in identifying strategies for provisioning PEP.
Limitations
Data limitations introduced bias and uncertainty to our estimates. For example, travel times from the MAP friction surface underestimated patient-reported travel times, with discrepancies between assigned transport speeds (from the Open Street Map user community, or country cluster data [18]) and realities of local travel. In Madagascar, the presence of paved roads does not necessarily reflect their quality or the modes of transport used. Moramanga ARMC patients reported various transport methods and highly variable travel times even within a single commune. While patient-reported travel times may lack precision from recall and estimation error, they likely better reflect lived experience; validated travel times [31] could improve estimates of spatial health inequities. Similarly, modeled estimates of population distribution [19] also introduce uncertainty. Our analysis of data from the Moramanga District indicate that variation at the sub-district level is high and impacts health seeking behavior. However, we lacked fine-scale data from other catchments for comparison. Additionally, we had to correct for underreporting and incomplete data; strengthening surveillance and routine data collection should improve understanding of health seeking behavior and access, and support monitoring and evaluation of PEP provisioning.
We assumed geographic access to PEP was the primary driver of health-seeking behavior, but socioeconomic status, education and awareness about rabies [27,32–34] all play a role. In Madagascar, where PEP is free-of-charge, the main cost to patients is transport and time lost. More remote communities tend to be of lower socioeconomic and educational status [2], so travel time may be a proxy for these correlated variables. We also assume that all ARMC reliably provision PEP, but a 2019 KAP survey reported clinics experiencing stock-outs [25]. Most ARMC charge fees (from 0.50 - 3.00 USD for consultations, wound treatment, etc [25]) which may also act as barriers. Significant overdispersion in the data that cannot be explained by travel times suggests that clinic-level variation (e.g. vaccine availability and charges) and regional differences (e.g. dog populations, outbreaks, awareness) further influence health-seeking behavior and vaccine demand.
Our burden estimates were most sensitive to assumptions about rabies exposure incidence, drawn from studies in the Moramanga District [15] and elsewhere [4]. As incidence of rabies exposures varies over time and space [35,36], we incorporated uncertainty into our estimates, but we did not find qualitative differences in the effects of travel times on rabies deaths. Our simplifying assumptions, regarding patient compliance, which is generally high in Madagascar [15], and on complete efficacy of PEP, are unlikely to greatly influence our burden estimates [28]. Likewise we do not account for differential risk for severely exposed patients not receiving Immunoglobulins (RIG), which is only available at IPM in Antananarivo, but recent studies show that even in the absence of RIG, PEP is extremely effective [4]. Although our estimates could be improved with better data on rabies incidence, health-seeking behavior, and PEP provisioning, predicting PEP impacts will remain challenging given the complex interactions between socioeconomic factors, access to and quality of care and human behavior, as illustrated by the complex case studies in Box 1.
Box 1: case studies of of health seeking behavior for PEP in Madagascar
Anosibe An’ala District (population ∼ 100,265), south of Moramanga, has moderate incidence of bite patients (∼ 54/100,000 persons) even though travel times often exceed 24 hours. While a road connects the main Anosibe An’ala commune to the Moramanga ARMC, it is only passable by large trucks during much of the rainy season, with speeds usually < 10km per hour. Over 9% of patients from Anosibe An’ala had been in close proximity or touched a person that died from rabies (four suspect human rabies deaths of patients who did not receive any PEP), whilst of patients with Category II and III exposures that were interviewed, 11/19 (58%) were bitten by probable rabid dogs. Given the high travel times (although underestimated by the friction surface) and incidence of reported rabies exposures and deaths, we predict a large but unobserved rabies burden in this remote community (∼6.02 deaths per year) and we ranked an ARMC provisioning PEP in Anosibe An’ala 28th for travel time reductions. Other remote communities likely experience similar high and unrecognized burden, but improved surveillance is necessary to identify such areas. Notably, bite patients in this district demonstrate willingness to travel for free PEP (in some cases walking 3 days to a clinic) with awareness of rabies risk. Community outreach and active surveillance in other remote areas could also greatly improve people’s awareness of risk and health seeking behavior.
Recently, a middle aged taxi driver died of rabies in suburban Antananarivo. The day after being bitten by an unknown dog, he reported to a clinic that referred him to the ARMC at IPM, approximately one hour’s drive from his home. His family urged him to get PEP, but he did not believe his risk was high and decided not to seek further care. He developed symptoms two weeks later and was confirmed as a rabies death by the National Rabies Reference Laboratory. Despite prompt reporting and referral, and socioeconomic indicators suggesting a high care-seeking probability, this person did not receive PEP. His story highlights the need for sensitization about rabies and how PEP at peripheral clinics (even in areas with reasonable access) could prevent additional deaths, but also how PEP alone is unlikely to prevent all rabies deaths.
Broader context
Recent studies have estimated access to health-seeking behavior and PEP completion and adherence, but not directly linked these estimates to burden [7,37]. Our approach for incorporating access to vaccines (echoing [38–42]) into burden estimation methods could guide provisioning of PEP to maximize impacts. This approach will have most value in settings with limited PEP access and poor health seeking, but will be less valuable where rabies exposures make up a small fraction of patients reporting for PEP e.g. [43,44]. Our revised estimate of rabies deaths in Madagascar using this approach was higher than previous [15], which assumed uniform reporting of 85%, but remained within the range of other empirical and modeling studies from low-income countries [26,27,45–47].
Our estimates of vial demand depend on use of the new abridged intradermal regimen [29], which has been adopted by the Ministry of Health in Madagascar. However, most ARMC staff were not aware of WHO classifications of exposure categories, and vaccination of Category I exposures (those not requiring PEP) remains common practice, comprising 20% of vial demand in Moromanga [15]. We predict that as clinics are expanded, throughput (daily patients reporting to a clinic) will decrease. This may make provisioning PEP more challenging and vial demand less predictable, leading to stock outs or wastage. Decentralized provisioning mechanisms, for example adopting routine childhood vaccine supply chains, or novel vaccine delivery methods such as drones [48], may mitigate these challenges. When nerve tissue vaccines were used in Madagascar, clinics requested vaccine upon demand and PEP access was more widespread, but provisioning the more expensive cell culture vaccines to all clinics became too costly [16]. Widespread vaccine provisioning is therefore feasible given Madagascar’s health infrastructure, if cost barriers are removed.
Gavi investment could greatly reduce the access and cost barriers to PEP [6,7,28,49]. Currently, each clinic in Madagascar serves an average catchment of 780,000 persons. Latin American countries, where significant progress has been made towards elimination, aim for one PEP clinic per 100,000 persons. In Madagascar this would require around 212 additional ARMC. We predict that Gavi investment would be highly cost-effective, greatly reducing deaths by expanding PEP supply to underserved areas.
However, our results suggest that PEP expansion alone cannot prevent the majority of rabies deaths, and even given maximal access, achieving ‘the last mile’, preventing deaths in the most remote populations, will require disproportionate resources [50]. To achieve ‘Zero by 30’, mass dog vaccination will be key to interrupting transmission, and eliminating deaths. Integrated Bite Case Management (IBCM) uses bite patient risk assessments to determine rabies exposure status, guide PEP administration, and trigger investigations of rabid animals, potentially identifying other exposed persons [15,51,52]. IBCM is one way to manage PEP effectively [43] and as it relies on exposed persons reporting to clinics, expanding PEP access could strengthen this surveillance framework. These same issues of access, however apply to both dog vaccination and surveillance, and understanding spatial heterogeneities will be critical to determining how control and prevention interventions can be best implemented [53,54].
Conclusion
Our study suggests that rabies deaths in Madagascar disproportionately occur in communities with the poorest access to PEP and that expanding PEP access should reduce deaths. Without data on rabies incidence and exposure risk, targeting PEP expansion to underserved areas is a strategic way to reduce rabies burden and provide equitable access, for example, by expanding provisioning to clinics serving populations that target an evidence-based travel time threshold or catchment size. Implementing outreach programs to raise awareness should further increase the efficacy of PEP expansion by improving care seeking. Better surveillance is also needed to understand the geographical distribution of rabies exposures and identify populations most at risk. Gavi investment could support countries to more equitably provision PEP and overcome barriers to access ([9], see Box 1 for case studies), but as PEP alone cannot prevent all rabies deaths, investment should be used to catalyse mass dog vaccination to interrupt transmission, and eventually eliminate rabies deaths.
Data Availability
All processed data, code, and outputs are archived (here)[zenodo] and maintained at github.com/mrajeev08/MadaAccess. Anonymized raw bite patient data and full data on geolocated clinics are available from IPM following insitutional data transfer protocols.
Acknowledgements
We thank all the clinicians and staff at the ARMC across the country. We are grateful to IPM and the Ministry of Public Health who collect and maintain data on PEP provisioning. In particular, we thank the GIS unit for assistance with spatial data, Michael Luciano Tantely for sharing driving time data, and Claire Leblanc, Rila Ratovoson, and Daouda Kassie for sharing results of their work in the Moramanga District. In addition, we thank Jean Hyacinthe Randrianarisoa, Ranaivoarimanana, Fierenantsoa Randriamahatana, Esther Noiarisaona, Cara Brook, Amy Winter, Christian Ranaivoson, John Friar, and Amy Wesolowski for assistance. This work was funded by grants from the Center for Health and Wellbeing and the Department of Ecology and Evolutionary Biology at Princeton University to CJEM and MR. MR is supported by an NSF Graduate Research Fellowship. KH is supported by the Wellcome Trust (207569/Z/17/Z).
Footnotes
↵* Joint senior authors.
References
