Consideration of Maintenance in Wine Fermentation Modeling

We show that a simple model with a maintenance term can satisfactorily reproduce the simulations of several existing models of wine fermentation from the literature, as well as experimental data. The maintenance describes a consumption of the nitrogen that is not entirely converted into biomass. We show also that considering a maintenance term in the model is equivalent to writing a model with a variable yield that can be estimated from data.


Introduction
The overall principle of wine fermentation consists of the conversion of sugar into ethanol by yeast. It has been observed for a long time that nitrogen consumed during the yeast growth also plays an important role. The fermentation can be indeed modeled by a two-step process in which the yeast first grows on nitrogen as a limiting resource and then degrades the non-limiting sugar into ethanol and carbon dioxide. However, experimental observations have shown that the consumed nitrogen was not entirely converted into biomass. Several mathematical models were proposed to take these characteristics into consideration. For instance, in [1,2], the biomass growth follows a logistic law whose carrying capacity depends on the initial quantity of nitrogen. In [3], a model that distinguishes part of nitrogen used for yeast growth from another part responsible of the synthesis of proteins (hexose transporters [4]) was developed. Both models were calibrated with different sets of experimental data and provide satisfactory fitting. However, both models present some drawbacks. The dependency of the dynamics on the initial condition of the first model makes it sensitive to the precise knowledge of the initial quantity of nitrogen (that needs to be "memorized" in the dynamical equations of the model). Moreover, it does not allow consideration of non-batch operations or continuous addition of nitrogen, such as in [5] for instance. The second model relies on the knowledge of the time-varying concentration of transporters, which is in general not easily accessible for experimental measurements, and several assumptions were necessary to estimate it from biomass measurements.
The objective of the present work is to propose a new model that reconciles both approaches in a single one.
The observation of the ratio of produced biomass over nitrogen consumption along the whole fermentation, determined on experimental database or numerical simulations of models [1,3], shows that this ratio is non-constant and depends on the initial quantities. This highlights that the conversion of nitrogen into biomass can be viewed as a variable yield process. The experimental evidence that nitrogen is not entirely converted into biomass therefore advocates for the consideration of a maintenance term in the modeling (see, for instance, [6]), without necessarily requiring a detailed representation of the internal mechanism or cells.
Indeed, different mechanisms in the internal functioning of the cells have been investigated in the literature, particularly the role of carbohydrate accumulation [7][8][9], which could explain that the growth dynamics of yeast in wine fermentation does not follow the classical mass-balanced models [10,11]. However, the measurements of these biochemical compounds is experimentally very difficult and is almost impossible in an industrial framework.
The rationale of the results presented here is to test if the introduction of a maintenance term (see [12][13][14] or [15][16][17]) can improve wine fermentation modeling. One of the original features of the proposed approach is to view nitrogen consumption as a global consumption for growth by considering a variable yield. This allows us to avoid to consider a specific structure to model the maintenance. Thus, the purpose of the present work is to investigate the ability of a simpler model with a maintenance term to reproduce and predict wine fermentation kinetics.
Here, we propose a new modeling approach based on a maintenance term (which gives rise to a variable yield), a feature that has not been yet considered in the wine fermentation literature, to the best of our knowledge.
It focuses mainly on the new modeling of the growth of yeast on nitrogen. This new model was validated using both data generated by existing models (Section 4) and experimental data (Section 5).

The Proposed Model
We denote by N, S, E, CO 2 and X the concentrations of (total) nitrogen, sugar, ethanol, dioxide carbon and biomass, respectively. For simplicity, we derive a model under isothermal conditions. For the first step N → X (yeast growth on nitrogen), we propose the following equations where Y is the growth yield, µ the Contois growth function and m a maintenance function, which is positive for N > 0 and X > 0. We choose here a ratio-dependent kinetics function µ N to reproduce the observation that the growth is slowing down under an excess of yeast, with a Contois expression as in [3]. In the literature, the maintenance m is often considered as constant [12,13], which was validated in continuous culture (chemostat). In general, continuous cultures are intended to be operated at a stationary phase, very differently to batch-operating mode. However, as already investigated in [17], maintenance terms have to depend on the level of available resources; say, R (N here). In particular, constant maintenance in a batch model would imply dR dt < 0 when the resource is exhausted, i.e., R = 0, and thus R could take unrealistic negative values, as underlined in [14]. In [15,16], the maintenance is directly related to the microbial activity, which is stopped in absence of nutrients. This is why we consider a maintenance function proportional to the growth activity, with a factor that might depend on the nitrogen concentration (one may expect that it decreases when the substrate N becomes rare) where α is a positive function equal to zero for N = 0. Then, one can consider the function y defined as follows (1) and (2) can be rewritten equivalently as where the function y is playing the role of a variable yield. Identifying the function m or the function y is thus formally equivalent. However, we shall see in the next section that identifying the function y instead of m presents some practical advantages. For the second step S → E + CO 2 , we follow the model proposed in the literature [3] where µ S is a Monod function and ν E a function inhibited by the ethanol Inhibition by the consumption of sugar S by ethanol E has been reported many times in the literature [18][19][20][21][22]. The constant yield of production k of CO 2 and consumption of S follows a mass balance assumption, verified experimentally [23], that can be determined using thermodynamics considerations [24]. Note that this model can be extended to anisothermal conditions, considering that the maximal specific rate parameters µ max N , µ max S and affinity constants K S , K E are temperature dependent, as in [3].

