INTRODUCTION
Campylobacter spp. continue to be a major problem in large parts of the world, including Denmark, being one of the most common causes of human bacterial gastroenteritis [1–3]. The species most frequently associated with human disease are Campylobacter jejuni and C. coli. The proportion between the species varies between countries [2]. In Denmark, most human campylobacteriosis cases are caused by C. jejuni (∼96%) [Reference Nielsen, Engberg and Madsen4].
Campylobacter spp. have been detected in many sources and are considered to be widespread in production animals and in the environment [Reference Miller, Mandrell, Ketley and Konkel5]. In many countries, travel is considered a major risk factor for acquiring campylobacteriosis [6]. In Denmark, travel is estimated to account for about one third of human campylobacteriosis cases [Reference Ethelberg7]. Broiler chicken meat is recognized as the largest single source of foodborne campylobacteriosis cases [6, Reference Wingstrand8]. Consequently, several countries have already established action plans against Campylobacter in the broiler production chain [9]. A number of initiatives have also been initiated in Denmark to reduce the Campylobacter burden in broiler production [Reference Rosenquist10].
The Danish strategy to control Campylobacter has had a positive effect; reducing the prevalence in broiler flocks and broiler meat. A small decrease in the number of human cases has also been observed. The explanation for the effect on human cases not being more significant is probably due to other factors counterbalancing the effect of the implemented interventions [Reference Rosenquist10]. In particular, imported broiler meat and sources of infection other than broilers are assumed to have influenced the limited effect on humans. Consequently, it could be interesting to evaluate what other sources could potentially add to the total number of human campylobacteriosis cases. With regard to the Danish Salmonella situation, the development of a source-attribution model has proved to be an important tool for risk managers in order to implement targeted interventions, resulting in a significant reduction in human salmonellosis [Reference Wegener11]. It would be of benefit if a similar tool could be developed for Campylobacter for prioritization of public health resources and source-specific implementation of control measures.
Several methods of source attribution are available [Reference Pires12]; however, some are more suitable to apply for Campylobacter than others. For example, analysis of outbreak data makes little sense, as recognized outbreaks caused by Campylobacter are considered to be rare in the EU [6]. Interesting results have been shown using the microbial subtyping approach; e.g. in New Zealand and England [Reference Wilson13, Reference Mullner14]. The models used were modifications of the original Danish Salmonella model and the newly developed Asymmetric Island (AI) model [Reference Wilson13, Reference Mullner14]. Both studies found poultry to be the main reservoir, and broilers to be the principal source of human campylobacteriosis caused by C. jejuni. Second to broilers came the cattle reservoir.
A challenge in source attribution for Campylobacter cases based on microbial subtyping is that individual subtypes appear to be widespread between sources. The discriminatory ability of the commonly used typing schemes does not at present allow for a very distinct separation between sources, which complicates the attribution of human cases to possible sources. To be able to attribute human cases to different sources, strain genetic diversity between source groups is essential. Future identification of more source-specific markers would be of great importance for optimizing the apportioning of cases.
With the perspective of the potential benefit from adoption of a Campylobacter source-attribution tool with regard to risk management, the aim of the present study was to assess the Danish Campylobacter situation using two different source-attribution approaches based on multilocus sequence typing (MLST) of C. jejuni isolates; each approach identifying the primary reservoirs of domestically acquired cases and cases with unknown travel history. Furthermore, the effect of adding an extra typing parameter (the flaA gene) for additional discriminatory power was explored.
MATERIALS AND METHODS
Isolates
C. jejuni isolates included in the study were collected in 2007 and 2008 and originated from various projects including the EU baseline study of broiler carcasses [15], the national Campylobacter surveillance for broiler meat (at slaughter plants and at retail), and the national surveillance for antimicrobial resistance (including various animal species) [16]. Some isolates were collected for this study only.
In total, 406 human isolates were collected mainly from three regions of Denmark: Northern Jutland, Funen, and Zealand. Of the human isolates, 246 were reported as domestic (i.e. acquired in Denmark), 109 were reported as travel related, and 51 had unknown travel history. Cases were reported as related to travel, if a person in the 7-day period prior to disease onset had stayed a minimum of one night in any country other than Denmark [Reference Ethelberg7]. All human isolates were from cases characterized as sporadic, i.e. not associated with any known outbreaks. Isolates from six putative sources were included: Danish broiler meat (185 isolates collected at retail), imported broiler meat (137 isolates collected at retail), turkey meat (96 isolates collected at retail), duck meat (70 isolates collected at retail), cattle (171 isolates from faeces), pig/pork (three isolates from caecum and one isolate collected at retail).
MLST, flaA typing and antimicrobial susceptibility testing
All isolates were characterized by MLST, flaA type and antimicrobial susceptibility.
MLST was performed according to the scheme described by Dingle et al. [Reference Dingle17], which sequenced seven housekeeping genes. For each house keeping gene, the different sequences were assigned as distinct alleles (by assignment of a single number) and, for each isolate, the alleles at each of the seven loci defined the allelic profile, which determined the sequence type (ST).
The flaA type was determined by sequencing of a short variable region of the flaA gene as described previously [Reference Meinersmann18]. The PubMLST database was used for identification of profiles (http://pubmlst.org/campylobacter/) [Reference Jolley and Maiden19].
Antimicrobial susceptibility was determined as described in DANMAP 2008 [16]. In brief, the antimicrobial susceptibility testing was performed by microbroth dilution minimum inhibitory concentration (MIC) using the Sensititre system (Trek Diagnostic Systems Ltd, UK). Inoculation and incubation procedures were in accordance with CLSI guidelines. MIC values were interpreted using EUCAST epidemiological cut-off values. Susceptibility was determined for the following antimicrobial agents: chloramphenicol, ciprofloxacin, erythromycin, gentamicin, nalidixic acid, streptomycin, and tetracycline. The result of the testing was coded as resistant (R) or sensitive (S). An isolate was categorized as resistant if it was resistant to at least one of the specified agents.
Analysis of molecular variance (AMOVA)
The genetic structure differentiation between groups was tested by AMOVA [Reference Excoffier, Smouse and Quattro20]. AMOVA is a method of testing population differentiation directly from molecular data, based on Euclidean distance metrics. Variance components are used to calculate statistics called Φ-statistics, summarizing the degree of differentiation between population groups. The genetic distance between a pair of isolates was defined as the number of loci, out of seven, at which they differed. AMOVA was performed using Arlequin software version 3.11. The number of permutations for significance was 999 and the level of significance was set to 0·05. The genetic diversity between groups (ΦGT) was tested by estimation of the pairwise difference. The lower the value of ΦGT, the lesser the variation between groups. A significant P value indicates statistical difference between groups.
Source-attribution modelling
Two models were used to attribute human cases: the AI model developed by Wilson et al. [Reference Wilson13] and a model modified after the Danish Salmonella attribution model [Reference Hald21]. The second model will henceforth be designated the CAMSA (Campylobacter source attribution) model.
The AI model apportions domestic cases to different putative sources defined by relatedness to groups comprising isolates collected from the respective source. This model is an evolutionary model taking into account mutation, recombination and migration rates. The model was used without modifications. Convergence of the model was monitored by multiple runs ensuring equal results; 100 000 iterations were run without thinning. The model was run using the freely available software iSource (downloaded from the website http://www.danielwilson.me.uk/software.html). The input for the model was the allelic profile of each isolate in the source groups. For attribution, domestic human cases and human cases with unknown travel history were used as input. In addition to running the model with the animal food sources, it was further explored by running the model with the inclusion of a ‘source’ called travel. Based on the 109 human cases categorized as travel related, an additional source was created. This was done in order to investigate the possibility of allowing cases with unknown travel history to be attributed to travel, and not making the assumption that these cases are all domestic or choosing not to attribute the cases at all.
The CAMSA model was adapted after the Danish Salmonella attribution model [Reference Hald21]. This model apportions human cases using a Bayesian framework. The modelling was based on the occurrence of types in included sources, combined with the amount of the food source available for consumption and two factors regarding (1) the type-specific ability to cause infection, and (2) the source-specific ability to serve as a vehicle for the types. The equation used to estimate the expected number of human cases was:
$$\lambda _{ij} = M_j p_{ij} q_i a_j, $$where λ ij is the expected number of cases of type i from source j, M j is the amount of source j available for consumption, p ij is the number of isolates of type i in source j (modification from the original model), q i is the ST-dependent factor, and a j the source-dependent factor. The equation represents a multi-parameter prior, where q i and a j are parameters of unknown value. These parameters were included as distributions; a hierarchical prior (modification from the original model) and a uniform prior, respectively. The use of a hierarchical prior was adapted after Mullner [Reference Müllner22], using a lognormal distribution N(0, τ). The prior distribution for τ was gamma (0·01, 0·01). The source-dependent factor, a j, was assumed equal for Danish-produced meat and imported chicken meat. A Markov Chain Monte Carlo simulation, specifically the Gibbs sampler, was applied to compute the posterior distributions for a j and q i. Five independent Markov chains of 40 000 iterations were run. Convergence was monitored using methods described previously [Reference Gelman and Rubin23]. The model was run in WinBUGS version 1.4. The input for the model was ST. The ten STs most frequently found in humans were further differentiated by supplementing the STs with information of antimicrobial susceptibility. Domestic human cases and human cases with unknown travel history were inputted for attribution. Cases categorized as related to travel were directly assigned as this.
The main difference between the two models is the principle behind attribution of cases and the factors which the models account for. The AI model estimates the probability of each source for each human isolate, accounting for evolutionary relationship, meaning that all human cases will be apportioned to sources. The CAMSA model estimates the expected number of cases per source based on comparison of the observed number of cases caused by a specific type with the occurrence of types in specified sources, weighted by amount of food source available for consumption and accounting for type- and source-specific ability to cause disease. The CAMSA model requires an exact type match between isolates from humans and sources. Therefore, human types that are not found in any source are referred to the group ‘unknown’.
Both models were run at two different discriminatory levels; level 1: differentiation based on the seven housekeeping genes (considered as the basic model), and level 2: differentiation based on the seven housekeeping genes + the flaA type.
RESULTS
The results of the AMOVA analysis (Table 1) showed that the genetic differentiation between source categories was statistically significant for almost all groups. The genetic difference was not significant between the group comprising domestic human cases and the group of human cases with unknown travel history. This indicates that human cases with unknown travel history were more similar to the group of domestic cases compared to the group of travel-related cases. The pig group was not statistically different from three of the groups, which was probably due to the very small number of isolates (n = 4).
Table 1. Genetic differentiation between groups
Values below the diagonal are the estimated pairwise differences between groups (ΦGT); values above the diagonal are the associated P values.
Non-significant Φ-statistics appear in bold.
The output from the basic AI model is presented in Table 2. The model attributed most cases to the broiler chicken reservoir; 52% [95% credibility interval (CrI) 37–67)] to Danish chicken and 17% (95% CrI 3–33) to imported chicken. The second most important reservoir in relation to human illness was cattle (17%, 95% CrI 7–28). The uncertainty around the attribution estimates ranged widely. This is influenced by the fact that most cases could not stringently be assigned to one particular source. In particular, the distinction between broiler chicken and cattle was vague (Fig. 1a). The inclusion of travel as a ‘source’ for cases with unknown travel history in the AI model only changed the attribution of cases slightly. Most cases were still attributed to broiler chicken and cattle (Table 2). Abouty 11% (95% CrI 1–29) of cases with unknown travel history were attributed to travel, i.e. 2% of all apportioned cases.
Fig. 1. Probability of each human case belonging to each of the included sources (results from the basic Asymmetric Island model). (a) Modelling based on multilocus sequence typing (MLST), (b) modelling based on MLST + flaA. The probability is depicted by colour coding: cattle (dark blue), Danish chicken (red), imported chicken (yellow), turkey (green), duck (cyan), pork (pink).
Table 2. Attribution of human cases with the Asymmetric Island model, excluding and including travel as a ‘source’
CrI, Credibility interval.
Values given are proportion of cases attributable to the specific source and corresponding uncertainty (95% CrI).
* Based on the separate models for domestic cases and cases without travel history.
For the CAMSA model, the ST-dependent factor (q i) was fairly equal between STs. Only the estimate for one ST (ST4811) tended to be higher than the rest indicating this type to result in relatively more human cases compared to the other STs. The food-related factor (a j) for cattle tended to be higher compared to the food-related factor for other sources. The output of the CAMSA model was very similar to the results of the AI model (Table 3). The primary reservoir being broiler chicken, comprising Danish chicken 38% (95% CrI 28–47) and imported chicken 14% (95% CrI 10–18). Sixteen percent of cases were attributed to cattle. About one fifth of cases could not be assigned to any of the sources.
Table 3. Proportion of cases attributable to the specific source and corresponding uncertainty (95% CrI)
CrI, Credibility interval.
Attribution of human cases (domestic and without travel history) with basic models; differentiation based on MLST (input 1) and differentiation based on MLST + typed flaA gene (input 2)
Both models attributed the majority of cases to the broiler chicken reservoir. Therefore Danish broiler chicken was found to be the largest contributor compared to imported chicken. The cattle reservoir was found to be the second most important in relation to human campylobacteriosis (Table 3).
Inclusion of flaA gene sequences caused different impacts on the attribution outputs from the two models (Table 3). For the AI model, the inclusion of the flaA gene boosted the proportion of cases attributed to chicken at the expense of the proportion attributable to cattle. The uncertainty around the attribution estimates still ranged widely, but diminished slightly following the inclusion of flaA. The inclusion of flaA resulted in a larger proportion of cases more strictly associated with chicken (Fig. 1b) compared to the basic model (Fig. 1a). For the CAMSA model a larger proportion of cases fell in the group ‘Unknown’ including a proportion from each group except travel. The proportion of cases attributed to cattle decreased slightly more compared to the other groups.
DISCUSSION
Both models agreed in recognizing broiler chicken as the primary source of human campylobacteriosis. This supports the original hypothesis of chicken being the most important single source of human campylobacteriosis. A higher proportion of cases were apportioned to Danish chicken compared to imported chicken. Several factors can explain this. First, the risk from Danish chicken might actually be higher compared to imported chicken. Second, the Danish broiler chicken reservoir comprises of more transmission routes compared to imported chicken, as the production is taking place in the country of concern. Besides meat, animal contact through, e.g. occupation, may also be a risk factor. Third, some of the cases assigned to Danish chicken might actually belong to another reservoir, as the STs comprising the Danish chicken reservoir are more closely related to the other sources than the STs comprising the source of imported chicken. A Danish case-control study found a population-attributable risk from fresh chilled chicken of 24% of domestically acquired cases [Reference Wingstrand8]. A Scientific Opinion from EFSA suggests that handling, preparation and consumption account for only 20–30% of human cases while 50–80% may be attributed to the broiler chicken reservoir as a whole [6]. This agrees very well with the Danish studies; attributing about 50–60% of cases to the broiler chicken reservoir as a whole with about half of these acquired from chilled chicken meat [Reference Wingstrand8]. It is not possible from the source attribution to estimate the proportion of cases caused by handling, preparation and consumption of Danish chicken meat because the isolates collected represent all transmission routes from the chicken reservoir to the consumer. This is in contrast to imported chicken meat, where transmissions routes prior to packaging are of no risk to the Danish population. However, we know that the proportion of Danish/imported meat available for sale is 60/40 (in 2008). Assuming no difference in ability of infection between types and combining source estimates and consumption data, we could infer, that human cases caused by Danish chicken meat would be about 21% of all cases.
Cattle were found to be the second most important source. High C. jejuni prevalence has been reported in cattle [24, 25]; however, very low occurrence has been found in beef [24, 25]. In light of the finding that cattle were the second most important source, other exposure routes as well as meat should also be considered. This would agree with the results from a Dutch comparative exposure assessment ranking farm animal contact higher compared to beef with regard to importance of transmission [Reference Evers26]. The Campylobacter occurrence in meat from ducks is high in Denmark [1, 27]. However, only a few cases were attributed to this reservoir by the models. The reasons for this observation is probably the way of handling and preparing this product and that the consumption of this type of meat is less compared to chicken meat. Traditionally, ducks are prepared whole and hours before garnishing, reducing the risk of cross-contamination considerably. The proportion of cases attributed to turkey, being smaller than broiler chicken and higher than duck, fits with the proportionate consumption of turkey meat being less than chicken and higher than duck.
The addition of travel-related cases as a ‘source’ in the AI model estimated ∼11% of cases with an unknown travel history to this group. This estimate corresponds to what was found with the CAMSA model and also the result of the AMOVA analysis; finding a greater similarity between cases without travel history and domestic cases compared to travel-related cases. The introduction of travel in the AI model is a way of coping with the fact that information about travel history might not always be available for all human cases. This scenario is probably the situation in many countries as it is not always feasible to obtain this information. As travel is considered an important risk factor for acquiring campylobacteriosis, the modelling could benefit from being able to handle this in case of missing travel information from some human cases.
The assumption that human cases, with travel history 7 days prior to onset of illness, have been infected abroad might not be true for every case. In Denmark, cases are reported as related to travel if a person on the day of, or 7 days prior to, illness onset has stayed for a minimum of one night in any country other than Denmark [Reference Ethelberg7]. However, there is still a possibility that the infection has been acquired in Denmark. Consequently, the estimate might be too high. In addition, the risk of acquiring campylobacteriosis varies between countries [Reference Ethelberg7, 28], which is not considered in the model at present.
The isolates representing each source were collected through different monitoring programmes and surveys, which varied in sampling size, sampling period, etc. To achieve the best possible attribution, the majority of the variation of STs within sources should be covered. Rarefaction analysis showed that this was not the case (with the exception of cattle) (data not shown). In order to show the full diversity of STs, a large number of additional samples would be required which would be extremely costly. We used national data collected within the project period only in order to limit potential biases due to both the possible occurrence of different STs in different countries and changes in ST frequencies over time. This, however, results in lower numbers of isolates available for modelling. A larger number of isolates could potentially influence the results, if this would change the mix of STs.
Thus far, this approach for Campylobacter has only been applied to studies with a ‘one nation approach’; however, studies are in the pipeline comparing STs from different countries and regions.
Genetic diversity was found between isolates collected from chicken meat at retail and live broilers (Table 1). This might reflect that we have not covered the whole genetic diversity or that a potential natural selection of types through the processing chain might occur. We chose to use isolates from retail chicken to represent the Danish broiler chicken reservoir. The use of retail isolates compared to isolates from live broilers might to some degree reduce the bias of environmental types in the Danish broiler chicken reservoir.
Data from other sources were sampled but not included in the study. Samples were collected from 125 pet dogs and cats; however, only one C. jejuni isolate was found. It was decided not to include this isolate. Furthermore, samples from fresh water streams, petting zoo goats, dairy cattle, and meat from lamb were collected but unfortunately lost. We did not collect samples from tap water, as this is not considered a source of sporadic cases in Denmark. The influence of these decisions and events may have resulted in less accuracy in assigning cases, as some cases might rightfully belong to the sources not included in the model. In theory, a larger part of the group ‘unknown’ from the CAMSA model may have been explained by the addition of more sources and the AI model might have attributed cases differently. Almost all tap water in Denmark is obtained from the groundwater reserve and is filtered during diffusing from the surface through the soil layers. This is considered to be a sufficient hygienic barrier and no further treatment is done. Water quality is monitored by extensive testing for indicator bacteria, but only sparse data exists for the actual occurrence of Campylobacter [Reference Jeppesen and Guldbæk29] and, furthermore, the relationship between Campylobacter and the indicator bacteria have not been confirmed. Waterborne outbreaks caused by Campylobacter have been observed in Denmark as result of contamination of the tap water supply and contamination of costal water after heavy rainfall [1]; however, as a source of sporadic cases, drinking water is disregarded. Puppies, compared to older dogs have been found to excrete C. jejuni in Denmark as well as in other countries [Reference Hald30–Reference Wieland32]. Therefore, young pets might pose a risk, which was not reflected by our sampling, as our sampling was not targeted towards puppies and kittens, but dogs and cats in general. Finally, a study by Evers et al. [Reference Evers26] demonstrated the potential importance of petting zoo animals. Further, there might remain other, unrecognized sources.
Hence, the inclusion of all relevant reservoirs in modelling is important for the ‘correct’ attribution of cases. Inclusion of more sources would probably not reduce the uncertainty in the attribution estimates, because of the great variation in STs within sources and the overlap in genotypes between sources would still exist and result in lack of clear separation. Even though MLST is highly discriminatory in typing C. jejuni, MLST cannot distinguish very clearly between animal reservoirs. Numerous types are detected in several reservoirs, i.e. ST21 and ST45. A more nuanced attribution of cases could be obtained in case of identification of more source-specific attributes. In addition to performing source attribution based on only the MLST data, and resistance profiles for the CAMSA model, the potential of adding another attribute, in this case the sequenced flaA gene, was explored. For the CAMSA model, this addition decreased the number of cases that the model was able to assign. This model seeks exact matches between human and source types, thus the added discriminatory power resulted in fewer matches between cases and sources. The proportion of cases attributed to the different sources decreased in approximately equal magnitude, suggesting that additional sampling of sources are needed to cover the large variation in STs. Only the proportion assigned to cattle decreased slightly more. For the AI model (Fig. 1a, b) the discriminatory power was increased by the addition of the flaA gene. Especially with regard to differentiation between Danish broiler chicken and cattle, the addition appeared valuable.
There are pros and cons in using either of the attribution models. Both models need a considerable amount of data. The AI model infers the apportioning of cases based on the available sources, accordingly producing estimates for only the implicated sources. This might skew the results if not all putative sources are represented in the data. The CAMSA model, on the other hand, has a category (unknown) for the cases that cannot be attributed to one of the sources in the model, however it requires an exact match of ST to assign cases to a source. Both ways of handling the data can be positive as well as negative. To attribute every case to a source might wrongly boost some categories if not all putative sources are included in the model; however, to only attribute cases with an exact ST match might be too stringent. For Campylobacter in particular, a close relatedness can be sufficient to reduce the large number of samples per source that should be collected to cover the great variation of types encompassed within the sources. Overall, the two models were considered to complement each other.
A prerequisite for using typing data for source attribution is that there is some degree of association between types and sources. A linkage between specific STs and source has been demonstrated [Reference Ogden33, Reference Sheppard34]. However, the biology of C. jejuni limits the chances for this linkage to be stable and universal. It has been shown that C. jejuni is genetically diverse, with a weakly clonal population structure, and that intra- and inter-species horizontal genetic exchange is common [Reference Dingle17]. Compared to Salmonella, the weaker population structure of Campylobacter will probably make it more difficult to obtain accurate results from source-attribution models that are primarily based on typing data obtained by any typing method.
The use of a reservoir model might produce different results than, e.g. a model attributing cases at the point of purchase or consumption (i.e. comparative exposure modelling). The application of a reservoir model enables evaluation of the problem at its root in contrast to source attribution at the point of purchase, which allocates cases to transmission routes. It could be of particular interest to try to combine the source attribution using microbial subtyping with comparative exposure assessment, thereby adding information regarding the potential transmission routes of Campylobacter.
Substantial work has been performed in New Zealand in the fight against Campylobacter. Included in this work was the use of source attribution for identification of specific sources for human campylobacteriosis (caused by C. jejuni). The initial source-attribution work found that 58–76% of human cases were caused by poultry sources [Reference Mullner14]. Similar results were found in England, estimating chicken as a source of human illness in 51–62% of cases [Reference Wilson13]. These findings are similar to the results from the present work on Campylobacter source attribution in Denmark. After implementation of a comprehensive intervention strategy, an evaluating source-attribution work revealed a decline in the number of cases attributed to poultry in New Zealand by 74% [Reference Sears35]. No marked decrease was observed in human cases in Denmark after the implementation of Campylobacter intervention strategies. The adaption of source attribution could be a valuable tool in assessing present and future intervention strategies. When comparing the results obtained in New Zealand, the difference in food trade patterns should be considered; New Zealand has no import of chicken meat compared to Denmark where ∼40% of the meat available for consumption is imported (in 2008). This of course affects the potential of national action plans to be effective.
Other studies have found that children aged <5 years are becoming ill from different sources compared to older population groups, and also living in rural compared to urban areas influences the epidemiology [Reference Mullner36, Reference Strachan37]. In future studies, we will look at the distribution in age groups and rural vs. urban areas.
Both models applied MLST data for source attribution and indicated that the major burden of human campylobacteriosis in Denmark originates from the broiler chicken reservoir. This was further emphasized by applying additional discriminatory power to the models by including flaA subtypes. Using source attribution as basis for national interventions for Campylobacter has proven effective in New Zealand. The results of the present study can be useful for future risk management decisions related to the control of Campylobacter in Denmark.
ACKNOWLEDGEMENTS
The Ministry of Food, Agriculture and Fisheries, The Danish AgriFish Agency funded the study. The authors gratefully acknowledge Sara Monteiro Pires, DTU for the introduction to the Danish Salmonella model and support in setting up the CAMSA model in WinBUGS and Jørgen Legind-Hansen for technical advice.
DECLARATION OF INTEREST
None.
