Dynamic Modeling of Carnobacterium maltaromaticum CNCM I-3298 Growth and Metabolite Production and Model-Based Process Optimization

Carnobacterium maltaromaticum is a species of lactic acid bacteria found in dairy, meat, and fish, with technological properties useful in food biopreservation and flavor development. In more recent years, it has also proven to be a key element of biological time–temperature integrators for tracking temperature variations experienced by perishable foods along the cold-chain. A dynamic model for the growth of C. maltaromaticum CNCM I-3298 and production of four metabolites (formic acid, acetic acid, lactic acid, and ethanol) from trehalose in batch culture was developed using the reaction scheme formalism. The dependence of the specific growth and production rates as well as the product inhibition parameters on the operating conditions were described by the response surface method. The parameters of the model were calibrated from eight experiments, covering a broad spectrum of culture conditions (temperatures between 20 and 37 °C; pH between 6.0 and 9.5). The model was validated against another set of eight independent experiments performed under different conditions selected in the same range. The model correctly predicted the growth kinetics of C. maltaromaticum CNCM I-3298 as well as the dynamics of the carbon source conversion, with a mean relative error of 10% for biomass and 14% for trehalose and the metabolites. The paper illustrates that the proposed model is a valuable tool for optimizing the culture of C. maltaromaticum CNCM I-3298 by determining operating conditions that favor the production of biomass or selected metabolites. Model-based optimization may thus reduce the number of experiments and substantially speed up the process development, with potential applications in food technology for producing starters and improving the yield and productivity of the fermentation of sugars into metabolites of industrial interest.


Introduction
Carnobacterium maltaromaticum is a psychotropic species of lactic acid bacteria widely found in food such as dairy products, fish, and meat. It is a Gram-positive, facultative anaerobic bacterium, able to grow at alkaline pH (up to 9.6) [1,2].
In the food industry, C. maltaromaticum has potential applications related to health protection and organoleptic properties. These include the biopreservation of food, by inhibiting the growth of foodborne pathogens such as Listeria sp. in cold conditions, and the development of flavor in ripened cheese varieties [2][3][4].
This lactic acid bacterium may also be used as a biological indicator in time-temperature integrators (TTI): 'smart-labels' that monitor the time-temperature history of chilled products throughout the cold-chain [5,6]. Concentrates of the strain CNCM I-3298 have been selected as inoculum for TopCryo ® labels, the only biological TTI that has been taken to market to date. A pH decline of the label medium, associated with bacterial growth and End products are shown in blue. ACK, acetate kinase; ADH, acetaldehyde dehydrogenase; LDH, lactate dehydrogenase; PFL, pyruvate formate lyase; PTA, phosphate acetyltransferase; PYK, pyruvate kinase; TreH, neutral trehalose. Adapted from [19][20][21][22][23].
For optimization purposes, modeling has proven to be a powerful tool, enabling the exploration of a wider range of operating conditions while minimizing cost, compared with the experimental approach [24][25][26][27][28][29]. To our knowledge, the only dynamic model dealing with C. maltaromaticum strains has been published by Ellouze et al. [6]. That research was oriented towards a biological TTI setting associated with a sausage-like packaging instead of a bioreactor and taking into account lactic acid as the single metabolite.
The aim of this study was thus to develop and validate a dynamic model predicting the impact of fermentation conditions (temperature and pH) on the growth and bioconversion fermentation dynamics of C. maltaromaticum CNCM I-3298 using trehalose as a carbon source and considering the four main identified metabolites: formic acid, acetic acid, lactic acid, and ethanol. This study was conducted as part of a research project on the production and conservation of C. maltaromaticum concentrates. In that context, the growth of C. maltaromaticum was tested in different sugars: glucose, maltose, mannitol, and trehalose, with similar growth rates. Trehalose was chosen in this study because this molecule is known for its ability to protect cells during bacterial stabilization processes (freeze-drying in particular). Therefore, the residual trehalose (not consumed during fermentation) could be used as cryoprotectant after production of bacterial concentrates.
The model development involved four major steps, presented in Section 3: derivation of the main governing equations based on the known mixed-acid fermentation pathway, mass balances, and kinetic rate expressions (Section 3.1); parameter identification for each fermentation experiment (Section 3.2); construction of response surfaces of the calibrated parameters as a function of temperature and pH (Section 3.3); and final validation of the complete model. The resulting model is shown to be a useful tool in determining the optimal conditions for producing bacterium concentrates in bioreactors and for assessing the productivity of the bioconversion fermentation of sugars into metabolites of potential industrial interest (Section 4.4).