Calibration of the Model
From model Equation (1), the parameters of the function µ N can be identified independently of the yield and maintenance terms. To validate the hypothesis of ratio dependency of the function µ N , one can first use experimental data to plot the slope of the logarithm of X versus the ratio r = N/X and check if it qualitatively follows a function of the form µ(r) = µ max N r K N + r A classical least-square method can be applied to fit parameters µ max N , K N on the data. Alternatively, one can plot the inverse of the slope of the logarithm of X versus the inverse of the ratio r to check if it qualitatively follows a linear dependency, as obtained from However, for the accurate identification of the parameters µ max N , K N , a linear regression on Equation (8) is expected to be less reliable than a non-linear least-square optimization of the solution X(·) of (3), because data might be too far to be uniformly distributed.
Note from Equations (1) and (2) that one has lim t→+∞ N(t) = 0 (because the derivative of N cannot vanish when N is not exhausted). In absence of the maintenance term m, one gets dX dt + Y dN dt = 0 which implies that one should have To test the validity of the model with maintenance, one can plot from experimental data the ratio for different values of N(0) to check that it is not constant. If this is the case, one can then look at identifying a non-constant function y. For this purpose, we write from Equations (3) and (4) and as t → N(t) is a monotone-decreasing function, one can make the change of variable n = N(t) in this last integral to obtain Therefore, if one fits a differential function f such that f (0) = 0 that satisfies for experimental data with different values of N(0), then one simply gets y = f . Let us underline that identifying the function y in this way can be achieved independently of the knowledge of the kinetics µ N , differently to the function m, which clearly presents some robustness advantages. Once the function µ N is identified, the maintenance function can then be determined as where Y = y(0) (to fulfill α(0) = 0). For model Equations (5) and (6), the coefficient k is kept from the literature, and the parameters β, µ max S , K S , K E are identified (with a least-square method) from experimental data of CO 2 production rate.

Validation of the Model on Synthetic Data
We have used synthetic data generated by models of the literature that were previously validated on experimental data [1,3] for a range of initial conditions and operating conditions.
Fitting comparisons of the proposed model with the different data sets are reported in Section 6.

Validation on Simulations of a Model with Transporter
We have considered the model with transporters developed in [3], which is more complex with two additional state variables: the concentrations of hexose transporters and the nitrogen dedicated to these transporters. Data were generated by simulating this model with the parameters given in [3] and operating conditions given in Table 1. This model explicitly distinguishes two forms of nitrogen, one available for the yeast N X and the other one N tr for the transporters. To compare with the variable N of our model, we have considered the total nitrogen N = N X + N tr .

Estimation of the Contois Function
We have used a non-linear least-square method based on a Newton algorithm with a finite difference approximation of the Jacobian matrix (function leastsq of scilab).  Table 2.

Estimation of the Variable Yield Function
On Figure 2, data X(T) − X(0) versus N(0) from the model with transporters were plotted for T = 350 h (we have checked that N is quasi-null at T and that X no longer increases after T). One can see that the points are aligned. However, the line that passes through these points does not touch 0, which is not possible for a constant yield (for a constant yield, the points have to be aligned on a line that passes through 0, because when N(0) = 0, there is no biomass production).
Then, we fitted a C 2 function f such that f (0) = 0 with the following expression whose parameters are given in Table 3. The calibration of the parameters a, b of the function f was performed with a linear regression (function reglin of scilab). Then, we obtain the variable yield function y as the C 1 function and the function α, which describes the maintenance as which are both depicted on Figure 3. Note that the model with transporters was validated only for N(0) in the interval [0.071, 0.57] g·L −1 , and that we have no a priori information about the behavior of the yield for values of N(0) smaller than 0.071 g·L −1 . The threshold parameter N † was simply chosen so that the simulations of the variables X and N of the model (3) and (4) were the closest to the ones of the transporter model.

