Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates

Food Microbiology xxx (2014) 1e9

Contents lists available at ScienceDirect

Food Microbiology journal homepage: www.elsevier.com/locate/fm

Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates
gis Pouillot*, Yuhuan Chen, Karin Hoelzer Re U.S. Food and Drug Administration / Center for Food Safety and Applied Nutrition (HFS-005), 5100 Paint Branch Parkway, College Park, MD 20740, USA

a b s t r a c t Keywords: Quantitative microbial risk assessment Bacterial concentration Bacterial number Modeling process pathway Simulation study

When developing quantitative risk assessment models, a fundamental consideration for risk assessors is to decide whether to evaluate changes in bacterial levels in terms of concentrations or in terms of bacterial numbers. Although modeling bacteria in terms of integer numbers may be regarded as a more intuitive and rigorous choice, modeling bacterial concentrations is more popular as it is generally less mathematically complex. We tested three different modeling approaches in a simulation study. The first approach considered bacterial concentrations; the second considered the number of bacteria in contaminated units, and the third considered the expected number of bacteria in contaminated units. Simulation results indicate that modeling concentrations tends to overestimate risk compared to modeling the number of bacteria. A sensitivity analysis using a regression tree suggests that processes which include drastic scenarios consisting of combinations of large bacterial inactivation followed by large bacterial growth frequently lead to a >10-fold overestimation of the average risk when modeling concentrations as opposed to bacterial numbers. Alternatively, the approach of modeling the expected number of bacteria in positive units generates results similar to the second method and is easier to use, thus potentially representing a promising compromise. Published by Elsevier Ltd.

1. Introduction Quantitative microbial risk assessment (QMRA) models typically evaluate the dynamics of bacterial1 populations in food along a food-production pathway to determine contamination at consumption, which is integrated with a dose-response curve to estimate risk to the consumer. Usually, the production pathway is described by using a number of discrete steps along the “farm-tofork” continuum, and these steps may be grouped into various distinct modules (Cassin et al., 1998). At each process step, bacteria may be introduced, can increase or decrease, or may be completely eliminated from certain food units. These steps can be modeled as a set of well-characterized basic processes that may impact bacterial prevalence and/or levels. Various different basic processes have been defined in the literature including, for example, growth, inactivation, food partitioning, mixing, removal, and cross-

* Corresponding author. Tel.: þ1 240 402 2123. E-mail address: [email protected] (R. Pouillot). 1 We used the term “bacteria” throughout this manuscript, but the study is generalizable to other microorganisms of interest.

contamination (Nauta, 2008). FDA-iRISK®, a web-based quantitative risk assessment tool, defines a set of slightly different basic processes: increase by growth, increase by addition, decrease, redistribution, food pooling, food partitioning and food evaporation (Chen et al., 2013). Because of the diversity of food production, processing, and handling situations that may occur, there is no standardized way to model these processes, even though some general modeling frameworks have been developed (FDA-iRISK, 2012; Nauta, 2008, 2005). When modeling these processes, a fundamental consideration for risk assessors is to decide whether to model changes in bacterial levels in terms of concentrations, defined on a continuous scale, usually expressed in cfu/g or in log cfu/g, or in terms of actual bacterial numbers, on a discrete scale, usually expressed in cfu per food unit. A generic distribution of concentrations was for example used as the starting step to characterize prevalence and levels of Listeria monocytogenes in 23 ready-to-eat foods in a QMRA conducted jointly by the U.S. Food and Drug Administration and the U.S. Department of Agriculture Food Safety and Inspection Service (FDA/FSIS, 2003). In some published QMRAs risk assessors chose to model changes in bacterial concentration, such as Escherichia coli O157:H7 concentrations in spinach (Danyluk and Schaffner, 2011)

http://dx.doi.org/10.1016/j.fm.2014.05.008 0740-0020/Published by Elsevier Ltd.

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

2

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

Table 1 Characteristics of the three evaluation methods. Characteristics Prevalence (P) Contamination

Definition Evaluation Definition

Unit System Evaluation Parameters (e.g. log-reduction, growth, redistribution factor, size of unit after partition, size of units after mixing) Integration of the model

Method #1

Method #2

Method #3

Prob (C > 0 cfu/g) Deterministic Concentration (C)

Prob (N  1 cfu/food unit) Deterministic Number per contaminated unit (N)

log10 cfu/g Rational Deterministic Variable

cfu in the food unit Natural  1 Stochastic Variable

Prob (N  1 cfu/food unit) Deterministic Expected number per contaminated unit (E[N]) cfu in the food unit Rational  1 Deterministic Variable

Monte-Carlo

Monte-Carlo

Monte-Carlo

and Salmonella concentrations in almonds (Lambertini et al., 2012). In other QMRAs, risk assessors have modeled the number of bacteria per food unit, such as Salmonella Enteritidis cells in shell eggs (Schroeder et al., 2006), while in yet other QMRAs both were considered, depending on the basic process being modeled (Rigaux et al., 2014). While modeling integer numbers of bacteria may be a more intuitive choice, modeling bacterial concentrations is generally less mathematically complex. Some recent efforts have strived to describe the way in which contamination is handled mathematically (FDA-iRISK, 2012; Nauta, 2008). The objective of this study was to evaluate the impact of specific modeling choices on risk estimates and model complexity in a QMRA. Specifically, we evaluated the impact of considering bacterial populations as concentration (per unit of volume or weight of food, e.g. cfu/g or log cfu/g) or as absolute numbers of bacteria per food unit (e.g. per batch, package, can or leaf, e.g. cfu/can). We report the results of a specific comparison among three different approaches to model changes in bacterial prevalence and levels in contaminated food products. One method tracks the bacterial concentrations; a second method models changes in the number of bacteria in contaminated food units, and a third, novel method, evaluates the expected number of bacteria in contaminated food units. To evaluate the potential impact of the choice in modeling approach on risk estimates, a simulation study was performed, evaluating randomly assembled, but plausible, sets of food pathways which included five distinct process steps. 2. Methods 2.1. Description of the modeling approaches We derived three distinct modeling approaches, designated henceforth as method #1, #2 and #3. 2.1.1. Method #1 This approach tracks bacterial concentrations in food units. At step j of the modeled food-production pathway, contamination of the food units is characterized by: i) prevalence Pj, defined as the proportion of food units with a mean bacterial concentration >0 cfu/g, and ii) mean bacterial level Cj in the contaminated food units, modeled as some rational number (e.g. 1.23, or 2.5; unit: log10 cfu/g). As an illustrative example, consider food units contaminated at a given step of the pathway with Pj ¼ 10% and Cj ¼ 1.22 log10 cfu/g (i.e., 6 cfu in 100 g). In this case, 90% of the food products have a concentration of exactly 0 cfu/g and the mean concentration in the other 10% units is 1.22 log10 cfu/g. Note that at a relatively low mean concentration, the actual number of bacteria in any given quantity of product may be 0 cfu even though the mean concentration is > 0 cfu/g. In this approach, changes in prevalence and concentration are evaluated deterministically for