Materials and Methods
Data used to calibrate and validate the model were partially reported in a previous study, in which a modified central composite experimental design was carried out to study the effect of operating conditions on the technological properties of C. maltaromaticum CNCM I-3298 [7]. Sixteen lab-scale fermentations (hereafter named F01 to F16) were performed using a wide range of regulated operating conditions ( Figure 2): temperature between 20 and 37 • C and pH between 6.0 and 9.5. and trehalose, with similar growth rates. Trehalose was chosen in this study because this molecule is known for its ability to protect cells during bacterial stabilization processes (freeze-drying in particular). Therefore, the residual trehalose (not consumed during fermentation) could be used as cryoprotectant after production of bacterial concentrates. The model development involved four major steps, presented in Section 3: derivation of the main governing equations based on the known mixed-acid fermentation pathway, mass balances, and kinetic rate expressions (Section 3.1); parameter identification for each fermentation experiment (Section 3.2); construction of response surfaces of the calibrated parameters as a function of temperature and pH (Section 3.3); and final validation of the complete model. The resulting model is shown to be a useful tool in determining the optimal conditions for producing bacterium concentrates in bioreactors and for assessing the productivity of the bioconversion fermentation of sugars into metabolites of potential industrial interest (Section 4.4).

Materials and Methods
Data used to calibrate and validate the model were partially reported in a previous study, in which a modified central composite experimental design was carried out to study the effect of operating conditions on the technological properties of C. maltaromaticum CNCM I-3298 [7]. Sixteen lab-scale fermentations (hereafter named F01 to F16) were performed using a wide range of regulated operating conditions ( Figure 2): temperature between 20 and 37 °C and pH between 6.0 and 9.5. Fermentation durations varied between 20 h and 45 h, and the initial conditions were: for biomass (X0) 0.077 molC·L −1 , trehalose (S0) between 0.091 mol·L −1 and 0.107 mol·L −1 , and medium volume (V0) 3.5 L.
The main fermentation settings and the kinetic measurements are reported below.  Fermentation durations varied between 20 h and 45 h, and the initial conditions were: for biomass (X 0 ) 0.077 mol C ·L −1 , trehalose (S 0 ) between 0.091 mol·L −1 and 0.107 mol·L −1 , and medium volume (V 0 ) 3.5 L.
The main fermentation settings and the kinetic measurements are reported below.
2.1. Fermentation 2.1.1. Culture Medium and Bacterial Strain The fermentation medium was composed of the following ingredients for 1 kg of final solution: 40 g of trehalose (Treha™; Tokyo Japan); 10 g of proteose peptone (Oxoid; Waltham, MA, USA); 5 g of yeast extract (Humeau; La-Chapelle-sur-Erdre, France); 5 g of Tween 80 (VWR; Leuven, Belgium); 0.41 g of MgSO 4 (Merck; Darmstadt, Germany); 0.056 g MnSO 4 (Merck; Darmstadt, Germany); and water to reach a total of 1 kg of solution. All medium components were sterilized together at 121 • C for 20 min. Fermentations were carried out on C. maltaromaticum CNCM I-3298 pre-cultures. Pre-cultures were prepared by inoculating 10 mL of sterilized fermentation medium with 100 µL of C. maltaromaticum CNCM I-3298 stock culture and were incubated for 13 to 16 h at 30 • C. An amount of 1 mL of the resulting culture was transferred into 50 mL of fresh medium and then incubated again for 11 h under the same conditions. The resulting culture was then used to inoculate the bioreactor. Inoculation was performed at an initial concentration of approximately 10 7 CFU mL −1 .

Bioreactor and Parameter Control
The bioreactor (Minifors, Infors HT, Bottmingen, Switzerland) had a total volume of 5 L and was equipped with a heat mantle and a cryostat for temperature control. It contained 3.5 L of fermentation medium, inoculated with an initial cell concentration of approximately 10 7 CFU·mL −1 . Initial pH was adjusted to the desired value with 5 M NaOH or 0.01 M H 2 SO 4 solutions. During fermentation, pH was controlled to the desired setpoint for each investigated condition ( Figure 2) by automatic addition of 5 M NaOH. Culture homogenization was performed with an agitation device set at 150 rpm. Temperature was set according to the investigated operating conditions mentioned above ( Figure 2).

Cell Growth
Cell growth was monitored using an infrared probe (Excell210, CellD, Roquemaure, France) continuously measuring absorbance at 880 nm and storing data every minute. The absorbance data were calibrated in dry weight. Dry cell weight was determined by filtering 10 mL of bacterial suspension (straight out of the bioreactor) through a 0.20 µm polyethersulfone membrane (Supor ® , PALL Biotech, Saint-Germain-en-Laye, France). The filter was then dried for 24 h at 80 • C. Measurements were obtained in triplicate. Mass concentrations were finally converted to mol C L −1 (carbon-mol of biomass per liter), assuming the simplified unit-carbon biomass formula CH 1.8 O 0.5 [30].

Total Acid Production
Total acid production was determined according to the volume of NaOH solution injected into the bioreactor to maintain a constant pH. The pH was regulated/controlled to set values using the IRIS NT V5 software (Infors, AG, Bottmingen, Switzerland).

