bioOED: Optimum Experiment Design for Microbial Inactivation

Alberto Garre, Jose Lucas Peñalver Soto, Pablo S. Fernandez, Jose A. Egea

2018-04-17

Introduction

The mathematical models commonly used for the description of microbial inactivation contain parameters whose exact value, at the time, cannot be analytically determined. Hence, they have to be estimated using experiments. This type of procedures do not provide an exact value of the model parameters. Instead, they estimate their probability distribution.

The level of uncertainty strongly depends on the experiment performed. For instance, Grijspeerd and Vanrrolleghem (1999) are able to half the estimated standard deviation of the model parameters for microbial growth using an optimum experiment design. This results in a more accurate description of the process.

The R-package bioOED implements functions for the design of optimum experiments, as well as for the calculation of local sensitivities, of microbial inactivation processes. This package can be downloaded form the Comprehensive R Archive Network (CRAN). Once installed, it can be loaded as

library(bioOED)

bioOED adds the following functions to the namespace:

The simulations of microbial inactivation performed by bioOED are based on the functions implemented in the bioinactivation (Garre et al., 2016) package. Hence, only the inactivation models implemented in this package (Bigelow, Peleg, Mafart and Geeraerd) are available in bioOED. The calculation of local sensitivities is based on the features implemented in the FME package (Soetaert et al., 2010).

Fisher Information Matrix

The Fisher Information Matrix (\(FIM\)) for a variable \(y\) which depends on a vector of parameters (\(\theta\)) at a series of time points (\(t_i\)) is defined by the following equation:

\[ FIM = \sum_i \left( \frac{\partial y}{\partial \theta} \left( t_i \right) \right)^T \cdot Q_i \cdot \left( \frac{\partial y}{\partial \theta} \left( t_i \right) \right) \]

where \(Q_i\) is a matrix of weights. In this document, \(Q_i\) will be considered as the identity matrix.

The FIM is relevant for the OED because its inverse is an estimator of the covariance matrix of the error. Hence, an OED can be made by optimizing certain measurements of the FIM. bioOED implements two of this criteria: D and modified E. The D-criterion consists on the maximization of the determinant of the FIM. It is equivalent to the minimization of the volume of the confidence ellipsoids. The modified E criterion is based on the minimization of the ratio between the maximum and minimum eigenvalue of the FIM (\(\mathrm{abs} \frac{\lambda_{max}}{\lambda_{min}}\)).

Calculation of local sensitivities

Local sensitivities (\(S_{ij}\)) play a key role on Optimum Experiment Design. These functions quantify that variation of the output variable (\(y_i\)) caused by a unitary variation of the model parameter (\(\theta_j\)), as described in the following equation. \[ S_{ij} = \frac{\partial y_i}{\partial \theta_j} \]

Typically, this derivative cannot be calculated in close form and one has to use numeric methods. bioOED uses the function sensFun() from the FME package for this purpouse, which approximates the local sensitivity functions using finite differences: \[ S_{ij}(t) \approx \frac{y \left( \theta + \Delta \theta \right) - y \left( \theta \right)} {\Delta \theta} \]

bioOED provides the function sensitivity_inactivation(), which allows to calculate estimate the local sensitivity functions. This function has 8 input arguments:

The inactivation model to use for the calculation is defined by the character argument inactivation_model. It must be consistent with the requirements of the function predict_inactivation from the bioinactivation package. Refer to the help package of this function for more information.

inactivation_model <- "Bigelow"

The environmental conditions are defined through the data frame temp_profile. Its columns, named time and temperature must provide discrete values of the temperature profile. Intermediate points will be calculated by linear interpolation.

temp_profile <- data.frame(time = c(0, 60),
                           temperature = c(30, 60))

The nominal values of the model parameters are provided by the arguments parms and parms_fix. Both arguments follow the same restriction as the equivalent ones for the predict_inactivation() function from the bioinactivation package. Note that local sensitivities will be calculated only for those parameters included in the parms argument.

parms_fix <- c(temp_ref = 57.5, N0 = 1e6)
parms <- c(D_R = 3.9,
           z = 4.2
           )

