The original motivation of SINDy is to discover governing equations for nonlinear dynamical systems, which is reconfigured here to apply to non-dynamical systems as follows: y = f ( x ) , (1) where x is the vector of state variables, and f denotes the nonlinear relationship between the input ( x ) and output variables ( y ). SINDy approximates f by a weighted linear combination of nonlinear terms, e.g., for the i t h output variable: y i = f i ( x ) ≈ ∑ k θ k ( x )   ξ k , i , (2) where θ k ( x ) and ξ k , i denote the k t h term and its weight, respectively. The above equations can be represented in a more succinct form as matrices, i.e., Y = Θ ( X ) Ξ , (3) where X = [ x 1   x 2 ⋯ x m ] , Y = [ y 1   y 2 ⋯ y n ] , Θ ( X ) is a library of candidate functions of X and the matrix of weights Ξ = [ ξ 1   ξ 2 ⋯ ξ m ] . In SINDy, the library Θ ( X ) is built by polynomial expansion of input variables X , i.e., Θ ( X ) = [ 1   X   X 2 ⋯ X d ⋯ ] where X d denotes a matrix with column vectors of all possible d t h degree monomials in the state variable x .

SINDy seeks a parsimonious model composed with minimal number of terms as possible without compromising model accuracy. Sparse regression methods such as Sequentially Thresholded Least Squares (STLS) and Least Absolute Shrinkage and Selection Operator (LASSO) are useful algorithms that can be used in SINDy for this purpose (Brunton et al., 2016). In this work, we employ STLS where Ξ in Eq. 3 retains the coefficients (weights) greater than the prescribed parameter λ (otherwise, zero weights are assigned), such that only the terms in the library with significant influence on the outputs are included in the final model structure. Here, λ is known as sparsity-promoting knob because the model sparsity increases with higher values of λ , while model accuracy may decrease.

We use SINDy to formulate microbial inactivation as functions of various process variables including temperature ( T ), pH, water activity ( a w ), NaCl content ( C N ), and phosphate level ( C P ), which are all known to significantly influence microbial growth rate (Juneja et al., 1995; Cerf et al., 1996). We employ D-value (i.e., the time for microbial population to shrink to 10% of initial level) as a standard measure for microbial inactivation, which is taken as our target variable to predict in applying SINDy. With a single target variable chosen, Eq. 3 is reduced to the following equation, i.e., y = Θ ( X ) ξ . (4)

While SINDy offers flexibility to pick any nonlinear terms for input and output variables, we determine the inclusion of their specific functional forms following the known mechanistic knowledge and characteristics of the system. Therefore, we used a vector of log ⁡ D as y (instead of a vector of D ) and determined X to be [ 1 / T , pH, a w , C N , C P ] (i.e., the use of 1/T, instead of T ). The rationale for our choice of functional forms of output and input variables are detailed in Section 3.1.

We tune the order of combination of primitive process variables and the sparsity index, λ , in stages. We first determine the maximum order of combination with λ = 0 (which will result in a non-parsimonious model), beyond which there are no significant improvements to model accuracy. Subsequently, by retaining the maximum polynomial order, we employ the maximum λ that does not significantly compromise the model accuracy. To facilitate determining optimal polynomial order and λ values for a balanced compromise between model accuracy and sparsity, we use an information-theoretic metric, Akaike Information Criterion (AIC) (Akaike, 1998). Specifically, we use the second-order information criterion that includes a correction term to alleviate the bias that may arise if the number of model parameters is large relative to the sample datapoints (Burnham and Anderson, 2002): A I C = n ⁡ ln ( M S E ) + 2 K + 2 K ( K + 1 ) n − K − 1 , (5) where MSE denotes mean squared error, n is the number of sample datapoints, K is the number of model parameters and the third term on the RHS corrects the bias where it tends to zero when n ≫ K . Generally, a model with the least AIC score is ideal as AIC penalizes the model based on the relative balance between error and complexity ( K ). The formulation above is often denoted as AICc in the literature. The methodical implementation of the general guidelines to develop microbial inactivation models is demonstrated in Section 3.2.

