Simultaneous inference of phylogenetic and transmission trees in infectious disease outbreaks

Whole-genome sequencing of pathogens from host samples becomes more and more routine during infectious disease outbreaks. These data provide information on possible transmission events which can be used for further epidemiologic analyses, such as identification of risk factors for infectivity and transmission. However, the relationship between transmission events and sequence data is obscured by uncertainty arising from four largely unobserved processes: transmission, case observation, within-host pathogen dynamics and mutation. To properly resolve transmission events, these processes need to be taken into account. Recent years have seen much progress in theory and method development, but existing applications make simplifying assumptions that often break up the dependency between the four processes, or are tailored to specific datasets with matching model assumptions and code. To obtain a method with wider applicability, we have developed a novel approach to reconstruct transmission trees with sequence data. Our approach combines elementary models for transmission, case observation, within-host pathogen dynamics, and mutation, under the assumption that the outbreak is over and all cases have been observed. We use Bayesian inference with MCMC for which we have designed novel proposal steps to efficiently traverse the posterior distribution, taking account of all unobserved processes at once. This allows for efficient sampling of transmission trees from the posterior distribution, and robust estimation of consensus transmission trees. We implemented the proposed method in a new R package phybreak. The method performs well in tests of both new and published simulated data. We apply the model to five datasets on densely sampled infectious disease outbreaks, covering a wide range of epidemiological settings. Using only sampling times and sequences as data, our analyses confirmed the original results or improved on them: the more realistic infection times place more confidence in the inferred transmission trees.


Introduction
As sequencing technology becomes easier and cheaper, detailed outbreak investigation increasingly involves the use of pathogen sequences from host samples [1]. These sequences can be used for studies ranging from virulence or resistance related to particular genes [1,2], to the interaction of epidemiological, immunological and evolutionary processes on the scale of populations [3,4]. If most or all hosts in an outbreak are sampled, it is also possible to use differences in nucleotides, i.e. single-nucleotide polymorphisms (SNPs), to resolve transmission clusters, individual transmission events, or complete transmission trees. With that information it becomes possible to identify high risk contacts and superspreaders, as well as characteristics of hosts or contacts that are associated with infectiousness and transmission [5,6]. Much progress has been made in recent years in theory and model development, but existing methods either include assumptions that do not take the full uncertainty in the evolutionary process into account [7,8], are designed for specific datasets, with fit-for-purpose code for data analysis [9][10][11], or make limiting assumptions about the relation between sampling times and infectivity [12]. An easily accessible method without these restrictions and with the flexibility to cover a wide range of infections is currently lacking, and would bring analysis of outbreak sequence data within reach of a much broader community.
The interest in easily applicable methods for sequence data analysis in outbreak settings is demonstrated by the community's widespread use of the Outbreaker package in R [8,[13][14][15]. However, the model in Outbreaker assumes that mutations occur at the time of transmission, which does not take the pathogen's in-host population dynamics into account, nor the fact that mutations occur within hosts. The publications by Didelot et al [7] and Ypma et al [11] revealed that within-host evolution is crucial to relate sequence data to transmission trees, as is illustrated in Fig 1A: there are four unobserved processes, i.e. the time between subsequent infections, the time between infection and sampling, the pathogen dynamics within hosts, and which mutations can accumulate [9,17], or absence of a clearly defined infection time [18]. To take the complete process into account, Didelot et al [7] and Numminen et al [10] took a twostep approach: first, phylogenetic trees were built, and second, these trees were used to infer transmission trees. Didelot et al [7] used the software BEAST [19,20] to make a timed phylogenetic tree, and used a Bayesian MCMC method to colour the branches such that changes in colour represent transmission events. Numminen et al [10] took the most parsimonious tree topology, and accounted for unobserved hosts by a sampling model (which is an additional complication). This two-step approach is likely to work better if the phylogenetic tree is properly resolved (unique sequences with many SNPs), but less so if there is uncertainty in the phylogenetic tree. However, also in that case construction of the phylogenetic tree is done without taking into account that only lineages in the same host can coalesce, and that these go through transmission bottlenecks during the whole outbreak. That is likely to result in incorrect branch lengths and consequently incorrect infection times.
Hall et al [12] implemented a method in BEAST for simultaneous inference of transmission and phylogenetic trees. BEAST allows for much flexibility when it comes to phylogeny and population dynamics reconstruction (for which it was originally developed [19,20]), e.g. by allowing variation in mutation rates between sites in the genome, between lineages, and in time. However, datasets of fully observed outbreaks often do not contain sufficient information for reliable inference: they typically cover only a few months up to at most several years (as in Didelot et al [7], with tuberculosis) and do not contain many SNPs (usually of the same order of magnitude as the number of samples). A more important limitation is that the transmission model implemented in BEAST is rather specific: it allows for transmission only during an infectious period informed by positive and negative samples, during which infectiousness is assumed to be constant. This may put prior constraints on the topology and order of events in the transmission and phylogenetic trees, which is undesirable if the primary aim is to reconstruct the transmission tree with little or no prior information about when hosts were infectious.
Previously, Ypma et al [11] had also developed a method for simultaneous inference of transmission and phylogenetic trees, albeit with rather specific assumptions on the within-host pathogen dynamics and the time and order of transmission events, and with no available implementation. However, their view on the phylogenetic and transmission trees was quite different. Instead of a phylogenetic tree with transmission events, they regarded it as a hierarchical tree. The top level is the transmission tree, with hosts having infected other hosts according to an epidemiological transmission model. The lower level consists of phylogenetic "mini-trees" within each host. A mini-tree describes the within-host micro-evolution. It is rooted at the infection time and has tips at transmission and sampling events; in its simplest form it is only a single branch from infection to sampling. The complete phylogenetic tree then consists of all these mini-trees, connected through the transmission tree. That description allowed them to develop new MCMC updating steps, some changing the transmission tree, some the phylogenetic mini-trees.
We built further on that principle to reconstruct the transmission trees of outbreaks, in a new model and estimation method. The method requires data on pathogen sequences and sampling times. The model includes all four underlying stochastic processes (Fig 1A), each described in a flexible and generic way, such that we avoid putting unnecessary prior constraints on the order of unobserved events (Fig 1B). This allows for application of the method to a wide range of infectious diseases, including new emerging infections where we have little quantitative information on the infection cycle. The method is implemented in R, in a package called phybreak. We illustrate the performance of the method by applying it to both new and previously published simulated datasets. We demonstrate the range of applicability by applying the model to five datasets on densely sampled infectious disease outbreaks, covering a wide range of epidemiological settings.

