Primary Model for Biomass Growth Prediction in Batch Fermentation

Predictive models may be considered a tool to ensure food quality as they provide insights that support decision making on the design of processes, such as fermentation. Objective: To formulate a mathematical model that describes the growth of lactic acid bacteria (LAB) in batch fermentation. Methodology: Based on real-life experimental data from eight LAB strains, we formulated a primary model in the form of a third-degree polynomial function that successfully describes the four phases observed in LAB growth, i.e., lag, exponential, stationary, and death. Our cubic mathematical model allows us to understand the fundamental nonlinear dynamics of LAB as well as its time-variant dependencies. Parameters of the model are written in terms of initial biomass, maximum biomass, maximum growth rate, and lag phase duration. Further, a statistical analysis was performed to compare our cubic primary model with the ones proposed by Gompertz, Baranyi, and Vázquez-Murado by computing the coefficient of determination ( R2 ) , the residual sum of squares (RSS), and the Akaike Information Criterion (AIC). Results: The average statistical results from the cubic model are as follows: R2 = 0.820 providing a better fit than the other three models, RSS = 0.658 and AIC = −6.499, where both values are lower than the other models considered in this study. Conclusion: The cubic primary model formulated in this work describes the behavior of biomass as it accurately represents the four phases of biomass growth in batch fermentation process.