We perform sensitivity analyses on our models as an alternative to arduous assessment of the relative effects of the process variables on microbial inactivation directly from highly distributed experimental data. As the models are linear combinations of nonlinear terms and the datasets used in this work span over a wide parameter space, the possibility of model forming stiff parameter dependency is high. Therefore, we employ a density-based global sensitivity analysis approach called PAWN (Pianosi and Wagener, 2015), instead of local sensitivity approach. Based on this approach, absolute deviation is calculated between an unconditional cumulative density function of model output, F ( y ) , where all input variables in the model are randomly sampled simultaneously over the whole parameter space, and conditional cumulative density functions, F ( y | X ¯ i , k ) , which are constructed by randomly sampling all but a single model variable of interest fixed at the k t h nominal value, X i = X ¯ i , k . The sensitivity index for i -th model variable, S i , is characterized as the maximum value across the distribution of absolute deviations collected for a range of k nominal values: S i = max X ¯ i , k [ K S ( X ¯ i , k ) ] ;   K S ( X ¯ i , k ) = max y | F ( y ) − F ( y | X ¯ i , k ) | . (6)

Here, K S is the Kolmogorov—Smirnov statistic, y is the model-estimated D-values and the variable X i is the i t h element of X = { T , p H , a w , C N , C P } .

The experimental datasets used for microbial inactivation modeling in this work are collated in Table 1. We chose datasets that are predominantly distinct in terms of microorganisms, media, process variables, and parameter space to demonstrate the tractability of our knowledge-informed data-driven pipeline for model development. We also found that the structure of the literature models developed from these datasets was under- or over-determined, rather than optimally determined. Consequently, the datasets and benchmark models we chose serve as an ideal testbed for evaluating the robustness of our approach.

Numerical codes were developed using MATLAB^® R2021a by adapting the prototype codes of SINDy provided in Brunton and Kutz (2019) and PAWN global sensitivity analysis given in Pianosi and Wagener (2015) and Pianosi et al. (2015).

Our data-driven modeling approach combines SINDy and global sensitivity analysis to identify model equations and key factors that govern microbial inactivation in food. As a main feature, users can define any functional forms of input variables ( T , p H , a w , C N , C P ) , i.e., g T ( T ) , g p H ( p H ) , g a w ( a w ) , g C N ( C N ) , g C P ( C P ), and output variables (i.e., g D ( D )). As explained below, we set g D ( D ) = log ⁡ D (instead of D ) and determined g T ( T ) = 1 / T based on the Arrhenius equation, while using first-order terms for the other input variables, i.e., g p H ( p H ) = p H , g a w ( a w ) = a w , g C N ( C N ) = C N , and g C P ( C P ) = C P . Subsequently, a library of input terms is generated through polynomial combinations of those basic input variables provided from the user. SINDy, then, identifies a sparse model by choosing a minimum number of input terms (included in the library) that is required to represent the output variable with an acceptable accuracy (cf. Section 2.1). The resulting equations derived by SINDy takes the form of a linear combination of nonlinear terms, and therefore, explicitly show the impacts of environmental variables on D-values. The impact of individual primitive input variables (not the combined terms) can be identified through PAWN global sensitivity analysis. The two complementary tools together identify key model equations and factors that govern D-value for a given pathogen or spoiler. The modeling workflow is illustrated in Figure 1. We term our approach knowledge-based data-driven modeling as we incorporate known insights of system characteristics (such as Arrhenius equation) as a key component to determine basic form of input and output variables as described in detail below.

The development of data-driven microbial inactivation model can be facilitated by known characteristics of the system. While first-order representation is a typical choice for input and output variables, it is possible to improve model performance by a more appropriate choice of their functional forms. For this purpose, we leverage mechanistic microbial growth models (Whiting, 1995) to inform our choice of functional forms for microbial inactivation dynamics as follows: d N d t = − k ( p ) N , (7) where N is population density and k ( > 0 ) is the deactivation rate constant, which is given as a function of a vector of environmental variables ( p ). If we maintain environmental variables constant over time, we can get the solution in an analytical form, i.e., N = N 0 ⁡ exp ( − k ( p ) t ) . (8)

By definition, the population density is N = 0.1 N 0 when t = D , i.e., 0.1 N 0 = N 0 ⁡ exp ( − k ( p ) D ) . (9)

Therefore, D-value is simply: D = − ln ( 0.1 ) k ( p ) . (10)

