bioinactivation: Software for modelling microbial inactivation

Introduction

R-package bioinactivation implements functions for the modelization of microbial inactivation. These functions are based on published inactivation models commonly used both in academy and industry. bioinactivation can be downloaded from CRAN. Once installed, the package is loaded with the command.

library(bioinactivation)

The bioinactivation package allows to predict the inactivation process of a population of microorganisms exposed to isothermal conditions (temperature remains constant with time) or non isothermal (temperature changes with time, dynamic conditions). Moreover, it can be used to adjust the parameters of different inactivation models to both isothermal and dynamic experimental data. All the functions contained in bioinactivation are written in R.

bioinactivation takes advantage of the capabilities of the packages deSolve (Soetaert, Petzoldt, and Setzer 2010a) and FME (Soetaert, Petzoldt, and Setzer 2010b) to make predictions of the inactivation process and to fit model parameters to experimental data. The inactivation models most commonly used in the industry have been implemented in this package:

Bigelow model (Bigelow, 1921),
Weibull-Peleg model (Peleg and Cole, 1998),
Weibull-Mafart model (Mafart et al., 2002),
Geeraerd’s model (Geeraerd et al., 2005).
Arrhenius model (Box, 1960).

The bioinactivation package provides five functions for the modelization of microbial inactivation:

predict_inactivation() predicts the inactivation process of a population for some given temperature conditions, isothermal or dynamic.
fit_isothermal_inactivation() adjusts the parameters of a selected inactivation model to a set of isothermal inactivation experiments.
fit_dynamic_inactivation() fits the parameters of a chosen inactivation model to a set of dynamic inactivation experiments using non-linear regression.
fit_inactivation_MCMC() adjusts an inactivation model to a set of dynamic inactivation experiments using a Markov Chain Monte Carlo algorithm.
predict_inactivation_MCMC() generates prediction intervals for isothermal or non-isothermal inactivation processes using a Monte Carlo method.

Moreover, two functions with meta information of the inactivation models implemented are included:

get_model_data(), provides meta information about the models used for prediction and adjustment of non-isothermal experiments.
get_isothermal_model_data(), provides meta information about the models for fitting of isothermal experiments.

Furthermore, four data sets mimicking the results of an isothermal inactivation experiment (isothermal_inactivation and laterosporus_iso) and a dynamic experiment (dynamic_inactivation and laterosporus_dyna) are included in the package.

bioinactivation does not require the user to introduce the data in any particular units system. The units system used for input is kept on the output. For this reason, units are not shown in any figure within this document.

Example data sets included in the package

The example data set of an isothermal experiment

The bioinactivation package includes the isothermal_inactivation data set, which imitates the data obtained during a set of isothermal inactivation experiments. This data set can be loaded as:

data(isothermal_inactivation)

This data set supplies experiments made at six different temperature levels, with a variable number of observations per experiment. In total, the data set contains 68 observations of 3 variables:

The first variable (time) provides the experiment time at which the observation was made.
The second column (temp) contains the temperature of the experiment at a given time.
The logarithmic difference between the initial concentration and the one at a time-temperature combination is given by the third column (log_diff).

The following plot illustrates the data contained in the data set.

library(ggplot2)
ggplot(data = isothermal_inactivation) +
    geom_point(aes(x = time, y = log_diff, col = as.factor(temp))) +
    ggtitle("Example dataset: isothermal_inactivation")

The example data set of a dynamic experiment

The dynamic_inactivation data set contained within the bioinactivation package provides a data set imitating the data obtained during a non-isothermal inactivation experiment. This data set can be loaded as:

data(dynamic_inactivation)

This example data set contains 3 variables and 19 observations. The first column, time, reflects the time point at which the observation was taken. The second and third ones provide the number of microorganism (N) and the temperature (temperature) recorded.

The following plots depicts the data contained in the data set.

ggplot(data = dynamic_inactivation) +
    geom_point(aes(x = time, y = log10(N))) +
    ggtitle("Example dataset: dynamic_inactivation. Observations")

ggplot(data = dynamic_inactivation) +
    geom_point(aes(x = time, y = temperature))  +
    ggtitle("Example dataset: dynamic_inactivation. Temperature profile")

The example data set of an isothermal study of Bacillus laterosporus

bioinactivation contains a data set with the isothermal inactivation observed in a population of B. laterosporus. This data set can be loaded as:

data(laterosporus_iso)

This data set contains 52 observations and 3 different variables:

The first variable (time) provides the experiment time at which the observation was made.
The second column (temp) contains the temperature of the experiment at a given time.
The logarithmic difference between the initial concentration and the one at a time-temperature combination is given by the third column (log_diff).

The following plot illustrates the data contained in the data set.

ggplot(data = laterosporus_iso) +
    geom_point(aes(x = time, y = log_diff, col = as.factor(temp))) +
    ggtitle("Example dataset: laterosporus_iso")

The example data set of a non-isothermal study of laterosporus

The data set laterosporus_dyna included in bioinactivation contains the results of a series of non-isothermals studies performed on a population of B. laterosporus. This data set can be loaded as:

data(laterosporus_dyna)

This data set contains 20 observations with 3 variables:

time includes the time spent betwen the beginning of the experiment and the measurement.
temp is the temperature recorded at the moment of the observation.
logN is the logarithmic number of microorganism observed.

The following plots depict the data within the data set.

ggplot(data = laterosporus_dyna) +
    geom_point(aes(x = time, y = logN)) +
    ggtitle("Example dataset: laterosporus_dyna. Observations")

ggplot(data = laterosporus_dyna) +
    geom_point(aes(x = time, y = temp))  +
    ggtitle("Example dataset: laterosporus_dyna. Temperature profile")

Inactivation models implemented

The inactivation models most commonly used in both academy and industry for the modelization of microbial inactivation have been implemented in the bioinactivation package:

Bigelow model (Bigelow, 1921),
Weibull-Peleg model (Peleg and Cole, 1998),
Weibull-Mafart model (Mafart et al., 2002),
Geeraerd model (Geeraerd et al., 2005).
Arrhenius model (Box, 1960).

A brief mathematical description of each one of them is presented through this section.

Bigelow model for isothermal data

Bigelow’s model (Bigelow, 1921) assumes that the inactivation process follows a first order kinetics. Therefore, a linear relation exists between the logarithm of the number of microorganisms (\(N\)) and the time (\(t\)).

\[ \mathrm{log}_{10}(N) = \mathrm{log}_{10}(N_0) - \frac{1}{D(T)}t \]

The parameter \(D(T)\) (also named D-value) equals the negative inverse of the slope of the line. Its relationship with temperature (\(T\)) is described by the equation

\[ \mathrm{log}_{10}D(T) = \mathrm{log}_{10}D_R + \frac{T_R-T}{z} \]

where \(T_R\) is a reference temperature and \(D_R\) the D-value at this temperature. The parameter \(z\) is usually named z-value and provides an indication of the sensitivity of the microorganism to temperature variations.

Bigelow model for nonisothermal data

Bigelow’s model for nonisothermal data is developed from the isothermal one by calculating the first derivative of its constitutive equation with respect to time.

\[ \frac{d\mathrm{log}_{10}(N)}{dt} = - \frac{1}{D(T)} \]

In this model, the evolution of \(D(T)\) with temperature is ruled by the expression used for the isothermal model:

\[ \mathrm{log}_{10}D(T) = \mathrm{log}_{10}D_R + \frac{T_R-T}{z} \]

Weibullian models

Weibullian models assume that the the inactivation process can be viewed as a failure phenomenon. According to this hypothesis, the time that each microorganism within the population can sustain some environmental conditions follows a Weibull-type probability distribution. Two models of this family have been implemented: Weibull-Peleg (Peleg and Cole, 1998) and Weibull-Mafart (Mafart et al., 2002).

Weibull-Peleg model for isothermal inactivation

The Weibullian model presented by Peleg and Cole (1998) for isothermal cases follows the equation:

\[ \log_{10}S(t) = - b \cdot t^n \]

where \(S = N/N_0\) represents the fraction of survivors and \(t\) the time, and \(b\) and \(n\) are the two parameters of the model.

In the Weibull-Peleg model the parameter \(n\) is assumed to be temperature independent, while the dependency of \(b\) with temperature is modeled as

\[ b(T) = \mathrm{ln} \big\{ 1 + \mathrm{exp} \big[ k_b(T - T_c) \big] \big\} \]

This equation is based on the definition of a critical temperature (\(T_c\)). For temperatures (\(T\)) much lower than \(T_c\), \(b(T) \simeq 0\) and the inactivation is insignificant. For temperatures close to \(T_c\), there is an exponential growth of \(b(T)\) with respect to the temperature. For temperatures much higher than the critical one, a linear relation exists between \(T\) and \(b(T)\) with a slope \(k_b\). This evolution is illustrated in the following plot:

Weibull-Peleg model for non-isothermal inactivation

