<figure><img src="image_1.png" /></figure>

HAL
open science

# Multipathogen quantitative risk assessment in raw milk soft cheese

Subhasish Basak, Laurent Guillier, Julien Bect, Janushan Christy, Fanny

Tenenhaus-Aziza, Emmanuel Vazquez

► To cite this version:

Subhasish Basak, Laurent Guillier, Julien Bect, Janushan Christy, Fanny Tenenhaus-Aziza, et al. Multipathogen quantitative risk assessment in raw milk soft cheese. Microbial Risk Analysis, 2024, 27-28, pp.100318. 10.1016/j.mran.2024.100318. anses-04823814

HAL Id: anses-04823814

https://anses.hal.science/anses-04823814v1

Submitted on 24 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Multipathogen quantitative risk assessment in raw milk soft cheese

Subhasish Basak<sup>a,b,*</sup>, Laurent Guillier<sup>a</sup>, Julien Bect<sup>b</sup>, Janushan Christy<sup>c</sup>, Fanny Tenenhaus-Aziza<sup>d</sup>, Emmanuel Vazquez<sup>b</sup>

## ARTICLE INFO

**Keywords:**

Raw milk soft cheese

Haemolytic uremic syndrome

Listeriosis

Salmonellosis

Shiga-toxin producing *Escherichia coli*

*Salmonella*

*Listeria monocytogenes*

Quantitative microbial risk assessment

## ABSTRACT

We propose a multipathogen Quantitative Microbiological Risk Assessment (QMRA) model to estimate the risk of foodborne illness from bacterial pathogens in raw milk soft cheese. Our work extends an existing QMRA model for pathogenic Shiga toxin-producing *Escherichia coli* (STEC) (Basak et al., under review; Perrin et al., 2014) by incorporating the effects of *Salmonella* and *Listeria monocytogenes*. This multipathogen model integrates microbial contamination of raw milk at the farm level, as well as the growth and survival of these bacteria during cheese fabrication, ripening, and storage. The public health impact of multipathogen risk associated with raw milk cheese consumption is assessed using Disability-Adjusted Life Years (DALYs). The model evaluates intervention strategies at both pre- and post-harvest stages to estimate intervention costs. Furthermore, it tests various scenarios of these strategies and optimizes intervention parameters to minimize multipathogen risk and associated costs. This article discusses challenges in QMRA model validation, emphasizes model limitations, and explores future perspectives for improvement.

## 1. Introduction

Microbiological food safety is a major challenge for the food sector (see, e.g., Plaza-Rodriguez et al., 2018). In this context, the microbiological food safety community—including food authorities, food industries, and food research institutes—have invested research efforts into the field of Quantitative Microbial Risk Assessment (QMRA). The aim is to establish risk-based control measures (see, e.g., Koutsoumanis and Aspridou, 2016). QMRA is a part of microbial risk analysis, which involves risk assessment, risk management, and risk communication (World Health Organization, 1997).

Microbial risk assessment enables the evaluation of the likelihood of illness caused by pathogenic microorganisms and environmental factors that impact microbial growth. According to Commission (1999), the framework for executing a QMRA for pathogens is built on several foundational components: hazard identification, exposure assessment, hazard characterization, risk characterization, and risk management options.

Hazard identification entails recognizing microbiological agents—such as bacteria, pathogens, and viruses—that exist in food and may

lead to adverse health effects. This identification process is initiated after the problem formulation (Ungaretti Haberbeck et al., 2018). It involves delineating a list of microbial pathogens associated with the risk assessment in question.

Exposure assessment identifies and characterizes the pathways of the exposure to the microbial hazards and quantifies the exposure levels to estimate the magnitude of intake via consumption of a given food commodity.

Hazard characterization then converts these levels of exposure into a probability, representing the likelihood of adverse health effects associated with the hazard.

Risk characterization consolidates this information, offering a quantitative estimation of the likelihood, along with its associated uncertainties, of experiencing known or potential adverse health effects in a specific population. This estimation draws on data from hazard identification, hazard characterization, and exposure assessment (Commission, 1999).

Generally, cheeses are considered safe and nutritious food, but food-borne illnesses related to cheese consumption can occur (see, e.g., Dubois-Brissonnet et al., 2022). In the microbial risk assessment

<footnote><sup>a</sup> Agence nationale de sécurité sanitaire de l'alimentation, de l'environnement et du travail (ANSES), 14 Rue Pierre et Marie Curie, Maisons-Alfort, 94700, France</footnote>

<footnote><sup>b</sup> Laboratoire des Signaux et Systèmes (L2S), Université Paris-Saclay, CNRS, CentraleSupélec, 3 Rue Joliot-Curie, Gif-sur-Yvette, 91192, France</footnote>

<footnote><sup>c</sup> Centre technique d'expertise agroalimentaire (ACTALIA), 419 Route des Champs-Laitiers, La Roche-sur-Foron, 74801, France</footnote>

<footnote><sup>d</sup> Centre national interprofessionnel de l'économie laitière (CNIEL), 42 Rue de Châteaudun, Paris, 75314, France</footnote>

<footnote>* Corresponding author.
E-mail address: subhasish.basak@centralesupelec.fr (S. Basak).</footnote>
<table><caption>**Table 1**
Notations and corresponding pathogen classification.</caption><thead><tr><th>Index notation</th><th>Pathogen membership set</th></tr></thead><tbody><tr><td>x</td><td>MPS-STEC, <i>Salmonella</i>, <i>Lm</i></td></tr><tr><td>x<sup>l</sup></td><td>MPS O157:H7, MPS nonO157:H7, HV <i>Salmonella</i>, <i>LV Salmonella</i>, <i>Lm</i></td></tr><tr><td>x<sup>d</sup></td><td>MPS O157:H7, MPS nonO157:H7, <i>Salmonella</i>, <i>Lm</i></td></tr><tr><td>x<sup>ll</sup></td><td>MPS-STEC, HV <i>Salmonella</i>, <i>LV Salmonella</i></td></tr></tbody></table>

literature, there exist a number of QMRA studies on the contamination of raw milk soft cheese (see, e.g., Sanaa et al., 2004; Tenenhaus-Aziza et al., 2014; Campagnollo et al., 2018; Ramos et al., 2021; Lindqvist et al., 2002), where the authors have proposed methods not only to compute bacterial prevalence and contamination at the time of consumption but also to identify major parameters contributing to the risk, using simulation studies (usually through “what-if” scenarios).

Among all food-borne pathogens, Shiga-toxin producing *Escherichia coli* (STEC) in soft cheese is a notable concern (Farrokh et al., 2013). The definition of the virulence potential of STEC is not straightforward (Lindqvist et al., 2023). Depending on the strains' characteristics and the exposed population, the symptoms can range from mild to severe illnesses, such as Haemolytic and Uremic Syndrome (HUS), which is a leading cause of renal failure in young children. The growth and survival of pathogenic STEC serotypes throughout different phases of cheese production, were studied by Maher et al. (2001) and Miszczycha et al. (2016). In this context, Perrin et al. (2014) proposed a stochastic QMRA model to assess the risk of HUS associated with the five Main Pathogenic Stereotypes of STEC (MPS-STEC) in raw milk soft cheeses, and explored the role of control measures for minimizing the risk of illness.

Building on the work of Perrin et al. (2014), we introduce a new farm-to-fork QMRA multipathogen model, integrating the impacts of three bacteria-specifically, STEC, non-typhoidal *Salmonella*, and *Listeria monocytogenes*-which can be potentially present in raw milk (see, e.g., Costanzo et al., 2020; Van et al., 2009; Sanaa et al., 2004). Like STEC, *Salmonella* and *Listeria monocytogenes* are hazards that, upon infecting humans, can trigger distinct health conditions. *Salmonella* can lead to a condition called salmonellosis, which is characterized by symptoms such as diarrhea, abdominal pain, fever, and vomiting. *Listeria monocytogenes*, on the other hand, can cause listeriosis, a more severe illness (see, e.g., Camargo et al., 2017; Leclercq et al., 2014), that can lead to fever, muscle aches, nausea, diarrhea, and in severe cases, can even result in meningitis or septicemia.

This new multipathogen QMRA model integrates the pre-harvest and post-harvest intervention steps in the cheese production processes, which are used as control measures. This integration will enable the evaluation of their impact on the risk of illness.

This article is structured as follows. Section 2 describes the model with the underlying assumptions and its components, along with their functionalities, developed according to the World Health Organization et al. (2021) framework. Section 3 elaborates on the implementation of the model, detailing the selection of model parameters and the mathematical and statistical techniques employed. Section 4 presents numerical results regarding bacterial prevalence and risk reduction, derived from exploring various intervention scenarios. Finally, Section 5 provides a discussion on the model's usage, applicability, and perspectives for future work.

## 2. Description of the proposed hierarchical model

### 2.1. Model overview

The multipathogen QMRA model builds upon the QMRA model for STEC proposed by Perrin et al. (2014) and the R implementation provided by Basak et al. (2024).

The QMRA model is used to assess the risk of microbial contamination by examining the evolution of pathogens throughout the entire cheese-making process, from the farm where milk is produced to the consumer's fork where it is consumed. This farm-to-fork model can be regarded as a stochastic simulator, a computational model used for simulating a complex system that incorporates inherent randomness. The model consists of two hierarchical components: a batch-level simulator and an output module. The batch-level simulator models all the various steps of the cheese manufacturing process, beginning with the collection of milk from a specific number of farms on a given day and continuing through the production process. The production process usually lasts up to 14 days until the cheese ripening step, followed by a cheese storage step until the 22nd day, after which the produced batch of cheese is sent to the market (see, e.g., Perrin et al., 2014). A typical batch of cheese usually contains 22,000-23,000 cheeses of 250 g, produced using a total volume of 50,000 l of raw milk, though it can vary depending on the cheese producer. The outputs of interest corresponding to a particular cheese batch is produced by the batch-level simulator, which are then used by the output module, that produces an estimate for the impact of the food-borne illnesses and the intervention costs.

*Notations and abbreviations* Throughout this article, we adopt a notation convention in which, depending on the variables used for the indices, it is implicitly assumed that these variables have values in different sets of pathogens, as summarized in Table 1.

Here, MPS-STEC denotes the Main Pathogenic Serotypes of Shiga Toxin Producing *Escherichia coli*, and it is classified into two subclasses based on serotype markers, namely O157:H7 and non O157:H7 (Perrin et al., 2014). For *Salmonella*, we consider two subclasses corresponding to their high and low virulent strains, respectively abbreviated by HV *Salmonella* and LV *Salmonella*. In mathematical notations, these are further compressed as HV-*Salmo* and LV-*Salmo*. For *Listeria monocytogenes*, the abbreviation *Lm* is used.

**Modules** The batch-level simulator has four modules: a farm module followed by a preharvest intervention step, a cheese production module, a consumer module, and a postharvest sampling module.

The sets of input parameters, as detailed in Section 3.2, for each of the modules, are aggregated in a vector denoted by $\theta = \{\theta^{\text{farm}}, \theta^{\text{cheese}}, \theta^{\text{con}}, \theta^{\text{post}}\}$, which forms the input parameter of the simulator.

The farm module, as explained in Section 2.3.1, models the collection of milk from different farms and outputs the initial pathogen concentration $Y_x^{\text{milk}}$ in the aggregated milk tank, for the pathogen $x$ (see Table 1). (Note that we use capital letters in the description of the model to designate random variables. This will be made explicit in the next sections.) The farm module also implements the preharvest intervention step, as explained in Section 2.5.1, which prevents contaminated farm milk from entering the production process.

The cheese module, as explained in Section 2.3.2, models the evolution of pathogens during the cheese production process and simulates the colony size $Y_{x^l}$ and the average number of colonies $\lambda_{x^l}^{\text{colony}}$ in a single cheese from that specific batch. Note that $Y_{x^l}$ and $\lambda_{x^l}^{\text{colony}}$ are random variables that depend on $Y_x^{\text{milk}}$ and $\theta^{\text{cheese}}$.

The consumer module, as described in Section 2.4.1, models the impact of the ingested dose of pathogen by the consumer and estimates the risk of illness.

The postharvest intervention step, as described in Section 2.5.2, implements the microbial cheese sampling plan and estimates the probability of detecting contamination in the produced batch of cheese.

Given the input vector $\theta$, the batch-level simulator yields the following outputs relevant to a particular batch: the milk loss per batch $M^{\text{batch}}$ due to the preharvest intervention step, the probability of rejecting the batch of cheese $p^{\text{batch}}$ due to the postharvest intervention step, and the batch risk $R_x^{\text{batch}}$, that is, the risk of encountering food-borne illnesses due to pathogen $x$, if a portion of 25 g of cheese is consumed from that specific batch of cheese. Again, note that $M^{\text{batch}}$,
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Batch level simulator: The four modules of the simulator, namely, farm, cheese, consumer and postharvest, are shown in pink boxes, along with their corresponding set of inputs in violet boxes. The outputs corresponding to a simulation of a single, namely batch are shown using orange boxes. The different modules are related hierarchically in the sense that the outputs of one model are used as inputs to other model(s) towards which it shows an arrow.</figcaption></figure>

<caption>**Table 2**</caption>

The estimated median values (Cassini et al., 2018) of years of life lived with disability (YLD) per 1000 cases, years of life lost (YLL) per 1000 cases, and DALY per 1 case, for STEC infections, listeriosis and salmonellosis. For MPS-STECA the DALY values are taken from ANSES (2020).

<table><thead><tr><th>Pathogen</th><th>YLL (1000 cases)</th><th>YLD (1000 cases)</th><th>DALY(1 case)<sub>x</sub></th></tr></thead><tbody><tr><td>MPS STEC infections</td><td>2700.0</td><td>1000.0</td><td>3.7</td></tr><tr><td>STEC infections</td><td>41.1</td><td>13.0</td><td>0.0541</td></tr><tr><td>Listeriosis</td><td>3300.0</td><td>400.0</td><td>3.7</td></tr><tr><td>Salmonellosis</td><td>15.0</td><td>4.0</td><td>0.019</td></tr></tbody></table>

$p^{\text{batch}}$ and $R_x^{\text{batch}}$ are random variables in what follows. When the batch level simulator is run once, it will simulate a sample value for these random variables, conditional on $\theta$.

Fig. 1 offers a schematic diagram of the batch-level simulator, encompassing several modules and simulating outputs corresponding to the production of a single batch of cheese. The output module, as explained in Section 2.6.1 and demonstrated by Fig. 9, is used to simulate several batches of cheese, that is, several sample values of $M^{\text{batch}}$, $p^{\text{batch}}$ and $R_x^{\text{batch}}$, and to estimate the final quantities of interest, e. g., the prevailing risk of illness, average milk loss, and average probability of rejecting a batch of cheese.

## 2.2. Hazard identification

The three pathogens identified as hazards in this study—the main pathogenic serotypes of Shiga-toxin producing *Escherichia coli* (also known as MPS-STECA), *Salmonella*, and *Lm*—can survive or grow during cheese making, particularly in raw milk soft cheeses (see, e.g., Costanzo et al., 2020). These bacterias are present in the intestines of lactating dairy animals, and can be transmitted through fecal matter to their udders, thereby contaminating milk during the milking process (Gopal et al., 2015).

**Main pathogenic serotypes of STEC** The five main pathogenic serotypes of STEC (MPS-STECA) identified thus far in Europe are O157:H7, O26: H11, O103:H2, O111:H8, and O145:H28. According to Panel et al (2020), EFSA (2017), 14 outbreaks involving STEC in milk, dairy, and cheese products were recorded between 2012 and 2017, affecting 775

individuals, with the primary causative agents being Shiga toxin-producing *E. coli* (STEC). In 2005, in France, an outbreak of Hemolytic Uremic Syndrome (HUS) was reported, linked to the contamination of raw milk soft cheese with *E. coli* O26 and O80, followed by another outbreak in 2009 among children up to 15 years (King et al., 2009) and in 2019 (Jones et al., 2019). However, other pathogenic serotypes have caused major outbreaks. In 2010, atypical STEC serotypes O104:H4 caused a large outbreak in Germany (see, e.g., Frank et al., 2011; EFSA, 2012). Recently, the O80 serotype has arisen in Europe (see, e.g., Bruyand et al., 2019).

**Listeria monocytogenes** The presence of *Lm* in raw milk and cheese has been extensively reported (see, e.g., Dalzini et al., 2016), with its widespread occurrence and potential for contamination possible at any stage of the production chain. Due to cold tolerance, i.e., the ability to grow at refrigeration temperatures as low as -1.5°C, and its capacity to form resilient biofilms resistant to sanitation, *Lm* exhibits the ability to persist and survive in various environments (see, e.g., McIntyre et al., 2015). According to Authority and Prevention (2022), there were 2183 confirmed invasive human cases of listeriosis in 2021. Cheese was estimated to be the origin of 1% of the listeriosis cases in Europe according to risk assessment models (Ricci et al., 2018).