Subsequently, applying logarithm to the equation above yields: log ⁡ D = C − log [ k ( p ) ] , (11) where C is a constant. Given that many prior Arrhenius-based models produce reasonable fit to growth data by relating the growth rate to various environmental variables as ln ⁡ k = f ( 1 / T , p H , a w , ⋯ ) (Whiting, 1995; Ross and Dalgaard, 2003), we similarly re-write Eq. 11 as: log ⁡ D = C + f ( 1 T , p H , a w , C N , C P ) . (12)

Consequently, the functional forms of output variable and input variables provided for SINDy implementations are y = log ⁡ D and X = [ 1 / T , p H , a w , C N , C P ] in Θ ( X ) , respectively, in reference to the generalized form in Eq. 4.

To substantiate our choice of functional forms for input and output variables in the preceding section, we compare the model performance per our approach (orange lines in Figure 2) against another base case with non-logarithmic D-values and non-reciprocal temperature and other process variables (blue lines in Figure 2). Here, complete non-parsimonious models are used to ensure fair comparison of the models without the influence of sparse regression. Our approach consistently performed better in terms of MSE calculated based on logarithmic D-values for both cases across different datasets and orders of polynomial combinations of input variables. For all the results that follow henceforth, our choice of the functional forms (i.e., orange lines in Figure 2) are adopted.

With log ⁡ D chosen as the target variable, we identify the optimal model structure that balances both accuracy and sparsity. We first determine the order of polynomial combination without accounting for model sparsity (i.e., with λ = 0 ) and subsequently choose the appropriate value of λ (now to promote sparsity). This two-step process is demonstrated through the case study of the dataset from Cerf et al. (1996) (Figures 3, 4). In doing this, we used three major criteria including AIC values, MSE, and the number of terms. The analysis based on the first criterion suggested us to choose the third order polynomial combination (Figure 3A) and λ = 0 (Figure 3B) where the AIC scores are minimal. In contrast, determination of the order of polynomial combination is not clear based on MSE because it keeps decreasing as the order increases (Figure 4A), highlighting the utility of information-theoretic criterion. While the third-order model with λ = 0 may be a desirable choice from a rigorous statistical point of view, we found that the increase of MSE is not significant up to λ = 1.36 (Figure 4B) where the number of terms can be further reduced from 20 to 17 (Figure 4C). MSE was significantly increased when λ > 1.36 without significantly reducing the number of terms, leading us to choose the third-order model with λ = 1.36 .

We also applied this stepwise model construction approach to the datasets from Juneja et al. (1995) (Supplementary Figures S1, S2) and Villa-Rojas et al. (2013) (Supplementary Figures S3, S4). The analysis for the dataset from Juneja et al. (1995) showed that AIC values have two local minima at the polynomial orders 2 and 4 (Supplementary Figure S1A). Through further checking with MSE and the number of terms, we chose the second-order model for better interpretability. After determining the polynomial order, we subsequently determined the optimal value of λ to be 0.24 (Supplementary Figure S1B). The changes of MSE and the number of terms as polynomial orders and λ values support our choice (Supplementary Figure S2). Lastly, the analysis of the dataset from Villa-Rojas et al. (2013) suggested the first-order model (Supplementary Figures S3, S4), which is because any further increase of model complexity would result in severe overfitting due to limited data points ( n = 16 ). In this case, we have not further reduced model complexity by fine tuning λ . The final model was therefore a simple equation with two input variables ( T and a w ).

Following the guidelines outlined in the preceding and Methods sections, we developed models for all experimental datasets considered in this work. The resulting model equations are summarized in Table 2. The individual models consist of varying number and degree of monomial terms as identified through our data-driven model development pipeline which are optimal to represent the output variable within the parameter space of the respective datasets. The identification of the optimal terms (especially higher-order terms) would not be possible with the previous approaches that rely on empirical choices of equation terms. The issues of model overfitting and uncertainties are also minimized with our stepwise approach in model design.

