Bioinactivation FE: A free web application for modelling isothermal and dynamic microbial inactivation

Food Research International 112 (2018) 353–360

Contents lists available at ScienceDirect

Food Research International journal homepage: www.elsevier.com/locate/foodres

Bioinactivation FE: A free web application for modelling isothermal and dynamic microbial inactivation Alberto Garrea, Marta Clemente-Carazoa, Pablo S. Fernándeza, Roland Lindqvistb, Jose A. Egeac,

T ⁎

a

Departamento de Ingeniería de Alimentos y del Equipamiento Agrícola, Instituto de Biotecnología Vegetal, Universidad Politécnica de Cartagena (ETSIA), Paseo Alfonso XIII, 48, 30203 Cartagena, Spain b National Food Agency, Box 622, SE 751 26 Uppsala, Sweden c Departamento de Matemática Aplicada y Estadística, Universidad Politécnica de Cartagena, Antiguo Hospital de Marina (ETSII), Av. Dr. Fleming S/N, 30202 Cartagena, Spain

A R T I C LE I N FO

A B S T R A C T

Keywords: Mathematical modelling Bioinactivation Free software Microbial inactivation Food safety Web application

Mathematical models developed in predictive microbiology are nowadays an essential tool for food scientists and researchers. However, advanced knowledge of scientific programming and mathematical modelling are often required in order to use them, especially in cases of modelling of dynamic and/or non-linear processes. This may be an obstacle for food scientists without such skills. Scientific software can help making these tools more accessible for scientists lacking a deep mathematical or computing background. Recently, the R package bioinactivation was published, including functions (model fitting and predictions) for modelling microbial inactivation under isothermal or dynamic conditions. It was uploaded to the Comprehensive R Archive Network (CRAN), but users need basic R programming knowledge in order to use it. Therefore, it was accompanied by Bioinactivation SE, a user-friendly web application including selected functions in the software for users without a programming background. In this work, a new web application, Bioinactivation FE, is presented. It is an extension of Bioinactivation SE which includes an interface to every function in the bioinactivation package: model fitting of isothermal and nonisothermal experiments, and generation of survivor curves and prediction intervals. Moreover, it includes several improvements in the user interface based on the users' feedback. The capabilities of the software are demonstrated through two case studies using data published in the scientific literature. In the first case study, the response of Escherichia coli to isothermal and non-isothermal treatments is compared, illustrating the presence of an induced thermal resistance. In the second, the effect of nanoemulsified D-limonene on the thermal resistance of Salmonella Senftenberg is quantified.

1. Introduction Several efforts have been dedicated to developing approaches for ensuring the safety and quality of food products. Various methods have been applied, advancing from qualitative towards quantitative approaches and risk assessments (Coleman & Marks, 1999; Lindqvist, Sylvén, & Vågsholm, 2002). These methods require an accurate description of microbial responses in the whole food chain, i.e. during harvesting, processing, distribution, storage, preparation and consumption. In this respect, predictive microbiology has proved its applicability to achieve this purpose. Predictive microbiology is based on the hypothesis that past experience can be used to predict the microbial response (Ross, 1996). Therefore, data generated under laboratory conditions can be used to

⁎

estimate the microbial concentration during the farm-to-fork life of the product using mathematical models (McMeekin, Olley, Ratkowsky, & Ross, 2002). The first scientific papers on predictive microbiology date from the early 1920s, when the inactivation kinetics of Clostridium botulinum spores during isothermal treatments was studied (Bigelow, 1921; Bigelow & Esty, 1920; Esty & Meyer, 1922). However, its applicability was quite limited until the 1980s, most likely due to the limited computational power available at the time. Nowadays, predictive microbiology models are basic tools for food safety studies and quantitative risk assessment (McMeekin, Mellefont, & Ross, 2007). However, the application of mathematical models is still, to some extent, limited. This is likely due to the absence of stand-alone software solutions, and the lack of free open databases and software (PlazaRodríguez et al., 2015). Furthermore, many food microbiologist lack

Corresponding author. E-mail addresses: [email protected] (P.S. Fernández), [email protected] (J.A. Egea).

https://doi.org/10.1016/j.foodres.2018.06.057 Received 22 May 2018; Received in revised form 21 June 2018; Accepted 24 June 2018 Available online 25 June 2018 0963-9969/ © 2018 Elsevier Ltd. All rights reserved.

Food Research International 112 (2018) 353–360

A. Garre et al.

Fig. 1. Overview of the layout of Bioinactivation FE as it is loaded by the web browser.

Thus, the R package bioinactivation includes all the functions usually required for mathematical modelling of microbial inactivation (model fitting and generation of predictions for isothermal and dynamic conditions), but there is no free user-friendly application including all these functions. Consequently, due to the positive feedback received from the users of bioinactivation, a novel web application, Bioinactivation FE, described in the present work has been developed. This application includes an interface to all functions in the bioinactivation package. In this article, the capabilities of Bioinactivation FE are showcased through two different case studies using data already published in the scientific literature. The first compares the inactivation patterns of Escherichia coli during isothermal and non-isothermal treatments, whereas the second assesses how different concentrations of nanoemulsified D-limonene affects the thermal inactivation of Salmonella Senftenberg. Besides, a general description of the functions implemented and the user interface is provided.