any of the basic processes, given a specified parameter (see next section and Table 1). As an illustrative example, consider bacterial growth of 3.2 log10 cfu in the food units. This would lead to no change in prevalence and an increase in the bacterial level in contaminated products Cjþ1 ¼ 1.22 þ 3.2 ¼ 1.98 log10 cfu/g (corresponding to 9509.4 cfu in a 100 g food unit). 2.1.2. Method #2 This approach tracks prevalence and number of bacteria in contaminated food units. The prevalence is the proportion of food units that contains 1 cfu. As an illustrative example, consider food units contaminated with prevalence Pj ¼ 10% and number of bacteria per contaminated unit nj  1, where nj is a natural number. In this case, 90% of the food units contain exactly 0 cfu while the remaining 10% of food units are contaminated with some natural number of bacteria (e.g., 1, 34, or 514; unit: cfu/food unit). For any given process with a given parameter, the prevalence is evaluated deterministically as the probability, at the end of the process, that the food unit contains  1 bacteria. The number of bacteria per contaminated unit is evaluated stochastically, using one random draw from a discrete distribution per iteration. As an illustrative example, assume one 100 g food unit contains nj ¼ 6 cfu (i.e., a concentration of 1.22 log10 cfu/g). In this case, bacterial growth of 3.2 log10 cfu would lead to no change in prevalence, while the number of bacteria in the contaminated units would be determined through a random draw from a negative binomial distribution used to characterize growth (See next section), leading for instance to a value of Njþ1 ¼ 9822 cfu in the unit for one iteration. The expected number of bacteria2 in positive units for this specific process would be 9509.4 cfu for a 100 g unit. 2.1.3. Method #3 The third approach considers the prevalence (same definition as in method #2) and the expected number of bacteria per contaminated unit (i.e., in food units containing  1 bacteria) at each step of the food pathway. The expected number of bacteria in contaminated units is a rational number 1 (e.g., 1.2 or 65.8; unit: cfu/food unit). For any given process defined by a given parameter, the prevalence is evaluated deterministically as the expected probability that the food unit contains 1 bacteria. The expected number of bacteria per contaminated unit is evaluated deterministically. As an illustrative example, consider food units contaminated with prevalence Pj ¼ 10% and E[Nj] ¼ 6 cfu/100 g food unit at step j. Bacterial growth of 3.2 log10 would lead to no change in prevalence and an expected number of bacteria of E [Njþ1] ¼ 6  103.2 ¼ 9509.4 cfu in the 100 g contaminated units.

2 The expected value E[X] of random variable x is the weighted average of all possible values of x weighted by the probability that x assumes it.

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

2.2. Derivation of functions to characterize basic processes for each modeling approach We derived functions to model changes in bacterial contamination in food units as a result of partitioning, mixing, growth, inactivation and partial redistribution (Table 2). A basic process is characterized by three parameters at its initial state e the food unit size (U0), the prevalence (P0) and the bacterial level (C0, N0 or E[N0], depending on the method) e and by one or two parameters characterizing the impact of the basic process (e.g., log10 cfu increase and maximal population density for bacterial growth, redistribution factor for partial redistribution, etc.). The functions are used to estimate the new unit size (U1), the new prevalence (P1) and the new bacterial level (C1, N1 or E[N1], depending on the method) after the basic process. The impacts of the basic processes on unit size, prevalence, and bacterial level were characterized using several implicit assumptions. For all of the basic processes, the assumption of homogeneous cell distribution applies, indicating that, if a food unit is contaminated, the bacterial cells are homogeneously distributed within the unit without clustering effect. Additional assumptions will be specified below where appropriate. 2.2.1. Partitioning In this basic process, each food unit is split, starting from a unit of size U0 and generating U0/U1 units of smaller size U1. The basic process may be characterized by a single value U1 or by a probabilistic distribution of values U1. Practical examples of partitioning include distributing cheese curd to cheese hoops or filling cheese curd into waxed wooden boxes during cheese making (Bemrah et al., 1998), cutting frozen minced fish blocks into fish sticks before battering and breading for frying, partitioning of pasteurized milk from a tank to bottles (Latorre et al., 2011), or division of a food into individual servings at the time of consumption (Nauta et al., 2009). Method #1: Assuming homogeneous contamination, the expected concentration (in log10 cfu/g) after the process, C1, remains unchanged. Similarly, given our definition of prevalence, each uncontaminated food unit will lead to U0/U1 uncontaminated units, and each contaminated units will lead to U0/U1 contaminated units leaving the overall prevalence unchanged. Method #2: The number of bacteria at the beginning of the process, N0, will be distributed among the newly generated U0/U1

3

