Development of a risk assessment model to predict the occurrence of late blowing defect in Gouda cheese and evaluate potential intervention strategies

Late blowing defect (LBD) is an important spoilage issue in semi-hard cheese, with the outgrowth of Clostridium tyrobutyricum spores during cheese aging considered to be the primary cause. Although previous studies have explored the microbial and physicochemical factors influencing the defect, a risk assessment tool that allows for improved and rational management of LBD is lacking. The purpose of this study was to develop a predictive model to estimate the probability of LBD in Gouda cheese and evaluate different intervention strategies. The spore concentration distribution of butyric acid bacteria (BAB) in bulk tank milk was obtained from 8 dairy farms over 12 mo. The concentration of C. tyrobutyricum from raw milk to the end of aging was simulated based on Gouda brined for 2 d in saturated brine at 8°C and aged at 13°C. Predicted C. tyrobutyricum concentrations during aging and estimated concentration thresholds in cheese at onset of LBD were used to predict product loss due to LBD during a simulated 1-yr production. With the estimated concentration thresholds in cheese ranging from 4.36 to 4.46 log most probable number (MPN)/ kg of cheese, the model predicted that 9.2% (±1.7%) of Gouda cheese showed LBD by d 60; cheeses predicted to show LBD at d 60 showed a mean pH of 5.39 and were produced with raw milk with a mean BAB spore count of 143 MPN/L. By d 90, 36.1% (±3.4%) of cheeses were predicted to show LBD, indicating that LBD typically manifests between d 60 and 90, which is consistent with observations from the literature and the cheese industry. Sensitivity analysis indicated that C. tyrobutyricum maximum growth rate as well as concentration threshold in cheese at onset of LBD are the most important variables, identifying key data needs for development of more accurate models. The implementation of microfiltration or bactofugation of raw milk (assumed to show 98% efficiency of spore removal) in our model prevented occurrence of LBD during the first 60 d of aging. Overall, our findings provide a framework for predicting the occurrence of LBD in Gouda as well as other cheeses and illustrate the value of developing digital tools for managing dairy product quality.


