Holistic Approach to the Uncertainty in Shelf Life Prediction of Frozen Foods at Dynamic Cold Chain Conditions

Systematic kinetic modeling is required to predict frozen systems behavior in cold dynamic conditions. A one-step procedure, where all data are used simultaneously in a non-linear algorithm, is implemented to estimate the kinetic parameters of both primary and secondary models. Compared to the traditional two-step methodology, more precise estimates are obtained, and the calculated parameter uncertainty can be introduced in realistic shelf life predictions, as a tool for cold chain optimization. Additionally, significant variability of the real distribution/storage conditions is recorded, and must be also incorporated in a kinetic prediction scheme. The applicability of the approach is theoretically demonstrated in an analysis of data on frozen green peas Vitamin C content, for the calculation of joint confidence intervals of kinetic parameters. A stochastic algorithm is implemented, through a double Monte Carlo scheme incorporating the temperature variability during distribution, drawn from cold chain databases. Assuming a distribution scenario of 130 days in the cold chain, 93 ± 110 days remaining shelf life was predicted compared to 180 days assumed based on the use by date. Overall, through the theoretical case study investigated, the uncertainty of models’ parameters and cold chain dynamics were incorporated into shelf life assessment, leading to more realistic predictions.


Introduction
One of the principal goals in food science and industrial practice is to develop ways to improve quality of foods, by means of controlling chemical, physical and microbiological changes during processing and storage [1]. In most cases, food quality is gradually deteriorating, as food matrices are physicochemically and biologically active systems [2]. Increased interest is focused on deriving kinetic information for different food products, so as to calculate change rates, and thus be able to estimate their shelf life. In order to be able to reliably assess product quality during processing and storage, appropriate mathematical equations and kinetic parameters are necessary. Further to the mean estimate of such parameters, complementary statistical evaluation methodology should also be employed, to allow for reliable assessment of the uncertainty of calculations in the conditions of the cold chain [3].
When reviewing literature on frozen food quality deterioration, most available experimental data are provided at a reference temperature (e.g., −18 • C) [4] and single measurements at selected, limited time points (e.g., after 6 and 12 months of storage) for different quality parameters of provide information on the combinations of model parameter values, that are encompassed by the joint confidence ellipsoid. This would allow for the exemption of combinations that fall outside the ellipsoid that denotes the joint confidence boundary [31][32][33]. In recent literature, there are scarce studies that have estimated and plotted the joint confidence intervals of kinetic parameters [25,31,[34][35][36][37][38]. Such approach would be particularly appropriate in order to account for the real uncertainty of model parameters, and thus proceed to realistic and reliable shelf life estimations, using Monte Carlo simulation. Although Monte Carlo techniques have been applied in recent literature for the probabilistic assessment of stochastic variability and uncertainty associated with various quality attributes of different food systems [16,18,[26][27][28][29], the statistical interrelation of the values applied in the iterative simulation process has not been considered.
Besides overlooking model parameter uncertainty (best expressed by the ±95% joint confidence plots), another problematic practice in shelf life determinations and predictions refers to assumed fixed temperature conditions at different stages of the cold chain, despite the well-known deviations that occur. Therefore, it would be important to also account for the real variability of cold chain conditions. The objective of this work was to holistically approach shelf life calculations through a double Monte Carlo technique applied to frozen food kinetics data, taking into account both kinetic parameter uncertainty and temperature variability at all stages of the actual cold chain. This will be demonstrated via case study computations encompassing the effect of different scenarios of cold chain control. The ultimate goal of this holistic approach is to obtain realistic shelf life predictions, examine the main factors affecting shelf life calculations through a sensitivity analysis and provide a tool for effective management and optimization of frozen food cold chain.

Basic Principles
According to the two-steps kinetic approach to model food quality degradation, the most representative quality attributes are carefully selected, their change is systematically assessed over time under isothermal storage conditions and an appropriate mathematical equation (primary model) is applied in order to describe this change. Subsequently, a secondary model is chosen to depict the effect of the most important factors (e.g., temperature), on the rate of changes. An alternative approach involves the implementation of a 1-step analysis, through a non linear fitting algorithm that uses an integrated unique equation, deriving more precise estimations [39].
Primary models describe the rate of loss of one or more quality factors, such as a nutrient or a pigment or the rate of the production of an undesirable compound, such as an off-flavor or discoloration (Equation (1)): where the quality index (A), that may represent a chemical, physical, microbiological or sensory attribute, is chosen to most representatively describe the quality degradation of the food system studied. The constant k is the reaction rate constant and n is the apparent order of the respective reactions. In order to calculate the values of k and n of Equation (1), the change of the experimentally measured values of (A) should be appropriately fitted with time; a mathematical function, frequently called "the Quality function" (Q(A)) vs. time is then derived [40], that allows for quality quantification. This function incorporates the dependency of quality change on several intrinsic and extrinsic factors, through their effect on the reaction rate constant, k (k = f (C i , E j ), with C i describing factors related to food composition and E j representing different environmental factors [41,42].
As far as secondary models are concerned, the Arrhenius equation (Equation (2)) is mostly applied to describe temperature dependence of quality changes of frozen foods [25]. Nonetheless, there have been serious arguments highlighting its restrictions in describing frozen matrix behavior in the temperature range near the glass transition temperature [43][44][45]. In [46,47], the pros and cons of the alternative secondary models often used for frozen foods, namely Arrhenius and WLF, are presented in detail, underlying their cautious applicability at the lower boundary of the temperature range investigated, since their logarithmic functions render them sensitive to rates.
where k ref is the rate constant at the reference temperature T ref (K), R: the universal gas constant and E a : the activation energy (J/mol or cal/mol).
In the case of frozen foods, the time of experimental measurements for a full kinetic study is often times impracticably long. In order to expedite the experimental procedure, without jeopardizing the accuracy of the results, the two-step kinetic methodology is implemented via the Accelerated Shelf Life Testing (ASLT) approach [48,49]. The value of the Q(A t ) quality function at time t, defined by Equation (1) in the case of isothermal conditions, is calculated, with T(t) describing the change of temperature as a function of time: where k eff is the value of the rate of the quality loss reaction at the effective temperature T eff . The T eff term represents the constant temperature that results in the same quality value as the variable temperature function over the same time period, which equals t tot .
If the T(t) distribution is discretized in small time increments t i of constant temperature T i (with t i = t tot ), and applying the Arrhenius equation as the secondary model, then Equation (3) can be alternatively written as Equation (4): From Equation (4), the value of k eff can be estimated and subsequently, from the Arrhenius model, the effective temperature T eff can be calculated.

Model Development and Determination of the Uncertainty of Kinetic Parameters
As a case study to test the proposed methodology, data from literature on the shelf life of frozen green peas was used [50], where Vitamin C loss was systematically measured at five constant temperatures. A first-order reaction order (primary model) was established for this chemical reaction (Equation (5)) and the Arrhenius equation was applied to efficiently quantify the temperature effect (secondary model), (Equation (2)): where C 0 the initial Vitamin C concentration (mg/100 g of food), k VitC is the reaction rate of the Vitamin C oxidation at a fixed temperature and at Equation (2) T ref for frozen foods equals to −20 • C. If a one step analysis is performed, a single equation (Equation (6)) that integrates both primary and secondary models (Equations (2) and (5)) is used and a non-linear regression is applied (SYSTAT 8.0).
then, one can estimate the mean value, as well as the 95% Confidence Intervals (C.I.) of the kinetic parameters, namely E a (in kJ/mol) and k ref (at −20 • C, in days −1 ). In order to estimate the Confidence Intervals of Equation (6) one needs to calculate the first derivative of the dependent variables. A model is linear in a parameter if the derivative with respect to that parameter is not a function of that parameter [25]. In the studied case, our model is non ϑE a = f(E a )); then the parameter cannot be solved for directly, but only through a non linear regression. The estimation of such non symmetric confidence intervals is not straightforward and the extent of their asymmetry depends on the nonlinearity of the function and the number of data [38]. As stressed out in [30], it is essential to describe the methodology used for confidence intervals calculation, in order to be able to interpret their meaning. In this work, MATLAB ® was used and appropriate code was written, in order to derive the asymptotic standard errors (SE), calculated by Equation (7) based on the t-parameter for a confidence level (±SE ·t(1-0.5a),n) ( [1,51]).
where J is the Jacobian matrix (Equation (8)), which is calculated by the partial derivatives of the model output with respect to the model parameters estimated at each measurement point, p is the number of estimated parameters (p = 2, in this work), n is the number of observations, the superscript T denotes the matrix/vector transpose operator and SSE represents the sum of standard errors.
The 95% C.I. are calculated based on Equations (7) and (8) [38], and the relative MATLAB command nlparci([kref,ea],res,j), where the res term, representing the Residuals, is calculated using Equation (9) at each experimental point, i: The confidence level selected actually expresses the probability that the confidence interval produced will contain the true parameter value. For example, a 95% confidence interval covers 95% of the normal curve and the probability of obtaining a value outside of this area is less than 5%. Therefore, it is assumed that the Arrhenius parameters E a and k ref can be described by a normal distribution curve, rather than a single value. Implementing this assumption, the variability calculated by the model of Equation (6) is incorporated within calculations of the shelf life of frozen green peas.
The next step in this analysis of the 1-step non-linear regression on the isothermal data, is to investigate the covariance of the simultaneously estimated parameters, by constructing the, joint confidence regions, according to the expression [52] (Equation (10)): where SSE(θ) is the least sum of squared differences, at optimal parameter values, and F is the upper 1 − a quantile for an F-distribution with p and n − p degrees of freedom. All combinations of kinetic parameters with sum of squares less than or equal to the calculated SSE(θ) values will be inside the joint confidence region. The confidence regions obtained with this method (Equation (9), called likelihood confidence regions, can be disjoint and unbounded [32]. All of our calculations were based on the iterative method of [53], which is described in [38].

Determination of the Variability of Temperature Conditions in the Cold Chain
When addressing the issue of shelf life determination, most studies on frozen foods refer to an average temperature throughout product handling in the cold chain. This oversimplistic assumption, besides being false, could possibly lead to unrealistic predictions of quality change and remaining shelf life, at any point of the cold chain. Consequently, it is crucial to account for the real temperature history of the product in the cold chain (including fluctuations and abusive conditions), in order to be able to predict in a more accurate way quality changes, at any point of the cold chain (at storage or distribution). Based on field studies of the real handling of frozen foods in the cold chain [54][55][56], a variable temperature environment was recorded, which often included stages of abusive storage or transport/transfer conditions. In order to estimate in a realistic way the loss of quality parameters at each stage, and predict accurately the remaining shelf life at the end of the cold chain, it is important to fully account for the effect of temperature history by incorporating the fluctuations occurring in the real cold, post-processing chain, and the conditions' variability within the model prediction algorithm. A very useful tool that provides plenty of data concerning frozen foods handling, in all distinct stages of storage/transport/distribution is the FRISBEE project cold chain database (www.frisbee-project.eu/coldchaindb, [55]) 1841 records of production warehouse, 53 records of distribution warehouse, 636 records of retail display and 354 records of domestic freezers were retrieved and statistically analyzed, in order to incorporate temperature variability within shelf life predictions, at any point of the cold chain. In Figure 1, temperature distribution of the three main stages of frozen food handling is depicted, based on the aforementioned statistical treatment of the database records.
of our calculations were based on the iterative method of [53], which is described in [38].

Determination of the Variability of Temperature Conditions in the Cold Chain
When addressing the issue of shelf life determination, most studies on frozen foods refer to an average temperature throughout product handling in the cold chain. This oversimplistic assumption, besides being false, could possibly lead to unrealistic predictions of quality change and remaining shelf life, at any point of the cold chain. Consequently, it is crucial to account for the real temperature history of the product in the cold chain (including fluctuations and abusive conditions), in order to be able to predict in a more accurate way quality changes, at any point of the cold chain (at storage or distribution). Based on field studies of the real handling of frozen foods in the cold chain [54][55][56], a variable temperature environment was recorded, which often included stages of abusive storage or transport/transfer conditions. In order to estimate in a realistic way the loss of quality parameters at each stage, and predict accurately the remaining shelf life at the end of the cold chain, it is important to fully account for the effect of temperature history by incorporating the fluctuations occurring in the real cold, post-processing chain, and the conditions' variability within the model prediction algorithm. A very useful tool that provides plenty of data concerning frozen foods handling, in all distinct stages of storage/transport/distribution is the FRISBEE project cold chain database (www.frisbee-project.eu/coldchaindb, [55]) 1841 records of production warehouse, 53 records of distribution warehouse, 636 records of retail display and 354 records of domestic freezers were retrieved and statistically analyzed, in order to incorporate temperature variability within shelf life predictions, at any point of the cold chain. In Figure 1, temperature distribution of the three main stages of frozen food handling is depicted, based on the aforementioned statistical treatment of the database records.

Shelf Life Assessment and Uncertainty Determination
The next step in the proposed methodology involves a Monte Carlo scheme. applied through a FORTRAN algorithm. This iterative algorithm is applied based on the previous assumption that Ea and kref variability are effectively represented by a normal distribution [14,25,57]. At each iteration, a random number is generated through an appropriate FORTRAN routine function, and a value is

Shelf Life Assessment and Uncertainty Determination
The next step in the proposed methodology involves a Monte Carlo scheme. applied through a FORTRAN algorithm. This iterative algorithm is applied based on the previous assumption that E a and k ref variability are effectively represented by a normal distribution [14,25,57]. At each iteration, a random number is generated through an appropriate FORTRAN routine function, and a value is assigned to E a and k ref (independently the one from the other). The exact parameter value assigned is based on the discretization of the corresponding normal distribution curve (Figure 2a), and thus the corresponding value frequency. The construction of such Gaussian distributions is based on the estimate of the mean value, and the ±95% C.I. of the kinetic parameters of Equation (6). assigned to Ea and kref (independently the one from the other). The exact parameter value assigned is based on the discretization of the corresponding normal distribution curve (Figure 2a), and thus the corresponding value frequency. The construction of such Gaussian distributions is based on the estimate of the mean value, and the ±95% C.I. of the kinetic parameters of Equation (6). Going further with analyzing kinetic data, one cannot overlook the potential correlation between the kinetic parameters (Ea and kref), which means that the confidence interval of one parameter depends on the value of the other parameter. Therefore, joint confidence regions were derived using MATLAB and Equation (9); this information was used in order to exclude some of the pairs of values of Ea-kref, obtained by the random algorithm of Monte Carlo.
A final step for improving shelf life assessment is based on the incorporation of temperature variability, as depicted in Figure 1, within the Monte Carlo scheme. Applying this nested iterative algorithm by means of a FORTAN code, at each iteration, a double scheme is applied: (1) as already discussed, a random number is generated and a value is assigned to Ea and kref based on the discretization of the corresponding normal distribution curves (2) pairs of Ea and kref that do not fall within the estimated joint confidence intervals are excluded from further analysis and (3) a nested Monte Carlo algorithm is applied, where temperatures at each of the three stages are randomly selected based on the discretization of the distribution curves of Figure 1. With all parameter value assigned, Vitamin C retention is then calculated, based on Equation (6), and the shelf life can be accordingly estimated.

Application of the Holistic Approach to Shelf Life Prediction in the Frozen Green Peas Cold Chain
In order to test and validate the abovementioned methodology, results from an isothermal study (storage at 5 sub-zero temperatures) of frozen green peas [45] were used and a global 1-step procedure was implemented, using Equation (6). Based on 50% Vitamin C loss, green peas shelf life Going further with analyzing kinetic data, one cannot overlook the potential correlation between the kinetic parameters (E a and k ref ), which means that the confidence interval of one parameter depends on the value of the other parameter. Therefore, joint confidence regions were derived using MATLAB and Equation (9); this information was used in order to exclude some of the pairs of values of E a -k ref, obtained by the random algorithm of Monte Carlo.
A final step for improving shelf life assessment is based on the incorporation of temperature variability, as depicted in Figure 1, within the Monte Carlo scheme. Applying this nested iterative algorithm by means of a FORTAN code, at each iteration, a double scheme is applied: (1) as already discussed, a random number is generated and a value is assigned to E a and k ref based on the discretization of the corresponding normal distribution curves (2) pairs of E a and k ref that do not fall within the estimated joint confidence intervals are excluded from further analysis and (3) a nested Monte Carlo algorithm is applied, where temperatures at each of the three stages are randomly selected based on the discretization of the distribution curves of Figure 1. With all parameter value assigned, Vitamin C retention is then calculated, based on Equation (6), and the shelf life can be accordingly estimated.

Application of the Holistic Approach to Shelf Life Prediction in the Frozen Green Peas Cold Chain
In order to test and validate the abovementioned methodology, results from an isothermal study (storage at 5 sub-zero temperatures) of frozen green peas [45] were used and a global 1-step procedure was implemented, using Equation (6). Based on 50% Vitamin C loss, green peas shelf life can be predicted at any arbitrary reference temperature, using Equations (2)  namely E a = 104.24 ± 11.34 kJ/mol and k ref = 0.00177 ± 0.000494 days −1 , are slightly different from those estimated via the 2-step analysis. Using a one step analysis, frozen green peas shelf life is also estimated at −18 • C at approximately 250 days. It can be also observed that when applying the two-step approach, the 95% C.I., (calculated via regression analysis) are usually wider than those calculated with a global-one step approach [22].
As discussed earlier, the next step involves the application of a Monte Carlo simulation scheme, assuming that E a and k ref variabilities are described by a normal distribution (Figure 2a). The Gaussian distributions illustrated in Figure 2a were constructed based on the estimate of the mean value, (for example E a = 104.24 kJ/mol) and its standard deviation, (in the case of E a , σ = 5.8 kJ/mol). The same procedure was followed to construct the corresponding distribution curve for k ref .
The Monte Carlo algorithm aimed at ascribing specific values to the E a and k ref parameters, based on their normal distribution curve, in order to estimate the shelf life (Equation (6)) at an arbitrarily chosen temperature of −18 • C, using the 50% Vitamin C loss as the acceptability limit. Results for Shelf Life (SL) calculation, including its uncertainty, are depicted in Figure 2b, expressed as a frequency curve with a mean value (SL estimate) ± 95% C.I., equal to 254.2 ± 29.9 days, giving a more realistic prediction than the single value estimation of 250 days, based on Arrhenius parameters' mean estimated values. The necessary number of Monte Carlo simulations for a specific kinetic model is not clearly defined in literature (ranging from some hundreds to few thousands), since there are other factors, such as the degrees of freedom, the range of inputs, possible parameter interactions, etc., that affect algorithm performance [58]. However, the number of 10 4 (applied in this study) seems to be a frequently used number of iterations [59,60] and an acceptable compromise between computing power/time and result accuracy.
The aforementioned procedure is in agreement with the approach described in [29], where the confidence intervals of E a and k ref parameters of the Arrhenius equation were estimated for degradation of cyanidins under dynamic conditions, by generating artificial data of the initial measurements, superimposing the experimental error. The use of a Gaussian distribution is a common practice in food engineering [18,61]; however, the same methodology can be implemented in case a different probability distribution describes better data variability [28,[62][63][64][65][66].
As a next step, the joint confidence intervals estimated by MATLAB and Equation (9)   In order to implement this methodology in the real cold chain, and make an attempt to expand it by introducing the recorded temperature variability, a 130-days distribution scenario of frozen green peas is assumed, including 60 days at the warehouse (Figure 1a), 40 days at the retail level ( Figure 1b) and 30 days at the domestic freezer (Figure 1c). The main goal is to predict the remaining shelf life of such a product prior to consumption, taking into account not only the estimated In order to implement this methodology in the real cold chain, and make an attempt to expand it by introducing the recorded temperature variability, a 130-days distribution scenario of frozen green peas is assumed, including 60 days at the warehouse (Figure 1a), 40 days at the retail level ( Figure 1b) and 30 days at the domestic freezer (Figure 1c). The main goal is to predict the remaining shelf life of such a product prior to consumption, taking into account not only the estimated uncertainty of kinetic parameters, but also temperature variability at each stage of the cold chain (Figure 1). When only average temperatures are introduced in Equation (6), and mean estimates of kinetic parameters E a -k ref are considered, the remaining Shelf Life after 130 days of handling is estimated at 106 days (at a 'reference' temperature of −18 • C), (black line, Figure 4a). Black line represents the SL distribution without considering Ea-kref correlation, blue line after excluding the Ea-kref pairs of Figure 3a.
In order to implement this methodology in the real cold chain, and make an attempt to expand it by introducing the recorded temperature variability, a 130-days distribution scenario of frozen green peas is assumed, including 60 days at the warehouse (Figure 1a), 40 days at the retail level ( Figure 1b) and 30 days at the domestic freezer (Figure 1c). The main goal is to predict the remaining shelf life of such a product prior to consumption, taking into account not only the estimated uncertainty of kinetic parameters, but also temperature variability at each stage of the cold chain ( Figure 1). When only average temperatures are introduced in Equation (6), and mean estimates of kinetic parameters Ea-kref are considered, the remaining Shelf Life after 130 days of handling is estimated at 106 days (at a 'reference' temperature of −18 °C), (black line, Figure 4a).  In Figure 4a, the degradation of Vitamin C is depicted vs. time, for different handling scenarios, throughout the frozen green peas cold chain, based on temperature data retrieved from Figure 1a-c. Applying a Monte Carlo iterative scheme, 10 4 different temperature cases (3 of the random scenarios are depicted in Figure 4a, using different colours) were ran and the remaining shelf life at the end of the 130 days cycle was calculated (Figure 4b). In this approach, it is important to bear in mind that  In Figure 4a, the degradation of Vitamin C is depicted vs. time, for different handling scenarios, throughout the frozen green peas cold chain, based on temperature data retrieved from Figure 1a-c. Applying a Monte Carlo iterative scheme, 10 4 different temperature cases (3 of the random scenarios are depicted in Figure 4a, using different colours) were ran and the remaining shelf life at the end of the 130 days cycle was calculated (Figure 4b). In this approach, it is important to bear in mind that kinetic parameters E a -k ref in Equation (6) were assumed constant, with fixed values equal to the mean estimates of regression analysis (Figure 2a).
Finally, in a further attempt to improve quality change predictions, and therefore obtain more realistic estimations of the remaining shelf life, kinetic parameter uncertainty was also incorporated in a double Monte Carlo scheme throughout the 130 days cycle (Figure 5), and thus all sources of 'error' (parameter uncertainty and temperature variability) were integrated in a holistic approach of kinetic data analysis. In this case, temperature variability is taken into account (by a Monte Carlo iterative algorithm using the temperature distributions of Figure 1a-c), and a second, 'nested' Monte Carlo scheme is used during each of the three distinct stages where a random pair of E a -k ref is selected considering not merely the distributions depicted in Figure 2a, but also the intercorrelation between these parameters, as depicted in the Joint Confidence Intervals (Figure 3a). Having incorporated within the model both temperature variability and parameter uncertainty, remaining shelf life predictions become significantly broader, as it is clearly depicted in Figure 6, where results deriving from different approaches are comparatively depicted. When all possible sources of error are incorporated, remaining shelf life is estimated at 93.4 ± 110.5 (days), vs. 106.0 ± 71.3 (days) when only temperature variability is considered, vs. 112.4 ± 26.0 when temperatures are assumed constant during the three stages (and equal to the mean estimate of each distribution of Figure 1), accounting only for parameter uncertainty.
All these predictions can also be compared to the mean value predicted (106 days) without considering any source of error.
Carlo scheme is used during each of the three distinct stages where a random pair of Ea-kref is selected considering not merely the distributions depicted in Figure 2a, but also the intercorrelation between these parameters, as depicted in the Joint Confidence Intervals (Figure 3a). Having incorporated within the model both temperature variability and parameter uncertainty, remaining shelf life predictions become significantly broader, as it is clearly depicted in Figure 6, where results deriving from different approaches are comparatively depicted. When all possible sources of error are incorporated, remaining shelf life is estimated at 93.4 ± 110.5 (days), vs. 106.0 ± 71.3 (days) when only temperature variability is considered, vs. 112.4 ± 26.0 when temperatures are assumed constant during the three stages (and equal to the mean estimate of each distribution of Figure 1), accounting only for parameter uncertainty. All these predictions can also be compared to the mean value predicted (106 days) without considering any source of error.

Effect of Parameter Uncertainty
In this case, it can also be observed that the prevailing source of error is mainly attributed to the temperature variability, and less to the kinetic parameter uncertainty; this could be justified by the relatively narrow ±95% confidence intervals of the kinetic parameters, Ea-kref, as obtained by the 1step non linear regression, based on the Vitamin C measurements. If an artificial error is introduced in raw data, and the derived ±95% confidence intervals (using MATLAB) are assumed broader (as depicted in the joint confidence region within Figure 7a), it can be observed that the contribution of

Effect of Parameter Uncertainty
In this case, it can also be observed that the prevailing source of error is mainly attributed to the temperature variability, and less to the kinetic parameter uncertainty; this could be justified by the relatively narrow ±95% confidence intervals of the kinetic parameters, E a -k ref , as obtained by the 1-step non linear regression, based on the Vitamin C measurements. If an artificial error is introduced in raw data, and the derived ±95% confidence intervals (using MATLAB) are assumed broader (as depicted in the joint confidence region within Figure 7a), it can be observed that the contribution of parameter uncertainty becomes more pronounced, than that depicted in Figure 6. Similarly, if initial measurements are improved so as to provide much narrower ±95% confidence intervals (as depicted in the joint confidence region within Figure 7b), one could conclude that parameter uncertainty introduces a relatively small error in comparison to the effect of the storage temperature variability that could practically be neglected.

Effect of Temperature Variability
Based on our results, temperature conditions throughout the current cold chain, deviating from the ideal ones, have an important impact on overall uncertainty, when assessing the quality status at any point of distribution. A 'what if' analysis was also performed regarding temperature variability (Figure 8), assuming a narrower distribution (red, narrower distribution within both Figure 8a,b), during storage at the domestic freezer (3rd stage), which is found to be the weakest link of the chain. As it can be seen in Figure 8, the remaining shelf life is estimated to be 103.9 ± 74.8 days (Figure 8a, narrower temperature conditions during home storage), and 113.1 ± 65.6 days (Figure 8b, assuming reference conditions during home storage), a quite narrower, improved distribution compared to the 93.4 ± 110.5 days, derived when temperatures during the third stage followed the broad, often abusive distribution of Figure 1c (based on real data).

Effect of Temperature Variability
Based on our results, temperature conditions throughout the current cold chain, deviating from the ideal ones, have an important impact on overall uncertainty, when assessing the quality status at any point of distribution. A 'what if' analysis was also performed regarding temperature variability (Figure 8), assuming a narrower distribution (red, narrower distribution within both Figure 8a,b), during storage at the domestic freezer (3rd stage), which is found to be the weakest link of the chain. As it can be seen in Figure 8, the remaining shelf life is estimated to be 103.9 ± 74.8 days (Figure 8a, narrower temperature conditions during home storage), and 113.1 ± 65.6 days (Figure 8b, assuming reference conditions during home storage), a quite narrower, improved distribution compared to the 93.4 ± 110.5 days, derived when temperatures during the third stage followed the broad, often abusive distribution of Figure 1c (based on real data).
during storage at the domestic freezer (3rd stage), which is found to be the weakest link of the chain. As it can be seen in Figure 8, the remaining shelf life is estimated to be 103.9 ± 74.8 days (Figure 8a, narrower temperature conditions during home storage), and 113.1 ± 65.6 days (Figure 8b, assuming reference conditions during home storage), a quite narrower, improved distribution compared to the 93.4 ± 110.5 days, derived when temperatures during the third stage followed the broad, often abusive distribution of Figure 1c (based on real data).

Conclusions
When revisiting recent literature on kinetics of food-related reactions, limited number of investigations have studied the significance of statistical aspects, such as the effect of kinetic parameter uncertainty. Additionally, few researches take into account quantitatively the significant temperature variability, recorded in all stages of the actual cold chain of frozen foods. Most physicochemical properties of raw food materials, even of the same species or cultivar, are found to assume values in a wide range, a variability that is difficult to be quantified and properly introduced within kinetic models. Therefore, next to the inconsistency of the initial quality, it is important to investigate in depth the significance and the impact of parameters uncertainty (as expressed by their 95% Confidence Intervals), as well as temperature conditions' variability, especially when the goal is to derive realistic model predictions. Reporting of parameter uncertainty is essential when kinetic results are presented. The use of computational tools such as Monte Carlo technique can provide a practical approach in providing an integrated picture of the impact of both experimental and statistical uncertainties as well as the potential conditions' dynamics.
In this work, by applying the proposed technique with the double Monte Carlo scheme on Vitamin C degradation of frozen green peas, a broader, more realistic prediction of the remaining shelf life was obtained based on the fact that the model incorporates more realistic conditions, as model validation was not implemented in this study. The aforementioned methodology allows also for a sensitivity analysis, revealing the importance of reducing the 95% Confidence Intervals of the kinetic parameters (possibly by optimizing the experimental design), and/or by improving temperature conditions of the cold chain, actions that would lead to narrower distributions of the remaining shelf life, with less unacceptable products.
As a final remark, it should be pointed out that, although the preceding analysis was exemplified by application on a first-order reaction (Vitamin C loss) and temperature dependence modelled by the Arrhenius equation, the same approach, using identical steps of methodology, can be implemented when an alternative primary or secondary model is used to describe kinetic data.
Author Contributions: M.G. and P.T. contributed equally to this research article. Conceptualization M.G. and P.T.; M.G. wrote the FORTRAN and the MATLAB code and run the simulations; both authors contributed in designing the structure, analyzing the data and writing the paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.