The set of inputs parameters of the farm module denoted by \(\theta^{\rm farm}\), are described in Table 1.

The concentration of STEC in farm milk is usually too low to be assessed quantitatively through microbiological methods (because of their limit of detection). For this reason, we rely on the "relative" approach proposed by Perrin et al. (2014), which uses collected data*4 on Escherichia coli (E. coli) to estimate the STEC concentration. This approach is based on the assumption that E. coli and STEC strains follow the same faecal routes in the cows.

For each farm \(i = 1, 2, ..., N^{\rm farms}\), the milk is collected into a bulk tank from \(N^{\rm cow}_i\)cows, where \(N^{\rm cow}_i\) is simulated from a cow distribution estimated from data (see Endnote *5). The STEC concentration in the bulk tank milk (BTM) corresponding to farm \(i\) is denoted by \(Y_{0, i}\) and is computed using a proportion rule as:

\(Y_{0, i} = Y^{\rm EC}_i \frac{\overline{F^{\rm STEC}_i}}{\overline{F^{\rm EC}_i}},\)

where \(Y^{\rm EC}_i\) is the concentration of E. coli (CFU/ml) in bulk milk tank, \(\overline{F^{\rm STEC}_i}\) is the average STEC concentration (CFU/gram) and \(\overline{F^{\rm EC}_{i}}\) is the average E. coli concentration in faecal matter (CFU/gram) from all the cows.

• Pre-harvest intervention*6

Bulk tank milk coming from each farm is tested for E. coli concentration and the farms with concentration \(Y_i^{\rm EC}\) higher than a certain threshold \(l^{\rm sorting}\) are rejected. Milk sorting is not performed for every batch of milk: instead, it is controlled by a parameter \(f^{\rm sorting}\), which denotes the frequency (in days) of testing the farms milk. Let \(S\) denote the set of farms that are qualified after milk sorting and let \(N^{\rm farms, sorted} = |S|\). The milk loss for a particular batch is given by \(M^{\rm batch} = \sum_i^{N^{\rm farms}}V_{i}\, {1}_{\{i \in S\}}\), where \(V_i = q_{\rm milk} N_i^{\rm cow}\) is the amount of milk produced by the \(i\)-th farm\(N_i^{\rm cow}\)*5. After milk sorting, milk from all the qualified bulk tanks is collected into a single aggregated milk tank. Hereafter, only the farms with accepted level of E. coli are considered.

For the \(i\)-th accepted farm, \(\overline{F^{\rm STEC}_i}\) is the average of the individual STEC concentrations \(F^{\rm STEC}_{i,\,j}\) in infected cows \(j = 1, 2, ..., k_{i}\), where \(k_i\) is the number of infected cows in the \(i\)-th farm simulated according to a binomial distribution:

\(k_{i} \sim {\rm Bin}(N^{\rm cow}_i, p_i) \), where \({\rm logit}(p_i) = u_i\) and \(u_i \sim \mathcal{N}(\mu_{u}, \tau_{u})\).

The number of cows infected with MPS-STEC in the \(i\)-th farm is also simulated according to a binomial distribution:

\(k^{\rm MPS}_{i} \sim {\rm Bin}(k_i, p^{\rm MPS}_{\rm STEC}).\)

Remark: If we are interested to compute the concentration of MPS-STEC in the aggregated milk, we use \(k_i^{\rm MPS}\) instead of \(k_i\)in the following computations:

\(F^{\rm STEC}_{i,\,j}\) is simulated according to a Weibull distribution:

\(\begin{equation}\begin{split}F^{\rm STEC}_{i,\,j} &\sim {\rm Weibull}(a^{\rm weibull}, b^{\rm weibull}),\\\overline{F^{\rm STEC}_i} &= \frac{1}{N_{i}^{\rm cow}}\sum_{j=1}^{k_{i}}F^{\rm STEC}_{i,\,j}.\end{split}\end{equation}\)

The concentration of E. coli (CFU/ml) in BTM \(Y^{\rm EC}_i\), is modelled by a lognormal distribution:

