# Revealing the Micro-scale Signature of Endemic Zoonotic Disease Transmission in an African Urban Setting.

Hervé Bourhy, Emmanuel Nakouné, Matthew Hall, Pierre Nouvellet, Anthony Lepelletier, Chiraz Talbi, Laurence Watier, Edward C Holmes, Simon Cauchemez, Philippe Lemey, et al.

► To cite this version:

Hervé Bourhy, Emmanuel Nakouné, Matthew Hall, Pierre Nouvellet, Anthony Lepelletier, et al. Revealing the Micro-scale Signature of Endemic Zoonotic Disease Transmission in an African Urban Setting. PLoS Pathogens, 2016, 12 (4), pp.e1005525. 10.1371/journal.ppat.1005525.t001. pasteur-01429570

HAL Id: pasteur-01429570

https://pasteur.hal.science/pasteur-01429570

Submitted on 8 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

CrossMark
click for updates

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Bourhy H, Nakouné E, Hall M, Nouvellet P, Lepelletier A, Talbi C, et al. (2016) Revealing the Micro-scale Signature of Endemic Zoonotic Disease Transmission in an African Urban Setting. PLoS Pathog 12(4): e1005525. doi:10.1371/journal.ppat.1005525

**Editor:** Colin Parrish, Cornell University, UNITED STATES

**Received:** September 19, 2015

**Accepted:** March 3, 2016

**Published:** April 8, 2016

**Copyright:** © 2016 Bourhy et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files. Sequence data are available from GenBank under accession numbers KT119634-KT119795.

**Funding:** This work was supported by European Union Seventh Framework Programme PREDEMICS (grant number 278433). ECH is supported by an NHMRC Australia Fellowship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

RESEARCH ARTICLE

# Revealing the Micro-scale Signature of Endemic Zoonotic Disease Transmission in an African Urban Setting

Hervé Bourhy¹*, Emmanuel Nakouné², Matthew Hall³,⁴, Pierre Nouvellet⁴,⁵, Anthony Lepelletier¹, Chiraz Talbi¹, Laurence Watier⁶,⁷, Edward C. Holmes⁸, Simon Cauchemez⁴,⁵,⁹, Philippe Lemey¹⁰, Christl A. Donnelly⁴,⁵, Andrew Rambaut³,¹¹

1 Institut Pasteur, Unit Lyssavirus Dynamics and Host Adaptation, WHO Collaborating Centre for Reference and Research on Rabies, Paris, France, 2 Institut Pasteur de Bangui, Bangui, République Centrafricaine, 3 Institute of Evolutionary Biology, University of Edinburgh, Ashworth Laboratories, Edinburgh, United Kingdom, 4 Department of Infectious Disease Epidemiology, Imperial College London, London, United Kingdom, 5 Medical Research Council Centre for Outbreak Analysis and Modelling, Department of Infectious Disease Epidemiology, Imperial College London, London, United Kingdom, 6 INSERM, UMR 1181 and Institut Pasteur, B2PHI, Paris, France, 7 Faculté de Médecine Paris Ile de France-Ouest, Université de Versailles-Saint-Quentin, Versailles, France, 8 Marie Bashir Institute for Infectious Diseases & Biosecurity, Charles Perkins Centre, School of Life and Environmental Sciences and Sydney Medical School, Sydney, Australia, 9 Mathematical Modelling of Infectious Diseases Unit, Institut Pasteur, Paris, France, 10 Rega Institute, KU Leuven, Leuven, Belgium, 11 Fogarty International Center, National Institutes of Health, Bethesda, Maryland, United States of America

* hbourhy@pasteur.fr; herve.bourhy@pasteur.fr

## Abstract

The development of novel approaches that combine epidemiological and genomic data provides new opportunities to reveal the spatiotemporal dynamics of infectious diseases and determine the processes responsible for their spread and maintenance. Taking advantage of detailed epidemiological time series and viral sequence data from more than 20 years reported by the National Reference Centre for Rabies of Bangui, the capital city of Central African Republic, we used a combination of mathematical modeling and phylogenetic analysis to determine the spatiotemporal dynamics of rabies in domestic dogs as well as the frequency of extinction and introduction events in an African city. We show that although dog rabies virus (RABV) appears to be endemic in Bangui, its epidemiology is in fact shaped by the regular extinction of local chains of transmission coupled with the introduction of new lineages, generating successive waves of spread. Notably, the effective reproduction number during each wave was rarely above the critical value of 1, such that rabies is not self-sustaining in Bangui. In turn, this suggests that rabies at local geographic scales is driven by human-mediated dispersal of RABV among sparsely connected peri-urban and rural areas as opposed to dispersion in a relatively large homogenous urban dog population. This combined epidemiological and genomic approach enables development of a comprehensive framework for understanding disease persistence and informing control measures, indicating that control measures are probably best targeted towards areas neighbouring the city that appear as the source of frequent incursions seeding outbreaks in Bangui.
**Competing Interests:** The authors have declared that no competing interests exist.

## Author Summary

Rabies is an endemic and zoonotic disease circulating in most African and Asian countries and is mainly transmitted by dogs. Despite its importance for public health, little is known about the spatiotemporal dynamics of rabies in dogs in an urban setting. Using spatial, temporal and virus genetic information, we reconstructed the fine-scale evolutionary and epidemiological dynamics of RABV among the dog population in a major African city—Bangui in the Central African Republic. Notably, we document that rabies is not self-sustaining in the city, and that numerous reintroductions of the virus occur in the city, generating a succession of epidemic waves. The data obtained here will inform the development of tailored control measures and avoid wasting resources, revealing that control measures are probably best targeted towards neighbouring geographic areas rather than to the city itself.

## Introduction

Understanding the processes responsible for the maintenance and extinction of infectious diseases in specific geographic locations is essential for the establishment of locally relevant control strategies. Resources for infectious disease surveillance are always in limited supply, especially those dedicated to zoonotic neglected tropical diseases [1]. As a consequence, a lack of comprehensive data on incidence and transmission chains often prevents the accurate inference of disease dynamics and patterns of transmission. Fortunately, the development of approaches combining both epidemiological and genomic data provides new opportunities to explore this problem at a fine scale offering a robust framework to account for the highly stochastic nature of both the disease transmission and observation processes.