units. Assuming homogeneous contamination in the unit U0, the number of bacteria present in the newly formed food unit of size U1, given that this food is contaminated, may be simulated using one random draw from a binomial distribution with size parameter N0 and a probability U1/U0, restricted on [1, N0]. A binomial distribution restricted on [1, N0] is named a positive binomial distribution (see Appendix for some characteristics of the positive binomial distribution). At low bacterial concentrations, the prevalence following this process may change, as some newly generated food units may contain exactly 0 bacteria. The probability of having exactly 0 bacteria in the new unit of size U1 is Prob(n ¼ 0) for n ~ Binomial(N0, U1/U0), that is ð1  U1 =U0 ÞN0 . The new expected prevalence is then P1 ¼ P0 ð1  ð1  U1 =U0 ÞN0 Þ. Method #3: The new prevalence P1 is derived as in method #2. The expected number of bacteria on a contaminated unit is the expected value of the binomial distribution with size parameter N0, a probability U1/U0 restricted on [1, N0], that is N0  U1 =U0  ð1  ð1  U1 =U0 ÞN0 Þ1 . 2.2.2. Mixing In this basic process, food units of size U0 are combined to form new units of larger size U1. Examples of mixing are mixing shredded lettuce on conveyor belts or in a flume tank (Danyluk and Schaffner, 2011) and the addition of milk from multiple bulk tanks to a tanker truck (Albert et al., 2005). U1 can be characterized by a single value or by a distribution of values. In the interest of parsimony, we here restrict U1 to be a multiple of U0, but this assumption can easily be relaxed. This process may be modeled in various ways, depending on the underlying assumptions regarding the distribution of the contamination from units to units (Nauta, 2005). In this study, we consider a simple case of mixing U1/U0 contaminated and uncontaminated units, where all contaminated units are similarly contaminated (i.e., same average contamination level for method #1, same number of bacteria for method #2, or same expected number of bacteria for method #3). The basic process of mixing has the same impact on prevalence, which is characterized by the same equation (Table 2), for all three modeling methods. Method #1: The prevalence will be impacted. The new proportion of uncontaminated food products is equal to the proportion of new food units of size U1 consisting only of uncontaminated units (i.e., with mean concentration equal to 0 cfu/g). The new prevalence is then given by 1  ð1  P0 ÞU1 =U0 . It takes (P0U0/U1)

Table 2 Basic process functions according to the approach. Basic process Partitioning from units of size U0 to units of size U1

Method #1 P1 ¼ P0 C1 ¼ C0

Method #2

Method #3 N0

P1 ¼ P0 ð1  ð1  U1 =U0 Þ Þ N1 ~ Binom(N0,U1/U0), with N1  1a

P1 ¼ P0 ð1  ð1  U1 =U0 ÞN0 Þ N U 0 1 N !

N1 ¼

0

U0

Mixing from units of size U0 to units of size U1

P1 ¼ 1  ð1  P0 ÞU1 =U0 C1 ¼ C0 þ log10(P0/P1)

Growth of G log10 cfu

P1 ¼ P0 C1 ¼ min(C0 þ G, Cmax)

Inactivation of D log10 cfu

P1 ¼ P0 C1 ¼ C0  D

Partial Redistribution over R units

minðM; P01 Þ

M0 ¼ P1 ¼ M'  P0 C1 ¼ C0  log10 ðM'Þ

1

U

1U1 0

P1 ¼ 1  ð1  P0 ÞU1 =U0 N1 ~ N0  n n ~ Binom(U1/U0,P0), with n  1a

P1 ¼ 1  ð1  P0 ÞU1 =U0 N1 ¼ 0 N0 P0 U1 1

P1 ¼ P0 N1 ~ N0 þ NegBin(N0,exp(G  log(10))) with N1 < Nmax Nmax ¼ U0  10Cmax P1 ¼ P0 ð1  ð1  10D ÞN0 Þ N1 ~ Binom(N0,10D),with N1  1a

P1 ¼ P0 N1 ¼ minðN0  10G ; U0  10Cmax Þ

M0 ¼ minðM; N0; P01 Þ  P1 ¼ M 0  P0  1  1   N1 ~ Binom N0 ; M 01 ;

N0   M 01 with N1  1a

B

U1

C

U0 @1ð1P0 ÞU0 A

P1 ¼ P0 ð1  ð1  10D ÞN0 Þ N1 ¼

N0 10D 1ð110D ÞN0

M0 ¼ minðM; N0; P01 Þ  N0  P1 ¼ M 0  P0  1  1  M 01 N1 ¼

N0 N M0 ð1ð1M01 Þ 0 Þ

Key: C: Concentration, P: Prevalence, U: Unit size. Index 0: pre-process. Index 1: post-process. a i.e. random draw from a positive binomial distribution.

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

4

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

units with a mean number of bacteria ðU0  10C0 Þ for the process to obtain P1 contaminated unit of size U1 with a mean number of bacteria 10C1 . The new average concentration is then C1 ¼ C0 þ log10(P0/P1). Method #2: Assuming a random distribution of a proportion of P0 contaminated units and (1  P0) uncontaminated units at the beginning of the process, n, the number of initially contaminated units contributing to the final contaminated unit follows a binomial distributions with size parameter U1/U0 and probability parameter P0 restricted on [1, U1/U0]. This distribution is positive binomial because at least one initial contaminated unit is needed to obtain a final contaminated unit. The number of bacteria in each ultimately contaminated unit may then be simulated (in our design) using N1 ~ N0  n where n ~ Binom(U1/U0,P0), n  1. The prevalence of contaminated units is modified accordingly to reflect the gathering of U1/U0 uncontaminated units. We have P1 ¼ 1  ð1  P0 ÞU1 =U0 . Method #3: The prevalence is evaluated similarly to method #2. The new expected number of bacteria in contaminated units is N0 times the expected values of n, that is N1 ¼ N0 P0 ðU1 =U0 Þð1  ð1  P0 ÞU1 =U0 Þ1 . 2.2.3. Growth Bacterial growth of G log10 occurs in this process. The process can be characterized by a single value or a distribution of values G. We assumed a maximal population concentration of Cmax log10 cfu/ g. The prevalence and the size of the unit remain unchanged for all three modeling methods. Method #1: The new concentration in log10 cfu/g is simply min(C0 þ G, Cmax). Method #2: During the exponential growth phase, given the assumptions of i) asexual replication (e.g., by cell division), ii) independence among individual bacterial cells, iii) multiplication being a Poisson process in time, and iv) complete uniformity among all individuals in the population, starting with N0 individuals (cfu), we have a Yule (1925) pure birth growth process. We have