the necessary skills in statistics, mathematical modelling and scientific programming required to develop and apply mathematical models (Granato, de Araújo Calado, & Jarvis, 2014; Granato, Nunes, & Barba, 2017; Nunes, Alvarenga, Santos, and Granato, 2015). A way to circumvent these limitations is the development of user-friendly software applications, so that mathematical models can be used and their results interpreted by users without advanced education in mathematical modelling and software development (Muñoz-Cuevas, Metris, & Baranyi, 2012). A recent review article described the software tools presented at the International Congress of Predictive Microbiology in Foods in 2013 (ICPMF8) software fair (Tenenhaus-Aziza & Ellouze, 2015). However, none of the freely available software included all the tools required to model microbial inactivation: fitting of isothermal and non-isothermal experiments, and generation of predictions under dynamic and static conditions (Garre, Fernández, Lindqvist, & Egea, 2017). Since the review by Tenenhaus-Aziza & Ellouze, two new free applications have been developed: IPMP Global Fit and bioinactivation. IPMP Global Fit (Huang, 2017) is a tool for curve fitting of isothermal growth and inactivation experiments but it does not generate predictions. It is freely available and its files are downloaded and executed locally. Bioinactivation (Garre et al., 2017) includes functions for the model fitting of both isothermal and non-isothermal experiments, as well as for generating predictions. It is available as an R package (bioinactivation), and uploaded to CRAN making it publicly available as an open software (https://CRAN.R-project.org/package=bioinactivation). The use of the bioinactivation package requires basic knowledge of R programming from the user. For that reason, a user-friendly web interface including only the functions for model fitting of iso- and non-isothermal experiments was developed: bioinactivation SE (https://opada-upct.shinyapps. io/bioinactivation_SE/).

2. Materials and methods 2.1. Bioinactivation FE Bioinactivation FE is a user-friendly interface to all the functions included in the bioinactivation package for R (Garre et al., 2017); i.e. simulation of survivor curves, generation of prediction intervals, and model fitting of isothermal and dynamic experiments. It has been built using the shiny (Chang, Cheng, Allaire, Xie, & McPherson, 2017) and shinydashboard (Chang & Ribeiro, 2017) packages. The combination of both packages allows the creation of a web page with input and output widgets able to run R code in the background (R Core Team, 2016). Therefore, users can access advanced mathematical operations implemented in R without any skills in this programming language. 354

Food Research International 112 (2018) 353–360

A. Garre et al.

2.1.2. The isothermal fitting module The isothermal fitting module includes functions to fit microbial inactivation models using the microbial counts observed in a set of isothermal inactivation experiments performed at different temperature levels. The user must input a table with the observed microbial density as log-fraction of survivors (log10S = log10N/N0), the elapsed time and the temperature of the experiment. This can be done using a text file, or by entering the data manually into the application. The model fitting is performed using Nonlinear Least Squares (Bates & Watts, 2007) following a one-step procedure, which has proved superior to the classical two-step algorithm (Fernández, Ocio, Fernández, Rodrigo, & Martinez, 1999). The user is free to select any inactivation model among those implemented in the bioinactivation package, as well as fixing any of the model parameters. Nevertheless, due to the fitting algorithm used, the user must provide initial guesses for unknown model parameters. This can be done using prior-knowledge or by exploring the experimental observations. The application generates a plot of model predictions and experimental data for each temperature. Moreover, a table with estimated values and standard deviations, as well as confidence intervals at the 95% confidence level is reported. Confidence intervals are calculated using the approximation to normality (Box, Hunter, & Hunter, 2005). In order to evaluate the goodness of the fit, the ME, RMSE, log-likelihood, AIC, AICc, BIC, Af and Bf of the fitted curve are reported, as well as the results of a Shapiro-Wilk normality test and a residual plot.

Bioinactivation FE is free of charge and accessible online (https:// opada-upct.shinyapps.io/bioinactivationFull/); the only requirement to access it is an internet connection and a web browser. Alternatively, the application can be compiled, so an internet connection is no longer a requirement. Fig. 1 illustrates the application as it is loaded by the Web browser. The body of the application includes all the widgets with which the user can interact, as well as the tables and figures for displaying results. It is interactive, in the sense that its contents are dynamically generated as the user interacts with the application. The side bar allows the user to choose any of the three modules of the application: predictions, isothermal fitting and non-isothermal fitting. It also includes a link to the article published by Garre, Fernández, Lindqvist & Egea (2017) describing the mathematical functions implemented in the bioinactivation package, as well as an About page with miscellaneous information regarding the application. Every module contains example data, so that the user can readily experiment with the application as it is loaded. In the following subsections, the functions included in each module are described. A more thorough description, including information and screenshots for every single widget, is provided in the user manual of Bioinactivation FE (Supplementary material A1). 2.1.1. The prediction module The prediction module includes functions to simulate the microbial response during isothermal or non-isothermal processing treatments. For that, the user must select an inactivation model among those implemented in the bioinactivation package (Bigelow, Peleg, Mafart and Geeraerd models) and to input the values of the model parameters. Values of model parameters characterizing the microbial thermal resistance can be obtained by model fitting of experimental data (see Sections 2.1.2 and 2.1.3), or from the published literature and databases. Furthermore, “hypothetical” values of the model parameters can be used to test “what-if” cases. The user also has to input the temperature profile of the treatment, providing discrete values for the temperature at different time points of the treatment. This information can be input either by uploading a text file or by manually typing the values in the application. Temperature values at time points not supplied by the user are estimated by linear interpolation. Note that isothermal profiles can be simulated by using the same temperature for every time point. Once this input data has been provided, bioinactivation solves the differential equation of the inactivation model using the LSODA algorithm (Hindmarsh, 1983), through the deSolve R package (Soetaert, Petzoldt, & Setzer, 2010). Bioinactivation FE generates a plot of the predicted microbial response, which can be exported as a text file for further processing using external software. In addition, it allows the user to compare model predictions with experimental data. Microbial counts potentially measured in the lab can be entered manually into the tool or uploaded via text files. Bioinactivation is not restrictive with respect to the system of units used for the analysis. It only requires consistency (e.g. using the same units for the parameter N0 and the experimental microbial counts, or for the D-value and the elapsed time). The application, then, generates a plot comparing the model predictions and the experimental data. This plot can be edited, so that it can be made publication-ready. Furthermore, the tool calculates the residuals of the prediction, and reports several statistical indexes measuring the goodness of the prediction: Mean Error (ME), Root Mean Squared Error (RMSE), log-likelihood, Akaike Intormation Criterion (AIC) (Akaike, 1974), Corrected AIC (AICc) (Hurvich & Tsai, 1993), Bayesian Information Criterion (BIC), and Accuracy (Af) and Bias factors (Bf). The two last indexes are broadly used in predictive microbiology for the evaluation of model predictions. Although they were both originally proposed by Ross (1996), a later article proposed a refined version (Baranyi, Pin, & Ross, 1999), which is the one implemented in Bioinactivation FE. Furthermore, the results of a Shapiro-Wilk normality test of the residuals are reported (Shapiro & Wilk, 1965).

2.1.3. The non-isothermal fitting module The non-isothermal fitting module allows the characterization of microbial inactivation using observations obtained during processing under dynamic temperature conditions. This module includes all the functionalities of Bioinactivation SE (Garre et al., 2017), with several improvements in the user interface based on the users' feedback. It requires the user to enter the observed microbial density, as well as the temperature profile of the treatment. This information can be introduced either manually or using a text file. The temperature at time points not provided by the user is estimated by linear interpolation of the data provided. Alternatively, by setting every temperature of the thermal profile to the same value, this module can be used to analyse individual isothermal treatments. The user can select which of all models implemented in the bioinactivation package to use. Two different fitting algorithms are implemented in the software: the Levenberg-Marquardt algorithm (Moré, 1978), a form of non-linear regression, and the adaptive Markov-Chain Monte Carlo method (MCMC) proposed by Haario, Laine, Mira & Saksman (2006), through the implementation included in the FME package (Soetaert & Petzoldt, 2010). The Levenberg-Marquardt algorithm is less computationally intensive, but can converge to sub-optimal solutions for the non-linear models used to describe microbial inactivation. The MCMC method is a stochastic method which may be better at dealing with non-linear problems, although the number of iterations required for its convergence must be evaluated by the user. Bioinactivation FE makes it very easy to change between both methods. It is, thus, recommended to perform the model fitting using both algorithms (as well as different number of iterations for the MCMC method) and to compare the results. Both the Levenberg-Marquardt and MCMC algorithms require the user to define initial guesses for model parameters, as well as upper and lower bounds. In order to support the user in choosing initial values, Bioinactivation FE reports a plot comparing the experimental observations against the model predictions for the initial guess, and the lower and upper bounds. The user is free to select which model parameters to fit to the data and which to set to fixed values. The application reports the estimated values and standard deviations of the model parameters, as well as confidence intervals. Note however, that the procedure followed for the calculation of the latter depends on the fitting algorithm selected. For the Levenberg355

Food Research International 112 (2018) 353–360

A. Garre et al.

The Mafart (Mafart, Couvert, Gaillard, & Leguerinel, 2002) model considers that the deviations from log-linearity are due to heterogeneity in the responses of separate microbial cells to the thermal stress. It is assumed that the time that a single cell can withstand a constant thermal stress follows a Weibull distribution, with shape factor, p, and scale factor,δ(T). In this model, the shape factor is considered temperature independent, whereas the scale factor is assumed to have an exponential relationship with temperature, in a similar way as in the Bigelow model. The Peleg (Peleg & Cole, 1998) model also assumes that the time required to inactivate an single cell under constant stress follows a Weibull distribution. However, it uses a different parametrization of the distribution, using the parameters n and b(T). Moreover, it assumes a log-logistic relationship between the parameter b(T) and temperature. In this secondary model, the inactivation rate is zero for temperatures lower than the critical temperature, Tc. For temperatures much higher than the critical one, a linear relationship between both variables is assumed, with slope kb. This model describes a super-linear transition between both regimes when the temperature is close to the critical one.

Marquardt, confidence intervals are calculated using the linear approximation to normality, whereas for the MCMC algorithm they are estimated as the credible intervals of the posterior probability distribution of the model parameters (Garre et al., 2017). Furthermore, a plot comparing model predictions and experimental observations is provided by the tool. The plot can be edited to make it publicationready. The goodness of the fit can be evaluated based on the reported ME, RMSE, log-likelihood, AIC, AICc, BIC, Af and Bf, the residual plot and the results of a Shapiro-Wilk test of normality. Moreover, estimates of parameter correlations are reported. If the MCMC algorithm was selected for model fitting, a trace plot of the Markov Chain is reported, allowing for a check of its stationarity. Based on the model fitted, the user can calculate a prediction interval for the microbial counts during a non-isothermal (or isothermal) treatment. This is calculated by applying bootstrap using the functions implemented in the bioinactivation package (Garre et al., 2017). The application does not check the convergence of the prediction interval, so the user should repeat the calculations using a different number of iterations to ensure convergence. The application generates a plot of the prediction interval, which can be edited to make it publication-ready. Both the calculation of the prediction interval and the model fitting using the MCMC algorithms use stochastic procedures which need the generation of random numbers. The application includes a button to set the internal seed of the pseudo-random number generator, providing reproducibility of the results.

2.3. Microbiological data The functions included in Bioinactivation FE are showcased in this article using two different case studies. The first study deals with the characterization of the thermal resistance of Escherichia coli during isothermal and non-isothermal treatments. In the second case study, the effect of adding a D-limonene nanoemulsion to the heating medium on the thermal resistance of Salmonella Senftenberg is investigated. The data for the cases studies are from Garre et al. (2018) and RosChumillas, Garre, Maté, Palop & Periago (2017), respectively. Consequently, the methodology followed is only briefly described here. The culture preparation and the evaluation of the cell viability were performed according to standard microbiological methods. The experiments were performed using E. coli type strain (CECT 515) and Salmonella Senftenberg CECT 4565, both supplied by the Spanish Type Culture Collection (CECT). The thermal resistance of both microorganisms in liquid media under controlled temperature (isothermal or dynamic) conditions was evaluated using a Mastia thermoresistometer (Conesa, Andreu, Fernández, Esnoz, & Palop, 2009). The thermal resistance of E. coli was evaluated using a set of isothermal experiments at four different temperature levels (52.5, 55, 57.5 and 60 °C), as well as several non-isothermal profiles of different shapes. In this case study, the complete set of isothermal profiles was used, as well as one of the non-isothermal profiles. The thermal resistance of S. Senftenberg in the presence of different concentrations of a D-limonene nanoemulsion (0.1, 0.5 and 1 mM) was evaluated using a set of isothermal inactivation experiments performed at 47.5, 50, 52.5 and 55 °C. The nanoemulsion was prepared according to Ros-Chumillas, Garre, Maté, Palop & Periago (2017). Moreover, experiments without any D-limonene in the heating media were kept as control samples.

2.2. Microbial inactivation models All microbial inactivation models included in the bioinactivation package can be accessed using Bioinactivation FE. This includes the Bigelow, Peleg, Mafart and Geeraerd models. Table 1 shows the differential and algebraic equations defining these models, where N stands for the microbial density at time t, usually expressed as CFU/ml. The Bigelow model considers that, under constant temperature, the microbial density follows a log-linear relationship with the elapsed time of the experiment (Bigelow, 1921). The rate of the inactivation is quantified through the D-value (D(T)), which stands for the time required to apply a thermal stress at a given temperature in order to reduce the microbial load by 90%. This models considers that the D-value follows an exponential relationship with temperature (Bigelow, 1921). The sensitivity of the bacterial cells to temperature (T) changes is quantified using the z-value (z), which equals the temperature increase required to cause a ten-fold reduction of the D-value. The Geeraerd model (Geeraerd, Herremans, & Van Impe, 2000) also considers that the microbial inactivation under a constant stress is loglinear. However, it includes additional terms to account for shoulder and tail effects. The shoulder is modelled under the assumption that a theoretical substance Cc exists that inhibits microbial inactivation. The tailing effect is modelled using the Verhulst log-logistic term (Verhulst, 1838), so that microbial numbers cannot become lower than the tail height, Nres.

3. Results and discussion

Table 1 Summary of the microbial inactivation models included in Bioinactivation FE. Primary model Bigelow

d log10 N dt

Peleg

d log10 N dt

Mafart

d log10 N dt

Geeraerd

dN dt dCc dt

=−

=−

1 D (T )

log10 D (T ) = log10 Dref −

( = −p ( )t

= −b (T )·n· −

log10 S b (T )

)

n−1 n

p 1 p−1 δ (T )

1 k (T ) 1 + Cc max

3.1. Case study a: thermal stress adaptation of Escherichia coli

Secondary model

(

1−

)

N

z

b(T) = ln (1 + ekb(T−Tc))

log10 δ (T ) = log10 δref − Nres N

The ability of bacterial strains to develop a stress resistance after surviving a mild stress has been reported by several authors. (Hassani, Mañas, Condón, & Pagán, 2006; Hill, Cotter, Sleator, & Gahan, 2002; Valdramidis, Geeraerd, Bernaerts, & Van Impe, 2006). The application of a sub-lethal stress (e.g. a heat shock) to a bacterial population can trigger an acclimation response of the bacterial cells, making them more resistant to subsequent stresses (Hassani, Manas, Pagán, & Condón, 2007; Hill et al., 2002; Lin & Chou, 2004; Richter, Haslbeck, & Buchner, 2010). This effect has also been observed during dynamic

T − Tref

kmax (T ) =

T − Tref z

T − Tref kref 10 z

= −kmax (T )·Cc

356

Food Research International 112 (2018) 353–360

A. Garre et al.

This overparameterization is confirmed by the values of AICc reported in Table 2. The AICc is lower for the Bigelow model (64.7) than the Mafart model (67), indicating that the inclusion of an extra parameter does not significantly improve the description of the data. On the other hand, the Peleg model assumes a different secondary model than the Bigelow and Mafart models. Hence, this model is different even when the shape factor of the Weibull distribution equals one. Indeed, according to the statistical indexes reported in Table 2, the Peleg model is best at describing the data set, with the lowest RMSE (0.38) and AICc (46.27). Using the feature of Bioinactivation FE to generate plots comparing experimental data and predictions of the fitted model, it can be illustrated how the Peleg (Fig. 2B) model is slightly better than the Bigelow model ((Fig. 2A) at describing inactivation at 57.5 °C. The model parameters in Table 2 were used to predict microbial inactivation during one of the non-isothermal inactivation profiles tested by Garre et al. (2018) using Bioinactivation FE. The application generates a plot comparing the model predictions against the experimental observations. Nevertheless, in order to compare the three curves, the model predictions were exported using the functions included in Bioinactivation FE and a plot comparing the predictions of the three models and the experimental data (Fig. 3) was generated using external software. All three models predict a microbial inactivation much higher than observed in the dynamic experiments which shows a dynamic effect possibly indicating the relevance of stress acclimation in this experiment, as reported by Garre et al. (2018). Bioinactivation FE calculates a ME of −3.32, −3.19 and − 3.96 log CFU/ml for the Bigelow, Mafart and Peleg models, respectively. Despite illustrating the dynamic effect (potential stress acclimation) in Fig. 3, its quantification is not possible. For instance, it would be hard to compare the magnitude of stress acclimation in two different treatments using only such plots. In order to do so, one can fit inactivation models to the non-isothermal data and do inference on the estimated model parameters. In order to illustrate this methodology, the nonisothermal inactivation data reported by Garre et al. (2018) was fitted using the non-linear regression algorithm included in Bioinactivation FE. The estimated parameter values and standard deviations, as well as the indexes describing the goodness of the fit, are reported in Table 3. Bioinactivation FE estimates model parameters with very high standard deviation for the Bigelow model, making parameters practically unidentifiable (Balsa-Canto, Alonso, & Banga, 2010). This is likely due to the high correlation between the D and z-value of the Bigelow model (ρ = − 1). Therefore, the D and z-values cannot be estimated simultaneously using this dataset. A way to circumvent this issue is by fixing one of the correlated parameters to a known value and fitting the other (Vilas et al., 2018). Bioinactivation FE allows the setting of any model parameter to a selected value. Hence, the fitting was repeated setting the z-value to the value estimated from the isothermal experiments (5.18 °C). As a result, the D-value has been estimated with an acceptable accuracy (29.55 min (1.60)) without a significant loss in accuracy, as shown by the moderate increase in RMSE (4.35 log CFU/ml fitting all the parameters, 4.51 fixing the z-value) and the better AICc (−38.53 versus 35.75). Moreover, using the same z-value as the one estimated from the isothermal experiments allows an easier comparison between the obtained results. From the isothermal experiments a Dvalue at 52.5 °C of 12.10 min was estimated. Thus, the dynamic effect, possibly stress acclimation, results in more than a two-fold increase in the D-value. Nevertheless, this conclusion is based on the hypothesis that the Bigelow model is able to describe the microbial inactivation in cases where the thermal acclimation is relevant. It is advisable to use models specifically developed to account for the stress acclimation for this type of profile (Garre et al., 2018; Stasiewicz, Marks, Orta-Ramirez, & Smith, 2008; Valdramidis, Geeraerd, & Van Impe, 2007). A similar issue is encountered for the Peleg model; the parameter k is not practically identifiable. Again, a high linear correlation is observed between the model parameter k and Tcrit (ρ = − 0.90). Hence,

mild heat treatments, where bacterial cells may undergo a mild heating phase before lethal temperatures are reached. Due to the physiological acclimation developed during the heating phase, the microbial inactivation achieved in the treatment will be lower than expected. This implies a potential risk for the consumer. The relevance of stress acclimation in a thermal treatment can be evidenced by modelling the microbial inactivation achieved during isothermal and mild dynamic experiments. During isothermal experiments, lethal temperatures are applied instantaneously to the microbial cells. Hence, they do not have time to develop any kind of stress acclimation prior to the treatment. As a consequence, model parameters (e.g. D and z-values) obtained under isothermal conditions are those representing the most sensitive, non-acclimated, bacterial response (Garre et al., 2018). These model parameters can, then, be used to predict the microbial density during a non-isothermal profile. If stress acclimation is relevant, the microbial density observed experimentally will be significantly higher than the one predicted by the model. This difference will increase with the degree of acclimation. The web interface of Bioinactivation FE includes functions for model fitting of isothermal experiments, as well as for the generation of predictions. It can, thus, be used to evaluate the relevance of stress acclimation of bacterial populations in a dynamic treatment. The isothermal inactivation data for E. coli reported by Garre et al. (2018) were used in this example. In that study, only the Bigelow model was fitted to the isothermal data. However, with Bioinactivation FE it is easy to switch between inactivation models, facilitating the comparison between models. Therefore, the four models included in Bioinactivation FE were fitted to the experimental data. Table 2 reports the model parameters estimated by the tool using a non-linear regression algorithm to isothermal data, as well as the statistical indexes describing the quality of the model fit. As expected, the D52.5 (12.10 min (sd = 0.56 min)) and zvalues (5.18 °C (0.12)) estimated for the Bigelow models are identical to those estimated by Garre et al. (2018). Garre et al. (2018) concluded that the Bigelow model was able to accurately describe the isothermal inactivation data; i.e. that it was loglinear. This hypothesis is confirmed from the results of the model fitting of the two Weibullian models (Mafart and Peleg). In this family of models, a shape factor of the Weibull distribution equal to one indicates that the microbial inactivation is log-linear. In both cases, the shape factor is not significantly different from one (p = .99 (0.08), n = 1 (0.07)). When p = 1, the Mafart model is equivalent to the Bigelow model. Indeed, Bioinactivation FE estimates the same value and standard deviation for the z-value of both models. Moreover, the estimated δ52.5 (11.96 min (1.46)) is not significantly different from the estimated D52.5 of the Bigelow model (12.10 min (0.56)). However, the estimated standard deviation is slightly higher for δ52.5 thanD52.5. This is likely due to the overparameterization of the Mafart model, which has one additional model parameter that does not significantly improve the model fit (Vilas, Arias-Mendez, Garcia, Alonso, & Balsa-Canto, 2018).

Table 2 Estimated values (and standard deviations) of the model parameters of each model implemented in Bioinactivation FE for the isothermal inactivation of Escherichia coli, as well as values of the statistical indexes describing the goodness of the model fit.

Bigelow Mafart

Peleg

Parameter

Estimate

ME

RMSE

AICc

Af

Bf

D52.5 (min) z (°C) δ52.5(min) z (°C) p (−) k (1/°C) n (−) Tcrit (°C)

12.10 (0.56) 5.18 (0.12) 11.96 (1.46) 5.18 (0.12) 0.99 (0.08) 0.58 (0.03) 1 (0.07) 56.95 (0.27)

0.02

0.49

64.7

3.08

1.04

0.02

0.49

67.00

3.08

1.05

0.02

0.38

46.27

2.42

1.04

357

Food Research International 112 (2018) 353–360

A. Garre et al.

Fig. 2. Comparison between the experimental results of isothermal inactivation of E. coli and the curve fitted using the Bigelow (A) and the Peleg (B) models.

the model fitting has been repeated fixing the value of k to that estimated from the isothermal experiments (0.58 1/°C). This allowed the accurate estimation of the model parameters. A shape factor equal to one is, then, estimated from the non-isothermal experiments, as well as from the isothermal experiments. Therefore, according to the Peleg model, the stress acclimation only has an effect on the critical temperature. A critical temperature of 59.47 °C (0.96) has been estimated, 2.5 °C higher than that estimated from isothermal experiments, indicating the presence of a stress acclimation.

3.2. Case study B: study of the combined effect of thermal stress and Dlimonene in the inactivation of Salmonella Senftenberg In the field of food microbiology, it is very common to evaluate the effect of an additive (e.g. an antimicrobial) on the thermal resistance of a bacterium based on parameter inference. This methodology implies fitting the same inactivation model to experiments realized under similar conditions and to compare the estimated model parameters (e.g. the D-values). Ros-Chumillas et al. (2017) used this technique to evaluate the effect that the inclusion of a D-limonene nano-emulsion in the heating medium had on the inactivation of S. Senftenberg. The data used in that research was introduced into Bioinactivation FE, and the isothermal inactivation data obtained at each D-limonene concentration

Fig. 3. Comparison of the microbial counts of E. coli obtained under non-isothermal conditions against the predictions of the Bigelow (−), Mafart (−) and Peleg (−) models based on isothermal conditions (Table 2). 358

Food Research International 112 (2018) 353–360

A. Garre et al.

respect to the experimental observations. This is not the case for the Mafart and Peleg model, whose ME remains between 0.01 and − 0.01 log CFU/ml for every case studied. Furthermore, the AICc calculated for the Bigelow model is much larger than that calculated for the Peleg and Mafart models. As a consequence, the survivor curves observed in these experiments significantly deviate from log-linearity. Between the Mafart and Peleg model, it is unclear which model is better at describing the dataset. For the control samples and those with a 0.05 mM D-limonene concentration, the Peleg model has a lower AICc than the Mafart model, whereas at concentrations of 0.1 mM the Mafart model is better, and at 0.01 mM they are similar. It is, thus, unclear whether the Mafart or Peleg model is the best for this dataset. The similarities between both models are reasonable, because both are based on the hypothesis that the thermal resistance of individual microbial cells in the population follow a Weibull probability distribution. The model parameters reported in Table 4 can be used to quantify the effect that the addition of the D-limonene nanoemulsion has on the thermal resistance of S. Senftenberg. The δ-value represents the time required to apply the stress in order to cause the first ten-fold reduction of the microbial population. The values of δ50 obtained with a D-limonene concentration of 0.01 mM, 0.05 mM and 0.1 mM were, respectively, 8.3, 52.7 and 116 times larger than those obtained for control samples. Moreover, the parameter Tcrit of the Peleg model indicates a temperature below which no thermal inactivation occurs. The addition of a concentration 0.01 mM, 0.05 mM and 0.1 mM D-limonene results in a reduction of Tcrit of 5, 7 and 10 °C with respect to the control samples. Therefore, according to both models, the addition of this antimicrobial induces a significant reduction in the thermal resistance of S. Senftenberg.

Table 3 Estimated values (and standard deviations) of the model parameters of each model implemented in Bioinactivation FE for the inactivation of Escherichia coli under non-isothermal conditions, as well as values of the statistical indexes describing the goodness of the model fit.

Bigelow

Bigelow (fixed zvalue)

Peleg

Peleg (fixed k)

Parameter

Estimate

ME

RMSE

AICc

D52.5 (min)

49.96 (67.62) 4.15 (2.09) 6.00 (0.19)

−0.02

4.35

−35.75

29.55 (1.60) 5.18 (*) 6.06 (0.14)

−0.02

4.51

−38.53

4.81 (12.10) 0.79 (0.10) 57.54 (0.16) 5.91 (0.18)

0.00

3.72

−32.21

0.58 (*) 1.00 (0.14) 59.47 (0.96) 6.05 (0.18)

−0.02

4.47

−35.87

z (°C) logN0 (log CFU/ ml) D52.5 (min) z (°C) logN0 (log CFU/ ml) k (1/°C) n (−) Tcrit (°C) logN0 (log CFU/ ml) k (1/°C) n (−) Tcrit (°C) logN0 (log CFU/ ml)

Table 4 Estimated values (and standard deviations) of the model parameters of each model implemented in Bioinactivation FE for the isothermal inactivation of Salmonella Senftenberg for different concentrations of nanoemulsified D-limonene. Control

0.01 mM

0.05 mM

0.1 mM

54.51 (3.87) 4.15 (0.18) 39.52 (4.92) 4.02 (0.21) 0.70 (0.07) 0.45 (0.03) 0.67 (0.06) 55.50 (0.41)

13.69 (1.10) 3.07 (0.09) 4.76 (1.21) 3.03 (0.12) 0.42 (0.06) 0.58 (0.07) 0.42 (0.06) 50.67 (0.35)

2.25 (0.15) 3.65 (0.12) 0.75 (0.18) 4.54 (0.49) 0.31 (0.06) 0.54 (0.08) 0.49 (0.06) 48.71 (0.27)

1.79 (0.08) 3.34 (0.07) 0.34 (0.09) 3.04 (0.08) 0.38 (0.04) 0.47 (0.05) 0.21 (0.03) 45.86 (0.50)

4. Conclusions Bigelow Mafart

Peleg

D50 (min) z (°C) δ50 (min) z (°C) p (−) k (1/°C) n (−) Tcrit (°C)

Bioinactivation FE is a user-friendly free web application for modelling microbial inactivation. It improves Bioinactivation SE with additional functions (model fitting of isothermal experiments and generation of predictions) and new developments in the user interface. The functions implemented in the tool have been showcased using two case studies. The first one illustrates how Bioinactivation FE can be used to compare the inactivation pattern of bacterial cells during isothermal and non-isothermal treatments. The functions included in the application can be used to characterize the microbial response observed during one type of treatment (isothermal or dynamic). Then, the model parameters obtained can be used to simulate the expected microbial concentration during different treatments (isothermal or dynamic) and compare it with experimental values. The second case study shows how the software can be used to analyse the effect that different concentrations of additives have on the thermal resistance of bacterial cells. Moreover, this case study illustrates how the statistical information provided by the tool can be used for discriminating between inactivation models. In conclusion, the functions included Bioinactivation FE make accessible advanced mathematical modelling techniques to food scientists without any expertise in scientific programming.

Table 5 Statistical indexes describing the quality of the model fit of Salmonella Senftenberg for different concentrations of nanoemulsified D-limonene.

AICc

Bf

ME

RMSE

Af

Bigelow Mafart Peleg Bigelow Mafart Peleg Bigelow Mafart Peleg Bigelow Mafart Peleg Bigelow Mafart Peleg

Control

0.01 mM

0.05 mM

0.1 mM

178.29 163.28 157.3 0.78 1.00 1.02 −0.11 0.00 0.01 0.59 0.54 0.53 3.92 3.5 3.37

82.73 43.13 43.13 0.61 1.00 1.00 −0.21 0.00 0.00 0.59 0.37 0.42 3.89 2.33 2.64

254.95 187.05 174.52 0.53 1.01 0.98 −0.28 0.00 −0.01 0.75 0.55 0.37 5.58 3.59 2.33

127.82 45.36 71.98 0.56 1.00 1.02 −0.25 0.00 0.01 0.69 0.34 0.52 4.9 2.17 3.34

Acknowledgements The financial support of this research work was provided by the Ministry of Economy and Competitiveness (MINECO) of the Spanish Government and European Regional Development Fund (ERDF) through projects AGL2013-48993-C2-1-R and DPI2011-28112-C04-04. Alberto Garre (BES-2014-070946) is grateful to the MINECO for awarding him a pre-doctoral grant.

was fitted using the one-step algorithm included in the tool. Table 4 reports the estimated parameter values and standard deviations, whereas Table 5 includes the values of the statistical indexes describing the quality of the model fit. The limitations of the Bigelow model to describe this data set are shown by the statistical indexes reported in Table 5. In every case, the ME is lower than one, indicating that the fitted curve is biased with

Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.foodres.2018.06.057. 359

Food Research International 112 (2018) 353–360

A. Garre et al.

References

Lin, Y.-D., & Chou, C.-C. (2004). Effect of heat shock on thermal tolerance and susceptibility of Listeria monocytogenes to other environmental stresses. Food Microbiol. 21(5), 605–610. http://dx.doi.org/10.1016/j.fm.2003.10.007. Lindqvist, R., Sylvén, S., & Vågsholm, I. (2002). Quantitative microbial risk assessment exemplified by Staphylococcus aureus in unripened cheese made from raw milk. Int. J. Food Microbiol. 78(1–2), 155–170. http://dx.doi.org/10.1016/S0168-1605(02) 00237-4. Mafart, P., Couvert, O., Gaillard, S., & Leguerinel, I. (2002). On calculating sterility in thermal preservation methods: application of the Weibull frequency distribution model. Int. J. Food Microbiol. 72(1–2), 107–113. http://dx.doi.org/10.1016/S01681605(01)00624-9. McMeekin, T. A., Mellefont, L. A., & Ross, T. (2007). Predictive microbiology: past, present and future. Modelling Microorganisms in Food (pp. 7–21). Elsevier. http://dx. doi.org/10.1533/9781845692940.1.7. McMeekin, T. A., Olley, J., Ratkowsky, D. A., & Ross, T. (2002). Predictive microbiology: towards the interface and beyond. Int. J. Food Microbiol. 73(2–3), 395–407. http://dx. doi.org/10.1016/S0168-1605(01)00663-8. Moré, J. J. (1978). The Levenberg-Marquardt algorithm: Implementation and theory. Numerical Analysis (pp. 105–116). Berlin, Heidelberg: Springer. http://dx.doi.org/10. 1007/BFb0067700. Muñoz-Cuevas, M., Metris, A., & Baranyi, J. (2012). Predictive modelling of Salmonella: From cell cycle measurements to e-models. Food Res. Int. 45(2), 852–862. http://dx. doi.org/10.1016/j.foodres.2011.04.033. Nunes, C. A., Alvarenga, V. O., de Souza Sant'Ana, A., Santos, J. S., & Granato, D. (2015). The use of statistical software in food science and technology: Advantages, limitations and misuses. Food Res. Int. 75, 270–280. http://dx.doi.org/10.1016/j.foodres. 2015.06.011. Peleg, M., & Cole, M. B. (1998). Reinterpretation of Microbial Survival Curves. Crit. Rev. Food Sci. Nutr. 38(5), 353–380. http://dx.doi.org/10.1080/10408699891274246. Plaza-Rodríguez, C., Thoens, C., Falenski, A., Weiser, A. A., Appel, B., Kaesbohrer, A., & Filter, M. (2015). A strategy to establish Food Safety Model Repositories. Int. J. Food Microbiol. 204, 81–90. http://dx.doi.org/10.1016/j.ijfoodmicro.2015.03.010. Richter, K., Haslbeck, M., & Buchner, J. (2010). The Heat Shock Response: Life on the Verge of Death. Mol. Cell, 40(2), 253–266. http://dx.doi.org/10.1016/j.molcel.2010. 10.006. Ros-Chumillas, M., Garre, A., Maté, J., Palop, A., & Periago, P. M. (2017). Nanoemulsified D-Limonene Reduces the Heat Resistance of Salmonella Senftenberg over 50 Times. Nanomaterials, 7(3), 65. http://dx.doi.org/10.3390/nano7030065. Ross, T. (1996). Indices for performance evaluation of predictive models in food microbiology. J. Appl. Microbiol. 81(5), 501–508. Shapiro, S. S., & Wilk, M. B. (1965). An analysis of variance test for normality (complete samples). Biometrika, 52(3–4), 591–611 https://doi.org/10.1093/biomet/52.3-4. 591. Soetaert, K., & Petzoldt, T. (2010). Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME. J. Stat. Softw. 33(3), http://dx.doi.org/10.18637/ jss.v033.i03. Soetaert, K., Petzoldt, T., & Setzer, R. W. (2010). Solving differential equations in R: Package deSolve. J. Stat. Softw. 33(9), http://dx.doi.org/10.18637/jss.v033.i09. Stasiewicz, M. J., Marks, B. P., Orta-Ramirez, A., & Smith, D. M. (2008). Modelling the effect of prior sublethal thermal history on the thermal inactivation rate of salmonella in ground Turkey. J. Food Prot. 71(2), 279–285. http://dx.doi.org/10.4315/0362028X-71.2.279. Tenenhaus-Aziza, F., & Ellouze, M. (2015). Software for predictive microbiology and risk assessment: A description and comparison of tools presented at the ICPMF8 Software Fair. Food Microbiology, 45(Part B), 290–299. http://dx.doi.org/10.1016/j.fm.2014. 06.026. Valdramidis, V. P., Geeraerd, A. H., Bernaerts, K., & Van Impe, J. F. (2006). Microbial dynamics versus mathematical model dynamics: The case of microbial heat resistance induction. Innovative Food Sci. Emerg. Technol. 7(1–2), 80–87. http://dx.doi.org/10. 1016/j.ifset.2005.09.005. Valdramidis, V. P., Geeraerd, A. H., & Van Impe, J. F. (2007). Stress-adaptive responses by heat under the microscope of predictive microbiology: Modelling the microbial heat resistance. J. Appl. Microbiol. 103(5), 1922–1930. http://dx.doi.org/10.1111/j.13652672.2007.03426.x. Verhulst, P.-F. (1838). Notice sur la loi que la population suit dans son accroissement. Correspondance mathématique et physique publiée par a. Quetelet. 10, 113–121. Vilas, C., Arias-Mendez, A., Garcia, M. R., Alonso, A. A., & Balsa-Canto, E. (2018). Towards predictive food process models: A protocol for parameter estimation. Crit. Rev. Food Sci. Nutr. 58(3), 436–449. http://dx.doi.org/10.1080/10408398.2016. 1186591.

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Autom. Control, 19(6), 716–723. Balsa-Canto, E., Alonso, A. A., & Banga, J. R. (2010). An iterative identification procedure for dynamic modelling of biochemical networks. BMC Syst. Biol. 4(1), 11. Baranyi, J., Pin, C., & Ross, T. (1999). Validating and comparing predictive models. Int. J. Food Microbiol. 48(3), 159–166. http://dx.doi.org/10.1016/S0168-1605(99) 00035-5. Bates, D. M., & Watts, D. G. (2007). Nonlinear Regression Analysis and Its Applications (1 edition). New York, NY: Wiley-Interscience. Bigelow, W. D. (1921). The logarithmic nature of thermal death time curves. J. Infect. Dis. 29(5), 528–536. Bigelow, W. D., & Esty, J. R. (1920). The thermal death point in relation to time of typical thermophilic organisms. J. Infect. Dis. 602–617. Box, G. E. P., Hunter, J. S., & Hunter, W. G. (2005). Statistics for Experimenters: Design, Innovation, and Discovery. Hoboken, N.J: Wiley-Blackwell. Chang, W., Cheng, J., Allaire, J. J., Xie, Y., & McPherson, J. (2017). Shiny: web application framework for R. Retrieved from https://CRAN.R-project.org/package=shiny. Chang, W., & Ribeiro, B. B. (2017). Shinydashboard: Create Dashboards with “Shiny”. Retrieved from https://CRAN.R-project.org/package=shinydashboard. Coleman, M. E., & Marks, H. M. (1999). Qualitative and quantitative risk assessment. Food Control, 10(4–5), 289–297. http://dx.doi.org/10.1016/S0956-7135(99)00052-3. Conesa, R., Andreu, S., Fernández, P. S., Esnoz, A., & Palop, A. (2009). Nonisothermal heat resistance determinations with the thermoresistometer Mastia. J. Appl. Microbiol. 107(2), 506–513. http://dx.doi.org/10.1111/j.1365-2672.2009.04236.x. Core Team, R. (2016). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from https://www.Rproject.org. Esty, J. R., & Meyer, K. F. (1922). The heat resistance of the spores of B. botulinus and Allied anaerobes. XI. J. Infect. Dis. 31(6), 650–664. Fernández, A., Ocio, M. J., Fernández, P. S., Rodrigo, M., & Martinez, A. (1999). Application of nonlinear regression analysis to the estimation of kinetic parameters for two enterotoxigenic strains ofBacillus cereus spores. Food Microbiol. 16(6), 607–613. http://dx.doi.org/10.1006/fmic.1999.0282. Garre, A., Fernández, P. S., Lindqvist, R., & Egea, J. A. (2017). Bioinactivation: Software for modelling dynamic microbial inactivation. Food Res. Int. 93, 66–74. http://dx.doi. org/10.1016/j.foodres.2017.01.012. Garre, A., Huertas, J. P., González-Tejedor, G. A., Fernández, P. S., Egea, J. A., Palop, A., & Esnoz, A. (2018). Mathematical quantification of the induced stress resistance of microbial populations during non-isothermal stresses. Int. J. Food Microbiol. 266, 133–141. http://dx.doi.org/10.1016/j.ijfoodmicro.2017.11.023. Geeraerd, A. H., Herremans, C. H., & Van Impe, J. F. (2000). Structural model requirements to describe microbial inactivation during a mild heat treatment. Int. J. Food Microbiol. 59(3), 185–209. http://dx.doi.org/10.1016/S0168-1605(00)00362-7. Granato, D., de Araújo Calado, V. M., & Jarvis, B. (2014). Observations on the use of statistical methods in Food Science and Technology. Food Res. Int. 55, 137–149. http://dx.doi.org/10.1016/j.foodres.2013.10.024. Granato, D., Nunes, D. S., & Barba, F. J. (2017). An integrated strategy between food chemistry, biology, nutrition, pharmacology, and statistics in the development of functional foods: A proposal. Trends Food Sci. Technol. 62, 13–22. http://dx.doi.org/ 10.1016/j.tifs.2016.12.010. Haario, H., Laine, M., Mira, A., & Saksman, E. (2006). DRAM: Efficient adaptive MCMC. Stat. Comput. 16(4), 339–354. http://dx.doi.org/10.1007/s11222-006-9438-0. Hassani, M., Manas, P., Pagán, R., & Condón, S. (2007). Effect of a previous heat shock on the thermal resistance of Listeria monocytogenes and Pseudomonas aeruginosa at different pHs. Int. J. Food Microbiol. 116(2), 228–238. Hassani, M., Mañas, P., Condón, S., & Pagán, R. (2006). Predicting heat inactivation of Staphylococcus aureus under nonisothermal treatments at different pH. Mol. Nutr. Food Res. 50(6), 572–580. http://dx.doi.org/10.1002/mnfr.200500171. Hill, C., Cotter, P. D., Sleator, R. D., & Gahan, C. G. M. (2002). Bacterial stress response in Listeria monocytogenes: jumping the hurdles imposed by minimal processing. Int. Dairy J. 12(2), 273–283. http://dx.doi.org/10.1016/S0958-6946(01)00125-X. Hindmarsh, A., Odepack, A., Stepleman, R. S., et al. (1983). Systematized collection of ODE Solvers. North-Holland, Amsterdam. vol. 1. North-Holland, Amsterdam (pp. 55– 64). IMACS Transactions on Scientific Computation (1, 55–64). Huang, L. (2017). IPMP Global Fit – A one-step direct data analysis tool for predictive microbiology. Int. J. Food Microbiol. 262(Supplement C), 38–48. http://dx.doi.org/ 10.1016/j.ijfoodmicro.2017.09.010. Hurvich, C. M., & Tsai, C.-L. (1993). A corrected Akaike information criterion for vector autoregressive model selection. J. Time Ser. Anal. 14(3), 271–279.

360