Detection probability models for bacteria, and how to obtain them from heterogeneous spiking data. An application to Bacillus anthracis

International Journal of Food Microbiology 241 (2017) 78–88

Contents lists available at ScienceDirect

International Journal of Food Microbiology journal homepage: www.elsevier.com/locate/ijfoodmicro

Detection probability models for bacteria, and how to obtain them from heterogeneous spiking data. An application to Bacillus anthracis Ronny Hedell a,b,⁎, Olga Stephansson c, Petter Mostad b, Mats Gunnar Andersson c a b c

Swedish National Forensic Centre (NFC), SE-581 94 Linköping, Sweden Department of Mathematical Sciences, Chalmers University of Technology and University of Gothenburg, SE-412 96, Gothenburg, Sweden National Veterinary Institute (SVA), SE-751 89, Uppsala, Sweden

a r t i c l e

i n f o

Article history: Received 28 December 2015 Received in revised form 26 August 2016 Accepted 7 October 2016 Available online 11 October 2016 Keywords: Sample analysis Pre-enrichment qPCR Sensitivity Bayesian modelling

a b s t r a c t Efficient and correct evaluation of sampling results with respect to hypotheses about the concentration or distribution of bacteria generally requires knowledge about the performance of the detection method. To assess the sensitivity of the detection method an experiment is usually performed where the target matrix is spiked (i.e. artificially contaminated) with different concentrations of the bacteria, followed by analyses of the samples using the pre-enrichment method and the analytical detection method of interest. For safety reasons or because of economic or time limits it is not always possible to perform exactly such an experiment, with the desired number of samples. In this paper, we show how heterogeneous data from diverse sources may be combined within a single model to obtain not only estimates of detection probabilities, but also, crucially, uncertainty estimates. We indicate how such results can then be used to obtain optimal conclusions about presence of bacteria, and illustrate how strongly the sampling results speak in favour of or against contamination. In our example, we consider the case when B. cereus is used as surrogate for B. anthracis, for safety reasons. The statistical modelling of the detection probabilities and of the growth characteristics of the bacteria types is based on data from four experiments where different matrices of food were spiked with B. anthracis or B. cereus and analysed using plate counts and qPCR. We show how flexible and complex Bayesian models, together with inference tools such as OpenBUGS, can be used to merge information about detection probability curves. Two different modelling approaches, differing in whether the pre-enrichment step and the PCR detection step are modelled separately or together, are applied. The relative importance on the detection curves for various existing data sets are evaluated and illustrated. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The correct interpretation of bacterial test results depends on exactly how tests have been performed (Dahms, 2004; Legan et al., 2001; Morzinski et al., 1996; Van Schothorst et al., 2009; Wehling et al., 2011). Specifically, the testing procedure must be described and preferably modelled statistically. Moreover, if we want to compute the probability of a bacteria being present given some test results, we need to specify a prior probability of presence, i.e., the probability before considering the test data. This is particularly evident if we want to compute the probability of bacterial presence given only negative test results. In fact, in addition to a prior probability of presence, we need to model how

⁎ Corresponding author at: Swedish National Forensic Centre (NFC), SE-581 94, Linköping, Sweden. E-mail addresses: [email protected] (R. Hedell), [email protected] (O. Stephansson), [email protected] (P. Mostad), [email protected] (M.G. Andersson).

http://dx.doi.org/10.1016/j.ijfoodmicro.2016.10.005 0168-1605/© 2016 Elsevier B.V. All rights reserved.

bacterial concentrations varies naturally in the tested material, how concentrations vary in the samples taken (which depends on sampling design), and how concentrations vary within each sample (which depends on sample treatment). Finally, we need to model the specificity of the testing method, and its sensitivity, i.e., the probability of detection at different concentrations. In this paper, we focus on models for sensitivity, and how such models can be obtained from various types of laboratory data. Specifically, we show how to establish models that, given the overall concentration c of a bacterium in a well mixed sample of some matrix, produce a probability distribution for the probability pc that the detection method in question will indeed detect the bacterium. We also consider the variability in the number of colony forming units (CFU) in the sub-samples analysed. Note that Section 2.3 includes an argument why a probability distribution for pc is preferable to a fixed value. Section 2.2.1 presents a basic logistic regression setup for sensitivity models. Section 2.2 then goes on to describe various model components that may extend this setup, depending on the types of available laboratory data.

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

Models for the probability of detection are optimally based on experimental data where the relevant matrix is spiked (i.e., artificially contaminated) with the relevant bacterium at various concentrations (Evers et al., 2010; Koyuncu et al., 2010) and subjected to the relevant detection method. However, for reasons such as cost or safety, one may want to combine this data with information from experiments where similar detection methods are applied to similar bacteria in similar matrices. When similarities are confirmed in the data analysis, the additional information can help reduce uncertainty and bias in the sensitivity model. Statistical methods for achieving such data merging are discussed in Section 2.2.2. Bacterial detection methods often consist of a pre-enrichment step followed by analysis using real time quantitative polymerase chain reaction (qPCR; Butler, 2010) or colony counting on blood agar plates (plate counts). When the available data consists of detections that are similar but not identical in their setup and parameters, one way to take these similarities into account is to use separate and explicit models for both the pre-enrichment phase and the subsequent bacterial detection. Even more fundamentally, data where repeated detections are applied to the same pre-enrichment broth should clearly not be treated as independent observations of detections. Sections 2.2.4–2.2.6 discuss these issues. Our main data example concerns obtaining detection probabilities for Bacillus anthracis, which is the cause of anthrax, a potentially lethal disease for both humans and animals. As the negative consequences of the presence of B. anthracis can be large, it is important to screen objects and areas for it whenever the prior probability of its presence is nonnegligible. Such cases may be related to feed and food quality, tracing or confirmation/falsification of an outbreak or bioterrorism. Incidences in the US 2001, where letters contaminated with anthrax were sent to several people, gave rise to a large scale sampling campaign that later highlighted the need for sound probabilistic methods for evaluating negative sampling results (GAO, 2006). Due to the dangerous nature of B. anthracis the less dangerous bacterium Bacillus cereus is sometimes used as positive control in spiking experiments (Fricker et al., 2011). As capacity of biosafety level 3 and level 4 laboratories is limited it might be tempting to perform spiking experiments using such a surrogate bacterium (in a biosafety level 2 lab) to avoid competing for capacity and to reduce the risk for crosscontamination between casework and experimental samples. In extraordinary events where quick results are crucial or in large investigations the laboratory could be forced to make a trade-off between number of casework and experimental samples analysed and the length of the pre-enrichment time. It is then important that the probability of detection at the shorter pre-enrichment time has been assessed, and that the sensitivity is still deemed acceptable by risk managers in relation to the scope of the investigation and the levels of contamination expected given presence. Quantification of the initial bacterial concentration might also be possible if samples for PCR are taken when the bacterial growth is still in exponential phase (Krämer et al., 2011). We consider models with varying pre-enrichment times, up to 24 h. In Section 2.3, we discuss how models for the probability of detection can be combined with models such as those mentioned in the introductory paragraph to obtain the probability of bacterial presence given test data. It is clear that the assumptions and models used, for example concerning the prior probability of contamination and the assumed variability of concentrations if the bacterium is present, influence results decisively. However, in this paper, possible approaches are only outlined. We analyse a quite heterogeneous data set containing detection results using both qPCR and plate counts, after varying amounts of preenrichment time of either B. anthracis or the surrogate bacterium B. cereus in five different matrices. The data set seems realistic in the sense that it includes data from experiments that where originally performed for other purposes. Section 2.1 presents the details of this data