\(\begin{equation}Y^{\rm EC}_i \sim {\rm Lognormal}(\alpha_{i}, \sigma_i).\label{eq:farm_alpha_sigma}\end{equation}\)

\(\overline{F^{\rm EC}_{i}}\) is the average of individual E. coli concentration in faecal matter for each cow, denoted by \(F^{\rm EC}_{i,\,j}\), \(j = 1, 2, ..., N^{\rm cow}_{i}\) cow. It is simulated according to a Lognormal distribution:

\(\begin{equation}\begin{split}&\log_{10}(F^{\rm EC}_{i,\,j}) \sim \mathcal{N}(\mu^{\rm ecoli}, \tau^{\rm ecoli}), \\&\overline{F^{\rm EC}_{i}} = \frac{1}{N^{\rm cow}_i}\sum_{j=1}^{N^{\rm cow}_i}F^{\rm EC}_{i,\,j}.\end{split}\end{equation}\)

The STEC concentration \(Y_0\) (in CFU/ml) in this aggregated milk tank is then given by:

\(\begin{equation}Y_0 = \sum_{i=1}^{N^{\rm farms}} \bigg (Y_{0, i} \cdot \frac{V_{i}{1}_{\{i \in S\}}}{\sum_{i=1}^{N^{\rm farms}}V_{i}{1}_{\{i \in S\}}} \bigg).\label{eq:farm_atm_STEC}\end{equation}\)

The quantities of interest from the farm module are \(Y_0\), the milk loss due to sorting, the \(M^{\rm batch}\) total milk utilised \M^{\rm batch} and the number of farms discarded (due to sorting), i.e. \(N^{\rm farms} - N^{\rm farms, sorted}\). Fig. 2 plots the histogram for the final STEC concentration in the aggregated milk tank without milk sorting testing.

The inputs of the cheese module are the initial concentration \(Y_0\) of STEC and the parameters described in Table 2, which are denoted as \(\theta^{\rm cheese}\).

The evolution of STEC involves six steps, with milk storage and moulding taking place during the liquid growth phase and draining and salting occurring during the solid growth phase. Ripening and cheese storage represent the (solid) decline phase of STEC.

In all the growth steps, the concentration \(y(t)\) for STEC at time \(t\), is modelled using an ordinary differential equation:

\(\begin{equation}\frac{dy}{dt} = \mu^{\rm max}(t) y(t) \bigg(1-\frac{y(t)}{y^{\rm max}}\bigg),\end{equation}\)

where \(\mu^{\rm max}(t)\) stands for the maximum growth rate (in \(\rm h^{-1}\)) and \(y^{\rm max}\) is a parameter that represents the hypothetical maximum population of STEC strain in milk or cheese. For each step, the physico-chemical parameters \(\{d, pH, T, a_w\}\) are computed from Table III in Perrin et al. (2014). The maximum growth rate \(\mu^{\rm max}(t)\) for a particular step, at time \(t\) is according to model 5 in Augustin et al. (2005). \(\mu^{\rm max}(t)\) is dependent on the values of physico-chemical parameters \(\{d, pH, T, a_{\rm w}\}\) at time \(t\) and the optimal growth rate parameter \(\mu^{\rm opt}\). The optimal growth rate is taken to be the average of the growth rates corresponding to different strains as provided in Table II in Perrin et al. (2014).

The milk storage step starts with initial concentration \(Y_0\) and outputs the final concentration \(Y^{\rm storage}\) after growth at the end of storage. The moulding step takes \(Y^{\rm storage}\) as input to model the corresponding growth of the bacteria during moulding and outputs the final concentration \(Y^{\rm molding}\). After moulding, the milk is converted into solid cheese and the STEC bacteria from colonies inside a single cheese made from a volume \(v^{\rm cheese}\) (ml) of milk. The number of colonies in a cheese is modelled as a Poisson variable \(N^{\rm colony}_s \sim \rm Poisson(\lambda^{\rm colony}_{\rm s})\), for each strain \(s\) of MPS-STEC \(s \in \{ {\rm O157}, \overline{\rm O157}\}\), with the following mean parameters:

\(\lambda^{\rm colony}_{\rm O157} = Y^{\rm molding} \cdot v^{\rm cheese}\cdot w^{\rm loss} \cdot p_{\rm O157} \cdot p^{\rm STEC}_{\rm MPS}\),

\(\lambda^{\rm colony}_{\overline{\rm O157}} =Y^{\rm molding} \cdot v^{\rm cheese} \cdot w^{\rm loss} \cdot(1-p_{\rm O157}) \cdot p^{\rm STEC}_{\rm MPS}\).

Remark: The factor \(p^{\rm STEC}_{\rm MPS}\) takes into account only the concentration of pathogenic STEC bacteria (MPS-STEC). Alternatively, this can be taken into account by directly computing the concentration of MPS-STEC in the aggregated milk tank as an output of the farm module using \(k_i^{\rm MPS}\) instead of \(k_i\). The FSKX implementation provides a flag flag_MPS = FALSE that allows the user to directly compute the concentration of MPS-STEC in the farm module using the random proportion of MPS-STEC infected cows in a farm.

Starting from the draining phase, the evolution of the size of colony, stemming for one immobilised STEC cell, is studied. The draining phase commences with an initial colony size of 1 CFU and growth continues until the salting phase. It is assumed that the evolution of each colony inside each cheese (weighing 250 g) is identical during these phases, since they have the same environmental conditions. The output of the draining and salting steps are called \(Y^{\rm draining}\) and \(Y^{\rm salting}\), respectively.

The growth in colony size stops after the salting phase. Then starts the decline phase, which is composed of two steps—namely, ripening and cheese storage. The ripening phase lasts until the 14^th day of cheese production and the cheese storage phase duration depends on the consumption time \(t^{\rm consum}\) of the cheese. The rate of decline is different for different strains—namely, MPS (\({\rm O157}, \overline{\rm O157}\)) and non-MPS STEC. The median number of bacteria per colony, for a particular strain \(s\) after decline, at the end of the ripening phase and at the time of consumption, is computed as:

\(Y^{\rm ripening}_s =Y^{\rm salting} \cdot 10^{- \rho_s \times (14 \times 24 - t)/24},\)

\(Y^{\rm consum}_s =Y^{\rm salting} \cdot 10^{- \rho_s \times (t^{\rm consum} \times 24 - t)/24}.\)

where \(t\) the time (in hours) taken until the salting step.

The expected number of colonies \(\lambda^{\rm colony}_s\) is adjusted taking the inactivation during the ripening and cheese storage phases. For instance if the median number of bacteria for a particular strain \(Y^{\rm consum}_s\) falls below 1, it signifies that the colony might have disappeared with probability \(Y^{\rm consum}_s\) and the corresponding expected number of colonies is obtained by multiplying \(\lambda^{\rm colony}_s\) by \(Y^{\rm consum}_s\). In case of \(Y^{\rm consum}_s > 1\), the expected number of colonies remains unchanged.

\(\lambda^{\rm colony}_s \leftarrow \lambda^{\rm colony}_s \cdot \min(1, Y_s^{\rm consum})\).

The outputs of interest for the cheese module are the average number of colonies \(\lambda^{\rm colony}_{s}\) and the median colony size at the time of consumption \(Y^{\rm consum}_s\), where the subscript s denotes the strain \(\{ {\rm O157},\overline{\rm O157}\}\). Fig. 3 shows the growth of STEC concentration (CFU/ml) starting from \(Y_0 = 10^{-3.8}\) CFU/ml, during the storage and moulding steps, for the baseline scenario (with no interventions and default values of the parameters). After the transition from liquid to solid state, the evolution of a single colony is shown in Fig. 4. The colonies grow during the draining and salting step and then decline during the ripening and cheese storage steps.