Introduction
Lactic Acid Bacteria (LAB) refer to a group of microorganisms that share certain morphological, physiological, and metabolic characteristics. They have the peculiarity of producing lactic acid from several carbohydrates through a process known as microbial fermentation. They are normally found in cultures such as milk, whey, and pickles [1,2]. The Danish microbiologist Orla-Jensen said that "True lactic acid bacteria are a large group of cocci and Gram positive, immobile, spore free bacilli, which produce lactic acid in the fermentation of sugar" [3]. Currently, LAB have a primary role in the food industry, as they are used to acidify and preserve food. Further, they contribute to texture, taste, smell, and scent development in all kinds of fermented foods [4,5].
In the industry, most fermentation is carried out through batch culture. This means that only a culture medium with the necessary nutrients is added to the fermenter, and this is incubated for the growth of bacteria and the production of its primary metabolite. During the time of incubation, nothing other than oxygen is added [6]. Hence, the bacterial growth follows an asymmetric cell division (ACD), which is a conserved mechanism evolved to generate cellular diversity. A key principle of ACD is the establishment of distinct sibling cell fates by mechanisms linked to mitosis [7]. ACD also occurs in the development and physiology of unicellular organisms ranging from bacterial species to yeasts and flagellates [8,9]. In these organisms, ACD underlies replicative aging as a means of maintaining the immortality of the mitotically proliferating population [10].
LAB are very demanding microorganisms and need a set of growing factors, including sugars such as carbohydrates, as well as amino acids and vitamins. Therefore, LAB can only grow in a medium that supplies these nutrients [3]. The interest in the physiology of LAB has been studied by its industrial importance and potential use of genetic engineering in strain optimization. For example, in the particular case of Lactococcus lactis spp. lactis, a minimal growth medium should contain glucose, acetate, vitamins, and amino acids [11]. Milk is the medium that contains all these nutrients. For this reason, LAB are used as starter cultures for the preparation and preservation of dairy products, such as acidified milk, yogurt, butter, cream, and cheese [1]. In batch fermentation, bacteria cannot grow exponentially indefinitely. The latter due to bacteria depleting nutrients as it grows, changes in the chemical composition of the medium (pH), and toxic compounds that are accumulated. The characteristic curve of bacterial growth is composed by four phases: (i) lag, (ii) exponential, (iii) stationary, and (iv) death. Lactic acid production occurs during the reproduction of LAB; the accumulation of this and other organic acids decrease the pH of the medium, which made culture conditions become more selective; hence, the more acid-tolerant bacterium will prevail.
The total number of cells in fermentation are usually represented by a logarithmic scale, and normally, at the end of this process, one can see values from 10 8 up to 10 9 CFU/mL (Cell-Forming Units per milliliter) [2]. In addition to nutritional requirements, temperature is one of the most important factors in LAB growth. Furthermore, there is an optimal temperature and other environmental conditions for which these microorganisms present a higher growth rate [4]. LAB are acid-tolerant and may grow either at pH values as low as 3.2 or as high as 9.6. Usually, most of LAB strains grow at pH values between 4 and 4.5 [1].
In microbiology, several mathematical models have been developed to predict the behavior of microorganisms on the influence of different factors in culture. This branch is known as predictive microbiology, and it aims to study the response of bacteria to environmental factors that can be controlled, such as temperature, pH, and water volume, among others. In this sense, a tool may be formulated to ensure both quality and safety of products to estimate their useful lifetime and to make decisions regarding the composition and design of fermented products [12]. Mathematical models that relate biomass production and time are called primary models, while those that describe the relationship between primary models parameters and environmental conditions are known as secondary models [13]. Some of these models are discussed below.
The Gompertz equation, formulated in 1825, is one of the most used in predictive microbiology. It is based on an exponential relationship between the growth rate and the density of a population. This equation was formulated to describe the law of human mortality. However, years later, it was adapted and re-parametrized for its use in microbiology [14]. Giraud et al. [15] observed that pH variations resulted in a decrease in the growth rate when measuring the fermentation performance of L. plantarum at a controlled pH between 4.0 and 8.0. Baranyi et al. [16] applied a non-autonomous differential equation to describe the dynamics of growing bacterial cultures. Based on more than 500 growth curves, the statistical properties of this equation were compared to the Gompertz approach, which is the most commonly used in food microbiology. After these results, Baranyi and Roberts [17] proposed a growth model where a single variable represents the physiological state of cells. The lag phase period is determined by the value of this variable at inoculation and by the post-inoculation environment. The model was able to describe bacterial growth in an environment where factors, such as temperature and pH, change over time. Nicolai et al. [18] constructed a dynamical model for the growth of LAB in vacuum-packed meat. The model was divided into two parts: one part describing the fermentation and the other describing the pH evolution in the liquid surface layer. These models were given as two differential equations and two algebraic equations, respectively. Passos et al. [19] developed an unstructured model to describe bacterial growth, substrate utilization, and lactic acid production by L. plantarum in cucumber juice. They also developed an equation that relates the specific mortality rate and the sodium chloride (NaCl) concentration. Drosinos et al. [20] formulated an empirical model to describe the growth and production of bacteriocin by Leuconostoc mesenteroides E131 under different conditions of pH and temperature. Further, a De Man, Rogosa, and Sharpe (MRS) broth was used as a growth medium in the fermenter. Vázquez and Murado [21] proposed a model based on the re-parametrized logistic equation to describe fermentation kinetics of lactic acid production by Lactococcus lactis and Pediococcus acidilactici in a batch system. Da Silva et al. [22] studied the growth of L. plantarum, W. viridescens, and L. sakei under different isothermal culture conditions at temperatures of 4,8,12,16,20, and 30 • C. They determined that LAB growth was strongly influenced by the culture temperature and created new models that allowed them to predict growth at temperatures ranging from 4 to 30 • C. Dalcanton et al. [23] built a response model of the growth rate of L. plantarum as a function of temperature, pH, and concentrations of NaCl and sodium lactate (NaC 3 H 5 O 3 ).
Despite the existence of several mathematical models in the literature that describe the behavior of LAB during the fermentation process in its three phases, i.e., lag, exponential, and stationary phases, they do not include the death phase in these models. In batch fermentation, after the stationary state, the depletion of nutrients in the culture medium occurs, and therefore, the death phase begins and should be considered in these models. By not taking the death phase into account, a problem arises when trying to fit real-life experimental data with the models usually used in predictive microbiology. Therefore, the objective of this research is to formulate a time-variant mathematical model that describes the growth of LAB, including its death phase in batch fermentation.
The remainder of this paper proceeds as follows. In Section 2, we present the real-life experimental data concerning biomass growth for eight LAB strains; we formulate our cubic mathematical and explain each parameter; and the most important kinetic models identified in the literature are explored. In Section 3, a statistical analysis is performed to establish which model better fits the experimental data, and these results are illustrated by means of several numerical simulations. Finally, discussions are described in Section 4, and conclusions are given in Section 5.