79

and experimental setup, and Section 3 presents the fitted statistical model for it. Finally, we use the data to exemplify our main idea: By comparing the detection probabilities obtained by using various subsets of the data, we show how the inclusion of heterogeneous data may help estimate and reduce uncertainty around detection probabilities, and how such detection curves may be used to assess how strongly the sampling results speak in favour of or against contamination. 2. Materials and methods 2.1. Experimental data The statistical modelling of the detection probabilities and of the growth characteristics of the bacterium types is based on data from four experiments where different matrices of food were spiked with B. anthracis or B. cereus; see Table 1. A table with more details about the spiking experiments is shown in Appendix B, Table B.1. The matrices are denoted A–E. 2.1.1. Preparation of samples The experiments were performed at different times during 2009– 2014. Samples from the same stored batches were used for matrix B, D and E in all experiments, while different batches were used for data set 1 and 2 for matrix A and C. For each data set a number of spiking replicates were made, in addition to negative controls. Sub-samples were taken from each tube of pre-enrichment broth at different times and analysed using qPCR and for some experiments also plate counts, as described below. The same strain of B. cereus was used in both experiments of data set 2 and 4. Two strains of B. anthracis were used; one in data set 1 and in replicate 1 of data set 3 and another in replicate 2 of data set 3. Initial data analysis didn't suggest any important differences between the two strains of B. anthracis. Therefore these are treated as one B. anthracis type in the statistical modelling. In order to minimize sampling effects 2 g aliquots of the matrix were individually added to the tubes prepared for each experiment and replicate. Subsequently 0.1 ml of the appropriately diluted spore suspension was added directly to the matrix. For data set 1 the spiking levels were different between the matrices with the target spiking levels varying between 0.9 and 155,000 CFU/g matrix. For data set 2 the spiking levels were different between replicates but the same for all matrices within the same spiking replicate. The target spiking levels were between 0.6 and 2300 CFU/g. For all replicates of data set 3 the two target spiking levels were 515,000 and 700,000 CFU/g. In data set 4 the target spiking levels were between 1.4 and 1400 CFU/g. Following absorption (approximately 5 min) 18 ml brain heart infusion broth (BHI) (SVA, art. no. 311001) was added. The total preenrichment volume was about 20 ml (including 2 g of the spiked matrix). Each tube was vortexed (2 min) and incubated, first for 15 min at 70 °C to trigger germination of the Bacillus spores and to reduce the amount of competing background flora followed by pre-enrichment at 37 °C. 2.1.2. Sub-sampling from pre-enrichment broth For the experiment of data set 1 matrices A–D were spiked. In total 5 (spiking levels) × 4 (matrices) × 3 (replicates) = 60 pre-enrichment broths were sampled at 6 and 24 h (plus the negative controls). Each such sample was analysed with qPCR, and one out of three spiking replicates was also analysed with a plate count. In data set 2, B. cereus was spiked at 4 different levels to the matrices A, B, C, D with 3 spiking replicates for A and C, and 2 spiking replicates for B and D. Each of these 40 pre-enrichment broths were sampled at three different times, and each sample was analysed with both qPCR analyses (1–3 repeats) and plate counts (0–2 plate series), with more analyses for matrix A and C and less for B and D. In data set 3, B. anthracis was spiked at 2 spiking levels in matrix E, with 3 spiking replicates. Each of the 6 pre-enrichment

80

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

Table 1 The data sets used. For each data set a number of spiking replicates were made. Sub-samples were taken from each tube of pre-enrichment broth at different times and analysed using qPCR or plate counts. For data set 1 and 2 the experimental setup was not equal in all replicates, with varying matrices, measurement types or points of time when measurements were taken. Data set, index

Bacillus type used

Spiked matrices

Number of spiking levelsa

Spiking replicates

Measurements taken at pre-enrichment times [h]

Measurements made for each spiking level, replicate, matrix and pre-enrichment time

1 2 3 4

anthracis cereus anthracis cereus

A, B, C, D A, B, C, D E E

5 4 2 4

3 2–3 3 1

6, 24 3–4, 5–6, 24 1.5, 3, 4.5, 6 1, 2, 3, 5

1 qPCR, 0–1 plate series 1–3 qPCR, 0–2 plate series 1 qPCR, 1 plate series 1 qPCR, 1 plate series

a

Plus 1–2 negative controls.

broths was sampled at 4 different times, and each such sample was analysed with both qPCR and plate counts. In data set 4, B. cereus was spiked at 4 spiking levels in matrix E, with 1 spiking replicate. Each of the 4 pre-enrichment broths was sampled at 4 different times, and each such sample was analysed with both qPCR and plate counts. 2.1.3. Analysis of samples The sub-samples taken from the pre-enrichment broths were analysed with blood agar plates (SVA, art. no. 341180) using 10-fold dilution series or with qPCR. In all experiments DNA extraction was performed using Biorobot EZ1 (Qiagen) (total volume 200 μl, including internal amplification control (IAC)) and EZ1 DNA tissue kit (Qiagen, art. no. 953034). For data set 1, 2 and 3 a multiplex qPCR analysis was made with Perfecta Multiplex qPCR Supermix with low ROX (VWR) and seal herpes virus type 1 (PhHV1) (Van Doornum et al., 2003) as IAC. In the experiment for data set 4 a SYBR-PCR assay (Applied Biosystems) was used. The qPCR detection results for data set 1 and 2 for 5–6 and 24 h preenrichment are shown in Appendix A, Fig. A.4. 2.2. Models for detection probabilities In the subsections below, we present and discuss a number of modelling fragments useful when putting together detection probability models based on various types of spiking data. 2.2.1. A logit model A single qPCR result yi (0 if bacterium not detected, 1 if detected) can be modelled as Bernoulli distributed (i.e. binomial with one trial): yi
Bernoulliðpi Þ

ð1Þ

A simple possible model for the detection probability pi at concentration ci is a logistic regression model1 logitðpi Þ ¼ a þ b  logðci Þ þ εi

ð2Þ

where
 ε i
N 0; σ 2ε

ð3Þ

is independent for each i. This over-dispersion term is included to reflect the fact that actual detection probabilities are unlikely to to fit exactly on a logit curve. If we have probability distributions reflecting the uncertainty in the parameters a, b and σ2ε after fitting this model using spiking data, we get a model for the uncertainty in pi given various values of ci. Specifically, one may set up a prior distribution for the unknown parameters, and use e.g. Markov Chain Monte Carlo sampling (MCMC) with a tool like OpenBUGS (Lunn et al., 2009) to obtain numerical representations of the joint posterior (Christensen et al., 2011; Gelman et al., 2003), given observations y1 , … , yn at concentrations c1 , … , cn in spiking experiments. The posteriors for a, b and σ2ε can 1

pi logitðpi Þ ¼ logð1−p Þ i