N1 eN0 þ NegBinðN0 ; expð  G  logð10ÞÞÞ where NegBin(n, p) is the negative binomial distribution with size parameter n and probability parameter p (Vose, 2008). N1 is limited to be below ðU0  10Cmax Þ. Method #3: The expected number of bacteria in a contaminated  unit is simply min N0  10G ; U0  10Cmax . 2.2.4. Inactivation Bacterial inactivation of D log10 occurs in this process. Inactivation may occur naturally such as during storage of low water activity food (Burnett et al., 2000; Edelson-Mammel et al., 2005) and lemon and lime juices (Enache et al., 2009); or obtained through treatments such as heating of juices (Enache et al., 2006), acidification of salad dressing (Beuchat et al., 2006), and high pressure processing of oysters (Calci et al., 2005). The process can again be characterized by a value or a distribution of values D. We assume that the inactivation process is applied homogeneously and independently on each bacterial cell contained in a given food unit. Method #1: The prevalence remains unchanged according to the definition of the prevalence in this method. The concentration is decreased by D, that is C1 ¼ C0  D. Method #2: The number of units with exactly 0 bacteria may possibly increase. Given a decrease of D log10, the probability of death for a bacterial cell is (1  10D). The probability that all bacteria are removed under the considered assumptions is then ð1  10D ÞN0 and the new prevalence is given by P1 ¼ P0 ð1  ð1  10D ÞN0 Þ. The number of bacteria on contaminated units is determined using one draw from a binomial

distribution with size parameter N0 and probability parameter 10D restricted on the domain [1, N0]. Note that, if initial bacterial contaminations are relatively small and D is large, as frequently observed in heating processes, most of the iterations will lead to a value of 1 cfu in the contaminated units. Method #3: The new prevalence is simulated similarly as in method #2. The expected number of bacteria in contaminated units is given by N0  10D  ð1  ð1  10D ÞN0 Þ1. 2.2.5. Partial redistribution Description of the process: In a partial redistribution, it is assumed that each food unit redistributes its bacteria amongst M > 1 food units. If M equals for example, 2.3, each contaminated unit will lead to (M  1) ¼ 1.3 newly generated contaminated units. This process allows simulation of cross-contamination among food units (Danyluk and Schaffner, 2011). An example of partial redistribution is the contact among contaminated and uncontaminated food units during processing and handling. It is characterized by a single value or a distribution of values M. Multiple other ways to simulate such cross-contamination have been described, depending on the exact process and the underlying hypotheses (Aziza et al., 2006; Hoelzer et al., 2012; Nauta, 2008). Method #1: Because the final prevalence cannot exceed 1, the actual redistribution factor M0 is then M0 ¼ minðM; P01 Þ and the final prevalence is M'  P0. The average concentration is reduced to C1 ¼ C0  log10(M'). Method #2: The number of bacteria on a contaminated unit of food after redistribution follows a binomial distribution with size N0 and probability M0 1, restricted on the [1, N0] domain. The   N0  , and the new probability of having 0 bacteria is 1  1  M 01   N0  0 01 prevalence is given by P1 ¼ M  P0  1  1  M . Method #3: The prevalence is evaluated similarly as in method #2. The expected contamination of positive units is the expected value of the previously derived binomial distribution, given by    N0 1 N0 M0 1  1  M 01 . 2.3. Dose-response model and expected number of cases An exposure assessment using method #1 will lead to an average concentration per serving. Method #2 will lead to a number of bacteria per contaminated serving, and method #3 to an expected number of bacteria per contaminated serving. The doseeresponse should be expressed in a way consistent with the way in which bacterial levels are expressed. Under the assumptions of an exponential dose-response model (Haas et al., 1999), the probability of infection following the ingestion of a contaminated product (C > 0 cfu/g) should be estimated as Pr(infjd, cont) ¼ 1  exp(r  E[d]) for methods #1, where E[d] is the expected dose estimated by using the concentration of bacteria under a Poisson distribution assumption and r is the exponential doseeresponse parameter. Pr(infjd, cont) ¼ 1  (1  r)d should be used with method #2, for an exact dose of d bacteria (Haas, 2002). Under the assumption of a positive Poisson distribution of bacteria in contaminated product (d  1) with mean E[d] in Method #3, the probability of infection following the ingestion of a contaminated product should be estimated as

Prðinfjd; contÞ ¼

1  expðr  E½dÞ : 1  expðE½dÞ

Similarly, if a beta-Poisson dose-response with parameters a and b is considered (Haas et al., 1999), the doseeresponse translates to Pr(infjd,cont) ¼ 1  (1 þ E[dose]/b)a for products in method #1,

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