Our focus is dog mediated rabies, one of the most neglected of all infectious diseases. Every year, nearly 15 million people receive post exposure prophylaxis and an estimated 59,000 people die from rabies [2], with 98% of these resulting from dog bites and with a burden disproportionately affecting poor rural communities and children. However, regional and national priorities for research and control are difficult to establish due to poor background data. A number of major challenges are apparent [3]. In particular, rabies surveillance systems are often non-existent or inadequate due to the lack of rabies awareness among health professionals and the target population. In addition, there is a lack of legislation for compulsory notification and an absence of laboratory confirmation leading to unreliable diagnoses. In addition, the long-term presence of rabies relies on the zoonotic source of transmission and domestic dog populations remain the principal target for prevention and interventions [4]. Evidently, uncertainty surrounding the transmission dynamics and maintenance of dog rabies must be addressed to plan, optimize and implement strategies. Intrinsic cycles and traveling waves linked to rabies-induced fluctuations in susceptible animal populations are observed in wild carnivores in Europe and North America [5,6]. In contrast, little is known about the spatial and temporal dynamics of rabies in domestic dogs at a large scale [7–10], and almost nothing at local scales, including the factors that lead to disease endemicity in an urban setting. Given the tight social and spatial structure of dog populations, the dynamics of rabies in dog populations are often considered in the context of meta-populations where endemic populations provide a source of infection for sink populations [11,12]. In such settings, targeted interventions would clearly benefit enormously from identifying these sources of infection.

Herein, we utilize a unique combination of epidemiological data and virus genome sequence data to determine the local population dynamics of a neglected tropical disease. We apply this
method to determine the local population dynamics of rabies, particularly its rate of importa-
tion and transmission potential at a highly localized scale in Bangui, the capital city of Central
African Republic (CAR), information that will be critical to its overall control.

Results

Periodicity of rabies in Bangui

Wavelet analysis and periodograms were performed on the data collected on dogs in Bangui from 1990 to 2012. Periods of 53.4 and 89.0 months were observed. No cases were detected over periods of 289 days, 470 days and 388 days, between 1996–1997, 2005–2006 and 2010–2012, respectively. Such long ‘epidemiologically silent’ intervals are intriguing and, given the mean incubation period plus infectious period of rabies in dogs (25.1 days) [7], suggest that the virus may not persist in the city between consecutive epidemiological periods.

Phylogeny and time-scale of RABV evolution

To further investigate the spatiotemporal dynamics of rabies in CAR, and to verify the likely periodic introduction of RABV in Bangui, rabies virus (RABV) isolates from Bangui collected from 2003 to 2012 were sequenced and analysed, together with those of viruses obtained earlier (in 1992) from Bangui and from neighboring or geographically close countries (Fig 1, Table A and Fig A in S1 Text). We first inferred the evolutionary history of rabies virus in this area using a set of 162 RABV sequences (5061 nt; encompassing the N, P, M, G and intergenic G-L region) (Fig 2). The tree topology is similar to earlier phylogenetic analyses, with sequences divided into two major clusters identified as the ‘Africa 1’ (AF1) clade which represents an African subset of the ‘Cosmopolitan’ clade, and an ‘Africa 2’ (AF2) clade [13–15]. Most of the

<figure><img src="image_4.png" /></figure>

**Fig 1. Sampling of RABV samples.** Map of the Central African Republic and neighbouring countries (a) showing the location (red points and squares) of the 162 isolates of RABV analyzed in this study. The cities in red indicate the presence of RABV belonging to the Africa 1 clade. Cities in blue represent those where RABV belonging to the Africa 2 clade were found. In Bangui, both clades were identified. Detailed map of Bangui (b), showing the location of the 122 samples analyzed in the context of the landscape. Number of specimens positive for rabies (c) by trimester during sentinel surveillance, January 1993 to March 2012.

doi:10.1371/journal.ppat.1005525.g001
isolates reported in Bangui belonged to the AF1 clade. Root-to-tip linear regression analysis of the 2003–2012 CAR sequences from those two clusters established that there was a strong temporal signal in these data, with correlation coefficients of 0.79 (AF1) and 0.71 (AF2) for the relationship between root-to-tip divergence and sampling date. We then estimated the time-scale of the evolutionary history of the 162 RABV isolates using a Bayesian Markov chain Monte Carlo approach [16]. The mean rate of nucleotide substitution estimated in this

<figure><img src="image_5.png" /><figcaption>**Fig 2.** MCC tree of 162 sequences from the Central African Republic and other locations in Africa estimated from 5000 nt of dog RABV genome.
Tips representing isolates from Bangui are coloured according to the selected subtypes of RABV; other tips are colored by location. Tip times are scaled to the date of sampling (years) and branches are estimated in time units as indicated by the time bar.</figcaption></figure>

doi:10.1371/journal.ppat.1005525.g002
manner, assuming an uncorrelated lognormal relaxed molecular clock model, was 5.9x10⁻⁴
substitutions per site, per year (95% highest posterior density (HPD) interval, 4.4–7.5x10⁻⁴
subs/site/year). These rates are within the range of those previously estimated for the N and G
genes in various dog RABV sample data sets in Africa [9,15] and worldwide [17].

Evidence of several independent introductions of RABV into Bangui

Clades in the RABV phylogeny were regarded as representing likely new introductions into Ban-
gui if it was possible to go back a time *L* from their root node without encountering either the
start of sampling in 2003, or a node that could not be ruled out as representing a transmission to
a dog that was sampled. *L* is intended to be a time period that is sufficiently long so that the prob-
ability was very low that the virus had been transmitted in Bangui, while sampling was going on,
but without any dogs in the transmission chain being sampled. In addition, we made the *a priori*
assumption that there were at least three introductions comprising AF2 and the two waves of
AF1 (2003–2005 and 2006–2012). This was because surveillance only began in 2003, such that
for larger values of *L* it was impossible to go back *L* years from the estimated times of the root
nodes of AF2 or the second AF1 wave without encountering the start of sampling, and hence our
algorithm alone could never determine that undetected passage for *L* years was unlikely. Not
making this assumption was deemed overly conservative due to the absence of any detected clini-
cal cases in the city for 15.5 months between the two waves. Our initial assumptions were there-
fore that there were at least three separate introductions, but for values of *L* of 1, 2, 3, 4 and 5
years the analysis identified 7, 6, 5, 5 and 5, respectively; even very conservative choices of *L* indi-
cated more than were obvious from the epidemiological data. We designate these introductions
as subtypes (ST): variants of RABV circulating during the time frame of the study.