then be used to obtain simulated values for detection probabilities in a case with real data. 2.2.2. Combining data from similar experiments In general, if detection methods s = 1 , … , S are applied to bacteria j = 1 , … , J in matrices k = 1 , … , K one could use independent parameters ajks and bjks in the model of Section 2.2.1. However, it may be quite biologically reasonable to assume that such parameters are related. For example, one may write ajks ¼ μ þ α j þ βk þ γs þ δjk þ ηjs þ πks þ λjks

ð4Þ

with the restrictions that the sum over the parameters in each group is zero, i.e., ∑jαj =∑kβk = … = ∑kλjks = ∑sλjks = 0. Note that if the last four set of parameters are all zero, i.e. δjk = ηjs = πks = λjks = 0, we get a linear model without interaction. If we regard the parameters in Eq. (4) as random variables, we get a model with a limited amount of interaction if the interaction parameters have distributions that constrain them to have much smaller absolute values than those of αj , βk and γs. If we model each of the parameters αj , … , λjks above as normally distributed with expectation zero and use hyperpriors for the variances of these distributions, we get a random effect model:

 α j ∼N 0; σ 2α ; …; λjks ∼N 0; σ 2λ

ð5Þ

Then the posterior variances may be quite different between variables, thus specifying how much interaction there seems to be in the data. Sum-to-zero constrains may be imposed after the MCMC-runs in order to re-center the parameter estimates (Kruschke, 2011). Given a particular data set, a Bayesian data analysis can proceed as follows: Several different models of the type above, where some or all of the parameters of Eq. (4) are included, are compared and a choice is made based on subject area knowledge and the fit to data. For the chosen model, the posterior distribution of the parameters may help transfer information between similar types of spiking experiments. 2.2.3. Uncertainty in spiking concentrations In many spiking experiments, and indeed in some of those of our main example, the initial target concentration is so low that it implies the presence of only a handful CFU's, or possibly no CFU's. For each spiking q = 1 , … , Q the number of CFU's in the pre-enrichment broth at time 0, denoted Mq, can be modelled as Poisson distributed:
 Mq
Poisson uq vs

ð6Þ

where uq is the concentration in the spiking tube containing the bacteria [CFU/ml] and vs is the spiking volume. Further we denote the concentration of CFU in the pre-enrichment broth [CFU/ml] at time 0 as cq0 and the total pre-enrichment volume as vp. Hence, cq0 ¼ Mq =vp

ð7Þ

When this model is combined with detection models like those in Section 2.2.1, the main effect is that the variability in the detection

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

probability increases for very low concentrations. It also allows for the possibility that the test volume by chance contains no CFU's. 2.2.4. Modelling pre-enrichment Bacterial detection methods generally include a pre-enrichment phase, lasting a time t. Let cqt be the concentration after preenrichment. According to Buchanan et al. (1997) a reasonable model connecting cqt with cq0 is a model where the growth of bacteria is characterized by a lag phase with no growth, then a growth phase with exponential growth, and finally a stationary phase where the concentration has reached its maximum. In our notation we write this as cqt ¼ cq0  w f ðt;

τL ; τG Þ

 ϕq

ð8Þ

where w is the growth factor per time unit and f(t,τL, τG) is a function of t, the lag time τL and the end time of the growth phase τG: f ðt; τL ; τ G Þ ¼

8 <

0; t−τ L ; : τG −τL ;

t ≤τL ; lag phase τL bt ≤τG ; growth phase t NτG ; stationary phase

ð9Þ

We also allow for independent over-dispersion factors ϕq for each spiking which may make the model more realistic, fitting the data better. This is an extension of the model of Buchanan et al. (1997). In practice we consider the logarithm of cqt. General models for microbial growth, of different complexity, have previously been discussed by several authors (Nauta, 2002; Peleg and Corradini, 2011; Vimont et al., 2006). Given the concentration cqt after the pre-enrichment process, one may then combine the model above with a detection model like the one in Section 2.2.1 for the actual detection data. An important benefit from using a separate model for the pre-enrichment is that experimental runs that share the same pre-enrichment process can be modelled using the same concentration cqt. We may also use the same postenrichment concentration cqt for experimental runs where this concentration is then tested with several different detection methods. In particular, we may include plate count data, as in the next section. The growth factor w in the model above may be extended to allow for multiple bacterium types, matrix types and interaction effects, similarly as in Section 2.2.2. If the spiking q is made with bacterium j and matrix k the new growth factor is denoted wjk: wjk ¼ h  g j  mk  θjk

ð10Þ

where h is the average growth factor and gj, mk and θjk are random effects related to the bacterium types, matrix types and interaction effects respectively:    


 log g j ∼N 0; σ 2g ; logðmk Þ∼N 0; σ 2m ; log θjk ∼N 0; σ 2θ

ð11Þ

2.2.5. Using plate-count data The model of Section 2.2.1 is best suited for detection methods that only produce a yes/no answer. An alternative, suitable for example when the data is a dilution series plate count, includes a Poisson distribution (Chen et al., 2003; Ishikuri and Hattori, 1985): zqtn



Poisson cqt vr φn

ð12Þ

where zqtn is the total number of colonies for a single plate count series at pre-enrichment time t and vr is the total volume of pre-enrichment broth taken to the plate series. Bacteria in suspension will show some degree of spatial aggregation and the difference in concentration between two sub-samples will be larger than suggested by the Poisson distribution, so we allow for an over-dispersion parameter φn for each

81

unique plate count series n:   logðφn Þ∼N 0; σ 2φ

ð13Þ

Inclusion of plate count data allows for explicit modelling of the bacterial growth. 2.2.6. Approaches to model choice for experimental data We provide two different examples/approaches in which we derive detection probabilities. In one approach the growth of bacteria during the pre-enrichment phase is modelled separately from a model for the probability of detection after pre-enrichment. In the other approach the pre-enrichment phase and detection phase are modelled jointly. The latter approach is denoted below as “approach 1” and the former as “approach 2”. For both approaches we use the spiking model described in Section 2.2.3. 2.2.6.1. Approach 1 - pre-enrichment step and PCR detection step modelled together. In this approach, the models are put together using the elements from Sections 2.2.2 and 2.2.3. For the model of Section 2.2.2, we consider two detection methods; 5–6 h pre-enrichment followed by qPCR and 24 h pre-enrichment followed by qPCR, with the data for 5– 6 h pre-enrichment considered having 6 h pre-enrichment. The models are fitted with the qPCR data of data set 1 and 2 with 5–6 and 24 h pre-enrichment. Hence this approach neither include any plate count data nor data set 3 and 4, as we use these for modelling of bacterial growth, which is not part of this approach. For multiple qPCR samples taken simultaneously from the same pre-enrichment broth only the first measurement is used in this approach, as the results cannot be treated as independent given only the initial concentration c0. 2.2.6.2. Approach 2 - pre-enrichment step and PCR detection step modelled separately. In the second approach, all data of the four datasets are used. Thus, we fit a model combining all the elements from Sections 2.2.1–2.2.5, so the pre-enrichment step and the PCR detection step are modelled separately. We model the difference in detection probability between B. anthracis and B. cereus via their growth characteristics; here as a possible difference in growth speed. Differences in growth speed between different Bacillus types has been observed e.g. by De Siano et al. (2006). Similar effects have been observed for other bacteria, e.g. between different serotypes of Salmonella by Singer et al. (2009). 2.3. Detection probabilities and likelihood ratios Consider the following scenario, in which n independent samples are taken from a batch containing a food matrix and tested for the presence of B. anthracis. The issue is whether the batch is free from contamination. Specifically, consider the hypotheses: H1. The batch contains B. anthracis, and the concentrations of CFU's c1 , … , cn in the n samples will be log-normally distributed with known parameters. H2. The batch does not contain B. anthracis. The parameters of the spatial model under H1 may be given by prior information or risk assessment. We let the first hypothesis refer to the scenario most obviously having legal consequences if being true, i.e. establishing the presence of harmful bacteria, similar to the practice in forensic science where the first hypothesis refers to the prosecutors side (insurance or criminal case) (Aitken and Taroni, 2004). Let each of y1 , … , yn be the sample results, so that yi is 1 or 0 depending on whether sample i has a positive or negative outcome. Assume the total number of detections y= y1 + … +yn is zero. If the risk of false positive results is negligible (so we set P(y = 0| H2) = 1) the likelihood ratio (LR)

