Three types of randomness were used to describe variability and uncertainty under the proposed bacterial inactivation model. Under the first type, which describes the variability in the initial bacterial count. The second type describes the uncertainty in an estimated parameter derived from the secondary models obtained from the bootstrap. The third type describes true randomness variability in bacterial reduction using Eq. 10. Using these approaches, the specific uncertainty in bacterial behavior was simulated as follows:

The above series of procedures was repeated 100 times to determine the variance, including variability and one type of uncertainty, in bacterial behavior. The R statistics (Ver. 3.5.1 for Mac OS X) statistical software was used to carry out all statistical analyses.

The spore survival kinetics under heating for each condition and the Weibull fitted survival curve based on 1,000 bootstrap replicates are shown together in Figure 4, in which the error bars indicate the standard deviations (from three observations) under the respective conditions. Based on the variation in the estimated Weibullian model, each condition has 1,000 replicates. The estimated Weibull model parameters (i.e., log₁₀δ andp) derived from the bootstrap were convergent, since all the R^-values of the parameters were equal to 1.0; they followed empirical distributions (Figure 5), with both decreasing as either the temperature increased, or the pH decreased. These distributions come from the variations in combinations of resampled observed values by the bootstrap. Therefore, it reflects the fitting uncertainty in Weibullian parameters and the thermotolerance variability of individual cells that arise from observations. The variation in p values was larger than that observed in log₁₀δ. Under some conditions in which the variation in p was particularly large, its distribution had two peaks (Figure 5). These indicate that it is difficult to detect the true one parameter from observed bacterial survival kinetics during thermal inactivation.

The secondary model was used to describe the relationship between the bootstrapped parameters and those for pH and temperature (Figure 6). In the figure, the points, error bars, and lines indicate the average values of the bootstrapped parameters, the standard deviations of the bootstrapped parameters, and the 1,000 replicates of the bootstrapped secondary model, respectively. It is seen that the secondary model has confirmed the results of previous studies (Mafart et al., 2002), in which the temperature dependencies of the δ parameter were shown to follow a conventional D-value form. The pH dependencies of δ also follow a conventional D-value form. Although it has been reported that p parameter does not show any temperature and/or pH dependency in most cases (Van Boekel, 2002), the p parameter values obtained in the present study show a temperature and pH dependency (Figures 6C,D). Similar to the present study, there are some reports that it was found that p decreased linearly as the temperature rose before reaching a constant level, suggesting that the results for p corresponded to a partially linear regression (Campanella and Peleg, 2001; Gomez et al., 2005).

Apart from temperature/pH dependency of the parameters (δand p), there is an interaction between the δ and p parameters, because the shapes depend on the both of parameters (Peleg, 2006). The developed 2DMC model enables to take into account the interaction of both the parameters, because it estimates the Weibull parameters using sets of secondary models of the corresponding δ and p parameters to resampled datasets. Comparing the fitted parameters’ interaction and the estimated parameters’ interaction by the bootstrapped secondary model sets, the bootstrapped secondary model sets explain linearly the interaction between δ and p parameter as shown in the right columns of Figure 5.

Comparing the variations in p-parameters and log₁₀δ derived from the bootstrapped parameters, the variations in the estimated p-parameter derived from the secondary model using bootstrapping were larger than those of log₁₀δ, possibly because they arose from small differences such as those between repetitions (Figure 6). These results confirm the difficulty to develop a secondary model of p from only one parameter per condition.

