An Agent-Based Model for Pathogen Persistence and Cross-Contamination Dynamics in a Food Facility.

We used an agent-based modeling (ABM) framework and developed a mathematical model to explain the complex dynamics of microbial persistence and spread within a food facility and to aid risk managers in identifying effective mitigation options. The model explicitly considered personal hygiene practices by food handlers as well as their activities and simulated a spatially explicit dynamic system representing complex interaction patterns among food handlers, facility environment, and foods. To demonstrate the utility of the model in a decision-making context, we created a hypothetical case study and used it to compare different risk mitigation strategies for reducing contamination and spread of Listeria monocytogenes in a food facility. Model results indicated that areas with no direct contact with foods (e.g., loading dock and restroom) can serve as contamination niches and recontaminate areas that have direct contact with food products. Furthermore, food handlers' behaviors, including, for example, hygiene and sanitation practices, can impact the persistence of microbial contamination in the facility environment and the spread of contamination to prepared foods. Using this case study, we also demonstrated benefits of an ABM framework for addressing food safety in a complex system in which emergent system-level responses are predicted using a bottom-up approach that observes individual agents (e.g., food handlers) and their behaviors. Our model can be applied to a wide variety of pathogens, food commodities, and activity patterns to evaluate efficacy of food-safety management practices and quantify contamination reductions associated with proposed mitigation strategies in food facilities. Published 2018. This article is a U.S. Government work and is in the public domain in the USA. Risk Analysis published by Wiley Periodicals, Inc. on behalf of Society for Risk Analysis.

INTRODUCTION

Persistence of microbial pathogens in the food‐facility environment and potential for cross‐contamination of food products during preparation have been implicated as factors contributing to several well‐documented outbreaks of foodborne illness (Brown, Hoover, Selman, Coleman, & Schurz Rogers, 2017; Ferreira, Wiedmann, Teixeira, & Stasiewicz, 2014; Heiman et al., 2016; Jackson et al., 2011; Munther & Wu, 2013). In its recommendations for food‐safety plans, based on requirements of the Food Safety Modernization Act, U.S. Food and Drug Administration draft guidance specifies that pathogen persistence and spread within a food facility are controlled and managed predominantly by measures that provide basic environmental and operational conditions needed to produce safe and wholesome foods. Such measures, also known as prerequisite programs (PRPs), are based on the principles of hazard analysis and critical control points (HACCP) and may include proper cleaning and disinfection as well as personal hygiene practices, among other practices (Gaze, 2015; Food and Drug Administration [FDA], 2018). While risk factors contributing to microbial persistence and spread in food facilities have been identified (Valero, Rodriguez, Posada‐Izquierdo, & Peres‐Rodriguez, 2016; Wirtanen & Salo, 2007), their potential impacts on the likelihood and levels of contamination in prepared foods are not well understood.

Mathematical models, when grounded in observational and experimental data, can provide numerous insights into the relationships between food‐safety risks and contamination events. Such models can further facilitate evaluation and exploration of different mitigation options, in terms of their potential impacts on reducing the risk of foodborne illnesses (e.g., Duret et al., 2017; Ivanek, Grohn, Wiedmann, & Wells, 2004; Mokhtari & Jaykus, 2009; Mokhtari, Oryang, Chen, Pouillot, & Van Doren, 2018; Pouillot et al., 2015; Schaffner, 2004). The purpose of the work described here was to develop a mathematical model to explain the complex dynamics of microbial persistence and spread within a food facility and to provide a framework to aid risk managers in identifying effective mitigation options. To this end, we developed the FDA's Quantitative Microbial Risk Assessment Model for Food Facilities, hereafter referred to as F2‐QMRA. This model can serve as a virtual laboratory to (i) identify “hot spots” of potential contamination in different areas of a food facility (e.g., food processing area, restroom, loading dock) that can act as contamination niches during facility operation, (ii) examine food‐handler activities and behaviors that can result in microbial contamination of foods, and (iii) evaluate the effectiveness of candidate mitigation options aimed at reducing the likelihood and level of microbial contamination in the facility environment and prepared foods.

Challenges inherent in modeling a dynamic system, such as operations in a food facility, include (i) tracking persistence and spread of microbial pathogens in tandem with movements of food handlers and foods over time, (ii) accounting for spatial heterogeneity of microbial contamination and presence of contamination‐harborage sites in the facility environment, (iii) quantifying likelihood of cross‐contamination events over time and space during facility operation, and (iv) capturing temporal and spatial relationships among multiple activities that may occur during facility operation (e.g., unloading food ingredients, using restrooms, preparing foods) and the impact of these activities on spread of microbial contamination in the facility environment and to prepared foods.