To determine a plausible value of *L*, we used estimates of the assumed detection probability *ρ* (range 10%-50%) to account for low surveillance efforts and combined them with data concerning the incubation plus infectious period reported for dogs in Africa [8]. From this, we were able to estimate the probability of the presence of a ST in Bangui without any dog in the transmission chain being sampled. The expected number of dogs in the chain is 14.37 in one year and 28.89 in two years. The probability of not detecting one of these dogs up after one year is 0.22 for a 10% sampling fraction, 0.0405 for 20% and 4.71x10⁻⁵ for 50%. At two years, the probability is less than 5% (0.0477) even for a 10% value of ρ. A value for *L* of two years therefore seems reasonable, from which we can infer at least 6 separate introductions of RABV occurred in Bangui during the time course of this study (2003–2012) (Fig 3). Importantly, this estimate may be conservative because sampling only started in 2003, and the posterior median TMRCA of all the isolates collected in Bangui during 2003–2004 occurs during 1998, such that it was impossible to determine whether these isolates were the result of a single introduction. In addition, we cannot exclude that the same ST was repeatedly introduced into the city during the time course of the study. The six different introductions were due to 4 different STs of AF1 and 2 different STs of AF2. One viral subtype was always predominant during each of the waves but the reported ST shifted from one wave to another (Fig 3). AF1 ST1 was predominant in 2003–2005 while AF1 ST2 became predominant in 2006–2010. In 2012, another new ST appeared in Bangui (AF2 ST2). Our findings are therefore consistent with a series of distinct clusters of RABV originating from independent introduction into Bangui, each of which then experienced relatively limited transmission in the city.

Phylogeography of RABV in Bangui

To determine the spread of rabies virus in the city, we used a combination of phylogenetic and
coalescent approaches to connect the historical process of dog rabies virus evolution in Bangui
<figure><img src="image_6.png" /><figcaption>Fig 3. Temporal distribution of RABV subtypes in Bangui as determined by the MCC analysis (number of specimens positive for rabies by months from January 2003 to March 2012 and same color code as in Fig 2). Black arrows indicates RABV introduction events in the city (first introduction of a new subtype circulating in Bangui). Grey areas indicated periods without any reported cases.</figcaption></figure>

doi:10.1371/journal.ppat.1005525.g003

together with epidemiological parameters. Our aim was to integrate spatial dynamics with tem-
poral inferences to follow the movement of RABV, and hence infected dogs, in the city. A Man-
tel r-test was performed on the pairwise genetic differences versus physical distances of the
samples, giving a significant (p = 0.023) but weak (r = 0.074) positive association between the
two differences. A continuous traits phylogeographic analysis was then performed in BEAST,
using the geographical coordinates of the largest data set of samples corresponding to AF1
ST4. The analysis of how the infected areas changes with time indicated that these viruses ini-
tially arrived from the Southwest and then eventually dispersed to the centre of the city (Fig B1
in S1 Text). The virus then largely spread along a Southeast Northwest axis in the city (Fig B2
in S1 Text). The estimated diffusion rate for the AF1 ST4 strain in Bangui obtained from the
data was 0.9 km per year (95% HPD interval: 0.65–1.2 km)

Temporal trends and role of importation of rabid dogs into Bangui

To determine the population dynamics of rabies in Bangui, the instantaneous effective repro-
duction number R<sub>t</sub> and the rate of introduction λ<sub>t</sub> were estimated based on the incident cases
from 1990 to 2012. Strikingly, we observed relatively little variation in R<sub>t</sub> over time, with
median monthly point estimates in the range 0.8 to 1.3 (Fig 4, Fig C in S1 Text), with R<sub>t</sub> rarely
significantly above the critical value of 1 (only during 12 weeks in 2009 according to our best
model). Substantial temporal variation in the 95% credible intervals for R<sub>t</sub> was observed and
likely reflects variation in the number of reported incident cases rather than dramatic changes
in the underlying transmission dynamics of the disease, i.e. since R<sub>t</sub> describes infectiousness of
a case, the number of incident cases in a time period determines the precision of the estimate
[18]. The best model (Table B in S1 Text) assumed a constant exogenous rate of introduction
λ<sub>t</sub> throughout the study period, with a median estimate of 0.13 (95% credible interval CI: 0.07,
0.27) rabid dogs imported to Bangui per week (i.e. ~7 per year) (Fig C, D and E in S1 Text).
Given the posterior distribution of R<sub>t</sub> and λ<sub>t</sub>, the simulated number of rabid dogs from the
model was in good agreement with the observed counts (Fig 4c). Based on the estimated
parameters, we were able to infer the likely proportion of dogs that became rabid as a result of
a local transmission event or were imported cases. Accordingly, from 1992 until 2012 we
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Estimation of the rate of introduction of rabid dogs in Bangui, λt (in rabid dogs imported to Bangui per week), the instantaneous effective reproduction number, Rt, and observed number of rabid dogs infected locally (red bar) or from outside the city (blue bar) and simulated number of rabid dogs from the model (black line: posterior median; grey area: 95% credible interval of the posterior distribution). If Rt, is above 1, a local self-sustaining epidemic in dogs may occur in the city.</figcaption></figure>

doi:10.1371/journal.ppat.1005525.g004