The model described by Peleg and Cole (1998) for dynamic inactivation has been implemented in the bioinactivation package. This model is based on the constitutive equation of the Weibull-Peleg model (Peleg and Cole, 1998) for isothermal inactivation:

\[ \log_{10}S(t) = -b(T) \cdot t^n \]

The inactivation rate for a constant temperature can be calculated as the first derivative with respect to time of this equation:

\[ \Bigg| \frac{d \log_{10}S(t)}{dt} \Bigg|_{T=const} = -b(T) \cdot n \cdot t^{n-1} \]

The time required to achieve any survival ratio, \(t^\star\), can be calculated from the model equation as:

\[ t^\star = \Bigg| \frac{-\mathrm{log}_{10}S(t)}{b(T)} \Bigg|^{1/n} \]

Substituting this expression into the rate equation, the differential constitutive equation of the model is obtained as:

\[ \Bigg| \frac{d \log_{10}S(t)}{dt} \Bigg| = -b(T) \cdot n \cdot \Bigg| \frac{ -\mathrm{log}_{10}S(t) }{b(T)} \Bigg|^\frac{n-1}{n} \]

where \(b(T)\) follows the same expression as for the isothermal model:

\[ b(T) = \mathrm{ln} \big\{ 1 + \mathrm{exp} \big[ k(T - T_c) \big] \big\} \]

Weibull-Mafart model for isothermal inactivation

The Weibullian model presented by Mafart et al. (2002) for isothermal cases is described by the equation:

\[ \mathrm{log}_{10}S(t) = - \Big( \frac{t}{\delta(T)} \Big)^\beta \]

where \(S = N/N_0\) represents the fraction of survivors and \(t\) is the time. \(\delta\) and \(\beta\) are, respectively, the rate and shape factors of the Weibull distribution.

In this model, \(\beta\) is considered to be temperature independent and \(\delta\) follows and exponential relation with temperature (\(T\)):

\[ \delta(T) = \frac {\delta_{ref}} {10^ {\big( T-T_{ref}/z \big)} } \]

where \(T_{ref}\) is a reference temperature, \(\delta_{ref}\) is the value of \(\delta\) at this reference temperature and \(z\) describes the sensitivity of the microorganism to temperature variations.

Weibull-Mafart model for nonisothermal inactivation

The dynamic Weibull-Mafart model is developed from the isothermal one. The isothermal rate of logarithmic reduction at each time is calculated as the first derivative of the constitutive equation with respect to time:

\[ \frac{d\mathrm{log}_{10}S(t)}{dt} = -\beta \frac{t^{\beta-1}}{\delta(T)^\beta} \]

As well as for the isothermal model, an exponential relation between \(\delta(T)\) and temperature is assumed:

\[ \delta(T) = \frac {\delta_{ref}} {10^ {\big( T-T_{ref}/z \big)} } \]

Geeraerd model

The inactivation model described by Geeraerd et al. (2005) is based on the assumption that the inactivation process follows a first order kinetics with rate \(k\).

\[ \frac {dN}{dt} = - k \cdot N \]

However, the parameter \(k\) is corrected in this model with two parameters to include for shoulder (\(\alpha\)) and tail (\(\gamma\)) effects:

\[ k = k_{max} \cdot \alpha \cdot \gamma \]

where \(k_{max}\) stands for the maximum inactivation rate of the microorganism at the current environmental conditions.

A secondary model equivalent to Bigelow’s is used to describe the evolution of \(k_{max}\) with temperature (\(T\)):

\[ k_{max} = \frac{1}{D(T)}\] \[ \mathrm{log}_{10}D(T) = \mathrm{log}_{10}D_R + \frac{T_R-T}{z} \]

where \(D(T)\), \(T_R\) and \(z\) are the D-value, reference temperature and z-value respectively.

The shoulder parameter is defined as

\[ \alpha = \frac {1}{1+C_c} \]

where \(C\) is a dummy component that limits the microbial inactivation. It is assumed that this component also follows a first order kinetics with rate \(k_{max}\):

\[ \frac {dC_c}{dt} = - k_{max} \cdot C_c \]

The tail parameter is defined by the equation

\[ \gamma = 1 - \frac{N_{res}}{N} \]

where \(N_{res}\) stands for the residual number of microorganism after the treatment.

Arrhenius inactivation model

Several versions of the Arrhenius model are available in the literature. In bioinactivation, a log-linear primary inactivation model:

\[ \frac{ dN}{dt} = -k(T)\cdot N \]

where \(k(T)\) represents the specific inactivation rate at temperature (\(T\)). The reparametrization proposed by Box (1960) was used as secondary model:

\[ k(T) = k_{ref} \exp{ \left[- \frac{E_a}{R} \left( \frac{1}{T} - \frac{1}{T_{ref}} \right) \right]}\]

where \(k_{ref}\) is the value of \(k\) at the reference temperature \(T_{ref}\). \(E_a\) is the activation energy and \(R = 8.31 J\cdot mol^{-1} K^{-1}\), is the universal constant gas.

Prediction of microbial inactivation

The bioinactivation package includes capabilities for the prediction of the inactivation process of a microorganism under either isothermal or dynamic conditions. This is accomplished through the predict_inactivation() function, which solves the differential equation of the model using the ode function from the deSolve package (Soetaert, Petzoldt, and Setzer 2010).

predict_inactivation() requires five input arguments:

simulation_model,
times,
parms,
temp_profile,
... and
tol0 (optional).

simulation_model is a string identifying the model to be used. A list of the available models can be obtained by calling the function get_model_data() without any argument.

get_model_data()

## [1] "Bigelow"   "Mafart"    "Metselaar" "Geeraerd"  "Peleg"     "Arrhenius"

In this example Geeraerd’s model will be used, so the value Geeraerd will be assigned.

example_model <- "Geeraerd"

times is a numeric vector which specifies the values of time when results will be returned. It does not define the points at which the numerical integration is done, so changing it does not affect the values calculated at individual time points. An arbitrary equally spaced vector will be used for this example.

times <- seq(0, 4.5, length=100)

parms is a named numeric vector specifying the values of the model parameters, as well as the initial values of the model variables. A call to the function get_model_data() with a valid model identifier as an input argument provides several data concerning this model as a list. Namely, its parameters and variables.

model_data <- get_model_data(example_model)
print(model_data$parameters)

## [1] "D_R"      "z"        "N_min"    "temp_ref"

print(model_data$variables)

## [1] "N"   "C_c"

The input argument parms must provide values for every degree of freedom of the model, i.e. model parameters and variables. The names of the model parameters must be identical to those provided by the function get_model_data(). On the other hand, the names assigned to the initial values of the model variables must be equal to the ones provided by get_model_data() plus a character "0". For instance, the initial value for the variable named N needs to be labeled "N0".

model_parms <- c(D_R = 3,
                 z = 10,
                 N_min = 1e2,
                 temp_ref = 100,
                 N0 = 1e5,
                 C_c0 = 1e1
                 )

The last input variable required for the function is a data frame defining a discrete temperature profile. Values of the temperature at time points not provided will be approximated by linear interpolation. In case a temperature value outside the time range specified is requested, the temperature of the closest extreme will be used.

Predictions for isothermal conditions can also be calculated using this function. In this case, a constant temperature profile has to be defined.

A temperature profile with a constant slope ranging between 70 and 120ºC will be used for the prediction.

temperature_profile <- data.frame(time = c(0, 4.5),
                                  temperature = c(70, 120))

plot(temperature_profile, type = "l")
title("Example temperature profile")

The inactivation models based on the Weibull distribution are singular for \(t=0\). For this reason, every value of time in times lower than the optional argument tol0 is changed for tol0. By default, tol0 = 1e-05.

Finally, the argument ... allows to pass adidtional arguments to the ode() function of the deSolve package.

The prediction is calculated by calling the function predict_inactivation() with the parameters described before. The optional argument will not be modified.

prediction_results <- predict_inactivation(example_model, times, model_parms, temperature_profile)

The value returned by the function is a list of class SimulInactivation. This list has four entries:

model, a string with the identifier of the model used for the simulation.
model_parameters, a named numeric list with the model parameters used.
temp_approximations, a list with two entries, each providing the functions used to approximate the temperature at each time point.
simulation a data frame with the results of the simulation.

The entry simulation provides a data frame with the calculated results.

head(prediction_results$simulation)

##         time         N       C_c     logN         S          logS
## 1 0.00001000 100000.00 10.000000 5.000000 1.0000000  0.000000e+00
## 2 0.04545455  99999.66  9.999630 4.999999 0.9999966 -1.460606e-06
## 3 0.09090909  99999.29  9.999214 4.999997 0.9999929 -3.101738e-06
## 4 0.13636364  99998.86  9.998746 4.999995 0.9999886 -4.945336e-06
## 5 0.18181818  99998.39  9.998222 4.999993 0.9999839 -7.013202e-06
## 6 0.22727273  99997.85  9.997633 4.999991 0.9999785 -9.337306e-06