82

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

concentrations c1j , … , cnj. Then we may approximate

for the two hypotheses becomes

LR ¼

P ðy ¼ 0jH1 Þ P ðy ¼ 0jH 1 Þ ¼ P ðy ¼ 0jH2 Þ 1

ð14Þ

To obtain the posterior odds of the hypotheses the prior odds is multiplied with the likelihood ratio (Gelman et al., 2003; Kass and Raftery, 1995). To compute the value of the LR above, a specific assumption like the one in H1 about the concentrations ci is needed. The value of the LR is also influenced by uncertainty in the detection probabilities. Intuitively, if there is such uncertainty, one may still be uncertain about the presence of the bacterium even if the expected detection probability is high and a large number of negative tests have been done. More precisely, if knowledge about pc is described with a probability distribution with expectation e and variance v, and if e is, say, 0.8, the probability of making, say, n = 10 independent negative observations will be only (1− 0.8)10 =0.0000001 if there is no variance, but it will grow quickly with increased variance v. As an illustration, assume for mathematical tractability that the uncertainty in pc is modelled with a Beta distribution with expectation e and variance v = e(1− e)/(u + 1) so that pc jc
Betaðue; uð1−eÞÞ

ð15Þ

We get that yjH 1
Beta‐binomialðn; ue; uð1−eÞÞ

ð16Þ

so that P(y = 0 | H1) can be computed mathematically. Fig. 1 shows how P(y = 0 | H1), i.e., the probability of no detections even if there is contamination, can increase dramatically with increased variance in pc. Analytic results like those above are in most cases not available, so in practice one may use simulation: For j = 1 , … , N let c1j , c2j , … , cnj be simulated values for the n sample concentrations (reflecting the uncertainty in these), and let p1j , p2j , … , pnj be corresponding simulated detection probabilities, reflecting the uncertainty in these given the

ð17Þ

2.3.1. Hypothetical sampling scenario The fitted detection probability models for our experimental data are used in a hypothetical sampling scenario. We assume a batch containing food matrix A is investigated for the presence of B. anthracis. From each of 10 randomly drawn samples a sub-sample of 2 g is analysed with a total of y = 0 detections. Let c1 , … , c10 be the number of CFU/g in each sample. We consider the same pair of hypotheses as in the previous section with log10(ci) ~ N(μ, 1) under H1 for all i and for a fixed value μ. Under H1 the number of CFU in each sub-sample is modelled as Poisson distributed; Po(2ck), while under H2 the number of CFU is zero for all sub-samples. Each individual detection result is modelled as an observation from a Bernoulli distribution. Assuming the risk for false positive results is negligible the likelihood ratio is equal to P(y= 0| H1)/1. We evaluate the likelihood ratio under both approach 1 and approach 2, comparing results when the fact finder has access to different subsets of the experimental data of data set 1–4. Either one has access to all the data with B. anthracis (including target matrix A) or to all the data with B. anthracis except for the target matrix A, but instead all the data with B. cereus (including the target matrix A). Hence, we see how the resulting LR changes when a surrogate bacterium is used instead of the target bacterium for the target matrix. 2.4. Software and computation of posterior distributions All models in Section 3.1, fitted to our data, were implemented in OpenBUGS version 3.2.3 (Lunn et al., 2009). In all simulations three MCMC chains with different starting values were used. The simulations were run long enough so that the trace plots did not exhibit any signs of non-convergence. R software version 3.2.2 (R Core Team, 2015) were used for data management, plotting and evaluation of posterior samples and all other computations.

0.15

Probability of n = 10 negative tests

N     1X P y1 ¼ 0jp1 j  …  P yn ¼ 0jpnj N j¼1 N     X 1 ¼ 1−p1 j  …  1−pnj N j¼1

P ðy ¼ 0jH1 Þ ≈

0.10

Below, we exemplify our ideas using two subsets of the data presented in Section 2.1. For each subset, a statistical model is put together from the model elements of Section 2.2. The fitted models for the data subsets are presented here in results Section 3.1 to emphasize how model choice is a process that depends on the actual data. For each of these two subsets of data, or approaches, Section 3.2 compares results for detection probabilities using smaller or larger parts of the subset. This illustrates how inclusion of various types of results from spiking experiments influences corresponding results for detection probabilities. In Section 3.3 use of detection curves is illustrated in the hypothetical sampling scenario.

0.00

0.05

Probability

3. Results

0.00

0.05

0.10

0.15

Variance in detection probability

Fig. 1. The probability of n = 10 negative results can increase quickly with the uncertainty in the detection probability, even if the expected detection probability is 0.8. Above we assume, for mathematical tractability, a Beta distribution for the detection probability pc.

3.1. Fitted model for the experimental data 3.1.1. Approach 1 - pre-enrichment step and PCR detection step modelled together We want the final model to be as simple as possible but still fitting the data well, so we started with a few simple candidate models of

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

83

parameters θjk were modelled as:

Eqs. (2)–(4); the most complex being   logit pijks ¼ μ s þ α j þ βk þ δjk þ bsj  logðc0 Þ þ εi

ð18Þ

Note that we use independent intercept parameters μs for each detection method (qPCR analysis at 6 and 24 h pre-enrichment) rather than using random effect parameters γs, as one can expect that the pre-enrichment times have a strong structural dependence on the detection curves. The different models, each one nested within the model depicted by Eq. (18), were compared via their log pseudo-marginal likelihood (LPML; Christensen et al., 2011), and checked with observed data using posterior predictive plots of the individual detection results (Gelman et al., 2003). The resulting model was   logit pijks ¼ μ s þ α j þ βk þ δjk þ bs  logðc0 Þ

ð19Þ

Notably, the model selection resulted in a model with interaction effects between bacteria type and matrix type δjk included, and with the slope parameter bs depending only on the detection method and with no over-dispersion parameters εi for the logit curves. The posterior predictive check of the final model did not suggest any overall misfit (see Appendix A Fig. A.1). The bacterium specific effects, the matrix specific effects and the interaction parameters were modelled as:
 α 1 ; α 2
N 0; σ 2α

ð20Þ

  β 1 ; …; β4
N 0; σ 2β

ð21Þ


 δ1;1 ; …; δ2;4
N 0; σ 2δ

ð22Þ