Substrate Consumption and Metabolite Production
Trehalose consumption and metabolite production were determined using highperformance liquid chromatography (HPLC, Waters Associates, Millipore; Molsheim, France). HPLC was performed on culture media samples of a few mL, aseptically retrieved from the bioreactor at different times during fermentation and filtered through 0.22 µm pores (Sartorius stedim, Biotech; Göttigen, Germany). Analyses were made using a cation exchange column (Aminex Ion Exclusion HPX-87 300 × 7.8 mm, Bio-Rad, Richmond, VA, HPLC analysis showed that C. maltaromaticum CNCM I-3298 produced not only lactic acid but also formic acid, acetic acid, and ethanol in variable proportions according to the fermentation conditions.

Dynamic Model
The mathematical model was a set of ordinary differential equations implemented in MATLAB R2018b (the MathWorks Inc. Natick, MA, USA). Model parameters and response surface coefficients were identified by nonlinear regression analysis using the Statistic and Machine Learning Toolbox of MATLAB.

Model Formulation
The dynamic model developed in this study combined biochemical knowledge about the metabolism of the selected bacterium and mass balances of the main compounds: substrate, biomass, and identified metabolites. Expressions of specific growth and metabolite production rates included substrate limitation, product inhibition phenomena, and time lags due to microbial metabolism adaptation [31]. The surface response method was used to express the empiric dependence of some model parameters on operating conditions. The model assumed the bioreactor was perfectly stirred and there were no differences between individual cells. It was thus unsegregated and zero-dimensional, predicting average spatial concentrations [32].
Seven  Figure 1) and the culture medium volume (V). This latter variable varied continuously with the addition of base (NaOH) for pH control but also changed in a discrete way due to periodic sampling for biological and chemical analysis.
Mass balances for the considered metabolites resulted in the following set of differential equations: Here, µ X is the specific growth rate (h −1 ); π F , π A , π L , and π E are the specific production rates of four metabolites (h −1 ); and Y X/S , Y F/S , Y A/S , Y L/S , and Y E/S are the yield of biomass and metabolites with respect to the substrate (mol.mol −1 ). Q is the experimentally measured rate of NaOH solution (L.h −1 ) added for pH control throughout fermentation.
In Equation (6), [A T ] is the total acid concentration, defined as the sum of formic, acetic, and lactic acid concentrations. These compounds are assumed to be mainly responsible for the pH change of the liquid medium. Specific growth and production rates were defined using the Monod law to account for substrate limitation, modified with product inhibition and enzymatic activation factors [33][34][35]: In these equations, I X and I m are inhibition factors that depend on the inhibitor concentration. They vary between 1 and 0. Inhibition increases with the inhibitor concentration, and its effect on the specific rate is maximal when the corresponding factor is 0. In this model, progressive inhibition factors of the following form were used [36,37]: K IX and K Im represent characteristic concentrations of the inhibitors (mol L −1 ) such that the corresponding rates (µ X and π m ) are reduced by a factor of 2 compared with the absence of inhibitor, n and p are shape factors, and C I is the concentration of the inhibitor. Since all the metabolites were produced in similar proportions and no biochemical knowledge about their relative inhibiting nature was available, C I was simply defined as the sum of the four metabolite concentrations: To illustrate the role of the shape factor n, Figure 3a depicts the evolution of I X with C I for different n values and a lag-time of 5 h. A more or less sharp change in the inhibition factor occurs around the characteristic inhibitor concentration, C I = K IX . The significance of the shape factor p is similar. Specific growth and production rates were defined using the Monod law to account for substrate limitation, modified with product inhibition and enzymatic activation factors [33][34][35]: In these equations, IX and Im are inhibition factors that depend on the inhibitor concentration. They vary between 1 and 0. Inhibition increases with the inhibitor concentration, and its effect on the specific rate is maximal when the corresponding factor is 0. In this model, progressive inhibition factors of the following form were used [36,37]: KIX and KIm represent characteristic concentrations of the inhibitors (mol L −1 ) such that the corresponding rates (µX and πm) are reduced by a factor of 2 compared with the absence of inhibitor, n and p are shape factors, and CI is the concentration of the inhibitor. Since all the metabolites were produced in similar proportions and no biochemical knowledge about their relative inhibiting nature was available, CI was simply defined as the sum of the four metabolite concentrations: To illustrate the role of the shape factor n, Figure 3a depicts the evolution of IX with CI for different n values and a lag-time of 5 h. A more or less sharp change in the inhibition factor occurs around the characteristic inhibitor concentration, CI = KIX. The significance of the shape factor p is similar.  The enzymatic adaptation factor E A is an empirical representation of the lag time, a period of adaptation to the culture environment where the microorganism produces new enzymatic machinery [38][39][40]. Based on the shape of experimental data, the following equation was proposed: where t lag (h) is the lag time experimentally observed. Figure 3b shows that E A is an increasing function of time, tending to 1 when t t lag . In analogy with n, r is a shape factor that describes the gradual transition from the lag phase to the active phase of growth. A higher value of r implies a steeper change of E A around t = t lag .
To illustrate the features of the proposed model, a representation of the dimensionless specific growth and production rates (µ/µ max and π/π max ) over time is depicted in Figure 4. The dynamic behavior of both variables is similar given the similarity of Equations (1)-(5). The specific rates achieve a maximum value in the active growth phase, and they are zero when t t lag and when the substrate is depleted. The shape of the curve is defined by three factors: in the increasing region (0 to 10 h in Figure 4), the dominant effect is enzyme activation E A (Equation (14)); in the slowly decreasing region (10 to 20 h), the rate is controlled by inhibition (Equation (11) or (12)), whereas in the sharply decreasing region (20 to 22 h) it is controlled by substrate limitation, corresponding to the Monod-like factor in Equation (9)   The enzymatic adaptation factor EA is an empirical representation of the lag time, a period of adaptation to the culture environment where the microorganism produces new enzymatic machinery [38][39][40]. Based on the shape of experimental data, the following equation was proposed: where tlag (h) is the lag time experimentally observed. Figure 3b shows that EA is an increasing function of time, tending to 1 when t ≫ tlag. In analogy with n, r is a shape factor that describes the gradual transition from the lag phase to the active phase of growth. A higher value of r implies a steeper change of EA around t = tlag.
To illustrate the features of the proposed model, a representation of the dimensionless specific growth and production rates (µ/µmax and π/πmax) over time is depicted in Figure 4. The dynamic behavior of both variables is similar given the similarity of Equations (1)-(5). The specific rates achieve a maximum value in the active growth phase, and they are zero when t ≪ tlag and when the substrate is depleted. The shape of the curve is defined by three factors: in the increasing region (0 to 10 h in Figure 4), the dominant effect is enzyme activation EA (Equation (14)); in the slowly decreasing region (10 to 20 h), the rate is controlled by inhibition (Equation (11) or (12)), whereas in the sharply decreasing region (20 to 22 h) it is controlled by substrate limitation, corresponding to the Monod-like factor in Equation (9) or (10).