inferred that among 390 rabid dogs observed, approximately 10% of them (~40) were intro-
duced in Bangui (blue bars in Fig 4, lower panel). The sensitivity analysis indicated that the
inferred pattern of temporal variation in both R_t and λ_t was robust to changes in the assumed
detection probability (Table C in S1 Text). Given an estimated median for R_t of 0.84 over the
whole study period, an outbreak seeded by a single case would end, on average, after a total of
6.3 cases (including the index case) with a variance linked to heterogeneity of transmission
[19,20]. While our analysis does not allow us to directly reconstruct the chains of infection,
given estimated values for the median R_t and λ_t, (and assuming a Poisson distribution of cluster
sizes, i.e. no heterogeneity in infectiousness), we expect the occurrence of large clusters; for
instance, over the whole study period we expect 3.2 clusters larger than 50 cases (10 observed
cases given the 20% best inferred sampling detection). Finally, between 2003 and 2012, our sta-
tistical analysis predicted that 10 observed cases had a probability greater than 90% of having
been imported. In contrast, the genetic analysis allowed us to determine that 6 new virus sub-
types were introduced during the same time period. Four importations (in late 2003, early
2006, late 2010 and early 2012) were identified by both methods. The known importations
(from genetic information) from mid-2006 and late 2007 were missed by the statistical analysis,
although not surprisingly as they arose during a large infection cluster. Finally, 6 likely
imported cases were missed by the genetic analysis, most likely reflecting multiple importations
of viruses with very similar sequences. In sum, our analyses confirms that, although rabies is
almost endemic in Bangui, its dynamics are in fact driven by a high rate of importation and
characterized by frequent cluster die-out.
Discussion

Phylogenetic and coalescent-based methods allow the explicit use of spatial information within a statistical framework [9,21,22]. When coupled with mathematical models for estimating time-varying reproduction numbers during epidemics this permits practical and operational questions to be asked. In the case of neglected and re-emerging zoonoses such as rabies, it is particularly important to understand the details of transmission dynamics to design efficient and sustainable strategies for dog rabies control and elimination. Here, we described transmission patterns of RABV infection in a dog population in Bangui, a relatively large African city with nearly 900,000 inhabitants. In this city, the domestic dog is considered to be the only reservoir of rabies virus—that is, the only set of epidemiologically connected animal populations maintaining rabies infection in a defined geographical area [23]—although actions taken to control the disease in dog populations are limited.

Our study demonstrates that the dog population of Bangui experiences the extinction of RABV transmission chains separating small outbreaks and that reflects a succession of epidemic waves, a pattern that is consistent with previous observations in Africa [7,8]. Importantly, when examined over the long term this gives the false impression of a stable endemic infection. Indeed, we noted that the instantaneous effective reproduction number $R_t$ was rarely significantly above the critical value of 1. Wide temporal variation in credible intervals for $R_t$ were observed, but are likely to reflect more the variation in incident cases reported rather than changes in the dynamics of the disease per se [18]. These estimates of $R_t$ are surprisingly low with median monthly point estimates (0.8 to 1.3) spanning the range of lower estimates of the basic reproduction number, $R_0$ (not instantaneous), obtained from compiled outbreak data from other localities [8]. These estimates of $R_t$ are consistent, although often lower than the $R_0$ estimated ($R_0$ = 1.1) for another African city, Nd'jamena in Chad [10]. It is interesting to note that the periods that correspond to clear outbreaks have $R_t$ above 1 (presumably corresponding to beginning of AF1 ST2 and of AF1 ST4) and that during the exponential growth of these outbreaks, $R_t$ corresponds to $R_0$. Although $R_t$ rarely exceeds 1, this is consistent with other endemic persistent infections (e.g. see variation in $R_t$ for polio, [24]).

The inference of successive small transmission chains is supported by our phylogenetic analysis that reveals the existence of repeated introductions of RABV from outside the city. Notably, the subtype of reported RABV circulating in the city shifted from one wave to the next, indicating a sudden occurrence and rapid disappearance for most of the ST identified. This periodic replacement of the subtypes argues against persistent diffusion of RABV through the local movement of rabid dogs in the city. Some of the subtypes infecting dogs in Bangui belonged to the Africa 2 clade, which was originally only found in the northern part of the country. Hence, this confirms the long-distance spread of RABV (several hundreds of km) likely associated with human-mediated displacement of dogs infected with rabies [9], and its role in repeated and successful re-introduction of RABV in rabies-free areas. In turn, our mathematical model supports a constant exogenous rate of introduction $\lambda_t$ throughout the study period with a median estimate at 0.13 (95% CrI: 0.07, 0.27) rabid dog imported to Bangui per week. This was substantiated by the phylogenetic analysis which allowed us to estimate that at least 6 different ST were imported in the city between 2003 and 2012. Further, the movement of infected dogs in the city as suggested by the rate of diffusion of one RABV ST in the city (0.9km/year) from the phylogeographic analysis appeared relatively low compared to other studies performed at a larger scale in Africa [8,9], indicating that the spatial spread of RABV within the urban dog population is limited. In sum, both the epidemiological and genomic data paint a consistent picture of serial importation of distinct viral clusters that experience relatively little evolution in the city.
Given that epidemiological surveillance of rabies in animals is not exhaustive, a substantial proportion of animal cases were likely never detected and therefore did not appear in our data set, in turn leading to an underestimation of the burden of rabies. However, this absence of information did not interfere with the reliability of the conclusions obtained, as these potential sampling biases were taken into account during the analyses. Further, estimates of R<sub>t</sub> over time are generally robust to underreporting as long as the proportion of underreported cases is constant over time. While assuming constant reporting was reasonable from field experience, it is likely that some variation occurred, particularly that reporting increased during periods of higher incidence as surveillance became more active. If such variations in reporting occurred, our estimate of R<sub>t</sub> during period of high incidence would be slightly overestimated, as the increase surveillance would reveal cases that would otherwise not be detected. However, we predict that any possible reporting bias would have a minimal impact on our estimate of importation as it relies more on period of low incidence. Together, our central conclusions (i.e. of disease dynamics driven by a high rate of importation and characterized by frequent cluster die-out) appear to be robust, and even conservative, to increased reporting during periods of high incidence.

The overall consistency of the results obtained by our two approaches, mathematical and phylogenetic modeling, is therefore reassuring and complements other approaches combining genetic and epidemiological data enabling the reconstruction of partially observed transmission trees and estimating the number of cases missing [25–29]. Further, our approach is designed to incorporate an exogenous force of infection and can be applicable in a wide range of settings including those where handling missing data is critical. We still have to determine how best to integrate commonly collected data types (incidence data, time series of reported cases and molecular sequence, etc.) to inform pathogen persistence and control. However, we have shown here that the combined use of mathematical and phylogenetic approaches has provided a new understanding of the transmission of RABV. This information could be further used to inform RABV prevention and response policies in defined areas in a way that maximizes cost-effectiveness and operational efficiency [1,30]. For example, it would be useful to identify areas where surveillance is critical and low reporting rates might leave epidemics undetected, and would aid the rapid implementation of public health strategies targeting identified geographical regions at risks.