The local sensitivities are estimated at points uniformly distributed between the maximum and minimum values of time provided in the temp_profile argument. The number of points is provided by the integer argument n_times (100 by default).

The local sensitivities can be scaled using the varscale and parscale input arguments (see the help page of sensFun from the FME package). By default, both arguments are set to 1 (i.e. no scaling).

The argument sensvar allows choosing whether the sensitivity of the microbial count (\(N\)) or its decimal logarithm (\(\log N\)) with respect to the model parameters will be calculated. By default, the sensitivity of \(\log N\) is calculated.

sensitivity <- sensitivity_inactivation(inactivation_model,
                                        parms, temp_profile,
                                        parms_fix)
## Warning: package 'bindrcpp' was built under R version 3.4.4

The sensitivity_inactivation() function provides a data frame of class "sensFun" (defined in the FME package). This data frame contains the values estimated of the local sensitivities at each time point with respect to each one of the model variables.

head(sensitivity)
##           x  var          D_R             z
## 1 0.0000100 logN 0.000000e+00  0.000000e+00
## 2 0.6060606 logN 0.000000e+00 -1.691768e-07
## 3 1.2121212 logN 2.277381e-08 -3.595008e-07
## 4 1.8181818 logN 4.554761e-08 -5.921189e-07
## 5 2.4242424 logN 4.554761e-08 -8.881784e-07
## 6 3.0303030 logN 9.109522e-08 -1.205385e-06

Moreover, it includes an S3 method for plot() which allows the visualization of the evolution of the local sensitivities through the experiment.

plot(sensitivity)

Structural identificability

Some mathematical models may be ill-defined in the sense that different combination of its model parameters may provide the same output. This is defined as structural identificability of the model. It can be calculated as the correlation between the local sensitivity functions with respect to the different functions.

The bioOED package implements the function calculate_pars_correlation() which calculates the correlation between the sensitivity function of the model, providing an indicator of its structural identificability. This function has 6 input arguments:

They are defined identically to those defined for sensitivity_inactivation() and will not be repeated here.

parms_fix <- c(temp_ref = 57.5)
parms <- c(delta_ref = 3.9, z = 4.2, p = 1, N0 = 1e6)
temp_profile <- data.frame(time = c(0, 60), temperature = c(30, 60))
pars_correlation <- calculate_pars_correlation("Mafart", parms,
                                               temp_profile, parms_fix)

calculate_pars_correlation() provides a matrix of class parCorrelation containing the correlation between the sensitivity functions of the model parameters considered.

print(pars_correlation)
##             delta_ref           z           p          N0
## delta_ref  1.00000000  0.03549925 -0.98982135 -0.07864353
## z          0.03549925  1.00000000 -0.17736312  0.44558859
## p         -0.98982135 -0.17736312  1.00000000  0.01399121
## N0        -0.07864353  0.44558859  0.01399121  1.00000000
## attr(,"class")
## [1] "parCorrelation" "matrix"

This object includes an S3 method for plot, allowing to visualize the results of the calculation.

plot(pars_correlation)

Optimization of the reference temperature

The reference temperature can have a significant influence in the structural identificability of the model (Dolan and Mishra, 2013). The function optimize_refTemp() is able to calculate the reference temperature which optimizes the structural identificability of the model for a given set of nominal parameters and environmental conditions. This function has 8 input arguments:

The optimization is performed using the Brent method implemented in the optim() function from the stats package. The numeric argument temp_ref0 provides an initial guess for the optimization. Moreover, lower and upper bounds for the optimization are provided by the numeric arguments lower and upper.

temp_ref0 <- 57
lower <- 50
upper <- 70

The inactivation model to use for the calculation is defined through the argument inactivation_model. This argument has to be compatible with the equivalent one for the predict_inactivation() function from the bioinactivation package.

inactivation_model <- "Mafart"

The conditions of the experiment are defined through the temp_profile argument. It must be a data frame with columns named time and temperature providing discrete values of the temperature profile. Intermediate values are calculated by linear interpolation.