**Salmonella** Salmonellosis is recognized as one of the most common bacterial food-borne illnesses in humans, with several outbreaks reported from *Salmonella* contamination in raw milk cheese (see, e.g., Ung et al., 2019; Robinson et al., 2020). *Salmonella* is commonly found in the intestinal tract of lactating animals, and milk contamination primarily occurs during the milking process (Ruzante et al., 2010). Among various serotypes, *Salmonella* Dublin, *Salmonella* Newport, and *Salmonella* Typhimurium are commonly linked to salmonellosis in both calves and adult cows, leading to varying degrees of illness. Additionally, dairy animals have been found to carry *Salmonella* serotypes such as *Salmonella* Cerro, *Salmonella* Kentucky, *Salmonella* Mbandaka, and *Salmonella* Montevideo without showing any symptoms, while shedding the bacteria in their feces (see, e.g., Bonifait et al., 2021; Van Kessel et al., 2012).
<caption>**Table 3**</caption>

The inputs for the Farm module, collectively represented as ($\theta^{\text{farm}}$), are categorized into two groups, visually separated by a dashed line. The first part includes inputs that the user can modify, while the second part comprises fixed parameters of the model.

<table><thead><tr><th>Symbol</th><th>Description</th><th>Values</th></tr></thead><tbody><tr><td>$N^{\text{farms}}$</td><td>Number of farms</td><td>31</td></tr><tr><td>$N_i^{\text{cow}}$</td><td>Number of cows in ith farm</td><td>Table 4</td></tr><tr><td>$q^{\text{milk}}$</td><td>Avg. quantity of milk from a cow</td><td>25 Liters</td></tr><tr><td>$f^{\text{sorting}}$</td><td>Milk testing frequency</td><td>10 days</td></tr><tr><td>$p^{\text{sorting}}$</td><td>Max. limit of E. coli conc.</td><td>50 CFU/mL</td></tr><tr><td>$\alpha_i, \sigma_i$</td><td>Hygiene parameter for ith farm</td><td>Section 3.1.1</td></tr><tr><td>$a^{\text{weibull}}, b^{\text{weibull}}$</td><td>Param. of dist. of STEC in feces</td><td>0.264, 16.288 (Perrin et al., 2014)</td></tr><tr><td>$\mu_{\text{ecoli}}, \tau_{\text{ecoli}}$</td><td>Mean &amp; sd of E. colist. in feces</td><td>6, 0.3 (Perrin et al., 2014)</td></tr><tr><td>$\mu_u, \tau_u$</td><td>Param. for estimating $p_{\text{STEC}}$</td><td>-0.927, 1.47411 (Perrin et al., 2014)</td></tr><tr><td>$\mu^{\text{Lm}}, \tau^{\text{Lm}}$</td><td>Parameters of Lmdist. in milk</td><td>-7.178, 0.552 (Sanaa et al., 2004)</td></tr><tr><td>$p_{\text{MPS-STEC}}$</td><td>Prop. of MPS-STEC carriers</td><td>0.025 (Perrin et al., 2014)</td></tr><tr><td>$p_{\text{Salmo}}$</td><td>Prop. of Salmonella infected cows</td><td>0.03 (Bonifait et al., 2021)</td></tr><tr><td>$p_{\text{HV-Salmo}}$</td><td>Prop. of HV Salmonella carriers</td><td>0.33 (Bonifait et al., 2021)</td></tr></tbody></table>

## 2.3. Exposure assessment

### 2.3.1. Farm module

At the farm level, milk is collected from a mixed herd of lactating animals, which can potentially includes both infected and non-infected cows, and stored in a bulk tank. During the milking process, pathogens potentially present in the fecal matter of infected cows can be transmitted into the bulk tank through their contaminated udders. The Bulk Tank Milk (BTM) is then collected from different farms and mixed into an Aggregated Milk Tank (ATM), which is used for cheese production. In the QMRA model, the farm module replicates this entire scenario of milk collection and simulates the concentration in CFU/ml (Colony forming units per mili liters) of all the pathogens in the ATM. The input parameters of the farm module are listed in Table 3.

The concentration $Y_{\text{Lm}}^{\text{milk}}$ of Lm in the ATM was studied by Sanaa et al. (2004). Following this work, we model $Y_{\text{Lm}}^{\text{milk}}$ as a log-normal random variable:

$$Y_{\text{Lm}}^{\text{milk}} | \theta^{\text{farm}} \sim \text{Lognormal}(\mu_{\text{Lm}}, \tau_{\text{Lm}}). \qquad (1)$$

Due to very low contamination level, insufficient knowledge and the unavailability of a reliable method to directly determine the pathogen concentration in milk for STEC and Salmonella, we employ the indirect approach of estimating concentration used by Perrin et al. (2014) for the QMRA model for STEC. This approach relies on collected data on E. coli concentration in farm milk, assuming that E. coli and other pathogen strains follow the same fecal routes.

Suppose that, for each farm indexed by $i = 1, 2, \dots, N^{\text{farms}}$, the milk is collected into a BTM from $N_i^{\text{cow}}$ cows; then the concentration of pathogen $x^{\ddagger}$ in the BTM corresponding to farm $i$ is denoted by $Y_{x^{\ddagger},i}^{\text{milk}}$ and is obtained as,

$$Y_{x^{\ddagger},i}^{\text{milk}} = Y_i^{\text{EC}} \cdot \frac{\overline{F_i^{x^{\ddagger}}}}{\overline{F_i^{\text{EC}}}}, \qquad (2)$$

where $Y_i^{\text{EC}}$ denotes the concentration (CFU/mL) of E. coli in BTM, $\overline{F_i^{x^{\ddagger}}}$ and $\overline{F_i^{\text{EC}}}$ respectively denote the average concentration (CFU/gram) in the fecal matter for the $x^{\ddagger}$ pathogen and E. coli, coming from all the cows.

To obtain $Y_i^{\text{EC}}$, $\overline{F_i^{x^{\ddagger}}}$ and $\overline{F_i^{\text{EC}}}$, the module first models the number of

infected cows using a binomial distribution $k_i^j | \theta^{\text{farm}} \sim \text{Binomial}(N_i^{\text{cow}}, p_j)$, in the ith farm, with the corresponding class probabilities $p_j$ of the two major pathogen classes $j \in \{\text{STEC}, \text{Salmonella}\}$. For STEC, the proportion $p_{\text{STEC}}$ of infected cows is estimated using a logit model (following Perrin et al., 2014). Next, the number of cows affected by the subclass serotypes is modeled, using the respective class probabilities:

$$k_i^{\text{MPS-STEC}} | \theta^{\text{farm}} \sim \text{Binomial}(k_i^{\text{STEC}}, p_{\text{MPS-STEC}}), \qquad (3)$$

$$k_i^{\text{HV-Salmo}} | \theta^{\text{farm}} \sim \text{Binomial}(k_i^{\text{Salmo}}, p_{\text{HV-Salmo}}). \qquad (4)$$

The number of Salmonella-infected cows carrying the low virulent serotype, is computed as $k_i^{\text{LV-Salmo}} = k_i^{\text{Salmo}} - k_i^{\text{HV-Salmo}}$. For the jth infected cow, $1 \le j \le k_i^{\ddagger}$, in the ith farm, $1 \le i \le N^{\text{farm}}$, the concentration in the fecal matter is modeled according to a Weibull distribution for MPS-STEC (Perrin et al., 2014),

$$F_{i,j}^{\text{MPS-STEC}} | \theta^{\text{farm}} \sim \text{Weibull}(a^{\text{weibull}}, b^{\text{weibull}}),$$

and the concentrations of the two serotypes HV-Salmonella, LV-Salmonella are modeled using a log-normal distribution (Bonifait et al., 2021),

$$\log_{10}(F_{i,j}^{\text{Salmo}}) | \theta^{\text{farm}} \sim \mathcal{N}(a^{\text{Salmo}}, b^{\text{Salmo}}).$$

For each farm, the average concentration of pathogen $x^{\ddagger}$ is

$$\overline{F_i^{\ddagger}} = \frac{1}{N_i^{\text{cow}}} \sum_{j=1}^{k_i^{\ddagger}} F_{i,j}^{x^{\ddagger}}.$$

The concentration of E. coli (CFU/mL) in a BTM $Y_i^{\text{EC}}$, is modeled by a log-normal distribution:

$$Y_i^{\text{EC}} | \theta^{\text{farm}} \sim \text{Lognormal}(\alpha_i, \sigma_i). \qquad (5)$$

The average $\overline{F_i^{\text{EC}}}$ of individual E. coli concentrations in fecal matter for each cow, denoted by $F_{i,j}^{\text{EC}}$, $j = 1, 2, \dots, N_i^{\text{cow}}$, is obtained using the model

$$\log10(F_{i,j}^{\text{EC}}) | \theta^{\text{farm}} \sim \mathcal{N}(\mu_{\text{ecoli}}, \tau_{\text{ecoli}}), \\ \overline{F_i^{\text{EC}}} = \frac{1}{N_i^{\text{cow}}} \sum_{j=1}^{N_i^{\text{cow}}} F_{i,j}^{\text{EC}}. \qquad (6)$$

Each of the BTMs is tested (also known as milk sorting, see Section 2.5.1) for E. coli concentration and accordingly accepted or rejected for cheese production. Let $S$ denote the set of farms that qualify after milk sorting and let $N_i^{\text{farms sorted}} = |S|$. After milk sorting, milk from all the qualified BTMs is collected into a single ATM. The final concentration (CFU/mL) of pathogen $Y_{x^{\ddagger},i}^{\text{milk}}$ in this ATM can be written as

$$Y_{x^{\ddagger},i}^{\text{milk}} = \sum_{i=1}^{N^{\text{farms}}} \left( Y_{x^{\ddagger},i}^{\text{milk}} \cdot \frac{V_i 1_{\{i \in S\}}}{\sum_{i=1}^{N^{\text{farms}}} V_i 1_{\{i \in S\}}} \right), \qquad (7)$$

where $V_i$ is the volume of milk in liters produced by the ith farm.

**Module outputs** The farm module models the concentrations (CFU/mL) of four different pathogenic serotypes in milk in the ATM, namely, MPS-STEC, HV Salmonella, LV Salmonella and Lm. It also yields the milk loss $M^{\text{batch}}$ (in Liters), due to the preharvest milk testing step, associated with the production of that particular batch of cheese. In addition it models the number of farms discarded due to milk testing $N^{\text{farms}} - N^{\text{farms sorted}}$, and the total volume of milk put in production, or in other words the volume of milk in the ATM $\sum_{i=1}^{N^{\text{farms}}} V_i 1_{\{i \in S\}}$.
<table><caption>**Table 5**
Cheese module inputs, collectively denoted as $\theta^{\text{cheese}}$.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Values</th></tr></thead><tbody><tr><td>$\theta(\mu^{\text{max}})$</td><td>Cardinal params. for $\mu^{\text{max}}$</td><td>Table 6</td></tr><tr><td>$\mu^{\text{opt}}_{\text{MPS-STEC}}$</td><td>Optimal growth rate for MPS-STECS</td><td>1.85 (Perrin et al., 2014)</td></tr><tr><td>$\mu^{\text{opt}}_{\text{Salmo}}$</td><td>Optimal growth rate for Salmonella</td><td>1.02 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$\mu^{\text{opt}}_{\text{Lm}}$</td><td>Optimal growth rate for Lm</td><td>0.55 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$y^{\text{max,milk}}$</td><td>Hypothetical max population in milk</td><td>10<sup>9</sup> CFU/mL (Perrin et al., 2014)</td></tr><tr><td>$y^{\text{max,cheese}}$</td><td>Hypothetical max population in cheese</td><td>10<sup>5</sup> CFU/g (Perrin et al., 2014)</td></tr><tr><td>$p_{\text{O157}}$</td><td>Class probability of MPS O157: H7</td><td>0.76 (Perrin et al., 2014)</td></tr><tr><td>$p_{\text{HV-Salmo}}$</td><td>Class probability of HV Salmonella</td><td>$Y_{\text{HV-Salmo}}^{\text{milk}} / (Y_{\text{HV-Salmo}}^{\text{milk}} + Y_{\text{LV-Salmo}}^{\text{milk}})$</td></tr><tr><td>$p_{\text{LV-Salmo}}$</td><td>Class probability of LV Salmonella</td><td>1 - $p_{\text{HV-Salmo}}$</td></tr><tr><td>$p_{\text{ colony HV-Salmo}}$</td><td>Consumption probability of HV Salmonella</td><td>$\lambda_{\text{HV-Salmo}}^{\text{colony}} / (\lambda_{\text{HV-Salmo}}^{\text{colony}} + \lambda_{\text{LV-Salmo}}^{\text{colony}})$</td></tr><tr><td>$p_{\text{ colony LV-Salmo}}$</td><td>Consumption probability of LV Salmonella</td><td>1 - $p_{\text{HV-Salmo}}^{\text{colony}}$</td></tr><tr><td>$\rho_{\text{O157}}$</td><td>MPS O157:H7 decline rate</td><td>0.14 (log<sub>10</sub> CFU/day) (Perrin et al., 2014)</td></tr><tr><td>$\rho_{\text{O157}}$</td><td>MPS nonO157:H7 decline rate</td><td>0.033 (log<sub>10</sub> CFU/day) (Perrin et al., 2014)</td></tr><tr><td>$\delta_{\text{core}}$</td><td>Core Salmonelladecline rate</td><td>1.4 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$\tau_{\text{core}}^{\text{d}}$</td><td>Inter cheese variability of Salmonella</td><td>0.021568496 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$\delta_{\text{rind}}$</td><td>Rind Salmonelladecline rate</td><td>3.1 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$\tau_{\text{rind}}^{\text{d}}$</td><td>Inter cheese variability of Salmonella</td><td>0.022373347 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$p_{\text{core}}$</td><td>Param. for core Salmonella decline</td><td>0.274 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$p_{\text{rind}}$</td><td>Param. for rind Salmonella decline</td><td>2.7 (Gonzales-Barron et al., 2022)</td></tr><tr><td>$\rho_{\text{Lm}}$</td><td>L. mono decline rate</td><td>Triangular(0.5,1,2) (Sanaa et al., 2004)</td></tr><tr><td>$p_{\text{Hlm}}$</td><td>pH for L. mono second growth phase</td><td>6 (ArtiSaneFood challenge tests)</td></tr><tr><td>$T_{\text{Lm}}$</td><td>Temp. for L. mono second growth phase</td><td>12.45 (ArtiSaneFood challenge tests)</td></tr><tr><td>$pH^{\text{storage}}$</td><td>Industrial params. for storage step</td><td>Table III in Perrin et al. (2014)</td></tr><tr><td>$d^{\text{storage}}$</td><td>Storage duration (hours)</td><td>Triangular(1,12,40) (Perrin et al., 2014)</td></tr><tr><td>$T^{\text{storage}}$</td><td>Storage temperature (°C)</td><td>Uniform(4,6) (Perrin et al., 2014)</td></tr><tr><td>$\{d, pH, T\}$</td><td>Physico-chemical parameters</td><td>Table III in Perrin et al. (2014)</td></tr><tr><td>$aw$</td><td>Parameter for water activity</td><td>0.99, Table III in Perrin et al. (2014)</td></tr><tr><td>$v^{\text{cheese}}$</td><td>Milk used in a single cheese</td><td>2200 mL (Basak et al., 2024)</td></tr><tr><td>$w^{\text{loss}}$</td><td>Proportion of water loss in molding</td><td>0.9 (Perrin et al., 2014)</td></tr><tr><td>$t^{\text{consum}}$</td><td>Consumption time</td><td>Triangular(22,30,60) (Basak et al., 2024)</td></tr></tbody></table>

### 2.3.2. Cheese module

The inputs for the cheese module are the initial pathogen concentrations from the farm module and a set of parameters denoted by $\theta^{\text{cheese}}$ detailed in Table 5.

After being collected in the ATM, the milk undergoes several processing steps, including milk storage, molding, draining, salting, ripening, and cheese storage. These steps can be categorized into two phases, namely the liquid phase and the solid phase, respectively illustrated by two schematic diagrams in Figs. 2 and 3.

During the premolding steps in the liquid phase, a growth in the bacteria concentration is observed, for all the three pathogens. At the end of the liquid phase, the pathogen cells are presumed to become immobilized within the cheese matrix, resulting in the formation of colonies. At this stage, the cheese module estimates the average number of colonies formed in a single cheese. After entering the solid phase, the cheese module models the evolution of the colonies starting from one single colony. Finally it estimates the colony sizes for different pathogens at the time of consumption. The rest of this section is organized in three different parts, corresponding to the evolution of pathogens during the different cheese processing steps.

