Monod model is insufficient to explain biomass growth in nitrogen-limited yeast fermentation

The yeast Saccharomyces cerevisiae is an essential microorganism in food biotechnology; particularly, in wine and beer making. During wine fermentation, yeasts transform sugars present in the grape juice into ethanol and carbon dioxide. The process occurs in batch conditions and is, for the most part, an anaerobic process. Previous studies linked limited-nitrogen conditions with problematic fermentations, with negative consequences for the performance of the process and the quality of the final product. It is, therefore, of the highest interest to anticipate such problems through mathematical models. Here we propose a model to explain fermentations under nitrogen-limited anaerobic conditions. We separated the biomass formation into two phases: growth and carbohydrate accumulation. Growth was modelled using the well-known Monod equation while carbohydrate accumulation was modelled by an empirical function, analogous to a proportional controller activated by the limitation of available nitrogen. We also proposed to formulate the fermentation rate as a function of the total protein content when relevant data are available. The final model was used to successfully explain experiments taken from the literature, performed under normal and nitrogen-limited conditions. Our results revealed that Monod model is insufficient to explain biomass formation kinetics in nitrogen-limited fermentations of S. cerevisiae. The goodness-of-fit of the herewith proposed model is superior to that of previously published models, offering the means to predict, and thus control fermentations. Importance Problematic fermentations still occur in the winemaking industrial practise. Problems include sluggish rates of fermentation, which have been linked to insufficient levels of assimilable nitrogen. Data and relevant models can help anticipate poor fermentation performance. In this work, we proposed a model to predict biomass growth and fermentation rate under nitrogen-limited conditions and tested its performance with previously published experimental data. Our results show that the well-known Monod equation does not suffice to explain biomass formation.


52
The yeast species Saccharomyces cerevisiae is the best-studied eukaryote. S. cerevisiae presents 53 unique characteristics such as its fermentation capacity and its ability to withstand adverse 54 conditions of osmolarity and low pH. These desirable properties have made of S. cerevisiae one 55 of the workhorses of bio-industries including food, beverage -especially wine-and biofuel 56 production industries (1). 57 In winemaking, the fermentation turns grape must into an alcoholic beverage. Wine production 58 occurs in a closed system (i.e. in batch conditions) and is, for the most part, an anaerobic process. 59 During fermentation, yeasts transform sugars present in the must into ethanol and carbon 60 dioxide. A complete alcoholic fermentation is achieved when the residual fermentable sugar is 61 less than 2 g/L. Despite improvements in fermentation control, problematic fermentations such 62 as, stuck -with a higher than desired sugar residual-and sluggish -unusually long-fermentations 63 still occur in real practice. sluggish fermentations (2). A minimum of 120-140 mg/L of assimilable nitrogen is required to 66 achieve a standard fermentation rate, while nitrogen content in grape juice may be as low as 60 67 mg/L (3). Therefore, it has become common practice to supplement nitrogen-deficient musts 68 with diammonium phosphate (4). Nevertheless, both excessive or insufficient nitrogen can lead 69 to the production of undesired metabolites affecting the organoleptic properties of wine (4, 5). 70 The need to predict and control wine fermentation has motivated a quest for mechanistic models 71 of alcoholic fermentation (see the reviews by (6, 7)). Previous studies proposed various 72 alternatives that differ on the way they explain biomass formation and the relation of cell mass 73 (or sometimes cell numbers) to fermentation rate. Dynamic kinetic models incorporate the 74 Monod equation to describe biomass growth (8-10). The estimation of the Monod parameters 75 requires nitrogen and sugars uptake data throughout the fermentation. Otherwise, lack of 76 identifiability of the parameters will limit the predictive capability of the model, and thus logistic 77 functions may be more suitable (11,12). Both formulations, Monod or logistic, achieve a 78 maximum biomass value which can not be exceeded. 79 However, (13) showed that, for nitrogen-limited anaerobic fermentations of S. cerevisiae, there is 80 a substantial increment in biomass content after the depletion of ammonia. Besides, 81 measurements of protein and messenger RNA (mRNA) contents revealed that their 82 concentrations stay stable while the concentration of carbohydrates increases. Similar effects 83 have been detected in baker's yeast exposed to nitrogen starvation but not in the carbon starved 84 cells; see (14) for a fed-batch aerobic example, or (15) for a chemostat anaerobic example. These 85 data would imply that, at least for ammonium deficient scenarios, the use of the Monod model 86 would be insufficient. 87 where X is the biomass (g/L). 132 To model growth associated with cell division we used a Monod type kinetics (25) where 133 nitrogen sources were considered the limiting nutrient: 134 where is the maximum specific growth rate, is extracellular YAN (g/L), is the 135 Monod equation parameter and lag is a function of time (t) representing the lag phase. The lag 136 phase is typically a period of adjustment in which a given microbial population adapts to a new 137 medium before it starts growing exponentially. To model the lag we used the model proposed by 138 (26): 139 where 0 is an estimated parameter and is the maximum growth rate associated with . 140 The secondary increase in biomass concentration ( ) is activated by the of decrease the 141 concentration of YAN: 142 includes a death rate proportional to the extracellular ethanol concentration. 157 2.2 Modeling fermentation rate and production of extracellular metabolites.