temp_profile <- data.frame(time = c(0, 60),
                           temperature = c(30, 60))

The nominal values of the model parameters are provided using the parms and parms_fix arguments. Both must be compatible with the requirements of the equivalent arguments for the predict_inactivation() function from the bioinactivation package. The parameters included in parms_fix will not be included in the calculation of the local sensitivities. Therefore, it is recommended to include in this argument those parameters which are known before performing the experiment. The D-value at the reference temperature (or \(\delta\)) and the z-value must be included in parms (i.e. cannot be fixed.)

Note that the D-value at the reference temperature (or \(\delta\)) depend on the reference temperature. Hence, the value provided in parms must be the one corresponding to the reference temperature provided in temp_ref0.

parms <- c(delta_ref = 3.9, z = 4.2, p = 1, N0 = 1e6)
parms_fix <- c()

The local sensitivities are estimated at a uniformly distributed time points. The number of time points used for the calculation is defined through the n_times argument. By default, it is set to 100.

optim_refTemp <- optimize_refTemp(temp_ref0, lower, upper,
                                  inactivation_model, parms, temp_profile,
                                  parms_fix)

optimize_refTemp() returns the object generated by the optim() function. The calculated optimum value of the reference temperature can be accessed in the par entry of this object.

print(optim_refTemp$par)
## [1] 57.77658

The convergence of the optimization algorithm should be checked in the convergence entry. If its value is different from 0, convergence was not achieved.

print(optim_refTemp$convergence)
## [1] 0

Optimum experiment design of microbial inactivation

The function inactivation_OED() is able to generate an optimum experiment design based on the FIM for a microbial inactivation experiment. This function has 10 input arguments:

The inactivation model to use for the calculation is defined by the input argument inactivation_model. It must be compatible with the predict_inactivation() function from the bioinactivation package.

inactivation_model <- "Mafart"

The nominal values of the model parameters are provided by the arguments pars and parms_fix. Note that local sensitivities will only be calculated for the model parameters included in pars. Hence, those included in parms_fix will not be considered for the OED. They also must be compatible with predict_inactivation().

parms_fix <- c(temp_ref = 57.5)
parms <- c(delta_ref = 3.9,
           z = 4.2,
           p = 1,
           N0 = 1e6
           )

The conditions of the experiment are defined through the input argument temp_profile. It must be a data frame which provides discrete values of time and temperature. Intermediate values will be calculated by linear interpolation.

temp_profile <- data.frame(time = c(0, 60), temperature = c(30, 60))

The number of measurements to be taken during the experiment are defined by the argument n_points. It must be an integer greater than 0.

n_points <- 5

The criteria for the OED is defined by the argument criteria. It must be a character with the value "D" or "E-mod". By default, the D-criterion is set.

The local sensitivity functions are calculated at a set of discrete uniformly distributed time points. The number of points to use is defined by the argument n_times (100 by default).

The variable to target for the OED is set by the argument sensvar. It must be a character equal to "logN" (default) or "N".

The OED can be generated using either a local or a global optimization algorithm. The selection is made using the argument optim_algorithm. In case it equals "global" (default) a global optimization algorithm is used. The local algorithm is used when this arguments equals "local".

The global optimization algorithm used is the MEIGO algorithm from the MEIGOR package (Egea et al., 2012). By default, a global solver with a maximum of 50000 function evaluations and printout on every step is selected. This can be changed through the opts_global argument. For information regarding the format of this argument, refer to the help page of the MEIGO() function.

opts_global <- list(maxeval=1000,  local_solver=0,
                    local_finish="DHC", local_iterprint=1)

global_OED <- inactivation_OED(inactivation_model, parms, temp_profile,
                               parms_fix, n_points, criteria = "D",
                               opts_global = opts_global)