Drawing upon these epidemiological and evolutionary data allows us to address more fun-
damental questions relating to the persistence of rabies and to the role of different dog popula-
tions. The size and structure of domestic dog populations in space and time can be extremely
variable and difficult to measure, but in general can be thought of as a spatial mosaic of large
and small populations that are occasionally connected. The respective roles of small metapopu-
lation-structured versus large and more homogeneous dog populations in the maintenance of a
zoonotic disease such as rabies are unknown [1]. The coupling of large and small dog popula-
tions and their spatial connectivity may therefore be the critical feature that generates waves
and persistence of RABV at a regional scale. As a consequence, the complex process of persis-
tence seems therefore to operate on multiple scales among which the metapopulation structure
of dogs could play a major role compared to densities [1,7,8,31]. The spatial distribution and
size of subpopulations can cause spatial hierarchies of metapopulation dynamics or travelling
waves of infection [32]. It is therefore essential to identify sub-populations important in
source-sink dynamics or bottleneck areas, each requiring tailored management plans [27,30].
More generally, our data demonstrate the need to understand hierarchical (‘core-satellite’) epi-
demic dynamics for the control of established and re-emerging diseases [23,33]. The explora-
tion of the potential ecological and epidemiological drivers behind these findings is a major
challenge for future research. From this perspective, the collection and study of rabies
specimens collected at a larger geographical scale including peri-urban and rural areas would be valuable. The role of metapopulations and in particular of the minimum viable metapopulation (MVM) size has also been addressed in rabies [34] and in other systems with short infectious period such as measles, in which large focal communities have been proposed to act as refuges that allow metapopulation persistence [35,36].

Phylogenetic and coalescent-based methods coupled with mathematical models also allow us to address practical and operational questions regarding rabies control, particularly the optimal strategies for vaccination and allocations in situations where resources are limited. The dynamics of rabies and the effectiveness of vaccination for reducing the occurrence of rabies are strongly influenced by the spatial distribution and the metapopulation structure of dog population [11,12]. According to our study, control measures are best directed toward preventing spill-over from neighboring populations and the introduction of novel RABVs into this target dog population, rather than only targeting transmission within the urban dog populations.

# Materials and Methods

## Epidemiological data collection and selection

The State health department of CAR conducted all observations and euthanization (if necessary) of animals suspected of rabies in Bangui, according to established national standardized protocols and WHO guidelines [37]. The Institut Pasteur of Bangui is the only laboratory in the country that is designated by the Ministry of Health of CAR to conduct the collection of specimens from dead animals suspected of rabies and to perform rabies laboratory diagnosis [38]. Rabies diagnosis was performed by direct immunofluorescence on smears of suspected animals, and confirmed in parallel using RT-PCR as previously described [37,39]. Confirmed animal cases from CAR were collated over the period 1990–2012. The data collected include the date of death, the species, the animal's city of origin and starting from 2003 the spatial coordinates for all the isolates originating from Bangui (Table A in S1 Text).

## Estimation of temporal trends of rabies

To determine periodic components of the monthly case time series, we performed a spectral analysis [40] on data collected from 1990 to 2012. We then used the number of dogs reported as rabid in Bangui and their date of submission to rabies diagnosis from 1992 to 2012 to infer the dynamics of rabies transmission. The instantaneous effective reproduction number $R_t$ is the average number of dogs that a single dog infected at time $t$ could be expected to infect should conditions remain unchanged. Simple methods were developed to estimate trends in $R_t$ from the epidemic curve of incident cases in the context of closed epidemics [18,41]. Here, we expanded these approaches for the situation where there is an exogenous force of infection. The number of incident cases at time $t$, $I_t$, was therefore modelled as a Poisson process with mean $\lambda_t + R_t \sum_{s=1}^{t} I_{t-s}$, where $\lambda_t$ was the rate of introduction of rabies into the city and $w_S$ was the distribution of the serial interval. The serial time interval, accounting for the incubation and infectious periods, was assumed to be constant over time and to follow a Gamma distribution parameterised using data presented in Hampson et al. [8] with a mean serial interval of 25 days and a variance of 171 days². Further, the number of observed rabid dogs at time $t$ followed a binomial distribution $Bin(I_t, \rho)$, where the detection probability $\rho$ was assumed to be constant over time. Given the data, the unit time of the inference model was the week; thus data on incidence were aggregated weekly (e.g. $I_t$ corresponds to the number of incident case during week $t$). Both $R_t$ and $\lambda_t$ were initially allowed to vary according to random walks characterised by
their respective variances parameters ($\sigma_{R_t}^2$ and $\sigma_{\lambda_t}^2$) (see [42, 43, 44] for details). Specifically, the value $R_t$ was modified using $R_t = R_{t-\tau}+\epsilon$ with $\epsilon$ a random variable following a normal distribution with mean zero and variance $\sigma_{R_t}^2$ and $\tau$ the time window at which $R_t$ was assumed to change (the variance $\sigma_{R_t}^2$ was increased during the 3 largest outbreaks to allow larger volatility in $R_t$, see Fig C, D and E in S1 Text). The same procedure was applied to $\lambda_t$. Given the limited data available, we explored the monthly variation of $R_t$ ($\tau = 1$ month) and 3 models for $\lambda_t$: i) a constant rate of rabies introduction throughout the period ($\lambda_t = \lambda$), ii) quarterly variation in $\lambda_t$ ($\tau = 3$ month), to explore potential seasonal trends and iii) annual variation in $\lambda_t$ ($\tau = 12$ month) to explore longer term trends. We estimated temporal trends in $R_t$ and $\lambda_t$ using a particle Monte Carlo Markov chain [42]. We report the posterior medians and 95% credible intervals of parameters. For each iteration of Markov chain we could extract a particular trajectory for the incidence (drawn randomly from the set of particles filtered at the end of the time series). The set of trajectories obtained constitutes the posterior incidence simulated. As all models had the same number of parameters, we compared them using an estimate of the mean deviance. A sensitivity analysis was performed to assess robustness of estimates to changes in the assumed detection probability $\rho$ (range 10%-50%). To further evaluate and validate our model, we linked our results to branching theory that gives predictions about the cluster size due to single introduction, i.e. the total of number of cases (including the index case) after an outbreak is seeded by a single case and end [45]. Finally, for each day 't', given estimation of $R_t$, $\lambda_t$ and the simulated trajectories, we could derive the force of infection due to either local transmission ($R_t \sum_{s=1}^t I_{t-s}w_s$) or imported cases ($\lambda_t$). This allowed us to determine the probability that an observed case was imported rather than due to a local transmission event, and we analyzed those results in relation to the available genetic information of introduction of new virus subtypes.

