In albgarre/biorisk: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(tidyverse)
library(biorisk)

Description of the case study

The biorisk package also includes functions for simulation and analysis of QMRA models that define the microbial concentration as a discrete variable. As an illustration, we illustrate the calculation of the percentage of milk servings that would cause illness if a tank of milk were contaminated by Listeria, a pathogen present in the dairy farm environment. A modular model is proposed, a first stage of packaging of contaminated milk from a tank in cartons followed by a period of refrigeration, and finally the ingestion of this milk. In all these steps, the bacterial concentration is described using discrete random variables (i.e., we cannot have 2.3 cells in a carton).

Model definition

Initial concentration

As an illustration, we define the initial concentration of microorganisms in the milk tank after heat treatment as a uniform distribution over -1 -- 0 log CFU/g.

logN0 <- Uniform$new("logN0")$ 
  map_input("min", 
            Constant$new("logN0_min", -1)  # with a constant min. value
  )$
  map_input("max", 
            Constant$new("logN0_max", 0)  # and a constant max. value
  )

Partition

In this case study, we model the number of bacterial cellsin the package as a discrete variable (e.g., we use units of CFU/package unit not CFU/ml). The first step we model is the packaging of the product from a tank to a milk carton. Because the volume of the tank is much larger than the volume of the milk carton (~100s of liters vs. 1 liter) we can assume that the former is infinite. Then, the sampling can be modelled as a Poisson process. The package includes the PartitionPoisson element for this.

package <- PartitionPoisson$new("Partition")$
  map_input("logN", logN0)$
  map_input("volume", 
            Constant$new("vol_package", 1000) # Volume of the milk carton (ml)
            )

Note that the output of this element is now in CFU. So, if we run the simulation we get how many Listeria cells are on each milk package.

package$simulate(5, seed = 123)
package$simulations

Growth during storage

The next stage is the microbial growth that takes place during storage of the product. We have a discrete number of cells, so we cannot (or rather, should not) use a continous growth model. Instead, we need to use a discrete growth model (technically, a Yule birth process). The package already includes the DiscreteGrowth element for doing this kind of calculation. Nonetheless, we first need to define the storage time. In this case, it will be a uniform distribution from 0 to 120 h.

stor_time <- Uniform$new("Storage time")$
  map_input("min", Constant$new("t_min", 50))$
  map_input("max", Constant$new("t_max", 150))

Then, we define the discrete growth model. For simplicity, we just give a constant value of 0.1 log CFU/h to $\mu$, although we could do something more complicated (e.g., a secondary model). For the initial concentration, we assign the the element describing the number of cells in the package. Note that, in this case, the input is labelled "N0", not "logN0" (as often in continuous elements). This is done to indicate that the input should be a discrete variable (a number of cells), not a microbial concentration (CFU/ml).

storage <- DiscreteGrowth$new("Storage")$
  map_input("N0", package)$
  map_input("time", stor_time)$
  map_input("mu", Constant$new("mu", .1))

If we run this part of the model, we see this output:

storage$simulate(1000, seed = 242)
storage$histogram() + scale_x_log10()

Because we are working on CFU/package, the histogram shows the number of Listeria cells on each milk carton. One important thing here is the stationary growth phase. Even if the model is discrete, the bacterial behavior is still affected by biological limitations, so we need to account for the carrying capacity of the media (and other reasons for the stationary phase).

However, we need to be a bit careful with the units here. If we assume $\log N_{max} = 8 log CFU/ml$, in a 1000 ml milk carton, we would have $N_{max} = 1e8 \cdot 1e3 = 1e11 CFU/package$. This can be implemented easily in biorisk using the DiscreteGrowthNmax element.

storage <- DiscreteGrowthNmax$new("Storage")$
  map_input("N0", package)$
  map_input("time", stor_time)$
  map_input("mu", Constant$new("mu", .1))$
  map_input("Nmax", Constant$new("Nmax", 
                                 1000*1e8
                                 ))

If we now run this part of the model...

storage$simulate(1000, seed = 421)
storage$histogram() + scale_x_log10()

We observe a very large bar on the right. This represents all the packages with a population on the stationary growth phase. This is often quite relevant for discrete growth, because the Yule process uses a negative binomial, which is right-tailed. As a result, the number of simulations above $N_{max}$ is often larger than calculated using a continuous model.

Consumer phase

Next, we need to simulate the number of cells that are consumed per serving. This is again a sampling process. However, here we cannot make the assumption of infinite volume (the glass will be, at least, a 10th of the package). Hence, we need to use a binomial process instead of a Poisson process. Therefore, we use PartitionBinomial, where we define the serving size as a uniform distribution between 100 and 500 ml.

dose <- PartitionBinomial$new("Dose")$
  map_input("N0", storage)$
  map_input("volume0", 
            Constant$new("volume0", 1000)
            )$
  map_input("volume1",
            Uniform$new("serving_volume")$
              map_input("min", Constant$new("min_vol", 100))$
              map_input("max", Constant$new("max_vol", 500))
            )

Dose-response & risk characterization

Finally, we introduce the dose-response model. We use the exponential dose response model with $r = 10^{-12}$, assuming constant pathogen-host survival probability.

Pill <- DoseResponse_Exponential$new("Pill")$
  map_input("r", Constant$new("r_dr", 1e-12))$
  map_input("dose", dose)

and calculate the number of cases assuming that each MC simulation represents one serving.

cases <- Pill2Cases_1$new("cases")$
  map_input("Pill", Pill)

Checks and model visualization

We can first check the consistency of the input types. Unlike in Case Study I, we do not get any warning in this case because the dose is modelled as a discrete variable.

cases$check_input_types(recursive = TRUE)

As always, it is important to plot the model to check all the dependencies are correct.

plot_model(cases)

Remember that we can always zoom-in on some parts of the model.

plot_model(dose)

Simulation

We can now simulate the model. We just need to call the simulate() method with the number of Monte Carlo iterations.

cases$simulate(1000, seed = 42124)

In this model, each MC simulation represents a single serving. That is, the output variable only takes binary values (0/1). This makes result interpretation a bit harder. For instance, we cannot look at quantiles because they will jump directly from 0 to 1 (although they can be 0.5 at a particular transition quantile):

quantile_table(cases, chosen = "cases", probs = c(.5, .9, .99))

Therefore, probably the best way to observe the output is by its mean:

cases$get_output() %>% mean()

That is, 1.6% of the servings will result in disease.

We can also look at the distribution of the simulations to evaluate the uncertainty of that estimate; e.g., using a histogram:

Pill$histogram()

albgarre/biorisk documentation built on Feb. 23, 2025, 5:19 a.m.