INTRODUCTION
Late blowing defect (LBD) is a major spoilage issue in semi-hard cheeses and is characterized by formation of cracks and slits in the cheese matrix as well as unpleasant rancid flavor due to unwanted butyric acid fermentation that produces gas and butyric acid. An anaerobic gram-positive spore-former, Clostridium tyrobutyricum is capable of butyric acid fermentation and is able to grow and produce gas at typical aging temperatures and salt concentrations used in production of cheeses such as Gouda and Emmental (Su and Ingham, 2000). These brined cheeses are particularly susceptible to LBD because the salt can take time to migrate to the interior of cheese matrix and is therefore unable to limit the growth of C. tyrobutyricum at initial phases of aging. For these reasons and due to the fact that it is more salt tolerant than other butyric acid bacteria (BAB), C. tyrobutyricum has been presented as the primary cause of semi-hard cheese LBD in multiple studies and book chapters (Klijn et al., 1995;Le Bourhis et al., 2007;Garde et al., 2013;D'Amico, 2014;Düsterhöft et al., 2017). Although some other clostridial species, such as Clostridium sporogenes and Clostridium beijerinckii, have also been shown to contribute to LBD (Le Bourhis et al., 2007;Morandi et al., 2021), C. tyrobutyricum is the Clostridium sp. most commonly isolated from cheese with LBD. Several different studies have been conducted that aimed to understand the effects of various factors on the occurrence of cheese LBD. These factors include cheese brining and aging conditions (Su and Ingham, 2000), the use of nisin-producing lactic acid bacteria (Toyoda and Nakajima, 1995;Rilla et al., 2003;Ávila et al., 2020), and C. tyrobutyricum concentration threshold in cheese at onset of LBD (Herman et al., 1997). Importantly, studies indicated that C. tyrobutyricum levels in cheese with LBD are difficult to define and can vary substantially depending on the characteristics of the cheese (e.g., shape, size, pH, structure; D'Amico, 2014) as well as the methods used to quantify C. tyrobutyricum, which have included quantifying spores (Stadhouders et al., 1985;Le Bourhis et al., 2007;Garde et al., 2011), viable cells (Gibbs and Freame, 1965), or metabolic products such as butyric acid (Bogovič Matijašić et al., 2007;Le Bourhis et al., 2007;Garde et al., 2011). Results from the few studies that have reported the concentration of viable C. tyrobutyricum cells at occurrence of LBD, however, suggest substantial variability and uncertainty with regard to the C. tyrobutyricum threshold associated with occurrence of LBD. For example, Herman et al. (1997) reported that a cell count of 0.4 most probable number (MPN)/g in Gouda cheese was sufficient to cause LBD, whereas Toyoda and Nakajima (1995) suggested, based on their growth data, that LBD did not occur until the cell concentration of C. tyrobutyricum exceeded 3.5 × 10 3 cfu/g of Gouda cheese. The high degree of variability and uncertainty poses a challenge in establishing an accurate C. tyrobutyricum concentration threshold in cheese at onset of LBD. As part of our current study, we calculated theoretical concentration thresholds as well as estimated the uncertainty of thresholds in specific cheese (i.e., Gouda) with defined size, pH, and aging temperature using a well-established growth model. Although previous studies (Toyoda and Nakajima, 1995;Herman et al., 1997;Su and Ingham, 2000;Rilla et al., 2003;Ávila et al., 2020) have made important contributions to investigating individual factors that facilitate the development of LBD, the overall high variability in factors such as concentration of C. tyrobutyricum spores in raw milk, cheese types, and cheese aging conditions requires a more sophisticated tool to account for all factors.
As a stochastic model that can simultaneously consider all possible combinations of model variables represented in a food system, Monte Carlo simulations can predict all possible outcomes associated with a defined problem (McNab, 1998). This modeling approach has been widely applied to predict spoilage in dairy products (Buehler et al., 2018a,b) and food safety risk in other food products (Pouillot et al., 2007;Pang et al., 2017;Zoellner et al., 2019). By simulating the growth, inactivation, transfer, or partition of microorganisms of interest over food processing or storage, this model can estimate the microbial concentrations in final food products and predict the associated risks. Compared with deterministic models that yield fixed prediction results, Monte Carlo simulation can predict a model outcome with some uncertainty measure, and thus is more appropriate for complex food systems, which exhibit high uncertainty. In addition, researchers can perform sensitivity analysis and what-if scenarios within this modeling framework (McNab, 1998), to evaluate the relative importance of different model variables and effectiveness of intervention strategies, respectively. All these characteristics make Monte Carlo simulations suitable for the purpose of risk management. Within the context of our study, a risk management tool will be beneficial for cheese manufacturers to predict and mitigate LBD in Gouda cheese. Because this defect typically occurs between 1 and 5 mo of cheese aging (Düsterhöft and Van Den Berg, 2007;D'Amico, 2014;Düsterhöft et al., 2017), cheese manufacturers can avoid significant economic loss by using this tool to make better decisions and avoid product loss after investing time and resources into the long aging process.

Model Overviews
The model developed here consists of 3 modules, including (i) a module that simulates transfer of spores from raw milk to cheese before the aging stage, (ii) a module that predicts the growth of C. tyrobutyricum during cheese aging, and (iii) a module that predicts the occurrence of LBD at a given aging time based on C. tyrobutyricum concentration in cheese ( Figure 1).
Module (i). Data on BAB spore concentrations in bulk tank raw milk were collected monthly from 8 dairy farms in New York State over a 12-mo period (S. Xin, C. Qian, S. I. Murphy, M. Wiedmann, and N. H. Martin, Cornell University, Ithaca, NY, unpublished data). Due to missing data for one farm (where data were available for only 2 mo), a total of 86 BAB spore concentration data points were available; spore counts below the detection limit (18 MPN/L; 11 samples) were considered as 4.5 MPN/L. Because neither a log normal, Poisson, nor a negative binomial distribution fit the raw data, a box-cox transformation was applied to identify the lambda value that would yield a normal distribution of the data; transformation of raw data using a power of 0.1 resulted in a well-fitted normal distribution with a mean of 1.53 (MPN/mL) 0.1 and a standard deviation of 0.2 (MPN/mL) 0.1 . This distribution was used to simulate the starting spore concentration for 10,000 bulk tanks of 900 L of raw milk; this simulation was performed by drawing a random variable from this normal distribution, followed by a reverse transformation (power of 10). For each bulk tank, we simulated a batch of curd, processed into 20 wheels of 4.5 kg (10 lbs) of Gouda cheese. We assumed that 60% of spores in the milk are retained in the cheese matrix and 40% are removed with the whey; this retention factor corresponds to a 6-fold higher spore concentration in cheese compared with the raw milk (assuming a 10% cheese yield; Morandi et al., 2021). The model subsequently randomly distributed all spores found in the curd produced from each simulated 900-L batch into 20 cheese wheels, using a Poisson distribution to account for variability between wheels (see Table 1 for details). Our model assumed salting at 8°C for 2 d; we assumed no germination or outgrowth of C. tyrobutyricum spores during this time. Germination of all spores present in a given cheese wheel was assumed to occur after salting and immediately at the start of aging at 13°C.
Overall, module (i) assumed that raw milk is the only contamination pathway through which BAB enter the cheese production system. Furthermore, we assumed that all BAB that are found in raw milk represent C. tyrobutyricum; this assumption was made due to (1) the lack of quantification of BAB count data at the species level and (2) lack of data on growth parameters for other BAB species.  For example, if n = 1.5 for 1 bulk tank milk, the initial spore count will have N 0 = 1.5 10 = 58 MPN/L (MPN = most probable number).  Spore concentration increased by 6-fold due to milk coagulation. 5 Spore concentration in cheese immediately after salting. The variability between cheese wheels is accounted for using Poisson distribution.

Module (ii).
The purpose of the second module is to predict the concentration of C. tyrobutyricum throughout Gouda cheese aging. Because insufficient data were available in the literature to estimate the duration of C. tyrobutyricum lag phase, a Baranyi growth model without lag phase (Equation [1]) was used to predict the concentration of C. tyrobutyricum in each simulated cheese wheel (Baranyi and Roberts, 1994): In this equation, μ is the bacterial growth rate (per hour) during aging; e is a mathematical constant (i.e., Euler's number); N (log cfu/kg) is the predicted microbial concentration in aged cheese at a given aging time t (hour); N 0 is the initial spore concentration in raw milk (MPN/L); N max (log cfu/kg) is the maximum population that C. tyrobutyricum vegetative cells can reach in Gouda cheese. A value of 4 × 10 7 cfu/kg was used as N max ; this is the concentration of viable C. tyrobutyricum cells that was reported to be present in 4-mo-old Gouda cheese with LBD (Toyoda and Nakajima, 1995). The maximum growth rate of C. tyrobutyricum in 200,000 simulated cheese wheels was calculated by adjusting the maximum growth rate under optimum conditions (μ max ) to the specific cheese aging conditions in a given wheel (Equation [2]); μ max (per hour) of C. tyrobutyricum in cheese was reverse-calculated using previously reported C. tyrobutyricum exponentialphase growth data in Gouda cheese collected over 2 to 3 mo (Toyoda and Nakajima, 1995). The gamma terms (Equations [3]-[5]; Zwietering et al., 1996) indicate the inhibitory effect by individual factors due to deviation from the optimal conditions and were determined based on the reported growth conditions (temperature of 15°C, pH 5.69, and water activity of 0.96) and growth characteristics (T min , T opt , pH min , pH max , pH opt , a w min ; Table 2):    Minimum growth water activity was calculated from 3.5% salt-in-cheese. 3 Cheese water activity is calculated from a standard 2.
0% salt-in-cheese for Gouda and assumed to be constant in the calculation of growth rate. [5] T aging , pH aging , and a w aging are, respectively, the temperature, pH, and water activity during aging. Given these conditions, the maximum growth rate under optimum conditions was determined as 0.76/h. Intrinsic and extrinsic physiochemical factors for each simulated cheese were selected based on the following assumptions. A typical aging temperature (13°C) for Gouda cheese was used for all cheese wheels (McSweeney, 2007;Wemmenhove et al., 2016;Düsterhöft et al., 2017). Each Gouda cheese wheel was assigned a pH value from a triangular distribution with a minimum of 5.1, maximum of 5.63, and mode of 5.40. This distribution of pH values was based on reported pH of young Gouda aged for less than 3 mo (Jo et al., 2018).
Although the pH of Gouda cheese can increase during aging due to proteolysis, this change was reported to not be more than 0.15 (Düsterhöft et al., 2017), and thus we assumed, for simplicity, no change of pH over time in our model. Water activity is heterogeneous across the cheese matrix and consistently decreases until the third week of aging (Wemmenhove et al., 2016), and sufficient data were not available to include water activity dynamics in the model at this point. Therefore, we assigned the typical water activity at d 60 for growth rate calculation to all cheese wheels across vats. This a w value was determined by a 2-step calculation. First, a standard 2% salt-in-cheese (i.e., 2% wt/wt) typical to young Gouda was converted to the water-phase salt as described by Guinee and Fox (2017). Then the water activity was determined from this calculated water-phase salt following Equation [6] (Wemmenhove et al., 2016): a w = 0.995 − 0.00721·water-phase salt. [6] With the calculated maximum growth rate (μ max ) for each cheese wheel and a fixed N max , we predicted the concentration of C. tyrobutyricum for each cheese wheel at any time point during the aging using Baranyi growth model without lag phase (Equation [1]). We specifically calculated the concentration of C. tyrobutyricum at d 30,60,90, and 120, as cheese manufacturers may want to perform monthly quality checks. Module (iii). The function of this module is to predict the proportions of cheese with LBD at a given aging time. Any simulated cheese wheels with a concentration of C. tyrobutyricum over the concentration threshold was considered as showing LBD. The concentration threshold in cheese was calculated using previously reported C. tyrobutyricum concentration data for experimental Valtellina Casera cheese with LBD (Morandi et al., 2021). Because no available data provide the precise time of appearance of LBD in naturally contaminated Gouda cheese, data for Valtellina Casera cheese were used instead; this cheese variety has moisture, size, aging regimen, and susceptibility to late blowing comparable to those of Gouda. Because LBD gradually occurred between d 60 and 70 for experimental cheese with various starting spore concentrations (Morandi et al., 2021), we assumed that cheese wheels with the lowest initial spore concentration after salting (≈1,300 MPN/kg) developed LBD latest (d 70), and that the cheese wheels with the highest spore concentration after salting (≈2,500 MPN/kg) developed LBD earliest (d 60). Using these 2 combinations of initial spore counts and aging times, we calculated the expected C. tyrobutyricum levels at d 60 and 70 (Equation [1]), using the growth conditions (temperature 13°C, pH 5.0, and water activity 0.968) reported for these cheeses ( Table 3). The calculated concentrations represent reasonable estimates for C. tyrobutyricum concentration thresholds in cheese at onset of LBD and were used as lower-bound and higher-bound thresholds, to account for uncertainty. Upper-and lower-bound concentration thresholds were used separately to predict whether a given cheese wheel showed LBD, using the C. tyrobutyricum concentrations obtained from module (ii).  The water activity is calculated from the mean salt-in-cheese between rind and core at d 60.
3 Clostridium tyrobutyricum concentration thresholds at onset of late blowing defect.

Sensitivity Analysis
Local sensitivity analysis was performed on 4 model parameters, which were (1) maximum growth rate (μ max ), (2) concentration factor, (3) C. tyrobutyricum concentration threshold at onset of LBD, and (4) cheese pH, to investigate their relative importance. Each parameter was independently adjusted to 4 biologically reasonable values, and the resulting proportional increase or decrease in cheese with LBD compared with the base model was recorded. Two best cases were defined as any changes in parameters that led to less occurrence of LBD, and 2 worst cases referred to the parameter changes that increased the percentage of late blown cheese. For μ max , the best cases were 40% and 20% decrease, whereas the worst cases were 20% and 40% increase in μ max . For the concentration factor, the best cases were −2 and −1, whereas the worst cases were +1 and +2.
For the concentration thresholds in cheese, both upper and lower thresholds were adjusted separately by −0.5 log and −0.25 log for the worst cases, and +0.25 log and +0.5 log for the best cases. For cheese pH, the best cases were determined by decreasing all cheese pH by 0.1 and 0.2, and the worst cases were calculated by increasing all cheese pH by 0.1 and 0.2.

What-If Scenarios
A total of 6 scenarios were proposed to evaluate the efficacy of different intervention strategies (Table 4), including (1) removal of spores in raw milk via microfiltration or bactofugation (Düsterhöft et al., 2017), (2) addition of sodium nitrate at 2.5 g per 100 L of milk (Düsterhöft et al., 2017), (3) addition of lysozyme at 2.5 g per 100 L of milk (Stadhouder, 1990b), (4) inoculation with bacteriocinogenic lactic acid bacteria to milk at 0.3% (Toyoda and Nakajima, 1995), (5) utilization of premium raw milk with reduced BAB levels, and (6) reduction in aging temperature to 12°C. Specific implementation of each intervention in the model is detailed in Table 4.

Model Validation
Model outcomes were compared with published data reported by Düsterhöft and Van Den Berg, (2007), D'Amico (2014), and Düsterhöft et al. (2017). Specifically, Düsterhöft et al. (2017) reported that late blowing caused by growth of C. tyrobutyricum typically occurred in 1-to 5-mo-old cheese. Similarly, D'Amico (2014) observed that LBD usually appears after 1 to 4 mo of aging. In Gouda-type cheeses, Düsterhöft and Van Den Berg (2007) reported that LBD typically occurred after 4 to 6 weeks, which the author referred to as prolonged aging. These findings confirmed that LBD can start as early as 1 mo and can last up to 5 mo of aging. The reported time frame for LBD was used to qualitatively validate our model result.
In addition, information obtained from one Gouda cheese producer (who requested anonymity) indicates that LBD in Gouda cheese typically occurs between 6 wk and 3 mo of aging, with 6 wk being the most common time at which LBD occurs when cheese is aged at temperatures between 10 and 12°C; however, this producer indicated that defect occurrence can be delayed beyond 3 mo if cheese is aged or stored at temperatures below 10°C. They estimated that aging Gouda at 10 to 12°C resulted in around 30% product loss due to LBD, whereas this loss was reduced to around 10 to 15% when cheese was aged at 8.5°C. Information from the same Gouda producer also indicated substantial reduction of occurrence of LBD when a combination of different protective cultures and lower aging temperature (8.5°C) was used in Gouda production. Limited availability of specific data limited our ability to establish information on the occurrence of LBD across different Gouda cheese producers.

Model Programming and Data Availability
The model was developed in the R programming language (R Core Team, 2019). The seed for Monte Carlo simulation was set to 1 for consistency. The raw data for initial spore concentration used in the second modules and R code for entire study can be accessed at https: / / github .com/ FSL -MQIP/ MC -Cheese -blowing.

Base Model
The base model reported in this study was designed to predict the occurrence of LBD in Gouda cheese (Figure 1). This model was used to simulate production of 10,000 cheese vats of 900 L each, with each cheese vat simulated to yield 20 4.5-kg (10-lb) wheels of Gouda. In practical terms, this could represent cheese produced over a year by approximately 30 small cheese makers (each producing about 330 cheese vats a year) with each sourcing their raw milk from a single farm. The model was designed by combining 3 consecutive modules. The first 2 modules tracked allocation of C. tyrobutyricum in Gouda cheese from raw milk to final cheese, whereas the third module determined the proportions of cheese developing LBD based on the C. 3 spore/mL Gouda cheese made from milk with lysozyme (2.5 g/100 L of milk) is safe from butyric acid fermentation provided that spore level in raw milk <0.3 per mL (Stadhouder, 1990b The lactic acid bacteria (LAB) in the original experiment was a nisin-producing starter, Lactococcus lactis ssp. lactis ATCC11454. It was inoculated at 0.3% and cultured at 30°C for 24 h.
tyrobutyricum concentration in cheese (Figure 1), using both upper and lower C. tyrobutyricum concentration thresholds in cheese expected at onset of LBD. Based on module (i), the initial concentrations of BAB spores in the simulated 10,000 raw milk vats ranged from 0 to 4,250 MPN/L with a median concentration of 69 MPN/L (equivalent to 1.84 log MPN/L); hence, all raw milk had BAB spore concentrations below 10,000 MPN/L (or 10 MPN/mL). A total of 98.9% of simulated raw milk vats had less than 1,000 MPN/L (1 MPN/mL). The subsequent growth module [module (ii)] predicted that during aging at 13°C, the median concentration (±SD) of C. tyrobutyricum increased to 3.68 (±0.67), 4.21 (±0.70), and 4.74 (±0.73) log MPN/ kg at d 60, 90, and 120 d, respectively (Figure 2). Using these data as input, module (iii) used upper and lower C. tyrobutyricum concentration thresholds to predict the proportion of cheese wheels expected to show LBD at d 60, 90, and 120. For d 60, use of the lower threshold predicted that 10.9% of wheels showed LBD, and use of the upper threshold predicted that 7.5% of wheels showed LBD, which can be summarized as a model prediction of 9.2% (±1.7) LBD for cheese aged for 60 d.
The corresponding values for d 90 and 120 were 36.1% (±3.4) LBD and 69.2% (±2.7) LBD, respectively (see Figure 3 and Table 5 for detailed results).  Model outputs also provided an opportunity to describe the characteristics of the cheeses that were predicted to show LBD at a given time point, such as initial spore count in raw milk and cheese pH (Table 5). For example, cheeses that developed LBD by d 60 based on the upper threshold showed a mean initial spore count in raw milk of 2.80 log MPN/L and a mean cheese pH of 5.43. In comparison, all cheeses that were predicted to show LBD by d 120 (which includes all cheese that showed LBD at d 60 or 90), also based on the upper threshold, showed a mean initial BAB spore count in raw milk of 2.12 log MPN/L and a mean cheese pH of 5.39 (Table 5). Although the mean BAB spore concentration for raw milk predicted to yield cheese with LBD varied widely (e.g., from 3.17 to 2.08 log MPN/L for d 30 and 120, based on the lower threshold), the mean cheese pH showed less variation across cheese late blown at different times, ranging from 5.39 for all cheese with LBD before d 120 to 5.43 for all cheese with LBD before d 30.
A comparison of cheeses predicted to show LBD and cheeses predicted not to show LBD at a given time point also provides information that can be used to identify raw milk batches associated with a higher risk of development of LBD. For example, based on the upper threshold, cheese that showed LBD at d 90 had mean BAB spore count for raw milk and pH of 2.42 log MPN/L and 5.40, respectively, but cheese that did not show LBD at this time point had mean BAB spore count for raw milk and pH of 1.53 log MPN/L and 5.36, respectively.
The results from module (iii) were qualitatively validated against the literature. The earliest occurrence of LBD predicted by the base model is at d 30, which coincides with literature observation and industry data that LBD may occur as early as 1 mo. In our model, most cheeses gradually developed LBD between 1 and 4 mo, which also fits the typical timeframe of LBD previously reported by Düsterhöft and Van Den Berg (2007), D'Amico (2014), and Düsterhöft et al. (2017).

Sensitivity Analysis
Sensitivity analysis was performed to assess the effects of 4 selected model parameters (maximum growth rate, concentration factor, C. tyrobutyricum concentration threshold in cheese, and cheese pH) on the proportional change of late blown cheese. The relative importance of each variable differed between d 60 and d 120 ( Figure  4A and B). At d 60, the concentration thresholds were identified as the most important. A decrease of thresholds by 0.5 log led to a 25.4 percentage point increase in the proportion of cheeses predicted as developing LBD, compared with the base model (i.e., 34.6% of cheese developed LBD at d 60 when thresholds were reduced by 0.5 log, whereas only 9.2% developed the defect at d 60 in the base model). The maximum growth rate ranked second in importance, as an increase of 40% on μ max could lead to 20.8 percentage point more cheeses with LBD. Although several intrinsic and extrinsic factors (e.g., temperature, pH, and salt and lactate concentrations) can influence the growth rate, this result indicates that naturally occurring differences in growth rate (e.g., between different strains of C. tyrobutyricum) can have a large effect on occurrence of LBD in Gouda.
The effects of concentration factor and pH were relatively minimal, resulting in only 4.9 and 5.7 percentage point increases in the number of late blown cheeses at the worst case ( Figure 4). Interestingly, the relative importance of maximum growth rate and concentration  percentage points more late blown cheese, compared with 18.8 percentage points more when the concentration thresholds were decreased by 0.5 log. The pH had more effect toward the end of aging, as an increase of pH by 0.2 can increase the occurrence of LBD by 13.6 percentage points at d 120, compared with only 4.9 percentage points at d 60. Increased metabolism (indicated by higher growth rate and gas production) at higher pH, until the optimal pH (approximately 5.75) is reached, is also supported by previous studies in liquid model systems (Huchet et al., 1995;Ruusunen et al., 2012;Silvetti et al., 2018).

What-If Scenarios
The potential effectiveness of different intervention strategies to reduce the development of LBD in Gouda cheese was evaluated independently through theoretical what-if scenarios (Table 4); for these analyses the average of the percent of LBD predicted by the lower  . Sensitivity analysis testing the relative importance of 4 variables (y-axis: tl = Clostridium tyrobutyricum concentration threshold, μ= maximum growth rate, pH = cheese pH, cf = concentration factor) for the change in proportions of late blown cheese at d 60 and 120 of aging (displayed on x-axis). Values on the x-axis were normalized such that a value of 0 represents 9.2% (d 60) and 69.2% (d 120) of cheese predicted to show LBD. Changes shown thus represent percentage points. The best cases for concentration thresholds are increasing thresholds by 0.5 log (high) and 0.25 log (low); worst cases are decreasing concentration thresholds by 0.5 log (high) and 0.25 log (low). Best cases for maximum growth rate are increasing the maximum growth rate by 40% (high) and 20% (low); worst cases are decreasing maximum growth rate by 40% (high) and 20% (low). Best cases for pH are decreasing pH by 0.2 (high) and 0.1 (low); worst cases are increasing pH by 0.2 (high) and 0.1 (low). Best cases for concentration factor are calculated by increasing the value by 2 (high) and 1 (low); worst cases are calculated by decreasing the value by 2 (high) and 1 (low). Baseline calculations were based on the lower C. tyrobutyricum concentration threshold. and upper concentration thresholds was used as the outcome. Microfiltration or bactofugation (each assumed to have a spore removal efficacy of 98%) were found to be most effective methods for reducing LBD; these approaches were predicted to reduce the occurrence of LBD from 69.2% to 0.6% at d 120 and completely prevent occurrence of the defect during 60 d of aging (0/200,000 wheels were predicted to show LBD). The second and third most effective strategies were the addition of nitrate at 2.5 g per 100 L of milk and lowering aging temperature from 13 to 12°C, which reduced predicted LBD prevalence to 1.2% and 6.3% at d 120, respectively. For the remaining 3 interventions, effectiveness differed substantially based on the aging time (see Table 4). The use of low-spore premium milk with arithmetic mean spore count of 25 MPN/L (compared with 143 MPN/L for the baseline model) was predicted to reduce LBD at d 60 to 1.3% (compared with a baseline of 9.2%), whereas the cumulative proportion of cheese predicted to show LBD at d 120 was 41.8% with and 69.2% without the intervention. Similarly, the inoculation of bacteriocinogenic lactic acid bacteria prevented the majority of LBD at d 60 (0.7% of cheeses were predicted to show LBD), whereas 12.4% of cheeses made by utilizing this intervention were predicted to develop LBD by d 120 (Table 4). The addition of lysozyme was predicted to prevent 87.8% of cheese from LBD throughout the late stage of aging, but 9% of cheese could still develop LBD in the first 2 mo of aging.

DISCUSSION
Here, we report on a one-dimensional MC model that can be used to predict the occurrence of LBD, which is a major quality issue in many semi-hard cheeses. Although we used Gouda cheese as a model, our work provides a framework for developing a digital tool that can aid cheese processors in making decisions necessary to mitigate the negative impact of LBD in different cheese types. Importantly, through a sensitivity analysis, the model used here was able to identify data needs for improved model performance. Although data on C. tyrobutyricum concentration thresholds at onset of LBD were identified as a particular data need, a general sparsity of data on C. tyrobutyricum growth in cheeses and limited data on BAB subtypes and their ability to cause LBD also represent limitations of the existing model. The predictive model also allowed for what-if analysis, which indicated that interventions that allow for substantial reduction in initial spore levels (e.g., bactofugation, microfiltration) have the greatest potential for reducing LBD.

Sensitivity Analysis
The sensitivity analysis performed here suggests that the concentration of C. tyrobutyricum in cheese that defines the concentration threshold at onset of LBD is the most important factor in predicting the proportion of cheese with LBD on d 60. This is important, as only very limited data are available on the concentration of vegetative C. tyrobutyricum cells at the onset of late blowing and because the level of C. tyrobutyricum to induce LBD will likely depend on several factors such as cheese type, shape, and structure, as well as pH, aging time, and temperature (D'Amico, 2014). Hence, this finding suggests that substantial additional data are needed to develop models for LBD in different cheese types. In addition, the few previous studies available reported a wide range in C. tyrobutyricum numbers present in cheeses when they developed LBD. For example, in naturally contaminated Gouda cheese, one study (Herman et al., 1997) reported C. tyrobutyricum vegetative cell concentrations as low as 0.4 MPN/g in cheese at onset of LBD, whereas another study (Toyoda and Nakajima, 1995) reported concentration of 3.5 × 10 3 cfu per gram of cheese. Although these data could indicate considerable variability of the concentration threshold, it is also possible that the MPN method used in the former study to quantity vegetative cells of C. tyrobutyricum in cheese is less reliable; for example, the nested PCR method that was performed in addition to MPN method was reported to have a high detection limit (Herman et al., 1997). Methodological problems in enumerating vegetative C. tyrobutyricum cells may also be responsible for the fact that additional studies (Stadhouders et al., 1985;Le Bourhis et al., 2007;Matijašić et al., 2007;Garde et al., 2011;Morandi et al., 2021) only reported spore concentrations at different stages of cheese aging and not concentrations of viable cells. Even within these studies, patterns of how spore concentrations change during aging differed widely. In one study that inoculated raw milk with C. tyrobutyricum to produce Emmental-type cheeses, the spore level increased with aging time, and at the onset of late blowing reached up to 4.7 × 10 5 spores per gram of cheese (Le Bourhis et al., 2007). However, in a different study where natural contamination of clostridia occurred in semi-hard Valtellina Casera cheese (Morandi et al., 2021), a decrease in spore concentration was observed over aging. At the day of LBD appearance, the spore concentration dropped to less than 0.5 MPN per gram of cheese. As these data indicated a high level of variability of the C. tyrobutyricum concentration threshold at onset of LBD, we elected to calculate the thresholds to use in our model instead of using literature data. Because the calculation of the thresholds was based on data for a single cheese type and aging condition, the variability in threshold was limited, as expected; using these calculated thresholds was appropriate, as our model currently only predicts LBD for one cheese type (Gouda). Overall, however, high variability in threshold levels reported from previous studies (Toyoda and Nakajima, 1995;Herman et al., 1997) along with our results suggest a need for additional research on C. tyrobutyricum concentration at onset of LBD, including characterization of uncertainty and variability (e.g., strain variability). These data will be particularly important as further work tries to develop models that can predict LBD across different cheese types.
After the concentration threshold, the maximum growth rate was identified as the second most important factor in our model. At worst scenario, 20% and 40% increases in the maximum growth rate led to 9.1 and 20.8 percentage points increases in occurrence of LBD, respectively. Although no previous studies have conducted sensitivity analysis for C. tyrobutyricum, the relative importance of the maximum growth rate is consistent with a risk assessment model for fluid milk spoilage due to sporeformers (Buehler et al., 2018a). In that study, growth rate was reported to be the most important factor in determining when sporeformers reached the concentration threshold of 20,000 cfu/mL. In another risk assessment model (Couvert et al., 2010) that estimated the risk of Listeria monocytogenes in cheese curds and smoked salmon, the maximum growth rate was also the most important factor tested in sensitivity analysis. Moreover, the importance of the maximum growth rate was reflected in an illustrative framework of food safety risk assessment (McNab, 1998) in which the maximum growth rate was the second most influential parameter in multiple model scenarios. Although our study did not have sufficient data to incorporate strain-specific growth characteristics for C. tyrobutyricum, previous studies (Ruusunen et al., 2012;Podrzaj et al., 2020) have indicated that C. tyrobutyricum strains can differ in favorable growth temperature, salt tolerance, and pH sensitivity. These phenotypical differences can significantly influence the ability of C. tyrobutyricum to grow in semi-hard cheese and cause LBD. This can lead to different abundance of strains in cheese as well. Thus, to refine the influence of growth rate, future studies should aim to identify the frequency and growth parameters of strains relevant to cheese late blowing.

What-If Scenarios
The what-if scenarios performed here suggest that implementation of microfiltration or bactofugation would have a profound impact in reducing LBD occurrence in Gouda cheese. With a raw milk supply with a mean spore count of 143 MPN/L, this intervention approach prevented LBD for the first 90 d of aging. This reduction translates to 36.1% of total cheese (72,200/200,000) saved from LBD during the first 90 d of aging. It should be noted that this result is specific to the efficiency of bactofugation we applied to our model (98% removal of anaerobic spores). The efficiency, however, is subject to change depending on the specific combination of centrifuge type, applied rate (i.e., speed), and temperature; 2 previous studies specifically reported that spore removal efficiency by bactofugation can range from 44% to 99.7% with a mode of 98% (Waes and Heddeghem, 1990;Gésan-Guiziou, 2010). Waes and Heddeghem (1990) also reported on results from various Dutch studies indicating a high effectiveness of bactofugation in preventing LBD. For example, in one study, a Jarlsberg cheese with no signs of LBD was produced after applying bactofugation to raw milk with a maximum of 1 anaerobic spore per mL (Strand, 1969, as cited in Waes and Heddeghem, 1990). The same chapter also mentioned studies by Bergère (1969) and Bottazzi (1982) that achieved similar inhibition of LBD in Gruyère and Grana cheeses by using bactofugation to lower the spore count in milk to less than 0.05 spores per milliliter. Similar to the spore counts reported in these studies, spore counts in our simulations were very low after bactofugation (with a median spore count of 0.01 spore per milliliter of milk). For this spore level, our model also predicted that bactofugation could minimize the LBD to below 1% of cheese after 4-mo aging at 13°C. Thus, the predictions from our model confirmed the previously reported effectiveness of bactofugation for preventing LBD in different cheese types (Waes and Heddeghem, 1990;Gésan-Guiziou, 2010). However, other studies have showed that bactofugation is not always able to eliminate LBD; for example, in the study of Berg et al. (1980), LBD still occurred after spore concentration in milk was reduced to between 0.25 and 1 per milliliter using bactofugation. This is in agreement with our results, which indicate that 9.2% (±1.7%) of cheese predicted to develop LBD by d 60 would have initial spore concentration between 0.25 and 4.3 spores per milliliter (Table 5). Our model predicted that, to prevent occurrence of LBD during 4 mo of aging, the concentration of spores in milk has to be at 0.01 spore per milliliter or lower. This concentration is substantially lower than the 0.25 to 1 spore per milliliter reported by Berg et al. (1980), which is probably the reason the authors suggested addition of sodium nitrate at 2.5 g per 100 mL of milk for complete prevention of LBD, even after bactofugation. Our model suggests that such addition is not needed when bactofugation can lower the median spore count to 0.01 spore per milliliter; this type of spore level should be achievable with good-quality raw milk and good efficiency of the bactofugation process. Although bactofugation and microfiltration were assumed to have the same effect in our model, in reality use of microfiltration to treat milk used for manufacture of cheeses susceptible to LBD is less practical, as microfiltration can only effectively be used on skim milk, which means that additional steps (e.g., cream separation, heat treatment of cream, and recombining) would need to be introduced as part of the intervention (Maubois and Schuck, 2005;Gésan-Guiziou, 2010).
Although less relevant in the United States, where use of nitrates in cheese is prohibited, adding sodium nitrate to milk was the second most effective intervention strategy demonstrated by our model (ranked after implementation of microfiltration or bactofugation), which is probably the reason why nitrates are still used as an intervention outside of the United States. With our model, using sodium nitrate at 2.5 g/100 mL was predicted to limit the proportions of late blown cheese to only 1.2% (2,400/200,000) over 4-mo aging. This coincides with a previous finding (Waes and Heddeghem, 1990) that adding 2.5 g of NaNO 3 was sufficient to prevent LBD when milk with 0.25 to 1 spore per milliliter was used for cheese manufacturing, a contamination level similar to that used in our model. To facilitate use of sodium nitrate for prevention of LBD, Stadhouder (1990a) previously provided a table of nitrate doses needed to inhibit LBD based on varying spore levels in raw milk used in Gouda production. This information, combined with our model, can be useful to help cheese makers, in countries where use of nitrate is allowed, make decisions on appropriate LBD interventions and also provides an opportunity for future integration of this information into an advanced version of our model.

Limitations
We have identified some limitations in our model. First, our model validation process is qualitative. Although our predictions obtained through what-if scenarios generally agree with literature reports, no previous studies to our knowledge include quantitative data (e.g., proportions of late blown cheese, distribution of C. tyrobutyricum concentration, and more) to statistically validate our model. In future studies, more comprehensive data sets that include parameters such as initial BAB spores, cheese aging conditions, and percentage of late blown cheese need to be procured either from pilot scale experiments or from industry for more rigorous validation. Second, in our secondary growth model, the cheese water activity used in calculation was simplified as the water activity that can be typically found in 60-d Gouda cheese. This simplification approach was selected because incorporating the spatial and temporal variation in cheese water activity is complex and beyond the scope of this study. Third, our model is a one-dimensional Monte Carlo simulation. In other microbial risk assessment models, second-order Monte Carlo simulations were applied, using a resampling method to address uncertainty in model parameters (Pouillot et al., 2004(Pouillot et al., , 2007. In the future, bootstrapping can be applied to our original BAB spore count data set, to obtain a distribution of hyperparameters in the central tendency and variance of spore count. In addition, although a dose-response relationship is typical in risk assessment models, our model also did not include such function to quantify the correlation between microbial concentration and risk of LBD, due to the lack of experimental data. Rather, we used upper and lower concentration thresholds to represent the uncertainty in model prediction. Future research is needed to quantify this dose-response function, to better characterize the risk.

CONCLUSIONS
Although cheese and other dairy manufacturers increasingly collect larger amounts of data, availability of data analytics and predictive tools for the dairy industry is still a bottleneck that negatively affects the ability of the dairy industry to make use of advancements in digital technologies. The model developed here not only represents a beta version of a decision-making tool that could be valuable for certain cheese makers, but also illustrates the types of digital tools that could be valuable for the dairy industry, particularly with enhanced capture of data (e.g., raw milk spore counts, water activity, and more). With integration of our model into R ShinyApp, a user-friendly interface (https: / / foodsafety .foodscience .cornell .edu/ digital -dairy/ cheese -blowing -model/ ), we provided an example of how tools developed in a sophisticated statistical environment (i.e., R) can be easily made available, enhancing deployability of digital tools developed in academia.