The prior distributions used for approach 1 are listed in Table 2. μ6 and μ24 are modelled as weakly informative; as normally distributed with zero mean and standard deviation 100. Realistically the slope parameters b6 , b24 cannot be too extreme and therefore their prior distributions were set as having limited variance with positive mean; N(1, 1). The results turned out to be relatively sensitive to the variance of the hyperprior distributions. Therefore these were set to be not too wide, otherwise unreasonably wide detection curves were obtained. Uniform distributions were applied with upper limits 10, 5 and 1 for σα, σβ and σδ respectively. 3.1.2. Approach 2 - pre-enrichment step and PCR detection step modelled separately For approach 2 we model the growth curve as in Section 2.2.4 but due to slow MCMC convergence we do not include any overdispersion effects ϕq and choose a slightly different expression for the growth factor: wjk ¼ g j  mk  θjk

ð23Þ

where gj were modelled as independent parameters (rather than as random effects). The matrix specific effects mk and the interaction Table 2 The prior distributions used in approach 1. Parameter(s)

Prior distribution

μ6 , μ24 b6 , b24 σα σβ σδ

N(0, 1002) N(1, 1) U(0, 10) U(0, 5) U(0, 1)


 logðm1 Þ; …; logðm5 Þ
N 0; σ 2m

ð24Þ




 log θ1;1 ; …; log θ2;5
N 0; σ 2θ

ð25Þ

As in Section 2.2.1 the individual qPCR results were modelled as Bernoulli variables, with the detection probabilities being dependent on the concentration cjkt of the target bacterium at time t:  
 logit pijkt ¼ a þ b  log cjkt þ εi

ð26Þ

Posterior predictive checks of the qPCR data indicated lack of fit between the model and the data if no over-dispersion effects were included (results not shown). As no obvious pattern were observed we choose to include independent over-dispersion parameters εi for each unique qPCR analysis i (where “qPCR analysis” means sub-sampling from the pre-enrichment broth, DNA extraction and amplification with qPCR):
 εi
N 0; σ 2ε

ð27Þ

Further, the plate counts were modelled as in Section 2.2.5. Again, posterior predictive checks indicated the need for over-dispersion effects φn for each unique plate count series n, modelled as independent parameters:   logðφn Þ∼N 0; σ 2φ

ð28Þ

The prior distributions used in approach 2 are listed in Table 3. The distributions were set as weakly informative; as uniformly distributed from 0 to 5 h for τL and 8 to 20 h for τG, independent half-normal distributions (normal distributions folded at zero taking only positive values) for all log(gj) with standard deviation 10 and N(0, 1002) for a. For the same reasons as in approach 1 the slope parameter b and the hyperprior distributions were set as having restricted variances; N(1, 1) for b and uniform distributions for the hyperprior parameters σm, σθ, σε and σφ with upper limits 2, 0.5, 2, and 5 respectively. Posterior predictive checks with the final models did not suggest any overall misfit (see Appendix A Fig. A.1). Under both approaches, 1 and 2, the posterior distribution of the spiking levels indicated that the number of spiked CFU were likely to be zero in a few cases (results not shown), explaining some of the negative sampling results. 3.2. Comparison of detection probability results based on subsets of experimental data Results for approach 1 and 2 fitted with different subsets of the experimental data are shown in Figs. 2–4 and A.2–A.4 and Table 4. In Fig. 2 curves for the probability of detection relating to approach 1 are shown. Note that the plots include an uncertainty interval for each spiking level. The solid line is based on all data (data set 1–2). The dotted Table 3 The prior distributions used in approach 2. Parameter(s)

Prior distribution

τL τG log(gj) a b σm σθ σε σφ

U(0, 5) U(8, 20) Half-N(0, 102) N(0, 1002) N(1, 1) U(0, 2) U(0, 0.5) U(0, 2) U(0, 5)

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

B,6h

1

2

3

4

5

0.8 P(detect)

0.0

0.0 0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

log10(cfu/g)

log10(cfu/g)

A , 24 h

B , 24 h

C , 24 h

D , 24 h

2

3

4

5

0

log10(cfu/g)

1

2

3

4

5

5

4

5

0.0

0.4

P(detect)

0.8 0.4

P(detect)

0.0

0.0

0.4

P(detect)

0.8 0.4

1

4

0.8

log10(cfu/g)

0.8

log10(cfu/g)

0.0 0

D,6h

0.4

P(detect)

0.8 0.0

0.4

P(detect)

0.8 0.4

P(detect)

0.0 0

P(detect)

C,6h 0.8

A,6h

0.4

84

0

1

log10(cfu/g)

2

3

4

5

0

1

log10(cfu/g)

2

3

log10(cfu/g)

Fig. 2. Probability of detection (50% intervals; interquartile range) for different spiking levels for 6 and 24 h pre-enrichment for matrix A–D. The target bacterium is B. anthracis. Approach 1 used (pre-enrichment step and PCR detection step modelled together). Solid line: Results based on all data (data set 1–2). Dotted line: Results based on all data (data set 1–2) but one matrix for B. anthracis (A–D respectively).

line is based on all data (data set 1–2) but one matrix for B. anthracis (A–D respectively). A similar figure for approach 2 based on data set 1–4 is shown in Fig. 3. In general, the intervals depicted by dotted lines (matrix A–D respectively for B. anthracis removed) are relatively similar to the intervals depicted by the solid lines (fitted to all data). I.e. extrapolating the information from all the experiments with the surrogate bacterium B. cereus together with the data for B. anthracis for all but one matrix gives quite similar results as using the full data. The results for approach 1 and 2 fitted to all B. anthracis data except for one matrix (A–D respectively) are shown in Figs. A.2 and A.3 respectively in Appendix A (solid lines). Using only the B. anthracis data to model the curve for the left out matrix for approach 1 (Fig. A.2, solid lines) gave somewhat wider intervals than when including also all the B. cereus data (Fig. A.2, dotted lines). For two of the matrices (A and C) the position of the curves is also somewhat different. I.e. by including the B. cereus data to the B. anthracis data both the position and the width of the curves are affected. With approach 2 the differences in interval widths are less but the position of the curves are still affected (Fig. A.3). Fig. 4 further exemplifies how the inclusion of heterogeneous data may help estimate and reduce uncertainty around detection

C,6h

3

4

5

0

1

2

3

4

5

0.4

0.8 1

2

3

4

5

0

1

2

3

log10(cfu/g)

A , 24 h

B , 24 h

C , 24 h

D , 24 h

3

4

5

0

1

2

3

log10(cfu/g)

4

5

5

4

5

0.0

0.4

P(detect)

0.8 0.4

P(detect)

0.0

0.0

0.4

P(detect)

0.8

2

log10(cfu/g)

4

0.8

log10(cfu/g)

0.4

1

P(detect) 0

log10(cfu/g)

0.0 0

0.0

0.0

0.0 2

log10(cfu/g)

0.8

1

D,6h

0.4

P(detect)

0.8 0.4

P(detect)

0.8 0.4

P(detect)

For the hypothetical sampling scenario the mean parameter μ is given value − 2, 0 or 2, approximately referring to low, medium and high expected detection probability for matrix A (see Fig. 2–3). The results are shown in Table 4. For this example, using the surrogate bacterium B. cereus instead of the target bacterium B. anthracis for the target matrix does not result in any substantial differences of the likelihood ratios for either approach. In all cases having no positive samples makes