## [1] "eSS R2014A - Enhanced Scatter Search"
## 
##  ------------------------------------------------------------------------------ 
##  essR - Enhanced Scatter Search in R 
## <c> IIM-CSIC, Vigo, Spain -  email: gingproc@iim.csic.es 
## ------------------------------------------------------------------------------ 
## 
## Number of diverse solutions automatically calculated: 50 
## Initial Pop: NFunEvals: 55 Bestf: -4.139976e-20 CPUTime: 0.16 Mean: -9.676263e-21 
## Iteration: 1 NFunEvals: 124 Bestf: -3.20892e-19 CPUTime: 0.38 Mean: -1.058837e-19 
## Iteration: 2 NFunEvals: 196 Bestf: -4.06661e-19 CPUTime: 0.54 Mean: -2.127713e-19 
## Iteration: 3 NFunEvals: 260 Bestf: -4.06661e-19 CPUTime: 0.69 Mean: -2.73448e-19 
## Iteration: 4 NFunEvals: 324 Bestf: -4.258508e-19 CPUTime: 0.83 Mean: -3.708941e-19 
## Iteration: 5 NFunEvals: 390 Bestf: -4.86536e-19 CPUTime: 1.13 Mean: -4.066797e-19 
## Iteration: 6 NFunEvals: 457 Bestf: -4.957436e-19 CPUTime: 1.3 Mean: -4.317302e-19 
## Iteration: 7 NFunEvals: 530 Bestf: -4.957436e-19 CPUTime: 1.47 Mean: -4.542365e-19 
## Iteration: 8 NFunEvals: 611 Bestf: -4.990607e-19 CPUTime: 1.67 Mean: -4.624459e-19 
## Iteration: 9 NFunEvals: 683 Bestf: -5.022175e-19 CPUTime: 1.83 Mean: -4.800706e-19 
## Iteration: 10 NFunEvals: 745 Bestf: -5.030958e-19 CPUTime: 1.99 Mean: -4.838313e-19 
## Iteration: 11 NFunEvals: 814 Bestf: -5.030958e-19 CPUTime: 2.14 Mean: -4.886249e-19 
## Iteration: 12 NFunEvals: 885 Bestf: -5.042043e-19 CPUTime: 2.31 Mean: -4.909638e-19 
## Iteration: 13 NFunEvals: 955 Bestf: -5.060717e-19 CPUTime: 2.47 Mean: -4.917669e-19 
## Iteration: 14 NFunEvals: 1033 Bestf: -5.062798e-19 CPUTime: 2.66 Mean: -4.926516e-19 
## 
## Final local refinement:  DHC 
## Initial point function value:  -5.062798e-19 
## NEvals: 10 Bestf: -5.063023e-19 
## NEvals: 13 Bestf: -5.064542e-19 
## NEvals: 25 Bestf: -5.066122e-19 
## NEvals: 32 Bestf: -5.071452e-19 
## NEvals: 38 Bestf: -5.087204e-19 
## NEvals: 51 Bestf: -5.089032e-19 
## NEvals: 53 Bestf: -5.089628e-19 
## NEvals: 56 Bestf: -5.089704e-19 
## NEvals: 65 Bestf: -5.090694e-19 
## NEvals: 70 Bestf: -5.090713e-19 
## NEvals: 76 Bestf: -5.091456e-19 
## NEvals: 89 Bestf: -5.092283e-19 
## NEvals: 91 Bestf: -5.092414e-19 
## NEvals: 92 Bestf: -5.092677e-19 
## NEvals: 94 Bestf: -5.093376e-19 
## NEvals: 105 Bestf: -5.093441e-19 
## NEvals: 108 Bestf: -5.093616e-19 
## NEvals: 115 Bestf: -5.093617e-19 
## NEvals: 121 Bestf: -5.093683e-19 
## NEvals: 126 Bestf: -5.093765e-19 
## NEvals: 140 Bestf: -5.093765e-19 
## NEvals: 142 Bestf: -5.093766e-19 
## NEvals: 148 Bestf: -5.09377e-19 
## NEvals: 157 Bestf: -5.093773e-19 
## NEvals: 159 Bestf: -5.093774e-19 
## NEvals: 165 Bestf: -5.093778e-19 
## NEvals: 171 Bestf: -5.093782e-19 
## NEvals: 182 Bestf: -5.093787e-19 
## NEvals: 184 Bestf: -5.093789e-19 
## NEvals: 187 Bestf: -5.093789e-19 
## NEvals: 202 Bestf: -5.09379e-19 
## NEvals: 210 Bestf: -5.09379e-19 
## NEvals: 222 Bestf: -5.09379e-19 
## NEvals: 224 Bestf: -5.09379e-19 
## NEvals: 225 Bestf: -5.09379e-19 
## NEvals: 226 Bestf: -5.09379e-19 
## NEvals: 237 Bestf: -5.09379e-19 
## NEvals: 248 Bestf: -5.09379e-19 
## NEvals: 251 Bestf: -5.09379e-19 
## NEvals: 258 Bestf: -5.09379e-19 
## NEvals: 263 Bestf: -5.09379e-19 
## NEvals: 272 Bestf: -5.09379e-19 
## NEvals: 275 Bestf: -5.09379e-19 
## NEvals: 278 Bestf: -5.09379e-19 
## NEvals: 292 Bestf: -5.09379e-19 
## NEvals: 296 Bestf: -5.09379e-19 
## NEvals: 307 Bestf: -5.09379e-19 
## Local solution function value: -5.09379e-19 
## Number of function evaluations in the local search: 321 
## CPU Time of the local search: 0.89 seconds 
## 
## Maximum number of function evaluations achieved 
## Best solution value -5.09379e-19 
## Decision vector 60 56.9697 50.90909 50.90909 3.030303 
## CPU time 3.55 
## Number of function evaluations 1034