Estimation of the Other Parameters and Comparison of the Models
For the model of the second step S → E + CO 2 , the stoichiometric parameter k was taken from the literature, while the other parameters β, µ max S , K S , K E were estimated with a least-square optimization on the CO 2 chronicles only (the CO 2 production rate being a variable that is usually measured in experiments), starting from values in [3]. Values are given in Table 4. Here, we also used a non-linear least-square method based on a Newton algorithm with a finite difference approximation of the Jacobian matrix (function leastsq of scilab). All data were re-normalized to 1 (i.e., for each variable, the figures were divided by the largest one).
Finally, we present on Figures 4-6 simulations of the new model for three largely different initial values of nitrogen from 0.170 g·L −1 to 0.567 g·L −1 . The evolution of the ethanol concentration E has not been reproduced as it is proportional to the CO 2 concentration.
These simulations show the ability of the new model to reproduce, with a single set of parameters, close simulations to the model with transporters, in terms of production of biomass and dioxide carbon, estimation of the peak of the CO 2 production rate and depletion of (total) nitrogen and sugar.

Validation on the SOFA Model
The model proposed in [1] does not explicitly consider transporters with an additional state variable as the previous model, and instead presents a more sophisticated expression of the dynamics that depend on the initial condition, with an additional latency term at the beginning of the simulations.
Differently to the previous model, which is built as a "mass-balanced" model, this one relies on an empirical dynamics of logistic shape for the biomass growth, with some parameters that depend on the initial concentration of nitrogen N(0), instead of the twodimensional model (3) and (4). Therefore, this is not a Markovian model. It has been validated on different operating conditions, and has been encoded into the SOFA software exploited for decision-making [2]. We launched simulations of this model for the same operating conditions than for the previous model (Table 1). Although simulations look qualitatively similar, they do not overlap, especially for the biomass chronicle. This could be explained by the fact that this model is intended to predict a number of cells and not a precise biomass (an average number of 4.15 × 10 9 cells for one g of biomass was used to have X expressed in g·L −1 as for the previous model). We proceeded to a new validation of our model on these data. Figure 7 shows that the data N X , dX dt X do not precisely follow the graph of a function (this is most probably due to the fact that the model is not Markovian). Indeed, this happens mainly for the large value N 0 of the initial nitrogen. We believe that this could be explained by the dynamics of the biomass X of this model, which is a logistic law with a carrying capacity given by an heuristic expression that depends on N 0 , and not dynamics coupled with the dynamics of N (indeed the interval of tested values of N 0 might be larger than the validity of this model). However, we have fitted the graph of a Contois function to these data with the parameters given in Table 5, which was able to satisfactorily reproduce the trajectories of the model for a large amplitude of values of N 0 , as we shall see later on.

Estimation of the Contois Function
As for the previous model, we used a non-linear least-square method based on a Newton algorithm with a finite difference approximation of the Jacobian matrix (function leastsq of scilab). As one can see in Table 5, the values of µ max N and K N are significantly larger and smaller, respectively, than in Table 2, which is consistent with the observation that this model predicts a faster convergence of the biomass to its maximal value, despite the latency term (compare Figures 4-6 with Figures 8-10).

Estimation of the Variable Yield Function
Data X(T) − X(0) from the simulation of the SOFA model were plotted on Figure 11 at T = 350 h, for different values of N(0) in the interval [0.071, 0.57] g·L −1 (here, we also checked that the fermentation was quasi-ended at T). One can see that the points follow an increasing concave curve and further increase very slowly, quite differently to the model with transporters (see Figure 2).
and parameters a, b, N † given in Table 6.
Parameters a and b were determined with a linear regression (function reglin of scilab). Then, we obtain the expression of the variable yield function as well as the function α whose graphs are drawn on Figure 12.