158
Regarding fermentation rate, we proposed two models with different underlying hypotheses: 159 • M : assuming that fermentation rate is proportional to the total protein content ( ), 160 • M : assuming that fermentation rate is proportional to the total biomass ( ) 161 The uptake of glucose reads as follows: 162 where is the expression governing glucose transport. A similar expression was also included 163 for fructose ( ) transport ( ): 164 Following (27), we used Micheaelis-Menten type kinetics kinetics coupled to ethanol inhibition 165 to model hexose transport: is the maximum rate of glucose transport, is the Michaelis-Menten constant 167 and ℎ ethanol inhibition. Ethanol has been reported as a non-competitive inhibitor (28) of 168 hexose transport. Here, we modeled its effect as follows (27): 169 where ℎ is the extracellular ethanol concentration and defines the strength of the inhibitory 170

effect. 171
The excretion of extracellular metabolites was assumed to be proportional to the consumption of 172 hexoses. In the case of ethanol, this is described by: 173 where / is the ethanol yield produced from glucose and fructose. 174 2.3 Calibration of candidate models.

175
The Table 1 summarizes the main characteristics of the three candidate models. Additionally, a 176 detailed description of models is given Supplementary Text S1. 177 The structural identifiability analysis of the models revealed that it is possible to uniquely 178 estimate parameter values provided the dynamics of biomass, glucose uptake, nitrogen sources 179 uptake and ethanol are measured throughout the fermentation. Remarkably, it is not necessary to 180 details in the Supplementary Text S1. 182 Candidate models were calibrated by data fitting using data-sets available in the literature. In 183 particular we used data-sets taken from (13) (Exp ) and (16) (Exp ). The data-set Exp included 184 one experiment with a nitrogen sufficient condition. The data-set Exp contained two 185 experiments, with high (Exp ↑ N) and low (Exp ↓ N) nitrogen concentrations (50 and 300mg 186 YAN). 187 We performed 10 runs of the parameter estimation problem for each case, so to guarantee 188 convergence to the best possible solution. Since we used a metaheuristic for the optimisation of 189 parameters, we obtained a distribution of values. We selected the best fit for ulterior analyses. 190 Bounds for the estimated parameters values are given in supplementary text S1. Particularly, 191 parameters related to biomass composition were chosen such that carbohydrate content does not 192 decrease after nitrogen depletion (i.e. > ). 193 The parameters recovered with the estimation procedure are given in Supplementary Text S1. show that all models could explain the data within a 9% normalised root mean square error. 197 However, the performance of different models differed substantially. The model proposed in this 198 work, M , which explains total biomass dynamics using two terms (due to nitrogen and 199 carbohydrates) produced NMRSE values that were 3.9 and 1.8 times lower than models M and 200 the corrected Akaike Information Criterion (AIC) to rule out any possible over-fitting. We 203 computed AIC scores excluding the residuals corresponding to protein, carbohydrates and 204 mRNAas these variables are not included in M and Mand considering the total number of 205 parameters of each model (including , , and ). Table 2   cases, the proposed model could recover the increase of biomass produced after YAN depletion. 222 The improvement was more notorious in the Exp data-set.   indicating that protein content is a helpful biomarker while modeling fermentation rate across 249 multiple experimental conditions. 250