B,6h

0.0 0

P(detect)

3.3. Hypothetical sampling scenario

0.8

A,6h

probabilities. The result is based on approach 1. The solid lines, based only on the B. cereus data for matrix C, predict the probability of detection for B. anthracis for the same matrix. The dotted lines on the other hand are based on all data (data set 1–2) but matrix C for B. anthracis. Hence, the inclusion of the extra datasets combined with our hierarchical model greatly affects the posterior distribution of detection probabilities. There were notable differences in detection curves between many of the matrices (see e.g. Fig. 2) and between bacteria types, with B. cereus having a higher probability of detection (Fig. A.4).

0

1

2

3

log10(cfu/g)

4

5

0

1

2

3

log10(cfu/g)

Fig. 3. Probability of detection (50% intervals; interquartile range) for different spiking levels for 6 and 24 h pre-enrichment for matrix A–D. The target bacterium is B. anthracis. Approach 2 used (pre-enrichment step and PCR detection step modelled separately). Solid line: Results based on all data (data set 1–4). Dotted line: Results based on all data (data set 1–4) but one matrix for B. anthracis (A–D respectively).

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

0.6 0.4 0.0

0.2

P(detect)

0.8

1.0

C,6h

0

1

2 3 log10(cfu/g)

4

5

Fig. 4. Probability of detection (95% intervals) for different spiking levels at 6 h preenrichment for matrix C. The target bacterium is B. anthracis. Approach 1 used (preenrichment step and PCR detection step modelled together). Solid line: Results based only on the B. cereus data for matrix C (data set 2). Dotted line: Results based on all data (data set 1–2) but matrix C for B. anthracis.

the likelihood for H1 higher than H2. For the highest concentration under H1 (i.e. log10(ci) ~ N(2, 1)) all results suggests strong evidence against contamination (LR b 0.001). For the lowest concentration under H 1 (log10 (ci ) ~ N(− 2, 1)) the evidence is more or less inconclusive about the contamination status (LR close to 1). With the concentrations log 10 (c i ) distributed according to N(0, 1), the evidence suggest substantial evidence against contamination only with 24 h pre-enrichment (LR close to 0.01). 4. Discussion Our data and statistical analysis indicated differences in detection probabilities between B. anthracis, B. cereus and also between several matrices (Fig. A.4). Thus, it was important to allow for possible differences between these in the statistical model when merging the available data. In general, when extrapolating information from a surrogate experiment to the case of interest it is important that the overall amount of data or expert knowledge is not too scarce. Otherwise very large uncertainties for detection curves may be obtained. Statistical decision theory may help in deciding whether a new spiking experiment with the bacterium and matrix of main interest should be performed instead. For evaluation of sampling results, the detection models must be combined with models for the bacterial concentration in the batch. For the sake of demonstration we have, in the working examples, assumed that the spore concentration in individual subsamples follows a Poisson log-normal distribution with standard deviation equal to 1 which may be considered a reasonable model for an organism unevenly distributed in a commodity (Van Schothorst et al., 2009). In actual casework it might be relevant to allow for other spatial structures (Jongenburger et al., 2012). Prior information, risk assessment, data

85

from other outbreaks or experimental data might guide what distributions are plausible and practically relevant assuming that the bacteria are present. Modelling of within-batch and between-batch variability was also discussed by Gonzales-Barron and Butler (2011). Further, models for the sampling design are also needed for the interpretation of evidence. We made no other assumption than simple random sampling of samples and sub-samples. In some cases other sampling strategies may be used, such as systematic sampling (Jongenburger et al., 2011). If the sample treatment differs from our setup this should also be accounted for in the models. For instance we assumed the samples taken are well mixed. For the specific data analysed in this paper, our models showed that the interaction between variation in the matrix and variation between the bacteria was non-ignorable, and thus that data from B. cereus experiments were somewhat less relevant for B. anthracis detection curves as one might have hoped. With our choice of priors the non-decline in uncertainty intervals was observed e.g. in Fig. A.3. As the experiments were performed at different occasions with several years in between it is possible that some of the interaction effects are due to experimentwise random effects. It would be interesting to study the differences between matrices and bacteria under even more similar conditions. Despite this, by using data from several matrices as well as data for B. cereus it was possible to mimic the detection curves for B. anthracis reasonably well, as shown in Figs. 2–3 and Table 4. There are some potential advantages with dividing the statistical modelling into a pre-enrichment step and a detection step. If the preenrichment protocol or the analytical detection protocol is changed in some way then parts of the model can be kept and do not have to be re-evaluated; e.g. if the PCR protocol is changed the posterior distributions of the pre-enrichment phase parameters could be kept as prior distributions for the new setup. If the knowledge in one part or the other of the model is good, a smaller experiment might be performed to assess the complete model rather than performing a more exhaustive experiment. It is also possible to use data from more types of experiments in order to further reduce the uncertainties in the detection probabilities, e.g. using data from quality assurance experiments specifically made to study the PCR performance (without pre-enrichment). A downside of using the divided model of approach 2 is that the computational burden increases a lot as more parameters have to be estimated, as well as consideration of more prior distributions. In our case the MCMC convergence was very much faster with approach 1 (pre-enrichment step and PCR detection step modelled together) than approach 2 (pre-enrichment step and PCR detection step modelled separately); minutes compared to several hours on a standard PC laptop and with less autocorrelations. Although we have not done so here, there is also some potential to utilize the cycle threshold (CT) value of the qPCR analysis, rather than just the detect/non-detect result. The relationship between the CT value and the original spiking level including 8 h pre-enrichment was studied by Krämer et al. (2011). Knutsson et al. (2002) considered a model relating the CT-values of the qPCR to the concentration of DNA of Salmonella. We implemented these models for B. anthracis and B. cereus, but we consider the modelling framework to be generally applicable to other bacteria types as well. With other bacteria it might be more relevant to consider also the difference between naturally and artificially contaminated samples.

Table 4 Likelihood ratios with approaches 1 and 2 fitted with different subsets of data and under different spatial models under hypothesis H1. The target matrix is A. Distribution of log10(ci) under H1

Pre-enrichment time

LR with approach 1 and all B. anthracis data

LR with approach 1 and all B. cereus data and all B. anthracis data except matrix A

LR with approach 2 and all B. anthracis data

LR with approach 2 and all B. cereus data and all B. anthracis data except matrix A

N(−2, 1) N(0, 1) N(2, 1) N(−2, 1) N(0, 1) N(2, 1)

6 6 6 24 24 24

0.98 0.38 b0.0002 0.78 0.01 b0.0002

0.88 0.23 0.0008 0.61 0.0072 b0.0002

0.99 0.54 0.0007 0.54 0.0013 b0.0002

0.90 0.19 0.0007 0.59 0.021 b0.0002

86

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

With the spore-form of Bacillus this difference is presumably less of a problem since the resting, metabolically inactive, spores found in nature and laboratory, can be expected to be quite similar with respect to germination rate. In contrast, the conditions in the commodity might have a significant impact on the metabolism of vegetative bacteria, something which may have to be accounted for in detection models. As illustrated in Section 3.3, a likelihood ratio, measuring the relative support for one hypothesis compared to the other, may be computed for the sampling results for each pair of hypotheses. Decision makers or risk assessors may then combine these with their own prior probabilities in order to assess the posterior odds' or posterior probabilities of the different hypotheses. Hence, our approach can be used to support informed decision making regarding possibly contaminated foods and environments, analogous to practice in forensic sciences (Aitken and Taroni, 2004). We have shown how flexible and complex statistical models can be used to merge information about detection probability curves from several diverse data sets. In many cases, experimental designs may deviate from an