Estimation of the Other Parameters and Comparison of the Models
For the second step, the same stochiometric parameter k was taken for the literature, and the other parameters β, µ max S , K S , K E were estimated with a least-square optimization on the CO 2 chronicles only, as for data generated by the model with transporters (see Table 7).  Figures 8-10 show the comparison between the SOFA model and our calibrated model for the same initial condition than for the former comparison with the model with transporters. Here also, we see that the proposed model reproduces quite faithfully the simulations of the SOFA model, with the advantage of being a simpler Markovian model. Indeed, the difference between the model with transporters and the SOFA model can be translated into different maintenance terms (see Figures 3 and 12): for large values of nitrogen, the model with transporters behaves like a model with a maintenance proportional to the growth, while the SOFA model amounts to have a strongly increasing maintenance. Recall that the simulations for the largest value of N(0) showed the most differences between these two models (for N(0) = 0.567 g·L −1 , the model with transporters predicts a biomass production of 5.11 g·L −1 , while the SOFA model predicts 3.88 g·L −1 ; see Figures 6 and 10). While the model with transporters was validated experimentally for N(0) in the interval [0.170, 0.567] g·L −1 , we believe the validation of the SOFA model for initial concentrations of nitrogen larger than 0.4 g·L −1 might need to be revisited (although our model once calibrated is able to reproduce the SOFA simulations).

Calibration of the Model on Real Data
We considered data from experiments conducted at SPO Lab (INRAE, Montpellier, France) in 2004, that were used to calibrate the model with transporters and the SOFA model (see [1,3]). The data consisted of a set of three experiments with the same operating conditions given in Table 1  We first calibrated a function f (·) to the data (N(0), X(T) − X(0)), with the same expression as in Section 4.2, to determine a yield function y(·) (see Figure 13), using a linear regression to estimate parameters a and b. Figure 13. Results of the fitting of the function f on the experimental data (left) and of the corresponding variable yield function y (right).
As we do not have measurements of N over the time, we cannot estimate the Contois parameters independently of the CO 2 measurements, as we did with the synthetic data. All the parameters of the model were fitted simultaneously with a least-square method (values are given in Table 8), except for the sugar-conversion yield, for which we have used the value of the literature k = 2.17, as before.
The non-linear least-square method uses a Newton algorithm with a finite difference approximation of the Jacobian matrix (function leastsq of scilab), and the data set was re-normalized to the maximal value of 1.  Figures 14-16 show the results of the fitting for the three experiments. One can appreciate the goodness of fit for a unique set of parameters. In particular, the production of biomass and CO 2 , as well as the height and date of the peak of dCO 2 /dt, are well predicted with this model.

Fitting Comparisons
For the calibration of the variable yield function on both synthetic and experimental data (Sections 4 and 5), we have used a linear regression (function reglin of scilab) for the determination of parameters a, b of the function f (for the model with transporter) and f (for the SOFA model and experimental data). The residual error is given in Table 9. This shows that the model with transporters behave very closely to a variable yield model. The fitting performances for the SOFA model and experimental data are more difficult to interpret, because the validity of the SOFA model for the large range of initial concentrations of nitrogen we considered is questionable, and the quantity of experimental data is quite poor compared to the synthetic data.
For the synthetic data, the calibration of the growth characteristics (parameters µ max N , K N of the Contois function) was performed first, independently of the CO 2 data. Then, parameters for the second step (parameters k, β, µ max S , K S , K E for the CO 2 production) were calibrated. In both cases, a non-linear least-square method based on a Newton algorithm with a finite difference approximation of the Jacobian matrix (function leastsq of scilab) was used. Table 10 shows a good fitting quality. We recall that for experimental data, we do not have measurement of N over time, so it was not possible to estimate the growth function independently of the CO 2 measurements. The estimation of all the parameters was made on the CO 2 measurements only. We have used the same non-linear least-square method, with data re-normalized to 1 (i.e., the figures were divided by the largest one), so that all points have equal weight in the criterion. The errors shows a good fitting of the CO 2 curves with the model with maintenance.

Conclusions
In this work, we demonstrated that the consideration of a maintenance term, or equivalently, a variable yield, in wine fermenting modeling can satisfactorily replace more sophisticated models with a simpler structure. Indeed, the effects of the underlying mechanisms of transporters or carbohydrate accumulation, which are difficult to capture experimentally, are somehow encoded into a maintenance term, and are translated into a variable yield between biomass and nitrogen. We showed that this variable yield, as a function of the nitrogen concentration, can be estimated from experimental data of biomass growth and nitrogen depletion, without the need to measure internal compounds. This consideration brings a flexibility to suit different kind of models or experimental data (once calibrated) with a single common structure, that could correspond to different operating conditions or hypotheses in wine fermentation. This new approach provides new perspectives of control of fermentation with nitrogen addition, based on a simple Markovian model, as well as model extensions with aromatic compounds [25] or multistrains [26].