Dynamic predictive microbiology models • predmicror

Overview

The dynamic modelling helpers in predmicror are designed for prediction and fitting when environmental conditions change over time. This is useful, for example, when temperature varies during storage, transport, processing, or abuse scenarios.

The basic workflow is:

define an environmental profile with dynamic_profile();
simulate microbial behaviour with predict_dynamic_growth() or predict_dynamic_inactivation();
fit selected dynamic parameters with fit_dynamic_growth() or fit_dynamic_inactivation();
inspect fitted values, residuals, metrics, and sensitivity coefficients.

A dynamic temperature profile

A dynamic profile stores time-varying environmental conditions. At present, the most common use case is temperature as a function of time.

profile <- dynamic_profile(
  time = c(0, 5, 10, 15, 20, 25, 30),
  temperature = c(10, 4, 2, 15, 15, 2, 10)
)

profile
#> Dynamic environmental profile
#>   points: 7
#>   time range: [0, 30]
#>   variables: temperature
#>   time temperature
#> 1    0          10
#> 2    5           4
#> 3   10           2
#> 4   15          15
#> 5   20          15
#> 6   25           2
#> 7   30          10

plot(profile$time, profile$temperature,
  type = "b",
  xlab = "Time",
  ylab = "Temperature"
)

Dynamic growth prediction

The Huang dynamic growth model can be combined with a secondary square-root model for the effect of temperature on the maximum specific growth rate.

growth_pred <- predict_dynamic_growth(
  profile = profile,
  model = "huang",
  secondary = "huang_sqrt",
  start = list(
    logN0 = 2,
    logNmax = 8.8,
    a = 0.0886,
    Tmin = 8.91,
    lag = 2
  ),
  times = seq(0, 30, by = 1)
)

head(growth_pred)
#> dynamic_growth prediction
#>   model: huang
#>   secondary: huang_sqrt
#>   rows: 6
#>   time response      lnN     logN temperature
#> 1    0 2.000000 4.605170 2.000000        10.0
#> 2    1 2.000002 4.605175 2.000002         8.8
#> 3    2 2.000002 4.605175 2.000002         7.6
#> 4    3 2.000002 4.605175 2.000002         6.4
#> 5    4 2.000002 4.605175 2.000002         5.2
#> 6    5 2.000002 4.605175 2.000002         4.0
tail(growth_pred)
#> dynamic_growth prediction
#>   model: huang
#>   secondary: huang_sqrt
#>   rows: 6
#>    time response      lnN     logN temperature
#> 26   25 2.352192 5.416123 2.352192         2.0
#> 27   26 2.352192 5.416123 2.352192         3.6
#> 28   27 2.352192 5.416123 2.352192         5.2
#> 29   28 2.352192 5.416123 2.352192         6.8
#> 30   29 2.352192 5.416123 2.352192         8.4
#> 31   30 2.353249 5.418557 2.353249        10.0

plot(growth_pred$time, growth_pred$response,
  type = "l",
  lwd = 2,
  xlab = "Time",
  ylab = "Predicted log10 count"
)

Predicted dynamic microbial growth curve under a changing temperature profile.

Dynamic growth fitting

Dynamic fitting estimates selected parameters by repeatedly solving the dynamic model and minimizing residual error at the observation times. In practice, it is often wise to estimate only a small number of parameters and keep others fixed, especially when the data set is small.

Here we simulate observations from the dynamic growth prediction and estimate the secondary-model parameter a.

obs_growth <- data.frame(
  time = growth_pred$time[seq(1, nrow(growth_pred), by = 3)],
  logN = growth_pred$response[seq(1, nrow(growth_pred), by = 3)]
)

growth_fit <- fit_dynamic_growth(
  data = obs_growth,
  profile = profile,
  time = "time",
  response = "logN",
  model = "huang",
  secondary = "huang_sqrt",
  start = list(
    logN0 = 2,
    logNmax = 8.8,
    a = 0.07,
    Tmin = 8.91,
    lag = 2
  ),
  estimate = "a",
  fixed = list(
    logN0 = 2,
    logNmax = 8.8,
    Tmin = 8.91,
    lag = 2
  )
)

growth_fit
#> predmicror dynamic fit
#>   type: dynamic_growth
#>   model: huang
#>   secondary: huang_sqrt
#>   response: logN (log10)
#>   SSE: 0.000000000012959
#>   convergence: 0
coef(growth_fit)
#>      logN0    logNmax          a       Tmin        lag 
#> 2.00000000 8.80000000 0.08859978 8.91000000 2.00000000
fit_metrics(growth_fit)
#>   model           type response response_scale  n p          SSE         RMSE
#> 1 huang dynamic_growth     logN          log10 11 1 1.295925e-11 1.085409e-06
#>            MAE         bias          RSE R2 adj_R2   logLik       AIC       BIC
#> 1 7.395177e-07 7.395177e-07 1.138387e-06  1      1 135.4608 -268.9215 -268.5236
#>   converged
#> 1      TRUE

plot(obs_growth$time, obs_growth$logN,
  pch = 16,
  xlab = "Time",
  ylab = "log10 count"
)
lines(growth_fit$prediction$time, growth_fit$prediction$response, lwd = 2)

Observed and fitted dynamic growth values.

The fitted object supports familiar helper methods:

head(fitted(growth_fit))
#> [1] 2.000000 2.000002 2.000002 2.000002 2.000002 2.048006
head(residuals(growth_fit))
#> [1] 0.000000e+00 9.289902e-12 9.289902e-12 9.289902e-12 9.289902e-12
#> [6] 2.370187e-07
head(predmicror_augment(growth_fit))
#>   time     logN  .fitted       .resid .model          .type
#> 1    0 2.000000 2.000000 0.000000e+00  huang dynamic_growth
#> 2    3 2.000002 2.000002 9.289902e-12  huang dynamic_growth
#> 3    6 2.000002 2.000002 9.289902e-12  huang dynamic_growth
#> 4    9 2.000002 2.000002 9.289902e-12  huang dynamic_growth
#> 5   12 2.000002 2.000002 9.289902e-12  huang dynamic_growth
#> 6   15 2.048006 2.048006 2.370187e-07  huang dynamic_growth