## Bayesian evolutionary analysis of spatiotemporal RABV dynamics

A total of 162 RABV isolates were analyzed from which 139 originated from CAR, including 123 isolates collected during the period 2003–2012, 2 in 1992 from Bangui, and 14 from other cities in CAR (Fig 1, Table A in S1 Text). To complete the data set, sequences of 25 isolates from neighboring countries (Democratic Republic of the Congo, Chad, Cameroon, Gabon and Zambia) were also determined. Total RNA from the original brain samples (when available) was extracted using Trizol reagent (Invitrogen) according to the manufacturer's instructions. To perform a reliable evolutionary analysis of dog RABV circulating in Bangui, we sequenced a large genomic area (42% of the viral genome) encompassing the N, P, M, G and intergenic G-L region, corresponding to positions 161 to 5221 of the RABV genome [46]. RNA extractions, RT-PCR and sequencing were performed as described previously [9, 14, 15, 17] and with additional primers (Table D in S1 Text). All 162 sequences were aligned using the MUSCLE [47] plugin for Geneious 5.6.6 (Biomatters). Unrooted maximum likelihood phylogenetic trees for the entire 162-sequence alignment and for the AF1 and AF2 clusters in the CAR data were estimated using RAxML [48] employing 100 bootstrap replicates and the GTR+GAMMA nucleotide substitution model. The AF1 and AF2 phylogenies were used to investigate the clock-like behaviour in the data using the Path-O-Gen program, which maximises the Pearson product-moment correlation coefficient between sampling date and root-to-tip divergence for each isolate. Time-scaled phylogenetic trees were estimated using BEAST v1.7.5 [16] using the SRD06 nucleotide substitution model [49,50], an uncorrelated relaxed molecular clock [51] and a GMRF Bayesian skyride tree prior [52]. The
analysis used multiple independent chains of 100 million MCMC states each, with samples taken sampling every 10,000 states. The convergence of parameter estimations was observed in all cases. The posterior set of phylogenetic trees from BEAST was used to identify separate introductions of RABV to Bangui; for full details see S1 Text. This procedure involves two numerical parameters, *L* and *P*. *L* is the length of time, in years, which is long enough that if rabies were being transmitted down a single chain while the population of infected dogs in the city was being sampled (with probability less than 1) upon death, we would expect to take at least one sample from a dog in the chain. *P* is a conservative estimate of the time from infection to death of an infected dog. In determining plausible values for *L* and *P*, we used previously-estimated figures for the distribution of incubation and infection periods of RABV in dogs [8] and the Welch-Satterthwaite approximation for the distribution of the sum of gamma-distributed random variables to estimate the distribution of the time from infection to death. The 99th percentile of this distribution corresponded to 98.5 days, and we used this value for *P*. We also used this distribution to calculate the probability of sampling no dogs in a transmission chain for various values of *L* and sampling proportions.

To analyse the diffusion of RABV in the city of Bangui, a phylogeographic analysis of the 88 sequences of AF1 ST4 (i.e. the most common ST) was then conducted, also in BEAST, using the same nucleotide substitution model and tree prior as described above, along with a continuous spatial diffusion model [22]. This analysis utilizes the latitude and longitude at which each virus sample was taken, and uses them, in combination with the phylogeny, to estimate the parameters of a process by which viral lineages diffuse across the landscape and infer the locations of internal nodes. Branch-specific variation in diffusion rates was assumed to be lognormally-distributed. The prior distribution on the precision matrix of the diffusion process was the Wishart distribution described in [22] and an Exp(3) distribution was used for the standard deviation of the lognormal distribution modelling branch-specific diffusion rates. The SPREAD utility [53] was used to identify regions of Bangui with a 95% posterior probability for the presence of viral lineages, at different points in time, and these were visualised using Google Earth.

**Supporting Information**

**S1 Text. Supporting information legends.** Table A: Relevant epidemiological information for all RABV isolates analyzed in this study. Text S1: Bayesian evolutionary analysis of spatiotemporal RABV dynamics. Table B: Summary of the fit for models where *R*<sub>t</sub> varied monthly while *λ*<sub>t</sub> was either kept constant varied every 3 months or yearly. Table C: Sensitivity analysis on the assumption on the surveillance effort *ρ*. Table D: List of primers used in this study. Figure A: Detailed figure showing the MCC tree of 162 sequences from the Central African Republic and other locations in Africa estimated from 5000 nt of dog RABV genome with names of the tips. Tips representing isolates from Bangui are coloured according to the selected subtypes of RABV; other tips are colored by location. Tip times are scaled to the date of sampling (years) and branches are estimated in time units as indicated by the time bar. Posterior clade probability values (>0.9) are shown for key nodes. Figure B: Diffusion of RABV in Bangui. Figure B1: 80% highest posterior density (HPD) regions (red areas) for areas of Bangui where the AF1 ST4 strain of RABV was present, at different times during the epidemic, taken from the phylogeographical analysis. Figure B2: The maximum clade credibility tree for the ST1 AF4 RABV sequences (2007–2010), superimposed on a map of Bangui. Figure C: Fitting results with quarterly variation of *λ*<sub>t</sub>. Figure D: Fitting results with yearly variation of *λ*<sub>t</sub>. Figure E: Fitting results with constant importation *λ*<sub>t</sub>.

(DOCX)
## Acknowledgments

The authors are grateful to Jean-Marc Reynes for providing the isolates collected from Cameroon during the period 2009–2011.

## Author Contributions