The dose \(\Gamma\), corresponding to the concentration of MPS-STEC per \(25\) g of cheese serving, is computed as:

\(\Gamma = \sum_s N^{\rm colony}_{\rm s, sample} Y^{\rm colony}_{s}\)

where \(N^{\rm colony}_{\rm s, sample}\) denotes the number of colonies in a cheese serving of weight \(wt^{ \rm serving}\) g, which is a Poisson random variable with mean \(\lambda^{\rm colony}_s \times \frac{wt^{ \rm serving}}{wt^{\rm cheese}}\). \(Y^{\rm colony}_s\) denotes the size of all colonies in a cheese, which follows a Log-Normal distribution:

\(Y_s^{\rm colony} = Y^{\rm consum}_s \cdot 10^{\epsilon_s},\)

where \(\epsilon_s \sim \mathcal{N}(0, \tau_{\epsilon_s})\) represents inter-batch variablity (i.e. variability between different cheeses in a single batch). Note that here we assume all the colonies inside a single cheese, for a particular batch and a given strain \(s\), have identical size \(Y^{\rm colony}_{s}\).

The average number of colonies \(\lambda_s^{\rm colony}\) depends on several random quantities, such as the initial STEC concentration \(Y_0\),\(\) \(t^{\rm storage}\), \(d^{\rm storage}\)and \(t^{\rm consum}\) (see Table 2). The median colony size \(Y_s^{\rm consum}\) depends on the random variable \(t^{\rm consum}\). In summary, the dose \(\Gamma\) is a combination of several Poisson and Log-Normal random variables corresponding to different strains and their parameters are dependent on the random quantities \(\xi^{\rm dose} = \{Y_0, t^{\rm storage}, d^{\rm storage}, t^{\rm consum}\}\). The probability of HUS (risk) associated with the ingestion of a dose \ (\gamma\), for age \(a\) is:

\(\begin{equation}P({\rm HUS}| a, \gamma) = 1-(1-r_{\rm a})^\gamma,\end{equation}\)

where \(r_{\rm a} = r_0 \times e^{-k a}\) for age \(a\) varying from \(1\) to \(15\). The risk of HUS for a particular batch of milk, conditional on \(\xi^{\rm dose}\), can be derived by integrating w.r.t \(\Gamma \vert \xi^{\rm dose}\):

\(R^{\rm batch}= \sum_{a=1}^{15}g(a)\int_0^{\infty}[1-(1-r_{\rm a})^\gamma]\,p(\gamma| \xi^{\rm dose})\,\rm{d}\gamma\) ,

where \ (p(\gamma|\xi^{\rm dose})\) is the conditional probability density function of \(\Gamma\) given \(\xi^{\rm dose}\) and \ (g(a)\) is the proportion of cheese consumed by the age group \(a\).

The inputs of the consumer module are the average number of colonies \(\lambda_s^{\rm colony}\), the median colony size \(Y_s^{\rm consum}\) and other parameters denoted by \(\theta^{\rm con}\), which are listed in Table 3.

For a given set of input parameters \(\xi^{\rm dose}\), the conditional risk \(R^{\rm batch} =P({\rm HUS}|\xi^{\rm dose})\) can be computed using simple Monte Carlo using i.i.d. samples from \(p(\Gamma|\xi^{\rm dose})\). The number of Monte Carlo samples is given by the numerical parameter \(N^{\rm dose}\). Note that the randomness in \(\Gamma|\xi^{\rm dose}\) is only due to the fact that the dose is a combination of several pairs of a Poisson and Lognormal random variables.

Alternatively, the batch risk can be computed by approximation of the integral \(\mathbb{E}[(1-r_a)^{\Gamma}]\). Since the two classes of MPS-STEC strains are independent, we can write the dose as a sum \(\Gamma = \Gamma_1 + \Gamma_2\) with:

\(\Gamma_s = Y^{\rm consum}_s \cdot N_s^{\rm colony} \cdot 10^{\tau_{\epsilon_s}\nu_s}= d_s \cdot N_s^{\rm colony}\cdot b_s^{\nu_s}\),