Dynamic inactivation prediction

Dynamic inactivation can be simulated with the Weibull-Peleg formulation. In the simplest case, the rate is kept constant.

inact_profile <- dynamic_profile(
  time = c(0, 10),
  temperature = c(60, 60)
)

inact_pred <- predict_dynamic_inactivation(
  profile = inact_profile,
  model = "weibull_peleg",
  secondary = "constant",
  start = list(
    logN0 = 7,
    b = 0.2,
    n = 1
  ),
  times = seq(0, 10, by = 1)
)

head(inact_pred)
#> dynamic_inactivation prediction
#>   model: weibull_peleg
#>   secondary: constant
#>   rows: 6
#>   time response logN log_survival temperature
#> 1    0      7.0  7.0          0.0          60
#> 2    1      6.8  6.8         -0.2          60
#> 3    2      6.6  6.6         -0.4          60
#> 4    3      6.4  6.4         -0.6          60
#> 5    4      6.2  6.2         -0.8          60
#> 6    5      6.0  6.0         -1.0          60
tail(inact_pred)
#> dynamic_inactivation prediction
#>   model: weibull_peleg
#>   secondary: constant
#>   rows: 6
#>    time response logN log_survival temperature
#> 6     5      6.0  6.0         -1.0          60
#> 7     6      5.8  5.8         -1.2          60
#> 8     7      5.6  5.6         -1.4          60
#> 9     8      5.4  5.4         -1.6          60
#> 10    9      5.2  5.2         -1.8          60
#> 11   10      5.0  5.0         -2.0          60

plot(inact_pred$time, inact_pred$response,
  type = "l",
  lwd = 2,
  xlab = "Time",
  ylab = "Predicted log10 count"
)

Predicted dynamic microbial inactivation curve.

Dynamic inactivation fitting

The same framework can fit selected inactivation parameters. Here we estimate the Weibull-Peleg rate parameter b.

obs_inact <- data.frame(
  time = c(0, 5, 10),
  logN = c(7, 6, 5)
)

inact_fit <- fit_dynamic_inactivation(
  data = obs_inact,
  profile = inact_profile,
  time = "time",
  response = "logN",
  model = "weibull_peleg",
  secondary = "constant",
  start = list(
    logN0 = 7,
    b = 0.1,
    n = 1
  ),
  estimate = "b",
  fixed = list(
    logN0 = 7,
    n = 1
  )
)

inact_fit
#> predmicror dynamic fit
#>   type: dynamic_inactivation
#>   model: weibull_peleg
#>   secondary: constant
#>   response: logN (log10)
#>   SSE: 0.000000000000000000000000051711
#>   convergence: 0
coef(inact_fit)
#> logN0     b     n 
#>   7.0   0.2   1.0
fit_metrics(inact_fit)
#>           model                 type response response_scale n p          SSE
#> 1 weibull_peleg dynamic_inactivation     logN          log10 3 1 5.171141e-26
#>           RMSE          MAE          bias          RSE R2 adj_R2   logLik
#> 1 1.312903e-13 1.030287e-13 -1.030287e-13 1.607971e-13  1      1 84.72728
#>         AIC      BIC converged
#> 1 -167.4546 -168.356      TRUE

plot(obs_inact$time, obs_inact$logN,
  pch = 16,
  xlab = "Time",
  ylab = "log10 count"
)
lines(inact_fit$prediction$time, inact_fit$prediction$response, lwd = 2)

Observed and fitted dynamic inactivation values.

Sensitivity analysis

Finite-difference sensitivity analysis helps identify which parameters most influence the model output. The scaled sensitivity coefficients are useful because they are expressed on the response scale.

sens <- dynamic_sensitivity(
  type = "growth",
  profile = profile,
  model = "huang",
  secondary = "huang_sqrt",
  start = list(
    logN0 = 2,
    logNmax = 8.8,
    a = 0.0886,
    Tmin = 8.91,
    lag = 2
  ),
  parameters = c("a", "Tmin"),
  times = seq(0, 30, by = 2)
)

head(sens)
#>   time parameter prediction  sensitivity scaled_sensitivity
#> 1    0         a   2.000000 0.000000e+00       0.000000e+00
#> 2    2         a   2.000002 4.249422e-05       3.764988e-06
#> 3    4         a   2.000002 4.249422e-05       3.764988e-06
#> 4    6         a   2.000002 4.249422e-05       3.764988e-06
#> 5    8         a   2.000002 4.249422e-05       3.764988e-06
#> 6   10         a   2.000002 4.249422e-05       3.764988e-06

plot(sens$time, sens$scaled_sensitivity,
  type = "n",
  xlab = "Time",
  ylab = "Scaled sensitivity"
)
for (p in unique(sens$parameter)) {
  z <- sens[sens$parameter == p, ]
  lines(z$time, z$scaled_sensitivity, lwd = 2)
}
legend("topleft", legend = unique(sens$parameter), lty = 1, lwd = 2, bty = "n")

Scaled sensitivity coefficients for selected dynamic growth parameters.

Practical notes

Dynamic fitting can be computationally more demanding than fitting explicit isothermal models because each objective-function evaluation requires solving the dynamic trajectory. For small data sets, avoid estimating too many parameters at once. Use fixed to keep parameters at known or literature-supported values and estimate only the parameters that are informed by the available observations.