Model Parameter Identification
The system of kinetic equations for a single fermentation experiment included 24 parameters: five yield coefficients, five inhibition parameters, five growth/production rates, five Monod-like saturation constants, three shape factors, and one lag time. Due to a limited number of experimental data and to facilitate the identification procedure, a single value was adopted for the inhibition parameter (KIm) and the Monod saturation

Model Parameter Identification
The system of kinetic equations for a single fermentation experiment included 24 parameters: five yield coefficients, five inhibition parameters, five growth/production rates, five Monod-like saturation constants, three shape factors, and one lag time. Due to a limited number of experimental data and to facilitate the identification procedure, a single value was adopted for the inhibition parameter (K Im ) and the Monod saturation constant (K Sm ) of the four identified metabolites. Moreover, 10 parameters were fixed for all experiments: the shape factors, the yield coefficients, and the Monod saturation constants (K SX and K Sm ). For each fermentation, lag time was determined by graphical readout. This simplification of fixing parameters independent of operating conditions is supported by two assumptions often used in the literature: (1) metabolite production yields are constant and therefore independent of culture conditions [41] and (2) the saturation constant of the Monod model depends only on the nature of the substrate [33,38], which was the same in all experiments of this study.
The remaining group of seven parameters (µ max,X , π max,F , π max,A , π max,L , π max,E , K IX , K Im ) were identified for each fermentation of the experimental design by nonlinear regression. Here, the Levenberg-Marquardt algorithm [42] was used to minimize the sum of squares of the errors between experimental and predicted concentrations. However, since the ranges and the number of measurements were slightly different among the metabolites, the values compared in the least squares function were normalized by dividing by their maximum value and were weighted by the relevant number of experimental measurements.
The quality of the model representation was quantified with two error indicators, defined as follows: Root mean square error: Relative mean error (as a percentage): where N is the number of available measurements, C model and C exp are respectively the values of the concentration variables calculated with the model and measured experimentally.

Response Surface Model for Parameter Dependence on Fermentation Conditions
Nonlinear regression was performed to model the relationship between the seven parameters of the dynamic model specific to each experiment and the fermentation operating conditions-namely, temperature (T) and pH. The regression model had a similar form for all parameters, the logarithm of the parameter being expressed as a second-order polynomial with interaction: The regression coefficients (β) for all parameters depending on operating conditions (µ max,X , π max,F , π max,A , π max,L , π max,E , K IX , K Im ) were simultaneously computed by leastsquares optimization based on all available concentration measurements. In this way, the accuracy and standard errors of the coefficients were statistically acceptable, due to a large number of degrees of freedom: several hundreds of concentration data were used to estimate 42 coefficients. Initial guesses for these coefficients were obtained using Equation (17), and parameter values were determined separately for each experiment.
In this procedure, two sets of data from the experimental design were defined as indicated in Figure 2: eight calibration experiments, located in extreme positions of the experimental domain, used simultaneously for coefficients (β) estimation, and eight validation experiments, only used a posteriori to verify the accuracy of the complete dynamic model.

