Ten predictive growth rate models, that included the effects of temperatures below 12°C, pH and a_w on the growth rate of either psychrotolerant or mesophilic-psychrotolerant intermediary thermotypes of B. cereus s.l., were extracted from the scientific literature. One model was from ComBase: A Combined Database For Predictive Microbiology (RRID:SCR_008181), which are partly based on the work by Sutherland et al. (1996), while of the remaining models eight came from Carlin et al. (2013) and one from Zwietering et al. (1996) (Table 5). These models were used to predict responses for different B. cereus s.l. isolates based on product characteristics and storage temperature as described in Table 2.

Carlin et al. (2013) developed cardinal parameter models including the effect of temperature (T), pH and a_w (Equations 4 and 5) on μ_max–values of B. cereus s.l. isolates from different phylogenetic panC groups. The eight models included in the present study (Table 5) had T_min–values from 1.4 to 9.1°C, pH_min–values from 4.59 to 4.96 and a_{w min}–values from 0.946 to 0.973 (Carlin et al., 2013).

For growth of naturally occurring B. cereus s.l. in milk, Zwietering et al. (1996) suggested a cardinal parameter model like Equation 4 and with terms for temperature, pH and a_w that were simpler than indicated by Equation 5. This model used T_min = 0.0°C; pH_min = 4.9 and a_{w min} = 0.95.

These 10 growth rate models were evaluated by using 21 growth/no-growth responses and corresponding product characteristics as determined in the present study (see Section 2.1 and Table 2). This screening of growth rate models was used to exclude the models with poor, or no potential of improved performance by product calibration and/or expansion with terms for interactions between T, pH and a_w.

The performance of growth rate models was evaluated by comparison of observed and predicted μ_max–values. Bias factor (B_f; Equation 6) and accuracy factor (A_f; Equation 7) values were calculated and compared with limits previously used for evaluating growth rate models for various bacteria: 0.95 < B_f < 1.11 indicate a good model performance, with B_f in the intervals of 1.11–1.43 or 0.87–0.95 corresponding to acceptable model performance and B_f < 0.87 or > 1.43 considered as unacceptable model performance (Mejlholm et al., 2010).

In addition, A_f > 1.5 indicate poor model precision or a systematic deviation between observed and predicted μ_max–values (Mejlholm and Dalgaard, 2013). Predicted and observed growth and no-growth responses were evaluated by calculating the percentage of samples that were correctly predicted. Incorrect predictions were considered as fail-safe (growth predicted when no-growth was observed) or fail-dangerous (no growth predicted when growth was observed). Criteria corresponding to good, acceptable and unacceptable model performance have not been established for the percentage of correct, fail-safe and fail-dangerous predictions. Nevertheless, when evaluating different models and using the same data set these indices allow the performance of models to be ranked. Larger validation studies found the better models to have >75% correct, < 15% fail-safe and < 10% fail-dangerous predictions (Koukou et al., 2022; Martinez-Rios et al., 2020; Mejlholm et al., 2010). Ideally, models should provide 100% correct, 0% fail-safe and 0% fail-dangerous predictions but when product characteristics are close to the growth boundary a small percentage of fail-safe and fail-dangerous predictions can be observed, even for precise models, due to for example variability in product characteristics. Therefore, it is particularly important to indicate if fail-dangerous predictions are close to the growth boundary and this can be done by using the ψ–value (see section 2.4, Equation 12) which has a value of 1.0 at the growth boundary (Mejlholm and Dalgaard, 2009).

Even though B_f –value, A_f –value and proportion of correct, fail-safe and fail-dangerous predictions of growth/no-growth responses are normally used as model performance indices, they were calculated here to be able to select the most promising models. Different approaches for model improvement were applied. For each of the selected literature models, 82 new models were developed in the following way. One model was developed by product calibrating the μ_opt–value using Equation 8 to create μ_opt-C. By keeping the μ_opt–value unchanged and expanding the model with the interaction term ξ (Equation 9) using three different values of n in Equation 11, for the effect of interaction between temperature, pH and a_w, resulted in 3³ = 27 different models. Another 2 × 27 = 54 models were developed, 27 by first product calibrating μ_opt and then expanding with the interaction term ξ and 27 by first expanding with the interaction term ξ and then product calibrating μ_opt.