We compare the performance of our models with existing models from the literature in Figure 5. Our models consistently perform better than the literature models across all datasets, particularly for the dataset from Juneja et al. (1995) that has considered additional process variables, i.e., C N and C P . It is certainly possible that C N and C P have significant interaction effects with other process variables and are critical to characterize microbial inactivation dynamics, which explains the enhanced accuracy that accompanies their inclusion in the model. Quantitative information of the models is tabulated in Table 3. In all cases, our models perform adequately with reasonable error measures, i.e., MSE ≤ O ( 10 − 2 ) and low AIC scores. In two of the cases (models for datasets from Juneja et al. (1995) and Villa-Rojas et al. (2013)), our models gained more than ten-fold increase in accuracy with fewer number of functional terms in the model structures as compared to the literature models. Moreover, we demonstrate the opportunity to further enhance the model accuracy over two-folds for dataset from Cerf et al. (1996) by considering a more complex model equation (greater number of functional terms) without the risk of overfitting as shown by the lower AIC score as compared to the literature model. By carefully adopting the stepwise model tuning scheme as described in Section 3.2, we were able to optimally tune the models to achieve better accuracy while minimizing the chances of overfitting as compared to existing literature models.

While our model governing equations offer good representations of experimental data, the individual effects of process variables remain elusive. In addition to data-driven modeling, we leverage global sensitivity analysis (Pianosi and Wagener, 2015) using the model-derived governing equations to identify key governing process variables and possibly divulge the interactions between them. Here, the sensitivities are evaluated for the entire parameter space encompassed by the respective datasets (Table 1). Our results demonstrate highly disparate sensitivity measures for the process variables across datasets as shown in Figure 6, positing that the relative effects of the process variables may be highly environment dependent. For example, T exhibited high sensitivity in the model for data from Cerf et al. (1996) but registered lower sensitivity than other process variables in the model for data from Juneja et al. (1995), although T is often viewed as the primary process variable in most microbial inactivation experiments. Conversely, a w consistently displayed relatively lower sensitivity measures than other process variables across the models. While condition-specificity is one of the plausible explanations for contrasting sensitivity measures, it should also be noted that the global sensitivity analysis strictly represents the individual influence of process variables on output and are indifferent to interactions effects between process variables. Therefore, there may be instances where the sensitivity measures of the individual process variables are insignificant, whilst considerable interactions effects with other process variables exist. For instance, T and pH registered low sensitivity in the model for data from Juneja et al. (1995) with statistically insignificant influence on the output variable (Figure 6), but the variables were repeatedly included in numerous terms in the governing equation identified through our knowledge-informed data driven approach (Table 2). This is expected given that SINDy retains the most influential terms irrespective of individual or interaction effects necessary for representation of output data.

To resolve this issue, we implement the combined use of data-driven approach and global sensitivity analysis that enables elucidation of the factors influencing the process dynamics from the contexts of governing equations, key process variables and critical model monomial terms. To this effect, we iteratively removed each term in the governing equations, and subsequently refitted the models with new MSEs as shown in Figure 7, whereby considerable increase in MSE (as compared to the original complete model, represented by horizontal dashed lines) indicates that the respective terms are essential to represent the output variable. In the model for data from Cerf et al. (1996), removal of terms containing p H and a w results in considerable increase in MSEs though the variables had relatively insignificant influences on the output in global sensitivity analysis. Clearly, p H and a w significantly influence the effects of T on microbial inactivation through strong interaction effects, but the reverse is not necessarily true. The finding here is reinforced by the outcome from SINDy, where the terms with p H and a w were unequivocally retained in the governing equations despite diminished main effects (individual effects), as the variables are still relevant to represent the output through their influence on T . Therefore, users can, hypothetically, fix the p H and a w at arbitrary optimal levels and tune only T as an alternative reduced-order experimental optimization.

Conversely, in the model for data from Juneja et al. (1995), removal of terms with C N and C P does not significantly increase the MSE despite their elevated global sensitivity. Therefore, the effects of C N and C P on microbial inactivation are possibly dependent on other process variables but they do not impose similar magnitudes of influence in return. Nevertheless, SINDy has retained the terms with C N and C P as their substantial individual effects are critical to represent the output when other process variables are invariant. For this system, a stepwise process optimization will work best where T and p H is independently tuned first followed by the optimization of C N and C P . Lastly, the model for data from Villa-Rojas et al. (2013) is fairly simple where both involved process variables possibly impose equivalent individual and interaction effects. While the global sensitivity analysis and reassessment of MSEs with iterative removal of terms offer additional insights that can aid process optimizations, we do not recommend the users to influence the model selection through these insights as they are model-derived contextual outcomes, and thus, the model should be thoroughly optimized beforehand, through the iterative feedback loop (Figure 1) as desired.