Model Parameter Identification
The values of the parameters that are independent of operating conditions, summarized in Table 1, were determined from the experimental data of experiment F10. This run was placed in a central position in the composite experimental design (T = 30 • C, pH = 8) (Figure 2). Monod saturation constants are usually difficult to determine from batch experiments because the number of measurements is typically very low in the substrate limitation zone. Saturation constants were thus fixed to a common value with a typical order of magnitude [43]. As for yields, they were found to differ from the theoretical ones defined through standard stoichiometric reactions of anabolism and catabolism. These differences can be due to other reactions involving the carbon substrate, whose products were not analytically measured and were not considered in the model. Table 1. Model parameters independent of operating conditions, determined from the experimental data of experiment F10 (T = 30 • C, pH = 8) with t lag = 10 h.

Parameter
Constant Value After fixing the parameters in Table 1 for the whole set of experiments, the group of seven adjustable parameters of the model (µ max,X , π max,F , π max,A , π max,L , π max,E , K IX , K Im ) were identified for each run by nonlinear regression.
The parameters obtained by this procedure are summarized in Table 2. Standard errors were computed from the variance-covariance matrix of the nonlinear optimization algorithm. These errors represented between 5% and 13% of the value of the identified parameters, a reasonable uncertainty level for a biological model. For the whole set of experiments, the prediction errors are reported in Appendix A Table A1. Except for some runs for variables S, F, and A, all RME were lower than 15%. Additionally, the average RMSE and RME values for each concentration were of the same order of magnitude as the experimental variability, here defined as the biological repeatability for run F01, for which three independent replicates were performed. These results validate the formulation and accuracy of the proposed model under the operating conditions included in the experimental design.
In the specific case of reference run F10, a comparison between the model simulation (using the corresponding parameters from Table 2) and experimental data is illustrated in Figure 5.  For the whole set of experiments, the prediction errors are reported in Appendix  Table A1. Except for some runs for variables S, F, and A, all RME were lower than 15%. Additionally, the average RMSE and RME values for each concentration were of the same order of magnitude as the experimental variability, here defined as the biological repeatability for run F01, for which three independent replicates were performed. These results validate the formulation and accuracy of the proposed model under the operating conditions included in the experimental design.
In the specific case of reference run F10, a comparison between the model simulation (using the corresponding parameters from Table 2) and experimental data is illustrated in Figure 5.  Three growth phases are apparent in Figure 5: a lag phase (phase 1, between 0 and 10 h); a phase of active growth, substrate consumption, and metabolite production (phase 2, between 10 h and 21 h); and a final phase where concentrations do not change over time, owing to the depletion of the carbon source or growth inhibition by metabolites (phase 3, after 21 h). Regarding culture volume evolution, as already mentioned, the discrete variations at regular intervals were due to sampling for analysis of the culture medium and the gradual increase was due to NaOH addition for pH control. One can also observe that the four metabolites were produced simultaneously, with no gap for the growth dynamics. The metabolites were thus primary end products generated during a single trophophase [44]. This justifies the choice of a global inhibitor concentration (Equation (13)), which included four correlated concentrations.
In consideration of the visual fit from Figure 5, the model representation is reasonably satisfactory. The most pronounced discrepancy between the model and experimental data appears for lactic acid, for which the model predicted a lower concentration before substrate depletion. This is related to a slightly underestimated yield factor Y L/S .