251
In this study, we proposed two candidate models to explain nitrogen-limited fermentations of S. 252 cerevisiae in industrially controlled conditions. We proposed two modifications to the standard 253 modelling approaches: 1) the biomass growth accounts for protein and carbohydrates, and 2) 254 fermentation rate is proportional to the total protein content. We reconciled candidate models to 255 previously published data. 256 Our results revealed that Monod model is insufficient to explain biomass formation kinetics in 257 nitrogen-limited fermentations of S. cerevisisae. This is so because, during exponential growth, 258 the biomass composition seems to remain unaltered with mRNA and protein comprising an 259 important percentage of the biomass. However, this is not the case during secondary growth 260 phase where fades and increases, as shown in (13)    qualitative point of view, this model corresponded to a significant improvement over past 300 models. Noticeably, the hexose transport functions considered the role of the initial nitrogen 301 concentration and temperature. Additionally, (16) showed that the protein content of biomass is 302 heavily dependent on the initial nitrogen concentration which supports the notion that hexose 303 transport is likely to be better described by total protein content rather than biomass per se. Our 304 model would support this hypothesis. Results obtained by M model are superior to those 305 obtained with the M model, at least, for those cases in which initial nitrogen is restricted. 306 Here we proposed an alternative modelling strategy which using a minimal number of 307 parameters, that can be estimated from easily measured data, provides a considerable published models. It should be noted that the model can be easily applied to controlled 310 fermentations -with a single starter-provided the grape must is well characterized, i.e. YAN and 311 carbon sources are measured beforehand and their values are followed throughout fermentation, 312 and biomass growth and desired products are monitored over time.

315
We considered data-sets from five different experiments as published in the literature. The first 316 experiment (Exp1) data set was obtained from (13). The authors conducted a nitrogen-limited 317 (165 mg NH 3 and 50g of glucose) batch fermentation while measuring biomass, protein content 318 of biomass along with several extracellular metabolites (glycerol, ethanol and ammonia). We 319 also considered four data-sets from (9). The authors performed various experiments at different 320 conditions. In this work we denote Exp2 as that performed at 11 C, with high-sugar and low-321 nitrogen; Exp3, the one performed at 15 C, with normal sugar and low-nitrogen; Exp4, the one 322 performed at 30 C with normal sugar and normal nitrogen; and Exp5, the one performed at 35 C 323 with high-sugar and low-nitrogen. However for some experimental set-ups protein was not measured. Therefore we performed a 335 structural identifiability analysis ((29, 30)) to assess if all parameters of the model can be 336 estimated even if protein measurements are not available. 337 We performed the analysis using GenSSI2 toolbox (31). The toolbox uses symbolic 338 manipulation to implement the so-called generating series of the model and to compute 339 identifiability tableaus. 340

Parameter estimation 341
The aim of parameter estimation is to compute the unknown parameters -growth related 342 constants and kinetic parameters -that minimize the weighted least squares function which 343 provides a measure of the distance between data and model predictions (32). In the absence of 344 error estimates for the experimental data a fixed value for each state, the maximum experimental 345 data, was assumed in order to normalize the least squares residuals: 346 In our particular case, we used five different experiments: the first, taken from (13), and second 350 to fifth, taken from (9). All experiments provided time series data of the biomass ( ), glucose 351 uptake ( ), amonia consumption ( ) and ethanol production ℎ. To avoid the risk of premature convergence by the optimization routine while searching for the 356 optimal parameter set, the parameter estimation procedure was repeated 10 times for each 357 candidate model starting from different initial guesses. The models were ranked by , with the best approximating model being the one with the 365 lowest value. The minimum AIC value (AIC ) was used to re-scale the Akaike's 366 information criterion. The re-scaled value = | − | was used to assess the 367 relative merit of each model. Models such as ≤ 2 have substantial support, models for which 368

Numerical methods and tools 370
The parameter estimation and model selection were implemented in the AMIGO2 toolbox (35). 371 The system of ODEs was compiled to speed up calculations and solved using the CVODES 372 solver (36), a variable-step, variable-order Adams-Bashforth-Moulton method. The parameter 373 estimation problem was solved using a hybrid meta-heuristic, the enhanced scatter search 374 method (eSS, (37)).