The first column contains the times at which results are output are given. They are equal to the values provided by the times input variable, although sorted. In case there were some repeated values in the input argument times, results are output only once.

The following columns provide the results obtained for the variables of the models. In this case, as Geeraerd’s model has been used, N and C_c. The last three columns provide the values of \(\mathrm{log}_{10}(N)\), \(S\) and \(\mathrm{log}_{10}(S)\).

The SimulInactivation object includes an S3 method for the plotting of its results. This method is able to generate the plot using ggplot2 (default) or base R. The ggplot2 graphic can be generated by calling the plot() function with the instance of SimulInactivation as argument. The function returns an instance of ggplot which can be represented by calling the function print().

p <- plot(prediction_results)
print(p)

Instead of directly generating the plot, the function returns the ggplot instance to ease further editting. For instance, the format and x-label can be easily modified.

library(ggplot2)
p + theme_light() + xlab("time (min)")

Note that if the returned object is not assigned to any variable, the plot is automatically printed.

On the other hand, the version of the plot in base R can be generated including the argument make_gg = FALSE in the call to plot().

plot(prediction_results, make_gg = FALSE,
     xlab = "Time (min)", ylab = "logN")

The implementation of the Geeraerd model allows to model inactivation process without shoulder and/or tail effects. For the model without shoulder, assign zero to the parameter C_c0:

parms_no_shoulder <- c(D_R = 3,
                       z = 10,
                       N_min = 100,
                       temp_ref = 100,
                       N0 = 100000,
                       C_c0 = 0
                       )


prediction_no_shoulder <- predict_inactivation(example_model, times, parms_no_shoulder, temperature_profile)

plot(prediction_no_shoulder)

Although the effect is not fully clear with the temperature profile used due to the increasing temperature.

The model without tail effect can be obtained by assigning zero to N_min.

parms_no_tail <- c(D_R = 3,
                   z = 10,
                   N_min = 0,
                   temp_ref = 100,
                   N0 = 100000,
                   C_c0 = 100
                   )


prediction_no_tail <- predict_inactivation(example_model, times, parms_no_tail, temperature_profile)

plot(prediction_no_tail)

Note that this operation causes some warnings. The differential equation of the system is solved for the variable N. For values of this variable very close to zero (of an order of magnitude lower than \(10^{-5}\)), the numerical error may cause the calculated value of N to be lower than zero. Afterwards, when the logarithmic count of microorganisms is calculated, warnings are generated for these values. However, this problem only occurs for values of N below the level of detection in normal laboratory conditions and are not relevant in most cases.

The prediction interval can be added to the plot by setting the parameter plot_temp to TRUE.

plot(prediction_results, plot_temp = TRUE)

## Warning: Removed 1 rows containing missing values (geom_path).

If this is the case, the temperature is plotted as a solid line and the prediction curve as a dashed line.

Fitting of isothermal data

The bioinactivation package includes functionality for the fitting of inactivation models to isothermal data. The function fit_isothermal_inactivation() makes use of the nls() function from the stats package to fit the model parameters using non-linear regression.

The fit_isothermal_inactivation() function requires the definition of five input parameters:

model_name,
death_data,
starting_point,
known_params and
adjust_log (optional).

The experimental data to fit is provided by the argument death_data, which is a data frame. It must contain at least three columns providing the times when the observations were made (time), the temperature of the experiment (temp) and the logarithmic difference observed with respect to the original population (log_diff). The example data set isothermal_inactivation, which already has this format, is used for this demonstration.

data(isothermal_inactivation)

The model to be used for the adjustment is defined by the character variable model_name. A list of the valid model keys for isothermal adjustment can be obtained by calling the function get_isothermal_model_data() without any arguments.

get_isothermal_model_data()

## [1] "Mafart"    "Peleg"     "Bigelow"   "Arrhenius" "Metselaar"

Note that Geeraerd’s model is not available for isothermal adjustment. The reason for this is the absence of a secondary model that explains the relation between \(N_{res}\) and temperature, making impossible the adjustment using non-linear regression. A possible alternative for the adjustment of Geeraerd’s model to individual isothermal experiments is the use of the function fit_dynamic_inactivation() with a constant temperature profile and fixing the value of \(z\) to an arbitrary value.

Bigelow model, whose key is "Bigelow", is used in this example.

model_name <- "Bigelow"

The names of the parameters of an inactivation model can be obtained by calling the function get_isothermal_model_data() with the model key as argument. Note that for the fitting of isothermal experiments the meta data of the models is gathered calling get_isothermal_model_data() while for making predictions and fitting dynamic experiments get_model_data() is called.

