*L. monocytogenes* in cold-smoked salmon"

In mc2d: Tools for Two-Dimensional Monte-Carlo Simulations

options("width"=100, "digits"=3)
set.seed(666)
knitr::opts_chunk$set(fig.path = "figs/")
knitr::opts_chunk$set(dpi = 96, fig.width = 6, fig.height = 4, dev = "pdf")

The objective of this case study is to assess the risk of invasive listeriosis from consumption
of cold-smoked salmon in France.
The process of interest lays from the end of the production line in the factory, when the
cold-smoked salmon is vacuum-packed, to the consumption.

The data and the model are adapted to illustrate the use of `mc2d`: the results will not and
*should not* be interpreted as an assessment of the actual risk of listeriosis from consumption
of cold-smoked salmon.
Interested readers could refer to \cite{POUILLOT-2007} and \cite{POUILLOT-2009} for a complete
risk assessment on that issue.

The model will be developed in a first section, without considering variability or uncertainty
(deterministic model). Variability will then be introduced in a second section, and a last
section will consider variability and a part of the data uncertainty.

# The Model

In this section, no variability nor uncertainty is considered.
We assess the final level of *L.\ monocytogenes* in the product, the exposure and the risk of
invasive listeriosis for an "average" individual of the "healthy" French
population.^[It makes little sense, but it will help us introducing smoothly the model.]

During the logistic, the retail and the home step, a bacterial growth is modeled considering
*i*) the fluctuating temperature during the various stages and; *ii*) the bacterial competition
with the food flora.
We use the models developed and/or used in \cite{POUILLOT-2007}. The data are adapted from
\cite{POUILLOT-2007} and \cite{DELIGNETTE-2006}:

- The DMS model predicts the bacterial growth during a stage of duration $d$,
  when the temperature is fluctuating, with an intra-stage average temperature $m_{T}$
  and an intra-stage standard deviation of the temperature $s_{T}$. It is written:
  \begin{equation}
  N_{1}=
  \min\left(
    N_{0} +
    \frac{\mu_{ref}}{\ln(10)} \times d \times
    \frac{\left(s_{T}^{2}+\left(m_{T}-T_{min}\right)^{2}\right)}{\left(T_{ref}-T_{min}\right)^{2}},
    \: N_{max}\right) \label{DMS}
  \end{equation}
  if $m_{T}>T_{min}$, with $N_{1}$ the $\log_{10}$ concentration of bacteria ($\log_{10}$ CFU/g)
  at the end of the stage, $N_{0}$ the $\log_{10}$ concentration at the beginning,
  $\mu_{ref}$ the specific growth rate (day$^{-1}$) at a reference temperature $T_{ref}$ (°C),
  $T_{min}$ the minimal growth temperature (°C), and $N_{max}$ the maximum achievable
  concentration ($\log_{10}$ CFU/g). If $m_{T} \leq T_{min}$, $N_{1}=N_{0}$.
- We will use $T_{ref}=25$°C. In this section $N_{max}=7.27\ \log_{10}$ CFU/g.
- The model for *L.\ monocytogenes* uses $\mu_{ref,Lm}=6.2\ \text{day}^{-1}$ and
  $T_{min,Lm}=-2.9$°C.
- The same model is used for the food flora, with $\mu_{ref,ff}=4.1\ \text{day}^{-1}$ and
  $T_{min,ff}=-4.5$°C.
- The growth model for bacterial competition considers the Jameson effect, i.e., the growth of
  *L.\ monocytogenes* and the food flora are stopped as soon as one population reaches $N_{max}$.

In practice, one evaluates $d_{Lm}$ and $d_{ff}$, the time needed for *L.\ monocytogenes*
or the food flora to reach $N_{max}$, and models growth during an effective duration of
$\min(d,d_{Lm},d_{ff})$.
The time needed to reach $N_{max}$ is evaluated by inverting \eqref{DMS}:
$$
d_{\left(N_{1}=N_{max}\right)}=
\left(N_{max}-N_{0}\right)\times
\frac{\ln(10)}{\mu_{ref}} \times
\frac{\left(T_{ref}-T_{min}\right)^{2}}{\left(s_{T}^{2}+\left(m_{T}-T_{min}\right)^{2}\right)}
$$