Experimental Data
Eight LAB strains were investigated in the present study; all of them isolated from autochthonous fermented milk. This process was performed based on certain desirable characteristics and against microbial spoilage [24]. LAB strains were incubated for 48 h at a temperature of 37 • C in reconstituted 10% milk powder to produce acidification. Biomass production was measured every 6 h by means of the Neubauer cell counting chamber. Two repetitions were performed for each strain, and the average was calculated as the final concentration value at each corresponding hour. The total biomass concentration was written on a logarithmic scale with units given in log 10 (CFU/mL). Furthermore, the experimental data were conditioned by a moving mean filter to minimize noisy measurements in the biomass concentration and to better illustrate LAB growth in each of its four phases, which we mentioned before as lag, exponential, stationary, and death. The results are illustrated in Figure 1, and the overall dynamics allow us to formulate a third-degree polynomial that will be discussed in the following section.

Mathematical Modelling: Cubic Polynomial Model
In order to formulate a mathematical model to describe the time-evolution for the eight different strains of LAB in the 48 h of the process (see Figure 1), we proposed the following third-degree polynomial where x β (t) describes biomass growth dynamics in log 10 (CFU/mL) as a function of time, x 0 is the initial biomass in each experiment, and coefficients a, b > 0 represent the relationships between maximum biomass, lag phase time, and maximum growth rate. Differential calculus concepts were applied to write coefficients as functions of the desired set of parameters. Now, let us compute the first derivative, as it represents growth rate as followṡ to find the maximum increase in biomass, which is denoted as x max , we calculate the local extremums as indicated below t(2a − 3bt) = 0, from the latter, it is evident that at t = 0, there is a minimum. Therefore, the time when the maximum biomass increase occurs is given at Then, by the substitution of Equation (3) into Equation (1), one can calculate x max as follows Now, by computing the second derivative of Equation (1), which is given beloẅ and the time to the inflection point is computed as follows from Equation (5) then, this value is evaluated in Equation (2) to calculate the maximum growth rate, which is denoted as µ max . This parameter represents the slope of the tangent at the inflection point, and its value is given by the next result The total amount of biomass at the inflection point, which is given as h, is obtained by evaluating Equation (6) into Equation (1), and the corresponding result is given below therefore, the inflection point is located at the following coordinate a 3b , Now, in order to find the lag time-period (L), the intersection of the tangent line on the inflection point at the t-axis was calculated, as illustrated in Figure 2. Hence, by applying the point-slope equation, we formulate the following and by isolating L from Equation (9), the following result is obtained then, replace T PI and h with Equations (6) and (8), respectively. Thus, L is rewritten as follows Equations (4) and (10) relate parameters of interest with coefficients a and b of Equation (1). Therefore, one can isolate a from Equation (10) as follows and substitute it into Equation (4) to find the value of b as indicated below then, Equation (12) should be replaced by Equation (11) to determine the final value of coefficient a; hence Values of a and b, given by Equations (12) and (13), respectively, are replaced into Equation (1) to formulate the model in terms of parameters x max , µ max and L. Therefore, our cubic polynomial model to describe the four phases of LAB growth is shown below

Primary Models
Numerical simulations were made to illustrate and analyze the time-evolution of our model and the other three chosen from the literature. These models are presented below. First, let us introduce the Gompertz model [14]: Then, the Baranyi model [17], which is given as follows: where A(t) and h 0 may be written as indicated below [25]: Now, the model from Vázquez-Murado [21]: .
In each model, x 0 represents the initial biomass, x max the maximum biomass, µ max the maximum growth rate, and L the lag phase time. Further, these three models were selected because they are among the most cited in predictive microbiology, and their parameters are the same as in our cubic primary model (14). Parameter values were calculated from the smoothed data of the eight LAB strains (see in Figure 1). Parameters x 0 and x max are, respectively, the initial and the maximum biomass values. The value of the parameter µ max is computed as shown below where S i+1 and S i are taken from the smoothed experimental data, i.e., a given value and its previous one. The latter is divided by the time interval of the measurements, which, in this particular case, ∆t = 6 h, and since we have nine measurements, then i = 0, . . . , 8. Furthermore, it should be noted that Equation (10) was used to calculate the Lag time period, L. The results concerning all necessary parameter values are shown in Table 1.  Figure 1.