model_data <- get_isothermal_model_data(model_name)
model_data$params

## [1] "z"        "D_R"      "temp_ref"

Therefore, the isothermal Bigelow model requires the definition of three parameters: z, D_ref and ref_temp. The fit_isothermal_inactivation() function allows to distinguish between known and unknown (i.e., to be adjusted) model parameters. The known model parameters are set by the input argument known_parameters, which must be a list. The reference temperature is usually fixed to a value close to the experimental ones, we will set it to 100ºC in this example.

known_params = list(temp_ref = 100)

According to the requirements of the nls function, initial points are required for the unknown model parameters. These initial values are provided by the input argument starting_point. This variable needs to be a list with its names identical to those retrieved using get_isothermal_model_data().

starting_point <- list(z = 10,
                       D_R = 1.5
                       )

Definition of unrealistic starting points might result in nls raising an error. Hence, it is important to define values relatively close to the solution. This can be done from experience, consulting the literature or testing isothermal predictions with the predict_inactivation() function.

The fit_isothermal_inactivation() functions can base the adjustment on the minimization of the error of either the logarithmic count (\(\mathrm{log}_{10}(N)\)) or the total number of microorganisms (\(N\)). This is controled by the boolean variable adjust_log. In case it is TRUE (default case), the adjustment is based on the error of the logarithmic count. Ohterwise, it is based on the error in the total number of microorganism. In this example, the default option will be used.

The adjustment is made by calling the fit_isothermal_inactivation() function:

iso_fit <- fit_isothermal_inactivation(model_name, isothermal_inactivation,  
                                       starting_point, known_params)

This function returns a list of class IsoFitInactivation. This object contains the results of the adjustment, as well as the statistical information of the adjustment process. This information is divided in four entries:

nls,
parameters,
model and
data.

The information of the adjustment process is provided under the entry nls, where the object returned by the nls function is saved. This argument provides extensive statistical information regarding the adjustment process. For more information, consult the help page for nls.

summary(iso_fit$nls)

## 
## Formula: log_diff ~ Bigelow_iso(time, temp, z, D_R, temp_ref)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## z     7.7605     0.2442   31.78   <2e-16 ***
## D_R   2.9805     0.1422   20.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4801 on 66 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 1.126e-06

vcov(iso_fit$nls)  # Calculates variance-covariance matrix {stats}

##               z         D_R
## z    0.05963167 -0.02944040
## D_R -0.02944040  0.02023196

confint(iso_fit$nls)  # Calculates confidence intervals {stats}

## Waiting for profiling to be done...

##         2.5%    97.5%
## z   7.285726 8.272271
## D_R 2.720058 3.296796

The entry parameters contains the values of both the parameters fixed through the input argument known_params and those adjusted by non linear regression. These values are saved as a list.

iso_fit$parameters

## $temp_ref
## [1] 100
## 
## $z
## [1] 7.760496
## 
## $D_R
## [1] 2.980496

The key of the model used for the adjustment is provided under the header model.

iso_fit$model

## [1] "Bigelow"

The data used for the adjustment is saved under the header data.

head(iso_fit$data)

##   time temp   log_diff
## 1    0  100  0.0000000
## 2    1  100 -0.2086205
## 3    2  100 -0.8500512
## 4    4  100 -0.8077948
## 5    6  100 -1.7080181
## 6    8  100 -2.0522733

The class IsoFitInactivation implements S3 methods for plot(). This method generates a comparison plot between experimental data and predicted values for every single experimental temperature. The title of each of the plots specifies the temperature of the experiments.

plot(iso_fit, ylab = "Number of logarithmic reductions",
     xlab = "Time (min)")

The plot can also be generated using ggplot2. To do that, the argument make_gg = TRUE must be added to the call to the plot() function.

plot(iso_fit, make_gg = TRUE) +
    ylab("Number of logarithmic reductions") +
    xlab("Time (min)")

Fitting of dynamic experiments using non-linear regression

The bioinactivation package implements the function fit_dynamic_inactivation(), which permits the adjustment of inactivation models to experimental data with time dependent temperature profiles. This is accomplished using the function modFit() of the FME package (Soetaert, Petzoldt, and Setzer 2010b), which adjusts the model using non-linear regression.

The function fit_dynamic_inactivation() requires the definition of eight input arguments:

experiment_data,
simulation_model,
temp_profile,
starting_points,
upper_bounds,
lower_bounds,
known_params,
...,
minimize_log (optional) and
tol0 (optional).