ideal design due to practical constraints. Still, if the statistical models are flexible enough they can include such data for the evaluation of evidence.

Acknowledgements The Swedish Laboratory for Food Safety and Preparedness (RUB) is acknowledged for laboratory work. Thanks to Anders Nordgaard for comments on the manuscript and for discussions. RH was partly financed by the Swedish Civil Contingencies Agency (MSB), project: MSB-SäkProv. Part of GA's work was supported by/executed in the framework of the EU project AniBioThreat (Grant Agreement: Home/2009/ISEC/AG/191) with financial support from the Prevention of and Fight against Crime Programme of the European Union, European Commission-Directorate General Home Affairs. This publication reflects the views only of the authors, and the European Commission cannot be held responsible for any use that may be made of the information contained therein.

Appendix A Figures for posterior predictive checks are shown in Fig. A.1. In Figs. A.2 and A.3 detections curves for approach 1 and 2 respectively, fitted with different parts of data, are presented. In Fig. A.4 are results for approach 1 presented, comparing results for B. anthracis and B. cereus.

0.8 0.6 0.0

0.2

0.4

P(ypred = y)

0.0

0.2

0.4

0.6

P(ypred = y)

0.8

1.0

Approach 2, post. pred. check

1.0

Approach 1, post. pred. check

0

50

100

150

200

0

50

100

qPCR analysis

150

200

250

300

qPCR analysis

Fig. A.1. Probability that the posterior predictive qPCR detection result (ypred) equals the observed qPCR detection result (y) for each qPCR analysis, sorted from lowest to highest agreement. Left: Approach 1 used (pre-enrichment step and PCR detection step modelled together). Results based on all data (data set 1–2). 95% of the qPCR analyses have above 50% agreement between predictions and observed data. Right: Approach 2 used (pre-enrichment step and PCR detection step modelled separately). Results based on all data (data set 1–4). 94.5% of the qPCR analyses have above 50% agreement between predictions and observed data.

B,6h

2

3

4

5

1

2

3

4

5

0.4

0.8 1

2

3

4

5

0

1

2

3

log10(cfu/g)

A , 24 h

B , 24 h

C , 24 h

D , 24 h

3

4

5

0

1

2

3

log10(cfu/g)

4

5

5

4

5

0.0

0.4

P(detect)

0.8 0.4

P(detect)

0.0

0.0

0.4

P(detect)

0.8

2

log10(cfu/g)

4

0.8

log10(cfu/g)

0.4

1

P(detect) 0

log10(cfu/g)

0.0 0

0.0

0.0 0

log10(cfu/g)

0.8

1

D,6h

0.4

P(detect)

0.8 0.0

0.4

P(detect)

0.8 0.4

P(detect)

0.0 0

P(detect)

C,6h 0.8

A,6h

0

1

2

3

log10(cfu/g)

4

5

0

1

2

3

log10(cfu/g)

Fig. A.2. Probability of detection (95% intervals) for different spiking levels for 6 and 24 h pre-enrichment for matrix A–D. The target bacterium is B. anthracis. Approach 1 used (pre-enrichment step and PCR detection step modelled together). Solid line: Results based on all B. anthracis data (data set 1) but one matrix (A–D respectively). Dotted line: Results based on all data (data set 1–2) but one matrix for B. anthracis (A–D respectively).

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

B,6h

3

4

5

0

1

2

3

4

5

0.4

0.8 1

2

3

4

5

0

1

2

3

log10(cfu/g)

A , 24 h

B , 24 h

C , 24 h

D , 24 h

3

4

5

0

1

2

3

4

5

5

4

5

0.0

0.4

P(detect)

0.8 0.4

P(detect)

0.0

0.0

0.4

P(detect)

0.8

2

log10(cfu/g)

4

0.8

log10(cfu/g)

0.4

1

P(detect) 0

log10(cfu/g)

0.0 0

0.0

0.0

0.0 2

log10(cfu/g)

0.8

1

D,6h

0.4

P(detect)

0.8 0.4

P(detect)

0.8 0.4

P(detect)

0.0 0

P(detect)

C,6h 0.8

A,6h

87

0

log10(cfu/g)

1

2

3

4

5

0

1

log10(cfu/g)

2

3

log10(cfu/g)

Fig. A.3. Probability of detection (95% intervals) for different spiking levels for 6 and 24 h pre-enrichment for matrix A–D. The target bacterium is B. anthracis. Approach 2 used (pre-enrichment step and PCR detection step modelled separately). Solid line: Results based on all B. anthracis data (data set 1, 3) but one matrix (A–D respectively). Dotted line: Results based on all data (data set 1–4) but one matrix for B. anthracis (A–D respectively).

0.6 0.4 0.0

0.2

P(detect)

0.8

1.0

C,6h

0

1

2 3 log10(cfu/g)

4

5

Fig. A.4. Probability of detection (50% intervals; interquartile range) for different spiking levels for 6 and 24 h pre-enrichment for matrix A–D. Approach 1 used (pre-enrichment step and PCR detection step modelled together). Results based on all data (data set 1–2). Solid line: results for B. anthracis. Dotted line: Results for B. cereus. Filled circles: qPCR detection results (data) for B. anthracis for the target spiking levels. Open circles: qPCR detection results (data) for B. cereus for the target spiking levels.

Appendix B Detailed descriptions of the spiking experiments are shown in Table B.1. Table B.1 The data sets used. For each data set a number of spiking replicates were made. Sub-samples were taken from each tube of pre-enrichment broth at different times and analysed using qPCR or plate counts. For data set 1 and 2 the experimental setup was not equal in all replicates, with varying matrices, measurement types or points of time when measurements were taken. An empty cell means that the value is the same as in the cell above. Data set, index

Bacillus type used

Spiked matrix

Spiking replicate, index

Spiking levelsa, [CFU/g]

Measurements taken at pre-enrichment times [h]

Measurements made for each spiking level, replicate, matrix and pre-enrichment time

1

anthracis

A

1

1.3, 13, 130, 1300, 13,000

6, 24

1 qPCR, 1 plate series

B

2, 3 1

6, 24 6, 24

1 qPCR 1 qPCR, 1 plate series

C

2, 3 1

6, 24 6, 24

1 qPCR 1 qPCR, 1 plate series

D

2, 3 1

6, 24 6, 24

1 qPCR 1 qPCR, 1 plate series

1.8, 18, 180, 1800, 18,000 15.5, 155, 1550, 15,500, 155,000 0.9, 9, 90, 900, 9000

(continued on next page)

88

R. Hedell et al. / International Journal of Food Microbiology 241 (2017) 78–88

Table B.1 (continued) Data set, index

Bacillus type used

Spiked matrix

Spiking replicate, index

Spiking levelsa, [CFU/g]

2

cereus

A

2, 3 1

2.3, 23, 230, 2300

2

0.6, 6, 60, 600

3