The function inactivation_OED returns a list of class OEDinactivation. It contains the following entries:

print(global_OED$optim_times)
## [1]  3.030303 50.909091 50.909091 56.969697 60.000000

The object OEDinactivation includes an S3 implementation which allows the visualization of the results.

plot(global_OED)

Optimum experiment design including penalty function

The function inactivation_OED can generate experiment design unrealizable, with measurements too close in time. This issue is circumvented with the function inactivation_OED_penalty , which implements a penalty function to penalize unfeasible solutions. This function has 11 input arguments, 10 of which are identical to those defined for * inactivation_OED()*:

The only argument added by this function is time_min. It is a numeric value defining the minimum feasible time between measurements.

time_min <- 8
global_OED_penalty <- inactivation_OED_penalty(inactivation_model, parms,
                                               temp_profile, parms_fix,
                                               n_points, time_min,
                                               criteria = "D", 
                                               opts_global = opts_global)
## [1] "eSS R2014A - Enhanced Scatter Search"
## 
##  ------------------------------------------------------------------------------ 
##  essR - Enhanced Scatter Search in R 
## <c> IIM-CSIC, Vigo, Spain -  email: gingproc@iim.csic.es 
## ------------------------------------------------------------------------------ 
## 
## Number of diverse solutions automatically calculated: 50 
## Initial Pop: NFunEvals: 55 Bestf: -1.3272e-25 CPUTime: 0.24 Mean: 1.354399e+15 
## Iteration: 1 NFunEvals: 125 Bestf: -5.129777e-21 CPUTime: 0.42 Mean: 6.55447e+14 
## Iteration: 2 NFunEvals: 184 Bestf: -5.129777e-21 CPUTime: 0.58 Mean: 6.179311e+14 
## Iteration: 3 NFunEvals: 245 Bestf: -5.129777e-21 CPUTime: 0.75 Mean: 2.797339e+14 
## Iteration: 4 NFunEvals: 308 Bestf: -8.06316e-21 CPUTime: 0.92 Mean: 1.29701e+14 
## Iteration: 5 NFunEvals: 366 Bestf: -8.06316e-21 CPUTime: 1.08 Mean: 1.29701e+14 
## Iteration: 6 NFunEvals: 423 Bestf: -1.215287e-20 CPUTime: 1.22 Mean: 1.29701e+14 
## Iteration: 7 NFunEvals: 485 Bestf: -1.215287e-20 CPUTime: 1.39 Mean: 1.29701e+14 
## Iteration: 8 NFunEvals: 547 Bestf: -1.215287e-20 CPUTime: 1.55 Mean: 1.263231e+14 
## Iteration: 9 NFunEvals: 605 Bestf: -1.469084e-20 CPUTime: 1.78 Mean: -5.037988e-21 
## Iteration: 10 NFunEvals: 661 Bestf: -1.469084e-20 CPUTime: 1.97 Mean: -5.037988e-21 
## Iteration: 11 NFunEvals: 718 Bestf: -1.469084e-20 CPUTime: 2.14 Mean: -5.210136e-21 
## Iteration: 12 NFunEvals: 775 Bestf: -1.469084e-20 CPUTime: 2.36 Mean: -5.260408e-21 
## Iteration: 13 NFunEvals: 833 Bestf: -1.469084e-20 CPUTime: 2.51 Mean: -5.576662e-21 
## Iteration: 14 NFunEvals: 891 Bestf: -1.469084e-20 CPUTime: 2.83 Mean: -6.384699e-21 
## Iteration: 15 NFunEvals: 951 Bestf: -1.593267e-20 CPUTime: 3.14 Mean: -7.048029e-21 
## Iteration: 16 NFunEvals: 1010 Bestf: -1.668121e-20 CPUTime: 3.31 Mean: -7.296826e-21 
## 
## Final local refinement:  DHC 
## Initial point function value:  -1.668121e-20 
## NEvals: 7 Bestf: -1.694677e-20 
## NEvals: 9 Bestf: -1.746196e-20 
## NEvals: 10 Bestf: -1.774195e-20 
## NEvals: 15 Bestf: -1.866978e-20 
## NEvals: 31 Bestf: -2.20885e-20 
## NEvals: 38 Bestf: -2.583065e-20 
## NEvals: 48 Bestf: -2.585651e-20 
## NEvals: 50 Bestf: -2.589773e-20 
## NEvals: 61 Bestf: -2.590447e-20 
## NEvals: 63 Bestf: -2.591076e-20 
## NEvals: 64 Bestf: -2.592105e-20 
## NEvals: 70 Bestf: -2.763439e-20 
## NEvals: 82 Bestf: -2.783898e-20 
## NEvals: 89 Bestf: -2.783929e-20 
## NEvals: 96 Bestf: -2.789056e-20 
## NEvals: 98 Bestf: -2.799323e-20 
## NEvals: 101 Bestf: -2.908723e-20 
## NEvals: 111 Bestf: -2.914176e-20 
## NEvals: 113 Bestf: -2.925098e-20 
## NEvals: 124 Bestf: -2.926465e-20 
## NEvals: 138 Bestf: -2.949977e-20 
## NEvals: 140 Bestf: -2.961426e-20 
## NEvals: 141 Bestf: -2.984391e-20 
## NEvals: 150 Bestf: -2.985528e-20 
## NEvals: 152 Bestf: -2.987803e-20 
## NEvals: 153 Bestf: -2.992356e-20 
## NEvals: 166 Bestf: -2.9938e-20 
## NEvals: 171 Bestf: -2.9938e-20 
## NEvals: 179 Bestf: -2.995447e-20 
## NEvals: 181 Bestf: -2.99585e-20 
## NEvals: 187 Bestf: -2.995851e-20 
## NEvals: 194 Bestf: -2.995952e-20 
## NEvals: 206 Bestf: -2.996081e-20 
## NEvals: 208 Bestf: -2.996138e-20 
## NEvals: 209 Bestf: -2.99625e-20 
## NEvals: 212 Bestf: -2.997393e-20 
## NEvals: 222 Bestf: -2.99745e-20 
## NEvals: 224 Bestf: -2.997562e-20 
## NEvals: 230 Bestf: -2.997562e-20 
## NEvals: 237 Bestf: -2.997591e-20 
## NEvals: 239 Bestf: -2.997647e-20 
## NEvals: 245 Bestf: -2.997647e-20 
## NEvals: 252 Bestf: -2.997661e-20 
## NEvals: 258 Bestf: -2.997661e-20 
## NEvals: 266 Bestf: -2.997685e-20 
## NEvals: 272 Bestf: -2.997685e-20 
## NEvals: 276 Bestf: -2.997686e-20 
## NEvals: 284 Bestf: -2.997742e-20 
## NEvals: 288 Bestf: -2.997742e-20 
## NEvals: 296 Bestf: -2.997753e-20 
## NEvals: 304 Bestf: -2.997753e-20 
## NEvals: 306 Bestf: -2.997754e-20 
## NEvals: 313 Bestf: -2.997782e-20 
## NEvals: 318 Bestf: -2.997782e-20 
## NEvals: 321 Bestf: -2.997782e-20 
## NEvals: 333 Bestf: -2.997782e-20 
## NEvals: 335 Bestf: -2.997782e-20 
## NEvals: 340 Bestf: -2.997785e-20 
## NEvals: 346 Bestf: -2.997785e-20 
## NEvals: 354 Bestf: -2.99779e-20 
## NEvals: 361 Bestf: -2.99779e-20 
## NEvals: 364 Bestf: -2.99779e-20 
## NEvals: 372 Bestf: -2.997795e-20 
## NEvals: 376 Bestf: -2.9978e-20 
## NEvals: 384 Bestf: -2.997802e-20 
## NEvals: 390 Bestf: -2.997802e-20 
## NEvals: 393 Bestf: -2.997802e-20 
## NEvals: 403 Bestf: -2.997805e-20 
## NEvals: 407 Bestf: -2.997805e-20 
## NEvals: 415 Bestf: -2.997805e-20 
## NEvals: 423 Bestf: -2.997805e-20 
## NEvals: 425 Bestf: -2.997805e-20 
## NEvals: 426 Bestf: -2.997805e-20 
## NEvals: 443 Bestf: -2.997805e-20 
## NEvals: 445 Bestf: -2.997806e-20 
## NEvals: 446 Bestf: -2.997806e-20 
## NEvals: 449 Bestf: -2.997806e-20 
## NEvals: 457 Bestf: -2.997806e-20 
## NEvals: 459 Bestf: -2.997806e-20 
## NEvals: 460 Bestf: -2.997806e-20 
## NEvals: 467 Bestf: -2.997806e-20 
## NEvals: 469 Bestf: -2.997806e-20 
## NEvals: 478 Bestf: -2.997806e-20 
## NEvals: 482 Bestf: -2.997806e-20 
## NEvals: 494 Bestf: -2.997806e-20 
## Local solution function value: -2.997806e-20 
## Number of function evaluations in the local search: 500 
## CPU Time of the local search: 1.67 seconds 
## 
## Maximum number of function evaluations achieved 
## Best solution value -2.997806e-20 
## Decision vector 44 52 10.99965 2.999651 60 
## CPU time 4.98 
## Number of function evaluations 1011