The experimental data to fit is provided by the variable experiment_data. It has to be a data frame with at least a column named time, which provides the times when the measurements were taken, and another one named N, with the number of microorganisms counted at each time. The function requires the initial number of microorganism \(N_0\) to be the same for each experiment. The data set dynamic_inactivation included in the bioinactivation package is used for this demonstration.

data(dynamic_inactivation)

The inactivation model is specified by the character variable simulation_model. A list of the available models for dynamic adjustment can be obtained by calling the function get_model_data() without any input argument.

get_model_data()

## [1] "Bigelow"   "Mafart"    "Metselaar" "Geeraerd"  "Peleg"     "Arrhenius"

The Weibull-Mafart model is used in this example.

simulation_model <- "Mafart"

The temperature profile of the experiment is given by the data frame temp_profile, which provides the values of the temperature at discrete time points. Temperatures at time points not included in temp_profile are calculated by linear interpolation. The dynamic_inactivation data set includes the temperature of each observation. Hence, the temperature profile defined by this observations is used for the adjustment.

dummy_temp <- data.frame(time = dynamic_inactivation$time,
                         temperature = dynamic_inactivation$temperature)

Although the temperature profile defines several different values of temperature for each time value, only the one calculated by linear interpolation will be used for the resolution of the differential equation.

The classification of every degree of freedom of the model (model parameters and initial values of the variables) as either known or to be adjusted is required for the adjustment. The degrees of freedom of the model selected can be obtained by calling get_model_data() with the model key as argument.

model_data <- get_model_data(simulation_model)
model_data$parameters

## [1] "delta_ref" "temp_ref"  "z"         "p"

model_data$variables

## [1] "N"

The degrees of freedom (parameters or inial values of variables) known are defined through the named list known_parameters. The names of known model parameters need to be identical to those provided by the function get_model_data(). The names of known initial values must be equal to the variable name provided by get_model_data() plus a caracter zero ("0"). In the specific casee of the initial count, if the parameters are called N0, the model will estimate the microbial count. If they are named logN0, the log-microbial count will be estimated. In this example, solely the parameter temp_ref will be considered as known. This parameter is usually considered fixed for microbial inactivation.

known_params = c(temp_ref = 90)

The rest of the degrees of freedom are adjusted to the experimental data. The fitting algorithm requires the definition of starting points for the adjustment of unknown parameters, as well as upper and lower bounds. This information is passed by the input arguments starting_points, upper_bounds and lower_bounds, which share format with the input argument known_params. In case of unbounded variables, +Inf or -Inf must be assigned to the desired entries of upper_bounds or lower_bounds.

starting_points <- c(delta_ref = 10,
                     p = 1,
                     z = 10,
                     logN0 = 5)
upper_bounds <- c(delta_ref = 20,
                  p = 2,
                  z = 20,
                  logN0 = Inf)

lower_bounds <- c(delta_ref = 5,
                  p = .5,
                  z = 5,
                  logN0 = 4)

Unrealistic initial points for the adjustment will cause the fitting algorithm to fail. Realistic initial points can be obtained from experience, literature or by testing different predictions using predict_inactivation().

It is advisable to use reasonable lower and upper bounds, so degrees of freedom do not become negative or unrealistically large.

The adjustment can be based on the minimization of the error of the total count of microorganism (\(N\)) or on the logarithmic count (\(\mathrm{log}_{10}N\)). This distinction is made by the boolean input argument minimize_log. In case it is TRUE, the error of the logarithmic count is minimized. By default, this parameter is set to TRUE. It will not be modified.

Further arguments can be passed to the modFit() function through the ... argument.

The adjustment is made by calling the function fit_dynamic_inactivation():

dynamic_fit <- fit_dynamic_inactivation(dynamic_inactivation, simulation_model, dummy_temp,  
                                        starting_points, upper_bounds, lower_bounds,  
                                        known_params)

The function returns a list of class FitInactivation with three entries:

fit_results,
data and
best_prediction.

fit_results provides the instance of modFit returned by the modFit() function. This variable contains broad information regarding the adjustment process. For more information, refer to the help page of the function.

summary(dynamic_fit$fit_results)

## 
## Parameters:
##            Estimate Std. Error t value Pr(>|t|)    
## delta_ref  9.286270   1.420749   6.536 1.34e-07 ***
## p          0.999935   0.007231 138.291  < 2e-16 ***
## z         10.777956   0.941328  11.450 1.48e-13 ***
## logN0      4.854266   0.131073  37.035  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5554 on 36 degrees of freedom
## 
## Parameter correlation:
##           delta_ref       p       z    logN0
## delta_ref   1.00000 -0.0971 -0.7540  0.06162
## p          -0.09710  1.0000  0.6869 -0.30026
## z          -0.75402  0.6869  1.0000 -0.38056
## logN0       0.06162 -0.3003 -0.3806  1.00000