2.2, 22, 220, 2200

1, 2, 3 1

515,000, 700,000 1.4, 14, 140, 1400

B C D A B C A C

3 4 a

anthracis cereus

D E E

Measurements taken at pre-enrichment times [h]

Measurements made for each spiking level, replicate, matrix and pre-enrichment time

6, 24 3, 24 5 3, 5, 24 4, 24 6 4 6, 24 3 5, 24 4 24 4, 6, 24 4, 6, 24 3 5, 24 3, 5, 24 1.5, 3, 4.5, 6 1, 2, 3, 5

1 qPCR 1 qPCR 1 qPCR, 1 plate series 1 qPCR 1 qPCR 1 qPCR, 1 plate series 1 qPCR 1 qPCR, 1 plate series 1 qPCR 1 qPCR, 2 plate series 1 qPCR 1 qPCR, 1 plate series 2 qPCR 3 qPCR 1 qPCR 1 qPCR, 2 plate series 1 qPCR 1 qPCR, 1 plate series 1 qPCR, 1 plate series

Plus 1–2 negative controls.

References Aitken, C., Taroni, F., 2004. Statistics and the Evaluation of Evidence for Forensic Scientists. John Wiley & Sons, Chichester. Buchanan, R.L., Whiting, R.C., Damert, W.C., 1997. When is simple good enough: a comparison of the Gompertz, Baranyi and tree-phase linear models for fitting bacterial growth curves. Food Microbiol. 14, 313–326. Butler, J.M., 2010. Fundamentals of Forensic DNA Typing. Elsevier Academic Press, San Diego. Chen, C.Y., Nace, G.W., Irwin, P.L., 2003. A 6 × 6 drop plate method for simultaneous colony counting and MPN enumeration of Campylobacter jejuni, Listeria monocytogenes, and Escherichia coli. J. Microbiol. Methods 55, 475–479. Christensen, R., Johnson, W., Branscum, A., Hanson, T.E., 2011. Bayesian Ideas and Data Analysis. CRC Press, Boca Raton. Dahms, S., 2004. Microbiological sampling plans - statistical aspects. Mitt. Lebensmittelunters. Hyg. 95, 32–44. De Siano, T., Padhi, S., Schaffner, D.W., Montville, T.J., 2006. Growth characteristics of virulent Bacillus anthracis and potential surrogate strains. J. Food Prot. 69, 1720–1723. Evers, E.G., Post, J., Putirulan, F.F., Van der Wal, F.J., 2010. Detection probability of Campylobacter. Food Control 21, 247–252. Fricker, M., Ågren, J., Segerman, B., Knutsson, R., Ehling-Schulz, M., 2011. Evaluation of Bacillus strains as model systems for the work on Bacillus anthracis spores. Int. J. Food Microbiol. 145, 129–136. GAO, 2006. United States Government Accountability Office, anthrax detection: agencies need to validate sampling activities in order to increase confidence in negative results. GAO-05-251. Gelman, A., Carlin, J.B., Stern, H.S., Rubin, D.B., 2003. Bayesian Data Analysis. second edition. Chapman and Hall, London. Gonzales-Barron, U., Butler, F., 2011. Characterisation of within-batch and between-batch variability in microbial counts in foods using Poisson-gamma and Poisson-lognormal regression models. Food Control 22, 1268–1278. Ishikuri, S., Hattori, T., 1985. Formation of bacterial colonies in successive time intervals. Appl. Environ. Microbiol. 49, 870–873. Jongenburger, I., Reij, M.W., Boer, E.P.J., Gorris, L.G.M., Zwietering, M.H., 2011. Random or systematic sampling to detect a localised microbial contamination within a batch of food. Food Control 22, 1448–1455. Jongenburger, I., Bassett, J., Jackson, T., Zwietering, M.H., Jewell, K., 2012. Impact of microbial distributions on food safety I. Factors influencing microbial distributions and modelling aspects. Food Control 26, 601–609. Kass, R.E., Raftery, A.E., 1995. Bayes factors. J. Am. Stat. Assoc. 90, 773–795.

Knutsson, R., Löfström, C., Grage, H., Hoorfar, J., Rådström, P., 2002. Modeling of 5′ nuclease real-time responses for optimization of a high-throughput enrichment PCR procedure for Salmonella enterica. J. Clin. Microbiol. 40, 52–60. Koyuncu, S., Andersson, M.G., Häggblom, P., 2010. Accuracy and sensitivity of commercial PCR-based methods for detection of Salmonella enterica in feed. Appl. Environ. Microbiol. 76, 2815–2822. Krämer, N., Löfström, C., Vigre, H., Hoorfar, J., Bunge, C., Malorny, B., 2011. A novel strategy to obtain quantitative data for modelling: combined enrichment and real-time PCR for enumeration of salmonellae from pig carcasses. Int. J. Food Microbiol. 145, 86–95. Kruschke, J.K., 2011. Doing Bayesian Data Analysis – a Tutorial with R and BUGS. Academic Press/Elsevier, Burlington, MA. Legan, J.D., Vandeven, M.H., Dahms, S., Cole, M.B., 2001. Determining the concentration of microorganisms controlled by attributes sampling plans. Food Control 12, 137–147. Lunn, D., Spiegelhalter, D., Thomas, A., Best, N., 2009. The BUGS project: evolution, critique, and future directions. Stat. Med. 28, 3049–3067. Morzinski, J.A., Jackson, P.J., Picard, R.R., Yoshida, T.M., 1996. Procedures for Collection and Handling of Forensic Samples Suspected of Containing Microbial Pathogens or Toxins. Los Alamos National Laboratory. Nauta, M.J., 2002. Modelling bacterial growth in quantitative microbiological risk assessment - is it possible? Int. J. Food Microbiol. 73, 297–304. Peleg, M., Corradini, M.G., 2011. Microbial growth curves: what the models tell us and what they cannot. Crit. Rev. Food Sci. Nutr. 51, 917–945. R Core Team, 2015. A language and environment for statistical computing. R Foundation for Statistical Computing (Vienna, Austria, http://www.R-project.org). Singer, R.S., Mayer, A.E., Hanson, T.E., Isaacson, R.E., 2009. Do microbial interactions and cultivation media decrease the accuracy of Salmonella surveillance systems and outbreak investigations. J. Food Prot. 72, 707–713. Van Doornum, G.J.J., Guldemeester, J., Osterhaus, A.D.M.E., Niesters, H.G.M., 2003. Diagnosing herpesvirus infections by real-time amplification and rapid culture. J. Clin. Microbiol. 41, 576–580. Van Schothorst, M., Zwietering, M.H., Ross, T., Buchanan, R.L., Cole, M.B., 2009. Relating microbiological criteria to food safety objectives and performance objectives. Food Control 20, 967–979. Vimont, A., Vernozy-Rozand, C., Montet, M.P., Lazizzera, C., Bavai, C., Delignette-Muller, M.-L., 2006. Modeling and predicting the simultaneous growth of Escherichia coli O157:H7 and ground beef background microflora for various enrichment protocols. Appl. Environ. Microbiol. 72, 261–268. Wehling, P., LaBudde, R.A., Brunelle, S.L., Nelson, M.T., 2011. Probability of detection (POD) as a statistical model for the validation of qualitative methods. J. AOAC Int. 94, 335–347.