Results
In order to validate our model, a statistical analysis was performed to calculate and compare the descriptive capacity of the cubic model (14) and models (15)- (17). The statistic used was the coefficient of determination R 2 , which measures the goodness of fit of a non-linear model according to experimental data. Its calculation is made based on the sum of squared residuals (RSS), and the explanatory capacity of a model is better the closer R 2 is to 1, see Section 11-8.2 [26]. Figures 3-6 illustrate the dynamics of the cubic model and its comparison with analyzed models from the literature concerning the experimental data of eight LAB strains, S i (t), i = 1, . . . , 8. One can see in Figure 3 that the cubic model describes the four expected growth phases in batch fermentation, while the Gompertz model, see Figure 4; the Baranyi model, see Figure 5; and Vázquez-Murado model, see Figure 6; only describe up to the stationary phase, which caused the RSS value to be higher for these models, and therefore, the results concerning R 2 and AIC indicate a poorer fit. The lag phase time was around 5 to 9 h. The maximum biomass growth was reached between 30 and 36 h for each of the strains, and from that time, the death phase begins.
According to our numerical simulations, one can observe that the proposed cubic model (14) and the other primary models (15)-(17) have a good adjustment in the lag phase, as is shown by each solution and the corresponding experimental data. Concerning the exponential phase, the Baranyi model is slightly above the others. However, this model has a lower growth rate, which causes a delay when reaching the maximum concentration, approximately at 48 h, while other models reach this value at around 30 or 36 h. Regarding the maximum biomass growth, the three primary models remain in the stationary phase, while the experimental data goes to the death phase. Therefore, the RSS of these models is higher than the cubic model, as this one better fits the observed data in this last phase. Table 2 shows the results of R 2 and RSS for the cubic, Gompertz, Baranyi, and Vázquez-Murado models with respect to the experimental data illustrated in Figure 1. Now, concerning strain 1, the cubic model presented an RSS of 0.843 and, consequently, a higher value of R 2 given by 0.805 compared to Gompertz, Baranyi and Vázquez-Murado models, which obtained an RSS equal to 1.431, 1.920, and 1.392 and an R 2 equal to 0.669, 0.556, and 0.678, respectively. For strain 2, the highest R 2 value was for the cubic model with 0.793, followed closely by the Vázquez-Murado and Gompertz models with 0.723 and 0.722, respectively, while the Baranyi model had a result of 0.641. In the particular case of strain 3, all models had a lower adjustment when comparing them with the other seven strains. However, the best fit was for the cubic model with an These higher R 2 values are due to fact that the experimental data for this strain have fewer outliers, as seen in Figure 1, and its death phase is lower compared to the other strains. Strain 7 had a better fit for the cubic model, with an R 2 of 0.792, than the Gompertz (0.590), Baranyi (0.491), and Vázquez-Murado (0.614) models. Finally, as in most cases, the cubic model had the best fit for strain 8 with an R 2 value of 0.898, while Gompertz and Vázquez-Murado obtained, respectively, 0.756 and 0.766. In general, the cubic model had higher values of R 2 than the rest of the models, except in strain 6, where the Vázquez-Murado model had the better fit to the experimental data. The latter was due to the data not yet presenting the death phase as the other strains in the 48 h of the experiment.
Furthermore, the Akaike Information Criterion (AIC) test was performed, which allows us to determine which model better fits the observed data. To calculate this value, the goodness of fit between predictions of models and experimental data is considered through the RSS while penalizing models that have a greater number of parameters due to these becoming more complex for its practice [27]. Table 3 shows results for every mathematical model. By taking the latter into account, the model that better fits the experimental data is the one with the lowest AIC value. For strain 1, the lowest AIC was for the cubic model with a value of −3.310, while Gompertz, Baranyi After analyzing the statistical results, it is evident that our formulated cubic primary model better fits the experimental data because it was the one that presented the lowest RSS in general and therefore had higher values for the R 2 in all LAB strains but number 6, which had a higher value with the Vázquez-Murado model (0.960). However, the result for the cubic model was of 0.912, which is still a good value concerning this parameter. Therefore, based on the results for the RSS, R 2 , and the AIC, we can establish that the cubic primary model (14) better represents the overall dynamics of the experimental data for the eight LAB strains under study in this research.