F2‐QMRA addresses these challenges by adopting an agent‐based modeling (ABM) framework that simulates behaviors and activities of food handlers in a food facility, while explicitly tracking spatial distribution of microbial contamination on individual objects in the facility environment. ABM is a computer simulation technique used to study the interactions among people, objects, places, and time. It employs a bottom‐up approach to modeling complex systems, starting from the individual agents involved (Abar, Theodoropoulos, Lemarinier, & O'Hare, 2017). ABM differs from traditional, regression‐based methods in that it allows for exploration of complex systems that display dependence of individuals and feedback loops in causal mechanisms (El‐Sayed, Scarborough, Seemann, & Galea, 2012). Compared with discrete‐event simulation models that represent a top‐down modeling approach with simulation events that typically are reactive and limited in capabilities (Duret et al., 2017; Pouillot et al., 2015), ABM is more flexible and provides a suitable framework for dynamic systems in which agents frequently interact with each other and whose behaviors can change (Chan, Son, & Macal, 2010). ABM has gained increasing attention over the past decade, with numerous applications in modeling risks associated with infectious disease transmission and recent applications in modeling microbial cross‐contamination (Carley et al., 2006; Mokhtari et al., 2018; Smolinski, Hamburg, & Lederberg, 2003; Venkatramanan et al., 2018).

To demonstrate the utility of F2‐QMRA, we created a hypothetical case study and used it to compare different risk mitigation options for reducing contamination and spread of Listeria monocytogenes in a food facility with four distinct rooms: food processing area, common area (office), restroom, and loading dock. The overall structure of F2‐QMRA is intended to be flexible, so that with minor modifications it can readily accommodate a variety of facilities with different designs and purposes (e.g., food service, retail) as well as activities.

METHODS

Overview of the Model

F2‐QMRA, including the ABM simulation framework, was written in the open‐source language R version 3.4.2 (the R code is available on request: FDAFoodSafetyRiskModel@fda.hhs.gov). We used an ABM simulation framework to model a food facility with food handlers represented as individual agents. Each food handler participated in selected daily activities (e.g., unloading ingredients in the loading dock, preparing foods in the food processing area, among other activities) and demonstrated certain behaviors (e.g., personal hygiene practices). Interactions among food handlers and between food handlers and facility environment then were used to produce emergent, facility‐level dynamics that could not be deduced or forecasted by the observation of each individual food handler.

F2‐QMRA allows users to explicitly define the facility layout, including individual rooms and areas. Each food facility was modeled as a network, with nodes representing distinct areas (e.g., food processing area, office, and restroom). For each node, room‐specific objects were explicitly defined to allow simulation of potential cross‐contamination events between food handlers and objects during facility operation (Fig. 1). Each object within a room was defined using a set of attributes (Table I) and further assigned to one of the four sampling zones, based on the principles of environmental monitoring (PEM) zoning concept (Almond Board of California, 2010; Brouillette et al., 2014; FDA, 2017). These sampling zones represent different levels of risk, in terms of likelihood of transferring contamination to food products once microbial contamination is introduced to contact surfaces (Table II).

Figure 1

Example layout for a food facility that includes four distinct areas: food processing area, common area (office), restroom, and loading dock.

Table I

List of Attributes Defined for Each Specific Type of Object in F2‐QMRAa

Object‐Specific Attributes	Definition
Name	Specific object name (e.g., slicer)
Surface area	Object surface area (ft2)
Fraction of the surface area representing holes and cracks	Represents the ability of the object to harbor pathogens and reintroduce microbial contamination to the environment during facility operation regardless of the surface sanitation events (see Section 2.2.2 for further discussion)
Total number of grids	Number of 4″ × 4″ grid cells defined on the object surface area representing potential swab sample locations (see Section 2.2.1 for further discussion)
Total number of grids for cross‐contamination	Number of grid cells on the object involved in a cross‐contamination eventb
Sanitation and disinfection technique	Sanitation or disinfection method defined for the object; options include no sanitation, using surface sanitizer, and using surface disinfectantc
Zone id	Represents sampling zone id (see Table II)

Each specific object is individually tracked within the ABM framework.

Number of grid cells involved in a cross‐contamination event is less than the total number of grids defined for a specific object (see Section 2.2.5 for further discussion).

Surface sanitizers were only applied to food‐contact surface areas (e.g., cutting board) and meant to reduce, not eliminate, the occurrence and growth of microbial pathogens. Surface disinfectants were applied to non–food‐contact surface areas that were frequently touched to eliminate microbial pathogens.

Table II

Definition of Environmental Sampling Zones in a Food Facility

Sampling Zones	Definition	Examples
Zone 1	Areas in the facility that are direct food product contact surfaces before the product is sealed in a package	Utensil, food‐contact surface, conveyor belt, slicer
Zone 2	Non–product‐contact areas in the facility that are closely adjacent to food product contact surfaces	Non–food‐contact surfaces, equipment framework
Zone 3	Non–product‐contact surfaces that are in the processing area but not adjacent to Zone 1 surfaces; Zone 3 surfaces, however, have the possibility of leading to food product cross contamination	Floors and walls in the food processing area
Zone 4	Areas remote from food product processing areas; Zone 4 areas, if not maintained in good hygienic condition, can lead to cross‐contamination of Zones 1, 2, and 3	Restroom environment (e.g., toilet, sink), office area, loading dock

In view of the numerous scenarios and spatial resolution of the model, the F2‐QMRA code was written to be launched on parallelized processors using a high‐performance computing cluster (Office of Management, Center for Food Safety and Applied Nutrition, FDA, College Park, MD). Nonetheless, the code can be run on a desktop. The model was vectorized to simultaneously simulate 10,000 independent food facilities with the same set of distinct rooms. We considered variability in initial (baseline) contamination assumptions, different contamination transfer rates among food handlers and food as well as non–food‐contact surface areas, and sanitation practices (e.g., wearing gloves, cleaning contact surfaces). Each food facility was assumed to operate 10 hours per day (between 8:00 a.m. and 6:00 p.m.) for 30 days.

Mathematical Descriptions

Spatial Distribution of Microbial Pathogens in Environment

Microbial pathogens may enter the food‐facility environment through raw food ingredients, food handlers, or mobile equipment, such as carts; through leaks and openings throughout the buildings; or through pests, and further spread within the facility environment (Den Aantrekker et al., 2003; Fredriksson‐Ahomaa, Korte, & Korkeala, 2000; Lawrence & Gilmour, 1995; Miettinen, Bjorkroth, & Korkeala, 1999; Norton et al., 2001; Thorberg & Engvall, 2001). To characterize the spatial distribution of microbial pathogens (e.g., L. monocytogenes) within the facility environment, we used the approach discussed in Mokhtari et al. (2018), and broke down the surface area for each object (e.g., floors, walls, and other objects in a room) to 4″ × 4″ grid cells representing potential swab‐sample locations. In each Monte Carlo realization of the model, the initial number of contaminated grid cells on each object (nc) was calculated as: where N is the number of grid cells defined for Object i (e.g., cutting board in the food processing area) and pc is prevalence of contamination, representing fraction of contaminated surface area, calculated for Object i using the environmental sampling data (Miconnet, Cornu, Beaufort, Rosso, & Denis, 2005). Assumptions regarding initial levels of contamination (cfu/+grid cell) were further assigned to the contaminated grid cells for each object (see Table VIII).

Table VIII

List of Model Inputs as Well as Discrete Levels Considered for the Unbalanced Factorial Design Experiments

			Discrete Levels
Model Inputs	Description	Units	Scenario (i)	Scenario (ii)	Scenario (iii)	Scenario (iv)
C0Z1	Initial contamination—individual objects in Zone 1a	log10 cfu/+grid	1,2,3,4,5,6	NA	NA	NA
C0Z2	Initial contamination—individual objects in Zone 2a	log10 cfu/+grid	NA	1,2,3,4,5,6	NA	NA
C0Z3	Initial contamination—individual objects in Zone 3a	log10 cfu/+grid	NA	NA	1,2,3,4,5,6	NA
C0Z4	Initial contamination—individual objects in Zone 4a	log10 cfu/+grid	NA	NA	NA	1,2,3,4,5,6
PSZ1	No. positive samples—individual objects in Zone 1b	Unitless	0,1,2,3,4,5	NA	NA	NA
PSZ2	No. positive samples—individual objects in Zone 2b	Unitless	NA	0,1,2,3,4,5	NA	NA
PSZ3	No. positive samples—individual objects in Zone 3b	Unitless	NA	NA	0,1,2,3,4,5	NA
PSZ4	No. positive samples—individual objects in Zone 4b	Unitless	NA	NA	NA	0,1,2,3,4,5
EffHP	Efficacy of hygiene practices	Unitless	0, 0.25, 0.50, 0.75, 1.00
EffSD,AS	Efficacy of sanitation and disinfection activities for pathogens on accessible surface areas	Unitless	0.50, 0.75, 1.00
EffSD,HC	Efficacy of sanitation and disinfection activities for pathogens in holes and cracks	Unitless	0, 0.25, 0.50
RHW	Pathogen reduction—handwashing	log10	0,1,2,3
RSH	Pathogen reduction—sanitizing shoes	log10	0,1,2,3,4,5
RS	Pathogen reduction—sanitizer	log10	1,2,3
RD	Pathogen reduction—disinfectant	log10	4,5,6
RFHC→AS	Pathogen release fraction during surface cleaning from in holes/cracks to more accessible surface areas	Unitless	10−6,10−5,10−4,10−3
HCHW	Hygiene compliance—handwashing	Unitless	0, 0.25, 0.50, 0.75, 1.00
HCSH	Hygiene compliance—shoes	Unitless	0, 0.25, 0.50, 0.75, 1.00
HCG	Hygiene compliance—gloving	Unitless	0, 0.25, 0.50, 0.75, 1.00
HCGC	Hygiene compliance—glove change	Unitless	0, 0.25, 0.50, 0.75, 1.00
SCSA	Sanitation compliance—surface areas	Unitless	0, 0.25, 0.50, 0.75, 1.00
HFGC	Hygiene frequency—glove change	No. prepared foods	5, 10, 20, 50, 100
SFFP	Sanitation frequency—food processing area	days	0.5,1,2,3,7
SFR	Sanitation frequency—restroom area	days	0.5,1,2,3,7
SFGF	Sanitation frequency—general facility area	days	0.5,1,2,3,7
AT	Ambient temperature	oC	15,20,25,30,35
aw	Water activity—individual objects in different sampling zonesc	Unitless	0.90,0.92,0.94,0.96,0.98,1.00
pH	pH—individual objects in different sampling zonesc	Unitless	4,5,7,9,10

Initial contamination levels (log10 cfu/+grid) were randomly assigned to individual objects in each sampling zone. For example, random contamination levels from the range of possible options (i.e., 1, 2, …, 6) were assigned to utensil, slicer, scale, sink, and food‐contact surface area in Scenario (i).

Number of positive swab samples (out of total number of samples of 10) were randomly assigned to individual objects in each sampling zone. For example, a random number of positive samples from the range of possible options (i.e., 0,1, 2, …, 5) were assigned to utensil, slicer, scale, sink, and food‐contact surface area in Scenario (i).

Each object in the facility environment (e.g., slicer) was assigned with an initial level of water activity as well as pH. These values were further used when adjusting the Listeria monocytogenes growth rate using Equations (14)–(18).

Persistence of Microbial Pathogens in Environment

Microbial pathogens can hide in harborage sites (contamination niches) that are difficult to clean and disinfect. Water and organic soil may stagnate in harborage sites that are difficult to clean and disinfect and difficult to rinse and where some disinfectant may thus still be present after rinsing. As bacteria are able to adapt to low disinfectant concentration, conditions are consequently met for bacterial growth (Carpentier & Cerf, 2011; FDA, 2018; Holah et al., 2012; Motarjemi & Lelieveld, 2014). These cells then can spread to nearby surroundings through a variety of activities, such as cleaning events (e.g., scrubbing the floor of the food processing area) or cross‐contamination events (e.g., transfer from surface to food during slicing). Presence of microbial pathogens in a contamination niche after a cleaning event depends on several factors, including the efficacy of cleaning, the efficiency of the sanitizer or disinfectant used during cleaning, and initial number of bacterial cells prior to the cleaning event (Carpentier & Cerf, 2011). Studies have shown that cleaning efficacy is lower on surviving bacterial cells attached to the surface than on recently attached cells (Marouani‐Gadri et al., 2010; Pan, Breidt, & Kathariou, 2006).

For an object with an ability to harbor microbial pathogens in holes and cracks (e.g., floor of the food processing area, slicer), we assumed that initial pathogen load (cfu) on a grid cell defined on the object's surface area (C0) can be divided into two compartments: (i) bacterial cells on more accessible parts of the grid's surface area (C0 − cfu) and (ii) bacterial cells that are harbored in holes and cracks and are difficult to clean (C0 − cfu), where . We also assumed that fresh cells added to the grid were initially on the more accessible surface areas (i.e., C0). Pathogen removal during the first cleaning event of grid i was modeled as: where R is the contamination reduction due to use of sanitizer or disinfectant during a cleaning event (log10), and Eff is the efficacy of sanitation/disinfection activity (values between 0 and 1) for pathogens on more accessible parts of the grid's surface area. This parameter represents the likelihood of achieving the expected R. A fraction of the residual pathogens on the accessible surface area after a cleaning event () was further added to the pathogen load in holes and cracks of the grid cell: where SA is surface area for holes and cracks in grid i, and SA is grid i surface area. A subsequent cleaning event not only reduces pathogen loads on the more accessible surface area of the grid (Equation (3)), but also can further reduce pathogen loads in the harborage sites, although with a lower cleaning efficacy (Marouani‐Gadri et al., 2010; Pan, Breidt, & Kathariou, 2006):where C0 is pathogen load in harborage sites on grid i, and Eff is the efficacy of sanitation/disinfection activity (values between 0 and 1) for pathogens in the harborage sites on grid i. We also assumed that a fraction of pathogen loads in the harborage site (RF) was released during the subsequent cleaning event (e.g., due to scrubbing, washing with low‐pressure hoses or high‐pressure jets) and added to the surface area that was more readily available for spreading to other objects in contact with the grid's surface area (e.g., subsequent contacts with food handlers’ hands or food servings). Therefore, contamination levels (cfu) on the more accessible surface area and in the holes and cracks of the grid i were updated as:

During a cleaning event, reduction of pathogens on accessible surface areas and in holes and cracks of contaminated objects could occur simultaneously with release of pathogens from harborage sites of contaminated objects to more accessible surface areas. We simplified the modeling approach by assuming that these two events occurred in order, as described in Equations (4)–(6) and Equations (7)–(10).

Growth of Microbial Pathogens in Contamination Niches

We assumed that, between two consecutive cleaning events, pathogen‐harborage sites on a grid cell would have optimal growth conditions; hence, we could anticipate increase in pathogen loads trapped in those areas. We used the following equations to calculate the microbial pathogen growth in contamination niches (Koseki & Isobe, 2005): where X(t) represents the log10 of the pathogen loads in holes and cracks (log10 C), q(t) is a dimensionless quantity related to the physiological state of the pathogens; μ is the specific growth rate; X represents maximum population density (MPD), calculated based on the maximum bacterial concentration on a contact surface (BC − log10 cfu/cm2) and the surface area for holes and cracks in grid cell i (SA); and X is the initial population size (log10 C). The lag time variable q evolves according to exponential growth, and q determines the duration of the lag time in microbial growth. We further described the effect of ambient temperature (T), pH, and water activity (aw) on μ of a microbial population using the cardinal secondary growth model as below (Rosso, Lobry, Bajard, & Flandrois, 1995): where is the optimal growth rate for pathogens and X, X, and X are minimum, maximum, and optimal values of . We substituted the model for into the coupled differential equations, then numerically solved the system of equations by the fourth‐order Runge‐Kutta method, to predict pathogen growth in holes and cracks of grid i during the time between two consecutive cleaning events. Initial values for specific growth parameters (i.e., aw, pH) were randomly selected from the range of possible values (see Table VIII) and assigned to individual objects that contained holes and cracks at time = 0. To simplify the modeling approach, we did not include changes in selected growth parameters during facility operation.

Inactivation of Microbial Pathogens on Contact Surfaces

The number of viable microbial pathogens generally decreases over time when the microorganism is in an environment that does not support growth. We assumed that more accessible parts of each contact surface did not support growth of pathogens and that, hence, pathogens on such areas may experience decline over time. However, when growth conditions were not met in pathogen‐harborage sites on a grid (e.g., , pathogens could decline over time. The rate of inactivation is dependent on a variety of factors, including pH, acidulant identity, acidulant concentration, water activity, concentration of antimicrobials, and temperature (Ahmad & Marth, 1989; Buchanan & Golden, 1995; Buchanan, Golden, & Philips, 1997; Buchanan, Golden, & Whiting, 1993; Buchanan, Golden, & Whiting, 1994; Cole, Jones, & Holyoak, 1990; El‐Shenawy & Marth, 1989; Parish & Higgins, 1989; Sorrells, Enigl, & Hatfield, 1989). Levels of microbial pathogens on more accessible parts of each contact surface and in contamination‐harborage sites (when growth conditions were not met) were updated as: where LR is the rate of inactivation (log10 cfu/hour) and Δt is the timestep (hours) between two consecutive events during which pathogen levels were updated (e.g., Δt = 1 hour).

Microbial Cross‐Contamination During Tactile Contact Events

We modified the approach discussed in Hoelzer et al. (2012) for modeling cross‐contamination between any two objects (e.g., A and B) as follows, to enable us to consider (i) the impact of spatial distribution of microbial contamination on objects involved in a tactile contact event and (ii) pathogen loads on surface areas that are more readily available for transferring to other objects during tactile contacts (CAS): where n and n are the numbers of grid cells on Object A and Object B, respectively, involved in the contact event; C0 is initial microbial contamination (cfu) on the readily accessible areas of grid i on Object A; and C0 is initial microbial contamination (cfu) on the readily accessible areas of grid j on Object B. TR and TR are contamination transfer rates (0 ≤ TR ≤ 1) from Object A to B and from Object B to A, respectively. Updated microbial contamination values on Objects A and B after a contact event (C and C, respectively) were randomly distributed among the more readily accessible surface areas of the selected grids involved in the contact event.

When an agent (i.e., food handler) was involved in a contact event, the model considered contamination level on the agent's hands (or gloves) when the receiving object in Equations (20) and (21) was another agent or an object other than the room floor. When the receiving object in these equations was the room floor, the model considered contamination level on the agent's shoes.

Dynamic Schedule of Events

Successful implementation of an ABM framework relies on proper scheduling of different events and real‐time activities (e.g., scheduling tactile contacts between a food handler and the slicer during a food‐preparation activity). Most food facilities operate in dynamic environments where unpredictable real‐time events may cause a change in prescheduled activity plans, and a previously feasible schedule may turn infeasible (Ouelhadj & Petrovic, 2008). We used the predictive‐reactive scheduling (PRS) approach to identify and prioritize a list of real‐time events and activities that may occur during the operation of a food facility. The PRS approach is the dynamic scheduling technique most commonly used in modeling complex systems, such as manufacturing operations (Aytug, Lawley, McKay, Mohan, & Uzsoy, 2005; Herroelen & Leus, 2005; Mehta & Uzsoy, 1999; Vieira, Herrman, & Lin, 2000; Vieira, Herrman, & Lin, 2003).

F2‐QMRA generated a daily dynamic schedule that included a list of all events and activities that could occur during a single day in a facility. In a dynamic schedule, an event may represent participation in a particular daily activity (e.g., morning food preparation, personal hygiene, or surface sanitation) or the occurrence of a tactile contact between two agents (e.g., contacts between food handlers in the loading dock), between an agent and an object (e.g., contact between a food handler and the slicer in the food processing area), or between two objects (e.g., contacts between a food serving and cutting board in the food processing area). Each event was defined by its occurrence time (e.g., 9:36 a.m.), its location (e.g., food processing area), and list of agents and/or objects involved in the event (e.g., a food handler and a slicer). The model used a predefined daily list of food handlers and facility‐operation activities to create the dynamic schedule, including all contact events. Table III shows examples of predefined daily activities, including activity priority, activity location, start and end times for each daily activity, and number of food handlers involved in each activity. Activity priorities were used for individual agents to select among the list of activities that could occur at the same time (e.g., morning idle time in the office area vs. unloading food ingredients in the loading dock), with higher likelihood of participation assigned to activities with higher priorities. For concurrent activities with similar priorities, an agent randomly selected among the list of available options (e.g., unloading food ingredients in the loading dock vs. sanitizing the food processing area).

Table III

Example List of Predefined Activities in F2‐QMRA

Activity Names	Activity Priority	Activity Location	Start Time	End Time	No. Agents Involved in the Activity
Morning arrivala	10	Office	7:00 a.m.	8:00 a.m.	All
Idle time in the morning	5	Office	8:00 a.m.	9:00 a.m.	All
Unloading food ingredients	10	Loading dock	8:30 a.m.	9:00 a.m.	2
Food processing area sanitation	10	Food processing area	8:30 a.m.	9:00 a.m.	2
Morning food preparation	5	Food processing area	9:00 a.m.	12:00 p.m.	All
Lunch and break	5	Office	12:00 p.m.	1:00 p.m.	All
Afternoon food preparation	5	Food processing area	1:00 p.m.	5:00 p.m.	All
Loading prepared foods to food trucks	10	Loading dock	4:00 p.m.	5:00 p.m.	2
Facility sanitation	10	All rooms	4:00 p.m.	5:00 p.m.	2
Idle time in the afternoon	5	Office	4:00 p.m.	6:00 p.m.	All
Restroom visitb	10	Restroom	8:00 a.m.	6:00 p.m.	All
Departurec	10	NA	5:00 p.m.	6:00 p.m.	All

Each food handler arrived at a random time between 7:00 a.m. and 8:00 a.m. in the morning.

Start times for restroom visit activities were randomly selected for food handlers between 8:00 a.m. and 6:00 p.m. Duration of the restroom visit activity was randomly selected for each event.

Each food handler departed at a random time between 5:00 p.m. and 6:00 p.m. in the afternoon.

Fig. 2 shows a simplified flow diagram illustrating how a dynamic schedule was created in F2‐QMRA. Each food handler arrived in the morning, at an arrival time randomly selected between 7:00 a.m. and 8:00 a.m. (bold box in Fig. 2). The Arrival Event for Agent i was added to the dynamic schedule. The model also generated Restroom Visit Events with random occurrence times between 8:00 a.m. and 6:00 p.m. for Agent i and added them to the schedule. When all agents arrived at the facility and corresponding arrival and restroom visit events were recorded, the model sorted the schedule and selected the first event that had the earliest occurrence time (e.g., Arrival Event for Food Handler #3 at 7:05 a.m.). If the selected event corresponded to initiation of a particular activity (e.g., arrival), the model added a new tactile contact event between the agent involved and the doorknob assigned to the activity location representing the agent entering the specific activity location (Table III). The model also went through the list of activity‐specific behavioral rules (Table IV) and added corresponding events to the schedule, when applicable (e.g., prior to entering the food processing area, each food handler might clean his/her shoes).

Figure 2

Dynamic scheduling component of F2‐QMRA.

Table IV

List of Activity‐Specific Behavioral Rules and Their Mathematical Formulations in F2‐QMRA

Activity	Applicable Behavioral Rules	Mathematical Formulations	Model Parameters
Food preparation	Food handler might sanitize his/her shoes prior to entering the food processing area	C_SH,i=C0_SH,i×(1−α)+CSH∗×α α∼Binomial(1,HC_SH) CSH∗∼Binomial(C0_SH,i,110^(R_SH×Eff_HP))	CSH,i: Contamination on shoes (cfu) for agent i after sanitation C0SH,i: Initial contamination on shoes (cfu) for agent i RSH: Pathogen reduction‐ sanitizing shoes (log10) EffHP: Efficacy of hygiene practices HCSH: Hygiene compliance‐shoes
	The food handler might wash his/her hands prior to preparing any food	C_H,i=C0_H,i×(1−α)+CH∗×α α∼Binomial(1,HC_HW) CH∗∼Binomial(C0_H,i,110^(R_HW×Eff_HP))	CH,i: Contamination on hands (cfu) for agent i after sanitation C0H,i: Initial contamination on hands (cfu) for agent i RHW: Pathogen reduction‐ handwashing (log10) EffHP: Efficacy of hygiene practices (values between 0 and 1) HCHW: Hygiene compliance‐handwashing
	The food handler might wear gloves during food preparation	{Binomial(1,HC_G)=0; agent i does not wear gloves Binomial(1,HC_G)=1; agent i weara gloves	HCG: Hygiene compliance‐gloving
	Food handler might change gloves after preparing certain number of food servings	{Binomial(1,HC_GC)=0; agent i does not wear gloves Binomial(1,HC_GC)=1; agent i weara gloves	HCGC: Hygiene compliance‐glove change
	Each food preparation event included random contacts with utensil, food‐contact surface area (cooking table), and food serving	See Section 2.2.5‐Equations (20) and (21)	NA
Unloading ingredients	If the food handler contacted a case (e.g., food ingredient package), a subsequent contact between the case and the cart was scheduled; the case was further moved to the food processing area	See Section 2.2.5‐Equations (20) and (21)	NA
Restroom visit	The food handler removed his/her gloves before visiting restroom	C_G,i=0	CG,i: Contamination on gloves for agent i
	The food handler made a contact with toilet seat	See Section 2.2.5‐Equations (20) and (21)	NA
	The food handler made a contact with the toilet flush	See Section 2.2.5‐Equations (20) and (21)	NA
	The food handler made a contact with handwashing sink	See Section 2.2.5‐Equations (20) and (21)	NA
	The food handler might wash his/her hands before leaving the restroom	See equations under Food Preparation/Handwashing activity	NA

C: Contamination on shoes (cfu) for agent i after sanitation

C0: Initial contamination on shoes (cfu) for agent i

R: Pathogen reduction‐ sanitizing shoes (log10)

Eff: Efficacy of hygiene practices

HC: Hygiene compliance‐shoes

C: Contamination on hands (cfu) for agent i after sanitation

C0: Initial contamination on hands (cfu) for agent i

R: Pathogen reduction‐ handwashing (log10)

Eff: Efficacy of hygiene practices (values between 0 and 1)

HC: Hygiene compliance‐handwashing

Each contact event potentially resulting in cross‐contamination was further executed by updating contamination levels on agents and/or objects involved in the event (Equations (20) and (21)). After executing a contact event, the model removed the event from the dynamic schedule and generated a list of new events representing potential future contact events among the agent involved in the original activity and other agents present in the room and/or objects assigned to the activity location. Times associated with these new contact events were generated using a Poisson process (Lawler, 2006): where t0 represents the time associated with the original contact event (e.g., time associated with contact with the office doorknob), α represents a random number generated between 0 and 1, and CR is the activity‐specific contact rate among the agent initially involved in the contact event and other agents or objects involved in the contact event (e.g., contact rate between a food handler and the slicer during the food‐preparation activity). F2‐QMRA further sorted the list of potential contact events and selected the one that occurred in the nearest future. If time associated with the selected contact event was within the ongoing activity timeline (e.g., between 9:00 a.m. to 12:00 p.m. during the morning food preparation activity), the new event was added to the dynamic schedule; otherwise, the agent was forced to leave the area at the end of the current activity time and select a new activity based on those listed in Table III. The model continued generating new events, executing current events in the dynamic schedule, updating the dynamic schedule, and moving agents between different rooms in the facility until the end of the operation time (e.g., day = 30; time = 6:00 p.m.).

Summary of the Model Inputs

F2‐QMRA includes inputs that require quantitative estimates. The data we used to populate selected model inputs related to the facility setup and tactile‐contact rates during facility operation are summarized in Tables V and VI, respectively. Transfer coefficients, used in Equations (20) and (21), varied from one transfer to another according to lognormal distributions. We used specific source‐recipient transfer‐rate data (e.g., transfer rate between stainless steel representing the slicer blade and the deli meat representing the food ingredient) generated in a systematic literature review conducted by Hoelzer et al. in this model (Hoelzer et al., 2012). Minimum, maximum, and optimal values of T, pH, aw, and for growth of L. monocytogenes (Equations (14), (15), (16), (17), (18)) are listed in Table VII (Augustin & Carlier, 2000; Augustin, Zuliani, Cornu, & Guillier, 2005; dos Reis‐Teixeira, Alves, & de Martinis, 2017; Le Marc et al., 2002). In the absence of surface‐specific data for modeling survival of L. monocytogenes on the more accessible parts of each contact surface and in contamination‐harborage sites, we used data generated for the survival of L. monocytogenes on a conveyor belt's material with or without antimicrobial additives, and at temperatures of 10 °C, 25 °C, and 37 °C (Chaitiemwong, Hazeleger, & Beumer, 2010). Data suggested large decreases in numbers of pathogens during the first 6 hours and lower decrease rates between 6 and 72 hours. Due to the dynamic nature of the facility environment, and to simplify the modeling approach, we calculated the total decline in the numbers of L. monocytogenes during 72 hours and estimated the hourly rate of reduction (LR‐log10 cfu/hour) as a function of ambient temperature (T—°C): For other model inputs, in which values were not specifically governed by the microbiology, chemistry, or physics of the situation (e.g., initial levels of contamination on different objects in rooms, surface‐sanitation frequency, surface‐sanitation efficiency, among others), we used discrete distributions representing values with equal probabilities spanning across the range of possible options, rather than choosing a specific set of baseline values. Using this approach, we were able to efficiently generate data on a wide range of what‐if scenarios during a single Monte Carlo simulation.

Table V

Room‐Specific Objects Defined for the Case Study Scenarios

Room‐Specific Objects	Zone id	No. Grids Defined for Each Objecta	No. Specific Object in Selected Roomb	Fraction of the Object Surface Area Representing Holes and Cracksb
Food processing area
Utensil	1	1	10	0
Slicer	1	30	3	1E‐1
Scale	1	40	1	0
Sink	1	200	2	0
Fridge	2	800	1	1E‐2
Trash can	3	100	1	0
Food‐contact surface	1	250	4	1E‐3
Non–food‐contact surface	2	50	10	1E‐2
Floor	3	4,000	1	1E‐3
Walls	3	8,000	4	0
Doorknob	3	1	1	0
Common area (office)
Chair	4	40	4	0
Desk	4	200	2	0
Floor	4	1,500	1	1E‐2
Walls	4	3,000	4	0
Doorknob	4	1	1	0
Restroom
Toilet	4	40	1	0
Sink	4	100	1	0
Flush	4	1	1	0
Floor	4	750	1	1E‐2
Walls	4	1,500	4	0
Doorknob	4	1	1	0
Loading dock
Cart	4	100	2	0
Ingredient package	3	50	20	0
Floor	4	3,000	1	1E‐1
Walls	4	6,000	3	0
Doorknob	4	1	1	0

Value of 1 indicates no spatial resolution was considered for the object.

Model assumptions—zero values represent no capacity to harbor microbial contamination.

Table VI

List of Random Tactile Contacts and Their Rates During Different Activities

	Number of Random Tactile Contactsa (Min, Most Likely, Max)
Activity‐Specific Objects	Food Handler	Ingredient/Food Packages
Food preparation (food processing area)
Utensil	ESb
Slicer	(5,8,10)
Scale	(2,3,5)
Sink	(2,3,5)
Fridge	(3,4,5)
Trash can	(2,3,5)	(0,1,1)d
Ingredient/food packages	(2,3,5)
Food‐contact surface	ESb
Non–food‐contact surface	(2,3,5)	(0,1,2)
Floor	(10,20,30)	(2,3,5)
Walls	(0,1,1)
Doorknob	EEc
Other food handlers	(1,2,3)
Arrival and idling (office)
Chair	(2,3,5)
Desk	(2,3,5)
Floor	(5,10,15)
Walls	(0,1,1)
Doorknob	EE c
Other food handlers	(3,4,5)
Restroom visit (restroom)
Sink	ESb
Toilet	ESb
Flush	ESb
Floor	(3,5,8)
Walls	(0,0,1)
Doorknob	EEc
Loading/unloading ingredients/prepared foods (loading dock)
Ingredient/food packages	(5,8,10)
Cart	(2,3,5)	(1,2,3)
Floor	(15,20,25)	(3,4,5)
Walls	(0,1,1)
Doorknob	EEc
Other food handlers	(2,3,5)

Numbers of random tactile contacts during different activities are basic model assumptions.

Event specific: (i) for preparing a food serving, series of random tactile contacts are created between the food handler, food, food‐contact surface, and utensil; and (ii) each restroom visit includes individual contacts between food handler and toilet, flush, and sink (in order).

Entry/exit: one doorknob contact is initiated upon entry and exit from each room.

Random contacts between food/ingredient packages and trash can represent accidental events based on the observed inspection data (data not shown here).

Table VII

Minimum, Maximum, and Optimal Growth Parameters for Listeria monocytogenes

Growth Parameters	Units	Minimum	Maximum	Optimal Value	Reference
Temperature (T)	°C	4.26	45.50	37	Le Marc et al. (2002)
pH	Unitless	4.71	9.61	7.10	Augustin and Carlier (2000); Augustin et al. (2005)
Water activity (aw)	Unitless	0.913	1.00	0.997	Augustin and Carlier (2000)
Growth rate (μ)	log10 cfu/hour	‐	‐	[0.57, 1.32]a	Augustin and Carlier (2000); Augustin et al. (2005)
Maximum bacterial concentration on a surface (BCmax)	Log10 cfu/cm2	‐	‐	[6, 8]b	dos Reis‐Teixeira et al. (2017)

Optimal growth rate for L. monocytogenes was modeled as a range of values based on the available data on wide variety of foods, including seafoods, dairy, and vegetables.

Maximum bacterial concentration on a contact surface was modeled as a range of values based on the available data.

Probabilistic Analysis Scenarios

F2‐QMRA includes inputs that represent variability or uncertainty in available data. For example, during a food‐preparation event, the number of tactile contacts between employee's hands (or gloves) with specific objects in the food processing area (e.g., food‐contact surface, slicer, and refrigerator, among others) was simulated using probability distributions representing variability in such events. Meanwhile, the prevalence of contamination or fraction of each contact‐surface area (e.g., floor of the food processing area) that was contaminated was represented with an uncertainty distribution based on the number of positive samples collected from that surface area (Equation (2)). Alternative approaches can be used to distinguish between variability and uncertainty in inputs and their contribution to the output variation (Cullen & Frey, 1999; Frey, Mokhtari, & Zheng, 2004). The choice of an appropriate approach should be made taking into account the assessment objectives, data quality objectives, and availability of data.

To demonstrate the model's capabilities in quantifying the spread of microbial contamination in a food facility, we chose a one‐dimensional probabilistic simulation. In this approach, we incorporated, but did not distinguish, variability and uncertainty in inputs, thereby exploring the impact of the possible/potential range of values for each model input (Mokhtari & Frey, 2005). We further developed a factorial design in which selected, predefined values were randomly assigned to different inputs and propagated through the model, using a Monte Carlo simulation. Each Monte Carlo simulation included 10,000 model realizations, or runs, with each realization representing a “treatment” in the factorial design, for which specific combinations of values were randomly assigned to different model inputs. For example, one realization could represent a factorial design treatment in which initial level of contamination on the slicer was assumed to be 1 − log10 cfu/+grid; number of positive samples taken from the slicer (out of 10 collected samples) was assumed to be 1; and the slicer was assumed to be cleaned every day, among other assumptions. A fully balanced factorial design would include all combinations of values assigned to all model inputs and was deemed unrealistic; instead, we created an unbalanced factorial design that included partial combinations of values assigned to selected model inputs, using our simulation sample size of 10,000 model realizations to reduce the required computational time and resources.

We further created four overarching scenarios involving initial levels of contamination on the surface areas of the objects located in each individual sampling zone (e.g., Scenario (i) represented initial contamination on objects assigned to Zone 1, while all other objects received no initial contamination). For each scenario, we generated four outputs across 30 days of model simulation: (i) daily pathogen prevalence in prepared foods (percent contaminated food servings with at least 1 cfu/25 g of sample as the level of detection for L. monocytogenes), (ii) daily average levels of L. monocytogenes in prepared foods that were positive (log10 cfu/+serving), (iii) daily average prevalence for objects assigned to different sampling zones, and (iv) daily average levels of L. monocytogenes on positive objects assigned to different sampling zones (log10 cfu). For the latter two outputs, contamination prevalence and level of individual objects assigned to different sampling zones were checked every hour during facility operation and average daily values were recorded at the end of each operation day. An object in a sampling zone (e.g., food‐contact surface in Zone 1) was considered positive for L. monocytogenes if one or more grid cells on the object was contaminated with the bacterium. By aggregating results across individual objects assigned to each sampling zone, average daily values of contamination prevalence and level were also calculated for each sampling zone. Full descriptions of selected model inputs and the discrete values assigned to the inputs for the four overarching scenarios are provided in Table VIII.

Global Model Sensitivity Analysis

Sensitivity analyses offer a vital tool for quantifying the effects of input parameters on model predictions, thereby providing insight into which of the model inputs may be most influential, and, hence, answer questions related to the most effective mitigation options for reducing risks (Cullen & Frey, 1999; Mokhtari & Frey, 2005; Saltelli, Chan, & Scott, 2000). The utility of sensitivity analysis can be maximized when the analysis considers the full range of values associated with model inputs by varying multiple inputs simultaneously. This approach, referred to as global sensitivity analysis, acknowledges potential interactions and nonadditive effects (Saltelli, Tarantola, & Chan, 1999). However, global sensitivity analysis has not been widely used because the probabilistic model outputs can be unwieldy, and methods of analyzing these data can be computationally intensive (Cariboni, Gatelli, Liska, & Saltelli, 2007; Fieberg & Jenkins, 2005; Naujokaitis‐Lewis, Curtis, Arcese, & Rosenfeld, 2009; Saltelli et al., 1999). To address these limitations, we used random forest as a global sensitivity analysis method to identify the most influential model inputs. Random forest is a nonparametric classification method that applies random subsets of the data to produce thousands of classification and regression trees (CART), which then are used to calculate classification error rates and to quantify model‐input importance values as a measure of sensitivity (Breiman, 2001; Gromping, 2009; Jones & Linder, 2015). Random forest has demonstrated several key advantages for sensitivity analysis of complex models (Harper, Stella, & Fremier, 2011): (i) it is accessible and easy to use because multiple packages are available in R and can be downloaded for free; (ii) sensitivity measures for model inputs include the effects of all higher‐order interactions and do not assume linearity; and (iii) generating robust global sensitivity results for models with a fairly large number of inputs requires a sample of hundreds of the model's Monte Carlo realizations, rather than a fully factorial set of several million realizations (Wagner, 1995).

The impact of model inputs on L. monocytogenes contamination prevalence (%) and levels (log10 cfu/+serving) in prepared foods was evaluated using permutation importance scores. The permutation importance scores in random forest reflect the decrease in prediction accuracy (represented by the % increase in the mean squared error [MSE]) resulting from permutation of a model input. If permuting the value of a model input does not affect (or increase) the prediction accuracy, then the input is not related to the model output. However, if permuting the value decreases the prediction accuracy of the model, the model input is related to the model output, and the larger the increase in MSE, the stronger the relationship (Genuer, Poggi, & Tuleau‐Malot, 2010; Jones & Linder, 2015).

RESULTS

Model Validation

Given the number of inputs in our model, a formal calibration (e.g., minimizing some objective function) or a validation currently is not possible. Nevertheless, we performed checks and controls to ensure the validity of the model's calculations. Specifically, the sum of the L. monocytogenes contamination levels in the facility environment, at any point in time, was equal to the sum of: (i) L. monocytogenes contamination introduced from outside sources (initial contamination on different objects assigned to one of the four environmental sampling zones); (ii) growth within the facility environment, including all contamination niches; and (iii) removal from the facility via exiting of prepared food and through sanitation, disinfection, and natural die‐off on contact‐surface areas. Verification of this mass balance additionally afforded a cross‐check of proper functioning of this model. The mass balance was controlled in all scenarios described in this study (data not shown).

Spread of L. monocytogenes Contamination in the Food‐Facility Environment

We investigated the spread of L. monocytogenes from objects in one sampling zone to objects in other zones, primarily arising from contacts with employees’ hands, gloves, or shoes, as well as movement of contaminated objects from one sampling zone to another (e.g., movement of contaminated food ingredient packages from the loading dock to the food processing area and subsequent contacts with non–food‐contact surfaces, floor, and trash can). Fig. 3 shows average contamination prevalence (%) on objects placed in different sampling zones on selected facility‐operation days when initial contamination on Day 1 started on objects in Zone 1 (Fig. 3(a)), Zone 2 (Fig. 3(b)), Zone 3 (Fig. 3(c)), and Zone 4 (Fig. 3(d)) and no additional contamination events took place after the Day 1 event. The whiskers in Fig. 3 represent the 1st and the 3rd quartiles and the median output values (horizontal bold line in each whisker); the ends of the whiskers (the vertical bars) represent the 2.5th and 97.5th percentile values. The isolated dots represent the outlier results. An outlier was any value that lay more than one and a half times the length of the whisker box from either side. Results showed that once contamination was introduced in a particular sampling zone (e.g., objects assigned to Zone 1), it reached objects placed in all other sampling zones; however, prevalence of contamination eventually reduced to values close to zero in most of the Monte Carlo realizations of the model (e.g., median average prevalence values in different sampling zones were close to zero), primarily due to sanitation/disinfection activities, pathogen die‐off, and removal from facility environment via exit of contaminated foods (e.g., enforcing pathogen mass balance in the facility environment). This finding is crucial, especially for scenarios in which initial contamination was added to objects that did not have direct contact with food servings and were not in close proximity to the food‐contact surface areas in the food processing environment. For example, once initial contamination was added to objects assigned to Zone 4 on Day 1 (e.g., objects in the loading dock and restroom environments), L. monocytogenes cells were transferred to objects assigned to Zone 1, with median contamination prevalence values of approximately 24%, 13%, and 2% on Days 1, 2, and 3, respectively. For the same scenario, median contamination prevalence for objects in Zone 1 was approximately 0% on Day 7 and beyond, although certain Monte Carlo realizations of the model resulted in prevalence values as high as 10% in Zone 1, even on Day 30 (Fig. 3(d)). Summary results also showed that L. monocytogenes persisted on objects for a longer time in the facility environment once the bacterium was initially introduced to Zone 4 (Fig. 3(d)) compared with other scenarios in which initial contamination was added to objects in Zones 1, 2, and 3 (Figs. 3(a)–3(c), respectively). For example, once contamination was initially introduced to Zone 4, the interquartile of the average daily contamination prevalence for objects placed in this zone ranged between 0% and 6% with a median value of 0.7% on Day 21 (Fig. 3(d)), while for the other three scenarios shown in Fig. 3, interquartile ranges of the average daily contamination prevalence for objects placed in those zones reduced to values very close to 0% within the first two weeks of facility operation. This was primarily due to abundance of contamination niches on objects assigned to Zone 4, represented by the fraction of holes and cracks on objects assigned to this zone as well as their surface areas (see Table V), and possibility of L. monocytogenes harboring in holes and cracks during routine sanitation and disinfection activities. For the same scenario, we also evaluated the impact of introducing additional contamination events after the initial Day 1 contamination event. This scenario resulted in persistence of L. monocytogenes in different sampling zones beyond the 30 days included in the baseline scenario with only one contamination event on Day 1 (data not shown).

Figure 3

Spread of L. monocytogenes contamination within food‐establishment environment. Average daily contamination prevalence (%) in different sampling zones is shown for selected operation days when initial contamination started from: (a) Zone 1, (b) Zone 2, (c) Zone 3, and (d) Zone 4.

We further investigated this scenario by identifying the contamination spots in the food processing area once L. monocytogenes was initially added to objects assigned to Zone 4 on Day 1. Figs. 4 and 5 show average L. monocytogenes contamination prevalence (%) and levels (log10 cfu/+object), respectively, for individual objects assigned to the food processing area on selected operation days when initial contamination started on Day 1 from Zone 4. Results showed that main modes of L. monocytogenes spread from Zone 4 to the food processing area included (i) movement of contaminated objects from the loading dock to the food processing area (i.e., contaminated food ingredient packages) and subsequent contacts with individual objects in the food processing area, such as the trash can, floor, and non–food‐contact surface areas; (ii) employees’ hands/gloves, in the case of the doorknob and equipment used during food‐preparation activities, such as food‐contact surfaces, slicer, and scale; and (iii) employees’ shoes, in the case of the floor. These transfer modes resulted in median prevalence values of approximately 100%, 75%, 74%, and 62% for trash can, floor, doorknob, and non–food‐contact surface areas, respectively, on Day 1 (Fig. 4(a)). Contamination prevalence values for objects in the food processing area declined in time, with approximately 0% median prevalence values on Day 7 and beyond for different objects (when no additional contamination events took place), although certain Monte Carlo realizations resulted in prevalence values as high as approximately 100%, 21%, 23%, and 37% for the floor, non–food‐contact surface areas, food‐contact surface areas, and the slicer, respectively, on Day 30 (Fig. 4(d)). Similarly, results showed that median contamination levels (log10 cfu) for contaminated objects in the food processing area declined with time, although for certain Monte Carlo realizations, the model predicted that some contaminated objects (e.g., slicers) had median contamination levels as high as 3‐log10 cfu after 30 days (Fig. 5(d)).

Figure 4

Spread of L. monocytogenes contamination from contaminated objects in Zone 4 to the food processing area. Average daily contamination prevalence (%) for objects in the food processing area is shown for selected operation days.

Figure 5

Spread of L. monocytogenes contamination from contaminated objects in Zone 4 to the food processing area. Average daily contamination levels for positive objects (log10 cfu) in the food processing area are shown for selected operation days.

Spread of L. monocytogenes Contamination from Objects in Different Sampling Zones to Prepared Foods

F2‐QMRA also predicted that prepared food servings could be contaminated with L. monocytogenes initiated from different sampling zones in the facility, although likelihood and levels of contamination in foods were substantially higher when initial contamination started from objects assigned to Zone 1 with direct contacts with food servings (e.g., food‐contact surfaces and slicer, among others). Fig. 6 shows the average daily contamination prevalence and level (log10 cfu/+serving) on selected facility‐operation days once microbial contamination had initiated from different sampling zones on Day 1. Results showed that once objects in Zone 1 were contaminated with L. monocytogenes, median value of contamination prevalence in prepared foods on Day 1 was approximately 50%, although median contamination prevalence values for prepared foods dropped to 9% and 1% on Day 2 and Day 3, respectively, and 0% on Day 7 and beyond. For the same scenario and for certain Monte Carlo realizations of the model, prevalence of L. monocytogenes contamination in prepared foods was approximately 40% on Day 30 (Fig. 6(a)). When L. monocytogenes contamination initiated from objects with no direct contacts with food servings (i.e., objects in Zones 2 to 4), contamination prevalence in prepared food servings was substantially lower, with median value of approximately 0% on Day 3 of facility operation and beyond (Fig. 6(a)).

Figure 6

Spread of L. monocytogenes contamination from different sampling zones with initial contamination to the prepared food products: (a) average daily contamination prevalence in prepared foods (%) and (b) average daily contamination levels in positive prepared foods (log10 cfu/+serving).

Similarly, Fig. 6(b) shows that when contamination initiated from objects with direct contact with food products (i.e., objects in Zone 1), average contamination levels in prepared foods that were positive (log10 cfu/+serving) were substantially higher, with a wider range of variability, compared with scenarios in which microbial contamination initiated from objects with no direct contacts with food servings (i.e., objects in Zones 2 to 4). Results also showed that levels of L. monocytogenes in positive foods generally declined during facility operation when no additional contamination events took place (after initial contamination on objects at the beginning of Day 1). For example, while contaminated objects in Zone 1 (e.g., slicer, food‐contact surface, among others) resulted in a median L. monocytogenes level of approximately 2‐log10 cfu in positive foods on Day 1, median values declined to approximately 0‐log10 cfu/+serving after 30 days.

Sensitivity Analysis

We further evaluated the impact of model inputs on prevalence and levels of L. monocytogenes in prepared foods during the 30 days of facility operation using the random forest method. Sensitivity analysis results showed that, in addition to the initial conditions of the food‐contact surfaces (i.e., both the initial contamination level and the number of positive swab samples), pathogen detachment factor (i.e., fraction of L. monocytogenes released from holes and cracks of different contact‐surface areas to more accessible surface areas during cleaning events), hygiene practices (e.g., sanitation/disinfection compliance, cleaning efficacy and efficiency, and sanitation frequency, for the food preparation area), and parameters supporting growth of L. monocytogenes on slicer (pH and aw) were among the top‐ranked inputs impacting the average prevalence of L. monocytogenes contamination in prepared foods (Fig. 7(a)). Average L. monocytogenes contamination levels in positive food products (log10 cfu/+serving) were primarily influenced by the initial contamination conditions (i.e., both initial L. monocytogenes levels and the number of positive swab samples) associated with objects with direct contact with food products (e.g., food‐contact surface, slicer, and utensil), although sanitation compliance, efficacy, and efficiency also influenced the L. monocytogenes levels in positive foods (Fig. 7(b)).

Figure 7

Top‐ranked model inputs impacting: (a) L. monocytogenes contamination prevalence (%) in prepared foods and (b) L. monocytogenes contamination levels in positive prepared foods (log10 cfu/+serving).

Mutual Impact of Selected Model Inputs on Prevalence and Levels of L. monocytogenes in Prepared Foods

We also investigated the interaction effects between the top two model inputs: sanitation/disinfection compliance and initial contamination on food‐contact surfaces (log10 cfu/+grid). Since food‐contact surfaces were placed in Zone 1, we focused on results obtained from Scenario (i), as described in Section 2.4, for which initial contamination levels were assigned to all objects in this zone (including food contact surfaces), while other objects in Zones 2–4 received no initial contamination. Figs. 8(a) and (b) show the mutual impacts of these two model inputs on average prevalence (%) and average levels (log10 cfu/+serving), respectively, of L. monocytogenes in prepared foods, on selected days of facility operations. Results showed that while compliance with sanitation/disinfection activities can reduce prevalence and levels of L. monocytogenes in prepared foods, full compliance can only lead to zero or very small average prevalence values when food‐contact‐surface areas were not initially contaminated with high levels of the bacterium (i.e., ≤2 log10 cfu/+grid). For example, Fig. 8(a) shows that the most likely value for average contamination prevalence in prepared foods was approximately 0% when facility employees fully complied with sanitation and disinfection activities during facility operation, although Monte Carlo realizations of the model resulted in average prevalence values as high as 12.5% in foods after 30 days. However, when food‐contact‐surface areas were highly contaminated (e.g., 6‐log10 cfu/+grid) on Day 1, full sanitation/disinfection compliance during facility operation timeline could only reduce the average prevalence of L. monocytogenes in prepared foods to 2.5%, with 50% of positive foods having average contamination levels as high as 0.3‐log10 cfu after 30 days (Fig. 8(b)).

Figure 8

Mutual impact of initial contamination on food‐contact surface areas (C0.FCS) and sanitation/disinfection compliance during selected facility operation days on: (a) average daily contamination prevalence in prepared foods (%) and (b) average daily contamination levels in positive foods (log10 cfu/+serving).

DISCUSSION

A limited number of studies have looked comprehensively at modeling persistence and spread of microbial pathogens in food facilities. For example, Schaffner (2004) proposed an approach to evaluating the possibility of modeling cross‐contamination of Listeria species, total L. monocytogenes, or specific L. monocytogenes strains using a probabilistic quantitative mathematical model. Schaffner's approach provided a starting framework for predictive modelers and scientists studying L. monocytogenes to begin research, together with an ultimate goal of understanding and controlling contamination in food‐processing plants. Schaffner's model did not incorporate a dynamic representation of cross‐contamination events, including impact of food handlers’ behavior, in terms of compliance with personal hygiene and environmental sanitation practices.

Ivanek et al. (2004) developed a plant‐specific compartmental mathematical model of L. monocytogenes cross‐contamination to describe the bacterium's transmission among food, food‐contact surfaces, employees’ gloves, and the environment in a smoked‐fish‐processing plant. Using a difference‐equation system based on the Reed‐Forest model, Ivanek et al. (2004) described only changes in contamination prevalence of food products and different contact surfaces during facility operation and did not quantify the level of L. monocytogenes contamination in their model.

Mokhtari and Jaykus (2009) developed a compartmental exposure model that characterized the dynamics of foodborne transmission of human norovirus in the retail‐food‐preparation environment. They used sensitivity analysis to identify model inputs that provided the greatest contribution to risk of virus contamination of prepared foods, while what‐if scenario analyses were applied to evaluate the relative efficacy of potential control strategies for reducing the likelihood of contamination during food preparation. While their exposure model tracked spread of contamination in a retail facility and to prepared foods in time, it followed each food handler individually, and, hence, did not consider possible interactions between food handlers that may influence the extent of microbial cross‐contamination in the facility environment. Furthermore, the model did not include presence of contamination niches in the facility environment.

Pouillot et al. (2015) developed a quantitative risk assessment model to simulate the behavior of retail employees in a delicatessen department and used it to track L. monocytogenes potentially present in such an environment and its foods. Using a discrete‐event simulation framework, that model assessed the risk of foodborne invasive listeriosis associated with current practices in retail delis and examined the impact on that risk from mitigations meant to reduce or prevent L. monocytogenes growth or contamination in RTE foods prepared in retail deli settings. The model demonstrated flexibility and granularity in modeling microbial cross‐contamination events; however, it did not incorporate a dynamic activity schedule beyond a chronological sequence of events that might occur in a deli department. Furthermore, while the model presented a simplified approach for modeling contamination niches associated with certain objects in the deli department (e.g., slicer, floor), it did not include spatial distribution of contamination on food‐ and non–food‐contact surfaces or pathogen persistence and growth in the facility environment, including contamination niches.

Duret et al. (2017) developed a quantitative risk assessment model using a discrete‐event simulation framework to quantify and study the risk associated with norovirus transmission to consumers through food contaminated by infected food employees in a retail food setting. The risk assessment focused on the impact of ill food workers with symptoms of diarrhea and vomiting, and potential control measures against transmission of norovirus to foods. The model examined the behavior of food employees regarding exclusion from work while ill and after symptom resolution, and preventive measures limiting food contamination during preparation. Similar to the model developed by Pouillot et al. (2015), the model described the series of consecutive and predefined activities undertaken by the food employees. Furthermore, while food‐employee behavior, in terms of compliance with hygiene practices (e.g., handwashing, gloving), was explicitly modeled, spatial distribution of norovirus particles in the facility environment, including food‐ and non–food‐contact‐surface areas, also were considered, using a simplified approach that assumed all surface contamination was available for potential transfer upon tactile contact.

Mokhtari et al. (2018) developed a probabilistic mathematical model for the postharvest processing of leafy greens, focusing on Escherichia coli O157:H7 contamination of fresh‐cut romaine lettuce as the case study. That model could support the investigation of cross‐contamination scenarios and evaluate and compare different risk mitigation options. They used an ABM framework to predict pathogen prevalence and levels in bags of fresh‐cut lettuce and quantify spread of E. coli O157:H7 from contaminated lettuce to surface areas of processing equipment. Mokhtari et al.’s model focused on a series of sequenced cross‐contamination events that could occur during the postharvest processing of romaine lettuce (e.g., manual trimming, shredding, washing, and bagging), without including potential dynamic events representing interactions between facility employees and processing equipment. Their model did not include growth and harborage of pathogens on processing equipment and in contamination niches.

The current model, F2‐QMRA, provides a comprehensive, flexible framework with a dynamic scheduling component that can accommodate a wide range of possible activities that could result in persistence and spread of microbial pathogens in the facility environment and in cross‐contamination of prepared foods. F2‐QMRA is able to not only keep track of contamination changes over time in different environmental sampling zones, but also to explicitly model spatial distribution of contamination in the facility environment, including contamination niches that may harbor microbial pathogens. These features are intended to enable more accurate attribution of contamination in different sampling zones to likelihood and levels of contamination in prepared foods, providing a better understanding of microbial persistence and spread dynamics, and to identify optimal risk mitigation options for controlling potential cross‐contamination of prepared foods. Furthermore, the spatial design of our model provides capacity to readily simulate different facility environment setups beyond the initial scope that included four rooms. We also demonstrated the ABM framework's capabilities to provide insights about a complex, dynamic system in which the emergent phenomena (e.g., spread of microbial pathogens within a facility environment and in prepared foods) could be better predicted in a bottom‐up approach starting from individual agents (e.g., food handlers) and simulating their behaviors and interactions with other agents and the facility environment.

The ABM framework is preferable to other simulation techniques, such as discrete‐event simulation, when agents have both dynamic relationships with other agents and their behaviors are unique. For example, our ABM framework can simply accommodate different personal hygiene practices for each of the four food handlers included in the model and has the possibility of adjusting behaviors upon observing certain trigger events (e.g., increased sanitation frequency once a positive sample is collected from a contact‐surface area). The ABM framework is also preferred when individual agents have spatial aspects to their behaviors. For example, our model can accommodate different hygiene practices for individual food handlers based on their current locations (e.g., a food handler may demonstrate higher likelihood of compliance with sanitation practices in the food processing area than in other rooms in a food facility).

In an ideal world, all parameters in a risk assessment model are expressed quantitatively and are based on peer‐reviewed scientific studies. Unfortunately, comprehensive data often are not available, and there are instances in which data are completely lacking. In other cases, data from reasonable surrogates may be used in place of agent‐specific data (Vose, 2000). Accordingly, some approximations and assumptions were necessary in our modeling approach, which could result in aggregation errors as a source of uncertainty. Aggregation errors can arise because of approximations or assumptions used in a risk assessment model to simulate the underlying dynamic system (e.g., spread of L. monocytogenes contamination in a food facility). Aggregation errors can be examined as a source of systematic error through which bias may be introduced in the model outputs. For example, in the absence of experimental data, we used a conceptual model of persistence/no persistence of L. monocytogenes pathogens in contamination‐harborage sites, as suggested by Carpentier and Cerf (2011). We further assumed that a certain fraction of the pathogens placed in contamination niches (RF) would be released to more accessible surface areas during each cleaning event. This assumption actually may be consistent with what food operators have been reporting regarding their inability to totally eradicate L. monocytogenes from the facility environment because of the contamination‐harborage sites (Ferreira et al., 2014; Malley, Butts, & Wiedmann, 2015; Ortiz et al., 2010; Tompkin, Scott, Bernard, Sveum, & Gombas, 1999). In the absence of available data, we further considered four possible values ranging between 10−6 and 10−3 (See Table VIII) that were randomly assigned to RF during the Monte Carlo simulation of the model. Considering this range of values, the global sensitivity analysis identified RF as the third most important model input impacting the prevalence of L. monocytogenes contamination in prepared foods (Fig. 7). Larger values of RF were associated with higher L. monocytogenes contamination prevalence values in prepared foods because more pathogens hidden inside contamination niches could become readily available to transfer to food servings during facility operation. RF also impacted the persistence of L. monocytogenes in the facility environment, especially in holes and cracks of objects assigned to Zone 4. Because of lower sanitation efficacy values assigned to pathogens hidden inside the contamination niches (see Table VIII), smaller values of RF were associated with longer survival of pathogens in the facility environment (data not shown). We also assumed certain values for the number of tactile contacts between food handlers and objects assigned to different sampling zones during selected activities (Table VI). Unfortunately, limited data are available for estimating common behaviors of food handlers during facility operation, such as the number of contacts among hands (or gloves), different surfaces, and foods. Most of the estimates used in this study were, therefore, based on our exploratory assumptions. We did not include our exploratory assumptions related to the number of tactile contacts in the global sensitivity analysis; however, these model inputs, if included in sensitivity analysis, could significantly impact the persistence of L. monocytogenes pathogens in the facility environment and the spread of contamination to prepared foods. Unlike RF that can be challenging to populate with real data, observational studies can be used as a reliable means of data collection regarding behaviors of food handlers during different food‐preparation activities (Clayton & Griffith, 2004; Green et al., 2006; Lubran et al., 2010).

F2‐QMRA is the first version of a highly flexible model that can be used to evaluate pathogen persistence and transfer in food facilities of varying designs and with varying operations. This study demonstrated that developing a flexible tool that can be used as a “virtual laboratory” to evaluate the impact of mitigation options in a food facility is possible and valuable. When combined with a global sensitivity analysis and what‐if scenario analysis, this model aided in our understanding of the complex dynamics of microbial pathogen spread in a food‐facility environment and in preliminary evaluation of candidate mitigation options. As new data become available, model inputs can be updated (Buchanan, Gorris, Hayman, Jackson, & Whiting, 2017). Further, additional transfer routes, such as airborne/aerosol transmission (Byrne, Lyng, Dunne, & Bolton, 2008; Vontayson, 2009), may be added as data become sufficient to characterize these processes. With further work, such as development of a user‐friendly graphical interface, this model also could be a valuable tool for industry and regulatory agencies as they seek to reduce the burden of foodborne disease associated with microbiological contamination.