Conceived and designed the experiments: HB EN MH PN SC PL CAD AR. Performed the experiments: HB MH PN AL CT LW SC. Analyzed the data: HB EN MH PN AL CT LW ECH SC PL CAD AR. Contributed reagents/materials/analysis tools: HB EN MH PN SC CAD AR. Wrote the paper: HB EN MH PN LW ECH SC PL CAD AR.

## References

1. Brunker K, Hampson K, Horton DL, Biek R (2012) Integrating the landscape epidemiology and genetics of RNA viruses: rabies in domestic dogs as a model. *Parasitology* 139: 1899–1913. doi: 10.1017/S003118201200090X PMID: 22814380

2. Hampson K, Coudeville L, Lembo T, Sambo M, Kieffer A, et al. (2015) Estimating the global burden of endemic canine rabies. *PLoS Negl Trop Dis* 9: e0003709. doi: 10.1371/journal.pntd.0003709 PMID: 25881058

3. Dodet B, Africa Rabies B (2009) The fight against rabies in Africa: From recognition to action. *Vaccine* 27: 5027–5032. doi: 10.1016/j.vaccine.2009.06.030 PMID: 19560430

4. Lembo T, Partners for Rabies P (2012) The blueprint for rabies prevention and control: a novel operational toolkit for rabies elimination. *PLoS Negl Trop Dis* 6: e1388. doi: 10.1371/journal.pntd.0001388 PMID: 22389727

5. Anderson RM, Jackson HC, May RM, Smith AM (1981) Population dynamics of fox rabies in Europe. *Nature* 289: 765–771. PMID: 7464941

6. Russell CA, Smith DL, Waller LA, Childs JE, Real LA (2004) A priori prediction of disease invasion dynamics in a novel environment. *Proc Biol Sci* 271: 21–25. PMID: 15002767

7. Hampson K, Dushoff J, Bingham J, Bruckner G, Ali YH, et al. (2007) Synchronous cycles of domestic dog rabies in sub-Saharan Africa and the impact of control efforts. *Proc Natl Acad Sci U S A* 104: 7717–7722. PMID: 17452645

8. Hampson K, Dushoff J, Cleaveland S, Haydon DT, Kaare M, et al. (2009) Transmission dynamics and prospects for the elimination of canine rabies. *PLoS Biol* 7: e53. doi: 10.1371/journal.pbio.1000053 PMID: 19278295

9. Talbi C, Lemey P, Suchard MA, Abdelatif E, Elharrak M, et al. (2010) Phylodynamics and human-mediated dispersal of a zoonotic virus. *PLoS Pathog* 6: e1001166. doi: 10.1371/journal.ppat.1001166 PMID: 21060816

10. Zinsstag J, Durr S, Penny MA, Mindekem R, Roth F, et al. (2009) Transmission dynamics and economics of rabies control in dogs and humans in an African city. *Proc Natl Acad Sci U S A* 106: 14996–15001. doi: 10.1073/pnas.0904740106 PMID: 19706492

11. Beyer HL, Hampson K, Lembo T, Cleaveland S, Kaare M, et al. (2011) Metapopulation dynamics of rabies and the efficacy of vaccination. *Proc Biol Sci* 278: 2182–2190. doi: 10.1098/rspb.2010.2312 PMID: 21159675

12. Beyer HL, Hampson K, Lembo T, Cleaveland S, Kaare M, et al. (2012) The implications of metapopulation dynamics on the design of vaccination campaigns. *Vaccine* 30: 1014–1022. doi: 10.1016/j.vaccine.2011.12.052 PMID: 22198516

13. Durr S, Naissengar S, Mindekem R, Diguimbye C, Niezgoda M, et al. (2008) Rabies diagnosis for developing countries. *PLoS Negl Trop Dis* 2: e206. doi: 10.1371/journal.pntd.0000206 PMID: 18365035

14. Kissi B, Tordo N, Bourhy H (1995) Genetic polymorphism in the rabies virus nucleoprotein gene. *Virology* 209: 526–537. PMID: 7778285

15. Talbi C, Holmes EC, de Benedictis P, Faye O, Nakoune E, et al. (2009) Evolutionary history and dynamics of dog rabies virus in western and central Africa. *J Gen Virol* 90: 783–791. doi: 10.1099/vir.0.007765-0 PMID: 19264663

16. Drummond AJ, Suchard MA, Xie D, Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST 1.7. *Mol Biol Evol* 29: 1969–1973. doi: 10.1093/molbev/mss075 PMID: 22367748
17. Bourhy H, Reynes JM, Dunham EJ, Dacheux L, Larrous F, et al. (2008) The origin and phylogeography of dog rabies virus. J Gen Virol 89: 2673-2681. doi: 10.1099/vir.0.2008/003913-0 PMID: 18931062

18. Cori A, Ferguson NM, Fraser C, Cauchemez S (2013) A new framework and software to estimate time-varying reproduction numbers during epidemics. Am J Epidemiol 178: 1505-1512. doi: 10.1093/aje/kwt133 PMID: 24043437

19. Blumberg S, Lloyd-Smith JO (2013) Comparing methods for estimating R0 from the size distribution of subcritical transmission chains. Epidemics 5: 131-145. doi: 10.1016/j.epidem.2013.05.002 PMID: 24021520

20. Garske T, Rhodes CJ (2008) The effect of superspreading on epidemic outbreak size distributions. J Theor Biol 253: 228-237. doi: 10.1016/j.jtbi.2008.02.038 PMID: 18423673

21. Lemey P, Rambaut A, Drummond AJ, Suchard MA (2009) Bayesian phylogeography finds its roots. PLoS Comput Biol 5: e1000520. doi: 10.1371/journal.pcbi.1000520 PMID: 19779555

22. Lemey P, Rambaut A, Welch JJ, Suchard MA (2010) Phylogeography takes a relaxed random walk in continuous space and time. Mol Biol Evol 27: 1877-1885. doi: 10.1093/molbev/msq067 PMID: 20203288

23. Haydon DT, Cleaveland S, Taylor LH, Laurenson MK (2002) Identifying reservoirs of infection: a conceptual and practical challenge. Emerg Infect Dis 8: 1468-1473. PMID: 12498665

24. Grassly NC, Fraser C (2006) Seasonal infectious disease epidemiology. Proceedings of the Royal Society B-Biological Sciences 273: 2541-2550.