Response Surface Model for Parameter Dependence on Fermentation Conditions
Model parameters were expressed as a function of temperature and pH, according to the surface model (Equation (17)). The values of the β regression coefficients were adjusted globally using the whole set of calibration data.
The resulting response surfaces for the seven model parameters are plotted in Figure 6. For the five kinetic parameters, (i.e., the maximum specific growth and production rates), the response surfaces have the same convex shape, with a well-defined maximum value at intermediate T and pH conditions. These maxima likely indicate the optimal temperatures and pH for cellular growth, as well as the enzymatic activity catalyzing each of the reactions, leading to the production of the different metabolites ( Figure 1).
Concerning the inhibition concentrations, the response surface for K Im has a concave shape with a local minimum, whereas that of K IX resembles a saddle surface. For this latter case, the surface shape indicates that for every pH there is a T where K IX is minimal, and for every T there is a pH where K IX is maximal. Both K Im and K IX represent the combined effect of several inhibiting metabolites (Equations (11)-(13)) with potentially different inhibition mechanisms. Three growth phases are apparent in Figure 5: a lag phase (phase 1, between 0 and 10 h); a phase of active growth, substrate consumption, and metabolite production (phase 2, between 10 h and 21 h); and a final phase where concentrations do not change over time, owing to the depletion of the carbon source or growth inhibition by metabolites (phase 3, after 21 h). Regarding culture volume evolution, as already mentioned, the discrete variations at regular intervals were due to sampling for analysis of the culture medium and the gradual increase was due to NaOH addition for pH control. One can also observe that the four metabolites were produced simultaneously, with no gap for the growth dynamics. The metabolites were thus primary end products generated during a single trophophase [44]. This justifies the choice of a global inhibitor concentration (Equation (13)), which included four correlated concentrations.
In consideration of the visual fit from Figure 5, the model representation is reasonably satisfactory. The most pronounced discrepancy between the model and experimental data appears for lactic acid, for which the model predicted a lower concentration before substrate depletion. This is related to a slightly underestimated yield factor YL/S.

Response Surface Model for Parameter Dependence on Fermentation Conditions
Model parameters were expressed as a function of temperature and pH, according to the surface model (Equation (17)). The values of the β regression coefficients were adjusted globally using the whole set of calibration data.
The resulting response surfaces for the seven model parameters are plotted in Figure  6. For the five kinetic parameters, (i.e., the maximum specific growth and production rates), the response surfaces have the same convex shape, with a well-defined maximum value at intermediate T and pH conditions. These maxima likely indicate the optimal temperatures and pH for cellular growth, as well as the enzymatic activity catalyzing each of the reactions, leading to the production of the different metabolites (Figure 1).
Concerning the inhibition concentrations, the response surface for KIm has a concave shape with a local minimum, whereas that of KIX resembles a saddle surface. For this latter case, the surface shape indicates that for every pH there is a T where KIX is minimal, and for every T there is a pH where KIX is maximal. Both KIm and KIX represent the combined effect of several inhibiting metabolites (Equations (11)-(13)) with potentially different inhibition mechanisms.   Table A2. All coefficients in Equation (17) for each model parameter were significantly different from zero at a 0.05 level. A comparison between the parameter values determined for each experiment (Section 4.1) and the parameter values computed with Equation (17) (from globally adjusted β coefficients) is depicted in Appendix A Figure A1. The goodness of the fit was assessed through the coefficient of determination, R 2 . This coefficient is higher than 0.89 for six out of seven model parameters, which is a high threshold for biological data. In the case of K Im , only 66% of the variance of this parameter was explained by variables T and pH. The remaining 34% could be associated with inherent experimental variability and factors not included in the model, for instance transient variability of the inhibition and kinetics parameters and actual dependence of the fixed parameters (Table 1) with T and pH [45]. From a more general point of view, differences from experimental data could be due to features that were not represented by the mathematical model, such as population segregation, internal pH variability, and concentration gradients in the culture medium [46,47].

Model Validation
The ability of the dynamic model including the parameters calculated from operating conditions (Equation (17)) to predict data of independent experiments was assessed with a set of validation experiments.
A comparison between the average RMSE values obtained in Section 4.1 (determined for each experiment) and Section 4.2 (calculated from operating conditions) for calibration and validation sets is depicted in Figure 7. In most cases, RMSE values were higher than the corresponding experimental variabilities, indicating that more complex models could capture additional phenomena not included in the present model, such as dependence of yields, saturation constants, or lag time (Table 1) on operating conditions. As one might expect, RMSE was generally lower for the calibration experiments than for the validation experiments, not used for parameter determination. However, the relative difference remained small (less than 30%), indicating a satisfactory ability of the developed model to predict time evolution of the considered biomass, substrate, and metabolites under new conditions within the explored experimental range. For completeness, the final values of the regression coefficients of Equation (17) for the seven adjustable parameters of the dynamic model are reported in Appendix Table  A2. All coefficients in Equation (17) for each model parameter were significantly different from zero at a 0.05 level. A comparison between the parameter values determined for each experiment (Section 4.1) and the parameter values computed with Equation (17) (from globally adjusted β coefficients) is depicted in Appendix Figure A1. The goodness of the fit was assessed through the coefficient of determination, R 2 . This coefficient is higher than 0.89 for six out of seven model parameters, which is a high threshold for biological data. In the case of KIm, only 66% of the variance of this parameter was explained by variables T and pH. The remaining 34% could be associated with inherent experimental variability and factors not included in the model, for instance transient variability of the inhibition and kinetics parameters and actual dependence of the fixed parameters (Table 1) with T and pH [45]. From a more general point of view, differences from experimental data could be due to features that were not represented by the mathematical model, such as population segregation, internal pH variability, and concentration gradients in the culture medium [46,47].