where \(\nu_s\) is a standard normal variable. Now the expectation can be written as \(\mathbb{E}[(1-r_a)^{\Gamma}] = \prod_s \mathbb{E}[\mathbb{E}[(1-r_a)^{\Gamma_s}|\nu_s]]\). Using the result in Endnote*2, the expectation can be further reduced to:

\(\mathbb{E}[(1-r_a)^{\Gamma}] = \exp(-\sum_s \lambda^{\rm colony}_s)\prod_s\mathbb{E}[c_{s,1}^{c_{s, 2}^{c_{s, 3}^{\nu}}}] ,\)

where \(c_{s, 1} = \exp(\lambda_s^{\rm colony}), c_{s, 2} = (1-r_a)^{d_s}\)and \(c_{s, 3} = b_s\). Since \(c_{s, 1} > 1, c_{s, 2} < 1\) and \(c_{s, 3} > 1\), the function

\((c_{s,1})^{(c_{s, 2})^{c_{s, 3}^{\nu_s}}}\)

is monotone (non-increasing) in \(\nu_s\). For such functions, deterministic quadrature methods (such as the trapezoidal rule) have convergence rate \(O(n^{-1})\), which is better than simple Monte Carlo (see, for example, Novak (1992)).

The choice of the method for computating the batch risk is determined through the parameter value \(N^{\rm dose}\). For a non-zero integral value, the simulator uses the simple Monte Carlo method and, if \(N^{\rm dose}\) is set to zero, it used the integral approximation method.

The current implementation also allows the user to compute the conditional batch risk \(\xi^{\rm dose}\setminus t^{\rm consum}\)averaged out with respect to the consumption time. This is done by using the flag flag_consum = TRUE that computes the integration of \(R^{\rm batch}\) with respect to \(t^{\rm consum}\). The integral \(I = \mathbb{E}_{t^{\rm consum}}[R^{\rm batch}(t^{\rm consum})]\) is approximated using a piecewise constant function on a regular grid \(\{t_1^{\rm consum}, t_2^{\rm consum}, \cdots, t_n^{\rm consum}\}\) of size \(n\), that spans the support of the triangular distribution of the variable \(t^{\rm consum}\). For each interval, the value of the piecewise constant function is taken equal to the value of \(R^{\rm batch}\) at the left endpoint of that interval and the approximation is computed as:

\(\hat{I} = \sum_{i=1}^{n-1}P[t_i^{\rm consum} \leq t^{\rm consum} \leq t_{i+1}^{\rm consum}] \cdot R^{\rm batch}(t_i^{\rm consum})\).

For STEC, due to the inactivation during the post-ripening storage phase, \(R^{\rm batch}\) turns out to be a decreasing function of \(t^{\rm consum}\), which slightly overestimates the risk with this approach. It is worth noting that the computation of \(R^{\rm batch}\) for a grid of values is not very expensive compared to the computation of the same for a single value of \(t^{\rm consum}\), since both of them require a single run to the cheese module that can compute all the necessary outputs, \(\lambda_s^{\rm colony}\)and \(Y^{\rm consum}_s\) corresponding to the consumption time points.

The output of the consumer module is the estimated batch risk of HUS. For a given set of input parameters, Fig. 6 plots the relative risk \(R^{\rm batch}/R^{\rm baseline}\) against \(log10(y_0)\). The baseline risk is computed with no intervention steps.

The inputs of the post-harvest module are initial STEC concentration in milk \(Y_0\), average number of colonies \(\lambda_s^{\rm colony}\) and the parameters denoted by \(\theta^{\rm post}\), listed in Table 6.