25. Jombart T, Cori A, Didelot X, Cauchemez S, Fraser C, et al. (2014) Bayesian reconstruction of disease outbreaks by combining epidemiologic and genomic data. PLoS Comput Biol 10: e1003457. doi: 10.1371/journal.pcbi.1003457 PMID: 24465202

26. Mollentze N, Nel LH, Townsend S, le Roux K, Hampson K, et al. (2014) A Bayesian approach for inferring the dynamics of partially observed endemic infectious diseases from space-time-genetic data. Proc Biol Sci 281: 20133251. doi: 10.1098/rspb.2013.3251 PMID: 24619442

27. Mollentze N, Weyer J, Markotter W, le Roux K, Nel LH (2013) Dog rabies in southern Africa: regional surveillance and phylogeographical analyses are an important component of control and elimination strategies. Virus Genes 47: 569-573. doi: 10.1007/s11262-013-0974-3 PMID: 23996607

28. Morelli MJ, Thebaud G, Chadoeuf J, King DP, Haydon DT, et al. (2012) A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data. PLoS Comput Biol 8: e1002768. doi: 10.1371/journal.pcbi.1002768 PMID: 23166481

29. Ypma RJ, Bataille AM, Stegeman A, Koch G, Wallinga J, et al. (2012) Unravelling transmission trees of infectious diseases by combining genetic and epidemiological data. Proc Biol Sci 279: 444-450. doi: 10.1098/rspb.2011.0913 PMID: 21733899

30. Gerardo-Giorda L, Puggioni G, Rudd RJ, Waller LA, Real LA (2013) Structuring targeted surveillance for monitoring disease emergence by mapping observational data onto ecological process. J R Soc Interface 10: 20130418. doi: 10.1098/rsif.2013.0418 PMID: 23864503

31. Morters MK, Restif O, Hampson K, Cleaveland S, Wood JL, et al. (2013) Evidence-based control of canine rabies: a critical review of population density reduction. J Anim Ecol 82: 6-14. doi: 10.1111/j.1365-2656.2012.02033.x PMID: 23004351

32. Grenfell B, Harwood J (1997) (Meta)population dynamics of infectious diseases. Trends Ecol Evol 12: 395-399. PMID: 21238122

33. Grenfell BT, Bjornstad ON, Kappey J (2001) Travelling waves and spatial hierarchies in measles epidemics. Nature 414: 716-723. PMID: 11742391

34. Viana M, Mancy R, Biek R, Cleaveland S, Cross PC, et al. (2014) Assembling evidence for identifying reservoirs of infection. Trends Ecol Evol 29: 270-279. doi: 10.1016/j.tree.2014.03.002 PMID: 24726345

35. Gunning CE, Wearing HJ (2013) Probabilistic measures of persistence and extinction in measles (meta)populations. Ecol Lett 16: 985-994. doi: 10.1111/ele.12124 PMID: 23782847

36. Hanski I, Moilanen A, Gylenberg M (1996) Minimum viable metapopulation size. The American Naturalist 147: 527-541.

37. WHO (2013) WHO Expert Consultation on Rabies: second report. Geneva: World Health Organization. 139 p.

38. Nakoune E, Digol M, Konamna X, Selekon B, Le Faou A (2012) New introduction and spread of rabies among dog population in Bangui. Acta Trop 123: 107-110. doi: 10.1016/j.actatropica.2012.04.005 PMID: 22569561
39. Dacheux L, Reynes JM, Buchy P, Sivuth O, Diop BM, et al. (2008) A reliable diagnosis of human rabies based on analysis of skin biopsy specimens. Clin Infect Dis 47: 1410-1417. doi: 10.1086/592969
PMID: 18937576

40. Chatfield C (1989) The analysis of time series: An introduction. In: Chapman & Hall L, editor. Texts in Statistical Science. London.

41. Fraser C (2007) Estimating individual and household reproduction numbers in an emerging epidemic. PLoS One 2: e758. PMID: 17712406

42. Andrieu C, Doucet A, Holenstein R (2010) Particle Markov chain Monte Carlo methods. JR Statist Soc B 72: 269-342.

43. Ionides EL, Breto C, King AA (2006) Inference for nonlinear dynamical systems. Proc Natl Acad Sci U S A 103: 18438-18443. PMID: 17121996

44. Rasmussen DA, Ratmann O, Koelle K (2011) Inference for nonlinear epidemiological models using genealogies and time series. PLoS Comput Biol 7: e1002136. doi: 10.1371/journal.pcbi.1002136
PMID: 21901082

45. Ferguson NM, Fraser C, Donnelly CA, Ghani AC, Anderson RM (2004) Public health. Public health risk from the avian H5N1 influenza epidemic. Science 304: 968-969. PMID: 15143265

46. Delmas O, Holmes EC, Talbi C, Larrous F, Dacheux L, et al. (2008) Genomic diversity and evolution of the lyssaviruses. PLoS ONE 3: e2057. doi: 10.1371/journal.pone.0002057 PMID: 18446239

47. Edgar RC (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res 32: 1792-1797. PMID: 15034147

48. Stamatakis A (2014) RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30: 1312-1313. doi: 10.1093/bioinformatics/btu033 PMID: 24451623

49. Hasegawa M, Kishino H, Yano T (1985) Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol 22: 160-174. PMID: 3934395

50. Shapiro B, Rambaut A, Drummond AJ (2006) Choosing appropriate substitution models for the phylogenetic analysis of protein-coding sequences. Mol Biol Evol 23: 7-9. PMID: 16177232

51. Drummond AJ, Ho SY, Phillips MJ, Rambaut A (2006) Relaxed phylogenetics and dating with confidence. PLoS Biol 4: e88. PMID: 16683862

52. Minin VN, Bloomquist EW, Suchard MA (2008) Smooth skyride through a rough skyline: Bayesian coalescent-based inference of population dynamics. Mol Biol Evol 25: 1459-1471. doi: 10.1093/molbev/msn090 PMID: 18408232

53. Bielejec F, Rambaut A, Suchard MA, Lemey P (2011) SPREAD: spatial phylogenetic reconstruction of evolutionary dynamics. Bioinformatics 27: 2910-2912. doi: 10.1093/bioinformatics/btr481 PMID: 21911333