Simulations based on 1,000 repetitions of the secondary models derived from the non-parametric bootstrap methods were used to model the variation in bacterial reduction behavior, including the variability and uncertainty in the fitting model estimation (Figure 7). In the figure, the temperature profile and dynamical kinetic model predictions based on prior research (Peleg and Penchina, 2000; Campanella and Peleg, 2001) are shown as red and blue lines, respectively, while the gray lines indicate the results of 100 runs of simulated spore counts under non-isothermal inactivation. It is seen that the reduction rate follows the heating temperature, increasing and decreasing as the latter increases and decreases, respectively. Moreover, the predicted survival cell counts can describe the tendency of the observed bacterial behavior (0.12 < RMSE < 0.55). However, there were differences between observed and estimated values in some conditions. As described below, the cause of these discrepancies would be considered to be another factor that cannot be expressed in the 2DMC model. At small bacteria counts, there is a large variation in inactivation times until an arbitrary bacterial count (Figure 7). The variability in the time needed to reduce the surviving cell count to a specific level is increased in the range below 2 log₁₀ CFU. The kinetic prediction passed through the center of the estimated range of survival spore count derived from the stochastic model, it is suggested that kinetic models tend to describe the average behavior of bacterial heterogeneity. The points in the figure indicate the observed surviving spore counts and are used to validate the 2DMC simulation model. Although the observed values in the region below 2 log₁₀ CFU are within the 2DMC-estimated range, those above 2 log₁₀ CFU are outside of the estimated range.

The 2DMC developed in this study gives the distribution derived from the bootstrap and reflects the stochastic process describing bacterial reduction. In this manner, it describes both the uncertainty in observed value and the variability in bacterial lethality. Although most previous stochastic models describing bacterial behavior were based on only variability (Corradini et al., 2010; Aspridou and Koutsoumanis, 2015; Abe et al., 2019; Koyama et al., 2019b) or uncertainty (Schaffner, 1994), some researchers have advocated using models that integrate variability and uncertainty estimation (Cassin et al., 1998; Nauta, 2000; Pouillot and Delignette-Muller, 2010; Abe et al., 2019). In this study, the variations of parameters were described with bootstrap while the randomness of dead cell counts were described with MCMC. The variations of parameters can describe the fitting uncertainty and the heterogeneity in individual cell thermotolerance; the randomness of dead cell counts can describe the true randomness in behaviors of bacteria population with exactly the same properties.

To further assess the 2DMC approach, we compared it to “only-bootstrap” and “only-MCMC” models, i.e., to a dynamic kinetic model derived from 1,000 secondary model sets using the bootstrap method and a dynamic stochastic model derived from a secondary model without the application of bootstrapping, respectively. Figure 7 shows the 99% confidence intervals for the time required to obtain a specific decrease in the number of surviving cells under a thermal condition (pH 5.4, heating temperature 85°C) derived from the only-bootstrap (blue region), only-MCMC (red region), and 2DMC models (gray region). We compared the three models in terms of their confidence intervals. At surviving cell counts larger than 2 log₁₀ CFU, the confidence intervals of the 2DMC model matched those of the only-bootstrap model, whereas those of the only-MCMC model were much smaller. At surviving cell counts of less than 2 log₁₀ CFU, the confidence intervals of the 2DMC were larger than those of the other models. Furthermore, the confidence intervals of the 2DMC model corresponded to the variations in the observed values (Figure 8). These results further validate the use of a 2DMC model that describes both uncertainty and variability and suggest that a stochastic model combining multiple Monte Carlo simulations, such as the one used in this study, can accurately estimate cell reduction behavior that includes both variability and uncertainty.

The 2DMC simulation model proposed in the present study (Figure 7) is based on the previously developed Weibullian dynamic models. One of the fundamental formulations used to develop the 2DMC, Eq. 6, is based on the previous kinetic Weibullian dynamic studies (Campanella and Peleg, 2001; Peleg and Normand, 2004; Corradini and Peleg, 2009). The applicability of the Weibullian dynamic kinetic model to other types of bacteria has been demonstrated in the previous studies, with the results of these studies confirming our description of bacterial reduction behavior under non-isothermal heating by an adaptation of conventional kinetic models. In addition to the conventional Weibullian dynamic models mentioned above, the proposed model enables to assess the variation in bacterial count during an actual thermal process. This suggests the possibility of developing a 2DMC model based on a database such as ComBase (Baranyi and Tamplin, 2004) for not only average behaviors but also variabilities and uncertainties of vegetative cells’ behaviors as well as spores, although the results would need to be validated. Generalizing the proposed model, the model can also be used to calculate the bacterial or sporular death behaviors including variation in real-time under existing or planned industrial thermal processes.