We assume that the colonies are homogeneously distributed inside a cheese, a colony contains at least one MPS-STEC cell and the test is accurate enough to detect a colony with a single MPS-STEC cell. We observe that probability of a sample unit testing positive is \(P^{\rm sample} = P[N^{\rm colony}_{\rm sample}>0]\), where \(N^{\rm colony}_{\rm sample}\) is the total number of colonies from all strains. Since the number of cells corresponding to two different strains grow independently, the total number of colonies is also a Poisson random variable, \(N^{\rm colony}_{\rm sample} \sim \rm Poisson((\lambda^{\rm colony}_{1} + \lambda^{\rm colony}_{2})\cdot \frac{m^{\rm sample}}{wt^{\rm cheese}})\). For given values of \(\xi^{\rm dose}\) using this Poisson distribution, we compute the probability of rejecting a particular batch \(\) \(P^{\rm reject} = 1-(1-P^{\rm sample})^{n^{\rm sample}}\)\(\).

The output of the post-harvest module is the probability of rejection. Fig. 5 plots the probability of rejecting \(P^{\rm reject}\) with fixed \(m^{\rm sample} = 25\) and \(n^{\rm sample} = 5\), for \(10^3\) batches corresponding values of \ (y_0\).

The inputs of the output module are the outputs of the farm, consumer and post-harvest module along with the parameters denoted by \(\theta^{\rm output}\), listed in Table 4.

The output module simulates several batches and produces estimates of the overall risk \(R^{\rm HUS}\), the probability of rejection \(P^{\rm avg}\) and the milk loss \(M^{\rm avg}\). We define the average milk loss \(M^{\rm avg} = \mathbb{E}[M^{\rm batch}]\), average batch rejection rate: \(\)

\(P^{\rm avg} = \mathbb{E}[P^{\rm batch} \cdot p^{\rm test}] = \int_{\xi^{\rm dose}}P^{\rm batch} \cdot p^{\rm test} \cdot p(\xi^{\rm dose})\, {\rm d}\xi^{\rm dose},\)

with \(p^{\rm test}\) being the proportion of batch tested and the average over batch risk: \(R^{\rm avg} = \mathbb{E}[R^{\rm batch} \cdot (1-P^{\rm batch} \cdot p^{\rm test})]\),

where the batch risk is set to zero for rejected batches.

The overall risk of HUS is conditional on the event that the batch actually goes into the market, i.e. not rejected. It is computed by dividing \(R^{\rm avg}\) by the probability that a produced batch is not rejected:

\(R^{\rm HUS} = \frac{\int_{\xi^{\rm dose}}R^{\rm batch} \cdot (1-P^{\rm batch} \cdot p^{\rm test}) \cdot p(\xi^{\rm dose}){\rm d}\xi^{\rm dose}}{1-P^{\rm avg}}\).

The quantities are estimated using simple Monte Carlo with sample size \(N^{\rm batch}\). We use i.i.d. samples \(\{\xi^{\rm dose}_{1}, \xi^{\rm dose}_{2}, \cdots , \xi^{\rm dose}_{N^{\rm batch}}\}\), drawn from the conditional distribution \(p(\xi^{\rm dose})\) to construct the unbiased simple Monte Carlo estimator:

\(\hat{P}^{\rm avg} = \frac{1}{N^{\rm batch}}\sum_{l=1}^{N^{\rm batch}}\hat{P}^{\rm batch}(\xi^{\rm dose}_{l}) \cdot p^{\rm test}\),

\(\hat{R}^{\rm avg} = \frac{1}{N^{\rm batch}}\sum_{l=1}^{N^{\rm batch}}\hat{R}^{\rm batch}(\xi^{\rm dose}_{l}) \cdot (1-\hat{P}^{\rm batch}(\xi^{\rm dose}_{l}) \cdot p^{\rm test})\).

The current implementation of the output module computes the relative risk of HUS with respect to a baseline scenario with no intervention steps. This quantity is obtained by dividing \(R^{\rm HUS}\) by the risk of HUS in the baseline scenario (the baseline risk is also estimated by simple Monte Carlo with sample size \(N^{\rm batch}\)).

The output module returns the relative risk of HUS \(R^{\rm HUS}\), average proportion of rejected batches \(P^{\rm avg}\) and average milk loss \(M^{\rm avg}\).