Model Validation
The ability of the dynamic model including the parameters calculated from operating conditions (Equation (17)) to predict data of independent experiments was assessed with a set of validation experiments.
A comparison between the average RMSE values obtained in Section 4.1 (determined for each experiment) and Section 4.2 (calculated from operating conditions) for calibration and validation sets is depicted in Figure 7. In most cases, RMSE values were higher than the corresponding experimental variabilities, indicating that more complex models could capture additional phenomena not included in the present model, such as dependence of yields, saturation constants, or lag time (Table 1) on operating conditions. As one might expect, RMSE was generally lower for the calibration experiments than for the validation experiments, not used for parameter determination. However, the relative difference remained small (less than 30%), indicating a satisfactory ability of the developed model to predict time evolution of the considered biomass, substrate, and metabolites under new conditions within the explored experimental range.  (Table 2) and the response surface models (Table A2 and Equation (17)).
It also appears in Figure 7 that average RMSE values with parameters given by the response surface model (Table A3) are about 50% higher than with parameters determined separately for each experiment (Table A1), for both calibration and validation sets. This result could be expected since in the global calibration step, data from eight independent experiments were combined as a whole for the least squares estimation, with a detrimental effect on the individual representation of each experiment. However, results   (Table 2) and the response surface models (Table A2 and Equation (17)).
It also appears in Figure 7 that average RMSE values with parameters given by the response surface model (Table A3) are about 50% higher than with parameters determined separately for each experiment (Table A1), for both calibration and validation sets. This result could be expected since in the global calibration step, data from eight independent experiments were combined as a whole for the least squares estimation, with a detrimental effect on the individual representation of each experiment. However, results with the parameters calculated from operating conditions are the most useful in engineering purposes since they enable a quick prediction of growth and metabolites production dynamics, based on the selected combination of temperature and pH.
In light of this quantitative analysis, the prediction accuracy of the empirical dynamic model coupled to the regression model may be considered satisfactory within the operating domain covered in this study.

Model-Based Optimization of Fermentation Operating Conditions for Industrial Use
Optimal conditions for growth and metabolite production of C. maltaromaticum calculated using the developed model are summarized in Table 3. Two optimization criteria were considered: final concentrations and final productivities calculated for a 99.9% substrate consumption.
For a detailed representation of the evolution of final concentrations and productivities for biomass and metabolites with temperature and pH, the reader is referred to Appendix A Figures A2 and A3. As a general trend, the highest productivities were obtained around 35 • C and pH 7.5, although the exact optimal conditions depended on the considered metabolite (Table 3). No general trend was readily apparent for the maximization of the final concentrations.  These data can be useful in optimizing industrial processes involving the growth of C. maltaromaticum cells in a trehalose-based substrate. A first application consists of producing C. maltaromaticum concentrates, regardless of metabolite production. In this case two conditions of cultivation appear advisable: 20 • C and pH 7.8 to maximize concentration (227 mmol C ·L −1 ) or 33.5 • C and pH 7.5 in order to maximize productivity (7.49 mmol C ·L −1 ·h −1 ) and thus the biomass production per unit of time, at the expense of a 12% reduction of the final biomass concentration (199 mmol C ·L −1 ).
A second application deals with the development and parametrization of timetemperature integrators (TTI), labels in which a pH decline, associated with acids synthesis, entails an irreversible color change from green to red. Modulating the acidifying activity of C. maltaromaticum thus allows a reliable shelf-life estimation of different food products. Long shelf-lives can be tracked using TTI composed of concentrates exhibiting low acidifying activities (minimal production of total acids), while short shelf-lives can be tracked using concentrates exhibiting high acidifying activities. In the scenario of maximizing acidifying activity, the production of total acids must be favored, and thus fermentation should be carried out under two possible conditions: 37.0 • C and pH 6.0 to maximize their final concentration (433 mmol·L −1 ) or 35.5 • C and pH 7.7 to maximize their productivity (14.90 mmol·L −1 ·h −1 ). Under these conditions, the biomass production decreases respectively by 20% and 4% with respect to its optimal values. If the objective is, on the contrary, to minimize acidifying activity, two conditions can be envisaged to favor the lowest production of total acids: 37.0 • C and pH 9.5 for a final concentration of 352 mmol L −1 or 25.0 • C and pH 9.5 for a final productivity of 4.53 mmol·L −1 ·h −1 . Under these conditions, the mean biomass production would decrease respectively by 48% and 77% with respect to the maximal values.
Data from Table 3 show that the conditions to minimize the total acids concentration (37.0 • C and pH 9.5) coincide with those to maximize the ethanol concentration (the nonacidifying metabolite, 139 mmol L −1 ) and lead to a lactic acid concentration close to its minimal value (128 mmol L −1 versus the minimum around 120 mmol L −1 ). Conversely, when the production of total acids is maximized, the lactic acid concentration is also maximal (296 mmol L −1 ) and that of ethanol is close to its minimum (72 mmol L −1 versus 68 mmol L −1 ).
Furthermore, it should be noted that the condition 27 • C and pH 7.6 leads both to a good biomass productivity (6.93 mmol·L −1 ·h −1 versus the maximum 7.49 mmol·L −1 ·h −1 ) and a low total acids concentration (372 mmol L −1 versus the minimum 352 mmol L −1 ). Cultivation under this condition turns out be advantageous to ally a high biomass production and a relatively low total acidification.