where E[dose] is the expected concentration of bacteria under a Poisson distribution assumption (Furumoto and Mickey's (1967) approximation), and to Pr(infjd,cont) ¼ 1  G(a þ b)G(b þ d)/ (G(b)G(a þ b þ d)) in method #2 where d is the actual number of bacteria in the product (Haas, 2002). Under the assumption of a positive Poisson distribution of bacteria in contaminated product (d  1) with mean E[d] in Method #3, the probability of infection following the ingestion of a contaminated product if a beta-Poisson dose-response with parameters a and b is considered should be estimated as

Prðinfjd; contÞ ¼

a  1  1 þ E½dose b 1  expðE½dÞ

:

The expected number of case X following the ingestion of S servings may be estimated as

X ¼ S  P  E½Prðinf jd; contÞ:

(1)

2.4. Simulation study In order to compare the three modeling approaches, we simulated 10,000 food pathways at random and evaluated the predicted risk of illness from each of these food pathways using the three methods. Each food pathway was simulated as a succession of the 5 basic processes described above (i.e., partitioning process, mixing, growth, inactivation and redistribution), utilizing each process exactly once but in a randomized order. The initial prevalence varied from simulation to simulation between 0 and 0.1. Each simulation started with 10,000 food units (each unit 1000 g). The initial bacterial concentrations followed a unit-to-unit triangular distribution of variability, with parameter min ¼ 0 þ x, mode ¼ 1 þ x, max ¼ 3 þ x log10 cfu/g, with x uniformly sampled on (5, 1) log10 cfu/g as simulation-to-simulation variability. This process allowed the simulation of various initial conditions ranging from low (down to 5 log10 cfu/g) to high (up to 4 log10 cfu/g) initial contamination levels, with a reasonable range of contamination variability between food units (3 log10 cfu/g). For methods #2 and #3, the initial prevalence and the initial concentration were translated as initial prevalence (probability that  1 bacteria is present in a unit) and number of bacteria per food unit, on the basis of a Poisson distribution. The five basic processes are modeled to vary from simulation to simulation, and from food unit to food unit within a given simulation as detailed in Table 3. The range of parameters used for each basic process represents a realistic set of conditions that may occur Table 3 Parameters from the simulation study. Characteristics

Simulation-to-simulation variability

Food-unit-to-food-unit variability in the basic process parameter

Initial Prevalence Initial level of contamination Partitioning

x ~ Uniform (0, 0.1) x ~ Uniform (5, 1) log10 cfu/g x ~ Uniform (1, 1000) g

Mixing

x ~ Uniform (1, 100) units

Growth

x ~ Uniform (0, 8) log10 cfu

Inactivation

x ~ Uniform (0, 8) log10 cfu

Partial Redistribution

x ~ Uniform (1, 5) new contaminated units

P0 ¼ x (no variability) C0 ¼ Triangular (x, xþ1, xþ3) log10 cfu/g U1 ¼ Uniform (min(U0, x), U0) g U1 ¼ Uniform (x5, xþ5)  U0 grams G ¼ Uniform (x1, xþ1) log10 cfu D ¼ Uniform (max(x1, 0), xþ1) log10 cfu M ¼ Uniform (max(x1), xþ1) new contaminated units

U0 unit size before the basic process.

5

in food pathways. The same set of parameters is used for all methods to avoid any difference linked to the variability in the input parameters. In each simulation of one food-production pathway, the mean risk per serving is estimated. Given the simulation of 10,000 pathways, 10,000 estimates of mean risk per serving are obtained for each of the modeling methods. This simulation study allowed for a comparison of the three modeling methods using a broad set of food pathways. As one illustrative example of doseeresponse, a beta-Poisson dose-response with parameters a ¼ 0.248 and b ¼ 48.80 was considered, adapted from Teunis et al. (2008) for E. coli O157:H7. This dose-response model is characteristic in its S shape of some other doseeresponse models commonly used in microbial risk assessment. The serving size was 100 g in all simulations. A sensitivity analysis was performed to characterize which simulated pathways and which sets of parameters led to large discrepancies among methods. For that purpose, classification trees (Frey et al., 2003; Ripley, 1996) were built to find factors that explain a >10-fold difference in the illness risk estimates between methods. The impact of the sequence of basic processes was explored by incorporating all 10 possible sequences for 2 given basic processes (e.g. “inactivation before growth”, “redistribution before mixing”, etc.) and all 36 possible sequences for 3 basic processes (e.g. “inactivation before growth before redistribution”, “inactivation before redistribution before growth”, “mixing before growth before inactivation”, etc.). All simulations and Monte-Carlo integrations were performed in R (© The R Core Team) using the mc2d package (Pouillot and DelignetteMuller, 2010) and the rpart package (Therneau et al., 2013). 3. Results Fig. 1a shows the scatter plot of the mean risk estimated based on 10,000 simulated food-production pathways, comparing method #1 vs. method #2. For most of the simulated pathways, the estimates based on the simple concentration method (method #1) are similar to or at least in the same order of magnitude as the estimates obtained using the more complex method #2, based on the simulation of the number of bacteria per food unit. Yet, for some simulations, method #1 provides an overestimation of the estimated risk. This overestimation exceeds 1 log10 in 720 out of 10,000 simulations (7.2%) and exceeds >3 log10, that is a 1000-fold overestimation, in 37 out of 10,000 simulations (0.37%). Fig. 1b compares method #2 and method #3, and suggests that these two methods lead to comparable estimates in all evaluated simulations, even though method #2 simulated additional variability (i.e., accounted for random variability of the stochastic processes) while method #3 did not include this variability (and thus is a simpler modeling choice). When the estimated risk is very low, method #3 may slightly overestimate the risk. The classification tree developed to characterize simulations for which method #1 overestimates the average risk by a factor of 10 or more (720 simulations out of 10,000) provides some insights as to when such overestimation may occur (Fig. 2). Large discrepancies are predominantly observed when the pathway includes considerable growth after a large inactivation event, and when the initial concentration was low. 4. Discussion Various approaches to develop QMRA and various models representing basic processes have been described in the literature. Approaches reflective of method #1, thus based on concentrations and used with or without separate estimation of the prevalence, are very popular because these types of equations are mathematically

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

6

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

Fig. 1. Log10 of the mean risk per serving estimated using method #1 (top) or method #3 (bottom) plotted against the log10 of the mean risk per serving estimated using method #2 from 10,000 random simulated pathways. On each graph, the dashed lines represent 1  5 log10 differences in the results.

relatively simple. In this method, an inactivation process leads to a simple subtraction; a growth process is a simple addition. In published QMRA framework proposals, it is frequently unclear if finite entities (i.e., bacteria) rather than concentrations should be used (Nauta, 2008). In published risk assessments, it is usually not explicitly specified if concentrations or numbers of bacteria are used. Because bacterial cells are finite entities, it is more rigorous to develop equations based on natural numbers rather than on rational numbers. Additionally, some basic processes such as cross-contamination are more easily handled using finite entities (Hoelzer et al., 2012). In this study, we considered method #2 as the ‘gold standard’ because modeling number of bacteria is more rigorous than modeling concentrations. We developed sets of equations associated with five basic processes typical in food production, processing and handling (Table 2), modeling either concentration (method #1), number of

bacteria (method #2) or expected number of bacteria in positive food units (method #3). The equations were based on probabilistic theory (Vose, 2008). The equations for methods #3 are new, while those for methods #1 and #2 are very similar to equations that have been reported previously (FDA-iRISK, 2012; Nauta, 2008). These equations can be used to model various food pathways. However, the user must be aware of the various assumptions underlying those equations, and evaluate whether key assumptions are met before incorporation in risk assessments. In favor of parsimony, we chose the simplest models, assuming, as an example, complete homogeneity of the contamination within a food unit at each step. However, there are multiple other approaches (FDA-iRISK, 2012; Nauta, 2008, 2005) to simulate the same basic processes, some of which may be more appropriate for modeling a given situation than the simplified approaches provided here. As an example, the mixing can be modeled in various ways depending on what units are actually mixed together and how variability is considered in the Monte-Carlo simulation. If the variability modeled in the simulation reflects only product-to-product variability within a process line, the mixing process could be simulated by randomly sampling values from different iterations of the Monte-Carlo simulation. In other situations, if other dimensions of variability (e.g., temporal, geographic) are modeled, this latter procedure may be flawed if mixing simply could not occur in reality, for example between simulated food units located in another place in the country, or generated during a different year. Similarly, clustering can have a very important impact in QMRAs but is rarely considered; we did not consider clustering in the equations provided here. Homogeneity of the contamination within a food unit is an assumption that might be unsubstantiated by the available data (ILSI, 2010). Similarly, bacterial growth leads, by definition, to spatial clusters around the original cells. This fact, however, is very rarely considered in risk assessments (Sanaa et al., 2004). The model and the derivation of the equations for the basic processes should then be developed on a case-by-case basis and the equations provided in Table 2 should be regarded as illustrative examples only. The equations for the various processes were developed so that the expected overall contamination (i.e. P  U  10C for method #1 or P  N for method #2 and method #3) was similar, at each step, across the three compared methods. This condition was necessary to permit comparisons across results. To test this condition, we confirmed that the average number of bacteria using the various methods were similar in all simulations. The impact of the chosen approach on the predicted risk of illness is clearly a function of the evaluated food pathway. Evaluation of one or few case studies would not have provided sufficient information to generalize the observations and to draw overarching conclusions. We chose to simulate a large number of plausible food pathways, using five basic processes (growth, inactivation, mixing, partitioning and redistribution) in a random order, with a random set of parameters for each of these processes. Classification trees were then a natural tool to handle the corresponding sensitivity analysis (Frey et al., 2003), notably for the evaluation of the impact of the order of the basic processes on the estimates. The results provide some insights as to the major factors leading to discrepancies between method #1 and method #2, which are notably a large bacterial inactivation step followed by a large bacterial growth step. While the expected number of bacteria in food units was similar for all methods, the mean risk per serving was somewhat different and, in some situations, largely overestimated by method #1 when compared to method #2. The classification tree helped to better characterize the most extreme situations, where method #1 overestimated the average risk per serving by a factor of more than 10

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

7

Fig. 2. Classification tree for 10,000 simulations of random pathways. “High” characterizes a simulation for which method #1 overestimates the mean risk per serving by a factor of 10 or more when compared to method #2. At each node, the left child is “yes” to the associated question on the root, and the right child is “no”. The classification is proposed at each leaf, and is boxed when “High” is the majority.

when compared to method #2. Note that the average risk per serving is proportional to the estimated number of cases (see Eq. (1)), so that an overestimation of the average risk by a factor of 10 would lead to an overestimation of the expected number of cases by a similar factor. The classification tree suggests that most discrepancies are observed when a large inactivation occurs in the process, followed by a large increase. Actually, in this situation, method #1 estimates a low concentration homogeneously across all previously contaminated food units and the prevalence is

unchanged; while method #2 estimates a reduction in prevalence, resulting in only a few highly contaminated food units (i.e., high risk from contaminated units in very rare situations). The resulting average risk, estimated using a classical sigmoidal doseeresponse function, leads to a higher estimated risk in method #1 compared to method #2. This point can be illustrated by a simple yet extreme, example (Fig. 3). Assume a hypothetical situation with a starting mean concentration of 100 cfu in a 100 g unit, i.e. a concentration of

Fig. 3. Extreme scenario illustrating method #1 and method #2 discrepancies.

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

8

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9

0 log10 cfu/g. An 8 log10 reduction (e.g. a heating process), followed by an 8 log10 growth (e.g. due to storage at temperature abuse on a product that supports bacterial growth) are applied. After the inactivation step, method #1 will consider a very low mean concentration in all contaminated food units. This mean level of concentration actually corresponds to the absence of bacterial contamination in most of the food units and some low numbers of bacteria in some very rare contaminated food units. However, method #1 does not specifically consider this variability and the following growth will be modeled from this very low concentration for all impacted food products. Overall, method #1 will lead to final estimates of 100% prevalence with 0 þ 8e8 ¼ 0 log10 cfu/g. Note that those final values are the same as in the initial step, as if the processes had no effects: if SR þ SI ¼ 0, where R is the total cumulative reduction of the hazard and I is the total increase of the hazard, the initial situation is the same as the final one using this method. The corresponding estimated probability of illness, using the E. coli O157:H7 doseeresponse model (Teunis et al., 2008), would be 0.24. Method #2 does consider the heterogeneity in food units, and the growth in food units with 0 bacteria after inactivation will have no impact. Here too, SR þ SI ¼ 0, but the final situation (with regard to the prevalence or the bacterial number) is not the same as the initial one using this method. The method would lead to a final prevalence of approximately 106, with an estimated concentration of 6 log10 cfu/g (N ¼ 108 cfu in 100 g) in positive food units, that is a very low prevalence, with a few highly contaminated units. The resulting risk will be 9.7  107. Most of the discrepancies between method #1 and method #2 are associated with this extreme variability of the contamination, which is explicitly modeled in method #2 but is averaged in method #1. Some selected food processes do include a high inactivation process step, followed by a potentially high level of bacterial growth. Examples include microbial food spoilage in canned food, L. monocytogenes in pasteurized fluid milk or in cheese made from pasteurized-milk, or Salmonella in low-moisture food that is added into a product that supports growth such as adding spices to a deli salad product that supports bacterial growth. Based on the extensive simulation results presented here, it appears that method #3 is comparable to method #2. In general, method #3 may thus be preferable due to the greater simplicity of the required modeling. Indeed, method #2 may in certain situations be complex if not impossible to derive. As an example, the simulation of a positive binomial distribution may be complicated, notably if the survival probability of one bacterial cell is very low, as observed in food processing (See Appendix). Alternatively, results from Fig. 1b suggest that method #3, estimating the expected number of bacteria in contaminated food units, is comparable to method #2. Method #3 is easier to implement, because no random draws are needed, and represent a potentially appealing alternative to model#1 and model #2 in situations where expectations can be readily evaluated. Additional studies, such as using more extensive consideration of basic processes and using different variability distribution or modeling concentration in cfu/g, should be developed to confirm broader applicability of this finding. 5. Conclusions In certain situations, the choice of how to model bacterial contamination in a QMRA can have a clear impact on modeling results and model complexity, even though in many of the cases studied here all three evaluated approaches generated highly similar results. A better understanding of those consequences will enable judicious selection of the most adequate modeling framework. Modeling prevalence and the number of bacteria in positive units, or prevalence and expected number of bacteria in positive

units, may be preferable to modeling concentrations defined on a continuous scale, notably in food pathways including large bacterial inactivation (e.g. pasteurization), followed by bacterial growth. Acknowledgments This work was supported by the U.S. Food and Drug Administration and, in part, by appointments to the Research Participation Program at the Center for Food Safety and Applied Nutrition administered by the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and the U.S. Food and Drug Administration. Appendix The positive binomial distribution is a binomial distribution restricted on [1, n], that is a binomial distribution with a probability of zeros being equal to 0. If x ~ PosBinom(n, p) with n the number of observation and p the probability, the density of this positive binomial distribution is

f ðx; n; pÞ ¼

8 > > > > < > > > > :

n

!

px ð1  pÞnx x 1  ð1  pÞn 0

x ¼ 1; …; n and 0 < p < 1 x¼0

we have: E[x] ¼ np/(1  (1  p)n) Given p(x, n, p) the probability mass function and F1(x; n, p) the quantile function of a binomial distribution Binom(n, p). A method to generate random values from the corresponding positive binomial is to use F1(u; n, p) where u is sampled from a uniform(p(0; n, p), 1). Simulation difficulties occur when n is small and p is close to 0, as observed in some inactivation processes, because of numerical precision issues. References Albert, I., Pouillot, R., Denis, J.B., 2005. Stochastically modeling Listeria monocytogenes growth in farm tank milk. Risk Anal. 25, 1171e1185. Aziza, F., Mettler, E., Daudin, J.J., Sanaa, M., 2006. Stochastic, compartmental, and dynamic modeling of cross-contamination during mechanical smearing of cheeses. Risk Anal. 26, 731e745. Bemrah, N., Sanaa, M., Cassin, M.H., Griffiths, M.W., Cerf, O., 1998. Quantitative risk assessment of human listeriosis from consumption of soft cheese made from raw milk. Prev. Vet. Med. 37, 129e145. Beuchat, L.R., Ryu, J.H., Adler, B.B., Harrison, M.D., 2006. Death of Salmonella, Escherichia coli O157 : H7, and Listeria monocytogenes in shelf-stable, dairybased, pourable salad dressings. J. Food Prot. 69, 801e814. Burnett, S.L., Gehm, E.R., Weissinger, W.R., Beuchat, L.R., 2000. Survival of Salmonella in peanut butter and peanut butter spread. J. Appl. Microbiol. 89, 472e477. Calci, K.R., Meade, G.K., Tezloff, R.C., Kingsley, D.H., 2005. High-pressure inactivation of hepatitis A virus within oysters. Appl. Environ. Microbiol. 71, 339e343. Cassin, M.H., Lammerding, A.M., Todd, E.C., Ross, W., McColl, R.S., 1998. Quantitative risk assessment for Escherichia coli O157:H7 in ground beef hamburgers. Int. J. Food Microbiol. 41, 21e44. Chen, Y., Dennis, S.B., Hartnett, E., Paoli, G., Pouillot, R., Ruthman, T., Wilson, M., 2013. FDA-iRISK-A comparative risk assessment system for evaluating and ranking food-hazard pairs: case studies on microbial hazards. J. Food Prot. 76, 376e385. Danyluk, M.D., Schaffner, D.W., 2011. Quantitative assessment of the microbial risk of leafy greens from farm to consumption: preliminary framework, data, and risk estimates. J. Food Prot. 74, 700e708. Edelson-Mammel, S.G., Porteous, M.K., Buchanan, R.L., 2005. Survival of Enterobacter sakazakii in a dehydrated powdered infant formula. J. Food Prot. 68, 1900e1902. Enache, E., Chen, Y.H., Awuah, G., Economides, A., Scott, V.N., 2006. Thermal resistance parameters for pathogens in white grape juice concentrate. J. Food Prot. 69, 564e569. Enache, E., Chen, Y.H., Elliott, P.H., 2009. Inactivation of Escherichia coli O157:H7 in single-strength lemon and lime juices. J. Food Prot. 72, 235e240.

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008

R. Pouillot et al. / Food Microbiology xxx (2014) 1e9 FDA-iRISK, 2012. FDA-iRISK(TM) Technical documentation (version 1.0). Retrieved May, 4th, 2014, from. http://irisk.foodrisk.org/SupportingDocumentation.aspx. FDA/FSIS, 2003. Quantitative Assessment of Relative Risk to Public Health from Foodborne Listeria monocytogenes Among Selected Categories of Ready-to-eat Foods. Food and Drug Administration, United States Department of Agriculture, Centers for Disease Control and Prevention, p. 541. http://www.fda.gov/ Food/FoodScienceResearch/RiskSafetyAssessment/ucm183966.htm. Frey, H.C., Mokhtari, A., Danish, T., 2003. Evaluation of Selected Sensitivity Analysis Methods Based upon Applications to Two Food Safety Process Risk Models. Office of Risk Assessment and Cost-Benefit Analysis, U.S. Department of Agriculture, Washington, DC. http://www.ce.ncsu.edu/risk/Phase2Final.pdf. Furumoto, W.A., Mickey, R., 1967. A mathematical model for the infectivity-dilution curve of tobacco mosaic virus: theoretical considerations. Virology 32, 216e223. Haas, C.N., 2002. Conditional dose-response relationships for microorganisms: development and application. Risk Anal. 22, 455e463. Haas, C.N., Rose, J.B., Gerba, C.P., 1999. Quantitative Microbial Risk Assessment. Wiley, New York. Hoelzer, K., Pouillot, R., Gallagher, D., Silverman, M.B., Kause, J., Dennis, S., 2012. Estimation of Listeria monocytogenes transfer coefficients and efficacy of bacterial removal through cleaning and sanitation. Int. J. Food Microbiol. 157, 267e277. ILSI, 2010. Impact of microbial distributions on food safety. In: ILSI Europe Risk Analysis in Food Microbiology Task Force. Report Series, Brussels, Belgium. http://www.ilsi.org/Europe/Publications/Microbial%20Distribution%202010.pdf. Lambertini, E., Danyluk, M.D., Schaffner, D.W., Winter, C.K., Harris, L.J., 2012. Risk of salmonellosis from consumption of almonds in the North American market. Food Res. Int. 45, 1166e1174. Latorre, A.A., Pradhan, A.K., Van Kessel, J.A.S., Karns, J.S., Boor, K.J., Rice, D.H., Mangione, K.J., Grohn, Y., Schukken, Y.H., 2011. Quantitative risk assessment of listeriosis due to consumption of raw milk. J. Food Prot. 74, 1268e1281. Nauta, M., 2008. The Modular Process Risk Model (MPRM): a structured approach to food chain exposure assessment. In: Schaffner, D.W. (Ed.), Microbial Risk Analysis of Foods. ASM Press, Washington, D.C., pp. 99e136.

9

Nauta, M., Hill, A., Rosenquist, H., Brynestad, S., Fetsch, A., van der Logt, P., Fazil, A., Christensen, B., Katsma, E., Borck, B., Havelaar, A., 2009. A comparison of risk assessments on Campylobacter in broiler meat. Int. J. Food Microbiol. 129, 107e123. Nauta, M.J., 2005. Microbiological risk assessment models for partitioning and mixing during food handling. Int. J. Food Microbiol. 100, 311e322. Pouillot, R., Delignette-Muller, M.L., 2010. Evaluating variability and uncertainty separately in microbial quantitative risk assessment using two R packages. Int. J. Food Microbiol. 142, 330e340. Rigaux, C., Andre, S., Albert, I., Carlin, F., 2014. Quantitative assessment of the risk of microbial spoilage in foods. Prediction of non-stability at 55 degrees C caused by Geobacillus stearothermophilus in canned green beans. Int. J. Food Microbiol. 171, 119e128. Ripley, B.D., 1996. Pattern Recognition and Neural Networks. Cambridge University Press, Cambridge. Sanaa, M., Coroller, L., Cerf, O., 2004. Risk assessment of listeriosis linked to the consumption of two soft cheeses made from raw milk: Camembert of Normandy and Brie of Meaux. Risk Anal. 24, 389e399. Schroeder, C.M., Latimer, H.K., Schlosser, W.D., Golden, N.J., Marks, H.M., Coleman, M.E., Hogue, A.T., Ebel, E.D., Quiring, N.M., Kadry, A.R., Kause, J., 2006. Overview and summary of the Food Safety and Inspection Service risk assessment for Salmonella enteritidis in shell eggs, October 2005. Foodborne Pathogens Dis. 3, 403e412. Teunis, P.F., Ogden, I.D., Strachan, N.J., 2008. Hierarchical dose response of E. coli O157:H7 from human outbreaks incorporating heterogeneity in exposure. Epidemiol. Infect. 136, 761e770. Therneau, T., Atkinson, B., Ripley, B., 2013. rpart: Recursive Partitioning. R package version 4.1-4. Retrieved Dec 24, 2013, from. http://CRAN.R-project.org/ package¼rpart. Vose, D., 2008. Risk Analysis: a Quantitative Guide, third ed. Wiley and Sons., Chichester, UK. Yule, G.U., 1925. The growth of population and the factors which control it. J. R. Stat. Soc. Ser. B 25, 1e58.

Please cite this article in press as: Pouillot, R., et al., Modeling number of bacteria per food unit in comparison to bacterial concentration in quantitative risk assessment: Impact on risk estimates, Food Microbiology (2014), http://dx.doi.org/10.1016/j.fm.2014.05.008