- Record: found
- Abstract: found
- Article: found

Author(s):
A. O'Hare ^{1} ,
R. J. Orton ^{1} ,
P. R. Bessell ^{1} ^{,} ^{2} ,
R. R. Kao ^{1}

Publication date PMC-release: 22 May 2014

Journal: Proceedings of the Royal Society B: Biological Sciences

Publisher: The Royal Society

Keywords: partial likelihood, ergodic, bootstrap, nonlinear dynamics

Fitting models with Bayesian likelihood-based parameter inference is becoming increasingly
important in infectious disease epidemiology. Detailed datasets present the opportunity
to identify subsets of these data that capture important characteristics of the underlying
epidemiology. One such dataset describes the epidemic of bovine tuberculosis (bTB)
in British cattle, which is also an important exemplar of a disease with a wildlife
reservoir (the Eurasian badger). Here, we evaluate a set of nested dynamic models
of bTB transmission, including individual- and herd-level transmission heterogeneity
and assuming minimal prior knowledge of the transmission and diagnostic test parameters.
We performed a likelihood-based bootstrapping operation on the model to infer parameters
based only on the recorded numbers of cattle testing positive for bTB at the start
of each herd outbreak considering high- and low-risk areas separately. Models without
herd heterogeneity are preferred in both areas though there is some evidence for super-spreading
cattle. Similar to previous studies, we found low test sensitivities and high within-herd
basic reproduction numbers (
*R*
_{0}), suggesting that there may be many unobserved infections in cattle, even though
the current testing regime is sufficient to control within-herd epidemics in most
cases. Compared with other, more data-heavy approaches, the summary data used in our
approach are easily collected, making our approach attractive for other systems.

- Record: found
- Abstract: found
- Article: not found

Tina Toni, David Welch, Natalja Strelkowa … (2009)

Approximate Bayesian computation (ABC) methods can be used to evaluate posterior distributions without having to calculate likelihoods. In this paper, we discuss and apply an ABC method based on sequential Monte Carlo (SMC) to estimate parameters of dynamical models. We show that ABC SMC provides information about the inferability of parameters and model sensitivity to changes in parameters, and tends to perform better than other ABC approaches. The algorithm is applied to several well-known biological systems, for which parameters and their credible intervals are inferred. Moreover, we develop ABC SMC as a tool for model selection; given a range of different mathematical descriptions, ABC SMC is able to choose the best model using the standard Bayesian model selection apparatus.

- Record: found
- Abstract: found
- Article: found

Roman Biek, Anthony O'Hare, David Wright … (2012)