**Growth phase** For all three pathogens, the growth phase occurs until the salting step. Starting from the farm module outputs, the cheese module models the growth of the pathogens over the different steps using an ordinary differential equation

$$ \begin{cases} \frac{dY_x}{dt} = \mu_x^{\text{max}}(t) \cdot Y_x(t) \cdot \left(1 - \frac{Y_x(t)}{y^{\text{max}}}\right), \\ Y_x(0) = Y_x^{\text{milk}} \end{cases} \quad (8) $$

where $\mu_x^{\text{max}}(t)$ stands for the maximum growth rate (in h<sup>-1</sup>) and $y^{\text{max}}$ is a parameter that represents the hypothetical maximum population of pathogen strains in milk or cheese. The maximum growth rate $\mu_x^{\text{max}}(t)$, as shown in Fig. 5a, is modeled against time $t$ according to Augustin et al. (2005) using the optimal growth rate parameter $\mu_x^{\text{opt}}$, several physico-chemical parameters {$d, pH, T, aw$}, and their nominal values, collectively denoted as $\theta(\mu^{\text{max}})$ and listed in Table 6. Fig. 4 shows the variation in physico-chemical parameters over time, as studied by Perrin et al. (2014).

**Remark 1.** For *Salmonella*, the overall initial concentration $Y_{\text{Salmo}}^{\text{milk}}$ is computed as the sum of the concentrations of its two subclass serotypes.

The concentrations during the liquid phase, specifically at the end of the storage step $Y_x^{\text{storage}}$ and molding step $Y_x^{\text{molding}}$, are computed using (8), as depicted in Fig. 5a.

Corresponding to the three major pathogens, we consider a total of five different subclass serotypes denoted by $x^i$.

At the end of the liquid phase, conditional on the concentration of the molding step $Y_x^{\text{molding}}$, the number of colonies for each pathogen subclass serotypes are modeled as a Poisson variable $N_{x^i}^{\text{colony}} | Y_x^{\text{molding}} \sim \text{Poisson}(\lambda_{x^i}^{\text{colony}})$, with the mean computed as

$$ \lambda_{x^i}^{\text{colony}} = Y_x^{\text{molding}} \cdot v^{\text{cheese}} \cdot w^{\text{loss}} \cdot p_{x^i}, \quad (9) $$

where for MPS-STECS and *Salmonella*, $p_{x^i}$ denotes the class probability of the corresponding subclass serotype, whereas $p_{\text{Lm}} = 1$ for Lm. The parameters $v^{\text{cheese}}$ and $w^{\text{loss}}$ denote the amount of milk used for a single cheese and the proportion of water loss during the molding step,

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Liquid phase steps: storage and molding, models the pathogen concentration, starting from initial concentration $Y_{\text{milk}}^{\text{molding}}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Solid phase steps: draining, salting, ripening and cheese storage, models the evolution of colony size, starting from 1 bacteria.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4. Dynamic physico-chemical parameters, namely pH and temperature, during milk storage, molding, draining and salting steps of cheese production, separated by blue dotted lines.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. The three pathogens STEC, *Salmonella* and *Lm* are shown in red, blue and black lines respectively. The evolution of bacteria is calculated starting from average values of initial concentrations, for a baseline scenario with no preharvest intervention step.</figcaption></figure>

<table><caption>**Table 6**
Cardinal parameters for the computation of μ<sup>max</sup> for MPS-STEC (Perrin et al., 2014), *Salmonella* (Gonzales-Barron et al., 2022) and *Lm* (Sanaa et al., 2004).</caption><thead><tr><th>Symbol</th><th>MPS-STEC</th><th>Salmonella</th><th>Lm</th></tr></thead><tbody><tr><td>T<sub>min</sub></td><td>5.5</td><td>3.4</td><td>-1.7</td></tr><tr><td>T<sub>opt</sub></td><td>40.6</td><td>38.5</td><td>37</td></tr><tr><td>T<sub>max</sub></td><td>48.1</td><td>46</td><td>45.5</td></tr><tr><td>pH<sub>min</sub></td><td>3.9</td><td>4</td><td>4.71</td></tr><tr><td>pH<sub>opt</sub></td><td>6.25</td><td>7</td><td>7.1</td></tr><tr><td>pH<sub>max</sub></td><td>14</td><td>9</td><td>9.61</td></tr><tr><td>aw<sub>min</sub></td><td>0.9533</td><td>0.94</td><td>0.913</td></tr><tr><td>aw<sub>opt</sub></td><td>0.999</td><td>0.99</td><td>0.997</td></tr></tbody></table>

respectively.

Starting from the draining step, with an initial size of 1 CFU, (8) models growth of the size of that single colony forming unit, until the salting step. The evolution of each colony inside a particular cheese (of weight 250 g) is assumed to be identical during the growth phase due to the same environmental conditions.

**Decline phase** At the end of the salting step, when the batch of cheese enters the ripening step, there is a decline in the population of the colonies. The ripening step lasts until the 14th day of the production (Perrin et al., 2014), followed by a cheese storage step until the 22-nd day (Basak et al., 2024), and after that the cheese is sent to the market. The cheese is consumed on the *t*<sup>consum</sup>th day, which is a triangular-distributed random variable with minimum possible value 22.

For MPS-STEC and *Salmonella* this decline phase is considered as a continuous process which lasts until the cheese is consumed, whereas for *Lm*, only a decline in the colony population is considered. Let *t*<sup>salting</sup> be

the time taken (in hours) until the salting step, then for any time point *t* > *t*<sup>salting</sup>, the cheese module computes four reference colony sizes *Y*<sub>x<sup>1</sup></sub>(*t*), for the four pathogenic serotypes denoted by *x*<sup>1</sup>. These reference colony sizes represent the size of the respective pathogen colonies inside a particular cheese in the batch, without taking into account the inter-cheese variability. Here, we assume that all colonies within a particular cheese have the same decline phase, indicating no intra-cheese variability.

Following Perrin et al. (2014), the decline phase for MPS-STEC is modeled differently for its different subclass serotypes, indexed by *s* ∈ {MPS O157:H7, MPS nonO157:H7}. For MPS-STEC the reference colony size at time point *t* (in days), is computed as *Y*<sub>s</sub>(*t*) = *Y*<sup>salting</sup><sub>MPS-STEC</sub>·10<sup>-ρ<sub>s</sub>·(*t*−*t*<sup>salting</sup>)/24</sup>, with the subclass serotypes distinguished by the decline rate parameter *ρ*<sub>s</sub>.

The decline phase for *Salmonella* is dependent on the position of the colonies in the cheese matrix. The decline is modeled using different set of parameters for the core region colonies and the rind region colonies (Gonzales-Barron et al., 2022), namely, {δ<sub>s</sub>, *p*<sub>s</sub>}, for *s* ∈ {core, rind}. The reference colony size at time *t*, for serotype *s* is computed as *Y*<sub>s</sub>(*t*) = *Y*<sup>salting</sup><sub>Salmo</sub>·10<sup>-((t−*t*<sup>salting</sup>)/24)/δ<sub>s</sub></sup>·*ρ*<sub>s</sub>. The reference colony size *Y*<sub>Salmo</sub>(*t*) for *Salmonella* is obtained by averaging over the core and rind region with respect to the proportion of colonies in those region *p*<sub>core,rind</sub>.

Following Sanaa et al. (2004), for *Lm*, there is a decrease in the population after the salting step during curd acidification in cheese vats and molds, and it is modeled using log apparent kills, *Y*<sup>post.salting</sup><sub>Lm</sub> = *Y*<sup>salting</sup><sub>Lm</sub>·10<sup>-ρ<sub>Lm</sub></sup>, with a triangular distributed decline parameter *ρ*<sub>Lm</sub>.

**Second growth phase** The second growth phase is only observed for *Lm* (Sanaa et al., 2004), and it initiates depending on the position of the colonies in the core and rind region. Challenge test data collected under
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Evolution of colony sizes starting from unit size, for MPS-STEC and *Salmonella*, during the draining, salting and ripening steps, separated by blue dotted lines.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Evolution of colony sizes starting from unit size, for *Lm*, during the draining, salting, ripening and cheese storage step, separated by blue dotted lines. The red and orange dotted lines respectively indicates the starting time of the second growth step, corresponding to the rind and core region colonies.</figcaption></figure>

<table><caption>**Table 7**
Inputs of the consumer module collectively denoted as $\theta^{con}$.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Values</th></tr></thead><tbody><tr><td>k, r<sub>0</sub></td><td>Param. in STEC dose-response</td><td>0.38, 1e-2.33</td></tr><tr><td>age<sup>max</sup></td><td>Maximum age group</td><td>15</td></tr><tr><td>w<sup>cheese</sup></td><td>weight of a single cheese</td><td>250 g</td></tr><tr><td>w<sup>serving</sup></td><td>weight of a single serving</td><td>25 g</td></tr><tr><td>g(age)</td><td>Prop. of cheese consumed per age</td><td>(Perrin et al., 2014)</td></tr><tr><td>μ<sub>O157</sub></td><td>Mean of ϵ<sub>O157</sub></td><td>0 (Basak et al., 2024)</td></tr><tr><td>τ<sub>O157</sub></td><td>SD. of ϵ<sub>O157</sub></td><td>0.000279659 (Perrin et al., 2014)</td></tr><tr><td>μ<sub>0157</sub></td><td>Mean of ϵ<sub>0157</sub></td><td>0 (Basak et al., 2024)</td></tr><tr><td>τ<sub>0157</sub></td><td>SD. of ϵ<sub>0157</sub></td><td>0.000065399 (Perrin et al., 2014)</td></tr><tr><td>α<sub>HV-Salmo</sub></td><td>HV Salmonella dose-response parameter</td><td>0.132 (Strickland et al., 2023)</td></tr><tr><td>α<sub>LV-Salmo</sub></td><td>LV Salmonella dose-response parameter</td><td>0.318 (Strickland et al., 2023)</td></tr><tr><td>β<sub>HV-Salmo</sub></td><td>HV Salmonella dose-response parameter</td><td>51.45 (Strickland et al., 2023)</td></tr><tr><td>β<sub>LV-Salmo</sub></td><td>LV Salmonella dose-response parameter</td><td>4729.9 (Strickland et al., 2023)</td></tr><tr><td>θ<sub>g</sub>, r</td><td>Lm dose-response parameters</td><td>(Ricci et al., 2018; Pouillot et al., 2015)</td></tr><tr><td>N<sup>dose</sup></td><td>Monte Carlo sample size</td><td>10000</td></tr></tbody></table>

the ArtiSaneFood project were used to identify the favorable environmental conditions, namely, the physico-chemical parameters $pH_{Lm}$, $T_{Lm}$ that initiate the second growth step. These parameters were used to estimate the maximum growth rate for *Lm*, and the colony size for the core and rind regions at time *t* was modeled separately using (8), starting

from an initial concentration $Y_{Lm}^{post salting}$. It was observed that colonies in the core region (comprising 90% of the total colonies) enter the second growth phase after 20 days of the salting step, while the colonies at the surface region (making up 10% of the total colonies) initiate their second growth phase after 7 days of the end of salting step, or equivalently, the start of the ripening step. The reference colony size $Y_{Lm}(t)$ is obtained by averaging over the core and rind regions, similarly as *Salmonella*.

The evolution of colonies for the three pathogens, during the solid phase, is depicted in Figs. 6a, b, and 7.

**Module outputs** For all the five subclass serotypes, the expected number of colonies $\lambda_{x^l}^{colony}$ at a particular time point *t* is adjusted with respect to the corresponding reference colony size of the subclass serotype (for example, $Y_{O157}(t)$ and $Y_{0157}(t)$ for MPS-STEC) or the reference colony size of the pathogen itself (for example, $Y_{Salmo}(t)$ and $Y_{Lm}(t)$ for *Salmonella* and *Lm* respectively), at time *t*. If the reference colony size $Y_{x^l}(t)$ falls below 1, we choose the assumption that the colonies have disappeared with a probability $Y_{x^l}(t)$, and the corresponding adjusted expected number of colonies is obtained as $\lambda_{x^l}^{colony}(t) = \lambda_{x^l}^{colony} \cdot Y_{x^l}(t)$. If $Y_{x^l}(t) > 1$, the expected number of colonies remains unchanged.

The outputs of the cheese module are the adjusted expected number of colonies $\lambda_{x^l}^{colony}(t)$, for pathogen $x^l$ and the corresponding reference size of the colonies $Y_{x^l}(t)$, for pathogen $x^l$. These quantities are computed at two specific time points, firstly, the time of testing samples from the cheese batches $t^{test}$, and secondly at the time of consumption $t^{consum}$.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Directed acyclic graph showing the dependence between the model variables. The dependent variables are positioned at the arrow tips, signifying their dependence on the specific variables, which are situated at the arrow ends. No arrows between any two variables signifies independence. This graph is used to simplify the conditional expectations defined later in this chapter, using the dependence relationship.</figcaption></figure>

## 2.4. Hazard characterization

### 2.4.1. Consumer module

The consumer module describes the nature and probability of adverse human effects as a function of viable pathogen numbers ingested, termed the dose. The consumer module uses the outputs of the cheese module, that is, the adjusted expected number of colonies $\lambda_{x^i}^{\text{colony}}(t^{\text{consum}})$ in a cheese and the reference size of one colony $Y_{x^i}^{\text{colony}}(t^{\text{consum}})$, at the time of consumption $t^{\text{consum}}$. In addition to the cheese module outputs, the other inputs of the consumer module, denoted by $\theta^{\text{con}}$, are listed in Table 7. The consumer module simulates the risk of the corresponding illness based on the cheese consumption behavior of people in different age groups.

**Inter-cheese variability** In the modeling of risk, considering inter-cheese variability is crucial as it accounts for the differences observed in the distribution of pathogenic colonies across different cheese samples. Specifically, even though the reference colony sizes $Y_{x^i}(t)$, as mentioned in Section 2.3.2, are defined without the inter-cheese variability, this variability is explicitly incorporated during risk computation.

For MPS-STEC, the final colony sizes for the two subclass serotypes, with the inclusion of the inter-cheese variability, follow a log-normal distribution (Perrin et al., 2014), described by $Y_s^{\text{colony}} = Y_s(t^{\text{consum}}) \cdot 10^{\delta_s}$, where $\delta_s \sim N(\mu_{\delta_s}, \tau_{\delta_s})$ represents the parameter for inter-cheese variability, with $s$ denoting the two subclass serotypes of MPS-STEC $s \in \{\text{MPS O157:H7}, \text{MPS nonO157:H7}\}$. While this assumption is made for MPS-STEC, we still posit that there is no intra-cheese variability, meaning all the colonies inside a specific cheese share an identical size $Y_s^{\text{colony}}$ for serotype $s$ (Perrin et al., 2014).

For *Salmonella*, the inter-cheese variability could be introduced through the variance of the decline step parameter $\delta_s$, for $s \in \{\text{core, rind}\}$, denoted by $\tau_s^\delta$ (Gonzales-Barron et al., 2022), as seen in Table 5. However, the current model assumes no inter-cheese variability for *Salmonella* and for *Lm* as well, resulting in their final colony size being equivalent to the reference colony size.

**Model for the dose** The risk of getting the food-borne illness is dependent on the number of cells of pathogen $x^i$ ingested by the consumer, which is referred to as the dose, denoted by $\Gamma_{x^i}$. The dose is

defined as the number of pathogenic cells present in a particular cheese serving of 25 g, which is obtained by multiplying the size of colonies $Y_{x^i}^{\text{colony}}$ with the number of colonies $N_{x^i}^{\text{colony-serving}}$ in a cheese serving, distributed as

$$N_{x^i}^{\text{colony-serving}} | \lambda_{x^i}^{\text{colony}}, t^{\text{consum}} \sim \text{Poisson} \left( \lambda_{x^i}^{\text{colony}} (t^{\text{consum}}) \frac{W^{\text{serving}}}{W^{\text{cheese}}} \right). \quad (10)$$

**Probabilities of illness** To obtain probabilities of illness, three different dose-response models were used corresponding to the three different pathogens. Following Perrin et al. (2014), an exponential dose-response model was used for MPS-STEC, that uses the combined dose $\Gamma_{\text{MPS-STEC}} = N_{\text{O157}}^{\text{colony-serving}} \cdot Y_{\text{O157}}^{\text{colony}} + N_{\text{O157}}^{\text{colony-serving}} \cdot Y_{\text{O157}}^{\text{colony}}$, which is the sum of the doses corresponding to the two different serotype classes of MPS-STEC. The probability of getting the HUS disease by consuming a serving of 25 g of cheese, conditional on Age = age and dose $\Gamma_{\text{MPS-STEC}} = \gamma_{\text{MPS-STEC}}$, can be written as