Outline of the method
The method infers infection times and infectors of all cases in an outbreak. The data consist of sampling times and sequences of all cases, where some of the sequences may be non-informative if no sequence is available. Using simple models for transmission, sampling, within-host dynamics and mutation, samples are taken from the posterior distributions of model parameters and transmission and phylogenetic trees, by a Markov-Chain Monte Carlo (MCMC) method. The main novelty of our method lies in the proposal steps for the phylogenetic and transmission trees that are used to generate the MCMC chain. It makes use of the hierarchical tree perspective, in which the phylogenetic tree is described as a collection of phylogenetic mini-trees (one for each host), connected through the transmission tree (see Methods for details).
The posterior samples are summarized by medians and credible intervals for parameters and infection times, and by consensus transmission trees. Consensus transmission trees are based on the posterior support for infectors of each host, defined as the proportion of posterior trees in which a particular infector infects a host. The Edmonds' consensus tree takes for each host the infector with highest support, and uses Edmonds' algorithm to resolve cycle and multiple index cases [21], whereas the Maximum Parent Credibility (MPC) tree is the one tree among the posterior trees with maximum product of supports [12].
The models and parameters used for inference are as follows: • transmission: assuming that all cases are sampled and the outbreak is over, the mean number of secondary infections must be 1. The transmission model therefore consists only of a Gamma distribution for the generation interval, i.e. the time interval between a primary and a secondary case. This transmission model is equivalent to a Kermack-McKendrick type renewal model with a generalized infectiousness function, with basic reproduction ratio R 0 = 1 in an infinite population [22]. The model contains two parameters: the shape a G , which we fixed at 3 during our analyses, and the mean m G , which is estimated and has a prior distribution with mean μ G and standard deviation σ G . In an uninformative analysis, μ G = 1, and σ G = 1.
• sampling: the sampling model consists of a Gamma distribution for the sampling interval, which is the interval between infection and sampling of a case. The model contains two parameters: the shape a S , which is fixed during the analysis, and the mean m S , which is estimated and has a prior distribution with mean μ S and standard deviation σ S . In an uninformative analysis, μ S = 1, and σ S = 1, but a S is chosen to reflect prior information (the coefficient of variation is a À 1=2 S ); in a naïve analysis we additionally set a S = 3.
• within-host dynamics: the within-host model determines the genetic relation between the tips of the within-host phylogenetic mini-tree (at sampling and transmission) through a coalescence model, assuming that samples are clonal lineages. The within-host model describes a linearly increasing pathogen population size w(τ) = rτ, at time τ since infection of a host. This within-host model results in a bottleneck at transmission of 1 lineage.The slope r has a Gamma distributed prior distribution with shape a r and rate b r . In an uninformative analysis, we used a r = b r = 3.
• The mutation model is a site-homogeneous Jukes-Cantor model, with per-site mutation rate μ. The prior distribution for log(μ) is uniform.