The header data contains a data frame with the experimental data used for the fit.

head(dynamic_fit$data)

##      time     logN
## 1 0.00001 5.025306
## 2 0.37000 5.110590
## 3 0.74000 5.037426
## 4 1.11000 4.832509
## 5 1.30000 4.255273
## 6 1.80000 2.602060

Under best_prediction, an instance SimulInactivation with the results of the best prediction is saved. For more information about this instance, refer to the documentation of predict_inactivation().

The parameters of the model calculated are saved under the entry model_parameters of the SimulInactivation instance.

dynamic_fit$best_prediction$model_parameters

## $delta_ref
## [1] 9.28627
## 
## $p
## [1] 0.9999346
## 
## $z
## [1] 10.77796
## 
## $temp_ref
## [1] 90
## 
## $N0
## [1] 71493.48

The FitInactivation object includes S3 methods for the creation of a comparison plot between the experimental results and the prediction made by the model parameters adjusted. Two different versions of this function are implemented, one in ggplot2 (default) and another one in base R. By default, the plot ggplot2 implementation is called, which returns an object of class ggplotwith the graphic generated. This allows for further editting:

p <- plot(dynamic_fit)
p + theme_light()

The temperature profile can be added to the plot by setting the argument plot_temp to TRUE.

plot(dynamic_fit, plot_temp = TRUE)

In this case, the temperature profile is plotted as a solid line and the survivor curve as a dashed line.

The base R graph can be generated including the argument make_gg = FALSE.

plot(dynamic_fit, make_gg = FALSE,
     xlab = "Time (min)", ylab = "logN")

Fitting of dynamic experiments using Markov Chain Monte Carlo methods

The function fit_inactivation_MCMC() is able to wrap the modMCMC() function from the FME package (Soetaert, Petzoldt, and Setzer 2010b). Therefore, inactivation models can be adjusted to dynamic experiments using a Markov Chain Monte Carlo algorithm.

The function fit_inactivation_MCMC() has ten input arguments:

experiment_data,
simulation_model,
temp_profile,
starting_points,
upper_bounds,
lower_bounds,
known_params,
...,
minimize_log and
tol0.

Every input argument except for ... are identical in meaning and format to those defined for fit_dynamic_inactivation(), enabling for easy switch between both adjustment methods. The only differents is that, in this case, the ... arguments allows to send further arguments to the modMCMC() function.

In order to highlight the similitudes between this function and fit_dynamic_inactivation(), the input arguments shared by both functions will be reused in this example. An additional argument will be passed to the modMCMC() function, the niter argument. This argument defines the number of iterations of the Markov Chain Monte Carlo algorithm for the adjustment.

set.seed(82619)

MCMC_fit <- fit_inactivation_MCMC(dynamic_inactivation, simulation_model, dummy_temp,  
                                        starting_points, upper_bounds, lower_bounds,  
                                        known_params, niter = 50)

## number of accepted runs: 9 out of 50 (18%)

fit_inactivation_MCMC() returns a list of class FitInactivationMCMC with three entries:

modMCMC,
best_prediction and
data.

modMCMC provides the instance of modMCMC returned by the function modMCMC(). This object includes broad information about the adjustment process. For more information, refer to the help of modMCMC().

summary(MCMC_fit$modMCMC)

##       delta_ref          p         z     logN0
## mean  8.2942017 1.04414259 11.538390 4.9446457
## sd    0.9378426 0.08807178  1.213006 0.1420907
## min   6.8341785 0.91978646 10.000000 4.6093818
## max  10.0000000 1.27024272 13.779134 5.2315942
## q025  7.4742757 0.96842000 10.606819 4.9310982
## q050  8.6457957 1.00000000 11.004503 4.9310982
## q075  9.0285319 1.08863485 12.245353 5.0059541

best_prediction contains a list of class SimulInactivation with the prediction of the adjusted parameters. This object has already been described in the description of the function predict_inactivation(). Finally, the data used for the fit is saved under data.

The FitInactivationMCMC instance includes a S3 implementation of plot() for comparing the prediction against the data. Two different versions of this function are implemented, one in ggplot2 (default) and another one in base R. The default version returns an instance of ggplot with the graphic generated, allowing further edition.

p <- plot(MCMC_fit)
p + xlab("Time (min)")