$$P_{\text{HUS}}(\text{age}, \gamma_{\text{MPS-STEC}}) = 1 - (1 - r_0 \cdot \exp(-k \cdot \text{age}))^{\gamma_{\text{MPS-STEC}}}. \quad (11)$$

For *Salmonella*, a Beta-Poisson dose-response model proposed by McCullough and Elsele (1951), World Health Organization (2002) is used, with model parameters (Strickland et al., 2023) dependent on the two subclass serotypes $s$, for $s \in \{\text{HV Salmonella, LV Salmonella}\}$,

$$P_{\text{Salmonellosis}}(\Gamma_s) = 1 - \left(1 + \frac{\Gamma_s}{\beta_s}\right)^{-\alpha_s}. \quad (12)$$

Conditional on dose $\Gamma_s = N_s^{\text{colony-serving}} \cdot Y_s^{\text{colony}}$, for each serotype $s$, the probability of getting salmonellosis from the consumption of a serving of 25 g of cheese, is a weighted average, using the respective consumption probabilities based on expected number of high and low virulent colonies in a cheese serving:

$$\begin{aligned} P_{\text{Salmonellosis}}(\Gamma_{\text{HV-Salmo}}, \Gamma_{\text{LV-Salmo}}) &= \\ P_{\text{Salmonellosis}}(\Gamma_{\text{HV-Salmo}}) \cdot P_{\text{HV-Salmo}}^{\text{colony}} + P_{\text{Salmonellosis}}(\Gamma_{\text{LV-Salmo}}) \cdot P_{\text{LV-Salmo}}^{\text{colony}} \end{aligned}$$

For *Lm*, the dose-response is adapted from Ricci et al. (2018), based on the Poisson model, which takes into account the variability in susceptibility across mutually exclusive population subgroups, as proposed by Pouillot et al. (2015). Conditional on dose $\Gamma_{Lm} = N_{Lm}^{\text{colony-serving}} \cdot Y_{Lm}^{\text{colony}}$, and population subgroup $\theta_g$, the probability of getting listeriosis, from the consumption of a serving of 25 g of cheese, is written as

$$P_{\text{Listeriosis}}(\theta_g, \Gamma_{\text{Lm}}) = \int_0^1 (1 - \exp(-r \cdot \Gamma_{\text{Lm}})) p(r; \theta_g) \, dr. \quad (13)$$

where $r$ is the probability of developing listeriosis from the ingestion of a bacterial cell in a given, specific serving, and $p(r; \theta_g)$ represents the remaining individual (within-group) susceptibility variability and strain virulence variability in $r$.

**Batch risk** The consumer module models the risk for each of the food-borne illnesses, associated with the consumption of a particular batch of cheese. For each pathogen $x$, the batch risk depends on a vector of stochastic internal variables of the QMRA model, denoted by $\Xi_x$, which characterizes a batch.

The vector $\Xi_x$ includes the initial concentration $Y_x^{\text{milk}}$ of the respective pathogen in the ATM, the stochastic parameters $\{d^{\text{storage}}, T^{\text{storage}}\}$, corresponding to the storage step, and the stochastic consumption time $t^{\text{consum}}$. Depending on the specific pathogen, $\Xi_x$ comprises of other stochastic variables. For *Salmonella*, $\Xi_x$ includes the rate parameters $\{\delta_{\text{core}}, \delta_{\text{rind}}\}$ and for *Lm* it includes the rate parameter $\rho_{Lm}$.

We define the batch risk $R_x^{\text{batch}}(\xi_x)$ as the probability of getting the particular illness from pathogen $x$, by consuming a portion of 25 g of cheese from a particular batch, characterized by the internal variables $\Xi_x = \xi_x$. Further we define two events $E_{\text{illness}(x)} \subset E_{\text{market}}$, where $E_{\text{market}}$ indicates that the particular batch goes to the market (i.e. is not rejected)
<table><caption>**Table 8**
Inputs for postharvest module collectively denoted as *o*post are shown in the first part of the table. In the second part the cost values of the intervention steps are shown.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Values</th></tr></thead><tbody><tr><td>n<sup>sample</sup></td><td>Number of test portions</td><td>5</td></tr><tr><td>m<sup>sample</sup></td><td>Mass of each test portion</td><td>25 gm</td></tr><tr><td>p<sup>test</sup></td><td>Prop. of batch tested</td><td>0.5</td></tr><tr><td>t<sup>test</sup></td><td>Time when batch is tested</td><td>14th Day</td></tr><tr><td>C<sub>milk</sub><sup>test</sup></td><td>Cost of testing farm milk</td><td>10 EUR</td></tr><tr><td>C<sub>milk</sub><sup>loss</sup></td><td>Cost of rejecting one Liter milk</td><td>0.2 EUR</td></tr><tr><td>C<sub>cheese</sub><sup>test</sup></td><td>Cost of testing one cheese sample</td><td>70 EUR</td></tr><tr><td>C<sub>cheese</sub><sup>loss</sup></td><td>Cost of rejecting one cheese</td><td>1.5 EUR</td></tr></tbody></table>

and $E_{illness(x)}$ stands for the event that the particular batch goes to the market and the consumer gets the illness from pathogen x. The mutual dependence between different model variables and these events are shown in Fig. 8. In this Directed Acyclic Graph, a specific dependent variable, conditioned by any of the variables placed at one of its arrow ends, is considered to be conditionally independent with respect to all the variables positioned at the subsequent arrow ends. From the definition of the batch risk we have

$$ 
\begin{aligned}
R_x^{\text{batch}}(\xi_x) &= P[E_{\text{illness}(x)} | \Xi_x = \xi_x, E_{\text{market}}] \\
&= \mathbb{E}[1_{E_{\text{illness}(x)}} | \Xi_x = \xi_x, 1_{E_{\text{market}}} = 1],
\end{aligned} 
\quad (14) $$

and from the definition of the probability of illness in (11)–(13), we have,

$$ 
\begin{aligned}
P_{\text{illness}(x)}(\text{age}, \gamma_x) &= P[E_{\text{illness}(x)} | \text{Age} = \text{age}, \Gamma_x = \gamma_x, E_{\text{market}}] \\
&= \mathbb{E}[1_{E_{\text{illness}(x)}} | \text{Age} = \text{age}, \Gamma_x = \gamma_x, 1_{E_{\text{market}}} = 1].
\end{aligned} 
\quad (15) $$

Since $1_{E_{illness(x)}}$ is conditionally independent of $\Xi_x$, given $\Gamma_x$, using (15) and the law of total expectation, we can write

$$ 
\begin{aligned}
\mathbb{E}[1_{E_{\text{illness}(x)}} | \Xi_x, 1_{E_{\text{market}}}] \\
&= \mathbb{E}[\mathbb{E}[1_{E_{\text{illness}(x)}} | \text{Age}, \Gamma_x, \Xi_x, 1_{E_{\text{market}}}] | \Xi_x, 1_{E_{\text{market}}}] \\
&= \mathbb{E}[\mathbb{E}[1_{E_{\text{illness}(x)}} | \text{Age}, \Gamma_x, 1_{E_{\text{market}}}] | \Xi_x, 1_{E_{\text{market}}}] \\
&= \mathbb{E}[P_{\text{illness}(x)}(\text{Age}, \Gamma_x) \cdot 1_{E_{\text{market}}} + 0 \cdot \bar{1}_{E_{\text{market}}} | \Xi_x, 1_{E_{\text{market}}}] .
\end{aligned} 
\quad (16) $$

Using (16) and the conditional independence with respect to $E_{market}$ given $\Xi_x$, the definition of $R_x^{\text{batch}}(\xi_x)$ in (14) boils down to

$$ 
\begin{aligned}
R_x^{\text{batch}}(\xi_x) &= \mathbb{E}[P_{\text{illness}(x)}(\text{Age}, \Gamma_x) | \Xi_x = \xi_x] \\
&= \sum_{\text{age}=1}^{\text{age max}} \int_{0}^{\infty} P_{\text{illness}(x)}(\text{age}, \gamma_x) \cdot p(\text{age}, \gamma_x | \Xi_x = \xi_x) d\gamma_x \\
&= \sum_{\text{age}=1}^{\text{age max}} g(\text{age}) \int_{0}^{\infty} P_{\text{illness}(x)}(\text{age}, \gamma_x) p(\gamma_x | \xi_x) d\gamma_x,
\end{aligned} 
\quad (17) $$

where the joint probability distribution of the Age and the dose $\Gamma_x|\Xi_x = \xi_x$ is $p(\text{age}, \gamma_x|\Xi_x = \xi_x) = g(\text{age}) \cdot p(\gamma_x|\xi_x)$, with $g(\text{age}) = P[\text{Age} = \text{age}|E_{\text{market}}]$ being the age distribution of cheese consumers.

For MPS-STEC, the parameter $g(\text{age})$ controls the proportion of cheese consumed by the age group (see, e.g., Perrin et al., 2014), and for *Salmonella*, the age parameter is considered to have no effects (Teunis et al., 2010; Teunis, 2022). The batch risk for *Lm* uses the model by Cadavez et al. (unpublished), where the population subgroups are based on different age groups, namely, $\{[1,4], [5,14], [15,24], [25,44], [45,64], [65,74], [75, +\infty]\}$, for males and females separately.

## 2.5. Risk management options

### 2.5.1. Preharvest: milk sorting

The preharvest intervention strategy, a.k.a milk sorting, is carried

out just before mixing the BTM from several farms into the ATM. In this step, the *E. coli* levels in BTM coming from the farms are tested and tanks with concentration above a certain threshold are rejected from the production chain. Each BTM is tested for *E. coli* concentration and the ith farm is rejected if $Y_i^{EC} > l^{\text{sorting}}$. The milk loss for a particular batch is given by $M^{\text{batch}} = \sum_{i=1}^{N^{\text{farms}}} V_i 1_{\{i \in S\}}$, where S denotes the set of farms that qualify after milk sorting. We note that the parameters controlling the preharvest intervention step are mainly the frequency of milk sorting $f^{\text{sorting}}$ (measured in days) and the threshold of milk sorting $l^{\text{sorting}}$ (measured in CFU/mL), as listed in Table 3. The milk sorting strategy is based on the hypothesis that *E. coli* and other pathogen strains follow the same fecal routes in the cows body (see, e.g. Perrin et al., 2014), as a results the *E. coli* concentration in the BTM, can be used as a measure of the farms hygiene conditions.

**Remark 2.** In the current implementation of the QMRA model this milk sorting strategy only affects the concentration of STEC, MPS-STECA and *Salmonella* in the ATM. For *Lm* since the concentration in the ATM is simulated directly, it remains unaffected by the milk sorting. Implementation of a more realistic and efficient preharvest intervention strategy is discussed in Section 5.1.

### 2.5.2. Postharvest: microbial cheese sampling

The postharvest intervention strategy, a.k.a cheese sampling, can be implemented at different stages of cheese production depending on the type of bacteria. Typically, for raw milk soft cheese, the sampling process is carried out at the end of the salting step, during the third day of production. However, some producers may choose to implement cheese sampling during cheese ripening, until the 14th day from the beginning of production. In our model, the cheese sampling time is determined by a parameter $l^{\text{test}}$, which is by default set at 14th day of production. Cheese sampling consists in inspecting a batch of cheese for pathogen contamination, and this is done with respect to small portions, called sampling units, taken out of the batch. Once a single sample unit is tested positive for any of the three major pathogens under consideration, the whole batch of cheese is not sent to the market. Let $n^{\text{sample}}$ and $m^{\text{sample}}$ respectively denote the number of sample units taken from a batch and the mass of each sample unit (usually fixed at 25 g). The probability for a sample unit being tested positive for the three main pathogens, is

$$ P_x^{\text{unit}}(\xi_x) = P(N_x^{\text{colony,sample}} > 0 | \Xi_x = \xi_x) $$

where conditionally on $\Xi_x = \xi_x$, $N_x^{\text{colony,sample}}$ is a Poisson distributed random variable, with expectation $(\lambda_x^{\text{colony}} n_x^{\text{sample}} / w_x^{\text{cheese}})$, and $\lambda_x^{\text{colony}}$ corresponds to the average number of colonies for pathogen x in the cheese sample, given $\Xi_x = \xi_x$.

For MPS-STECA and *Salmonella*, the average number of colonies $\lambda_x^{\text{colony}}$ is obtained as the sum of the average number of colonies of their respective subclass serotypes. The probability of detecting a colony in any of the sample units is $P_x^{\text{sample}}(\xi_x) = 1 - (1 - P_x^{\text{unit}}(\xi_x))^{n^{\text{sample}}}$. Assuming the detection events of the three types of pathogens are independent, the probability of rejecting the batch of cheese, given $\Xi_x = \xi_x$, where $\Xi_x = (\Xi_x)_x$, is obtained using the formula for the union of three independent events:

$$ P^{\text{batch}}(\xi) = \sum_{x} P_{x}^{\text{sample}}(\xi_x) - \sum_{x<\xi} P_{x}^{\text{sample}}(\xi_x) \cdot P_{x'}^{\text{sample}}(\xi_x) + \prod_{x} P_{x}^{\text{sample}}(\xi_x). \quad (18) $$

The parameters for the postharvest step are listed in Table 8.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Output module denoted by the green box is the next hierarchical module of the batch level simulator. It simulates N<sup>batch</sup> independent batches to estimate the average risk of illness R<sub>x</sub> from pathogen x, average milk loss M<sup>avg</sup> and average probability P<sup>avg</sup> of rejecting a cheese batch.</figcaption></figure>

## 2.6. Risk characterization

### 2.6.1. Output module

The output module, demonstrated in Fig. 9, is outside the batch level simulator and its purpose is to obtain estimations of the expectation of the quantities of interest for a given set of $\theta = \{\theta^{\text{farm}}, \theta^{\text{cheese}}, \theta^{\text{con}}, \theta^{\text{post}}\}$. As will be explained in Section 3, the estimations will be carried out using Monte Carlo simulations. The output module uses three batch level simulator outputs, namely, the batch risk $R_x^{\text{batch}}(\xi_x)$, the probability of rejecting a batch $P_x^{\text{batch}}(\xi)$ and the milk loss $M_x^{\text{batch}}$.

The average milk loss is defined as $M_x^{\text{avg}} = \mathbb{E}[M_x^{\text{batch}}]$.

The average batch rejection rate, or equivalently the average probability of a batch not going to the market is defined as $P_x^{\text{avg}} = P[\bar{E}_{\text{market}}]$. Further we have $\mathbb{E}[1_{\bar{E}_{\text{market}}} | \Xi = \xi] = p^{\text{test}} p_x^{\text{batch}}(\xi)$, where $p^{\text{test}}$ is the proportion of cheese batches tested. Then using the law of total expectation we can write

$$ 
\begin{aligned}
P_x^{\text{avg}} &= \mathbb{E}[\mathbb{E}[1_{\bar{E}_{\text{market}}} | \Xi]] \\
&= \mathbb{E}[p^{\text{test}} p_x^{\text{batch}}(\Xi)],
\end{aligned}
\quad (19) $$

The average risk of illness R<sub>x</sub>, from pathogen x, is the conditional probability P[E<sub>illness(x)</sub> | E<sub>market</sub>], of getting the illness from consuming a portion of 25g of cheese from a batch of cheese, that was not rejected, produced with input parameters $\theta$. Recalling $\Xi = (\Xi_x)_x$ and using the law of total expectation twice, we can write

$$ 
\begin{aligned}
P[E_{\text{illness}(x)}] &= \mathbb{E}[1_{E_{\text{illness}(x)}}] \\
&= \mathbb{E}[\mathbb{E}[1_{E_{\text{illness}(x)}} | \Xi_x, 1_{E_{\text{market}}}]] \\
&= \mathbb{E}[R_x^{\text{batch}}(\Xi_x) \cdot 1_{E_{\text{market}}} + 0 \cdot 1_{\bar{E}_{\text{market}}}] \\
&= \mathbb{E}[R_x^{\text{batch}}(\Xi_x) \cdot \mathbb{E}[1_{E_{\text{market}}} | \Xi]] \\
&= \mathbb{E}[R_x^{\text{batch}}(\Xi_x) \cdot (1 - p^{\text{test}} p_x^{\text{batch}}(\Xi))].
\end{aligned}
\quad (20) $$

Now using (20) and recalling $E_{\text{illness}(x)} \subset E_{\text{market}}$, the definition of R<sub>x</sub> boils down to,

$$ 
\begin{aligned}
R_x &= \frac{P[E_{\text{illness}(x)}]}{P[E_{\text{market}}]} \\
&= \frac{\mathbb{E}[R_x^{\text{batch}}(\Xi_x) \cdot (1 - p^{\text{test}} p_x^{\text{batch}}(\Xi))]}{1 - P^{\text{avg}}}. 
\end{aligned}
\quad (21) $$

**Module outputs** The three quantities of interest produces by the output module are, the average milk loss M<sup>avg</sup>, average batch rejection rate P<sup>avg</sup> and the average risk of illness R<sub>x</sub> from pathogen x. These quantities of interest characterizes a cheese production process with input parameters $\theta$.

### 2.6.2. DALY: disability-adjusted life years

A key feature of our multipathogen QMRA model is its ability to assess the impact of pathogens on public health due to the consumption of raw milk soft cheese. Within all microbial hazards potentially linked to soft cheese (Dubois-Brissonnet et al., 2022), we focus on the most severe pathogens, and also to evaluate their collective impact on public health, specifically in terms of the years of healthy life lost due to the illnesses resulting from cheese consumption.

To assess the impact on public health from the potential illnesses induced by cheese consumption, we adopt an approach based on DALYs (Murray and Lopez, 1997), an acronym for Disability-Adjusted Life Years. DALYs combine the years of life lived with a disability (YLD) and the years of life lost (YLL) due to premature death from the illness.

The DALY takes into account a spectrum of symptoms, which include death, loss of kidney function, prolonged or short-term hospitalization, bloody diarrhea, diarrhea, among others. Furthermore, the age of the afflicted person impacts DALYs, with a child's death yielding higher DALYs compared to the death of an older individual.

**DALY metric for cheese portions** First we define the expected DALY due to illness(x) caused by the consumption of a cheese portion by a consumer of Age = age,

$$ \bar{\text{DALY}}_{\text{portion illness}(x)}(\text{age}) = \mathbb{E}[\text{DALY}_x | 1_{E_{\text{illnessx}}} = 1, \text{Age} = \text{age}], \quad (22) $$

where DALY<sub>x</sub> is a random variable, denoting the DALY caused by the consumption of pathogen x. Further we define the age specific average risk R<sub>x</sub>(age), from pathogen x, which is the conditional probability P[E<sub>illness(x)</sub> | Age = age, 1<sub>E<sub>market</sub></sub> = 1], for consumers with a particular age, of getting the illness from consuming a portion of 25g of cheese from a batch of cheese, produced with input parameters $\theta$ and given that the particular batch actually goes into the market (not rejected). Following similar derivation steps of R<sub>x</sub>, we can write

$$ R_x(\text{age}) = \frac{\mathbb{E}[R_x^{\text{batch}}(\Xi_x, \text{age}) \cdot (1 - p^{\text{test}} p_x^{\text{batch}}(\Xi)) | \text{Age} = \text{age}]}{1 - P^{\text{avg}}}, \quad (23) $$

where the age specific batch risk $R_x^{\text{batch}}(\text{age})$ is given by

$$ R_x^{\text{batch}}(\xi_x, \text{age}) = \int_0^\infty P_{\text{illness}(x)}(\gamma_x, \text{age}) \cdot p(\gamma_x | \Xi_x = \xi_x) \, d\gamma_x. \quad (24) $$

Then using the law of total expectation, the definition of overall risk in (21) can be expressed as $R_x = \mathbb{E}[R_x(\text{age})]$.

Now for the assessment of impact on public health, we are interested in the average DALY due to the consumption of cheese portions contaminated with pathogen x, from a batch of cheese that went to the market. This metric is denoted by $\bar{\text{DALY}}_{\text{portion},x}$, is derived using the law of total expectations.
<figure><img src="image_13.png" /><figcaption>Fig. 10. Estimated values of hygiene parameters, namely α and σ², using the data collected during 2019 – 2022 from all the 97 farms (which had no missing data, out of 104 farms) under three different cheese producers in France.</figcaption></figure>

$$
\begin{align*}
\overline{\text{DALY}}_{\text{portion},x} &= \mathbb{E}[\text{DALY}_x | 1_{E_{\text{market}}} = 1] \\
&= \mathbb{E}[\mathbb{E}[\text{DALY}_x | \text{Age}, 1_{E_{\text{illness}(x)}}, 1_{E_{\text{market}}}] | 1_{E_{\text{market}}} = 1] \\
&= \mathbb{E}[\mathbb{E}[\text{DALY}_x | \text{Age}, 1_{E_{\text{illness}(x)}}] | 1_{E_{\text{market}}} = 1] \\
&= \mathbb{E}[\overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}} (\text{Age}) \cdot 1_{E_{\text{illness}(x)}} + 0 \cdot \overline{1}_{E_{\text{illness}(x)}} | 1_{E_{\text{market}}} = 1] \\
&= \mathbb{E}[\overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}} (\text{Age}) \cdot 1_{E_{\text{illness}(x)}} | 1_{E_{\text{market}}} = 1]. \tag{25}
\end{align*}
$$

The expectation in (25) can be expressed as a weighted sum of
$\overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age})$, with respect to $R_x(\text{age})$ and $g(\text{age})$,

$$
\begin{align}
\overline{\text{DALY}}_{\text{portion},x} &= \sum_{\text{age}=1}^{\text{age}^{\text{max}}} \overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age}) \cdot \nonumber \\
&\quad P[E_{\text{illness}(x)} | \text{Age} = \text{age}, 1_{E_{\text{market}}} = 1] \cdot P[\text{Age} = \text{age} | 1_{E_{\text{market}}} = 1] \nonumber \\
&= \sum_{\text{age}=1}^{\text{age}^{\text{max}}} \overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age}) \cdot R_x(\text{age}) \cdot g(\text{age}). \tag{26}
\end{align}
$$

The expectation in (25) can be expressed equivalently as

$$
\begin{align}
\overline{\text{DALY}}_{\text{portion},x} &= \sum_{\text{age}=1}^{\text{age}^{\text{max}}} \overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age}) \cdot \nonumber \\
&\quad P[\text{Age} = \text{age} | E_{\text{illness}(x)}] \cdot P[E_{\text{illness}(x)} | E_{\text{market}}] \nonumber \\
&= \sum_{\text{age}=1}^{\text{age}^{\text{max}}} \overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age}) \cdot \tilde{g}_{\theta}(\text{age}) \cdot R_x \nonumber \\
&= R_x \cdot \sum_{\text{age}=1}^{\text{age}^{\text{max}}} \overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age}) \cdot \tilde{g}_{\theta}(\text{age}), \tag{27}
\end{align}
$$

where $\tilde{g}_{\theta}(\text{age}) = P[\text{Age} = \text{age} | E_{\text{illness}(x)}]$ is the age distribution for the cases of illness caused by pathogen x, given the simulator inputs $\theta$.

DALY estimation We note that the estimation of the average DALY metric $\overline{\text{DALY}}_{\text{portion},x}$ using (26) or (27) includes the age specific DALY values $\overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age})$, which are not directly available in the literature. However, (27) can be simplified further given certain assumptions.

The age distribution $\tilde{g}_{\theta}(\text{age})$ for the illness($x$), is dependent on the QMRA simulator inputs $\theta$, through the corresponding dose-response model and is directly related to the proportion of cheese consumption $g(\text{age})$ for different age group. Acknowledging that our dose-response models heavily rely on epidemiological data primarily sourced from food-borne illness outbreaks, we tentatively posit that the age distribution $\tilde{g}_{\theta}(\text{age})$ specific to each illness remains relatively impervious to the influences of simulator inputs $\theta$, for reasonable variations around their nominal (baseline) values. Under this assumption, the age distribution generated by our QMRA raw milk cheese simulator can be assumed to be closely aligned with the global pattern. Although for the purposes of this study we make the above assumption, it is important to acknowledge

that this hypothesis is made for simplification and modeling conve-
nience. The validity of this hypothesis in the real-world context remains
an open question and deserves further investigation.

For the calculation of the DALY metric, we reside to the study by
Cassini et al. (2018), which is based on the Burden of Communicable
Diseases in Europe (BCoDE) project (Kretzschmar et al., 2012; Mangen
et al., 2013), focused on the EU/EEA population between 2009 and
2013. In this study, for each of the concerned illnesses, a model was
created using the BCoDE toolkit (European, 2019). Within each model,
age group-specific and sex-specific annual case numbers, multiplication
factors to account for underestimation, and population data were input
into the software. Finally a Monte Carlo simulation was performed with
10,000 iterations for each illness, to estimate the median DALYs per
case. Table 2 shows the estimated median DALY corresponding to a
singular case, represented as DALY(1 case)$_x$ along with the two com-
ponents, namely, YLL and YLD for 1000 cases, for each associated
illness. Given our assumptions, the average DALY metric in (27) can be
simplified by substituting $\sum_{\text{age}=1}^{\text{age}^{\text{max}}} \overline{\text{DALY}}_{\text{portion illness}(x)}^{\text{portion}}(\text{age}) \cdot \tilde{g}_{\theta}(\text{age})$ by
DALY(1 case)$_x$ from Table 2,

$$
\overline{\text{DALY}}_{\text{portion},x} = \text{DALY}(1 \text{ case})_x \times R_x. \quad (28)
$$

It is important to highlight that, in our study, the values presented in Table 2 are used as best available approximations from the literature. This choice is motivated by the absence of pertinent analyses from epidemiological studies in the existing literature. Moreover, it serves the purpose of streamlining the computations within our study. However, it is crucial to recognize that utilizing global DALY values for assessing the health impacts of cheese consumption in France is not ideal. This approach relies on substantial assumptions, as elaborated earlier. Consequently, the computation of the DALY metric in our study does not fully capture the precise and prevailing impact on public health attributable to cheese consumption. Therefore, these results should be interpreted judiciously. Recognizing the necessity for greater accuracy, the computation of a more refined DALY metric incorporating relevant epidemiological studies is deferred as a prospective avenue for future research.

Moreover, our model considers the dose-response from Perrin et al. (2014), which estimates the risk of HUS, but the DALY values that could be found in the literature (see, e.g., Table 2) are either based on all symptoms of STEC or MPS-STEC. The investigation conducted by Lindqvist et al. (2023) identifies more predominant strains of MPS-STEC that specifically contribute to HUS, which is our primary focus. Therefore, a more comprehensive epidemiological study is essential to accurately estimate DALYs attributed solely to HUS. However, in this study we use the DALYs corresponding to MPS-STEC to compute the metric in (28).

The combined influence of the three pathogens, corresponding to the
aggregated impact on public health from the consumption of raw milk
cheese, is defined as $\overline{\text{DALY}}_{\text{portion}} = \sum_{x} \overline{\text{DALY}}_{\text{portion},x}$, ignoring the con-
sequences of concurrent instances of two or more illnesses.
<figure><img src="image_14.png" /><figcaption>Fig. 11. Scatterplot of the estimated values of $\alpha$ and $\sigma^2$.</figcaption></figure>

# 3. Model implementation

## 3.1. Mathematical and statistical techniques

### 3.1.1. Estimation of farm hygiene parameters

The hygiene level of each of the farms is characterized by two parameters, namely $\alpha$ and $\sigma$, which control the distribution of E. coli in the farm milk. In Perrin et al. (2014), the authors proposed a hierarchical Poisson mixed model to express the relationship between these parameters and the daily E. coli concentration $x_d$ in BTM, that can be written as

$$ 
\begin{aligned}
x_d | \Lambda_d &\sim \text{Poisson}(\Lambda_d) \\
\log(\Lambda_d) &= \alpha + \varepsilon_d, \text{ where } \varepsilon_d \sim \mathcal{N}(0, \sigma^2),
\end{aligned} 
\quad (29) 
$$

where $\Lambda_d$ denotes the average E. coli concentration in BTM of the farm on day $d$.

We propose a Bayesian approach based on Markov Chain Monte Carlo (MCMC) sampling to estimate the parameters of this model, separately for each farm. E. coli test data was collected by CNIEL and ACTALIA from three French cheese producers, covering a total of 104 farms. The collected data ranges from 10 CFU/ml (which corresponds to the detection limit of the method) to either 150 or 300 CFU/ml (upper quantification limit) depending on the producers.

Fig. 10 plots two histograms of the values of $\alpha$ and $\sigma^2$ respectively, as estimated for all the farms. Fig. 11 displays the pairs of values for the two hygiene parameters across different farms, showing a tendency for farms with higher $\alpha$ values to also have higher $\sigma^2$ values.

### 3.1.2. Computing $R^{\text{batch}}(\xi_x)$ using simple Monte Carlo

The computation of batch risk for Lm is addressed by the JEMRA (Cadavez et al., unpublished) toolbox in R, whereas for MPS-STECE and Salmonella one possible approach is based on the simple Monte Carlo integration method. It involves the computation of the integral in (17), with respect to the dose $\Gamma_x$ conditional on the stochastic parameters $\xi_x$. This integral can be approximated using $N^{\text{dose}}$ i.i.d samples $\{\gamma_{1,x}, \gamma_{2,x}, \dots, \gamma_{N^{\text{dose}},x}\}$ from the conditional distribution $p(\Gamma_x | \Xi_x = \xi_x)$, of the random variable dose $\Gamma_x$, given the vector of stochastic internal variables $\xi_x$.

The approximated batch risk conditional on $\xi_x$ is obtained as,

$$ \hat{R}_x^{\text{batch}}(\xi_x) = \sum_{\text{age}=1}^{15} g(\text{age}) \frac{1}{N^{\text{dose}}} \sum_{i=1}^{N^{\text{dose}}} P_{\text{illness}(x)}(\text{age}, \gamma_{i,x}), \quad (30) $$

where $P_{\text{illness}(x)}(\text{age}, \gamma_{i,x})$ is the probability of getting the illness from pathogen $x$, for the consumers of a particular age group, by consuming a portion of 25g of cheese with dose $\gamma_{i,x}$, coming from a batch associated with internal variables $\xi_x$.

### 3.1.3. Computing $R^{\text{batch}}(\xi_{\text{MPS-STECE}})$ using integral approximation

For MPS-STECE, the batch risk in (17) can be expressed as

$$ R_{\text{MPS-STECE}}^{\text{batch}}(\xi_{\text{MPS-STECE}}) = \sum_{\text{age}=1}^{\text{age max}} g(\text{age}) \mathbb{E}[(1 - r_a)^{\Gamma_{\text{MPS-STECE}}}], \quad (31) $$

where $r_a = r_0 \cdot \exp(k \cdot \text{age})$. The random variable dose $\Gamma_{\text{MPS-STECE}}$ can be decomposed additively with respect to the two classes of MPS-STECE strains, $\Gamma_{\text{MPS-STECE}} = \sum_s \Gamma_s$ for $s \in \{\text{MPS O157:H7}, \text{MPS nonO157:H7}\}$, with

$$ 
\begin{aligned}
\Gamma_s &= Y_s^{\text{consum}} N_s^{\text{colony}} 10^{r_s \epsilon} \\
&= d_s N_s^{\text{colony}} b_s^{\epsilon},
\end{aligned} 
\quad (32) 
$$

where $d_s = Y_s^{\text{consum}}$, $b_s = 10^{r_s \epsilon}$, and $\epsilon$ is a standard normal variable. Now assuming the probability of getting HUS from each of the strains of MPS-STECE is independent, the expectation in (31) can be decomposed using the law of total expectations,

$$ 
\begin{aligned}
\mathbb{E}[(1 - r_a)^{\Gamma_{\text{MPS-STECE}}}] &= \prod_s \mathbb{E} \left[ \mathbb{E}[(1 - r_a)^{\Gamma_s} | \epsilon] \right] \\
&= \prod_s \mathbb{E} \left[ \mathbb{E} \left[ (1 - r_a)^{d_s N_s^{\text{colony}} b_s^{\epsilon}} | \epsilon \right] \right].
\end{aligned} 
\quad (33) 
$$

The inner expectations $\mathbb{E}[(1 - r_a)^{\Gamma_s} | \epsilon]$, can be analytically derived,

$$ 
\begin{aligned}
\mathbb{E}[(1 - r_a)^{N_s^{\text{colony}} d_s b_s^{\epsilon}} | \epsilon] &= \sum_{N_s^{\text{colony}}=0}^{\infty} (1 - r_a)^{N_s^{\text{colony}} d_s b_s^{\epsilon}} \exp(-\lambda_s^{\text{colony}}) \frac{(\lambda_s^{\text{colony}})^{N_s^{\text{colony}}}}{N_s^{\text{colony}}!} \\
&= \sum_{N_s^{\text{colony}}=0}^{\infty} \frac{\{ (1 - r_a)^{d_s b_s^{\epsilon}} \lambda_s^{\text{colony}} \}^{N_s^{\text{colony}}} \exp\left(-\left\{ (1 - r_a)^{d_s b_s^{\epsilon}} \lambda_s^{\text{colony}} \right\}\right)}{N_s^{\text{colony}}! \exp\left(-\left\{ (1 - r_a)^{d_s b_s^{\epsilon}} \lambda_s^{\text{colony}} \right\}\right)} \exp(-\lambda_s^{\text{colony}}) \\
&= \exp\left(-\lambda_s^{\text{colony}} \{ (1 - r_a)^{d_s b_s^{\epsilon}} \} \right).
\end{aligned} 
\quad (34) 
$$