Analysis of newly simulated datasets
We generated 25 new simulated datasets of 50 cases with the above model, which we modified by taking a population of 86 individuals and a basic reproduction number R 0 = 1.5 (instead of an infinite population with R 0 = 1). Parameters were a G = a S = 10, m G = m S = r = 1, μ = 10 −4 and sequence length 10 4 , resulting in 1 genome-wide mutation per mean generation interval of one year. Table 1 shows some summary measures on performance of the method (see S1 Results for additional measures and results for more simulations). A 5,000 cycle burn-in followed by sampling a single chain of 25,000 MCMC cycles took about 30 minutes on a 2.6 GHz CPU (Linux). Four sets of results are shown, all with an uninformative prior for μ: one with all parameters other than μ fixed at their correct value, and three with uninformative priors for m G and r, and different levels of prior knowledge on m S : informative with correct mean,  Table 1 shows the results on tree inference. Infection times (using all MCMC samples) are well recovered if the mean sampling interval does not have a strong incorrect prior: coverage of 95% credible intervals is good, and medians may only be slightly positively biased (later than true infection time) if uninformative priors are used. For this simulation scenario, consensus transmission trees contained almost 70% (35 out of 50) correct infectors, as determined by counting infectors and resolving multiple index cases and cycles in the tree (Edmonds' method [21]) and slightly fewer when choosing the maximum parent credibility (MPC) tree [12] among the 50,000 posterior trees. Infectors with high support are more likely correct: 84% (28 out of 33) are correct if the support is above 50%, and 96% (15.2 out of 15.8) are correct if the support is above 80%. These numbers are similar in smaller outbreaks (S1 Results). If sampling and generation interval distributions are wider, the sampling times contain less information on the order of infection, which reduces the accuracy of transmission tree reconstruction (S1 Results). Using prior information on the mean sampling interval did not improve on this, but if prior information is incorrect, fewer hosts have a strongly supported infector, which makes inference more uncertain. In conclusion, the method is fast and efficient if applied to simulated data fitting the model. In that case, no informative priors are needed for transmission tree inference, though correct estimation of the infection time is aided by some information.
For comparison, we analysed the same datasets with the Outbreaker package in R [8], which uses the assumption of mutation at transmission, and with the TransPhylo package [7,24], which requires input of a phylogenetic tree that we created in BEAST v2 [19] with a constant population coalescent model and Jukes-Cantor substitution model. Both Outbreaker and TransPhylo require input of a generation and sampling interval distribution, for which we supplied the distributions used to simulate the data. Thus, the results are best compared to the results of the reference scenario of our model (Table 1, left-most column). The numbers of correctly identified infectors (Edmonds' consensus tree [21]) were smaller with both alternative methods: in the 25 outbreaks of Table 1 (50 cases, a G = a S = 10), Outbreaker identifies on average 27.5 out of 50 infectors correctly, TransPhylo 32.2, and phybreak 34.9. Also in smaller outbreaks or with different generation and sampling interval distributions, phybreak identified 8-22% more infectors correctly (S1 Results).
We also analysed the simulated results with 20% of the cases removed from the dataset, to assess performance if outbreaks are not completely observed. Table 2 shows the results with reference (parameters fixed and correct) and uninformative analyses, in comparison with the reference scenario and all data observed. With some of the cases removed, some of the remaining cases did not have their infector in the dataset anymore; these cases are referred to as orphans in Table 2. Infection time estimation was less accurate, with only 85% of credible interval containing the correct value, and more infection times estimated too early in the outbreak. Surprisingly, this was not only the case with orphans, for which this may have been expected with their infector not present in the data. It turns out that infectors are correctly identified about 20% less accurately, for all threshold levels of support. However, when correcting for presence of the infector in the data, infectors are identified with the same accuracy as in the complete dataset. We also checked how frequently the identified infector of orphans was in fact an earlier ancestor in the transmission tree, i.e. the infector's infector in most cases. It turned out that ancestors were often identified as infector, but not as accurately as the true infector identification in complete datasets ( Table 2).

Analysis of previously published simulated data
We applied the method to previously published outbreak simulations [12]. Briefly, a spatial susceptible-exposed-infectious-recovered (SEIR) model was simulated in a population of 50 farms, with a latent period (exposed) of two days and a random infectious period with mean 10 days and standard deviation 1 day, at the end of which the farm was sampled. Two mutation rates were used with an HKY substitution model, either Slow Clock or Fast Clock, equivalent to 1 or 50 genome-wide mutations per generation interval of one week, respectively. Table 3 shows performance of the method with naïve and informative prior information on the sampling interval distribution (see S1 Results for uninformative). Effective sample sizes of parameters and phylogenetic trees are a bit smaller than with analysis of the new simulations, but in most cases still good for infection times, whereas sampling of infectors was excellent. The low variance of the sampling interval distribution caused some problems in efficient sampling of m S because of its high correlation with the associated infection times, but it also caused problems in the burn-in phase if inference starts with parameter values far from their actual values (not shown). This was especially the case in the uninformative Slow Clock analyses, resulting in unreliable estimates of the mean sampling interval and infection times (S1 Results). With the Fast Clock analyses there were no such problems, as long as the full set of proposal paths in the MCMC chain was used (see Methods for details). Posterior median mutation rates are slightly higher than used for simulation, which could be due to different rates for transition and transversion in the simulation model [12]. Consensus trees with uninformative and informative settings were almost as good as in the original publication [12], in which spatial data were used and in which it was known that there was a latent period and that hosts could not transmit after sampling. In the Slow Clock simulations about 62% of infectors were correct, and in the Fast Clock simulations about 92%. Infection times were also well recovered in most cases, but not in the uninformative Slow Clock MPC, maximum parent credibility. a a S = 3, μ S = 1, σ S = 1 b a S = 144, μ S = 12, σ S = 1 c path-distance approximate ESS [23] https://doi.org/10.1371/journal.pcbi.1005495.t003 Simultaneous inference of phylogenetic and transmission trees in infectious disease outbreaks analysis (S1 Results). In the naïve analyses, the Slow Clock consensus trees were only slightly less good (but not the infection times), whereas the Fast Clock consensus trees became much worse, with only 65% of infectors correct. In conclusion, the method performs well if applied to data simulated with a very different model. Good prior information on the sampling interval can significantly improve both MCMC mixing and transmission tree inference, especially if the genetic data contain many SNPs.
The results for the four smaller datasets are shown in Table 4, which shows that mixing of the MCMC chains was generally good. The Mtb data were analysed with naïve prior information, which resulted in a median sampling interval of 419 days (similar to estimated incubation times [31]), a median generation interval of 107 days, and a mutation rate equivalent to 0.3-1.3 mutations per genome per year, as estimated before [32,33]. The Edmonds' consensus transmission tree (Fig 2A) shows low support for most infectors, which is a reflection of the low number of SNPs, but also of the long sampling interval relative to the generation interval, which makes the sampling time less Simultaneous inference of phylogenetic and transmission trees in infectious disease outbreaks informative of the order of infection. However, the same index case K02 and three clusters as identified in Didelot et al [7] are distinguished: one starting with K22, one with K35, and the remaining cases starting with K16 or infected by the index case. The main difference compared to the original analysis lies in the shape of the phylogenetic tree and the estimated infection times ( Fig 3A). Whereas the topology is very similar, the timing of the branching events is different: in the original tree, internal branches decrease in length when going from root to tips. That shape is consistent with a coalescent tree based on a single panmictic population but also reflects the fact that three mutations separate the two clades after the root node, whereas the posterior median genome-wide mutation rate is estimated at 0.48 per year (mutation rate × sequence length). By taking into account the fact that coalescent events take place within individual hosts, our analysis shows branch lengths that are more balanced in length across the  [7]; (B) MRSA data [25]; (C) FMD2001 data [26]; (D) FMD2007 data [27]. https://doi.org/10.1371/journal.pcbi.1005495.g002 Simultaneous inference of phylogenetic and transmission trees in infectious disease outbreaks tree. Importantly, this results in a more recent dating of root of the tree: early 2008 ( Fig 3A) instead of early 2004 [7]. The MRSA data were analysed with an informative prior for the mean sampling interval m S and a shape parameter a S based on data on the intervals between hospitalisation and the first positive sample. The estimated mutation rate is similar to literature estimates [34,35], but the posterior median m S of 31 days is considerably higher than the prior mean of 15 days (Table 4). This may be explained by the two health-care workers (HCW_A and HCW_B), which have very long posterior sampling intervals that were not part of the data informing the prior (Edmonds' consensus tree, Fig 2B). In contrast with the original analysis, we now identify a transmission tree rather than only a phylogenetic tree, resulting in the observation that the two health-care workers may not have infected any patient in spite of their long infection-to-sampling interval. Almost all transmission events with low support (<20%) involved unsequenced  [7]; (B) MRSA data [25]; (C) FMD2001 data [9,26]; (D) FMD2007 data [9,27]. hosts. Two of them were identified as possible infector (P5 and P7), in the initial stage of the outbreak, when only few samples were sequenced. This indicates that some unsequenced hosts may have played a role in transmission, but that it is not clear which. Finally, a major difference between our results and those in the original publication is the shape of the phylogenetic tree and the dating of the tree root: around 1st January ( Fig 3B) instead of 1st September the year before [25].
Analysis of the FMD2001 and FMD2007 datasets resulted in posterior sampling intervals with means of 14 and 10 days, respectively, close to the 8.5 days estimated from epidemic data [36]. Generation intervals were about the same (Table 4). Both datasets contained more SNPs than the Mtb and MRSA data, with unique sequences for each host and higher mutation rates, similar to published rates in FMD outbreak clusters [37]. This resulted in equal Edmonds' and MPC consensus transmission trees, and much higher support for most infectors (Figs 2C, 2D, 3C and 3D). Our transmission tree is almost identical to the one from Ypma et al [11], who used a closely related method but did not allow for transmission after sampling. When comparing to the analysis of these data by Morelli et al [9], the topologies of the phylogenetic trees (Fig 3C and 3D) match the topologies of the genetic networks (Fig S18 in [9]), but the transmission trees are quite different. The main differences are that in the FMD2001 outbreak, they identify farm B as the infector of C, E, K, L, O, and P; and in the FMD2007 outbreak, they have IP4b infecting IP3b, IP3c, IP6b, IP7, and IP8. Differences are likely the result of their use of the Simultaneous inference of phylogenetic and transmission trees in infectious disease outbreaks spatial data [9]. Use of additional data is expected to improve inference, although their estimates of infection-to-sampling intervals (about 30 days) were unrealistically long. The H7N7 dataset was analysed with the sequences of the three genes HA, NA, and PB2 separately, and combined; with informative priors for both the mean sampling and mean generation intervals. Five parallel chains were run, and mixing was generally good (Table 5); it took about 7 hours on a 2.6GHz CPU to obtain 25,000 unthinned samples in a single chain. Analysis of the three genes combined resulted in a posterior median m S of 8.4 days, slightly longer than the 7 days on which the informative prior was based [38], and longer than in the separate analyses. The mean generation time was slightly shorter than the prior mean: 4.6 days with all genes. We also calculated the parsimony scores of the posterior sampled trees, defined as the minimum numbers of mutations on the trees that can explain the sequence data [39], and compared these with the numbers of SNPs in the data (Table 5). It appeared that with the genes separately analysed, parsimony scores were 6-13% higher than the numbers of SNPs, indicating some homoplasy in the phylogenetic trees (this was not seen with any of the other datasets). The analysis of all genes together resulted in parsimony scores of 18% higher than the number of SNPs. The estimated mutation rates are among the highest estimates for Avian Influenza Virus, as already noted before in earlier analyses of the same data [28,40]. Fig 4 shows the Edmonds' consensus tree in generations of infected premises, indicating locations and inferred infection days (full details in S1 Results). Without the use of location data, there is a large Limburg cluster, a Central cluster including two sampled Limburg cases, and a small Limburg cluster of three cases with an exceptionally long generation time (asterisk in Fig 4). A closer look at the sequences makes clear that the first of these cases (L24/34) has 3 SNPs different from assigned infector G4/11, and 4 SNPs different from cases in the large Limburg cluster. Using geographic data as in earlier analyses [12,30] will probably place these cases within that cluster.

Discussion
We developed a new method to reconstruct outbreaks of infectious diseases with pathogen sequence data from all cases in an outbreak. Our aim was to have an easily accessible and widely applicable method. For ease of access, we developed efficient MCMC updating steps which we implemented in a new R package, phybreak. We tested the method on newly simulated data, previously published simulated data, and published datasets. Our model is fast: 25,000 iterations took roughly 30 minutes with the Mtb and MRSA datasets of about 30 hosts, and 7 hours with the full three-genes H7N7 dataset in 241 hosts. Two chains with 50,000 posterior samples proved sufficient (measured by ESS) for tree inference (infectors and infection times) and most model parameters with most simulated and published data. The package contains functions to enter the data, to run the MCMC chain, and to analyse the output, e.g. by making consensus trees and plotting these (as in Figs 2 and 3).
Analysis of simulated datasets showed that the sampling times play an important role in transmission tree reconstruction. Firstly, the use of prior information on the sampling interval distribution (shape parameter as well as mean) greatly improves mixing of the MCMC chain (Tables 1 and 3). Secondly, the use of (correct) prior information on the sampling interval distribution can significantly improve infection time estimation as well as transmission tree reconstruction (Table 3). Thirdly, the extent to which sampling times are correlated with infection times determines how well the method is capable of reconstructing transmission trees, which appears from the fact that outbreaks are less well reconstructed with wider sampling interval distributions (Table 1 vs S1 Results) and the low support for the posterior infectors in the Mtb analysis, where sampling intervals were much longer than generation intervals. Therefore, it is advisable to use prior information on sampling intervals in the analysis (if available), and also to base conclusions not only on the summary transmission tree, but also on the posterior support of links in that tree.
We tested the method on five published datasets, with outbreaks of viral and bacterial infections, and in diverse settings of open and closed populations, in human and veterinary context. The method performed well on all datasets in terms of MCMC chain mixing and tree reconstruction. With naive priors on mean sampling intervals and mutation rates, we obtained estimates that were all very accurate when compared to literature, and the inferred transmission trees seemed as good, or even better when considering estimated infection times. With two datasets (MRSA and H7N7) we included some prior information on sampling and/or generation intervals, which mainly affected the inferred infection times, but not so much the transmission trees. It is possible that not all cases have been observed in these outbreaks, especially in the Mtb and MRSA outbreaks, an assumption nevertheless made by our model. If not too many cases are missing, the analyses of simulations show that this does not disturb identification of infector-host pairs that are in the data. It will only limitedly affect identification of transmission clusters, because if a host's true infector is not in the data, the true infector's infector is often selected as the most likely infector. Only some of the infection times may have been estimated too early.
For wide applicability, we kept the underlying model simple without putting prior constraints on the order of unobserved events such as infection and coalescence times. Four submodels with only one or two parameters each were used for sampling, transmission, withinhost pathogen dynamics, and nucleotide substitution. The sampling model, a gamma distribution for the interval between infection and sampling, has a direct link to inferred infection times, and is the model for which it is most likely that prior information is available from epidemiological data in the same or other outbreaks. We used simulated data to study the effect of uninformative or incorrect prior information on shape parameter a S and mean m S . It appears that an incorrect a S or an incorrect informative prior for m S does reduce accuracy of inferred infection times. However, consensus trees are hardly affected, at least if the number of SNPs is in the order of the number of hosts as we saw in the actual datasets (Table 1 and  Table 2 Slow Clock). Only the precision of consensus trees is reduced, i.e. there are fewer inferred infectors with high support. Results with the Fast Clock simulations did show a significant reduction in consensus tree accuracy. In that case, there are so many SNPs that the phylogenetic tree topology and times of coalescent nodes are almost fixed; then, too much variance in sampling intervals (low a S ) results in many incorrect placements of infection events on that tree. Possibly, with so many SNPs it could be more efficient to first make the phylogenetic tree, and then use that tree to infer transmission events [7,10], but it is questionable whether genome-wide mutation rates are ever so high that this becomes a real issue [41].
The submodel for transmission is relevant for transmission tree inference in describing the times between subsequent infection events. Transmission is modelled as a homogeneous branching process, implicitly assuming that there was a small outbreak in a large population, with a reproduction number (mean number of secondary cases per primary case) of 1. If all, or almost all, infectors are in the data, the generation interval distribution reflects the course of infectiousness, separating the cases in time along the tree. This interpretation may be obscured with many unobserved cases, as in the absence of the actual infector, the method often identifies an earlier ancestor in the transmission tree as infector ( Table 2). Apart from not taking heterogeneity across hosts into account (an extension we wish to leave for future development, see below), the current model also neglects the possibility that susceptibles can have contact with several infecteds in a smaller population or more structured contact network. That could be modelled by a force of infection, which would more realistically describe contraction of the generation interval during the peak of the outbreak, and provide estimates for relevant quantities such as reproduction ratios [6]. However, it requires information about uninfected susceptibles in the same population and a more complicated transmission model, which is a significant disadvantage when it comes to general applicability, one of our primary aims. More importantly, for transmission tree inference it does not seem to be a problem: the analyses of the published simulations were almost as accurate as in the original publication [12], and these simulations were in very small populations with almost all hosts infected.
The role of the within-host model is to integrate over all possible phylogenetic mini-trees and mutation events within the hosts, and through that, to obtain a posterior distribution of all transmission trees consistent with the (genetic) data. For this, we used a coalescence model based on a linearly growing within-host population, combined with a Jukes-Cantor substitution model. These models contain each only one parameter, but we think that-as long as only few mutations occur in each host, as in our own simulations, the published Slow Clock simulations, and most datasets-for most applications more complex models are not needed for the following reasons. First, the gross structure of the phylogenetic tree topology and branch lengths result from transmission and sampling models, and only the finer within-host details are determined by the within-host model. With only few mutations within each host, precise mini-tree inference is not possible, and for our aim of inferring transmission trees, unnecessary. Second, and confirming this imprecise mini-tree inference, most tree proposal steps include simulation of the within-host phylogenetic mini-trees, resulting in good mixing of transmission and phylogenetic tree topologies. The fact that proposing from the prior distribution works so well indicates that the sequence data do not contain much information on within-host branch lengths. Third, if there are few SNPs, the posterior contains almost only phylogenetic trees with fewest mutations (maximum parsimony). It is therefore not likely that tree inference will improve with more general substitution models. Fourth, inference of transmission trees and infection times appears not to be biased if the underlying simulation model was more realistic ( Table 2). If data do contain many SNPs, as in the Fast Clock simulations, more detailed and realistic models for within-host pathogen growth and nucleotide substitution do probably improve inference, especially on the phylogenetic tree. Nonetheless, even then our method was still capable of correctly inferring the infection times and transmission trees with almost the same accuracy as in the original publication.
With two exceptions, the parsimony scores of posterior tree samples were always equal to the number of SNPs in the datasets (the minimum possible). The first exception is the set of Fast Clock published simulations, which had so many SNPs that many of the same mutations had occurred in parallel. The second exception is the H7N7 dataset. In that case, the analyses of the three genes separately resulted in parsimony scores with 6-12 (6%-13%) more mutations than the number of SNPs, whereas the analysis of all genes together resulted in a parsimony score of 313 (median) to explain only 257 SNPs, a surplus of 56 mutations (18%). The results for separate genes could indicate positive selection, confirming the analysis by Bataille et al [28], who even identified specific mutations that had occurred multiple times. The even higher discrepancy for the combined analysis is suggestive of reassortment events, also recognised by Bataille et al [28].
The proposed method and implementation opens perspectives for further extending the methodology to reconstruct phylogenetic and transmission trees from pathogen sequence data. One possible set of extensions arises from changes to the models embedded in our method, to include additional aspects of outbreak dynamics. For instance, the generation time distribution (infectiousness curve) could be made dependent on the sampling interval, which may be relevant for the MRSA outbreak analysis in which the two health-care workers may have transmitted the bacterium until late after infection. This dependence is implicit in methods in which transmission is modelled more mechanistically (e.g. [11,12,16]), but we chose not to do that to keep the model more generic. Another important extension would be to relax the assumption of a complete bottleneck at transmission; the bottleneck may be larger in reality [42,43] and it has previously been relaxed by looking at transmission pairs [44] or modelling it as separate transmission events [18], but not yet in a timed transmission tree. In our model, this would mean that a host can carry multiple phylogenetic mini-trees, rooted at the same infection time to the same infector. A third extension would be to include the possibility of reassortment of genes within a host, primarily motivated by the results of the H7N7 analysis. This may be done by modelling the coalescent process within hosts, the phylogenetic minitrees, differently for different genes, but constrained by a single transmission tree. Finally, it would be possible to allow for multiple index cases, which may play a role in open populations with possible re-introductions (as in the MRSA setting), or when only a subset of a large epidemic is analysed (the FMD2001 dataset). This is implemented in models using genetic models based on pairwise genetic distances [8,16] and with a model assuming coalescence at transmission [45], but is considered a major challenge with a within-host coalescent model [46]. Multiple index cases could also reflect unobserved hosts in the outbreak itself, recently addressed by Didelot et al [24] in their two-step approach of first inferring a phylogenetic and then a transmission tree.
A second type of extension would stem from incorporating additional data. An example is the use of data that make particular transmissions more or less likely, such as contact tracing data, or censoring times for infection times per host or transmission times between sets of hosts, motivated by the MRSA dataset in which admission and discharge days are known for each patient. Sampling of infection times and infectors could be constrained by these additional data (as in [12,30]) and could then become less dependent on the sampling times and sampling interval distribution, as in the current implementation. Another example is the use of spatial data in combination with a spatial transmission kernel, so that the likelihood of infectors includes a distance-dependency, a possible extension motivated by the FMD and H7N7 analyses (as in [9,30]). A third example is the use of host characteristics to model infectivity as a function of covariates. With the MRSA data, it would then be possible to test for increased infectivity of the health-care workers, or to test for differences in transmissibility in the three wards. In general, the use of additional host data would make dealing with hosts for which a sequence is not available less problematic: the method currently can include these hosts, but without additional data their role is unclear and they are often placed at the end of transmission chains in consensus trees (Fig 2B, Fig 3).

Methods Data
We developed our model for fully observed outbreaks of size n hosts. Data consist of the sampling times S and DNA sequences G, which means that for each host i we know the time of sampling or diagnosis S i and the sequence G i associated with the sampling time. Hosts without available sequence are given a sequence with noninformative nucleotides (only 'n').
We illustrate the method with the following five datasets from earlier publications (all in S1 Data): 1. Tuberculosis (Mycobacterium tuberculosis, Mtb). This dataset was analysed by Didelot et al [7]. It consists of 33 Mtb cases in a population of drug users (approximate population size 400), with samples collected in a 2.5 years time frame. The 4.4 Mbp long sequences contained 20 SNPs. Analysis of this dataset tests the performance of this method in an outbreak with relatively few cases in a large population.
2. Methicillin-resistant Staphylococcus aureus (MRSA). This is the dataset from Nubel et al [25], with 36 MRSA cases in a neonatal ICU sampled within a time period of 7 months. Sampling dates were available for all cases, but sequences only for 28 cases, revealing 26 SNPs in the non-repetitive core genome of 2.7 Mbp. Analysis of this dataset tests for the performance of this method in an outbreak in a small population, including cases without sequence.
3. Foot-and-mouth disease (FMD2001). This is the dataset from Cottam et al [26] also analysed by several others [9,11], with 15 infected premises within a time period of 2 months. Sequences were available for all cases, with 85 SNPs among 8196 nucleotides. Analysis of this dataset and the next tests for the performance of this method in a small completely sampled outbreak in a large population and allows comparison of the estimated transmission tree to earlier results.
4. Foot-and-mouth disease (FMD2007). This is the dataset from Cottam et al [27], also analysed by Morelli [9], with 11 infected premises within a time period of 2 months. Sequences were available for all cases, with 27 SNPs among 8176 nucleotides 5. H7N7 avian influenza (H7N7). This dataset has been analysed by several authors [12,[28][29][30], and consists of 241 poultry farms in a time period of about 2.5 months. Sequences of the HA, NA, and PB2 genes were available on GISAID for 228 farms, with associated sampling dates. The total number of SNPs was 257 in 5541 nucleotides. For the 13 unsampled farms we used the culling date minus 2 days as the observation day (the mean sampling-toculling interval was 2.4 days in the 228 sampled farms). We analysed the data for the three genes separately, and combined (concatenated). To inform a prior distribution for the interval from infection to sampling, we used estimated infection times from Boender et al [38]. Analysis of this dataset tests for the performance of this method in a large outbreak, including cases without sequence.

The model and likelihood
The model describes the spread of an infectious pathogen in a population through contact transmission, the dynamics of the pathogen within the infected hosts, and mutation in the DNA or RNA of that pathogen. Furthermore, it describes how these dynamics are observed through sampling of pathogens in infected hosts. We infer the transmission tree and parameters describing the relevant processes by a Bayesian analysis, using Markov-Chain Monte Carlo (MCMC) to obtain samples from the posterior distributions of model parameters and transmission trees (infectors and infection times of all cases). We first introduce the models and likelihood functions; then we explain how we update the transmission trees and parameters in the MCMC chain. The posterior distribution is given by We now introduce the four models, the associated likelihoods, and prior distributions for associated parameters.
Transmission. We assume that the outbreak started with a single case. Each case produced secondary cases at random generation intervals after their own infection (Gamma distribution with shape a G and mean m G ). We consider that all untimed transmission tree topologies are equally likely, so that the probability of the transmission tree only depends on its timing. The outbreak is described by the vectors I and M with infection times I i and infectors M i for all numbered cases i. The infector of the index case is 0. The likelihood is the product of probability densities (d Gða G ;m G Þ ðÁÞ) of all generation times in the outbreak: Sampling. We assume that all cases are observed and sampled once at random times after they were infected, according to a Gamma distribution with shape a S and mean m S . Transmission and sampling are independent, so transmission can take place after sampling, and a case can be sampled earlier than its infector. The likelihood is the product of probability densities of all sampling intervals in the outbreak: Within-host dynamics. The main function of the within-host model is to allow for a stochastic coalescent process within the host. Each host i harbours its own phylogenetic mini-tree P i , with the tips being the transmission and sampling events, and the root being the time of infection, before the first coalescent node. Thus, samples are assumed to be clonal lineages. In its simplest form, a phylogenetic mini-tree consists of a single branch from infection time to sampling time. The likelihood is the product of all likelihoods per host: in which r is the parameter describing the within-host dynamics (see below). The likelihoods per mini-tree are dependent on all infection times and infectors, because these determine the transmission times with host i as infector. Going backwards in time, coalescence between any pair of lineages within a host takes place at rate 1/w(τ,r), where w(τ,r) = rτ denotes the linearly increasing within-host pathogen population size at (forward) time τ since infection of the host. With this particular function coalescent nodes tend to be close to the transmission events if r is small, whereas they tend to be soon after infection of the infector if r is large. This function also naturally results in only one lineage at the time of infection (complete transmission bottleneck), as the coalescence rate goes to infinity near the time of infection.
In the complete phylogenetic tree P, three types of nodes x are distinguished (Fig 1A): nodes x = 1. . .n are the sampling nodes of the corresponding hosts i = 1. . .n, i.e. the tips of the tree at which sampling took place; nodes x = n+1. . .2n-1 are the coalescent nodes; nodes x = 2n. . .3n-1 are the transmission nodes, i.e. the points in the tree at which a lineage goes from one host to the next. By h x we identify the host in which node x resides; for transmission nodes it identifies the primary host (infector). The mini-tree P i is the set of nodes within host i, and τ x is the time of node x since infection of host h x , so τ i is the time of sampling. Let L i (τ) denote the number of lineages in host i at time τ since infection: 0. The first term is the probability to have no coalescent events during the intervals in which there are two or more lineages, the second term is the product of coalescent rates at the coalescent nodes. Mutation. We use a single fixed mutation rate μ for all sites, with mutation resulting in any of the four nucleotides with equal probability (Jukes-Cantor). This parameterisation means that the effective rate of nucleotide change is 0.75μ. Given the phylogenetic tree, this results in the likelihood: Here, we multiply over all coalescent and transmission nodes x, which occur at time t x and have parent node v x ; I mut indicates if a mutation occurred on the branch between x and ν x , and N indicates if a branch ends with a tip with noninformative nucleotide ('n' in the sequence). The likelihood is calculated using Felsenstein's pruning algorithm [47]. Prior distributions. Here we describe our general choice of prior distributions, not the particular parameterization in our analyses (Section Evaluating the method). We chose fixed values for a G and a S , the shape parameters of generation and sampling intervals. For their means m G and m S , we used prior distributions with means μ G and μ S and standard deviations σ G and σ S , which are translated into Gamma-distributed priors for rate parameters b G = a G /m G and b S = a S /m S , distributed as Γ(a 0,G ,b 0,G ) and Γ(a 0,S ,b 0,S ) (see S1 Methods). For the slope r of the within-host growth model, we chose a Gamma-distributed prior with shape and rate a r and b r . We chose log(μ) to have a uniform (improper) prior distribution, equivalent to Pr(μ) / 1/μ.

Inference method
We use Bayesian statistics to infer transmission trees and estimate the model parameters from the data, and MCMC methods to obtain samples from the posterior distribution. The procedure is implemented as a package in R (phybreak), which can be downloaded from GitHub (www.github.com/donkeyshot/phybreak) and is available on CRAN (cran.r-project.org/web/ packages/phybreak/index.html). The package also contains functions to simulate data, and to summarize the MCMC output.
The main novelty of our method lies in the proposal steps for the phylogenetic and transmission trees, used to generate the MCMC chain. It makes use of the hierarchical tree perspective, in which the phylogenetic tree is described as a collection of phylogenetic mini-trees (one for each host), connected through the transmission tree. Most proposals are done by taking one host, changing its position in the transmission tree, and simulating the phylogenetic minitrees in the hosts involved in that change. In a second type of proposal, the transmission tree is changed while keeping the phylogenetic tree intact. A third type of proposal only resimulates the within-host mini-tree topology, keeping the transmission tree and coalescent times intact.
Initialization of the MCMC chain requires initial values for the six model parameters (a G , m G , a S , m S , r, and μ). The transmission tree is initialized by generating an infection time for each host (sampling day minus random sampling interval). The first infected host is the index case, and for the remaining hosts an infector is randomly chosen, weighed by the density of the generation time distribution. Finally, the phylogenetic mini-trees in each host are simulated according to the coalescent model and combined with one another to create a complete phylogenetic tree.
Each MCMC iteration cycle starts with updates of the transmission and phylogenetic trees, followed by updates of the model parameters. To start with the latter, the parameters m S and m G are directly sampled from their posterior distribution given the current infection times and transmission tree (Gibbs update). This is done by sampling the rate parameters b S and b G , which were given conjugate prior distributions (see above). If T S ¼ X S i À I i is the sum of n sampling intervals in the tree, a 0,S and b 0,S are the shape and rate of the prior distribution for b S , then a new posterior value is drawn as Per MCMC iteration cycle, n proposals are done with each host as a focal host once, in random order. Each proposal starts by taking a focal host i, drawing a sampling interval T $ G 2 3 a S ; m S À Á from a Gamma distribution with shape parameter 2 3 a S and mean m S , and calculating a preliminary proposal for the infection time I i ' = S i −T. The 2 3 is chosen to have a slightly higher variance for the proposal than in the likelihood, to improve mixing. Based on this preliminary proposal, the topology of the transmission tree is changed (see below), and in most cases the phylogenetic tree as well (80% probability). However, we also allowed for proposal steps without changing the phylogenetic tree (20% probability); this greatly improves mixing of the MCMC chain if there are many SNPs, which more or less fixes the phylogenetic tree topology. The 80%-20% distribution for the two types of proposal was not optimized but chosen such that mixing of the phylogenetic tree is only limitedly less efficient than without the second type of proposal (keeping the phylogenetic tree fixed). It is possible to include a third type of proposal, in which only the topology of the phylogenetic mini-tree of focal host i is resimulated, to improve phylogenetic tree mixing if there are many SNPs. In the default setting in the phybreak package, 80% of proposals are of the first type and 20% of the second type, but the user is free to change these percentages. We used these settings in all our analyses, expect for the Fast clock simulations, where we used a 75%-20%-5% distribution.
Proposals for changes in transmission and phylogenetic trees. Here we describe how changes in the transmission and phylogenetic trees are proposed for six different situations, based on the preliminary proposal for the infection time I i ' and on whether the index case is involved. Fig 5 shows the proposed changes. More detail on the proposal distribution and calculation of acceptance probability is given in the S1 Methods.
A. The focal host i is index case, and the preliminary I i ' is before the first transmission event.
In that case, the infection time of host i becomes I i ', and no topological changes are made in the transmission tree (Fig 5A).
B. The focal host i is index case, and the preliminary I i ' is after the first transmission event, but before host i's second transmission event, if there is any. In that case, the infection time of host i becomes I i ', and host i's first infectee becomes index case, transmitting to i (Fig 5B).
C. The focal host i is index case, and the preliminary I i ' is after host i's second transmission event, if there is any. In that case, the infection times of host i and its first infectee are switched, and host i's first infectee becomes index case. They may or may not exchange infectees, with 50% probability (Fig 5C).
D. The focal host i is not index case, and the preliminary I i ' is before infection of the index case. In that case, the infection time of host i becomes I i ', and host i becomes index case, transmitting to the original index case (Fig 5D).
E. The focal host i is not index case, and the preliminary I i ' is after infection of the index case, but before host i's first transmission event. In that case, the infection time of host i becomes I i ', and a new infector is proposed from all hosts infected before I i ' (Fig 5E).
F. The focal host i is not index case, and the preliminary I i ' is after host i's first transmission event. In that case, the infection times of host i and its first infectee are switched, as well as their position in the transmission tree. They may or may not exchange infectees, with 50% probability (Fig 5F).
Each change in the transmission tree is followed by proposing new phylogenetic mini-trees for all hosts involved, i.e. if their infection time was changed or transmission nodes were added or removed (grey hosts in Fig 5).
Proposals for changes in the transmission tree only. Here we describe how changes in the transmission tree are proposed without changing the phylogenetic tree, based on the preliminary I i ' and on whether the index case is involved. Fig 6 shows the proposed changes. More detail on the proposal distribution and calculation of acceptance probability is given in the S1 Methods.
G. The focal host i is the index case. If the preliminary I i ' is before the first coalescence node, the infection time of host i becomes I i ', and no changes are made in the transmission and phylogenetic trees. If the preliminary I i 'is after the first coalescence node, the proposal is rejected.
H. The focal host i is not the index case, and the preliminary I i 'is after the most recent common ancestor (MRCA) of the samples of host i and his infector j, which is a coalescent node in infector j. In that case, the infection time of host i becomes I i ', and infectees may move from host i to infector j or vice versa (Fig 6A).
I. The focal host i is not the index case, but his infector j is the index case, and the preliminary I i ' is before the MRCA of the samples of host i and his infector j. In that case, an infection time I j ' is proposed for the infector j. If I j ' is after the MRCA, the infection time of the infector j becomes I j ', and the infection time of host i becomes the original infection time of his infector j. Infectees may move from host i to infector j or vice versa (Fig 6B). If I j ' is before the MRCA, the proposal is rejected. In panels A, B, D, and E, the initial situation is at the top, and the proposal below. In panels C and F, the initial situation is in the middle, and two alternative proposal above and below. Every panel J. The focal host i is not the index case, and neither is his infector j, and the preliminary I i ' is before the MRCA of the samples of host i and his infector j, but after the MRCA of the samples of host i and infector j's infector. In that case, an infection time I j ' is proposed for the infector j. If I j ' is after the MRCA, the infection time of the infector j becomes I j ', and the infection time of host i becomes I i '. Infectees may move between host i, infector j, and infector j's infector (Fig 6C). If I j ' is before the MRCA of host i and infector j, or I i ' is before the MRCA of host i and infector j's infector, the proposal is rejected.

Proposal for changing the phylogenetic mini-tree only
In a single proposal path K, only a new topology of the phylogenetic mini-tree of focal host i is proposed; the coalescent times are kept unchanged.

Evaluating the method
We took three approaches to evaluate our method: analysis of newly simulated data, analysis of published simulated data [12], and analysis of published observed data. When not specified, the following parameter settings and priors were used: shape parameters for sampling and generation interval distributions a S = a G = 3, uninformative priors for mean sampling and generation intervals with μ S = μ G = 1 and σ S = σ G = 1, and a prior for within-host growth parameter r with a r = b r = 3. The prior for log(μ) (mutation rate) is always uniform. Analyses were done by two MCMC chains, in each taking 25,000 samples (25,000 MCMC cycles). Burn-ins were different: 5000 MCMC cycles for the newly simulated data, 25,000 for the published simulated data [12], and 5000 for the observed data. With the H7N7 data, five MCMC chains were run, with a burn-in of 5000 samples, followed by 25,000 samples. Burn-in lengths of simulated data were based on visual inspection of convergence for two datasets, and then choosing a burn-in period at least 10 times longer than necessary for all the other simulations, followed by comparing ESS and infector sampling in the parallel chains. The Slow clock published simulations had not all converged in the uninformative analysis (S1 Results). For the published data, all chains were inspected visually to confirm convergence.
Analysis of newly simulated data. Four outbreak scenarios were simulated, each replicated 25 times: outbreak sizes of 20 and 50 cases, each with a G = a S = 3, resulting in overlapping generations and cases sampled earlier than their infector, or a G = a S = 10, resulting in more discrete generations and cases mostly sampled in order of infection. Further, the mean generation and sampling intervals were m G = m S = 1 year, the mutation rate μ = 10 −4 per year in a DNA sequence with 10 4 sites resulting in a genome-wide mutation rate of 1 per year and a number of SNPs in the same order of magnitude as the outbreak size. For the within-host model we used r = 1 per year.
The transmission trees were simulated assuming populations of size 35 or 86 individuals and R 0 = 1.5, corresponding to expected final outbreak sizes of about 20 and 50 [22], respectively. Simulations started with one infected individual. All individuals were assumed to be equally infectious, resulting in a Poisson-distributed number of contacts at times since shows an outbreak with four hosts, with red arrows indicating transmission: the purple host is the focal host, with the purple arrow indicating the proposal for the new infection time I i '; filled hosts have a new phylogenetic mini-tree proposed; greyed-out hosts do not play a role in the proposal. (A) the focal host is the index case, and I i ' is before the first transmission event; (B) the focal host is the index case, and I i ' is after the first, but before the second secondary case; (C) the focal host is the index case and I i ' is after his second secondary case; (D) the focal host is not the index case and I i ' is before infection of the index case; (E) the focal host is not the index case and I i ' is before his first secondary case; (F) the focal host is not the index case and I i ' is after his first secondary case.
https://doi.org/10.1371/journal.pcbi.1005495.g005 Given the infection times, sampling times were drawn, and phylogenetic mini-trees were simulated for each host. These were combined into one phylogenetic tree on which random mutation events were placed according to a Poisson process with rate 1. Each mutation event was randomly assigned to one site, and generated one of the four nucleotides with equal probabilities (reducing the effective mutation rate by 25%). By giving the root an arbitrary sequence, the sampled sequences were obtained by following the paths from root to sample and changing the nucleotides at the mutation events.
The simulated data (sampling times and sequences) were analysed with four sets of parameter settings: • Reference: a G = 3, all other parameters at simulation value (except for μ); • Informative Correct: a S at simulation value, informative prior for m S with μ S = 1 and σ S = 0.1; • Uninformative: a S at simulation value; • Informative Wrong: a S at simulation value, informative prior for m S with μ S = 2 and σ S = 0.1.
In addition, the data were analysed with the Outbreaker [8] and TransPhylo [7,24] packages in R, with the correct generation and sampling intervals, in Outbreaker discretized in steps of 0.1 year. Outbreaker analyses consisted of a 200,000 iterations burn-in, followed by 2000 samples with thinning interval of 500. TransPhylo analyses consisted of 500,000 iterations burn-in, followed by 500,000 samples without thinning; the required phylogenetic trees were maximum clade credibility (mcc) trees obtained in BEAST v2 [19], assuming a constant coalescent population model. Both the Outbreaker and TransPhylo outputs were summarized by Edmonds' consensus transmission trees (see below).
Analysis of published simulated data. We used two sets of 25 simulated outbreaks, identified as Fast clock and Slow clock in the original paper [12], in which full details on the simulations can be found. Briefly summarizing some characteristics, 50 hosts were placed on a grid and a spatial transmission model was run, with exponential transmission kernel. Outbreaks with fewer than 45 cases were discarded. An SEIR (susceptible-exposed-infectious-removed) transmission model was used, with fixed latent period of 2 days and normally distributed infectious period (mean(sd) of 10(1) days). Sampling occurred at the time of removal. Phylogenetic mini-trees were simulated using a logistic within-host growth model w(τ) = 0.1(1 + e 6 )/ (1 + e 6-1.5τ ), starting at w(0) = 0.1, then growing to w(4) = 20.2 and going to w(1) = 40.4. Sequences were generated with a 14,000 base pair genome and a mutation rate of 10 −5 per site per day (Slow clock) or 5Á10 −4 per site per day (Fast clock), using a Hasegawa-Kishino-Yano (HKY) substitution model. The Slow Clock resulted in a mean number of mutations of 0.14 per part of an outbreak, with red arrows indicating transmission to depicted or undepicted hosts. Only in panel B host I must be the index case. The purple host is the focal host, with the dark purple arrow indicating the proposal for the new infection time I i '; the light purple arrow in panels B and C indicate the proposal for the new infection time I j ' of the focal host's infector. The grey parts of the phylogenetic tree are moved between the hosts. (H) the focal host is not the index case, and I i ' is after MRCA I,II of the focal host and his infector; (I) the focal host is not the index case, and I i ' is before MRCA I,II of the focal host and his infector (the index case), and I j ' is after MRCA I,II ; (J) the focal host is not the index case, and I i ' is before MRCA II,III of the focal host and his infector, but after the MRCA I,III of the focal host and his infector's infector; also, I j ' is after MRCA II,III .
https://doi.org/10.1371/journal.pcbi.1005495.g006 day, or 0.98 per mean generation time of 7 days (latent period plus half infectious period), equivalent to the rate used in the new simulations; the Fast Clock was 50 times as fast.
The simulated data (sampling times and sequences, not locations and removal times) were analysed with three levels of prior knowledge on the sampling interval distribution: • Naive: default settings; • Uninformative: a S = 144 (coefficient of variation of 0.083, as in the simulation); • Informative: a S = 144, an informative prior for m S (μ S = 12, σ S = 1).
Analysis of published datasets. The published Mtb, FMD2001, and FMD2007 datasets were analysed with default settings. The MRSA data contained information on times between hospital entry and first positive sample for 32 patients. Because of their mean and standard deviation of 20 days, we analysed these data with different prior information on the sampling interval only: a S = 1, μ S = 15, σ S = 5. For the H7N7 outbreak data, infection times of the flocks had been estimated [38], from which the mean and standard deviation of the sampling interval was calculated (7.0 and 2.2 days). We used this to inform the sampling intervals with: a S = 10, μ S = 7, σ S = 0.5. Because transmission after culling is not possible, we also used a weak informative prior for the mean generation interval: μ G = 5, σ G = 1.
Performance and outcome measures. The aim of the method is to reconstruct outbreaks in terms of infection times of all hosts and the transmission tree. This requires good mixing of the MCMC chain, especially of infection times and infectors, and a useful method to summarize all sampled transmission trees into a consensus tree.
To test for good mixing, we used effective sample sizes (ESS, calculated with the coda package in R) to evaluate mixing of the parameters and infection times. There are no strict thresholds, but in BEAST, an ESS < 100 is considered too low, whereas an ESS > 200 is considered sufficient [48]. Phylogenetic tree topology mixing was evaluated by the approximate topological ESS [23], available through the rwty package in R. Mixing of the transmission tree topology (infector per host) was evaluated as follows. To test for 200 independently sampled infectors per host, the chains were thinned by 250, giving 100 sampled infectors per chain. Then two Fisher's exact tests were done for each host. The first test was to compare the posterior frequency distributions of infectors between the two chains, with a two-row contingency table, entry {i, j} counting how often infector j occurred in chain i (100 infectors per chain in total). The second test was to assess independency of subsequent samples within the chains, i.e. absence of autocorrelation, with a contingency table in which entry {i, j} counts how often infector i was followed by infector j in the chains (198 pairs of infectors). We used the proportion of successful tests (i.e. P > 0.05) as a measure of mixing, expecting 95% successful tests with good mixing.
Two methods were used to make consensus transmission tree topologies (who infected whom), both based on the frequencies of infectors for each host among all posterior trees. The support of host j being the infector of host i is defined as the proportion of posterior trees in which host i infected host j. The first consensus tree is the maximum parent credibility (MPC) tree [12], which is the tree among all posterior trees that has the highest product of infector supports. The second consensus tree is the tree constructed using an adaptation of Edmonds' algorithm, which starts by taking the infector with highest support for each host, and resolves cycles if there are any [21]. Whereas in the original algorithm the sum of weights between nodes is minimized conditional on the absence of cycles, we maximize the sum of supports. Because the algorithm requires prior choice of an index case, we adapted it by repeating the algorithm for all supported index cases, and selecting the tree with highest sum of posterior supports.
Posterior infection times were summarized either outside the context of a consensus tree, i.e. based on all MCMC samples, or for a particular consensus tree, i.e. for each host based only on those samples in which the infector was the consensus infector. The latter is to improve consistency between topology and infection times, although even then consistency is not guaranteed. For plotting transmission trees only, we used the Edmonds' consensus tree; for plotting transmission and phylogenetic trees together, we used the MPC consensus tree, which comes with a consistent phylogenetic tree because it is one of the sampled trees.