Product calibration of the selected Carlin et al. (2013) and Zwietering et al. (1996) models was performed by dividing the original μ_opt –values for each model with the B_f –value determined for the specific model (Equation 8).

where μ_opt-C is the maximum specific growth rate (h⁻¹) after product calibration and at the optimum growth temperature as suggested by Koukou et al. (2021).

ξ in Equation 9 described the effect of interactions between the environmental factors and its effect was modeled as previously reported by using the Le Marc approach (Le Marc et al., 2002). The value of ξ was between 0 and 1 and calculated according to Equations 10–12.

where ξ(ϕ(T, pH, a_w) is the term describing the effects of interactions between environmental factors on μ_max. For temperature, pH and a_w the contribution of each of these terms in Equation 9 to the interaction term (ξ, Equation 10) was calculated by using Equation 11 and Equation 12. Le Marc et al. (2002) applying a value of 2.0 for n in Equation 11, however, in the present study the effect of using values of 1, 2 or 3 was evaluated as described below.

where e_i represents the environmental factors and ϕ_e the contribution of each environmental term to the effect of interactions between the factors.

The ψ–value provides a measure of how far a specific set of environmental factors is from the growth boundary (Mejlholm and Dalgaard, 2009) and a ψ–value higher than 1.0 indicated no growth (Equation 10).

To find the most promising interaction terms for the effect of T, pH and a_w in the better performing models, all the 27 combinations, which resulted from using values of 1, 2, or 3 for n in Equation 11 when used for T, pH or a_w, respectively, were tested. For each model, the B_f – and A_f –values as well as the percentages of correct, fail-safe and fail-dangerous predictions of growth/no-growth responses were calculated using the dataset in Table 2. This approach was performed both before and after product calibration of the models. Models with good or acceptable B_f – and A_f –values (see section 2.3), ≥ 75% correct and ≤ 5% fail-dangerous predictions for growth/no-growth responses were selected for further evaluation.

All new models, constructed as described in Section 2.4 and fulfilling the acceptability criteria for performance, were selected as promising models. Performance of these promising models were evaluated using the independent data reported in Table 3. B_f – and A_f –values (Equations 6 and 7) and proportion of correct, fail-safe and fail-dangerous predictions of growth/no-growth responses were calculated and used as model performance indices with the purpose of selecting two of the models for further evaluation using growth/no-growth responses reported in the scientific literature.

A total of 33 kinetic responses, for more than 12 different isolates of B. cereus s.l., in a range of single component and composite starchy foods were extracted from five available studies and 10 ComBase records. The studies with single component starchy foods included mashed potatoes from powder, cooked rice, cooked noodles, sliced bread and potato purée whereas the composite starchy foods included meat loaf, composite fried rice meal, pizza, meat lasagna, cottage pie and vegetable pie (Table 6). Exclusively, responses reported for storage temperatures of max 12°C were studied and exclusively for psychrotolerant or mesophilic-psychrotolerant intermediary thermotypes, i.e., strains belonging to the phylogenetic (panC) groups II, IV, V, and VI. When phylogenetic groups were not reported, then the thermo-type of strains was considered psychrotolerant when growth was observed below 10°C or mesophilic-psychrotolerant intermediary when growth was observed at 10°C but not below (Table 6). Product characteristics (pH, NaCl/a_w) and storage temperature were recorded for kinetic responses extracted from literature (Table 6). When no information regarding NaCl/a_w or pH was provided, then an average value was assumed from reported values for a similar type of food. Other environmental factors, including organic acids, were not mentioned for any of the eight studies analyzed and, therefore, assumed not to be present. When a_w was not reported, it was estimated using concentrations of NaCl and moisture to determine % WPS (Equation 13) and converting this to a_w using Equation 14 (Resnik and Chirife, 1988; Ross and Dalgaard, 2004) shown below.

where WPS is water phase salt.

Evaluation of the two most promising models were performed using the scientific literature data shown in Table 6 and the composite starchy foods from the present study reported in Table 4. For these evaluations, B_f – and A_f –values (Equations 6 and 7) and proportion of correct, fail-safe and fail-dangerous predictions of growth/no-growth responses were calculated and used as model performance indices with the purpose of selecting the best model to use for predicting safe shelf-lives.

The initial 21 challenge tests, which included nine B. cereus s.l. isolates, five starchy food products, storage temperatures from 6.6 to 11.7°C, pH 4.8–7.8 and % WPS of 0.02–9.0 with a_w–values of 0.935–0.999, resulted in nine tests with no-growth responses and 12 tests with average μ_max –values of 0.016–0.200 h⁻¹ (Table 2). When compared to the experimental data, two models were excluded from further evaluation due to observed growth at temperatures ≤7.7°C (Table 2), which is lower than the T_min–value of strains F4430/73 (T_min = 9.1°C) and ATCC 14579 (T_min = 7.8°C) (Carlin et al., 2013). Another two of the 10 studied models from the literature were excluded from further evaluation due to observed growth in bulgur adjusted to a_w of 0.957 and in pasta adjusted to a_w of 0.962 (Table 2), which were below the a_{w min} –values of 0.964 and 0.973 for the strains KBAB4 and ADRIA I21, respectively (Carlin et al., 2013). The model from ComBase was also excluded from further study. This model had a proportion of correct growth/no-growth predictions of 65% (Table 5) being lower than 75% which is aimed for in validation studies. Moreover, as this model is not a cardinal parameter-type model, and its model parameter values are not known (ComBase, 2024), it was not possible to expand the model with an interaction term to improve the percentage of correct and fail-safe predictions.

The five remaining models were the models for group II strains RIVM BC120 and NVH 0862–00, group V strains F2769/77 and NVH 141 (Carlin et al., 2013) and the model in Zwietering et al. (1996). They are all cardinal parameter models and had no fail-dangerous predictions of the growth/no-growth responses (Table 5). Two of these models, i.e., the Carlin et al. (2013) model for group II strain RIVM BC120 and the model in Zwietering et al. (1996), both had B_f – and A_f –values that were well above 1.0 and close to each other (Table 5). This does not necessarily disqualify these models from further studies as such a situation has previously been solved by product calibration of the model where the μ_opt–value is calibrated to include the effect of specific foods (Dalgaard and Mejlholm, 2019; Koukou et al., 2021; Rosso et al., 1996). High B_f – and/or A_f –values have also been linked to evaluation of predictive growth rate models without an interaction term (Mejlholm and Dalgaard, 2009). When approaching the growth boundary, growth rates are often reduced due to the interaction between environmental factors and if not accounted for in a predictive model, then growth rates can be over-predicted resulting in increased B_f – and A_f –values and a high proportion of fail-safe predictions of the growth/no-growth responses (Mejlholm and Dalgaard, 2009). Neither the Carlin et al. (2013) models nor the Zwietering et al. (1996) model included an interaction term, suggesting expansion of these models with this term could be an option to decrease the number of fail-safe and increase the number of correct predictions of the growth/no-growth responses. Importantly, Carlin et al. (2013) and Zwietering et al. (1996) models were developed in markedly different ways. Carlin et al. (2013) developed models for growth rates of individual isolates in BHI broth whereas Zwietering et al. (1996) estimated cardinal parameter values from data for growth of naturally occurring B. cereus in milk.

Of the overall 410 models (82 modifications of each of the five original models) that were tested using the growth responses in Table 2, nine models complied with the criteria of having a good or an acceptable B_f –value (0.87 ≤ B_f ≤ 1.43), an acceptable A_f –value (A_f ≤ 1.5) and at the same time resulting in ≤5% fail-dangerous and ≥ 75% correct predictions of growth/no-growth responses (Table 7). All nine of the best performing models were derived from the two original models with the lowest T_min–values, i.e., 1.4°C for the group II strain RIVM BC120 model (Carlin et al., 2013) and 0.0°C for the model in Zwietering et al. (1996) (Table 5). As growth responses at low temperatures was predicted this is probably to be expected. Importantly, all these nine models included a term for the growth inhibiting effect of interactions between the environmental factors, temperatures, pH and a_w (Table 7). This confirmed the importance of taking the effect of interactions into account when growth responses are predicted as previously observed, e.g., for mesophilic B. cereus (Le Marc et al., 2021), L. monocytogenes (Augustin and Carlier, 2000; Le Marc et al., 2002; Mejlholm and Dalgaard, 2009) and non-proteolytic Clostridium botulinum (Koukou et al., 2021).

The μ_opt–value was calibrated for eight of the nine best performing models (Table 7) indicating that growth rates of B. cereus s.l. in starchy foods differ from growth rates in the BHI broth or milk used for development of the original models (Carlin et al., 2013; Zwietering et al., 1996). Only, one of nine best performing models included the original μ_opt–value from milk (Table 7; Zwietering et al., 1996).

As changing the μ_opt –value does not affect the cardinal parameter values, which define the growth/no-growth conditions, product calibration as the sole approach for updating models (i.e., approach i, see section 2.4) could not improve the predictions of the growth/no-growth responses. However, product calibration in combination with expanding the original models with an interaction term (i.e., approaches iii and iv, see section 2.4) was very effective, as resulting in acceptable performance of eight of nine models in Table 7 which were updated in this way. Half of these eight models were a result of first calibrating the μ_opt –value of the original models and then expanding the model with an interaction term, while the other four models were expanded with an interaction term before calibration of μ_opt –values (Table 7). Hence, no clear picture on best practice could be deduced indicating that it could be depending on the specific model.

The six models, originating from Carlin et al. (2013) group II strain RIVM BC120 model, kept 0% fail-dangerous predictions after updating, suggesting the best model for the purpose of predicting the growth/no-growth response of B. cereus s.l. in starchy foods would be found among these six. However, the 5% fail-dangerous predictions of growth/no-growth responses obtained for the updated Zwietering et al. (1996) models, do not necessarily disqualify these models from further studies. The result of 5% fail-dangerous predictions of no-growth, when growth was observed, corresponded to one of the challenge tests, i.e., Exp. no. 20 (Table 2). In this challenge, with cooked bulgur inoculated with group VI strain ADRIA I21, the ψ–value was determined to be 1.03 (Table 7) with the updated models when the average measurements of storage temperature of 11.7°C, pH of 6.8 and a_w of 0.958 (Table 2) were applied. As a value of 1.0 indicates the growth boundary (Le Marc et al., 2024; Mejlholm and Dalgaard, 2009), the value 1.03 predicted no-growth. As the a_w –value of this challenge is close to the a_{w min} –value of 0.950 for the model even small deviations in the a_w measurement can change the ψ–value from predicting no-growth to predicting growth, e.g., using the highest a_w –value of 0.959, measured in this case, would have resulted in a ψ–value of 0.96, thereby, predicting a growth response in this sample. Therefore, none of the three models originating from the Zwietering et al. (1996) were disqualified and all nine updated models were studied further as promising candidates being evaluated using the independent experimental data shown in Table 3.

The 21 performed challenge tests, used as the independent growth/no-growth responses (Table 3) in the evaluation, included the same nine B. cereus s.l. isolates and the same five starchy food products as used for the updating of literature models. Regarding storage temperature, pH, % WPS and measured a_w, all levels were within the ranges used for updating the models (Tables 2, 3). Challenge tests resulted in ten no-growth responses and 11 growth responses with average μ_max–values of 0.011–0.128 h⁻¹ (Table 3).

Four of the nine studied models performed with B_f–values between 0.87 and 1.01 as well as A_f–values between 1.15 and 1.32 indicating acceptable to good performance for prediction of μ_max–values (h⁻¹) (models in bold, Table 8). Of these four models, two stood out with better results for the prediction of the growth/no-growth responses and resulted in more than 75% correct, less than 25% fail-safe and no fail-dangerous predictions (models with *, Table 8). This was an improvement of the number of correct predictions of 14 percentage points for both the Carlin et al. (2013) model for group II strain RIVM BC120 as well as for the model in Zwietering et al. (1996) (Tables 5, 8). Figure 1 compares observed and predicted μ_max–values for these two models and illustrates, that predictions obtained using the updated Zwietering et al. (1996) model, on average were less biased with equal number of data points scattered around the line of perfect match, though with two results positioned further above the line as indicated with square symbols in Figure 1 (□). These two growth responses were both from challenge tests with the panC group VI strain ADRIA I21 in samples with no added NaCl (a_w, 0.999) having pH–values (6.2 and 6.6) around the optimal of 6.4 for this strain (Carlin et al., 2013) and stored at 7.7 and 11.5°C, respectively (Table 3). However, as shown in Figure 1 (■), these two growth responses did not deviate markedly from the predicted μ_max–value when using the updated Carlin et al. (2013) model for panC group II strain RIVM BC120. The difference appeared to be related to the environmental term for a_w, CM₁(a_w), where the predicted values for the updated Carlin et al. (2013) model for panC group II strain RIVM BC120 were lower compared to values obtained using the updated Zwietering et al. (1996) model (results not shown). The cardinal parameter a_{w opt} has the value of 1.0 in the Zwietering et al. (1996) model meaning that a_{w max} becomes irrelevant. With a_{w opt} and a_{w max} of 0.997 and 1.0, respectively, in the Carlin et al. (2013) model for panC group II strain RIVM BC120, the consequence is that a_w–values of more than 0.9985 will result in a shift toward lower CM₁(a_w) and lower predicted μ_max–values for this model compared to the Zwietering et al. (1996) model.

Nevertheless, both models performed within an acceptable range for starchy foods, and both were additionally evaluated using growth responses partly reported in the scientific literature (Table 6) (n = 33), partly generated in the present study (Table 4) (n = 8).

The dataset extracted from the scientific literature consisted of 33 growth/no-growth responses, leading to the acquisition of 21 μ_max –values and 12 no-growth responses to be included in the evaluation (Table 6). Different sub-datasets of growth/no-growth responses were created for the evaluation based on data in Table 6; one for each of the two thermotypes psychrotolerant (n = 13) and intermediary (n = 18), one for single component starchy foods (n = 16), one for composite starchy foods (n = 17), one for meat loaf (n = 6) and a final sub-dataset excluding the data from meat loaf (n = 27) (Table 9). Using the 21 μ_max–values, the updated model from Carlin et al. (2013) for panC group II strain RIVM BC120 performed better than the updated Zwietering et al. (1996) model with B_f – and A_f –values closer to the acceptance criteria (Table 9). This applied regardless of sub-dataset suggesting a systematic difference resulting in generally higher μ_max predictions for the updated Zwietering et al. (1996) model. When looking closer into this difference, it revealed that the predicted ψ–values were lower for 27 out of the 33 challenge tests when applying the updated Zwietering et al. (1996) model (data not shown) indicating a lower dampening effect of the interaction term compared to the updated Carlin et al. (2013) model for panC group II strain RIVM BC120. On the other hand, the higher dampening effect, seen for the updated Carlin et al. (2013) model for panC group II strain RIVM BC120, resulted in two predictions of fail-dangerous growth responses (Table 9). When looking closer into these two fail-dangerous predictions, both cases appeared to concern pizzas, one with pH 5.1 having a_w –value of 0.983 (ComBase, 2023b) and another with pH 5.1 and a_w–value of 0.994 (ComBase, 2023a), and both had been stored at 10°C (Table 6). These pH, a_w and temperature conditions resulted in predicted ψ–values of 1.04 and 1.02 (Table 9). So, both were very close to the growth boundary at 1.0, which means that even small uncertainties in the product characteristics or in the storage temperature could change the prediction from a no-growth response to a growth response. For these two specific observations, e.g., a change in pH–value to 5.13 or a change in storage temperature to 10.4°C, would change the ψ–values to become less than 1.0, moving these fail-dangerous no-growth responses to correct growth responses. These relatively small changes are within the uncertainties that would be expected for pH and temperature measurements when conducting challenge tests (Tables 2, 3). Therefore, care should be taken when disqualifying models exclusively based on data where uncertainties for intrinsic and extrinsic factors are not reported. Taking this into consideration, the results in Table 9 pointed at the updated Carlin et al. (2013) model for panC group II strain RIVM BC120 as the less biased and more accurate of the two models for predicting growth of B. cereus s.l. The model performed with an overall acceptable B_f –value of 1.34, an A_f –value of 1.57 close to being acceptable and with 70% correct predictions of growth/no-growth responses, classifying the model as generally fail-safe for foods containing starchy ingredients and stored at max 12°C.

As shown in Table 9, the B_f – and A_f –values, obtained using the updated Carlin et al. (2013) model for panC group II strain RIVM BC120, were better for the sub-dataset of challenge tests using intermediary thermotypes (panC groups IV and V) and for the sub-dataset of challenge tests involving single component starchy foods. For both sub-datasets, the evaluation was based on n = 7 growth responses (Table 9) but only for one growth response, these two sub-datasets overlapped, i.e., only one μ_max observation was found for an intermediary thermotype (B4ac-1) in a single component starchy food (mashed potatoes from powder) (Mahakarnchanakul and Beuchat, 1999; Table 6), whereas six were found for intermediary thermotypes in composite starchy foods and six for psychrotolerant thermotypes in single component starchy foods. This indicated that n = 13 (n = 6 + 6 + 1) of the n = 19 μ_max predictions, used in total for this evaluation, actually were less biased (B_f, 1.14) and more accurate (A_f, 1.43) than the overall averages (B_f/A_f, 1.34/1.57) (Table 9). Consequently, the remaining n = 6 (n = 19–13) μ_max predictions represented the combination of psychrotolerant thermotypes in composite starchy foods. Applying this sub-dataset, which turned out to be the six challenges conducted for meat loaf, resulted in B_f – and A_f –values of 1.94, which were much higher than the overall averages (Table 9), indicating that growth was strongly over-predicted by the updated Carlin et al. (2013) model for panC group II strain RIVM BC120. This over-prediction was unexpected, as other studies reported faster growth of B. cereus s.l. when animal proteins were available in the substrate as compared to cereal proteins (Ellouze et al., 2021; Kang et al., 2010; Morita and Woodburn, 1977). The additional challenge tests (Table 4) using products with animal or vegetable proteins, and some of the strains as were used for updating the models, were, therefore, included in the present study to investigate this matter. For both of the updated models, μ_max –values were strongly under-predicted with unacceptable B_f –values below 0.7 (Table 9). Of the six observed growth responses, four were even below the lower acceptable A_f –limit meaning that the observed μ_max –values were more than 1.5-fold higher than predictions (results not shown). Interestingly, three of these four low-scoring A_f challenges tests were from composite starchy foods rich in animal proteins and the remaining contained vegetable protein from split peas (Table 4). This confirms previous findings of growth rates of B. cereus s.l. in carbohydrate-rich foods being lower than in protein-rich foods, such as meat patties and tofu (Kang et al., 2010). Taken together this means that the updated Carlin et al. (2013) growth rate model for panC group II strain RIVM BC120 should not be used for composite protein-rich foods, as the growth rate might be under-predicted creating unsafe situations.

Knowing the time to reach a critical concentration of, e.g., 10⁵ cfu/g of B. cereus s.l. is an important input when deciding on the safe shelf-life for ready-to-eat or ready-to-cook chilled foods. The updated Carlin et al. (2013) model for panC group II strain RIVM BC120 (Table 9) can support this decision for foods consisting mainly of starchy ingredients, if the initial concentration (N₀) and the lag time are known (Equation 3). With μ_max –values predicted by the best performing model, lag times can be determined from the relative lag time (RLT) as Lag time = RLT × ln(2)/μ_max. RLT is often a constant (Mellefont and Ross, 2003; Ross, 1999) and lag time has been calculated in this way for different pathogens and foods (Dalgaard and Mejlholm, 2019). In the present study, RLT was estimated using data from all the challenge tests showing growth in single starchy foods after a statistically significant lag time (n = 58) (Supplementary Tables S3, S4). The median RLT-value of 7.2 (95%-CI, 1.6–40) was selected as a representative value for the predicted examples (results not shown). Product examples were chosen based on known N₀ of B. cereus s.l., i.e., concentrations measured close to the production time, as well as measured product characteristics (Table 10). The predictions in Table 10 demonstrated that keeping the storage temperature at max. 5°C was by far the most effective way of achieving a long safe shelf-life, i.e., at least 38 days. At this temperature, the time to reach the critical concentration of B. cereus s.l. was less affected by lowering pH or a_w than seen at the higher temperatures. Under storage at 8°C, the lowering of pH from 6.5 to 5.8 and a_w from 0.996 to 0.990 increased the time to reach the critical concentration of B. cereus s.l. with approx. 1.5–fold. The effect of having a low initial cell concentration can be seen when comparing rice with N₀ of 0.1 cfu/g to pasta with N₀ of 3 cfu/g. This showed that an approx. 10-fold lower N₀ resulted in 1, 2 and > 11 days longer time to reach the critical concentration of B. cereus s.l. at 10, 8 and 5°C, respectively (Table 10).