Substituting (34) in (33), we get,

$$ \mathbb{E}[(1 - r_a)^{\Gamma_{\text{MPS-STECE}}}] = \exp\left(-\left(\lambda_{0157:H7}^{\text{colony}} + \lambda_{0157:H7}^{\text{colony}}\right)\right) \prod_s \mathbb{E}\left[c_{s,1}^{c_{s,1}^{\text{colony}}}\right], \quad (35) $$

where $c_{s,1} = \exp(\lambda_s^{\text{colony}})$, $c_{s,2} = (1 - r_a)^{d_s}$ and $c_{s,3} = b_s$.

Given that $c_{s,1} > 1$, $c_{s,2} < 1$, and $c_{s,3} > 1$, the function $c_{s,1}^{c_{s,2}^{\text{colony}}}$ is monotonically non-increasing with respect to $\epsilon$. Such functions can be integrated using deterministic quadrature methods (e.g., the trapezoidal rule), which offer a better convergence rate compared to simple Monte Carlo (see, for example, Basak et al., 2022).

### 3.1.4. Computation of simulator outputs

The final outputs of the simulator as mentioned in Section 2.6.1 are estimated using simple Monte Carlo method, by simulating $N^{\text{batch}}$ independent batches. The average milk loss $M^{\text{avg}} = \mathbb{E}[M^{\text{batch}}]$, is approximated using the average of $M^{\text{batch}}$ from each batch,

$$ \hat{M}^{\text{avg}} = \frac{1}{N^{\text{batch}}} \sum_{i=1}^{N^{\text{batch}}} M_i^{\text{batch}} \quad (36) $$

Each of the simulated batches yields $\{\xi_1, \xi_2, \dots, \xi_{N^{\text{batch}}}\}$, the set of internal stochastic parameters which are used to construct the unbiased estimated of the quantities of interest:
<caption>**Table 4**</caption>

Distribution of number of cows: data provided by CNIEL and ACTALIA, collected from 31 producers of milk. This empirical distribution is used to simulate the number of cows in farm in the QMRA model.

<table><thead><tr><th>Cows</th><th>5-20</th><th>20-40</th><th>40-60</th><th>60-80</th><th>80-100</th><th>100-120</th><th>Total</th></tr></thead><tbody><tr><th>Farms</th><td>1</td><td>2</td><td>18</td><td>8</td><td>1</td><td>1</td><td>31</td></tr></tbody></table>

$$ \hat{P}^{\text{avg}} = \frac{1}{N^{\text{batch}}} \sum_{l=1}^{N^{\text{batch}}} \hat{P}^{\text{batch}}(\xi_l) p^{\text{test}}. \qquad (37) $$

$$ \hat{R}_x = \frac{1}{N^{\text{batch}}(1 - \hat{P}^{\text{avg}})} \sum_{l=1}^{N^{\text{batch}}} \hat{R}_x^{\text{batch}}(\xi_l) (1 - \hat{P}^{\text{batch}}(\xi_l) p^{\text{test}}) \qquad (38) $$

### 3.1.5. Computation of intervention cost

One of the several objectives of the ArtiSaneFood project includes the implementation of intervention strategies for controlling the risk of food-borne illness. In this case the cost of the intervention strategies plays an important role. As detailed in Sections 2.5.1 and 2.5.2, the cheese production process involves two types of intervention strategies: one involving the testing of farm milk, and the other involving the testing of cheese after production. In addition to the expenses related to conducting the microbiological tests, these intervention steps involve the rejection of both farm milk and cheese batches. The QMRA simulator produces two outputs concerning the average loss of milk per batch $M^{\text{avg}}$, and the average probability of rejecting a batch of cheese $P^{\text{avg}}$. Thus the total cost of the intervention steps (C), can be decomposed into two part corresponding to the preharvest and postharvest step, as $C = C^{\text{pre}} + C^{\text{post}}$.

$$ C^{\text{pre}} = N^{\text{farm}} C_{\text{test}}^{\text{milk}} (1/f^{\text{sorting}}) + M^{\text{avg}} C_{\text{loss}}^{\text{milk}} \qquad (39) $$

$$ C^{\text{post}} = n^{\text{sample}} C_{\text{test}}^{\text{cheese}} p^{\text{test}} + p^{\text{avg}} p^{\text{test}} C_{\text{loss}}^{\text{cheese}} N^{\text{cheese}} $$

The parameters $\{C_{\text{test}}^{\text{milk}}, C_{\text{loss}}^{\text{milk}}, C_{\text{test}}^{\text{cheese}}, C_{\text{loss}}^{\text{cheese}}\}$ corresponding to different intervention costs, are provided French cheese producer representatives and $N^{\text{cheese}}$ denotes the average number of cheeses present in a batch. All the parameters in (39) are described in Table 8.

<caption>**Table 9**</caption>

Baseline results of QRMA model: Summary statistics (mean, median, standard deviation and quantiles of order 5% and 95%) of the simulated concentration (in $\log_{10}$ scale) for STEC, MPS-STECA, *Salmonella* and *Lm*, in ATM with no milk testing intervention step.

<table><thead><tr><th>Symbol</th><th>Mean</th><th>Median</th><th>SD</th><th>q<sub>0.05</sub></th><th>q<sub>0.95</sub></th></tr></thead><tbody><tr><td>Y<sub>STEC</sub><sup>milk</sup></td><td>-3.44</td><td>-3.5</td><td>0.58</td><td>-4.28</td><td>-2.39</td></tr><tr><td>Y<sub>MPS-STECA</sub><sup>milk</sup></td><td>-5.91</td><td>-5.94</td><td>1.0</td><td>-7.51</td><td>-4.24</td></tr><tr><td>Y<sub>LV-Salmo</sub><sup>milk</sup></td><td>-5.78</td><td>-5.86</td><td>0.94</td><td>-7.21</td><td>-4.12</td></tr><tr><td>Y<sub>LV-Salmo</sub><sup>milk</sup></td><td>-5.2</td><td>-5.28</td><td>0.82</td><td>-6.41</td><td>-3.72</td></tr><tr><td>Y<sub>Lm</sub><sup>milk</sup></td><td>-3.12</td><td>-3.12</td><td>0.24</td><td>-3.5</td><td>-2.73</td></tr></tbody></table>

### 3.2. Data and parameters

This subsection lists all the parameters of the multipathogen model with their corresponding values are the references. Table 3 lists the farm module parameters $\theta^{\text{farm}}$. Table 4 shows the frequency distribution of number of cows in different farms, used to simulate the number of cows in different farms. Table 6 lists the physico-chemical parameters required to compute the maximum growth rate $\mu^{\text{max}}$ for different pathogens, and Table 5 lists the other cheese module parameters, which are collectively denoted as $\theta^{\text{cheese}}$. Tables 7 and 8 respectively lists the consumer module $\theta^{\text{con}}$ and the postharvest module parameters $\theta^{\text{post}}$.

## 4. Model evaluation

### 4.1. Batch level outputs

As mentioned in Section 2.6.1, the batch level simulator produces three outputs, namely, the batch risk $R_x^{\text{batch}}(\xi_x)$, the probability of detecting contamination $P_x^{\text{sample}}(\xi_x)$ while testing and the milk loss per due to preharvest testing $M^{\text{batch}}$. These quantities are computed conditionally on the stochastic parameters $\xi_x$ of the model. To compute the relative batch risk for each pathogen, the risk associated with the consumption of a particular batch of cheese, denoted by $R_x^{\text{batch}}(\xi_x)$, is divided by a baseline risk value $R_x^{\text{baseline}}$. $R_x^{\text{baseline}}$ is the average risk of

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Relative batch risk against initial concentration of pathogen in milk $Y_x^{\text{milk}}$, with fixed values of other internal variables in $\xi_x$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Probability of detecting contamination at the time of testing, against initial concentration of pathogen in milk $Y_x^{\text{milk}}$, with fixed values of other internal variables in $\xi_x$.</figcaption></figure>
<caption>**Table 10**</caption>

Baseline prevalence summary statistics (in percentage) computed at time of consumption, for cheese servings (25 g) with no intervention steps.

<table><caption>**Table 11**</caption><thead><tr><th>Symbol</th><th>Mean</th><th>Median</th><th>SD</th><th>q<sub>0.05</sub></th><th>q<sub>0.95</sub></th></tr></thead><tbody><tr><td>MPS-STEC</td><td>1.97</td><td>0.16</td><td>8.21</td><td>0</td><td>7.43</td></tr><tr><td>HV <i>Salmonella</i></td><td>0.37</td><td>0.02</td><td>3.22</td><td>0</td><td>0.91</td></tr><tr><td>LV <i>Salmonella</i></td><td>0.75</td><td>0.06</td><td>4.67</td><td>0</td><td>2.27</td></tr><tr><td>Lm</td><td>39.47</td><td>37.03</td><td>21.97</td><td>8.37</td><td>80.03</td></tr></tbody></table>

Expected probability $E[P_x^{sample}(\Xi_x)]$ of detecting at least one colony in any of $n^{sample} = 5$ sample units of mass $m^{sample} = 25$ g.

<table><thead><tr><th>MPS-STEC</th><th>Salmonella</th><th>Lm</th></tr></thead><tbody><tr><td>0.14</td><td>0.1</td><td>0.5</td></tr></tbody></table>

illness from pathogen $x$ in a baseline scenario which corresponds to a specific situation with no preharvest or postharvest intervention steps. The relative batch risk against the initial concentration of pathogen in milk is plotted in Fig. 12. 1000 independent batches are simulated with fixed values of the internal stochastic variables, with $d^{storage} = 12$ hours, $T^{storage} = 5$ degree Celsius, $t^{consum} = 30$ days, $\rho_{Lm} = 1$, $\delta_{core} = 1.4$ and $\delta_{rind} = 3.1$. This shows the monotonically increasing relationship with respect to $Y_x^{milk}$. For *Salmonella* the computation of risk still includes some extent of randomness due to the use of Monte Carlo method as explained in Section 3.1.2.

In the same framework, the second output of interest the probability of detecting contamination $P_x^{sample}(\xi_x)$, i.e. detecting at least one colony at the time of testing $n^{sample}$ sample units of mass $m^{sample}$ g, is plotted in Fig. 13, as a function of initial pathogen concentration in milk. Primarily $Y_x^{milk}$ is the major influencing factor for determining the probability of detecting contamination as well. While MPS-STEC and *Salmonella* shows similar behaviour, due to high initial concentration, the probability of detection of contamination is higher for $Lm$.

## 4.2. Concentration and prevalence

The model was used to compute the prevalence of the three pathogens during different stages of cheese production. More precisely, the concentration of the pathogens in the ATM and the corresponding summary statistics over different batches were simulated, as listed in Table 9, in presence of no preharvest intervention. Based on the assumptions of the farm module stated in Section 2.3.1, a baseline scenario with no milk testing step was simulated 10,000 times to obtain these prevalence values.

**Remark 3.** It is to be noted that these figures does not represent the

actual scenario of farm milk contamination in France, since these values are based on a hypothetical baseline scenario with no intervention steps, which is not the case in reality. The baseline scenario is used in this study as a reference to compare the impact of the intervention strategies.

The prevalence of the pathogens in the cheese (in a standard serving of 25 g) was also simulated at the time of consumption. The prevalence is defined as the probability of observing at least one colony of the particular pathogen in the food-item. The prevalence is computed using the expected number of colonies given by (9), adjusted with respect to the size of the colonies at the time of consumption. In a baseline scenario 10,000 batches were simulated with different consumption time, to compute the average prevalence in a cheese serving, as listed in Table 10.

Evidently the high prevalence of $Lm$ can be attributed to two major reason, firstly the high concentration in ATM and as shown in Table 9. Secondly, as described in Section 2.3.2, $Lm$ has a secondary growth phase during the ripening step of cheese production, which continues until the time of consumption. The prevalence for all the three pathogens directly affects the postharvest sampling plan.

Consider a postharvest sampling plan with $n^{sample} = 5$ sample units of mass $m^{sample} = 25$ g, for each of the pathogens, such that the sample units are taken at the end of the ripening phase, i.e. at the 14th day of production. Table 11 lists the expected probability $E[P_x^{sample}(\Xi_x)]$ of detecting a colony of the respective pathogen, in any one of the sample units. This is estimated using the average probability of rejection over 10,000 simulated batches with no preharvest intervention step.

Clearly the high prevalence of $Lm$ has a significant effect on the final output of the postharvest sampling plan, i.e. the probability of rejecting the batch $P^{batch}(\xi)$ for detecting at least one of the pathogen colonies, in any one of the sample units. A particular batch will have a high probability of getting rejecting due to the high prevalence of $Lm$. Implementation of a more realistic postharvest sampling plan addressing thisissue is discussed in Section 5.1.

## 4.3. Impact of intervention

As discussed in Sections 2.5.1 and 2.5.2, there exists two types of intervention steps in the cheese production process, and we are interested in the finding the optimal values of the intervention parameters. The impact of the two types of intervention was studied qualitatively, using a series of different intervention scenarios defined using different combinations of intervention parameter values. More specifically, the effect of different intervention scenarios on the relative batch risk, for the preharvest step on and on the probability of batch rejection, for the postharvest step, were monitored.

Fig. 14 shows the impact of different scenarios of the preharvest intervention steps on the relative batch risk, only for two pathogens,

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Effect of preharvest intervention: The empirical cumulative distribution functions of relative batch risk, corresponding to different intervention scenarios (varying sorting limit $l$ and probability of batch testing $p$) and the baseline scenario. The leftmost curve corresponds to the most effective intervention strategy.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Effect of postharvest intervention: The empirical cumulative distribution functions of probability of detecting a colony in any of the sample units, corresponding to different intervention scenarios with varying *n*<sup>sample</sup>. The leftmost curve corresponds to the intervention strategy with highest probability of rejecting a cheese batch.</figcaption></figure>

<table><caption>**Table 12**

Values of intervention parameters to construct different intervention scenarios.</caption><thead><tr><th>Parameter</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td>f<sup>sorting</sup></td><td>1, 2, 10</td><td>Days</td></tr><tr><td>l<sup>sorting</sup></td><td>10, 20, 30, 50, 100, 200</td><td>CFU</td></tr><tr><td>p<sup>test</sup></td><td>0.1, 0.3, 0.5</td><td>Proportion</td></tr><tr><td>n<sup>sample</sup></td><td>1, 5, 10, 15</td><td>Sample units</td></tr></tbody></table>

MPS-STEC and *Salmonella*. Four simple scenarios were considered, implementing four different preharvest intervention strategies, with different values of the milk sorting limit *l*<sup>sorting</sup> = {10,50}, and different probability of testing a particular batch of milk *p* = 1 / *f*<sup>sorting</sup> = {0.5,1}. Corresponding to each of the intervention scenario, 1000 independent batches were simulated to compute the empirical cumulative distribution functions (ECDF) of the relative batch risk. Evidently the most strict scenario with a sorting limit at 10 CFU/ml and with a probability of testing each batch has the leftmost ECDF curve for the relative batch risk, as shown by the orange curve in Fig. 14.

**Remark 4.** Here we study the effect of preharvest intervention only on MPS-STEC and *Salmonella* since the proposed QMRA model uses directly the concentration of *Lm* in ATM, which remains unaffected by the preharvest intervention step (see Section 5.1).

Fig. 15 shows the impact of postharvest intervention step on the probability of detecting a colony in any of the sample units, if it is tested separately for each of the pathogens. The ECDFs are used to study qualitatively, the effect of number of cheese samples a.k.a sample units

taken for the postharvest testing of cheese. The number of sample units considered, were *n*<sup>sample</sup> = {5, 10, 15}, each of mass *m*<sup>sample</sup> = 25g, drawn at the end of ripening step. For each of the postharvest scenarios, 1000 independent batches were simulated, and the probability of detecting a colony in any of the sample units *P*<sub>*x*</sub><sup>sample</sup>(*ξ*<sub>*x*</sub>) was computed separately for three different pathogens. Certainly, a higher number of sample units tend to produce chance of detecting contamination and, in turn, a higher probability of rejecting the batch. For MPS-STEC and *Salmonella*, the batch rejection probabilities were obtained in presence

<table><caption>**Table A1**

Allocation of computational time (using a single core) over different modules of the QMRA simulator, corresponding to the simulation of single batch.</caption><thead><tr><th>Module</th><th>Farm</th><th>Cheese</th><th>Consumer</th></tr></thead><tbody><tr><td>Time spent</td><td>0.7%</td><td>96.4%</td><td>2.9%</td></tr></tbody></table>

<table><caption>**Table A2**