Conclusions
The dynamic model developed in this study is able to predict with satisfactory accuracy the growth of C. maltaromaticum CNCM I-3298 (average error of 10%) as well as the conversion of trehalose into four primary metabolites (average error of 14%) under a wide range of conditions of temperature and pH. The interpolation capability of the model was verified with a set of eight independent validation experiments, for which the average relative error was 13%.
This model constitutes a useful tool for optimizing C. maltaromaticum cultures. Based on two easily controllable parameters, pH and temperature, it could be implemented in industrial applications of food technology to define optimal growth and metabolite production conditions with various objectives, such as the maximization of biomass for production of bacterium concentrates or the maximization or minimization of the acidifying activity. A typical operating condition for this bacterium could be, for instance, 30.0 • C and pH 7.0. If the goal is to produce bacterium concentrates, to maximize final biomass concentration, our results suggest that a quite different condition should be selected (20.0 • C and pH 7.8), while for maximum biomass productivity, 33.5 • C and pH 7.5 is most appropriate. Such results are quite difficult to anticipate from the qualitative knowledge of the bacterium alone, and a large number of time-consuming experiments would be required to locate these optimal conditions experimentally without constructing a dynamic model of the process.
The effort of developing the model is especially cost effective when a variety of scenarios are explored. If the goal is to develop time-temperature integrators (TTI) to track the cold-chain of food products, a set of labels with specific shelf-lives has to be designed for various target products. The range of desired shelf-lives can be as large as 1 to 30 days, requiring very different TTI designs. In a traditional approach, for each desired shelf-life duration, a range of factors such as the initial bacterium concentration and the buffer properties of the medium have to be explored in a series of relatively time-consuming experiments. In such an environment, temperature varies in an arbitrary but known way, and pH depends on the produced acids. The presented dynamic model can be extended to predict the moment when a specific amount of acids is produced, corresponding to the pH-induced color change of the TTI label and hence to the desired shelf-life. Modelbased design of the TTI labels is expected to be faster and more accurate than a trial and error procedure.
On a more fundamental level, further work is required to incorporate the effect of other culture parameters, such as aeration, nutrient concentrations, or the use of a different carbon source, which may modify growth kinetics and metabolite production. Additionally, it would be relevant to deepen the understanding of inhibition mechanisms of the metabolites to give more biological significance to the associated parameters in the model.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to an ongoing research project.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Concentration of species i (substrate, metabolite, biomass) in the culture medium (in the system of differential equations) I m Production inhibition factor of metabolite m I X Biomass growth inhibition factor K Im (mol·L −1 )

Nomenclature
Concentration for 50% production rate inhibition of metabolite m K IX (mol C ·L −1 ) Concentration for 50% growth rate inhibition of biomass K Sm (mol·L −1 ) Concentration of production rate saturation of metabolite m K SX (mol C ·L −1 ) Concentration of biomass growth rate saturation mol C Carbon-mol of biomass n Shape factor of the growth inhibition function p Shape factor of the production inhibition function pH Potential of hydrogen Q (L·h −1 ) Rate of base addition for pH control R Shape factor of the enzymatic activation function T (K) Temperature TTI Time-temperature indicator RMSE Root-mean square error RME Relative mean error SE Standard error t (h) Time t lag (h) Lag time V (L) Culture medium volume Y i/S (mol·mol −1 ) Yield of product i on substrate S Y X/S (mol C ·mol −1 ) Biomass yield on substrate S µ X (h −1 ) Specific growth rate µ max,X (h −1 ) Maximal specific growth rate π m (h −1 ) Specific production rate of metabolite m π max,m (h −1 ) Maximum specific production rate of metabolite m Appendix A.
Appendix A.1. Model Parameter Identification Table A1. Residual modelling error with model parameters determined for each experiment and summarized in Table 2.   Figure A1. Comparison between model parameters determined for each experiment (Table 2) and parameters computed with the response surface models (Equation (17) using coefficients in Table A2). Figure A1. Comparison between model parameters determined for each experiment (Table 2) and parameters computed with the response surface models (Equation (17) using coefficients in Table A2).  Figure A2. Evolution of final concentrations (left) and batch-average productivities (right) with temperature and pH for biomass, formic acid, and acetic acid. Figure A2. Evolution of final concentrations (left) and batch-average productivities (right) with temperature and pH for biomass, formic acid, and acetic acid.