There is a possibility of application of the Markov property for other dynamic model than thermal inactivation. The proposed 2DMC model performs MC simulation by applying the Markov property, i.e. by using only the current state to predict the next state, with any other information about the past states disregarded (Durrett, 2011). This simulation approach, referred to as MCMC, is generally used for the simulation of stochastic processes with probability densities known up to a constant of proportionality (Geyer, 1992). In this paper, the Markov property is introduced in Eq. 10. Unlike the conventional stochastic Weibullian dynamic model (Corradini et al., 2010), which does not consider the influence of the randomness of instantaneous survival ratio, the proposed stochastic dynamic model takes both S_(t) and the instantaneous count of surviving bacteria N_(t) into account. From Eq. 5, as adapted from previous work (Mattick et al., 2001; Peleg et al., 2005), the randomness of the instantaneous survival ratio depends on the instantaneous values of S_(t). Considering the dependency to instantaneous S_(t), it is important to take an idea of the Markov property into account for not only in Weibull models but in any type of reduction model. In addition to reduction behavior, formulations based on the Markov property such as Eq. 10 can describe growth behavior because growth behaviors depend on the instantaneous bacterial counts. Thus, by applying MCMC it is possible to create multiple stochastic prediction models even under dynamic environments.

The proposed dynamic model, however, might have some limitations in its ability to predict variations in the reduction behavior of arbitrary bacterial types. The proposed model ignores several factors, including the interactive influence between bacteria and thermal unevenness. As cell (and in this particular case spore) interactions are not taken into account in the developed modeling approach, it is possible that the proposed model cannot be used to predict bacterial population dynamics involving such interactions. As Eq. 10 describes change over an infinitesimal interval under any type of bacterial reduction model, it can be applied via Eq. 10 to express average bacterial behavior through a kinetic model. In this study, very small (100 μL) bacterial suspensions were tested to avoid uneven temperature conditions across the samples. Under conditions of actual thermal inactivation in food or beverages, however, there will be significant temperature unevenness as a result of heating, which should therefore be taken into account through the use of a temperature-predictive physical model established via the modeling of heat transfer with a fine mesh. By doing so, the problem of kinetic modeling of the interactive dependence among bacteria can be solved.

Another critical limitation with the proposed model is that it cannot completely express uncertainty and variability. The proposed 2DMC method attempts to use the bootstrap method to express uncertainty in the Weibull model parameter estimations, although it does not account for uncertainty related to experimental procedures (e.g., plate counting), model selection and other sources. The same applies for variability sources not explicitly taken into account in the present modeling approach such as intra-species (i.e., strain-to-strain) variability, other than true randomness and thermotolerance variability within a bacterial strain. The differences seen in Figure 7 between the estimated and observed surviving cell counts in the range above 2 log₁₀ CFU demonstrate its failure to obtain accurate parameter values. This discrepancy arises from various factors in the creation of the primary and secondary models, and the development of a more robust dynamic model requires the selection of a more appropriate combination of these models. The proposed model is based on a conventional kinetic dynamic model, and the use of inappropriate secondary or primary models can therefore result in inappropriate dynamic predictions. In particular, the bootstrap process through which the model estimates the parameter distribution involves the selection of a most suitable value for each iteration via the least square method. However, the standard error of fitting using the least squares method was not considered in this study. In future applications of the proposed method, this could be addressed by resampling the prediction parameters based on a normal distribution with a standard error equivalent to that of the fitting to widen the parameter distribution and obtain a more accurate parameter value. In addition, Bayesian analysis might be useful in describing the uncertainties in the primary parameters. Given that the observed values used for validation will themselves have inherent and unavoidable uncertainties, full validation of the stochastic process model will require a method that can completely account for uncertainty in observed values.