Allocation of computational time among different steps of cheese production and their duration. The steps namely, storage, premolding, draining and salting represents the total time required for all the three pathogens and shows the combined duration for the three pathogens as well. The second growth step only concerns *Lm* and shows the combined duration for the core and rind region.</caption><thead><tr><th>Step</th><th>2nd growth (L. mono)</th><th>Storage</th><th>Premolding</th><th>Draining</th><th>Salting</th></tr></thead><tbody><tr><td>Time spent</td><td>34.9%</td><td>20.9%</td><td>4.9%</td><td>31.4%</td><td>7.9%</td></tr><tr><td>Duration (in h)</td><td>1136</td><td>36</td><td>9</td><td>51</td><td>13.5</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Figures showing the objective space consisting of the cost (in Euros) of intervention and the DALY per one portion (in μ := 1e-6 scale), respectively for a two pathogen (MPS-STEC and *Salmonella*) and three pathogen (MPS-STEC, *Salmonella* and *Lm*) framework. The scatter plots were obtained with 216 different scenarios made up of different combinations of four intervention parameters {*f*<sup>sorting</sup>, *l*<sup>sorting</sup>, *p*<sup>test</sup>, *n*<sup>sample</sup>}, among which the parameters, *f*<sup>sorting</sup> and *p*<sup>test</sup> are indicated using different symbols and colors, respectively.</figcaption></figure>
of no preharvest intervention, and for *Lm* it accounts for the particular preharvest intervention strategy as explained in Sanaa *et al.* (2004). Among the three pathogens, the ECDFs of *Lm* show higher chances of detecting contamination at testing, compared to the other two pathogens.

## 4.4. Objectives of interest

The previous subsection studies qualitatively the effect of intervention parameters on the batch level outputs, namely, the relative batch risk and the probability of detecting contamination. However, one of the goals of the ArtiSaneFood project was to find the optimal values of the intervention parameters that minimizes the two main objectives of the QMRA model, namely, the DALY per one portion $\bar{DALY}_{portion}$ and the total cost of intervention C. These two objectives are chosen in order to take into account the impact of intervention strategies both on the cheese consumers and the cheese producers. As demonstrated Section 4.3, strict intervention strategies can be helpful in reducing the risk of illness for the consumers, but at the same time it increases the probability of rejecting the milk and cheese batches. This trade off gives rise to a multiobjective optimization problem of a stochastic and computationally expensive simulator (see Appendix A), which is addressed using Bayesian approaches of design and analysis of computer experiments (see, e.g., Barracosa et al., 2021; Basak et al., 2023). The problem of optimization is out of the scope of this article, however this subsection provides a visualization of the existing trade-off of the two objectives of interest.

To formalize the optimization problem, four intervention parameters corresponding to the preharvest and postharvest intervention steps, are considered. They are respectively, the frequency of testing the farm milk $f_{sorting}$, the threshold for milk testing $f_{sorting}$, the probability of testing a cheese batch $p^{test}$ and the number of cheese sample units tested $n^{sample}$. The two objectives that are considered for the multiobjective optimization problem, are the the DALY per one portion of cheese $\bar{DALY}_{portion}$, caused by any of the pathogen and the total cost of intervention C, due to the loss of milk, cheese and analytical costs of testing. The computation of cost follows the computation of individual costs for different intervention steps, given by (39).

To visualize the trade-off between the two objectives $\bar{DALY}_{portion}$ and C, a series of different intervention scenarios are considered, using all possible combinations of values of the intervention parameters given in Table 12. A total of 216 different scenarios are constructed with these values and the QMRA model was evaluated for each of them with a Monte Carlo batch size 5000. Fig. 16, demonstrates the relation between the two objectives, using scatter plots on the objective space. Colors are used to indicate the values of the cheese testing parameter $p^{test}$ and different symbols are used to indicate the values of the milk testing parameter $f_{sorting}$. The left figure illustrates the scenario within a dual-pathogen framework, focusing solely on MPS-STEC and *Salmonella*, with both objectives considering only the impact of these two pathogens. The figure on the right incorporates the effects of all the three pathogens. The color transition from blue to red, signifying a stricter postharvest intervention, clearly shows an increase in the cost, however it is not so effective in reducing the the other objective DALY. This behaviour can attributed to the adapted strategy for rejecting a particular batch of cheese, in the postharvest sampling scheme, as described in Section 2.5.2. Unlike MPS-STEC, for the other two pathogens the rejection rule is not based on their corresponding high virulent or highly pathogenic strains. As a result due high prevalence, as shown in Table 10, pathogens like *Lm* are easily detectable when tested, but rather a smaller proportion of the entire population is actually pathogenic, which impacts the DALY. This explains why the effect of postharvest intervention on the DALY, is even less or almost null in the three-pathogen framework. On the other hand, the preharvest intervention parameter $f_{sorting}$ as indicated by the symbols, shows an increasing

impact on the cost and decreasing impact on the DALY metric. This explains the effectiveness of the preharvest intervention scheme on the two objectives of interest, though in the present model the concentration of *L. mono* remains unaffected by this intervention scheme.

# 5. Discussions and perspectives

## 5.1. Model calibration

Model calibration or anchoring is a technique of improving the models efficiency, to be more compatible with observed data. The proposed model is intended to replicate the real life scenario of fabrication of a batch of cheese, however in reality not all practical gestures are possibly reproduced in the model. In this subsection, we list a few perspectives to be adapted in the current model, which are left as future work.

**Preharvest intervention strategy** The current model implements a simple preharvest intervention strategy based on the concentration of *E. coli* in the farms BTM. In other words the acceptance or rejection of farms is based on the *E. coli* concentrations which is assumed to be positively correlated with the concentrations of other pathogenic serotypes, as explained in Section 2.5.1. Although the farms BTM is tested for all the three types of pathogenic contamination, only the *E. coli* concentration can be measured quantitatively. Whereas for the other pathogens, due to low concentrations, only their presence can be detected and the strains can be identified. However in reality preharvest intervention strategies based on contamination of farm milk by *Salmonella* and *Lm* are practiced. To enhance the preharvest intervention strategy further, we can utilize this additional information, which is considered as a future perspective. Moreover, in the proposed QMRA model the concentration of *Lm* in the ATM, that is, the milk to be used for cheese production, is computed directly using (1). This was motivated by the work of Sanaa *et al.* (2004), who performed a Monte Carlo simulation study, based on the data collected in the years 2000–2001, from respectively 347 and 79 farms, respectively, for Camembert and Brie in France. To compute the concentration of *Lm* in the milk used for cheese making, the authors followed the particular milk sorting strategy, for accepting the milk coming from different farms. The preharvest intervention step as explained in Section 2.5.1, is thus not applicable on the concentration of *Lm* and thereafter on the risk values, as simulated using the current version of our proposed model. The implementation of a more realistic and impactful preharvest step is left as a future perspective (see Section 5.1).

**Computation of the DALY metric** The computation of the DALY metric $\bar{DALY}_{portion,x}$ as defined previously, is based on a series of assumptions which allows us to use the available DALY estimates from the literature (Cassini et al., 2018). This approach was adopted to simplify computations and it is to be noted that, the DALY metric thereby estimated is an approximate which requires careful interpretation and usage. A potential future perspective in this context includes using appropriate epidemiological studies to estimate the DALYs related to the consumption of raw milk cheese.

**Postharvest intervention strategy** As mentioned earlier the construction of this model was motivated by an optimization problem to help cheese producers. The optimization was applied to four parameters of the model but other parameters could be of interest to cheese producers. In the context of own-checks, the choice of the sampling date for batches of cheese could prove decisive. In the current model, the date is set at 14 days which represents a certain compromise among different pathogens. The probability of detecting batches at the end of maturation increases for *Lm* but decreases for *Salmonella* and STEC. The other option would be to use two sampling dates optimized for each category of pathogens.

**Pathogen classification based on virulence** Pathogen classification based on virulence refers to categorizing microorganisms, into different groups or classes according to their ability to cause disease. Food-borne
outbreaks tend to result from more virulent serotypes of the pathogens, which are commonly implicated in human illness. In the proposed model, high and low-virulence *Salmonella* (see, e.g., Strickland et al., 2023) serotype proportions were incorporated into each process module to account for their different impact and exposure. For STEC the highly pathogenic and non pathogenic classification, was based on serotypes markers as proposed in Perrin et al. (2014). However, in Auvray et al. (2023) a more robust classification of pathogenic *E. coli* strains was proposed based on their potential virulence. The proposed model did not incorporate this new classification; however, it presents a promising avenue for enhancing the model in the future. For *Lm* there exists classification proportion of high, low and medium virulent serotypes for cold smoked salmon (see, e.g., Fritsch et al., 2018), but as far as our knowledge is concerned, these proportions remain unknown within the context of raw milk soft cheese in the existing literature. Developing a postharvest sampling plan centered around high-virulent serotypes could offer enhanced efficiency from an industrial perspective and provide a rational basis for decision-making.

Second growth step for L. mono As shown in Fig. 7, Lm has no decline phase after the salting step, and the colonies experience a second growth phase depending on the environmental conditions. The second growth phase is modeled using (8), which involves the maximum growth rate parameter $\mu_{Lm}^{\max}(t)$, which depends on the physico-chemical parameters namely the pH, temperature and water activity. The dynamics of these environmental parameters are taken from Perrin et al. (2014), as shown in Fig. 5a, and the available information extends only up to the salting step. According to challenge test data, it has been observed that the second growth step is activated when the pH increases to 6, from 4.52 starting from the end of salting step. Depending on the location of the colonies in the cheese matrix, the pH becomes favorable approximately after the 7th and 20th day of production, respectively for the rind and core region colonies. For simplicity, during the second growth phase, the proposed model uses a fixed value of the maximal growth rate $\mu_{Lm}^{\max}(t)$, which is computed using pH 6, temperature 12.45 degree Celsius and water activity parameter 0.99. However a more realistic implementation can be achieved, by using a more dynamic approach to model the environmental physico-chemical parameters.

## 5.2. Model validation and applicability

This work presents a QMRA model that offers a scientific approach to simulate the real-life scenarios encountered during the production of raw milk soft cheese. The model builds upon the previous work available in the literature, on the QMRA modeling of raw milk soft cheese (see, e. g. Basak et al., 2024; Perrin et al., 2014; Sanaa et al., 2004; Tenenhaus-Aziza et al., 2014) as well as expert opinions from ANSES, CNIEL, ACTALIA, and L2S. The primary goal of this type of model is to study the impacts of different process intervention parameters, in order to implement intervention strategies and make recommendations to cheese producers. However before deploying into industrial applications, careful precautions should be taken in terms of model quality assurance, that includes model verification, validation and calibration. Model verification includes checking the software code used to implement the model and providing a proper documentation. Future directions for this research work would involve the publication of the multipathogen model in the FSX format, thereby promoting open access and facilitating easy and quick peer comparisons. The next crucial step is model validation which ensures the accuracy and reliability of the model's predictions. This involves assessing whether the QMRA model accurately reflects the real-world conditions and produces results that are consistent with observed data. Given the complexity of the QMRA model, the validation step involved a detailed study and precise calibration of various components within the separate modules that replicate different phases of cheese production. As previously mentioned, the proposed QMRA model is based on existing models from the literature

and has been updated with data and expertise collected under the ArtiSaneFood project. The outputs from different components of the QMRA model were compared and calibrated with previous QMRA works and published reports on contamination and outbreaks. For STEC the prevalence rates and baseline risk are compared to Perrin et al. (2014), for *Salmonella* the contamination rates were compared to reports published by the Fédération National des éleveurs de Chèvres (FNEC) in France, and for *Lm* the references from Sanaa et al. (2004), Food and Administration et al. (2012) were compared. The comparisons assured similar orders of magnitude for the simulated outputs of the QMRA model, validating it and making it ready for simulation studies. However, it should be noted that, given the complexity of the model, the outputs can be sensitive to certain input parameters. Therefore, direct comparison of the results with the literature might not be ideal when data-driven parameter values are updated and when the model is enhanced with new functionalities. Despite of model quality assurance, the reliability and applicability of the model still remain subject to ongoing evaluation and refinement. According to World Health Organization et al. (2021) "Models are always incomplete representations of the system they are intended to model, but they can still be useful." Hence, it is essential to note that the outputs obtained using the simulator, such as the batch risk, loss of milk and proportion of rejecting cheese batches are just the estimates of a hypothetical scenario simulated with a state-of-the-art QMRA model. Depending on situations and model inputs these output can be significantly different from the actual prevalence of observed in reality. Nevertheless, the proposed multipathogen model continues to serve as a valuable tool for evaluating the efficacy of intervention strategies and aiding cheese producers in their decision-making processes.

## Funding

This work is a part of the European project ArtiSaneFood funded by the PRIMA program (grant number: ANR-18-PRIM-0015) of the European Union. This project aims on microbial safety of artisanal fermented food produced in the Mediterranean region.

## CRediT authorship contribution statement

**Subhasish Basak:** Writing – review & editing, Methodology. **Laurient Guillier:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Julien Bect:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Janushan Christy:** Writing – review & editing, Methodology, Data curation. **Fanny Tenenhaus-Aziza:** Writing – review & editing, Supervision, Project administration. **Emmanuel Vazquez:** Writing – review & editing, Supervision, Methodology, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

The code will be made available on GitHub.

## Acknowledgements

The authors would like to express their gratitude to PRIMA and ANR for providing financial support for this research work. Special thanks are extended to the partners of the French consortium of the ArtiSaneFood project, for their valuable contributions and insights during the implementation of intervention steps and collection of field data.
# Appendix A. Programming tools and computational time

The multipathogen model is implemented in R and for STEC an implementation of the single pathogen version of the QMRA model is made available in Food Safety Knowledge Markup Language (FSK-ML) format to facilitate its reuse (Basak et al., 2024). This open format is based on predefined terms, metadata and controlled vocabulary to harmonize annotations of risk assessment models (see, e.g., Ungaretti Haberbeck et al., 2018).

The current implementation of the multipathogen model takes around ~ 4.5 s, to simulate one batch. The Table A.13 lists the proportion of computational time (on a single core) spent while running the simulator. Evidently the cheese module takes more than 96% of the computational time, which is dedicated to the time required for solving the ordinary differential equations that models the growth of the pathogens, as shown in (8).

The current implementation in R, uses the ode45 function provided by the pracma (Borchers, 2022) package, that implements the Dormand-Prince (4,5) method. A detailed profiling on the computational time on the cheese module is shown in Table A.14, which precises the proportion of time taken by each of the cheese production steps and their duration in the production process. The number of function evaluations performed by the ODE solver, depends on the duration of the cheese production step, the behavior of the ODE's solution within that duration, and the desired level of accuracy specified by the user. Table A.14 shows the decomposition of the total time consumed by the ODE solver, while simulating a single batch. Except the second growth step of *Lm*, the other cheese production steps runs the ODE solver separately for three pathogens. The time taken is proportional to the total duration of the corresponding production step, however the second growth step spends relatively less time than others, despite of having a long duration. This can be explained by the behaviour of the function, or equivalently the maximum growth rate $\mu_{Lm}^{\text{max}}$ which is considered to be a constant (see, e.g., Section 5.1) for that particular step.

## References

* ANSES, 2020. Méthodologie de hiérarchisation des dangers biologiques et chimiques dans les aliments.

* Augustin, J.C., Zuliani, V., Cornu, M., Guillier, L., 2005. Growth rate and growth probability of Listeria monocytogenes in dairy, meat and seafood products in suboptimal conditions. J. Appl. Microbiol. 99, 1019–1042. https://doi.org/10.1111/j.1365-2672.2005.02710.x.

* Authority, E.F.S., Prevention, E., 2022. The European union one health 2021 zoonoses report. EFSA J. 20 (12), e07666.

* Auvray, F., Cointe, A., Desvaux, M., Jones, G., Kurkdjian, P.M., Oswald, E., Kooh, P., Guillier, L., Audiat-Perrin, F., 2023. Avis de l'Anses relatif à la définition des souches pathogènes d'Escherichia coli productrices de shigatoxines. Technical Report Saisine n°2020-SA-0095. Anses. Citation suggérée : Anses (2023). Avis relatif à la définition des souches pathogènes d'Escherichia coli productrices de shigatoxines (saisine n°2020-SA-0095). Maisons-Alfort : Anses, 63 p.

* Barracosa, B., Bect, J., Dutrieux Baraffe, H., Morin, J., Fournel, J., Vazquez, E., 2021. Extension of the Pareto active learning method to multi-objective optimization for stochastic simulators. SIAM Conference on Computational Science and Engineering (CSE21).

* Basak, S., Bect, J., Vazquez, E., 2022. Integration of bounded monotone functions: revisiting the nonsequential case, with a focus on unbiased Monte Carlo (randomized) methods. 53èmes Journées de Statistique de la SFdS.Lyon, France