Introduction The application of whole genome sequencing (WGS) technology to infectious bacterial diseases has resulted in unprecedented advances in our ability to resolve epidemic data at the global scale [1], [2], provided new insights into within-host replication processes [3], and been used to corroborate the importance of exhaustively identified transmission chains and social drivers of transmission [4], [5]. However, evidence of its value when observing fine scale transmission dynamics in a partially sampled population is thus far limited to some tantalising observations [1] with as yet no quantitative evaluation of the underlying contact processes using nonlinear mathematical models. Such evaluations are particularly important where the sampling is biased, such as when the epidemiology involves an unobserved ‘reservoir’ host. Interpretation of sequence data at this scale is further complicated by the lag between the time of transmission and the time of sample collection, which can be considerable, especially for pathogens with extended latent periods. Comparing genetic data with mathematical models based on epidemiological contact data should allow us to develop more general transmission principles, and to compare and contrast the types of information that these different data sources provide. Combining mathematical models and pathogen sequence data has been an area of increased attention in the epidemiology of fast-evolving RNA viruses [6]–[9], but is as yet largely unexplored for bacteria, particularly for TB and other slow growing mycobacteria, where transmission intervals and routes tend to be more uncertain and evolutionary rates are slower, warranting more extensive sequence information. While this presents a unique set of challenges, WGS now offers promising research solutions to this problem. Mycobacterium bovis is the causative agent of bovine tuberculosis (bTB), an important disease of both livestock and wildlife with severe socio-economic consequences and impacts on animal health. Historically, it is believed to have been a major contributor to human TB cases worldwide, and it remains a zoonotic concern in both developed and developing countries [10], [11]. While most countries that employ well developed test and slaughter programs have eliminated bTB from their livestock, the control of M. bovis has proven problematic in Britain and Ireland, with the Eurasian badger (Meles meles) implicated as an important wildlife reservoir of M. bovis [12]. However, despite considerable research efforts, the role of badgers in the transmission of M. bovis remains controversial both on scientific and socio-political grounds [13], [14]. Genotyping of M. bovis from cattle and badgers based on spoligotyping and Variable Number Tandem Repeat (VNTR) typing has provided considerable insight into the epidemiology of M. bovis. In particular, the spatial clustering of genotypes in isolates from British and Irish cattle is indicative of a locally driven transmission process [15], [16]. However, neither marker has the resolution to identify fine scale transmission patterns such as occurs at the individual herd level. Links between cattle and badgers have been identified via analysis of genotype associations and statistical analyses of observed outbreak data [17], however direct evidence of transmission chains linking the two hosts at a local farm scale remains lacking. Here, we exploit the resolution of WGS to address these questions, using samples from badgers and cattle collected from a group of neighbouring farms in Northern Ireland (NI) with a decade-long history of repeated bTB outbreaks. The epidemic of bTB in Northern Irish cattle is particularly well described: annual ‘tuberculin’ testing of the entire cattle population creates a uniform sampling framework, and the majority of M. bovis isolates are now genotyped [16]. Between 2003 and 2010, this amounted to 10596 isolates from cattle, which either had a positive tuberculin test (a “reactor”) or were identified by post mortem testing of non-reactors. These data revealed 193 VNTR types within NI, with three types accounting for over 50% of sampled bacteria. M. bovis isolates from badgers are encountered at a lower frequency via an ongoing road traffic accident (RTA) survey, but their VNTR types show strong spatial associations on a regional scale with those found in cattle (R. Skuce, unpublished data). Complementing this extensive information on the pathogen are detailed demographic and network data: cattle locations and movements are recorded on an individual, daily basis, and retained on a centrally held database that can be cross-referenced with the M. bovis sample data. This combination of bacterial testing and cattle life history data provides an ideal test bed for analysing the transmission patterns reflected in WGS data. In this study, we used whole genome sequences of M. bovis and mathematical modelling to analyse the transmission between and within cattle farms and the potential role of badgers in this system. A spatial cluster of five farms with recently recorded bTB breakdowns (i.e. herds with one or more reactors, or identified abattoir cases) due to VNTR type 10 was identified. This genotype is a single locus variant of VNTR type 1 (the second most common type, representing 19.2% of identified genotypes from 2003 to 2008) that is as yet of low prevalence. Type 10's relatively recent emergence means that an identifiable common source is more likely than with more broadly distributed, longer established genotypes. VNTR type 10 breakdowns had a median duration of eight months and a median of four cases each. Median breakdown size across Northern Ireland during the study period was two cases, with a median duration of seven months. Our sequence data are derived from M. bovis isolates from 26 cattle from the years 1999 to 2010, and 5 isolates from 4 road-kill badgers (including two from the same badger) collected from within the farm cluster between 2004 and 2007 (Figure 1). 10.1371/journal.ppat.1003008.g001 Figure 1 Main holdings associated with herds in the dataset and badger locations by year. Herd locations indicate centroids of main holdings, and do not include isolated land parcels or rented land. Based on WGS and anonymised cattle data for these 31 isolates, our study aimed to i) determine the amount of genomic divergence among M. bovis isolates within and among herd outbreaks, including subsequent outbreaks affecting the same farm; ii) compare this to genomic diversity seen among badger isolates sampled from the same spatial area and time period as cattle; iii) test whether genetic distances among isolates correlate with either spatial distance among farms or the probability of past movement events connecting them; iv) compare the genetic results to independently obtained contact structures identified by fitting network and transmission models to the herd history data. Results WGS from the 31 bTB isolates revealed a total of 39 polymorphic sites, of which 7 were shared amongst two or more isolates (Table S2 in Text S1). One of these shared polymorphisms was identified as potentially unreliable and thus excluded from the final data set (see Materials and Methods). Because M. bovis is believed to be clonal [15] these genome-wide data can be combined to establish the phylogenetic relationships among isolates (Figures 2 & S1). Using the known sampling dates, we found an increase in genetic divergence from the tree root through time, consistent with a molecular clock (Figure S2). The estimated evolutionary rate was 3.40 (CI: 0.87–5.93)×10−8 substitutions per site per year (equivalent to 0.147 substitutions per genome per year), about an order of magnitude lower than within-host mutation rates recently estimated for M. tuberculosis using WGS [3]. 10.1371/journal.ppat.1003008.g002 Figure 2 Maximum likelihood network of M. bovis genomes with tips arranged according to sampling date. Position of other nodes is simply shown for convenience and does not reflect known branch times. Black circles represent single nucleotide polymorphisms separating sequences, dashed lines indicates branches without mutational events. The size of the circle proportional to the number of isolates sharing the same sequence. The inferred phylogeny revealed that most outbreaks involved genetically distinct isolates (Figure 2) demonstrating the suitability of WGS to track M. bovis spread at the herd level. The more extensively sequenced outbreaks were dominated by a single common genotype (see below). Repeated outbreaks within the same herd tended to involve closely related isolates falling into the same genetic lineage (e.g. Herd 1: 1999/2004; Herd 5: 2008/2010). The exception was Herd 3, which fit this pattern from 1999 to 2004 and from 2007 to 2010 but with distinct lineages causing outbreaks during those two periods. Though few badger isolates were available, their genetic distances from cattle isolates were small and comparable to those observed among cattle isolates; two of the five M. bovis genomes sequenced from badgers were genetically indistinguishable from those seen in cattle during the same year. No M. bovis sample taken from a badger was more than four mutational steps away from the nearest cattle isolate. Interestingly, the two isolates from different tissues of the same badger were separated by five mutational steps, as great a genetic distance as found across different cattle in serial outbreaks within the same cattle herd. This suggests either multiple infections of the same animal or a long-term infection that had allowed for within-host divergence to evolve. In contrast, isolates collected during a series of four outbreaks in Herds 3 and 5 between 2007 and 2010 showed little divergence, with 9 out of the 20 genomes being indistinguishable based on our data, despite the considerable timeframe involved. Outbreaks could be epidemiologically linked through local transmission, which may for example involve badgers, unobserved infections in cattle, environmental contamination or contiguous contact between farms. Alternatively, transmission may be due to the network of livestock movements that connects the five farms through shared links to other farms. We compared the number of mutations separating isolates to two proxy measures of epidemiological distance, 1) the Cartesian distance between the main holding locations of cattle herds, and 2) one minus the relative probability that herds had been connected by movement of infected cattle, called here a “network separation” (see Materials and Methods for our definition of this). As there is a strong genetic and spatial autocorrelation for samples obtained from single herds, we compared only across herds, while noting that this does not completely eliminate the dependence problem. We found a positive relationship between the pair-wise Cartesian distance among farms and the smallest genetic distance among their bTB isolates (Figure 3A), consistent with at least one of the possible local transmission mechanisms being important. As would be expected, where different lineages are considered as separate outbreaks (e.g. the two lineages associated with Herd 3), there was a much poorer correlation (R2 = 0.444) for the earlier outbreaks compared to the later ones. To construct the cattle network, life history data were extracted for all cattle in the database from 01/01/1990 to 31/12/2010 as well as all farms in NI recording a breakdown due to VNTR type 10 between 1999 and 2010, amounting to 58 herds, 3321 cattle and 14258 recorded individual tuberculin test results from 29/07/1993 to 09/12/2010, plus records of post-mortem examination for tuberculous lesions and subsequent confirmatory tests including histology and culture. The movement network included movements amongst all herds in this group of 58 farms. The resulting network was highly connected, including our five farm cluster (Figure S3). The only two farms directly connected through recorded cattle movements were Herds 3 and 5 but this reflected a single animal and, in the absence of further indirect contact, corresponded to a high network separation. In contrast, other farm pairs were much more strongly connected through indirect contact, generally involving one or two intermediary farms. However, these reduced network separations correlated only weakly with genetic distance (Figure 3B). 10.1371/journal.ppat.1003008.g003 Figure 3 Genetic versus spatial and network distances. On both axes, all values are scaled to the maximum value for herd-to-herd interactions. On the x-axes, the minimum number of SNPs differentiating isolates from the two herds (X-Y). In panel (A) above, spatial distance versus genetic distance between herds (black squares). On the y-axis, the cartesian distance between the main holdings of two herds (X-Y) showing a high level of correlation with genetic distance (R2 = 0.720). For reference, the equivalent data for the badger isolates (not fitted) are shown as unfilled circles and diamonds for badger-badger and badger-herd relationships, respectively. Panel (B) below, network separation versus genetic distance between herds. On the y-axis, the network separation defined as (1−pij ), where pij is the probability that herds i and j and linked via cattle movements through the network, considering all possible pathways through any herd from which the same genotype of M. bovis has been isolated, and panimal is the per animal probability of contact that best explains the genetic distance data. The best fit value (panimal = 1.35×10−3) shows a poor correlation with genetic distance (R2 = 0.094). Because of the indeterminate and potentially long timescales for M. bovis transmission, the relationship between dates of sampling and dates of transmission are uncertain. Fitting a mathematical model of transmission to the observed life history data allowed us to compare likely epidemiological processes to the observed tree structure, in order to analyse dynamics at the within-herd, individual animal level. The long period from 2007/8 to 2010 during which no reactors were identified in Herds 3 and 5, combined with the large scale depopulation that occurs following identification of a breakdown (all test reactors are slaughtered and herds are tested every two months until two consecutive clear whole herd tests are achieved), would suggest multiple introductions of infection into these herds. In contrast, the ‘flat’ genetic structure for the same outbreak, with no obvious chains, would suggest one outbreak from a single source. We therefore asked whether the recorded cattle life history and movement data, interpreted in the context of a simple nonlinear model, are consistent with the observed genetic signature. We constructed a compartmental model of transmission, assuming all animals were in one of four states: susceptible to disease (S), exposed (E), potentially test positive but not yet infectious (T) or infectious (I) [18]. Exploiting the explicit animal life histories in all five herds (Figure 4), plus import/export data from the remaining 53 herds in the group, state probability distributions for each animal were generated based on the known dates of bTB detection [6], [19]. Mechanisms other than direct transmission among identified reactors were accounted for with a single external force of infection function ‘ ’, corresponding to the suite of possible local forces as described above, and a fitted force of infection due to latent or hidden infections, where each non-reactor is assumed to contribute to the overall force of infection in proportion to the calculated probability that itself had been infected. The overall force of infection was calculated assuming homogeneous, density dependent mixing within each herd and with contacts between herds via recorded cattle movements. A best-fit model was determined using a likelihood-based Markov Chain Monte Carlo approach (further details in the Materials and Methods section below and in Figures S4 to S6, including posterior distributions for all parameters) [20]. 10.1371/journal.ppat.1003008.g004 Figure 4 Life histories of all cattle with from which Mycobacterium bovis samples of VNTR type 10 were obtained. Showing all individuals residing within the five study herds at some time from 1994 to 2010. Cattle residence times indicated by the length of the horizontal bars (each bar representing a single animal). In black, all cattle from which sequences are derived (herd indicated by surrounding type). Test dates on which one herd received a whole herd test are indicated by vertical dashed lines. Herd colours correspond to colours in Figure 1 (1 – pink, 2 – purple, 3 – blue, 4 – orange, 5 – red). The considerable overlap in cattle life histories (Figure 4) could have potentially allowed for long within-herd epidemics. However, the probabilities of a transmission chain linking cattle-based isolates suggest that while cattle-to-cattle transmission can likely explain some maintenance of individual genetic lineages (at its peak, the cattle to cattle force of infection was calculated to be roughly 50× that from hidden sources), new outbreaks are usually better explained by an unspecified ‘reservoir’ (Figure 5). In particular, despite similar overlaps in reactor life histories across the entire recorded history of Herd 3, and similar observed genetic distances, the earlier outbreak is epidemiologically distinct from the later outbreak in 2007 to 2010, with separate introductions suggested for 1999 and 2004 unless at least one particular individual had an uncharacteristically long infectious period. In contrast, the 2007 to 2010 outbreaks in Herds 3 and 5 are well supported by cattle-to-cattle transmission alone. Herds 3 and 5 are the only ones recorded as having directly traded with each other, however our herd-to-herd analysis suggests that this link was unlikely to have connected the outbreaks in the two herds (Figure 3). All other herds appear only linked by the external force of infection, which was of the same order of magnitude as the force of infection due to a single infectious animal (see parameter estimates in Figure S5). Correlations between genetic distances and network distances are poor (Spearman rank correlation coefficient r = −0.326) with any correlation due to between- rather than within-herd differences, indicating that the WGS data are not of sufficient resolution to make inferences on transmission chains at the within-herd scale. 10.1371/journal.ppat.1003008.g005 Figure 5 Probabilities of pairwise transmission pathways amongst infected cattle with sequenced isolates. The weighted, directed network shows the probability that a transmission path exists between cattle with sequenced isolates that does not pass through other sequenced isolates. Infection events poorly explained by transmission amongst reactor cattle are therefore more likely to be caused by a ‘reservoir’, which potentially encompasses infected badgers, between-herd interactions, latent infections, or environmental contamination. Sequences belonging to the same herd are surrounded by dashed outlines. Discussion Measuring genetic variation at the whole genome scale enabled us to genetically distinguish most isolates of M. bovis. This is particularly notable given the small spatial extent of the study cluster, with no two farms being more than 5 km apart. Compared to traditional typing methods, for which the same genotype may be distributed over many hundred square kilometres and only broad associations can be rigorously defined [21], WGS affords a level of resolution for epidemiological studies previously limited to rapidly evolving RNA viruses [6]–[8]. In addition to most isolates and outbreaks being genetically unique we found that subsequent outbreaks on the same farm tended to involve the same genetic lineage previously detected in that location. This indicates that lineages are commonly able to persist locally within the direct environment of a farm, although the mechanisms for this are not yet clear. Results of our mathematical models based on individual cattle histories indicate that persistence on a farm is overall poorly explained by ongoing transmission within herds. In the cases of Herds 1 and 3 for example, the model identified independent introductions for subsequent outbreaks (Fig. 5), despite the fact that these serial outbreaks involved the same genetic lineage (Fig. 2). Based on these findings, the detected infections (including unsequenced reactors) are insufficient in explaining local persistence on farms, instead suggesting a number of possible alternative mechanisms, such as re-introduction of the same lineage from neighbouring herds, environmental persistence, or alternative hosts. In contrast, and despite the relative simplicity of the modelling approach used, the persistence of the outbreaks in Herds 3 and 5 from 2007 to 2010 are consistent with lineage persistence resulting from extensive within-herd transmission, despite the multiple clear whole herd tests that would have occurred in between dates for reactors. While forward simulations were used to corroborate the robustness of our modelling approach, any extrapolation of our results for bTB epidemiology in general must be viewed with caution, both because of the small size of the dataset and because some of the modelling assumptions (in particular the assumption of explicitly time dependent generation times, see supplementary information) were not intended to be mechanistic descriptions of the underlying transmission process. Nevertheless, the fact that such a parsimonious model identifies a cattle-only contact structure largely consistent with the observed phylogeny generates confidence in our results. In addition to local persistence, we also found evidence for the introduction of new genetic lineages onto farms and our analyses allow us to partially resolve how these introductions occur. Though cattle movements are a known risk for between-herd spread of bTB in Britain [22], [23], they do not appear to be important for the events observed here, as the probabilities of transmission amongst herds involved in the extensive network of all observed VNTR10 outbreaks only poorly predict the genetic divergence amongst isolates. In contrast, we found Cartesian distance to be a good predictor of genetic distance among isolates at a very fine scale. Though the small sample size means that inferences regarding between-herd contacts should be viewed with caution, these results are consistent with the relatively low importance accorded to movements compared to local risk factors in bTB endemic areas that was previously observed at a national scale in GB [23], and suggests that a more extensive analysis of the balance between local spatial and network processes would be merited. As it stands, the most parsimonious explanation for these outbreaks involve a local transmission process that could be due to a number of causes. A non-exhaustive list of these includes both cross-boundary contact or unrecorded local movements between herds and transmission from a common badger reservoir (where the interaction is spatially localised, consistent with the badger's stable social structure and strong territoriality [24]), or a combination of these factors. While our sample size for badgers is low, the badger-derived sequences are remarkably similar to those in cattle, demonstrating very recent cross-over events between the two populations, or alternatively recent infections from a common source, such as the environment [25]. The demonstration of a high M. bovis diversity in a single badger suggests either a lengthy infection in that badger, or multiple exposures to different sources of infection. Although our current estimate for the rate at which M. bovis genomes evolve must be considered preliminary, it is considerably slower than the rate observed in M. tuberculosis [3]. Should it be confirmed, this has obvious implications for the level of temporal resolution that WGS can provide for unravelling epidemiological dynamics for bTB. In our current data, we were unable to genetically resolve relationships among multiple isolates stemming from the same outbreak for example and saw serial outbreaks commonly involving identical genotypes. This is corroborated by the poor correlation between the genetic distances and our estimates of the within-herd contact structure. However, apart from limiting opportunities for molecular epidemiological inference, these observations may also hold clues with respect to M. bovis biology and transmission. A recent study conducting experimental infections with M. tuberculosis in primates, found mutation rates to be equivalent during latent and active infections and proposed oxidative damage as a potential mechanism [3]. If this is relevant to M. bovis, one could hypothesise that the slower rates of evolution seen here at the population level, could be caused by the bacterium spending extended periods outside the host, in the environment. Future studies and analyses are needed to obtain more accurate estimates for the genomic rate of evolution in M. bovis and to test for potential rate heterogeneity and its underlying factors. While cattle movements and long-term, hidden persistence within herds have both been shown to contribute significantly to herd breakdowns [22], [23], [26], these previous analyses were aimed at the identification of statistical associations; here we have shown that WGS data are able to identify local interactions as the principle culprit in specific herds. This makes WGS both a valuable tool for forensic epidemiology, and an aid in the construction of improved mathematical and statistical models of disease dynamics. In contrast, the poor correlation between network and genetic distance at the within-herd level suggests important limits to the resolution that WGS can provide for this system. The local effects identified here may be due to the local infected badger population, but are also consistent with local herd- to-herd spread. Our simplified modelling approach was chosen to maximise the use of available epidemiological contact data, but at the expense of a more detailed exploration of the possible hypotheses regarding the sources of transmission. However, it is likely that WGS based on more extensive sampling will allow for more sophisticated approaches, that could be used to directly estimate the role of badgers in the maintenance of bTB in British and Irish cattle. While insights into particular disease problems will depend on many factors we cannot consider here, our study supports the proposition that WGS data alone can provide insight into the impacts of unobserved populations on observed epidemics even in the absence of detailed transmission chain information, for M. bovis, other members of the M. tuberculosis complex, and other pathogens involving reservoir hosts. Materials and Methods VNTR typing M. bovis is a member of the closely related Mycobacterium tuberculosis complex, which consists of several species with a shared ancestry [27] but which have evolved marked but not absolute host preferences [28], [29]. On a global scale, the M. tuberculosis complex can now be subdivided into discrete lineages, which show strong phylogeographical localisation to regions [30], [31]. VNTR profiling is a genotyping technique based on determining the copy number of a series of short, simple DNA repeats, originally identified by genome analysis [32]. However, while mutations in VNTR type have been observed within the timescale of observation at the regional level, in most cases, probable transmission events are associated with the same VNTR-type, therefore requiring finer resolution typing to order the members of these groups. Herd-level M. bovis genotyping has been performed by the Agri-Food and Biosciences Institute, Belfast, UK since 2003, as follows. The first (disclosing) M. bovis isolate from all bovine TB herd incidents was subjected to genotyping (eight-VNTRs and spoligotyping convention). Heat-inactivated M. bovis cell lysates were used directly as PCR-ready templates. VNTR profiling, spoligotyping, nomenclature, reference strains and quality control were as described [32]. The inferred tandem repeat copy number at each VNTR locus was used to produce a concatenated multi-locus VNTR profile (a string of integers), which was then simplified to a number indicating the prevalence of that profile. Genotype 001 (SB0140), with a spoligotype of SB0140, was the most prevalent in Northern Ireland when surveyed in 1999 to 2003 [32]. Spoligotypes were named according to an agreed international convention (www.mbovis.org). Datasets Anonymised records of cattle tuberculin tests, farm locations and cattle movements among them were drawn from the Animal and Public Health Information System (APHIS), a database containing details of all cattle in Northern Ireland that is administered by the Department of Agriculture and Rural Development [33]. The locations (main farm building) of 58 herds that had a breakdown with VNTR type 10 between 1999 and 2010 were extracted from APHIS. A cluster of five of these herds was selected based on their spatial clustering and the proximity of badger-derived isolates of type 10. The life histories of all 3299 cattle that passed through them since 1995 were compiled, comprising birth and death dates, and the dates of movements into and out of the cluster herds. These animals underwent routine bTB skin testing every year and the lifetime test history of each animal was extracted, containing the dates and results of all tuberculin tests (a total of 14258 individual test results). Results of post-mortem examination for tuberculous lesions and any subsequent confirmatory tests (laboratory based histology, culture and VNTR typing) were also incorporated. Cattle movement was investigated by extracting a broader dataset, comprising all herds that animals passing through the 58 VNTR 10 herds had also visited (a total of 14 096 herds, excluding livestock markets). All cattle movements among the expanded set between 1992 and 2010 were extracted, a total of 5,875,510 individual animal movements. DNA extraction and sequence analysis M. bovis was isolated and confirmed from suspect bovine granulomatous tissue using standard protocols. Confirmed cultures were grown to single colonies on LJ slopes and single colonies were amplified for DNA extraction using the standard CTAB and solvent extraction protocol [34]. Extracted DNA was sequenced at the Sir Henry Wellcome Functional Genomics Facility at the University of Glasgow using an Illumina Genome Analyser IIx. Pair-end reads of 70 bp in length, separated by an average of about 350 bp, were trimmed from both ends based on phred quality scores so as to result in an error rate of 0.001 or less for each base call in the remaining sequence. Reads were mapped to a published UK reference genome (AF2122/97) [35] using the Geneious assembler under the “medium-low sensitivity” option, allowing for a maximum of 10% gaps and mismatches per read [36]. The reference sequence belongs to the same spoligotype (SB0140) as VNTR type 10 and shares identical repeat numbers with it for four out of eight loci used for typing. Mapping resulted in greater than 99% genome coverage with at least 1× and an average read depth of 60–112× for all isolates (see Table S1 in Text S1 for full details). Consensus sequences were generated from the mapped contig based on the quality score sum for each position. A cumulative quality score threshold of 60 (corresponding to an error probability of 1 in 1,000,000) was applied to each position to ensure that accuracy of the final consensus sequence was dependent on both quality and read depth, rather than read depth alone. Below this threshold, the consensus base call was scored as unknown (“N”). Alignment of consensus sequences was carried out using Mauve [37], as implemented within Geneious, assuming collinear genomes and with automatic calculation of seed weight and of the minimum Locally Collinear Blocks (LCB) score. Regions that were difficult to align or which contained >3 consecutive columns of unknown bases or gaps were removed from the final alignment. Similarly, sites that were polymorphic solely due to one or more sequences having ambiguity base calls were removed; this was the only context in which ambiguities were observed. The final alignment, which still represented 99.2% of the reference genome, thus only contained dimorphic single site polymorphisms situated within otherwise invariable regions. After stripping identical sites, a total of 39 SNPs were identified (38 substitutions, 1 deletion, Table S2 in Text S1), of which seven were shared between two or more sequences. All SNPs were examined to confirm their validity before further analysis. Of particular concern was the potential inclusion of spurious SNPs associated with repeat regions in the genome for which mapping may be unreliable. While four of the SNPs were found to fall either in or close to potentially problematic regions, the reliability of the mapping and SNP calling could be confirmed in all four cases (see Supplementary Materials). All SNP calls were supported by at least 38× coverage, with high consistency among reads (usually>95%). The only exception to this was a SNP in position 221927 (G−>A), for which consensus calling was ambiguous in one of the four isolates in which it occurred (Herd5_E_2010, 92×, A: 64%, G: 36%, Table S2 in Text S1). Preliminary analyses further showed that the phylogenetic information provided by this site was in conflict with that of other informative positions (which were in complete agreement). Because these observations raised doubts about the reliability of scoring this SNP as well as about the information it provided, the site was removed from the data set. The final data set was used to generate a maximum likelihood tree using phyml [38] under a Jukes-Cantor model using a heuristic search and the reference genome for outgroup rooting. All sequence data generated for this project are available from the European Nucleotide Archive (http://www.ebi.ac.uk/ena/) under accession number ERP001418. Modelling of transmission chains Using the anonymised life history data of the herds containing genotyped isolates, we adapted a method previously used to study foot and mouth disease virus epidemiology and transmission phylodynamics at the between-farm scale [6], [19]. All cattle were assumed to be in one of four states that reflect the critical points in the cattle infection process [18], [39]–[41]: susceptible (S), exposed (E), where an animal is infected but neither tests positive nor infects others, potentially test positive (T), where the animal can test positive but is not yet infectious and infectious (I), with a constant transition rate from the E to T states and from the T to I states. States T or I can be truncated by a positive test. Knowing the dates each reactor (an animal that tests positive) was detected, we use the transition rates to determine the set of distribution functions for each animal describing the probability of being in each of the four infection states at any time in its life. We used these distributions to determine the likelihood of the observed pattern of reactors in the cluster of herds. This key simplification (conditioning on the observed test data) is less rigorous than an approach where the likelihood is integrated forward across all possible dates of infection and the observed data is treated as a random variable. However the latter approach would have been computationally prohibitive for our system. The simplification implies that mean generation times are consistent over the observed timeframe, even though they would naturally be expected to contract over the course of an epidemic [42]. The assumption is expected to be reasonable if the disease is in an endemic state, or an endemic state is rapidly reached after the initial introduction, where ‘rapid’ is relative to the observed timescale – i.e. the distribution across states is ‘quasi-stationary’, with E and T proportional to the reciprocal of their respective transition rates, i.e. and respectively. Though observation of our incidence data suggested this to be likely, to test this assumption, we allowed the proportion in state I at the time of detection to vary over the course of the epidemic, with probability that it is infectious (i.e. in the I state) and for being in the T state, where t is the time since the genotype was initially detected. Here, we have assumed that has the form , where a and b are fitted constants and is a probability with range [0…1]. We assumed that all cattle are equally infectious if in state ‘I’ but that an infected animal does not contribute to further infection after it is detected; most reactor cattle are removed from the herd within days of being identified. All infection is caused by either infectious cattle or by a reservoir (either undetected infected cattle within the herd or external, probably local, factors). At this population scale, we assumed that the probability of a false positive is negligible. To account for untested and undetected contributions to the epidemic, we assumed all test-negative/untested cattle to have a fitted probability of being infectious. The date(s) of an animal “i” having a positive test is denoted as , at which it may be infectious with a probability , where is the time the outbreak was initially detected in the herd and t is equal to at the point of evaluation. We assumed that a reactor in the infectious state became infectious (i.e. moved from the test sensitive state) within a maximum of days prior to the positive test. is the maximum infectious lifetime of the animal prior to detection. Because a uniform test is applied to all cattle over 6 months of age on an annual basis, the date of first test after infection would fall in the range of 0 to 1 year; therefore we assume that the duration of the effective infectious period is uniform. This has the (conservative) effect of minimising the potential duration of cattle-only transmission chains since it does not allow for long infectious periods. The probability that a reactor was infectious at any time prior to the date it tested positive is: All other model states were assumed to be exponentially distributed. The probability that an individual is test sensitive at a given time t is therefore where is the transition rate from the T to the I states in the forward SETI model. The parameter determines the duration of the state T and therefore the time when the backwards T to E transition occurs, conditional on the backwards transition from I to T and the time of detection . The first term accounts for the animal being in the potentially test positive state at and the second term accounts for the animal being infectious at . Here, is the increase in that occurs due to transition from at time t′. Similarly the probabilities of being in the exposed and susceptible states are where is the transition rate from the E to the T states in the forward-in-time model. The forward-in-time transition from the S to the E state is determined by and is estimated via the calculated force of infection based on the probabilities of being in the I state. It is incorporated directly into the calculation of the likelihood as described below. A similar approach was used for animals identified by a post mortem test in which case is the date the animal was post mortem examined. We assumed further that post mortem testing only identifies the same infection classes as the standard tuberculin test (i.e. only T or I class individuals are detected). We calculated the probability that an animal became infected by considering the forces of infection on it during its lifetime. The life history of each animal that resided in the cluster between January 1990 and 31 December 2010 was converted into a lookup table of the herds in which it resided in fortnightly time steps. If an animal moved between herds in a time step it was treated as belonging to both herds since it can contribute to an outbreak in both herds during that step. Animals born into a herd are added to the lookup table for that time step and conversely death results in an animal being removed from the lookup table – movements back to the herd cluster via other herds are also retained. Using this lookup table we were able to calculate the probability an animal became infected in each period. We calculate the probability that an animal would have become infected at any time t, , given that is was a known reactor based on the force of infection calculated from the model. If an animal was a reactor, then the incremental probability that an animal, i, was infected by all sources over the time interval t, t+dt is therefore given as (1) calculated using the explicit infection histories. If the animal never tested positive then (2) i.e. we consider the probability that an animal was infected despite never testing positive. Young calves destined for slaughter at a young age are test exempt and are considered to be a negligible risk due to their short lifespan; they are therefore excluded from the analysis. The parameter is the transmission rate from a contact with an infected animal in the same herd at the same time (the summation is over animals in the same herd at t). The reservoir term, , can be written as i.e. a combination of any external forces of infection, , (e.g. infected wildlife) that is represented by a fixed rate and , the fraction of herd at t that never tested positive, for which the contribution of latent infections is summarized by the probability Platent . We calculated the probability each animal was infected during its lifetime by integrating over the animals' entire life (approximated by a summation over fortnightly time steps) and thus calculate the likelihood for the model as: Here, is the fitted sensitivity of the test for bTB infection. Consistent with the low number of reactors in officially TB free countries such as Scotland, the test specificity is assumed to be effectively 100%. Using L, we used the Metropolis-Hastings algorithm to derive the posterior distributions over the parameter space defined by . The length of burn-in for the MCMC simulations was checked using multiple MCMC chains with dispersed starting points, perturbing a subset of the parameters at each step in the chain (further descriptions in the supplementary information, and Figures S4 and S7 for the evolution of the MCMC chains and the posterior distribution of the parameters). Prior distributions were chosen to be uniform, with ranges , , , , , , , and , based on the parameters estimates previously reported in various field and experimental studies [43]–[45]. A running sum of the probability of a transmission event ( ) was retained, ignoring potential transmissions that were incompatible with the phylogenetic tree (for example between animals with isolates in different clades of the tree, excluding these animals from being linked by transmission). We used a modified Dijkstra's algorithm [46] to identify all possible infection paths between cattle for which we have isolates. Each link “i” in a defined chain has a probability pi of being associated with transmission, so that the probability of that transmission chain occurring is given by the product of the pi's. Then the probability that a path will exist between two individuals A and B is simply the probability that at least one of the identified possible chains will connect the two, and therefore The calculation included all chains that do not pass through another sequenced isolate, and for computational convenience, we have limited the calculation to chains of less than three individual intermediate cattle. In order to test our assumptions, we ran forward simulations of a SETI model, based on the most likely parameter values of the posterior distributions from the outbreak data analysis. Applying our estimation approach to the simulated outbreak data confirmed that all the input parameter values fell within the 95% credible intervals of their corresponding posterior distributions derived from the original analysis of the data (not shown). Identifying the best-fit network distance between herds Similar to the approach described above, all paths amongst Herds 1 to 5 that do not pass through the remaining herds were calculated. The link strength probability between Herds C and D was taken to be where the product considers any direct path (only between Herds 3 and 5), all chains of path length 2 (i.e. with one intermediate), or where these do not exist, chains of path length 3. The parameter nlink is the number of cattle moved between pairs of premises in the chains. The ‘network separation’, defined as {1−PCD }, was used as a measure of within-herd network distance, and this was plotted against the minimum genetic distance between the two herds. The value of panimal (i.e. the relative weighting assigned to each animal moved, where 0

- Record: found
- Abstract: found
- Article: found

Introduction The number of cattle herds in Great Britain (GB) placed under movement restrictions due to the suspected presence of bovine tuberculosis (bTB) has progressively increased over the past 25 years [1]. This increase in the rate of so-called “breakdown” herds is despite an intensive and costly test-and-slaughter control program [2]. Recent studies have focused on estimating the contributions of cattle movements and wildlife transmission to incidence [1],[3]–[8] as measured by the rate of new breakdowns. However, less attention has been paid to quantifying the dynamics of transmission within herds, even though this is arguably the most data-rich unit within the wider ecology of M. bovis. Previous history of disease within a herd is an important predictor of breakdown [5], [8], [9], with 38% of herds that clear movement restrictions experiencing a recurrent incident within 24 months [10]. This high rate of recurrence suggests that infection may be persisting within herds in the face of repeated testing. In GB and internationally, detection and clearance of herds is dependent on variants of the imperfect tuberculin skin test. In GB and Ireland this takes the form of a single intra-dermal comparative cervical tuberculin (SICCT) [11] test. Infection missed by SICCT testing is likely to be contributing to recurrence within herds. However, the lack of a gold-standard diagnostic test for bovine tuberculosis means that the efficiency of the SICCT test is poorly characterized and the contribution of missing infection to recurrence cannot be easily assessed. Furthermore, reactivity to the SICCT test is dependent on the time-from-infection [12], most often characterized [13], [14] as an occult period where animals are infected but not yet detectable. As a consequence, the efficiency of testing within a herd depends not only on the characteristics of the diagnostic test, but also on the competing timescales of transmission and the frequency of testing. Within-herd models of bTB have been developed to address this problem with a view to informing government policy [13], [14]. However, extant models have been based on ad-hoc parameterizations informed by disparate experimental studies and expert opinion. There exists considerable uncertainty in the assumed values of key parameters, in particular the occult period, the scaling of transmission rates with herd size [15], [16] and the duration of latency between infection and infectiousness. To address this uncertainty, we propose a novel basis for parameterization of within-herd models using measures of stochastic persistence [17] as a metric for Approximate Bayesian Computation (ABC) [18]–[19]. Persistence measures have proven to be a powerful probe on which to parameterize models of childhood infectious diseases [20]–[23]. A successful approach has been to assume that the intrinsic rate of transmission within a population is rapid compared to the combined rate of transmission from sources extrinsic to the local population. This time-scale separation allows us to model local populations independently. Extrinsic routes of transmission are modeled through a generalized infectious pressure [20]–[22]. Comparatively less theoretical attention has been paid to modeling the persistence of managed endemic diseases. For chronic diseases, such as bTB, demographic turnover of the population is the only natural mechanism acting to clear infection from populations. In this context persistence can be used as an indirect measure of the efficiency of diagnostic testing. In this study we model the within-herd persistence of bTB as a balance between three key processes: the infectious pressure acting to introduce infection into the herd from extrinsic sources, the rate of cattle-to-cattle transmission within the herd and the rate of removal of infection through testing and demographic turnover. Herds are considered as isolated populations loosely connected to a reservoir of infection modeled as an infectious pressure. We are therefore not concerned with modeling the routes of introduction to the herd – which may be through movements of infected animals or contact with wildlife reservoir populations. Instead we focus on the processes of transmission within a herd with relation to the detection and resolution of breakdowns. We do so using two mechanistic models of within-herd transmission that we parameterize using routinely collected epidemiological data. We finally apply our parameterized models to estimate the hidden burden of infection and its implications for control of bTB in Great Britain. Results Measures of persistence of bTB The probability of extinction within epidemic models is dependent on the past history of infection within the population [24]. Alternative empirical measures of persistence, that capture different aspects of the transmission dynamics, can be constructed depending on how we condition on the past history of infection [21]. For bTB, infection missed during a given test is likely to contribute to the probability of the herd failing subsequent tests. Contingent on the natural time-scale of transmission and the scheduling of testing, missing infection may act to prolong the duration of breakdowns and/or increase the probability of a recurrent breakdown. We therefore quantify within-herd persistence through two competing measures related to the duration and the rate of recurrence of breakdowns. The duration of breakdowns is captured by the probability that breakdowns are prolonged [25], defined as lasting longer than 240 days. Recurrence is captured through the probability of a breakdown recurring [10] within a fixed time horizon of 6, 12 and 24 months after the end of a breakdown. Study population Bovine tuberculosis is a statutory infectious disease. Incidence and testing data are routinely collected by the Animal Health and Veterinary Laboratories Agency (AHVLA) and collated within the VetNet database. It is not feasible or desirable to model the full complexity of the British testing regime within a herd level model. In addition to the schedule of statutory surveillance testing, there exists a diverse range of auxiliary tests including pre-movement testing, contiguous tests and tests trigged by epidemiological investigations of “at risk” premises. As a consequence there is considerable variability in the duration of time that infection can spread unobserved within herds before detection. In an attempt to control for this uncertainty we restrict our current analysis to new breakdowns with start dates between 2003–2005 that were detected through routine surveillance, either through the slaughterhouse or by the detection of reactors at a routine or whole herd test (tests classified as ‘VE-WHT’, ‘VE-WHT2’,‘VE-RHT’ or ‘VE-SLH’). We choose to restrict our study to the period 2003–2005 due to systematic changes to the testing system surrounding the 2001 foot-and-mouth disease epidemic [5] and a later increase in use of the gamma-interferon test [26]. Although discretionary use of the gamma interferon test increased after the end of our study period, this does not appear to have impacted upon persistence and our model still provides an equally good fit to target measures taken from more recent data (Figures S7 , S11 ). The cessation of testing during the 2001 foot-and-mouth epidemic artificially increased the duration of time that herds were kept under movement restrictions, delayed the scheduling of routine surveillance tests and was associated with an increase in incidence and spread of bTB to new areas [5]. Given the slow rate of transmission of bTB, the perturbative effect of this disruption to testing is likely to have continued for many years. However, our interest is primarily in the sequence of transmission and testing after the disclosure of a breakdown. Disruptions to testing prior to the disclosure of infection will increase the duration of time that infection was able to transmit silently within our study herds. However, as this variation is included within the empirical testing distributions used within our model (Figure S1 , Dataset S1 ) we do not expect it to affect our results. Of 10,174 breakdowns recorded within our study period, 3,456 (34%) breakdowns match our criteria for inclusion. Restricting our analyses to this sub-population has an important advantage. The scheduling of surveillance tests in GB is based on the local incidence (Figure 1 ), that determines the so-called parish testing interval (PTI). The duration of time that infection may have remained undisclosed within our sub-population of herds is strongly constrained by the herd's local PTI. After detection of infection, regardless of the disclosing test type, all herds must undergo the same statutory sequence of testing. We therefore do not believe that our inclusion criteria impacts upon the generality of our inference for rates of within-herd transmission and the efficiency of surveillance. 10.1371/journal.pcbi.1002730.g001 Figure 1 Parish Testing Intervals (PTI) for bTB in GB (2003–2005). Routine surveillance for bTB in GB is based upon the regular (SICCT) testing of herds at a frequency determined by the local incidence of affected premises. Panel (A) maps the shortest recorded PTI for each parish over our study period of 2003–2005. High incidence areas are spatially clustered with the greatest incidence, and thus intensity of testing, in the south-west of England and south Wales. PTI therefore offers a crude categorization of herds according to epidemiological risk, past history of testing and to a lesser extent, geographical location. In contrast to surveillance testing, the sequence of tests (B) following a breakdown are dependent only on the outcome of tests on the affected premises. A failed surveillance test leads to a sequence of short interval tests (SI) at intervals of at least 60 days. Confirmation of infection, through isolation of M. bovis or evidence of visible lesions at slaughter, leads to tests being re-interpreted under a “severe” interpretation of the SICCT, until one test is passed; testing then continues until an additional test at the standard interpretation is passed. After a breakdown is cleared follow-up tests are scheduled at intervals of at least 6 and 12 months, after which testing frequency reverts to the local parish testing frequency. Of the remaining breakdowns the majority are either recurrent breakdowns (2,102; 21%) initiated by a follow-up ‘VE-6M’ or ‘VE-12M’ test or breakdowns that started with a so-called “inconclusive” reactor (2,032; 20%). Inconclusive reactors (IRs) demonstrate a response to the SICCT that is close to the cut-off value defining a reactor. IRs do not necessarily trigger a breakdown but require the animal, rather than the herd, to be retested at an interval of 60 days. The population of IRs will be composed of both false reactor and truly infected animals and cannot be rationally treated within our model framework, requiring us to omit these herds from our analysis. The remaining 25% of breakdowns were initiated through a mixture of contiguous testing of affected premises and contact tracing. The persistence of bTB has previously been demonstrated to scale with herd size [27]; a known risk factor for both prolongation [25] of breakdowns in GB and recurrence in Irish herds [28]. We extend these analyses to quantify the relationship of our two persistence measures with herd size. The size of a cattle herd varies dynamically, even over the course of a breakdown. We therefore define herd size as the maximum herd size over the breakdown. The distribution of herd sizes is right skewed with 90% of breakdowns having herd sizes less than or equal to 360 cattle (Figure S4). The scarcity of herd sizes larger than this limits our ability to measure the relationship with persistence [27]. We therefore finally restrict our study population to breakdowns with a herd size of less than or equal to 360 cattle, leaving us with a final study population of 3,094 breakdowns. These 3,094 breakdowns were then binned into 6 groups with histogram mid-points of [30,90,150,210,270,330]. Patterns of persistence of bTB in Great Britain We further stratify these herds by the parish testing interval (PTI) and confirmation status of breakdowns to produce empirical distributions of persistence (Figure 2 , Figures S6 , S10 ). This classification is motivated by the systematic differences in testing for herds within different PTIs and after confirmation. Confirmed (recently re-classified as Officially TB-Free Status Withdrawn or OTF-Withdrawn) breakdowns are required to pass an additional clear test at a more strict (severe) interpretation of the SICCT test (Figure 1B ). The severe interpretation increases the sensitivity of the SICCT test at the expense of reducing the specificity. Confirmation is triggered by the discovery of reactor animals with visible lesions and/or culture of M. bovis. 10.1371/journal.pcbi.1002730.g002 Figure 2 The persistence of bTB in GB herds (2003–2005). The within-herd persistence of bTB in GB as measured by the probability of all GB breakdowns from our study population being prolonged (duration of greater than 240 days, top panel) or recurrent within 24 months (middle panel). The relationship of each measure is plotted against herd size, with breakdowns further stratified by parish testing interval (PTI 1, 2 and 4, left to right) and confirmation status (unconfirmed breakdowns: lime green, circles, confirmed breakdowns: magenta squares). Uncertainty in each (mean) target observation (thick lines) is illustrated by an envelope (thin lines) of ±1.96 standard errors around the mean. Predictive distributions from our within-herd (SORI) model for each of these measures are plotted as shaded density strips where the intensity of color is proportional to the probability density at that point [35]. The proportion of prolonged and recurrent breakdowns both scale with herd size, but demonstrate distinct relationships with respect to both confirmation status and the local background risk of infection as measured by PTI. These empirical relationships are consistent with previous analyses suggesting that confirmation is associated with an increase in the duration of breakdowns [25], but has negligible impact on recurrence [10]. In contrast to the consistency of the duration of breakdowns across all areas there is a marked increase in the rate of recurrence with local risk as measured by PTI. These differential relationships of persistence with herd size are the basis on which we set out to infer the within-herd transmission dynamics of bTB. Modelling the persistence of bTB We consider the persistence of bTB to be a product of the non-linear interaction of both the disease and testing dynamics. Heuristically, our model can therefore be considered as having two interacting dynamic components: an epidemic model that describes transmission within and into the herd and a testing model that models the sequence of tests and removal of reactors. We estimate the parameters of our model (Table 1 ) directly from our target measures of persistence using a sequential Monte Carlo implementation of Approximate Bayesian Computation (ABC-SMC) [18], [19]. This framework is designed to bypass the calculation of computationally infeasible likelihood functions and instead generates distributions of parameters for which model outputs are consistent with the data according to a set of pre-defined goodness-of-fit metrics. 10.1371/journal.pcbi.1002730.t001 Table 1 Model parameters. Parameter Description Standard Sensitivity. Probability of positive test from R,I compartments at standard definition. Standard Specificity. Probability of negative test from all compartments at standard definition.(1 – probability of a false positive ) Severe Sensitivity. Probability of positive test from R,I compartments at severe definition Severe Specificity. Probability of negative test from all compartments at severe definition.(1 – probability of false positive ) Confirmation. Probability of confirmation of breakdown per reactor based on slaughter-house inspection and culture ( only applies to reactors from O,R,I compartments) Slaughterhouse routine inspection. Probability of detecting an infected animal (O,R,I compartment) at slaughter under routine inspection Occult Period. Mean length of time that animals are undetectable (occult) to SICCT Reactive Period. Mean length of time between infection and animals becoming infectious Transmission parameter associated with density dependence (rate per day, dimensions change with q) Transmission parameter measuring the strength of density dependence (range 0–1) Transmission parameter measuring infectious pressure per susceptible per year in PTI 1 Transmission parameter measuring infectious pressure per susceptible per year in PTI 2 Transmission parameter measuring infectious pressure per susceptible per year in PTI 4 Constant equal to mid-point of range of herd sizes within study population. Used to transform density dependence of force of infection. The period of latency between infection and infectiousness is a key epidemiological parameter that sets the time-scale between subsequent epidemic generations. Given the chronic, progressive nature of bTB, models have conventionally assumed long epidemiological latent periods of ∼6–20 months [13], [14]. However, animal challenge studies have suggested that bacterial shedding, and therefore transmission, may occur over shorter time-scales of ∼30 days [29]. In order to explore these two scenarios of latency, we fitted two models using vague (uniform) prior assumptions (Table 2 ), one based on the conventional model structure assumed for bTB [13], [14] (SORI, Table 3 ) and an alternative model allowing for “early” infectiousness (SOR, Table 4). In the SORI model cattle are assumed to be infectious (I) only after passing through two latent stages: an occult stage (O) where animals are infected, but not detectable by SICCT testing and a reactive stage (R) where animals are ‘reactive’ to the SICCT test but not yet infectious. The SOR model decouples this relationship between epidemiological latency and reactivity to the SICCT test. Animals are assumed to be potentially infectious from both the occult (O) and reactive (R) classes, dispensing with the infectious class (I). 10.1371/journal.pcbi.1002730.t002 Table 2 Prior assumptions. Parameter Prior Constraints Initial sampling distribution Uniform [0.05,1] Uniform [0.05,1-0.9997] Uniform [0.4,1] Uniform [0,1-0.9990] Uniform [0.0,0.5] Uniform [0.0,0.5] Uniform [0.0,1.5] Uniform [0.0,0.35] Uniform [0.0,1.5] Uniform [0,2.0] Uniform [0,1] Uniform [0,3e-4] Uniform [0, 3e-4] Uniform [0, 3e-5] 10.1371/journal.pcbi.1002730.t003 Table 3 Events defining SORI stochastic epidemic model. Event Effect Probability per unit time Move susceptible animal onto herd Remove animal from herd Transmission Emergence (Occult) Emergence (Reactive) 10.1371/journal.pcbi.1002730.t004 Table 4 Events defining SOR stochastic epidemic model. Event Effect Probability per unit time Move susceptible animal into herd Remove animal from herd Transmission Emergence Both models provide comparable fits to the empirical target distributions (Figure 2 , Figures S6 , S10 ) despite remarkably different estimates for the duration of the occult period and the epidemiological period of latency between infection and infectiousness (Figures S5 , S9 ). The occult period is estimated to last ∼1.8 days (0.0–7.7 days, 95% CI) for the SOR model and ∼275 days (24–517, 95% CI) for the SORI model. Whereas infectiousness is implicitly assumed to begin with infection under the SOR model, the period of epidemiological latency estimated under the SORI model is more than twice that of previously assumed values [13], [14] with a point estimate of 406 days (116–827, 95% CI). The occult period estimated from the SORI model is an order of magnitude larger than the range of 8–65 days observed in animal challenge studies [12], [30], [31]. Placing a more informative prior (uniform on the range 0–128 days) on the occult period for the SORI model has no appreciable impact on the fidelity of the model fit and reduces the estimated occult period to a median point estimate of ∼28 days (1–119 days, 95% CI). We therefore select the SORI model fit with the more informative prior for comparison with the alternative SOR model within this paper. Both models estimate that the rate of cattle-to-cattle transmission within a herd increases, non-linearly, with herd size. The potential for transmission within a herd can be characterized by the basic reproductive ratio R0, defined as the expected number of secondary cases within a herd of size N on the introduction of a single infectious individual. Within the range of our study population our (median) point estimate of R0 from the SORI model increases from 1.5 (0.26–4.9; 95% CI) in a herd of size 30 up to 4.9 (0.99–14.0; 95% CI) in a herd of 400 cattle (Figure 3A ). Estimates from the SOR model are smaller, but with overlapping credible intervals, increasing from 0.52 (0.1–1.6, 95% CI) in a herd of size 30 up to 3.6 (0.73–8.85, 95% CI) in a herd of size 400 (Figure 3A ). As a consequence, both models predict that the efficiency of control will also scale with herd size. 10.1371/journal.pcbi.1002730.g003 Figure 3 Herd level measures of efficiency of transmission and clearance of infection for SORI and SOR models. (A) Predictive distribution for the within-herd reproduction ratio (R0), plotted as a shaded density strip where the intensity of shading is proportional to the probability that R0 takes a given value [35]. We calculate where measures the strength of density dependence, is the estimated transmission parameter and the per capita rate of turnover of the herd sampled from an empirical distribution and . Both models estimate that transmission is non-linearly density dependent with point estimates for the invasion threshold (R0 = 1) of 71 (8–674,95% CI) for the SORI model and 12 (1–833, 95% CI) for the SOR model. (B) The effective sensitivity of the SICCT test within our study population of herds measured under the standard (solid line) and severe (dashed line) interpretations and relative to the gold standard of confirmation with visible lesions (dotted line). (C,D) The infectious burden remaining after movement restrictions are lifted can be characterized in terms of the number of herds with at least one infectious animal remaining when movement restrictions are lifted (C) and the expected distribution of animals on those herds (D). Predictive distributions for (B,C,D) are calculated by simulating from the empirical distribution of herds taken from our study population, representing the national distribution of herds. The SOR and SORI models provide contrasting estimates of the efficiency of SICCT testing in Great Britain. The SOR model estimates a true median SICCT test sensitivity of 66% (52–80%, 95% CI) at the standard interpretation, rising to 72% (56–88%, 95% CI) under the severe interpretation. Estimates of true sensitivity from the more traditional SORI model are far lower, at 36% (24–51%, 95% CI) for the standard interpretation rising to 48% (34–69%, 95% CI) for the severe interpretation (Figure 3B ). However, these model estimates are relative to the true infection status of animals. Given the lack of a gold standard diagnostic test for bTB, such information is not available in real populations. Rather, out of necessity, sensitivity of diagnostic tests for bTB is routinely measured relative to the presence of visible lesions and/or culture. The limitations of such comparative measures of sensitivity are aptly illustrated by comparison of our estimates from the SOR and SORI models. Despite the differences in the true sensitivity between the two models, the effective test sensitivities relative to visible lesions are indistinguishable and considerably higher than the true values (Figure 3B ). As a consequence the diverse sensitivity estimates from the SOR and SORI models are nonetheless both consistent with published estimates of the sensitivity of SICCT relative to visible lesions of up to 93.5% at the severe interpretation [11]. In order to quantify the efficiency of control we introduce a new measure - the infectious burden. We define infectious burden as the probability that at least one infected animal remains within a herd after movement restrictions are lifted. By simulating our within-herd models using the distribution of herd sizes from our study population we can generate predictive distributions for the infectious burden at the national level (Figure 3 C,D ). Once more, the SOR and SORI models provide contrasting views of the efficiency of control. Under the SOR model we estimate that 8% (3–17%; 95% CI) of breakdowns will have an infectious burden when they clear restrictions, with a median of 1 (1–3; 95% CI) infectious animal remaining in these herds. Under the SORI model this estimate increases to 21% (12–33%; 95% CI) of breakdowns with an infectious burden when they clear restrictions, with a median of 1 (1–4; 95% CI) infectious animal remaining in these herds. Implications for control We apply our fitted models to predict how different herd-level interventions may affect the resolution of breakdowns (Figure 4 ). Specifically we consider two treatments: application of a ‘perfect test’ and eliminating the extrinsic infectious pressure through ‘perfect isolation’. In the SOR model, recurrence is driven almost completely by re-introduction, with ‘perfect isolation’ having the potential to eliminate recurrence completely in some herds (Figure 4B ). Perfect isolation is predicted to be less effective at reducing recurrence in herds with a low extrinsic rate of re-introduction (i.e. small herds in low incidence areas). Although these herds are predicted to have a lower infectious burden (Figures S8 , S12 ), when they do experience a recurrent breakdown it is more likely to be caused by infection missed by SICCT testing than re-introduction. 10.1371/journal.pcbi.1002730.g004 Figure 4 Impact of herd-level interventions on probability of recurrence within 24 months. Change in the probability of a herd experiencing a recurrent breakdown after application of a ‘perfect’ test (left column) or perfect isolation (right column). The perfect test is assumed to have 100% sensitivity and specificity and no occult period. Perfect isolation corresponds to setting the extrinsic infectious pressure to zero at the end of a breakdown ( ). Plotted values correspond to the average % difference in the probability of recurrence relative to the fitted SORI (Panel A) and SOR models (Panel B). Separate series are plotted for herds in PTI 1 (lime green circles), 2 (magenta squares) & 4 (sky blue diamonds). Predictive distributions illustrating the variability in these point (mean) estimates are presented in supplementary figure S13 . Under the SORI model there is a similar relationship in the response to ‘perfect isolation’, except that a greater proportion of recurrence is attributable to persistence of infection. At the national level, averaging over our study population of herds once again, we estimate that 50% (33–67, 95% CI) of recurrent breakdowns are attributable to persistence within the SORI model, compared to 24% (11–42, 95% CI) under the SOR model. However, eliminating this hidden burden of infection is not sufficient to eliminate recurrence if the extrinsic infectious pressure acting on herds is not simultaneously addressed. Under both models a ‘perfect test’ with 100% sensitivity, specificity and no occult period fails to improve the probability of recurrence in high incidence areas (Figure 4 ). Although the perfect test reduces the duration of breakdowns, it can also detect infection within herds more quickly. In the short term, rates of recurrence will therefore increase in high incidence areas and such a perfect test would only be of benefit for low incidence (PTI 4) breakdowns where the rate of re-introduction is sufficiently low. This counter-intuitive result demonstrates an important limitation of our approach. Our herd-level model does not distinguish between movements to slaughter or to other herds, so the infectious burden output from our model may potentially be contributing to the extrinsic rate of transmission that drives recurrence in our herd-level model. Both of our within-herd models can equally well fit the empirical patterns of persistence of bTB despite very different predictions for the level of the infectious burden. However, such a difference would place very different weights on the importance of cattle movements in network models of herd-to-herd transmission. Recent analyses of the between-herd transmission of the disease in GB have simplified, or ignored these within-herd dynamics of transmission [6], [7]. Discussion A fundamental challenge in epidemiological modeling concerns identifying the appropriate level of model complexity required to understand the dynamics of transmission and form a rational basis for policy development. Tuberculosis has been described as an infectious disease with a period of latency ranging from one day to a lifetime [32]. However, this uncertainty surrounding the progression of disease in individuals is rarely considered as a part of epidemiological modeling studies. In this study we have demonstrated how assumptions concerning the relative timing of infectiousness and reactivity to tuberculin profoundly impacts upon the estimated efficiency of SICCT testing. Both the SOR and SORI models are equally well supported by the population level data used in this study, despite very different estimates for the efficiency of testing. This suggests that persistence measures alone are insufficient to distinguish the true burden of infection and points to experimental studies that could resolve this uncertainty. Neither model identifies, without the support of informative priors, an occult period within the range observed from animal challenge studies [12], [30], [31]. However, if there is a relationship between infectious dose and the duration of latency, estimates from challenge studies must also be treated with caution. Given the importance of this parameter in determining the hidden burden of infection, further research is required to clarify the relationship between infectiousness and sensitivity to diagnostic tests. Our modeling suggests that transmission of bTB ‘early’ in infection necessitates a lower level of persistence of infection than predicted by traditional (SORI) transmission models. However, evidence for such ‘early’ transmission comes from animal challenge studies [29], [33] and has not been verified under natural transmission conditions. Our modeling emphasizes the critical importance of understanding how the pattern of bacterial shedding in naturally infected animals changes over time. Both models estimate that the rate of cattle-to-cattle transmission in GB herds is non-linearly density dependent. This result has immediate importance for the formulation of bTB policy at the herd level, suggesting that additional controls may need to be targeted towards larger herds. Our models suggest that the key to addressing the ongoing spread of bTB lies with reducing the rate of transmission into herds. The central question remains as to whether this requires management of the reservoir of infection in wildlife populations, or simply improved surveillance and diagnostic testing to reduce the movement of undisclosed infection between herds. We have shown that stochastic persistence measures can provide insights into the efficiency of control measures for managed populations. However, the interpretation of these patterns of persistence requires a modeling approach that simultaneously accounts for the dynamics of control as well as the intrinsic dynamics of disease transmission. In the case of bTB, the dynamics of infection at the individual level have a profound impact on the estimated burden of infection missed by testing. It is therefore imperative to improve our understanding of the, still mysterious, life history of infection of bTB in individual cattle. Methods Epidemic model Within-herd transmission of bTB is modeled using the standard compartmental approach where animals are classified only by their epidemiological status. We consider two alternative models (SORI and SOR) corresponding to different assumptions concerning the relationship between latency, transmission and reactivity to the SICCT skin test. In the traditional SORI model for bTB, occult (O) and reactive (R) animals are infected but not yet infectious, differing only in their response to the SICCT test. Infectious animals (I) are assumed to be both infectious and detectable by the SICCT test with the same efficiency as reactive animals. In the SOR model occult (O) and reactive (R) animals are both assumed to be potentially infectious eliminating the need for the infectious class (I). Both epidemic models are implemented as stochastic Markov chains in continuous time and can be defined by the allowed transitions between the four state variables: Susceptible (S), Occult (O), Reactive (R) and Infectious (I) (Tables 3 , 4 ). A per capita turnover rate μ is sampled from an empirical distribution for each simulation (Figure S3 ). A constant target herd size ( ) is maintained by balancing a constant per-capita removal rate ( ) with a fixed import rate of . Herd size therefore fluctuates, with an instantaneous herd size of . The extrinsic infectious pressure ( ) is the only parameter to vary with the parish testing interval (PTI) taking unique values for PTI 1, 2 and 4 ( , , ). Testing model The sequence of tests before, during and after a breakdown is simulated by a model where the timing of tests and number of animals to be tested changes dynamically according to the state-variables of the epidemic model and the outcome of individual animal tests. Model simulations are initialised with the entire herd in the susceptible compartment (S,O,R,I) = (N, 0, 0, 0). The model is then simulated forward, piecewise, between the dynamically scheduled tests before, during and for 5 years following the end of the first breakdown, or until a recurrent breakdown is triggered. The sequence of decisions following the outcome of herd tests is summarized in Figure 1B . Simulations begin with the herd undergoing routine surveillance through slaughterhouse inspection and whole herd tests (classified as RHT or WHT) at 1, 2 or 4 yearly intervals (described below). Detection of a reactor animal triggers a breakdown. The herd then enters a sequence of short interval tests (SIT). Unconfirmed breakdowns end after a single clear test at the standard interpretation, while confirmed breakdowns must clear two tests – one at severe interpretation and the second at standard interpretation. Two follow-up tests, one six months after the end of a breakdown (VE-6M) and one 12 months later (VE-12M) are then scheduled. The time between all tests associated with a breakdown (SIT, VE-6M, VE-12M) are sampled from empirical distributions (Figure S1 , Dataset S1 ). The duration of time between routine tests is also sampled from an empirical distribution (with separate distributions for PTI 1, 2 and 4) to account for the additional variation in the time to detection that is a consequence of delays in testing and the transition of herds between different parish testing intervals (Figure S1 , Dataset S1 ). Breakdowns are triggered by the detection of a reactor, either due to the presence of infected animals in the herd or the generation of a false positive test result. Nominally, we simulate the full sequence of tests until either of these events occurs with the proportion of false-positive breakdowns determined by the relative values of the specificity ( ) and the infectious pressure ( ). In practice, and to increase the speed of simulations, this can be pre-calculated by explicitly calculating the probability of a false positive breakdown occurring between periods where there are no infectious animals within the herd. Breakdowns can also be triggered by routine slaughterhouse surveillance that is modeled as a fixed probability ( ) that removals from the O, R and I compartments will be detected. Breakdowns triggered by slaughterhouse surveillance are treated as confirmed breakdowns, with the first whole herd test carried out under the severe interpretation. Simulating herd tests The application of herd tests in GB can be modeled by simulating three basic processes 1) the number of animals to test 2) the number of reactor animals detected at the standard and severe interpretations of the skin test and 3) the confirmation process. For tests associated directly with a breakdown (SIT, VE-6M, VE-12M) the whole herd is tested. However, there is more variation in the type of test, and numbers of animals tested, in PTI 2 and 4 herds. We simulate this process by choosing the test type – either a whole herd test or a routine herd test – at random according to the proportion of tests recorded within the parish testing interval of the simulated herd (Figure S2 , Dataset S1 ). Whole herd tests specify that all bovines older than 6 weeks should be tested. We simulate this requirement by approximating the instantaneous proportion of the herd ineligible for testing to be (6/52) μ/N, where μ is the per-capita turnover of the herd. The number of animals tested with a WHT (X) is then sampled from a binomial distribution: There is greater variability in which non-breeding animals are tested during a routine herd test (RHT), and therefore in the proportion of the herd tested. In order to account for this we sample the proportion from a Cauchy distribution fitted to the empirical distribution from VetNet data by maximum likelihood (Figure S2 ) with scale parameter 0.0932 and shift parameter 0.494 (to 3 s.f.). The number of animals tested with a RHT (X) is then sampled from a binomial distribution as before with: The outcome of diagnostic tests within our model is determined by the set of parameters defining the sensitivity and specificity of the SICCT test at both the standard and severe interpretations (Table 1 ). In the field, the classification of reactors is based on cut-off values for the difference in reaction between avian and bovine tuberculin, with the cut-off value for a reactor changing with the severe and standard interpretation. As a consequence, tests can, and are, re-interpreted at the severe interpretation following the confirmation of a reactor animal after slaughter. In order to model the process of confirmation in a consistent fashion we must simulate the test outcome for each individual animal in the herd separately to ensure that the number of reactors at the severe interpretation is strictly greater than or equal to the number at the standard interpretation. Given X animals to test we sample them randomly (and uniformly) from each of the model compartments (S,O,R,I) to generate the number of animals from each compartment that are tested (XS,XO,XR,XI). For each (XS,XO,XR,XI) we sample a uniform random number and use the value to simulate the number of reactor animals at the standard (Standard Reactors) and severe (Severe Reactors) interpretations: For each XS: if ( ): Standard Reactors +1; if ( ) : Severe Reactors +1 For each Xo: if ( ) : Standard Reactors +1; Z+1 if ( ) : Severe Reactors +1 For each XR: if ( ) : Standard Reactors +1; Z+1 if ( ) : Severe Reactors +1 For each XI: if ( ) : Standard Reactors +1; Z+1 if ( ) : Severe Reactors +1 We must also keep track of the number of true reactors (Z) in order to simulate the number of confirmed reactors (C): Provided that the breakdown has not been previously confirmed and C = 0, then all reactors at the standard interpretation are removed from the herd. Otherwise, if the number of confirmed reactors , then the breakdown status is switched to confirmed (requiring an additional severe interpretation clear test to move back to the standard interpretation) and the reactors at the severe interpretation (including those from the current test) are removed. ABC-SMC implementation We use the ABC-SMC algorithm described in Toni et al. [19]. In essence the method replaces the calculation of a likelihood function with an approximation based on matching model simulations to the observed data using a set of goodness-of-fit metrics (in this case corresponding to a set of key epidemiological target measures). These target measures can be combined to produce a single metric (described in the next section). For a given set of parameters, a series of stochastic simulations from the model are produced using the algorithm described in the previous section. A simulation is said to “match” the observed data if the corresponding (simulated) metric value lies below a pre-determined threshold. A Monte Carlo estimate of the probability of matching can then be used as an approximation to the likelihood in an SMC framework [18], [19]. As the tolerance applied to the metric is reduced, the approximate (posterior) distribution should in principle converge towards the (true) posterior distribution. The algorithm begins by generating 10,000 particles - each “particle” corresponding to a set of model parameters – from a set of uniform proposal distributions (Table 2 ). For each particle we produce a binary Monte Carlo estimate for the probability of matching. This information is combined with the prior distributions for the parameters to produce a set of weights across the whole population of particles. The algorithm then proceeds through a series of repeated steps, whereby the population of particles is re-sampled from the previous weighted population and then each particle perturbed according to a set perturbation kernel. A new set of weights is then generated in a similar manner to before. The tolerance controlling the matching is reduced at each step until the predictive distributions from the simulated model generate an acceptable agreement with the target epidemiological measures. All parameters are log transformed and the perturbation kernel is uniform for each parameter (on this scale) , where is the range of the marginal distribution for parameter from the previous SMC round. At each successive round the threshold was reduced semi-automatically to the median value of the metric from the previous round. Heuristically the ABC-SMC procedure can be thought of as using goodness-of-fit criteria to inform the shape of the approximate posterior distributions, rather than the likelihood function. Prior assumptions Uniform prior distributions are applied to individual parameters and combinations of parameters to constrain their values to biologically relevant ranges. Probabilities are constrained to be in the interval [0,1] and rates are constrained to be positive. The sensitivity and specificity of the SICCT test are constrained to increase and decrease respectively under the severe interpretation and the probability of confirmation of reactors under routine slaughterhouse surveillance is assumed to be less than the probability of confirmation of reactors. All prior assumptions are intended to be uninformative, apart from the occult period for the SORI model and the upper bound (0.0003) placed on pFP. This value, equivalent to a lower bound on the specificity of the skin test at the standard definition of 99.97%, was obtained by calculating the value of pFP required to explain all of the unconfirmed breakdowns within VetNet data given the total number of animal tests. Target epidemiological measures Table 5 summarizes the target measures used to build our ABC metric. For each measure there are 6 empirical targets, by 2 values of confirmation status, by 3 values of PTI to give 36 target distributions/probabilities for each epidemiological measure. 10.1371/journal.pcbi.1002730.t005 Table 5 Epidemiological target measures for ABC. Description Type of Measure Number of bins per target distribution Weighting ( ) Breakdown Length Distribution (Days) [100,200,300,400,500,1000,2000] 7 1/7 Reactors at first test Distribution (Reactors) [1,2,3,4,5,10,47] 7 1/7 Reactors at VE-6M Distribution (Reactors) [1,2,3,4,5,10,47] 7 1/7 Reactors at VE-12M Distribution (Reactors) [1,2,3,4,5,10,47] 7 1/7 Total reactors removed within breakdown (until movement restrictions are lifted) Distribution (Reactors) [2,4,6,8,10,12,14,16,18,20,47] 11 1/11 Probability of recurrence within 6 months Probability 1 1 Probability of recurrence within 12 months Probability 1 1 Probability of recurrence within 24 months Probability 1 1 All of the target epidemiological measures for our final ABC-SMC scheme can be expressed either as probabilities or as (binned) probability distributions. This motivated the choice of an ABC metric based on the relative entropy, also known as the Kullback-Leibler divergence [34], that measures the distance between a proposed probability distribution (p) and a reference distribution (q): Two properties of the relative entropy should be noted: firstly, the relative entropy is asymmetric to the choice of reference distribution, with . We choose to remove the potential ambiguity introduced through the choice of reference distribution by using a symmetrical metric ( ): Secondly the relative entropy is undefined if any of the elements of the reference distribution qi = 0. We numerically approximate the distribution of each of our target measures through histograms, with bin-sizes chosen to capture the range of observed values within VetNet data. Where we are free to choose appropriate bin sizes for the empirical distributions (q) such that we avoid any empty bins, we cannot ensure the same for the proposed distributions (p) generated from model simulations. To ensure that our metric is always defined, we add 1 to every bin of our empirical and simulated histograms. For each proposed set of parameters (particle) we simulated a fixed number of realizations of the model (500) at the midpoint of each of the 6 herd-size histogram bins [30,90,150,210,270,330] for PTI 1,2 and 4 and generate a set of j proposed distributions ( ) for each corresponding target measure ( ). We use a superscript to indicate the jth distribution with elements indexed by i. We calculate the distance between the proposed and target distributions using the symmetric metric ( ) introduced above. The maximum value of this metric will increase with the number of histogram bins associated with that target distribution. In order to ensure that the overall metric places a more equal weight on each epidemiological measure we weight the contributions from each target measure proportionally to the total number of bins forming that measure ( ) (Table 5 ): Supporting Information Dataset S1 Target and auxiliary distributions necessary for simulation and parameterization of models using ABC-SMC. (ZIP) Click here for additional data file. Figure S1 Distribution of times between scheduled whole herd tests in GB (2003–2005). Farmers are responsible for scheduling tests as close as possible to the statutory intervals. Historically, this has lead to variation in the time between tests that we quantify for our study period (2003–2005). The frequency of routine herd tests is determined by the current parish-testing interval (PTI) for a herd. However, the time since the previous whole herd test is also determined by the historical testing intervals for the parish and other epidemiological factors. In practice the time since the previous surveillance test for breakdown herds in PTI 1, 2 and 4 is distributed with the greatest variation seen in PTI 2 (left). Short interval tests (SIT) must be carried out at least 60 days after the last whole herd test leading to a skewed distribution where test intervals are more likely to be late than early (middle). Likewise the follow up tests after a breakdown, that must be scheduled at intervals of at least 6 and 12 months respectively (VE-6 M, VE-12 M) are skewed to be late (right). (PDF) Click here for additional data file. Figure S2 Proportion of animals tested in routine surveillance tests (2003–2005). The number of animals within a herd that are tested during routine surveillance varies depending on the demographic structure of the herd and the perceived epidemiological risk. PTI 1 herds should receive a whole herd test (WHT) where all bovines older than 6 weeks are tested. In PTI 2, 3 and 4 a routine herd test (RHT) may be carried out where there is greater discretion as to which animals are tested based on perceived epidemiological risk. As a consequence the proportion of animals tested is smaller and more variable for RHTs (Right) as compared to WHTs (Middle). The proportion of breakdowns reported as being disclosed by WHTs and RHTs also varies by PTI (Left table). The small proportion of breakdowns in PTI 1 initiated by a RHT, despite WHTs being mandated in these herds, stem from herds whose PTI was updated retrospectively after disclosure. Likewise despite the majority of tests in PTI 4 being RHTs, a higher proportion of breakdown herds are initiated by WHTs, applied when there was a perceived higher risk or consequence of infection (e.g. in dealer or milk retailer herds). (PDF) Click here for additional data file. Figure S3 Annual per bovine rate of turnover in breakdown herds. We define turnover as the time-averaged rate at which animals move into and are removed from a herd. A representative distribution of turnover rates for breakdown herds was calculated from CTS data from 1st January 2003 – 1st January 2005 for all breakdown herds with start dates in 2004. Note that since the CTS data only records movements at the holding (CPH) rather than the herd (CPHH) level, this is an indirect measure corresponding to the annual per capita rate of movement of bovines through the CPH associated with a breakdown herd. Turnover was calculated three ways: using movements into a CPH (“On Movements”, red line), movements out of a CPH (“Off Movements”, green line) and the combined rate of both types of movements (black line, points). The median number of “Off” movements is slightly smaller with than for on-movements consistent with the increase in herd size nationally over the period. The “All movements” estimate is used as the empirical distribution for the within-herd model. (PDF) Click here for additional data file. Figure S4 Distribution of Breakdown Herd Sizes. Smoothed density of the (maximum) herd size during a breakdown for our study population. There is some variation in the size of herds with PTI, with a longer tail of herds beyond our cutoff value of 360 (vertical dashed line) for PTI 4 herds. (PDF) Click here for additional data file. Figure S5 Estimated Parameter Distributions for SORI model. Distributions of parameters consistent with the persistence measures and reactor distributions estimated from VetNet data (Figure S6 ). The severe values of sensitivity and specificity are constrained to be greater than and less than their respective values at the standard interpretation. Likewise the probability of infected animals being detected by routine slaughterhouse surveillance is assumed to be less than or equal to the probability of confirmation. All parameters are constrained to be positive, with probabilities and the density dependent parameter q further constrained to be less than or equal to 1. (PDF) Click here for additional data file. Figure S6 Target measures for the within-herd persistence of bTB (2003–2005) and predictive distributions for SORI model. The within-herd persistence of bTB in GB as measured by the probability of breakdowns being prolonged (duration of greater than 240 days) or recurrent within a 6, 12 and 24 month horizon. The relationship of each target measure is plotted against herd size, with breakdowns further stratified by parish testing interval (PTI 1 top row, PTI 2 middle row, PTI 4 bottom row) and confirmation status (confirmed breakdowns: green, circles, unconfirmed breakdowns: magenta squares). Target measures are calculated from breakdowns trigged within 2003–2005 by a routine surveillance test (VE-WHT, VE-WHT2, VE-RHT, VE-SLH). The probability of confirmation varies between PTI, as does the proportion of confirmed breakdowns initiated by a slaughterhouse case (white diamonds) and the mean number of reactors reported at the disclosing test. Uncertainty in each (mean) target observation (thick lines) is illustrated by an envelope (thin lines) of ±1.96 standard errors around the mean. Predictive distributions for each of these target measures from the finalized within-herd transmission model are plotted as shaded density strips where the intensity of color is proportional to the probability density at that point. (PDF) Click here for additional data file. Figure S7 Target measures for the within-herd persistence of bTB (2006–2008) and predictive distributions for SORI model. The within-herd persistence of bTB in GB as measured by the probability of breakdowns being prolonged (duration of greater than 240 days) or recurrent within a 6, 12 and 24 month horizon. The relationship of each target measure is plotted against herd size, with breakdowns further stratified by parish testing interval (PTI 1 top row, PTI 2 middle row, PTI 4 bottom row) and confirmation status (confirmed breakdowns: green, circles, unconfirmed breakdowns: magenta squares). Target measures are calculated from breakdowns trigged within 2003–2005 by a routine surveillance test (VE-WHT, VE-WHT2, VE-RHT, VE-SLH). The probability of confirmation varies between PTI, as does the proportion of confirmed breakdowns initiated by a slaughterhouse case (white diamonds) and the mean number of reactors reported at the disclosing test. Uncertainty in each (mean) target observation (thick lines) is illustrated by an envelope (thin lines) of ±1.96 standard errors around the mean. Predictive distributions for each of these target measures from the finalized within-herd transmission model are plotted as shaded density strips where the intensity of color is proportional to the probability density at that point. (PDF) Click here for additional data file. Figure S8 Burden remaining after resolution of a breakdown using SORI model. redictive distributions for the probability of at least one infectious bovine remaining within a herd after a breakdown is resolved as a function of herd size, classified by PTI (1,2,4 left to right) and confirmation status (Green circles confirmed, magenta squares unconfirmed). Predictive distributions are plotted as shaded density strips where the intensity of shading is proportional to the probability density at that point. Solid lines, and points for each herd-size category, indicate the median of the predictive distribution to aid comparison between confirmed and unconfirmed breakdowns. (PDF) Click here for additional data file. Figure S9 Estimated Parameter Distributions from SOR model. Distributions of parameters consistent with the persistence measures and reactor distributions estimated from VetNet data (Figure S7 ). The severe values of sensitivity and specificity are constrained to be greater than and less than their respective values at the standard interpretation. Likewise the probability of infected animals being detected by routine slaughterhouse surveillance is assumed to be less than or equal to the probability of confirmation. All parameters are constrained to be positive, with probabilities and the density dependent parameter q further constrained to be less than or equal to 1. (PDF) Click here for additional data file. Figure S10 Target measures for the within-herd persistence of bTB (2003–2005) and predictive distributions for within-herd transmission model using SOR model. The within-herd persistence of bTB in GB as measured by the probability of breakdowns being prolonged (duration of greater than 240 days) or recurrent within a 6, 12 and 24 month horizon. The relationship of each target measure is plotted against herd size, with breakdowns further stratified by parish testing interval (PTI 1 top row, PTI 2 middle row, PTI 4 bottom row) and confirmation status (confirmed breakdowns: green, circles, unconfirmed breakdowns: magenta squares). Target measures are calculated from breakdowns trigged within 2003–2005 by a routine surveillance test (VE-WHT, VE-WHT2, VE-RHT, VE-SLH). The probability of confirmation varies between PTI, as does the proportion of confirmed breakdowns initiated by a slaughterhouse case (white diamonds) and the mean number of reactors reported at the disclosing test. Uncertainty in each (mean) target observation (thick lines) is illustrated by an envelope (thin lines) of ±1.96 standard errors around the mean. Predictive distributions for each of these target measures from the finalized within-herd transmission model are plotted as shaded density strips where the intensity of color is proportional to the probability density at that point. (PDF) Click here for additional data file. Figure S11 Target measures for the within-herd persistence of bTB (2006–2008) and predictive distributions for within-herd transmission model using SOR model. The within-herd persistence of bTB in GB as measured by the probability of breakdowns being prolonged (duration of greater than 240 days) or recurrent within a 6, 12 and 24 month horizon. The relationship of each target measure is plotted against herd size, with breakdowns further stratified by parish testing interval (PTI 1 top row, PTI 2 middle row, PTI 4 bottom row) and confirmation status (confirmed breakdowns: green, circles, unconfirmed breakdowns: magenta squares). Target measures are calculated from breakdowns trigged within 2003–2005 by a routine surveillance test (VE-WHT, VE-WHT2, VE-RHT, VE-SLH). The probability of confirmation varies between PTI, as does the proportion of confirmed breakdowns initiated by a slaughterhouse case (white diamonds) and the mean number of reactors reported at the disclosing test. Uncertainty in each (mean) target observation (thick lines) is illustrated by an envelope (thin lines) of ±1.96 standard errors around the mean. Predictive distributions for each of these target measures from the finalized within-herd transmission model are plotted as shaded density strips where the intensity of color is proportional to the probability density at that point. (PDF) Click here for additional data file. Figure S12 Burden remaining after resolution of a breakdown using SOR model. Predictive distributions for the probability of at least one infectious bovine remaining within a herd after a breakdown is resolved as a function of herd size, classified by PTI (1,2,4 left to right) and confirmation status (Green circles confirmed, magenta squares unconfirmed). Predictive distributions are plotted as shaded density strips where the intensity of shading is proportional to the probability density at that point. Solid lines, and points for each herd-size category indicate the median of predictive distribution to aid comparison between confirmed and unconfirmed breakdowns. (PDF) Click here for additional data file. Figure S13 Impact of herd-level interventions on probability of recurrence within 24 months. Change in the probability of a herd experiencing a recurrent breakdown after application of a ‘perfect’ test (left column) or perfect isolation (right column). The perfect test is assumed to have 100% sensitivity and specificity and no occult period. Perfect isolation corresponds to setting the extrinsic infectious pressure to zero at the end of a breakdown ( ). Plotted values correspond to the average % difference in the probability of recurrence relative to the fitted SORI (Panel A) and SOR models (Panel B). Separate series are plotted for herds in PTI 1 (lime green circles), 2 (magenta squares) & 4 (sky blue diamonds). Predictive distributions are plotted as shaded density strips where the intensity of shading is proportional to the probability density at that point, with the mean of the predictive distribution plotted as a solid line. (PDF) Click here for additional data file.

Journal ID (nlm-ta): Proc Biol Sci

Journal ID (iso-abbrev): Proc. Biol. Sci

Journal ID (publisher-id): RSPB

Journal ID (hwp): royprsb

Title:
Proceedings of the Royal Society B: Biological Sciences

Publisher:
The Royal Society

ISSN
(Print):
0962-8452

ISSN
(Electronic):
1471-2954

Publication date
(Print):
22
May
2014

Publication date PMC-release: 22
May
2014

Volume: 281

Issue: 1783

e-mail:
rowland.kao@
123456glasgow.ac.uk

DOI: 10.1098/rspb.2014.0248

PMC ID: 3996616

PubMed ID: 24718762

Copyright statement:

License:

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.

Subject:
1001

Subject:
44

Subject:
87

Subject:
Research Articles

cover-date May 22, 2014

Data availability:

ScienceOpen disciplines: Life sciences