The other assumptions are:

- A cold-smoked salmon package is homogeneously contaminated with *L.\ monocytogenes*
  at the end of production at a level of $0.1$ CFU/g.
- The food flora level at the end of production is $10^{2.78}$ CFU/g.
- The time-temperature profile is:
    - 1.1 days at an average temperature of 3.2°C (logistic step), with an intra-stage
      SD of 2.1°C.
    - 4.7 days at 5.5°C at retail, with intra-stage SD of 1.0°C.
    - 4.3 days at 8.2°C in the consumer's home, with intra-stage SD of 2.0°C.
- A healthy, non-elderly, non-pregnant individual eats 35 g of this product.
- The individual dose-response model is a one-hit model
  $\Pr(\text{Illness} \mid D)=1-(1-r)^{D}$
  with $r=4.7\times10^{-14}$ for this healthy sub-population. The mean-population
  dose-response (Poisson-distributed dose with mean $D$) is
  $\Pr(\text{Illness} \mid D)=1-\exp(-r\times D)$.

The question is: *What is the risk for this "average" individual?*

```r
Nmax <- 7.3
murefLm <- 6.2; TminLm <- -2.9
murefFF <- 4.1; TminFF <- -4.5
Lm0 <- log10(1); FF0 <- 2.78
d1 <- 1.1; mT1 <- 3.2; sdT1 <- 2.1
d2 <- 4.7; mT2 <- 5.5; sdT2 <- 1.0
d3 <- 4.3; mT3 <- 8.2; sdT3 <- 2.0
conso <- 35
r <- 4.7e-14

modGrowth <- function(duration, mTemp, sdTemp,
                      N0Lm, murefLm, TminLm,
                      N0FF, murefFF, TminFF,
                      Nmax, Tref = 25) {
  N0Lm <- pmin(N0Lm, Nmax)
  N0FF <- pmin(N0FF, Nmax)
  dLm <- (Nmax - N0Lm) * log(10) / murefLm * (Tref - TminLm)^2 /
         (sdTemp^2 + (mTemp - TminLm)^2)
  dLm <- ifelse(mTemp < TminLm & N0Lm != Nmax, Inf, dLm)
  dFF <- (Nmax - N0FF) * log(10) / murefFF * (Tref - TminFF)^2 /
         (sdTemp^2 + (mTemp - TminFF)^2)
  dFF <- ifelse(mTemp < TminFF & N0FF != Nmax, Inf, dFF)
  realDuration <- pmin(duration, dLm, dFF)
  xLm <- N0Lm + (mTemp > TminLm) * murefLm / log(10) *
         (sdTemp^2 + (mTemp - TminLm)^2) / (Tref - TminLm)^2 * realDuration
  xFF <- N0FF + (mTemp > TminFF) * murefFF / log(10) *
         (sdTemp^2 + (mTemp - TminFF)^2) / (Tref - TminFF)^2 * realDuration
  return(list(xLm = xLm, xFF = xFF))
}

x1 <- modGrowth(d1, mT1, sdT1, Lm0, murefLm, TminLm, FF0, murefFF, TminFF, Nmax)
x2 <- modGrowth(d2, mT2, sdT2, x1$xLm, murefLm, TminLm, x1$xFF, murefFF, TminFF, Nmax)
x3 <- modGrowth(d3, mT3, sdT3, x2$xLm, murefLm, TminLm, x2$xFF, murefFF, TminFF, Nmax)
x3
conta <- 10^x3$xLm; conta
expo <- conso * conta; expo
risk <- 1 - (1 - r)^expo; risk

modGrowth is a convenient function for the growth model. Within this function dLm is the time needed for L.\ monocytogenes to reach Nmax, dFF for the food flora, and realDuration is the effective duration of growth during the stage. Note that:

• This function is "vectorized" --- it can take vectors for any parameter and returns a vector. pmin (parallel minimum) is used instead of min (which returns the minimum of all values), and ifelse instead of if.
• It handles all edge cases that can occur in a Monte-Carlo simulation, such as $N_{0} \geq N_{max}$ or $m_{T} \leq T_{min}$, for either bacterial population.

x1, x2 and x3 are the bacterial concentrations at the end of the logistic, retail and home steps, respectively.

Including Variability

We now specify variability distributions for some inputs, following \cite{DELIGNETTE-2006} and \cite{POUILLOT-2007}.

library(fitdistrplus)
library(mc2d)
ndvar(10001)

Specifying Variability Distributions

Initial Contamination

For the initial contamination levels in L.\ monocytogenes, we have 62 enumeration data from a representative sample of positive packages of cold-smoked salmon: 43 samples < $0.2$ CFU/g, 7 at $0.2$ CFU/g, 4 at $0.4$ CFU/g, 2 at $0.6$ CFU/g, and values of 0.3, 1.0, 1.6, 2.4, 5.4 and 7.0 CFU/g \cite{POUILLOT-2007}. We use fitdistrplus to fit a normal distribution to the $\log_{10}$ values, accounting for censoring. The initial concentrations are then modelled through a truncated normal^[so that at least one CFU is included in one 100 g package] on $[-2,\infty)\ \log_{10}$ CFU/g.

For the food flora, we use the distribution from \cite{DELIGNETTE-2006}: $N_{0ff} \sim N(2.78, 1.14)$.

dataC <- data.frame(
  left  = c(rep(NA, 43), rep(.2, 7), .3, rep(.4, 4), 1, 1.6, .6, .6, 2.4, 5.4, 7),
  right = c(rep(0.2, 43), rep(.2, 7), .3, rep(.4, 4), 1, 1.6, .6, .6, 2.4, 5.4, 7)
)
fit <- fitdistcens(log10(dataC), "norm")
fit
Lm0V <- mcstoc(rnorm, mean = fit$est["mean"], sd = fit$est["sd"], rtrunc = TRUE, linf = -2)
FF0V <- mcstoc(rnorm, mean = 2.78, sd = 1.14)

By default, the type of variability modelled is "V" (variability).

Growth Parameters

Distributions derived from \cite{DELIGNETTE-2006}:

• $N_{max} \sim N(7.27,\ 0.86)\ \log_{10}$ CFU/g.
• $\mu_{ref,Lm} \sim N(6.24,\ 0.75)$ day$^{-1}$, truncated on $[0,\infty)$. $T_{min,Lm} \sim N(-2.86,\ 1.93)$°C.
• $\mu_{ref,ff} \sim N(4.12,\ 1.97)$ day$^{-1}$, truncated on $[0,\infty)$. $T_{min,ff} \sim N(-4.52,\ 7.6)$°C.

NmaxV    <- mcstoc(rnorm, mean = 7.27,  sd = 0.86)
murefLmV <- mcstoc(rnorm, mean = 6.24,  sd = 0.75, rtrunc = TRUE, linf = 0)
TminLmV  <- mcstoc(rnorm, mean = -2.86, sd = 1.93)
murefFFV <- mcstoc(rnorm, mean = 4.12,  sd = 1.97, rtrunc = TRUE, linf = 0)
TminFFV  <- mcstoc(rnorm, mean = -4.52, sd = 7.66)

Time-Temperature Profiles

The time-temperature profiles are modelled using Table \ref{tab:4} (adapted from \cite{POUILLOT-2007}).^[$\Gamma$ denotes the Gamma distribution parameterized as $\Gamma(shape,\ scale)$. Exponential($x$) denotes the exponential distribution with mean $x$.] We assume a shelf life of 28 days with $d_{1}+d_{2}+d_{3}\leq28$.^[See the code for how to model this using truncated distributions.]

\begin{table} \caption{\label{tab:4}Time-Temperature Profiles} \centering{}\begin{tabular}{|c|c|c|c|} \hline Stage & Mean Temperature (°C) & Intra-Stage SD of T (°C) & Time (days)\tabularnewline \hline\hline logistic & $N(3.2,\ 2.2)$ trunc.\ on $[-3;25]$ & $\Gamma(1.16,\ 4.61)$ & Exp(1.1)\tabularnewline \hline retail & $N(5.5,\ 2.2)$ trunc.\ on $[-3;25]$ & $\Gamma(0.65,\ 2.09)$ & Exp(4.7)\tabularnewline \hline consumer & $N(8.2,\ 3.8)$ trunc.\ on $[-3;25]$ & $\Gamma(0.35,\ 19.7)$ & Exp(4.3)\tabularnewline \hline \end{tabular} \end{table}

d1V   <- mcstoc(rexp,   rate = 1/1.1)
mT1V  <- mcstoc(rnorm,  mean = 3.2, sd = 2.2, rtrunc = TRUE, linf = -3, lsup = 25)
sdT1V <- sqrt(mcstoc(rgamma, shape = 1.16, scale = 4.61))

d2V   <- mcstoc(rexp,   rate = 1/4.7, rtrunc = TRUE, lsup = 28 - d1V)
mT2V  <- mcstoc(rnorm,  mean = 5.5, sd = 2.2, rtrunc = TRUE, linf = -3, lsup = 25)
sdT2V <- sqrt(mcstoc(rgamma, shape = 0.65, scale = 2.09))

d3V   <- mcstoc(rexp,   rate = 1/4.3, rtrunc = TRUE, lsup = 28 - (d1V + d2V))
mT3V  <- mcstoc(rnorm,  mean = 8.2, sd = 3.8, rtrunc = TRUE, linf = -3, lsup = 25)
sdT3V <- sqrt(mcstoc(rgamma, shape = 0.35, scale = 19.7))

Serving Size

From observed data, a discrete empirical distribution is used \cite{POUILLOT-2007}: values $V = {10, 12, 19, 20, 30, 34, 40, 50, 60, 67.5, 80, 100, 250}$ g, observed $F = {11, 1, 1, 29, 12, 1, 41, 4, 4, 1, 4, 1, 1}$ times.

consoV <- mcstoc(rempiricalD,
  values = c(10, 12, 19, 20, 30, 34, 40, 50, 60, 67.5, 80, 100, 250),
  prob   = c(11,  1,  1, 29, 12,  1, 41,  4,  4,    1,  4,   1,   1))

Applying the Model

r <- mcdata(4.7e-14, type = "0")
x1V <- modGrowth(d1V, mT1V, sdT1V, Lm0V, murefLmV, TminLmV, FF0V, murefFFV, TminFFV, NmaxV)
x2V <- modGrowth(d2V, mT2V, sdT2V,
                 x1V$xLm, murefLmV, TminLmV, x1V$xFF, murefFFV, TminFFV, NmaxV)
x3V <- modGrowth(d3V, mT3V, sdT3V,
                 x2V$xLm, murefLmV, TminLmV, x2V$xFF, murefFFV, TminFFV, NmaxV)

contaV <- 10^x3V$xLm
expoV  <- consoV * contaV
riskV  <- 1 - exp(-r * expoV)
Lm1 <- mc(Lm0V, FF0V, NmaxV, murefLmV, TminLmV, murefFFV, TminFFV,
          d1V, mT1V, sdT1V, d2V, mT2V, sdT2V, d3V, mT3V, sdT3V,
          consoV, r, contaV, expoV, riskV)
Lm1
sLm1 <- mc(contaV = Lm1$contaV, expoV = Lm1$expoV, riskV = Lm1$riskV)
summary(sLm1, probs = c(0, 0.5, 0.75, 0.95, 1))

Lm1 contains all parameters and outputs. sLm1 extracts a short summary.

Final Estimate

If 6.5% of cold-smoked salmon packages are contaminated, 49,090,000 French people are in the "non-susceptible" population, and they consume smoked salmon 6.4 times per year on average, the expected number of listeriosis cases is:

meanRisk  <- mcapply(riskV, "var", mean)
expectedN <- round(0.065 * unmc(meanRisk) * 6.4 * 49090000)
expectedN

Including (a Part of the) Uncertainty

We include both variability and uncertainty, focusing on the uncertainty in initial contamination, growth parameters, and prevalence.

Specifying Uncertainty

Initial Contamination

The uncertainty around initial L.\ monocytogenes contamination is modelled using a bootstrap procedure with bootdistcens from fitdistrplus.

ndunc(101)
bootLm0 <- bootdistcens(fit, niter = ndunc())
MLm0    <- mcdata(bootLm0$est$mean, type = "U")
SLm0    <- mcdata(bootLm0$est$sd,   type = "U")
Lm0VU   <- mcstoc(rnorm, type = "VU", mean = MLm0, sd = SLm0, rtrunc = TRUE, linf = -2)

For the food flora, uncertain hyperparameters $M_{N0ff}$ and $\sigma_{N0ff}$ govern the variable parameter $N_{0ff}$ \cite{DELIGNETTE-2006}: $$ N_{0ff} \sim N(M_{N0ff},\sigma_{N0ff}), \quad M_{N0ff} \sim N(2.78, 0.265), \quad \ln(\sigma_{N0ff}) \sim N(0.114, 0.172) $$

MLm0FF <- mcstoc(rnorm,  type = "U", mean = 2.78,  sd = 0.265)
SLm0FF <- mcstoc(rlnorm, type = "U", meanlog = 0.114, sdlog = 0.172)
FF0VU  <- mcstoc(rnorm,  type = "VU", mean = MLm0FF, sd = SLm0FF)

Growth Parameters

Uncertainty around $\mu_{ref,Lm}$, $T_{min,Lm}$, $\mu_{ref,ff}$, $T_{min,ff}$ and $N_{max}$ is modelled via hyperparameters \cite{DELIGNETTE-2006}.^[Note: there was a typo in \cite{DELIGNETTE-2006} that led to an error in \cite{POUILLOT-2007}: the standard-error for $\ln(\sigma_{\mu ref,Lm})$ is $1.03$, not $-1.03$. We use the correct value here.]

\begin{eqnarray} \mu_{ref,Lm} & \sim & N(M_{\mu ref,Lm},\sigma_{\mu ref,Lm}), \quad M_{\mu ref,Lm} \sim \Gamma(69.7,\ 0.0896), \quad \ln(\sigma_{\mu ref,Lm}) \sim N(1.03,\ 0.191)\ T_{min,Lm} & \sim & N(M_{Tmin,Lm},\sigma_{Tmin,Lm}), \quad M_{Tmin,Lm} \sim N(-2.86,\ 0.459), \quad \ln(\sigma_{Tmin,Lm}) \sim N(0.638,\ 0.208)\ \mu_{ref,ff} & \sim & N(M_{\mu ref,ff},\sigma_{\mu ref,ff}), \quad M_{\mu ref,ff} \sim \Gamma(32.5,\ 0.127), \quad \ln(\sigma_{\mu ref,ff}) \sim N(-0.656,\ 0.221)\ T_{min,ff} & \sim & N(M_{Tmin,ff},\sigma_{Tmin,ff}), \quad M_{Tmin,ff} \sim N(-4.52,\ 1.23), \quad \ln(\sigma_{Tmin,ff}) \sim N(2.00,\ 0.257)\ N_{max} & \sim & N(M_{Nmax},\sigma_{Nmax}), \quad M_{Nmax} \sim N(7.27,\ 0.276), \quad \ln(\sigma_{Nmax}) \sim N(-0.172,\ 0.218) \end{eqnarray}

with $\mu_{ref}>0$ and $T_{min}<25$.

MmurefLm  <- mcstoc(rgamma, type = "U", shape = 69.7,  scale = 0.0896)
SmurefLm  <- mcstoc(rlnorm, type = "U", meanlog = 1.03,   sdlog = 0.191)
murefLmVU <- mcstoc(rnorm,  type = "VU", mean = MmurefLm, sd = SmurefLm, rtrunc = TRUE, linf = 0)

MTminLm   <- mcstoc(rnorm,  type = "U", mean = -2.86, sd = 0.459)
STminLm   <- mcstoc(rlnorm, type = "U", meanlog = 0.638,  sdlog = 0.208)
TminLmVU  <- mcstoc(rnorm,  type = "VU", mean = MTminLm,  sd = STminLm,  rtrunc = TRUE, lsup = 25)

MmurefFF  <- mcstoc(rgamma, type = "U", shape = 32.5,  scale = 0.127)
SmurefFF  <- mcstoc(rlnorm, type = "U", meanlog = -0.656, sdlog = 0.221)
murefFFVU <- mcstoc(rnorm,  type = "VU", mean = MmurefFF, sd = SmurefFF, rtrunc = TRUE, linf = 0)

MTminFF   <- mcstoc(rnorm,  type = "U", mean = -4.52, sd = 1.23)
STminFF   <- mcstoc(rlnorm, type = "U", meanlog = 2.00,   sdlog = 0.257)
TminFFVU  <- mcstoc(rnorm,  type = "VU", mean = MTminFF,  sd = STminFF,  rtrunc = TRUE, lsup = 25)

MNmax     <- mcstoc(rnorm,  type = "U", mean = 7.27,  sd = 0.276)
SNmax     <- mcstoc(rlnorm, type = "U", meanlog = -0.172, sdlog = 0.218)
NmaxVU    <- mcstoc(rnorm,  type = "VU", mean = MNmax,    sd = SNmax)

Prevalence

The prevalence (6.5%) was estimated from 41 positive packages out of 626 tested \cite{POUILLOT-2007}. Assuming 100% sensitivity and specificity, uncertainty around the true prevalence is modelled with a Beta$(1,1)$ prior:

prevU <- mcstoc(rbeta, type = "U", shape1 = 41 + 1, shape2 = 626 - 41 + 1)

Applying the Model

x1VU <- modGrowth(d1V, mT1V, sdT1V,
                  Lm0VU, murefLmVU, TminLmVU, FF0VU, murefFFVU, TminFFVU, NmaxVU)
x2VU <- modGrowth(d2V, mT2V, sdT2V,
                  x1VU$xLm, murefLmVU, TminLmVU, x1VU$xFF, murefFFVU, TminFFVU, NmaxVU)
x3VU <- modGrowth(d3V, mT3V, sdT3V,
                  x2VU$xLm, murefLmVU, TminLmVU, x2VU$xFF, murefFFVU, TminFFVU, NmaxVU)

contaVU <- 10^x3VU$xLm
expoVU  <- consoV * contaVU
riskVU  <- 1 - exp(-r * expoVU)
Lm2 <- mc(Lm0VU, FF0VU, NmaxVU, murefLmVU, TminLmVU, murefFFVU, TminFFVU,
          d1V, mT1V, sdT1V, d2V, mT2V, sdT2V, d3V, mT3V, sdT3V,
          consoV, r, contaVU, expoVU, riskVU)
Lm2
sLm2 <- mc(contaVU = Lm2$contaVU, expoVU = Lm2$expoVU, riskVU = Lm2$riskVU)
summary(sLm2, probs = c(0, 0.5, 0.75, 0.95, 1))

The summary provides the mean, SD, minimum, median ... and a 95% credible interval. The point estimate is the median of the 101 values in the uncertainty dimension; the credible interval runs between the 2.5th and 97.5th percentiles of the uncertainty dimension.

Final Estimate

meanRiskU  <- mcapply(riskVU, "var", mean)
expectedNU <- round(prevU * meanRiskU * 6.4 * 49090000)
summary(expectedNU)

This estimate reflects uncertainty from initial contamination, bacterial growth parameters, and sampling of positive packages. Many other uncertainties (notably around the dose-response model) are not considered here; see \cite{POUILLOT-2007,POUILLOT-2009} for a complete analysis.

The tornado chart in the variability dimension suggests a large impact from the growth rate of L.\ monocytogenes, the storage duration during the consumer step, and the initial contamination level. The tornado in the uncertainty dimension highlights the impact of uncertainty around $N_{max}$ on the mean risk.

torn    <- tornado(Lm2); torn
tornunc <- tornadounc(Lm2, quant = .975); tornunc

plot(torn)
plot(tornunc, stat = "mean risk")

plot(torn)

plot(tornunc, stat = "mean risk")

As a conclusion, this example illustrates how predictive growth models can be implemented within mc2d...

mc2d documentation built on Oct. 23, 2026, 1:06 a.m.