* Basak, S., Bect, J., Vazquez, E., 2023. Bayesian multi-objective optimization for stochastic simulators. MASCOT-NUM 2023.

* Basak, S., Christy, J., Guillier, L., Audiat-Perrin, F., Sanaa, M., Tenenhaus-Aziza, F., Bect, J., Vazquez, E., 2024. Quantitative risk assessment of haemolytic and uremic syndrome (HUS) from consumption of raw milk soft cheese. Food Ecol. Syst. Modell. J. 5, e109502.

* Bonifait, L., Thépault, A., Baugé, L., Rouxel, S., Le Gall, F., Chemaly, M., 2021. Occurrence of Salmonella in the cattle production in France. Microorganisms 9 (4), 872.

* Borchers, H. W., 2022. pracma: Practical Numerical Math Functions. R package version 2.4.2. https://CRAN.R-project.org/package=pracma.

* Bruyand, M., Mariani-Kurkdjian, P., Le Hello, S., King, L.A., Van Cauteren, D., Lefèvre, S., Gouali, M., Jourdan-da Silva, N., Mailles, A., Donguy, M.P., et al., 2019. Paediatric haemolytic uremic syndrome related to shiga toxin-producing Escherichia coli, an overview of 10 years of surveillance in France, 2007 to 2016. Eurosurveillance 24 (8), 1800068.

* Cadavez, V., Pouillot, R., Guillier, L., Gonzales-Barron, U., Sanaa, M., unpublished. JEMRA QMRA models for Listeria monocytogenes. R package version 0.2.0. https://github.com/vcadavez/JEMRA/.

* Camargo, A.C., Woodward, J.J., Call, D.R., Nero, L.A., 2017. Listeria monocytogenes in food-processing facilities, food contamination, and human listeriosis: the Brazilian Scenario. Foodborne Pathog. Dis. 14 (11), 623–636.

* Campagnollo, F.B., Gonzales-Barron, U., Pilao Cadavez, V.A., Sant'Ana, A.S., Schaffner, D.W., 2018. Quantitative risk assessment of Listeria monocytogenes in traditional minas cheeses: the cases of artisanal semi-hard and fresh soft cheeses. Food Control 92, 370–379.

* Cassini, A., Colzani, E., Piini, A., Mangen, M.J.J., Plass, D., McDonald, S.A., Maringhini, G., van Lier, A., Haagsma, J.A., Havelaar, A.H., et al., 2018. Impact of infectious diseases on population health using incidence-based disability-adjusted life years (DALYs): results from the burden of communicable diseases in Europe Study, European Union and european economic area countries, 2009 to 2013. Eurosurveillance 23 (16), 17–00454.

* Commission, C. A., 1999. Principles and guidelines for the conduct of microbiological risk assessment. CAC/GL-30.

* Costanzo, N., Ceniti, C., Santoro, A., Clausi, M.T., Casalinuovo, F., et al., 2020. Foodborne pathogen assessment in raw milk cheeses. Int. J. Food Sci. 2020.

* Dalzini, E., Bernini, V., Bertasi, B., Daminelli, P., Losio, M.N., Varisco, G., 2016. Survey of prevalence and seasonal variability of Listeria monocytogenes in raw cow milk from Northern Italy. Food Control 60, 466–470.

* Dubois-Brissonnet, F., Auvray, F., Cerf, O., Delbès, C., Fravalo, P., Gautier, M., Jourdan Da-Silva, N., Korsak Koulagenko, N., Mathieu, F., Oswald, E., Spinnlier, H.E., Talon, R., Bergis, H., Bavourthe, C., Duret, S., Ramarao, N., Guillier, L., Audiat-Perrin, F., Chaix, E., Kooh, P., Sanaa, M., 2022. Avis de l'Anses relatif aux modalités de maîtrise du risque lié à la présence de dangers microbiologiques dans les fromages et autres produits laitiers fabriqués à partir de lait cru. Technical Report Saisine n° 2019-SA-0033. Anses. Citation suggérée : Anses (2022). Avis relatif aux modalités de maîtrise du risque lié à la présence de dangers microbiologiques dans les fromages et autres produits laitiers fabriqués à partir de lait cru (saisine 2019-SA-0033). Maisons-Alfort: Anses, 126 p.

* EFSA, 2012. The European Union summary report on trends and sources of zoonoses, zoonotic agents and food-borne outbreaks in 2010. EFSA J. 10 (3), 2597.

* EFSA, 2017. The European Union summary report on trends and sources of zoonoses, zoonotic agents and food-borne outbreaks in 2016. EFSA Journal 15 (12).

* Farrokh, C., Jordan, K., Auvray, F., Glass, K., Oppegaard, H., Raynaud, S., Thevenot, D., Condon, R., De Reu, K., Govaris, A., et al., 2013. Review of shiga-toxin-producing Escherichia coli (STEC) and their significance in dairy production. Int. J. Food Microbiol. 162 (2), 190–212.

* Food and Drug Administration, et al., 2012. Joint FDA/health Canada quantitative assessment of the risk of listeriosis from soft ripened cheese consumption in the United States and Canada. Center for Food Safety and Applied Nutrition. US Department of Health and Human Services, p. 2014. Retrieved July 15

* Frank, C., Werber, D., Cramer, J.P., Askar, M., Faber, M., an der Heiden, M., Bernard, H., Fruth, A., Prager, R., Spode, A., et al., 2011. Epidemic profile of shiga-toxin-producing Escherichia coli O104: H4 outbreak in Germany. N. Engl. J. Med. 365 (19), 1771–1780.

* Fritsch, L., Guillier, L., Augustin, J.C., 2018. Next generation quantitative microbiological risk assessment: refinement of the cold smoked salmon-related listeriosis risk model by integrating genomic data. Microbial Risk Anal. 10, 20–27.

* Gonzales-Barron, U., Cadavez, V., Valero, A., Skandamis, P., Kintzios, S., de Cesare, A., Manfreda, G., Tenenhaus-Aziza, F., Guillier, L., Boudhrioua, N., Achemchem, F., 2022. Report on the First Predictive Dynamic Models of the Viability of Pathogens along Processing of Mediterranean Artisanal Fermented Foods and Report on the Optimised Process Variables to Enhance their Microbiological Safety. 10.5281/zenodo.8118475.

* Gopal, N., Hill, C., Ross, P.R., Beresford, T.P., Fenelon, M.A., Cotter, P.D., 2015. The prevalence and control of bacillus and related spore-forming bacteria in the dairy industry. Front. Microbiol. 6, 1418.

* Jones, G., Lefèvre, S., Donguy, M.P., Nisavanh, A., Terpant, G., Fougère, E., Vaissière, E., Guinard, A., Mailles, A., de Valk, H., et al., 2019. Outbreak of shiga toxin-producing Escherichia coli (STEC) O26 paediatric haemolytic uremic syndrome (HUS) cases associated with the consumption of soft raw cow's milk cheeses, France, March to May 2019. Eurosurveillance 24 (22), 1900305.

* King, L.A., Espié, E., Haeghebaert, S., Grimont, F., Mariani-Kurkdjian, P., Filliol-Toutain, I., Bingen, E., Weill, F.X., Loirat, C., De Valk, H., et al., 2009. Surveillance du syndrome hémolytique et urémique chez les enfants de 15 ans et moins en France, 1996–2007. Bull. Epidémiol. Hebd. 14, 125–128.

* Koutsoumanis, K.P., Aspridou, Z., 2016. Moving towards a risk-based food safety management. Curr. Opin. Food Sci. 12, 36–41.

* Kretzschmar, M., Mangen, M.J.J., Pinheiro, P., Jahn, B., Fevre, E.M., Longhi, S., Lai, T., Havelaar, A.H., Stein, C., Cassini, A., et al., 2012. New methodology for estimating the burden of infectious diseases in Europe. PLoS Med. 9 (4), e1001205.
Leclercq, A., Charlier, C., Lecuit, M., 2014. Global burden of listeriosis: the tip of the
iceberg. Lancet Infect. Dis. 14 (11), 1027-1028.

Lindqvist, R., Flink, C., Lindblad, M., 2023. Classification and ranking of shigatoxin-
producing Escherichia coli (STEC) genotypes detected in food based on potential
public health impact using clinical data. Microbial Risk Anal. 23, 100246.

Lindqvist, R., Sylvén, S., Vågsholm, I., 2002. Quantitative microbial risk assessment exemplified by staphylococcus aureus in unripened cheese made from raw milk. Int. J. Food Microbiol. 78 (1-2), 155-170.

Maher, M.M., Jordan, K.N., Upton, M.E., Coffey, A., 2001. Growth and survival of E. coli O157: H7 during the manufacture and ripening of a smear-ripened cheese produced from raw milk. J. Appl. Microbiol. 90 (2), 201-207.

Mangen, M.J.J., Plass, D., Havelaar, A.H., Gibbons, C.L., Cassini, A., Mühlberger, N., van Lier, A., Haagsma, J.A., Brooke, R.J., Lai, T., et al., 2013. The pathogen-and incidence-based DALY approach: an appropriated methodology for estimating the burden of infectious diseases. PloS one 8 (11), e79740.

European centre for disease prevention and control, 2019. ECDC BCoDE toolkit,1.7 Solna. https://ecdc.europa.eu/en/toolkit-application-calculate-dalys

McCullough, N.B., Elsele, C.W., 1951. Experimental human salmonellosis: I. Pathogenicity of strains of Salmonella meleagridis and Salmonella anatum obtained from spray-dried whole egg. J. Infect. Dis. 278-289.

McIntyre, L., Wilcott, L., Naus, M., et al., 2015. Listeriosis outbreaks in British Columbia, Canada, caused by soft ripened cheese contaminated from environmental sources. BioMed Res. Int. 2015.

Miszczyza, S.D., Bel, N., Gay-Perret, P., Michel, V., Montel, M.C., Sergentet-Thevenot, D., 2016. Behavior of different shiga toxin-producing Escherichia coli serotypes (O26: H11, O103: H2, O145: H28, O157: H7) during the manufacture, ripening, and storage of a white mold cheese. J. Dairy Sci. 99 (7), 5224-5229.

Murray, C.J.L., Lopez, A.D., 1997. Global mortality, disability, and the contribution of risk factors: Global burden of disease study. Lancet 349 (9063), 1436-1442.

Panel, E.B., Koutsoumanis, K., Allende, A., Alvarez-Ordonez, A., Bover-Cid, S., Chemaly, M., Davies, R., De Cesare, A., Herman, L., Hilbert, F., Lindqvist, R., Nauta, M., Peixe, L., Ru, G., Simmons, M., Skandamis, P., Suffredini, E., Jenkins, C., Monteiro, P.S., Morabito, S., Niskanen, T., Scheutz, F., da Silva Felicio, M.T., Messens, W., Bolton, D., 2020. Pathogenicity assessment of shiga toxin-producing Escherichia coli (STEC) and the public health risk posed by contamination of food with STEC. EFSA J. 18 (1), e05967. https://doi.org/10.2903/j.efsa.2020.5967.

Perrin, F., Tenenhaus-Aziza, F., Michel, V., Miszczyza, S., Bel, N., Sanaa, M., 2014. Quantitative risk assessment of haemolytic and uremic syndrome linked to O157:h7 and non-O157:h7 shiga-toxin producing Escherichia coli strains in raw milk soft cheeses. Risk Anal. 35 (1), 109-128.

Plaza-Rodriguez, C., Ungaretti Haberbeck, L., Desvignes, V., Dalgaard, P., Sanaa, M., Nauta, M., Filter, M., Guillier, L., 2018. Towards transparent and consistent exchange of knowledge for improved microbiological food safety. Curr. Opin. Food Sci. 19, 129-137. https://doi.org/10.1016/j.cofs.2017.12.002.Food Chemistry and Biochemistry * Food Bioprocessing. https://www.sciencedirect.com/science/artic
le/pii/S2214799317301029

Pouillot, R., Hoelzer, K., Chen, Y., Dennis, S.B., 2015. Listeria monocytogenes dose response revisited-incorporating adjustments for variability in strain virulence and host susceptibility. Risk Anal. 35 (1), 90-108.

Ramos, G.L.P.A., Nascimento, J.S., Margalho, L.P., Duarte, M.C.K.H., Esmerino, E.A., Freitas, M.Q., Cruz, A.G., Sant'Ana, A.S., 2021. Quantitative microbiological risk assessment in dairy products: concepts and applications. Trends Food Sci. Technol. 111, 610-616.

Ricci, A., Allende, A., Bolton, D., Chemaly, M., Davies, R., Fernandez Escamez, P.S., Girones, R., Herman, L., Koutsoumanis, K., Norrung, B., Robertson, L., Ru, G.,

Sanaa, M., Simmons, M., Skandamis, P., Snary, E., Speybroeck, N., Ter Kuile, B., Threlfall, J., Wahlstrom, H., Takkinen, J., Wagner, M., Arcella, D., Da Silva Felicio, M.T., Georgiadis, M., Messens, W., Lindqvist, R., (BIOHAZ), E.P.o.B.H., 2018. Listeria monocytogenes contamination of ready-to-eat foods and the risk for human health in the EU. EFSA J. 16 (1) https://doi.org/10.2903/j.efsa.2018.5134.
https://efsa.onlinelibrary.wiley.com/doi/abs/10.2903/j.efsa.2018.5134

Robinson, E., Travanut, M., Fabre, L., Larréché, S., Ramelli, L., Pascal, L., Guinard, A., Vincent, N., Calba, C., Meurice, L., et al., 2020. Outbreak of Salmonella newport associated with internationally distributed raw goats' milk cheese, France, 2018. Epidemiol. Infect. 148, e180.

Ruzante, J.M., Lombard, J.E., Wagner, B., Fossler, C.P., Karns, J.S., Van Kessel, J.A.S., Gardner, I.A., 2010. Factors associated with Salmonella presence in environmental samples and bulk tank milk from US dairies. Zoonoses Public Health 57 (7-8), e217-e225.

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. Int. J. 24 (2), 389-399.

Strickland, A.J., Sampedro, F., Hedberg, C.W., 2023. Quantitative risk assessment of Salmonella in ground beef products and the resulting impact of risk mitigation strategies on public health. J. Food Prot. 86 (6).

Tenenhaus-Aziza, F., Daudin, J.J., Maffre, A., Sanaa, M., 2014. Risk-based approach for microbiological food safety management in the dairy industry: the case of Listeria monocytogenes in soft cheese made from pasteurized milk. Risk Anal. 34 (1), 56-74.

Teunis, P.F.M., 2022. Dose response for Salmonella typhimurium and enteritidis and other nontyphoid enteric Salmonella. Epidemics 41, 100653.

Teunis, P.F.M., Kasuga, F., Fazil, A., Ogden, I.D., Rotariu, O., Strachan, N.J.C., 2010. Dose-response modeling of Salmonella using outbreak data. Int. J. Food Microbiol. 144 (2), 243-249.

Ung, A., Baidjoe, A.Y., Van Cauteren, D., Fawal, N., Fabre, L., Guerrisi, C., Danis, K., Morand, A., Donguy, M.P., Lucas, E., et al., 2019. Disentangling a complex nationwide Salmonella dublin outbreak associated with raw-milk cheese consumption, France, 2015 to 2016. Eurosurveillance 24 (3), 1700703.

Ungaretti Haberbeck, H., Plaza-Rodriguez, C., Desvignes, V., Dalgaard, P., Sanaa, M., Guillier, L., Nauta, M., Filter, M., 2018. Harmonized terms, concepts and metadata for microbiological risk assessment models: the basis for knowledge integration and exchange. Microbial Risk Anal. 10, 3-12. https://doi.org/10.1016/j.mran.2018.06.001 SPECIAL ISSUE on 10th International Conference on Predictive Modelling in Food: Interdisciplinary Approaches and Decision-Making Tools in Microbial Risk Analysis

Van, C.D., Jourdan-da Silva, N., Weill, F.X., King, L., Brisabois, A., Delmas, G., Vaillant, V., De Valk, H., 2009. Outbreak of Salmonella enterica serotype muenster infections associated with goat's cheese, France, March 2008. Eurosurveillance 14 (31), 19290.

Van Kessel, J.A.S., Karns, J.S., Wolfgang, D.R., Hovingh, E., Schukken, Y.H., 2012. Dynamics of Salmonella serotype shifts in an endemically infected dairy herd. Foodborne Pathog. Dis. 9 (4), 319-324.

World Health Organization, 1997. Risk Management and Food Safety: Report of a Joint FAO/WHO Consultation, Rome, Italy, 27 to 31 January 1997. Food & Agriculture Org.

World Health Organization, 2002. Risk Assessments of Salmonella in Eggs and Broiler
Chickens, Vol. 2. Food & Agriculture Org.

World Health Organization, et al., 2021. Microbiological Risk Assessment-Guidance for food, Vol. 36. Food & Agriculture Org.