Discussion
In predictive microbiology, primary models based on sigmoidal functions are usually used. Among the most common models, we have the Gompertz [14], the Baranyi [17], and the Vázquez-Murado [21] models. Nonetheless, the latter only describe the first three growth phases, i.e., lag phase, exponential growth, and the stationary state. According to Buchanan et al. [28] and Garre et al. [29], these phases are sufficient to fit experimental data for growth models. However, according to Mandigan et al. [6], in batch fermentation, the bacterial population presents a death phase due to nutrient depletion in the culture medium. Therefore, when adjusting these models to observed data in this kind of fermentation, values of greater discrepancy are found in the death phase between the experimental data and approximated values. Chatterjee et al. [30] incorporated the death phase into the Gompertz model to better describe the behavior of E. coli and S. aureus. Studies by Bevilacqua et al. [31] were focused on demonstrating the importance of the death phase of microorganisms and its impact on the shelf life of food. Solano et al. [32] evaluated the capacity of five models, including the Gompertz, Baranyi, and Vázquez-Murado, to predict acid production in lactic fermentation of fishing products, and they found that models with the lowest residual variance were those of Gompertz and Baranyi. Further, this study shows that the Vázquez-Murado model failed to give an adequate adjustment for lactic fermentation. In the same way, Zwietering et al. [33] managed to describe the behavior of L. plantarum in an MRS medium with the Gompertz model, while the Vázquez-Murado logistic model was not suitable to accurately approximate the experimental data. In contrast, Kedia et al. [34] obtained a good fit with the Vázquez-Murado model for L. reuteri, L. plantarum, and L. acidophilus in oat fractions. Da Silva et al. [22] studied the growth of L. plantarum, W. viridescens, and L. sakei in vacuum-packed meats, and, as well as Baty and Delignette [35] who studied different growth models for various LAB, all authors in these two works found that the Baranyi model has a slightly better fit than the Gompertz model.
According to the R 2 and AIC values, it can be established that the cubic model has a better fit to the experimental data with average values of R 2 = 0.820 and AIC = −6.499. There was not enough difference between the average values for the Gompertz model with R 2 = 0.693 and AIC = −2.248, while the Vázquez-Murado had average results of R 2 = 0.695 and AIC = −2.366. Regarding the Baranyi model, the average results were as follows R 2 = 0.618 and AIC = 0.785. Further, the overall adjustment to the experimental data for the four primary models under study is illustrated in  In this research, the four primary models presented a good adjustment in the lag phase (L) and reached the maximum growth (x max ) registered in the experimental data. However, the maximum growth rate value (µ max ) had a greater impact on the adjustment of the exponential growth phase, in which the cubic model was the one that had the better fit. Nonetheless, in order to get a good fit with primary models, it is necessary to measure the largest amount of experimental data with the greatest possible continuity, which will ensure a better prediction by the models. In the batch fermentation process performed for this work, the eight LAB strains showed outliers, which could be due to the method applied for the measurements. Therefore, the experimental data had to be smoothed to calculate all required parameters for each mathematical model. It is important to consider that outliers are present in all real-life data measurements, such as biomass growth.
In the cubic model, the parameter µ max also influences the death phase of the biomass. The latter because a higher growth rate implies that the maximum biomass concentration will be reached, and due to the inherit dynamics of this kind of function, after reaching this maximum point, the biomass begins to decline. It is important to consider that if this mathematical model is simulated at time values further than 70 h, biomass will take negative values, which is not biologically feasible in real-life scenarios. However, this can be neglected because the fermentation process to produce acidified milk generally takes around 48 h, and it is not necessary to predict the behavior of this variable beyond this time, as this is the final desired product. Predictive models aim to ensure food quality and provide a tool that supports decision-making on fermentation design. Therefore, based on the results and as it is illustrated in Figures 3-6, a decision can be made to perform a fermentation process under the same culture conditions only up to 36 h, i.e., when the maximum biomass growth is reached for this type of fermented milk.

Conclusions
A time-variant mathematical model given in the form of a third-order polynomial was formulated to describe the growth dynamics over a period of 48 h for eight LAB strains. Our so-called cubic primary model (Equation (14)) was constructed with the parameters that are frequently used in other primary models, i.e, biomass initial concentration, (x 0 ); lag phase time, (L); maximum growth rate, (µ max ); and maximum growth, (x max ). Furthermore, both the numerical simulations and the statistical results concerning RSS, R 2 , and the AIC support our statement that the cubic model x β (t) accurately represents the four biomass growth phases in the batch fermentation process when comparing it with three other models that are among the most cited in predictive microbiology literature.