The object returned by inactivation_OED_penalty() is identical to the one returned by inactivation_OED(). The newly calculated optimum design can be checked as before:

print(global_OED_penalty$optim_times)
## [1]  2.999651 10.999651 44.000000 52.000000 60.000000

Again, the results can be easily plotted.

plot(global_OED_penalty)

References

Alberto Garre, Pablo S. Fernandez and Jose A. Egea (2016). bioinactivation: Simulation of Dynamic Microbial Inactivation. R package version 1.1.2. https://CRAN.R-project.org/package=bioinactivation

Karline Soetaert and Thomas Petzoldt (2010). Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME. Journal of Statistical Software, 33(3), 1-28. URL http://www.jstatsoft.org/v33/i03/.

K. Grijspeerdt and P. Vanrolleghem (1999). Estimating the parameters of the Baranyi model for bacterial growth. Food Microbiology, 16(6), 593-605.

Kirk D. Dolan and Dharmendra K. Mishra (2013). Parameter Estimation in Food science. Annual revie of Food Science and Tehnology, 4, 401-422.

Jose A. Egea, David Enriques, Alexandre Fdez. Villaverde and Thomas Cokelaer (2012). MEIGOR: MEIGO - Metaheuristics for bioinformatics global optimization. R package version 1.0.0.